A periodic configuration of the Kepler-25 planetary system?
Cezary Migaszewski Krzysztof Gozdziewski




firstpage–lastpage 2018
[
    [
    Accepted .... Received ...; in original form ...
====================================================



[figure]labelfont=bf,font=small



We study proximity of the Kepler-25 planetary system to a periodic configuration, which is known to be the final state of a system that undergoes smooth migration resulting from the planet-disc interaction. We show that the system is close to the periodic configuration of 2:1 mean motion resonance (MMR) what indicates that its past migration was neither disturbed significantly by turbulence in the disc nor the orbits were perturbed by planetesimals that left after the disc dispersal. We show that, because of the TTV model degeneracy, a periodic configuration is difficult to be found when the standard modelling of the transit timing variations (TTVs) is used. The TTV signal of a periodic configuration (with anti-aligned apsidal lines) may be misinterpreted as an aligned non-resonant system. We demonstrate that the standard MCMC modelling of the Kepler-25 TTVs is very sensitive to an a priori information on the eccentricities (prior probability distributions). Wide priors (of the order of the ones typically used in the literature) result in favouring the aligned non-resonant configurations with small planets' masses and moderate eccentricities, while for the narrower priors the most likely are the anti-aligned resonant systems with larger masses and very low eccentricities. 




Planetary systems – planets and satellites: dynamical evolution and stability – planet-disc interactions






§ INTRODUCTION


^⋆Email: migaszewski@umk.pl (CM), chris@umk.pl (KG)

It is known that a periodic configuration (in a reference frame co-rotating with the apsidal lines) is one of the possible outcomes of the smooth disc-induced migration of a two-planet system <cit.>. The migration can be alternately convergent or divergent during the whole lifetime of the disc. Depending on particular history of the migration, the final period ratio P_2/P_1 (where P_1, P_2 are the orbital periods of the inner and the outer planet, respectively) may be close to a nominal value of a particular mean motion resonance (e.g., 2/1, 3/2, 4/3, etc.) or shifted away from such a value. The proximity of the final system to the periodic configuration is, however, not related to the final period ratio. The system may have P_2/P_1 very distant from the resonant value but still be a periodic configuration, i.e., be resonant in terms of the resonant angles librations.

Nevertheless, as mentioned above, the periodic configuration is not the only possible result of the migration. The system can deviate from periodicity if the migration is too rapid (i.e., the evolution is non-adiabatic) or if the system passes through the resonance. The latter may happen if the migration is divergent or if the resonance capture is only temporary, i.e., librations around the periodic configuration (an equilibrium in the averaged model of the resonance) are overstable <cit.>. It is also possible that the amplitude of the librations saturates at a non-zero value, resulting in a resonant system that is shifted away from the periodic configuration <cit.>. Although, in general, the migration does not necessarily lead to the periodic configuration, a system that is close-to-periodic has been likely formed on the way of migration.

This work is related to the problem of explaining the observed period ratio distribution of the KEPLER systems <cit.>. Only small fraction of multi-planet systems have P_2/P_1 close to resonant values. There are several explanations of how a given resonant pair moved away from the resonance. They may be divided into two groups. In the scenarios from the first group the system is being moved away from the resonance because of perturbations resulting from remnant planetesimals <cit.> or the interaction with a turbulent disc <cit.>. In the scenarios from the second group the migration is smooth (i.e., not disturbed by the forces mentioned above) but not necessarily convergent during the whole evolution of the system. The divergence of the migration may result from particular physical conditions in an evolving disc <cit.> or from the tidal star-planet interaction <cit.>. The latter mechanism acts, however, only for very short-period planets (P ∼ 1day). Another mechanism, that can be counted to this group has been recently proposed <cit.>. They show that the equilibrium values of the period ratio of two short-period planets may differ significantly from the nominal resonant values if a disc in which they migrated has a small aspect ratio and is significantly flared.
<cit.> studied the Kepler-25 system with two transiting planets <cit.> in this context, and showed that its period ratio of 2.039 may result from the migration in the disc of properties suggested by <cit.>.

The first group scenarios result in not only the P_2/P_1 deviation from the resonant value but also in the deviation of the system from the periodic configuration. As discussed in <cit.> stochastic forces acting on a planet affects both the semi-major axis a and the eccentricity e in similar magnitude. That means that if the Kepler-25 with P_2/P_1 ≈ 2.039 was shifted away from the 2:1 MMR by such forces, the variation of the eccentricities should be ≳ 0.01, which is far from the periodic system, as e_1 ∼ 10^-3 and e_2 ∼ 10^-4 for the periodic configuration of P_2/P_1 = 2.039 and the planets' masses in super-Earths regime <cit.>. [We use the indices of 1 and 2 for the inner and the outer planet, instead of b and c.] That means that the Kepler-25 system is a good tester of the migration as a formation mechanism of configurations with P_2/P_1 relatively close to (but not exactly at) the resonant values.

This work is organized as follows. In Section 2 we study branches of periodic configurations as a function of the planets' masses m_1, m_2 (for the inner and the outer planets, respectively) and the period ratio. As the eccentricities e_1, e_2 depend on m_1, m_2, P_2/P_1 and the TTV amplitudes depend on e_1, e_2, we show that it is possible to constrain the masses of the Kepler-25 system, when assuming a periodic configuration of this system. The periodic configuration fitting procedure is presented in Section 3. In the next section we try to verify whether or not a system which is far from a periodic configuration may be misinterpreted as a periodic system, when applying the procedure explained in Section 3. On the other hand, in Section 5 we try to find out if a periodic configuration may be misinterpreted as far form periodic when applying the standard TTV fitting procedure. The last section is devoted to summary and conclusions.




§ MIGRATION AND PERIODIC CONFIGURATIONS


Families of periodic orbits of two-planet systems has been widely studied for different MMRs and in a wide range of planets' masses <cit.>. Their connection with the migration has been also pointed out <cit.>. Father in this work we show the connection between the periodic configurations and the TTVs.

When using the averaging approach to the resonant two-planet system <cit.>, a periodic configuration corresponds to an equilibrium in a reference frame co-rotating with the apsidal lines. For a stable equilibrium (which is the case we are interested in) the eccentricities, semi-major axes, the difference of the longitudes of pericentres (Δϖ = ϖ_1 - ϖ_2, where ϖ_i is the longitude of the i-th planet's pericentre) and the resonant angles are constant. The individual ϖ_i, on contrary, vary linearly in time. 
The so-called free eccentricities of a system in an equilibrium equal zero and the apsidal lines rotate with a period that equals the so-called super-period <cit.>, T_q = |q/P_1 - (q+1)/P_2|^-1, for the (q+1):q resonance.

From the observational point of view, a uniformly rotating orbit of given (and fixed) a and e should lead to a periodic signal in TTV, as an actual orientation of the orbit determines whether the transit occurs earlier or later than it stems from the Keplerian motion of the planet. As the true anomaly ν may be expressed for e ≪ 1 by the mean anomaly  as ν≈ + 2 e sin <cit.>, the TTV should be sinusoidal (for low e) with a semi-amplitude A that depends on e through the relation A/P = (2 e)/(2π). If the orbits of a periodic two-planet configuration are anti-aligned, Δϖ = π, the TTV signals for the planets should be in anti-phase. If the orbits are aligned, Δϖ = 0, the signals are in phase, while for Δϖ different from 0 or π <cit.> the difference in phases of the signals are neither 0 nor π.



As the amplitudes of TTV depend on the eccentricities, a natural question arises on how the eccentricities of periodic configurations depend on planets' masses. Figure <ref> presents branches (families) of periodic orbits obtained for different m_1, m_2 in certain range of P_2/P_1 > 2. The procedure of finding periodic orbits is the following <cit.>. For given masses we search for such p⃗ = (a_1cos_1, a_1sin_1, a_2cos_2, a_2sin_2, e_1cosΔϖ, e_1sinΔϖ, e_2) that satisfies δ = 0, where δ≡p⃗(t = T) - p⃗(t = 0), and T is the period.

For 2:1 MMR the period T corresponds to one revolution of the outer planet and two revolutions of the inner planet. The stable equilibrium (when we use the averaged model of the resonance) of 2:1 MMR in a regime of small eccentricities exist for Δϖ = π and values of the resonant angles depend on P_2/P_1. For P_2/P_1 > 2 the angles defined as ϕ_1 ≡λ_1 - 2 λ_2 + ϖ_1 and ϕ_2 ≡λ_1 - 2 λ_2 + ϖ_2 (where λ_i ≡_i + ϖ_i is the mean longitude of i-th planet) equal 0 and π, respectively, while for P_2/P_1 < 2, ϕ_1 = π and ϕ_2 = 0. The Kepler-25 system has P_2/P_1 = 2.039, thus we will consider only the first case. Naturally, saying that the equilibrium of the averaged system corresponds to Δϖ = π and ϕ_1 = 0 means that for the unaveraged system Δϖ and ϕ_1 oscillate around those values. The amplitudes of the oscillations depend on the distance of P_2/P_1 to the resonant values of particular resonance <cit.>. The closer P_2/P_1 is to the nominal value of given MMR, the smaller are the amplitudes.

When searching for a periodic configuration one needs to set initial values of the mean anomalies, such that Δϖ = π and ϕ_1 = 0. There are two combinations of (_1, _2) that satisfy the condition, i.e., (0, 0) and (0, π). The results presented in Fig. <ref> were obtained for the former choice of the mean anomalies. When the angles _1, _2, Δϖ are set, we search for p⃗(t=0) that gives δ=0 in the parameter space of (P_2/P_1, e_1, e_2). As it is known <cit.> there exists a curve satisfying δ=0. In the regime of (P_2/P_1, e_1, e_2), that we are interested in, the curve can be parametrized with P_2/P_1, i.e., e_1 = e_1 (P_2/P_1) and e_2 = e_2 (P_2/P_1). Therefore, in order to find a family of periodic configurations one needs to search for δ(e_1, e_2) = 0 for a series of values of P_2/P_1 in a given range. The equation δ(e_1, e_2) = 0 is being solved numerically with a help of the Powell's method combined with the golden section <cit.>.

Figure <ref>a presents the families of periodic orbits obtained for P_2/P_1 ∈ (2, 2.09] and for various m_1, m_2 with their ratio kept constant, m_1/m_2 = 1. Each curve, presented in the (P_2/P_1, e_1)-plane corresponds to different value of m_1 + m_2. When the masses are lower, e_1 (and similarly e_2) are lower for given P_2/P_1. The period ratio value of Kepler-25, P_2/P_1 = 2.039, is marked with large gray dots. 
Panel (b) of Fig. <ref> presents the branches of periodic configurations in the eccentricities plane. Here, the sum of the masses is fixed, while the ratio varies. Clearly, when m_1/m_2 increases (for a given P_2/P_1) e_1 decreases, while e_2 increases.

As A_i ∝ e_i, by repeating the analysis for different m_1, m_2 one can find an approximate relation between the amplitudes A_1, A_2 and the masses, i.e., A_i = A_i (m_1, m_2; P_2/P_1). Roughly speaking, A_1 ∝ m_2 and A_2 ∝ m_1 in a regime of small masses. Therefore, it should be possible, in principle, to constrain the planets' masses by fitting the model of a periodic configuration to the TTV data (naturally, if the real system is close-to-periodic, what we demonstrate later in this work). The fitting procedure is presented in the next section.

In contrast with the amplitudes dependence on the planets' masses, a period of the rotation of the system as a whole (that equals to the (O-C)-signal period, ) does not depend on the masses significantly (in a small mass regime). Figure <ref>c presents  as a function of P_2/P_1 for periodic configurations found for the masses in a range of [5,30]. We note that the  dependence on the period ratio agrees with the analytic model of <cit.>, i.e., the super-period, that equals |1/P_1 - 2/P_2|^-1 for the 2:1 MMR (the analytic prediction is marked with a dashed curve in Fig. <ref>c). The Kepler-25 system, whose position at the plot is marked with a circle symbol, lies exactly at the curve corresponding to the families of periodic configurations. Nevertheless, the agreement between  of Kepler-25 and the super-period does not mean that the system has to be a periodic configuration, as the super-period corresponds to the variation of the longitude of conjunction <cit.>.




§ MODELLING TTV WITH PERIODIC CONFIGURATIONS


We use the TTV data from the catalogue of <cit.> and search for the best-fitting parameters of the two-planet model in terms of the minimum of the standard Χ function. What differs the procedure from the standard fitting approach is that here the orbital parameters (a_1, a_2, e_1, e_2, ϖ_1, ϖ_2, _1, _2) given at time t_0 are not free parameters. Each configuration for which Χ is being evaluated has to be periodic. Therefore, the optimization occurs at a hyper-surface embedded in the parameters space.

As we already showed in the previous section, the periodic configurations form a family parametrized by the period ratio [for given initial values of (Δϖ, _1, _2); we will use (π, 0, 0)] and the planets' masses. Particular values of the eccentricities as well as the angles (Δϖ, _1, _2) are functions of time, i.e., they depend on the phase of the periodic evolution, namely τ∈ [0, T). Although the dynamics of the system in scalable in a sense of physical dimensions as well as it is rotation invariant, when modelling the observations one needs to find an appropriate scale (that is given by initial P_1) and orientation. For the TTV analysis the only Euler angle that needs to be fitted is the angle that measures the rotation in the orbital plane (as we assume that the inclination I=π/2 and the TTV signal is invariant with respect to the rotation in the sky plane; formally we put Ω=0 for both longitudes of the ascending nodes). We denote the Euler angle by ϖ_0, and for a given configuration tested in the fitting procedure this value is being added to both longitudes of pericentre. 

Finally, the complete set of free parameters of the model x⃗ = (m_1, m_2, P_2/P_1, P_1, τ, ϖ_0). Direct minimizing Χ = Χ(x⃗) would be too long, as finding a periodic configuration for a given (m_1, m_2, P_2/P_1) occurs in a numerical process of solving δ(e_1, e_2) = 0. In order to make the procedure work efficiently, we treat the first two parameters (m_1, m_2) as fixed in a given fitting process, while the fitting is being repeated for subsequent points at the (m_1, m_2)-plane, taken from a grid, namely m_i ∈ [5, 20], with an increment of 0.75.

The period ratio (that has an osculating sense) is being chosen before a given TTV fitting process in such a way that the ratio of the periods of the transit times series {n, t_n} agrees with the observational value. Speaking in more details, for a given initial set of osculating Keplerian elements, the equations of motion are solved numerically, the transit times for both the planets are found and the linear model of TTs, i.e., T_n = T_0 + n ⟨ P ⟩ is fitted to the {n, t_n} series. Next, the value of ⟨ P_2 ⟩ / ⟨ P_1 ⟩ is compared to the observational value of the Kepler-25 system, i.e., ≈ 2.039. The scale of the system is being found in the same way, i.e., ⟨ P_1 ⟩ needs to equal the observational value of ≈ 6.2385d.




Therefore, for given (m_1, m_2) the Χ is being minimized in the two-paremeter space, i.e., Χ = Χ(τ, ϖ_0). We use the same numerical optimization scheme as in searching for the periodic orbits, described in the previous section.  Figure <ref> presents a scan of Χ computed at a grid of masses, as described above. 

There is a clear minimum of Χ around m_1 = 10.8 and m_2 = 14.5. Those values agree with the results in <cit.>, where they used an analytic model of a near-resonant system and obtained m_1 = (9.0 ± 2.6) and m_2 = (14.3 ± 4.5). Formal confidence levels plotted with white curves indicate that the masses are constrained within a ∼ 1 uncertainty. The best-fitting model, that is presented together with the observations in Fig. <ref>, reconstruct the data satisfactorily well. Nevertheless, it is not perfect, as Χ≈ 1.59 is greater than 1 (what could be, in principle, explained by underestimated measurement uncertainties). The parameters of the best-fitting system are listed in Table <ref>. Formal 1-σ uncertainties are very small. One should keep in mind, though, that the orbital parameters are not free parameters of the model in the common sense. The results listed in the table has to be understood as the most likely values under the assumption that the configuration is periodic. The standard fitting procedure, in which the parameters are free, will be applied to the data further in this work.

The parameters uncertainties listed in Tab. <ref> were computed with the uncertainty of the stellar mass was accounted for. The mass of the star is computed from the surface gravity (logg) and the radius of Kepler-25 listed in the catalogue in <cit.>. We obtain m_0 = (1.19 ± 0.05). The stellar mass uncertainty enlarges uncertainties of the planets' masses (slightly) and of the semi-major axes (significantly). For the completeness of the errors estimates Tab. <ref> lists also the planet-to-star mass ratio as well as the Keplerian periods, computed from the semi-major axes with a help of the Kepler's third law.









§ CAN A NON-RESONANT SYSTEM MIMIC A PERIODIC CONFIGURATION?


As we already mentioned, a periodic configuration is characterized by the (O-C)-signals in anti-phase. The observed system fulfils the criterion. Lets assume a sinusoidal model of TTV, i.e., (O-C)_i = A_i sin[(2π/) t + Φ_i], where  is the period and Φ_i is the phase of the TTV signal of the i-th planet. By fitting the model to the data one obtains A_1 = (3.8 ± 0.4)min, A_2 = (1.6 ± 0.2)min, = (325 ± 5)d and ≡ |(Φ_1 - Φ_2) - π| = (5 ± 9)deg. The latter quantity measures the deviation of the signals from the anti-phase. [Naturally, (Φ_1 - Φ_2) - π is kept in the range of (-π, +π).] For the periodic configuration = 0.

One may ask what is the value of  for a system far from periodic. In general, we expect that  is a function of initial orbital elements. Figure <ref>a presents a scan of this quantity in a plane of the eccentricities for a fixed period ratio and four representative pairs of the resonant angles (ϕ_1, ϕ_2). Each quarter of the plane represents a different combination of the angles, i.e., (0,0), (0,π), (π,0)and (π,π) (lets enumerate those quarters by I, II, III and IV, respectively). That means that there are two quarters with Δϖ=0 and two with Δϖ=π. 




A given representative combination of (ϕ_1, ϕ_2) can be achieved by two different combinations of (_1, _2, ϖ_1, ϖ_2). The choices for particular quarters are dictated by the continuity requirement at the x and y axes of the plane, where x = e_1 cosϕ_1 and y = e_2 cosϕ_2. We chose for quarter I (_1, _2, ϖ_1, ϖ_2) = (0, 0, π, π), for quarter II (0, π, π, 0), for quarter III (π, 0, 0, π) and for quarter IV (π, π, 0, 0). The period ratio equals 2.039 and the planets' masses equal the best-fitting values of m_1 = 10.8 and m_2 = 14.5. The white symbol points the position of the periodic configuration. The values of Δ_a-ph(x,y) (that are coded with a shade of gray) are computed by integrating the N-body equations of motion for a particular initial system, computing the series of TTs (in a window of 4.3yr, as it is for the Kepler-25 system) and computing the phases of TTV signals for both planets. 

There is clearly a line corresponding to Δ_a-ph=0 (given by y = 2.78722 x - 0.00412482) going through the periodic system, what means that other configurations with different (x,y) also produce TTV signals in anti-phase. By looking at this kind of plot alone, one cannot distinguish between the periodic and non-periodic (non-resonant) systems lying at the black line in Fig. <ref>a. Dashed lines indicate a value of 14 degrees, what is the maximal Δ_a-ph for the observed system.

Remaining panels of Fig. <ref> show scans of other quantities. Panel (b) presents the results for . The period of the TTV signal does not change significantly over the plane. For the whole plane  agrees with the value of the Kepler-25 system, i.e., (325 ± 5)days. Panels (c) and (d) illustrate the results for the semi-amplitudes of the TTVs. The line for which Δ_a-ph=0 corresponds to minima of A_1 and A_2. At this line both A_1 and A_2 equal the values of the best-fitting model. In panel (d) an additional information is given. Dashed areas indicate regions of oscillations of Δϖ and ϕ_1. The latter encompasses a small area around the periodic configuration. The region in which the second resonant angle ϕ_2 librates is smaller than the size of the white dot pointing the position of the periodic system. The areas of Δϖ-oscillations are wide and exist for both oscillations centres 0 and π. They correspond to non-resonant dynamics and represent the two modes of secular oscillations <cit.>. That means that a non-resonant system can mimic the periodic configuration.





 §.§ Probability that a non-resonant system mimics a periodic configuration


A periodic configuration does not evolve in long time-scales (apart from a uniform rotation of the system as a whole), therefore, regardless the epoch in which we observe the system the (O-C)-diagram looks the same. On contrary, a non-periodic configuration evolves in the secular time-scale. In order to illustrate that we chose an initial system that lies in a vicinity of the (Δ_a-ph=0)-line, i.e., e_1 = 0.006, e_2 = 0.013 and all the angles equal 0 (the point lies in quarter I). The evolution of this example configuration is shown in Fig. <ref>. The eccentricities as well as Δϖ vary in ∼ 680yr secular time-scale. The resonant angle ϕ_1 rotates much faster. 



The evolution of the system results in the variation of , what is illustrated in
Fig. <ref>a (note the wider time-window with respect to Fig. <ref>). We observe the variation in ∼ 680yr period (reaching ∼ 0) as well as a longer-period modulation (∼ 5kyr), that corresponds to the rotation of the system as a whole (see Fig. <ref>b for the evolution of individual values of ϖ_1, ϖ_2). We observed that the long-period modulation of  has minimal amplitude for ϖ_1 ≈ϖ_2 ≈ -π/2. Therefore, a system that is not periodic may look like one (in a sense of ∼ 0) for particular orientations of the orbits as well as for particular phases of the secular evolution.

A shaded area in Fig. <ref>a denotes the 14degree limit. There are epochs in the evolution of this system in which the TTV signals can be even in phase, i.e., Δ_a-ph = π. As a result, only for some part of the time the system, that is not periodic, looks like a periodic configuration. One can compute the probability that Δ_a-ph<14deg (denoted with p_1) by dividing the amount of time in which Δ_a-ph<14deg by the whole time of the integration. 
Because  depends on the two characteristics, i.e., on the phase in the secular modulation of the eccentricities and the spatial orientation of the system, instead of integrating a given initial configuration for a very long time, we integrate the system for the time that equals the secular period of ∼ 680yr and rotate the configuration within the whole range of 360degrees. The probability is computed by evaluating  every 20yr of the integration and by rotating the system at each epoch with an increment of 20deg. That makes 34 × 18 = 612 (O-C) diagrams for each initial configuration to be tested. A number of diagrams for which <14degrees divided by 612 defines p_1.  
For the system considered p_1 ≈ 40 per cent. In general, the probability p_1 is a function of (x,y).




The (x,y)-scan of the probability p_1 is presented in Fig. <ref>a. Naturally, for the periodic system and in its vicinity p_1 = 100per cent. The region of relatively high p forms an X-shape structure at the plane. That stems from the fact that a given point at the plane has three counterparts, as a given configuration, in general, intersects the representative plane in four points during the evolution. The positions of the intersections for the example configuration with e_1 = 0.006 and e_2 = 0.013 (ϕ_1 = ϕ_2 = 0) are marked with black circles. 

One can see that for higher eccentricities, especially in quarters II and III, p_1 is relatively high. We computed the scan in a wider range of e < 0.1, and observed that p_1 can reach even 100per cent. Panels (b), (c) and (d) of Fig. <ref> show that not only  depends on (x,y), also  and the amplitudes A_1, A_2 are functions of the eccentricities and the resonant angles. Moreover, those quantities depend on the phase in the secular evolution as well as the orientation of the system. Therefore, a given system with high p_1 can be, in general, characterized by , A_1, A_2 very different from the values of the observed system. Such a system cannot be considered as consistent with the observations. 

Another characteristic of the (O-C)-signal that needs to be considered is whether or not there is a second mode in it. Figure <ref> illustrates the Lomb-Scargle periodograms of the (O-C)-signals of the Kepler-25 system. There is no secondary peak higher than ∼ 0.27 and ∼ 0.33 of the primary peaks heights (for planets b and c, respectively).
Therefore, we require that for the synthetic systems considered as consistent with the observations the secondary-to-primary peak ratio is below the limits given above.

Figure <ref> presents the (O-C)-signals as well as the Lomb-Scargle periodograms for two configurations that have ∼ 0 (see the caption of this figure for the parameters). The system with e_1 = 0.006 has a uni-modal (O-C)-diagram, as it should be for a close-to-periodic configuration. The system with higher eccentricities is characterized with bi-modal (O-C)-signals. Therefore, this kind of configuration is not consistent with the observations.







In order to incorporate those criteria into the probability that a given configuration from the (x, y)-plane can be the real configuration of Kepler-25, we defined another quantity, p_2, which is a probability that a given system has < 14deg and also that A_1 ∈ (2.6, 5.0)min, A_2 = (1.0, 2.2)min, = (320, 330)d and that the second peaks in the periodograms of both signals are smaller than the values given above. The ranges of the allowed amplitudes correspond to the 3-σ confidence levels for the values of the observed signal. We widened the range because the amplitudes depend on the planets' masses, that are fixed for the plane at the values constrained with the periodic configuration assumption, that is not fulfilled for the whole plane.

Figure <ref>(b) shows the scan of p_2 in the same manner as in panel (a). Clearly, the only configurations that might correspond to the observed system are confined to the X-shape structure around the periodic configuration. Still, the probability that a non-resonant low-eccentric configuration lying in a vicinity of the line given by y = 2.78722 x - 0.00412482 can mimic the periodic configuration is relatively high. We computed a 1-dimensional scan of p_2 along the ( = 0)-line (not shown) and found that for e_1 ≳ 0.01 (and e_2 ≳ 0.027) the probability is already very low, p_2 ≲ 5per cent, while for e_1 ≳ 0.02 (and e_2 ≳ 0.054) p_2 = 0.



§ THE MCMC ANALYSIS


Before going to the MCMC (Markov Chain Monte Carlo) analysis of the TTVs we consider an influence of the third non-transiting planet in the system on the TTVs of the two transiting planets. <cit.> measured the radial velocities (RV) of the Kepler-25 and constrained the masses of planets b and c to be m_1 = (9.6 ± 4.2) and m_2 = (24.6 ± 5.7), that is roughly consistent with the model of the periodic configuration, although the mass of planet c is slightly smaller in our work. They also found that apart from the two planets discovered by the  mission <cit.>, there is a weak RV signal of the third, more distant companion.
The amplitude of the RV variation is relatively small and the precision of the measurements was not good enough to constrain the parameters of the third planet satisfactorily. There are two possible periods of planet d reported, i.e., (123 ± 2) and (93 ± 2) days. The mass of the planet m_3 = (89.9 ± 13.7) and the eccentricity e_3 = 0.18 (the uncertainty of e_3 was not given). There is a possibility, though, that the third planet affects noticeably the TTVs of the two transiting planets b and c. We tried to verify the possibility by performing a following test.

An expected result of an existence of the third planet in relatively wide orbit is an additional signal in (O-C) of the period equal to the orbital period of this planet and its harmonics <cit.>. The perturbing planet should not change the phase nor the period of the (O-C) of the resonant pair. In order to verify that and to check the amplitudes of the additional (O-C)-modulation due to the outermost companion, we 
added the third planet into the model, constructed the (O-C)-diagrams for planets b and c and compared the diagrams with the ones obtained for the unperturbed two-planet system. At first we chose the most likely values of m_3 and e_3. For planet b the difference is practically none, while for planet c the additional signal amplitude ≲ 0.1min, depending on the assumed period of the third planet (123 or 93 days) and a given angles (_3, ϖ_3) = (0, 0), (0, π), (π, 0) or (π, π). The A_2 variation is well below the TTVs uncertainties, moreover neither the period nor the phase of the (O-C)-signal are altered due to the third planet.
Next, we increased the mass of the outermost companion to 89.9 + 13.7 = 103.6 and tried to find at what value of e_3 (when increasing above the most likely value of 0.18) the amplitude of the additional (O-C)-modulation equals ∼ 0.5min, that is the noise level for (O-C) of planet c (see the Lomb-Scargle periodogram illustrated in Fig. <ref>). We found that for P_3 = 123days, the limiting e_3 = 0.4, while for P_3 = 93days, the limiting e_3 = 0.3. For eccentricities higher than the limiting values the third planet would produce an additional signal in (O-C)-diagram of planet c that is detectable. As we do not observe any additional periodicities in the (O-C)-diagrams, we conclude that the eccentricity e_3 is below the limiting values we found and the long-period companion does not need to be incorporated into the model.




 §.§ The standard TTV modelling and the periodic configuration


At present, the Bayesian inference is a de facto standard  for the analysis of the  light-curves and the TTV measurements. A crucial step in this approach is to define correct priors for determining the posterior distribution of model parameters. This is particularly important for interpretation of TTV models which
tend to degenerate solutions characterised by strongly aligned orbits with moderate and large eccentricities <cit.>. Therefore, the TTV fitting must be monitored whether or not its results depend on the adopted eccentricity priors, to avoid drawing incorrect conclusions on the orbital
archituecture of the studied planetary system. 

<cit.> assumed the Rayleigh distribution for the prior of the eccentricities and chose two different values of the Rayleigh parameter σ_e = 0.1 (wide prior, i.e., weak constraints on the eccentricities) and σ_e = 0.02 (narrow prior, i.e., stronger constraints). The posterior distributions of Δϖ for a few systems obtained for the two different values of σ_e differ significantly one from another. For the wide prior the most likely are the aligned configurations, i.e., with Δϖ∼ 0. When the prior is narrower, the peaks of the posterior distributions of Δϖ move towards ±π, although for the value they used, 0.02, the maxima of the posterior distribution are shifted by only 20-50degrees with respect to 0, depending on the system. On may expect, that for even lower σ_e, the peaks could move to the ±π positions. 

We showed previously (Fig. <ref>) that the aligned non-resonant configurations can mimic the periodic configuration. The (O-C)-signals of such qualitatively different systems look the same, what may bring problems when the direct fitting approach is applied. The eccentricities of the periodic configuration for P_2/P_1 = 2.039 are very small, ∼ 0.0001 - 0.001, therefore, finding such a system without strong a priori constraints on the eccentricities may be difficult. We presume that non-periodic (aligned) configurations are favoured in the fitting procedure (with wide priors put on the eccentricities), regardless the real architecture of the system. 

In this section we aim to verify this hypothesis by running extensive MCMC sampling of the posterior distribution. We used the TTV model described, for instance, in our previous papers <cit.>. Similarly, we imposed the Gaussian priors with the mean equal to 0 and a prescribed standard deviation σ set on the Poincaré parameters, X ≡ e cosϖ and Y ≡ e sinϖ, rather than on the eccentricities. The Poincaré elements are the free parameters of the TTV model and encode the eccentricity and the longitude of pericentre, respectively. We sampled the posterior distribution with the emcee package by <cit.>, choosing up to 256,000 iterations made with 1024 “walkers” in a small hypercube in the parameter space. 

We did a series of the MCMC experiments for 
σ∈ [0.00033, 0.1].
The results, selected for representative values of σ=0.033,0.0033,0.001, are illustrated in Fig. <ref>. Each row presents the results for one of the values of σ (given in the middle-column panels). The left-hand column shows the posterior probability distribution at the plane of (ϖ_1, Δϖ), the middle column is for the (m_1, m_2)-plane, while the right-hand column – for (e_1, e_2)-plane.
These results confirm our predictions. Indeed, for high σ aligned configurations are preferred, while lower σ leads to anti-aligned orbits. Although for high σ the anti-aligned systems are also allowed, they are less likely than the aligned ones. Another general observation is an existence of a correlation between σ and the masses and eccentricities, i.e., higher σ means lower masses and higher eccentricities.

The top row of Fig. <ref> illustrates the results for σ = 0.033, which may be representative for high σ, as for even higher values the qualitative picture is the same. The posterior probability distribution is bi-modal, as both aligned and anti-aligned systems are possible. The aligned orbits are, as already mentioned, more likely. For higher σ the disproportion between the modes is even more significant (the corresponding plots are not shown). The bi-modality is seen also at the eccentricities plane (right-hand column), although it is not that clear as for the (ϖ_1, Δϖ)-distribution. The most likely masses are very low, i.e., m_1 ≈ 1, m_2 ≈ 4.5, which gives low densities for both the planets, ∼ 0.3g/cm^3.



The second row from the top of Fig. <ref> shows the posterior probability distribution obtained for σ = 0.0033. The anti-aligned mode becomes more significant than the aligned one. The masses increase and the eccentricities decrease when compared to the previous example. The masses increases by a factor of 6 and 2.5, for the inner and the outer planet, respectively. That leads to the densities, respectively, ∼ 1.8g/cm^3 and ∼ 0.68g/cm^3. Clearly, different prior information on the eccentricities leads to qualitatively different both dynamical structure of the system and possible physical composition of the planets. 



The bottom row corresponds to the lower value of σ = 0.001. For this prior there exists only the anti-aligned mode (understood as a maximum of the posterior probability distribution). Predicted masses are higher with respect to the previous case, and the eccentricities are lower. For lower σ = 0.00033 (not shown) the resulting eccentricities are even lower (e_1 ≈ e_2 ≈ 0.0003) and the masses are higher (m_1 ≈ m_2 ≈ 16). Due to very low eccentricities the apsidal lines are poorly constrained, although, similarly to the previous case, the anti-aligned system is the most likely. The masses obtained for σ = 0.033 and σ = 0.00033 differ by a factor of ∼ 4 and ≳ 10 for the inner and the outer planets, respectively. The differences in the eccentricities are even larger, i.e., up to two orders of magnitude.

While different σ gives different masses, eccentricities and apsidal lines, the best-fitting configurations in terms of the highest posterior probabilities have different values of the standard Χ, for different σ. From a set of values we chose, the lowest Χ was obtained for σ = 0.001, i.e., Χ = 1.532. The parameters of that configuration are listed in Tab. <ref>. The quality of the fit is slightly better when compared with the parameters obtained for the periodic configuration (see Tab. <ref>). Nevertheless, the masses are both in agreement between the models, so do the eccentricities. The longitudes of pericentre as well as the mean anomalies differ significantly (probably due to very small eccentricities), however, the mean longitudes are in perfect agreement, as they should be.

The mass and radius of the outer planet are very similar to the values of Uranus, only the radius is ∼ 10% larger. The inner planet's mass and radius suggest the composition based mainly on water <cit.>, although an existence of water at such small distance from the star could be problematic. We stress again that the masses, eccentricities as well as the relative orientation of the apsidal lines depend strongly on an a priori information on the eccentricities. Such different orbital and physical parameters bring very different boundary conditions for the planetary systems formation theories, both in the aspect of the orbital characteristics as well as the internal structure of the planets.




§ CONCLUSIONS


We showed that a periodic configuration of two planets, that is a natural outcome of the migration, can be a good model of the TTV of the Kepler-25 planetary system. The period ratio of the system P_2/P_1 = 2.039 is significantly shifted from the nominal value of the 2:1 mean motion resonance, what may suggest a non-resonant nature of the system. We demonstrated that an anti-aligned resonant system produces the same (O-C)-diagrams as an aligned non-resonant configuration, however, the latter needs to fulfil certain criteria, like the orientation of the apsidal lines, or the relation between the eccentricities, in order to mimic the TTV signal of the periodic configuration. Due to low eccentricities of the resonant system of this value of P_2/P_1 (e_1 ∼ 0.0015,  e_2 ∼ 0.0002), as well as a degeneracy of the model mentioned above, finding a resonant configuration would be difficult (even if it was the true configuration of this system), without an a priori information on the eccentricities. 

We studied the probability that a non-resonant configuration mimics a periodic system. Although, such a non-resonant configuration can explain the TTV produced by the periodic system, the probability is lower when the configuration is further from the periodic system. We conclude that the real architecture of the Kepler-25 planetary system is very likely resonant in terms of librating resonant angles[We do not check the dynamical neighbourhood of the system, that could show whether or not the system lies in a region separated by a separatrix from the rest, non-resonant part of the phase space.], not only because of the probability test, but also because such a configuration is a natural outcome of the disc-induced migration, both convergent and divergent, what is believed to act a crucial role in the formation of the planetary systems.

We illustrated the dependence of the final orbital structure of the system as well as planets' masses on the assumed a priori information on the eccentricities. Wide prior probability distributions lead to aligned orbits and small masses, while narrower priors result in anti-aligned orbits and larger masses. 
The mass-eccentricity anti-correlation shown here in the series of MCMC experiments (Fig. <ref>) results from the degeneracy of the TTV signals that was illustrated with a help of an analytic model of a near-resonant two-planet system in <cit.>. On the other hand, the dependence of Δϖ (aligned/anti-aligned orbits) on the priors <cit.> can be understood from Fig. <ref>, i.e., the systems that lie along the line with = 0 have the same TTV amplitudes and the period. Both the degeneracies make the standard MCMC modelling a challenging task and additional knowledge on the formation mechanisms (and their expected outcomes) may be very useful when setting the eccentricity priors, that are crucial for determining the systems' parameters.

We argue that aligned configurations with relatively high eccentricities that seem to be common among analysed  systems may be artefacts and we believe that other systems with clear TTV signals should be verified in terms of their closeness to the periodic configurations, as true nature of the systems is essential for our understanding of the planetary systems formation, both their orbital configurations as well as physical compositions of the planets.



§ ACKNOWLEDGEMENTS


We would like to thank the anonymous reviewer for comments and suggestions that helped to improve the manuscript.
This work was partially supported by Polish National Science Centre MAESTRO grant DEC-2012/06/A/ST9/00276. 
K.G. thanks the staff of the Poznań Supercomputer and Network Centre (PCSS, Poland) for
the generous and continuous support, and for providing computing
resources (grant No. 313).


mn2e


