© ESO, 2017

1. Introduction

From the distribution of orbital period ratios of Kepler multiplanet systems, two features stand out. On the one hand, mean-motion resonances (MMRs) seem scarce and most of the systems appear nonresonant. This may either indicate that planets formed in situ (e.g., Hansen & Murray 2012; Petrovich et al. 2013; Chatterjee & Tan 2014) or that planetary migration did not necessary result in resonance trapping (e.g., Rein 2012; Baruteau & Papaloizou 2013). A second characteristic is the existence of a statistically significant number of planetary pairs close to resonances (particularly the 2:1 and 3:2 MMRs), but with orbital period ratios that are higher than the nominal value. Although such an offset is expected from resonant dynamics, particularly for low eccentricities, the observed values are much higher than obtained from simulations of resonance trapping. It is generally believed that these systems are near-resonant and located outside the libration domain.

Let us assume two planets of masses m₁ and m₂ orbiting a star m₀, orbital periods P₁<P₂, and in the vicinity of a first-order (p + 1) /p MMR. We can then define the so-called resonance offset as ${\mathrm{\Delta }}_{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathit{p}}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{2}}}{{\mathit{P}}_{\mathrm{1}}}\mathrm{-}\frac{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{p}}\mathit{,}$(1)whose value indicates the distance from exact resonance. Although the expected value depends on the planetary masses and eccentricities, N-body and hydrodynamical simulations of resonant capture in the 2:1 MMR generally lead to Δ_{2 / 1} ≃ 10^-3 (e.g., Silburt & Rein 2015) and similar values for the 3:2 resonance. However, a significant proportion of Kepler systems exhibit much higher values (i.e., Δ_{2 / 1} ~ 0.04 and Δ_{3 / 2} ~ 0.02, see Fig. 1), apparently placing them outside the resonance domain.

The origin of these near-resonant systems is a dilemma. Divergent migration through tidal evolution (e.g., Lithwick & Wu 2012; Batygin & Morbidelli 2013; Delisle & Laskar 2014) can reproduce some of the dynamical properties of these systems, although the magnitude of the tidal parameters required in many cases are not consistent with current models (Lee et al. 2013). Planetary migration in a turbulent disk (Paardekooper et al. 2013) can, under certain circumstances, lead to near-resonant configurations, although it is not clear whether this mechanism can in fact reproduce the observed near-resonant distribution (e.g., Quillen et al. 2013).

More recently, Migaszewski (2015, 2016) studied the effect of more complex disk models on migration of planetary pairs, including photo-evaporation and different opacity laws. They found that several regions of divergent migration appear in the disk that may lead to departure from exact resonance and even ejection from the commensurabilities. However, they found little correlation between their model and the observed distribution of bodies.

Fig. 1 Distribution of orbital period ratios of adjacent pair of planets as a function of the orbital period of the inner planet in the vicinity of the 2:1 MMR (top) and 3:2 MMR (bottom). Data obtained from exoplanet.eu. The red circles identify confirmed planets detected by transits or TTV, while the blue circles correspond to recently validated Kepler planets proposed by Morton et al. (2016). Planets discovered by any other method are depicted with open circles. In the top panel, open green circles denote the GJ876 system, while the dashed black circle groups the resonant (or near-resonant) members of HD 73526, HD 82943, HD 155358, and HD 92788.

The origin of asymmetric distributions of Kepler planets has been studied by Xie (2014) both analytically and numerically, suggesting that the offset could be explained regardless of dissipation. However, they conclude that dissipation or other mechanisms may play a more important role in the 2:1 than in the 3:2 MMR.

Perhaps a more basic question is how do we know these systems are actually nonresonant. Calculations of the libration width require information of the planetary masses, a luxury that is usually unavailable in Kepler systems. Estimations of the eccentricities are also important since the libration widths are a strong function of this element. However, the observed resonance offset is believed to be larger than the libration domain, especially for super-Earths and Neptune-size bodies.

In this work we present a simple analytical and N-body study of the dynamics of two-planet systems leading to capture in the 2:1 and 3:2 MMRs in a laminar disk. We show that the observed distribution of resonance offsets amongst the exoplanetary population is consistent with disk-planet interactions in a flared disk. Evidence from our simulations seems to indicate that the main difference between the warm and hot planetary systems is their distance from the central star, and many of the Kepler near-resonant planets may actually be inside the libration domain of the 2:1 and 3:2 commensurabilities.

The paper is structured as follows. In Sect. 2 we review the distribution of all confirmed planets in the vicinity of the 2:1 and 3:2 MMR and analyze the observed offset as function of the orbital period. Section 3 discusses how the resonance offset depends on the planetary and disk parameters and how high values may be attained without extreme eccentricity damping. We also analyze two different analytical prescriptions for migration. Section 4 presents a series of applications of our analytical model for flat and flared disks and different mass ratios of the planets. In particular, we search for values of the disk and planetary parameters that lead to offsets similar to those observed in real systems. Applications to several concrete cases are discussed in Sect. 5, while conclusions close the paper in Sect. 6.

2. Population near the 2:1 and 3:2 MMRs

Figure 1 shows the distribution of orbital period ratios P_{(i + 1)}/P_i for adjacent planets in multiplanetary systems, as function of the orbital period P_i of the inner member of each pair. The top plot shows the vicinity of the 2:1 MMR while the bottom panel presents results for the 3:2 commensurability. In both cases, red circles correspond to planetary systems detected by transits or transit time variations (TTV) and are mostly fruits of the Kepler mission. The recently validated systems proposed by Morton et al. (2016) are indicated separately in blue, since they are not actually confirmed planets and still subject to uncertainties. Finally, bodies discovered by all other methods, particularly with radial velocity, are identified by open black circles.

Both distributions show distinct features. Among the “non-transit” population, the 2:1 MMR hosts several well-known examples of resonant pairs and even multiple resonances. The HR 8799 planets are located very far from the central star (orbital periods in excess of 10⁴ days) and appear to be locked in a Laplace-type resonance. The orbital period ratio of the outermost planets are 2.0009 and 1.9995, respectively, placing them extremely close to the location of exact resonance (e.g., Goździewski & Migaszewski 2014.)

Closer to the star, and encompassed by an open dotted circle, lies a group of planetary systems very close to the 2:1 resonance and detected through radial velocity surveys. Perhaps the best known examples are HD 82943, whose two confirmed planets have a period ratio of ~2.017 (Baluev & Beaugé 2014), and HD 73526 with a value of 2.006 (Wittenmyer et al. 2014). Uncertainties in the orbital fits make it difficult to establish whether all other cases are truly resonant (e.g., HD 128311; see McArthur et al. 2014; Rein 2015), although stability considerations usually require libration of the resonant angles.

The last example of RV systems in the vicinity of the 2:1 MMR is GJ876 with three planets locked in a Laplace resonance (e.g., Rivera et al. 2010; Baluev 2011). Both adjacent pairs are identified in the upper panel by broad open green circles. Although the system also houses a much smaller inner planet, it is located very close to the star and with negligible gravitational interactions with the other members. The resonance offset for GJ876 is larger than observed in other cases with orbital period ratios of ~2.02 and ~2.04 (Nelson et al. 2016).

The GJ876 system has been shown to be significantly chaotic (Martí et al. 2013), although stable for time spans comparable to the age of the system. Batygin et al. (2015) argued that the strong chaoticity is incompatible with a smooth migration in a laminar disk, proposing instead that resonance capture occurred in a turbulent disk characterized by low-amplitude stochastic perturbations. However, it is unclear why turbulence would have affected this system and not the other members of the resonant population. On the other hand, Martí et al. (2016) showed that the Laplace resonance for the GJ876 also contains an inner region of very low chaoticity. Even though the current orbital fits place the system in the more stochastic outer resonant domain, it is possible that future observations may change this picture and be more consistent with the regular inner zone.

The distribution of transit and TTV systems around the 2:1 MMR (Fig. 1, filled red circles) shows two main features that have been the focus of several works over recent years. First, there is a significant excess of planets with positive values of Δ_{2 / 1} with a median close to 0.04 for P₂/P₁ ∈ [ 2.0,2.1 ]. Since the planetary radii of these planets are consistent with masses in the range of super-Earths and Neptunes, it is generally believed that these systems are outside the librational domain and nonresonant in nature. It is indeed curious to see that more distant systems, with typically higher masses, have smaller offsets even though their librational domain is significantly larger.

A second feature in the distribution of transit and TTV systems is a significant increase in the value of Δ_{2 / 1} for planets closer to the star. Delisle & Laskar (2014) also found a larger offset closer to the central star by analyzing two subsets: the first with P₁ ≤ 15 days and a second subset with P₁ ≥ 15 days. These authors attributed this dichotomy to tidal effects, which should be important for orbital periods below (say) P₁ ≤ 10−20 days. Figure 1 shows similar results, but goes one step further. The observed offset may indicate a possible smooth trend in which Δ_{2 / 1} ~ 0.05 for P₁ ~ 1 day, down to Δ_{2 / 1} ~ 0.01 for orbital periods P₁ ~ 10³ days. Such a trend could suggest that the distribution of resonance offsets is primarily dependent on the distance from the star and not so much on the migration type. The case of GJ876 appears to comply with this idea. Although detected through radial velocity and consisting of giant planets that should have experienced Type-II migration, their value of Δ_{2 / 1} is similar to that of Kepler systems and not to other, more distant, resonant cases comprised of giant bodies.

The bottom panel of Fig. 1 shows results for the vicinity of the 3:2 MMR. Presently there is only one confirmed resonant RV system (i.e., HD 45364; Correia et al. 2009), with P₁ ~ 226 days and an offset of Δ_{3 / 2} ~ 0.01. Thus, at least in the case of radial velocity surveys, the resonant population of the 3:2 is much less abundant than observed in the 2:1.

The 3:2 transit and TTV systems also show significant differences with respect to the 2:1 commensurability. On the one hand, the resonant offset is noticeable smaller with a median of Δ_{3 / 2} ~ 0.01, as calculated for systems with P₂/P₁ ∈ [ 1.5,1.6 ]. This value is on the same order as that observed for HD 45364, indicating a shallower dependence with the orbital period. In fact, as pointed out in Wang & Ji (2014), the distribution of Kepler systems is consistent with an accumulation close to exact resonance with small offsets, and no significant gap is observed for low negative values of Δ_{3 / 2}. A second interesting feature is an apparent pile-up of (near)-resonant systems with orbital periods P₁ ~ 10 days, which is also not observed in the case of the 2:1 resonance.

In this paper we focus our attention on Kepler and other close-in systems detected mainly through transits and TTV. Although in most cases there is little information on the planetary masses, it is believed they are relatively small bodies (Earth up to Neptune mass range) and would have experienced planetary migration without having carved a significant gap in the disk. We then use analytical models for Type-I migration in an attempt to understand the observed distribution and general trend.

3. Resonance offset

3.1. Definition and analytical modeling

In principle, there are two ways to explain the (near)-resonant structure of the exoplanetary systems and its dependence on the distance to the central star. In the first explanation, the planets were initially trapped inside the resonance (exhibiting low values of their offset) but later were driven out of the libration domain through tidal effects. This appears reasonable for some systems in the 2:1 resonance, although it does not explain why the offset observed in the 3:2 commensurability is much smaller. Moreover, tidal effects do not seem to explain the trend observed in the 2:1 MMR for more distant planets (e.g., Lee et al. 2013).

A second possibility is that resonance trapping, due to disk-planet interactions, led to different values of Δ_{(p + 1) /p} in different parts of the disk. Thus, a planetary migration that was halted far from the central star would be characterized by small offsets, while this value increased from smaller orbital distances. This is the mechanism that we analyze in the present work. However, before simulating the migration process, it is important to review the dynamical structure of planetary resonances and analyze the types of motion expected for different values of Δ_{(p + 1) /p}.

Fig. 2 Dynamical maps of max(Δe) for the 2:1 (left) and 3:2 MMR (right) for fictitious two-planet systems with m1 = 0.05mJup and m2 = 0.10mJup orbiting a central star with mass m0 = 1m⊙. The orbit of the outer planet was initially circular with a2 = 1 AU, and all angular variables where chosen equal to zero. Total integration time was 103 yr. The green lines indicate the location of the zero-amplitude ACR-type librational solutions, estimated from the simple analytical model (see Eq. (14)). The white dots in the left panel are the result of three N-body simulations of resonance trapping in the 2:1 resonance. The broad black curve in the right plot corresponds to the Marchal-Bozis stability limit. See text for details.

Figure 2 shows two dynamical maps: the left panel for the 2:1 resonance and the right panel for the 3:2 commensurability. In each case we integrated a series of two-planet systems with initial conditions in a grid defined in the (P₂/P₁,e₁) plane. The masses of the planets were chosen equal to m₁ = 0.05m_Jup and m₂ = 0.10m_Jup, and the central star assumed of solar mass. We specifically chose a mass ratio m₂/m₁> 1 to guarantee symmetric fixed points for the resonant angles (e.g., Beaugé et al. 2006; Michtchenko et al. 2008). The orbit of the outer planet was initially considered circular with a₂ = 1 AU and all angles equal to zero.

The color code corresponds to the maximum values of | e₁(t)−e₁(t = 0) | (denoted here as max(Δe)) attained during a 10³ year integration time span. Darkened (lighter) tones are associated with small (large) variations in the eccentricity of the inner planet. For the 3:2 resonance, the broad black curve indicates the stability limit as calculated using the Marchal-Bozis criterion (Marchal & Bozis 1982). Although the max(Δe) indicator does not measure chaotic motion, it is an important tool with which to probe the structure of resonances and help identify the locus of stationary solutions, known as the ACR solutions (see Beaugé et al. 2003, 2006), and the separatrix delimiting the librational from the circulation domains (e.g., Ramos et al. 2015).

In the same figure, the green lines show the approximate location of the family of zero-amplitude ACR solutions (Beaugé et al. 2003, 2006), which are characterized by the simultaneous libration of both resonant angles $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{-}\mathit{p}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{-}\mathit{\varpi }1\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{\theta }}_{\mathrm{2}}\mathrm{=}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{-}\mathit{p}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{2}}\mathit{,}\end{array}& & \end{array}$(2)where λ_i and ϖ_i are the mean longitudes and longitudes of pericenter of each planet. Using a simple Sessin-type approximation for the resonant Hamiltonian (e.g., Sessin & Ferraz-Mello 1984), which neglects the secular perturbations and is valid only for very low eccentricities, we can estimate the value of the resonance offset as function of e₁, along with a relationship between the eccentricities of both planets. These give ${\mathrm{\Delta }}_{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathit{p}}\mathrm{=}{\mathit{C}}_{\mathrm{1}}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\frac{\mathrm{1}}{{\mathit{e}}_{\mathrm{1}}}\mathrm{;}{\mathit{e}}_{\mathrm{2}}\mathrm{=}{\mathit{C}}_{\mathrm{2}}\frac{{\mathit{m}}_{\mathrm{1}}}{{\mathit{m}}_{\mathrm{2}}}{\mathit{e}}_{\mathrm{1}}$(3)(see Lee 2004), where the coefficient C_i depends solely on the semimajor axes ratio a₁/a₂. For the MMRs under consideration, they take the following values: C₁ ≃ 1.5,C₂ ≃ 0.29, for the 2:1 resonance; and C₁ ≃ 1.2,C₂ ≃ 1 for the 3:2 case.

For given resonances and planetary masses, the value of the offset is inversely proportional to the equilibrium eccentricities; thus, very low values of e_i are necessary to obtain a significant deviation from an exact resonance. As an example, reproducing a resonant offset of Δ_{2 / 1} ≃ 0.05 would require the eccentricities to remain as low as e₁ ~ 10^-3−10^-4 during the final stages of planetary migration.

If the disk-driven planetary migration is sufficiently slow and smooth, and if the orbits remain almost circular, we expect the orbital evolution to follow the pericentric branch into the librational domain and exhibit low-amplitude oscillations of the resonant angles. In such an ideal scenario, the final eccentricities and resonant offset Δ_{(p + 1) /p} depends on the relative strength between the eccentricity damping and orbital migration timescales (Lee & Peale 2002; Beaugé et al. 2006), usually denoted by τ_ei and τ_ai, respectively. Thus, in the adiabatic limit, the final outcome of a resonance trapping depends fundamentally on the ratios , which we denote here as the K-factors.

Papaloizou & Szuszkiewicz (2005) deduced a closed relation between the equilibrium eccentricities and the orbital and damping timescales, which, after some simple algebraic manipulations, acquires the form ${\mathit{e}}_{\mathrm{2}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{𝒦}}_{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{D}\mathrm{)}\mathrm{+}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{𝒦}}_{\mathrm{1}}\mathit{B}\left(\frac{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{2}}}}{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{1}}}}\right)\mathrm{=}\frac{\mathit{D}}{\mathrm{2}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{2}}}}{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{1}}}}\right)\mathit{,}$(4)where the new coefficients $\mathit{D}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}}{\left(\mathrm{1}\mathrm{+}\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathrm{2}}}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}\right)}^{-1}\mathrm{;}\mathit{B}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{2}}{\mathit{a}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{2}}{\mathit{n}}_{\mathrm{1}}{\mathit{a}}_{\mathrm{1}}}\mathrm{+}\mathit{D}$(5)depend solely on the resonance and the mass ratio of the planets. Using Eqs. (3), (4), we finally obtain an expression for the resonance offset in terms of the migration timescales as ${\mathrm{\Delta }}_{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathit{p}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}}{\mathit{D}}{\left({\mathit{C}}_{\mathrm{1}}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\frac{\left[\mathrm{(}\mathrm{1}\mathrm{-}\mathit{D}\mathrm{)}{\mathrm{𝒦}}_{\mathrm{2}}{\mathrm{(}{\mathit{C}}_{\mathrm{2}}{\frac{{\mathit{m}}_{\mathrm{1}}}{{\mathit{m}}_{\mathrm{2}}}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{+}\mathit{B}{\mathrm{𝒦}}_{\mathrm{1}}\left(\frac{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{2}}}}{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{1}}}}\right)\right]}{\mathrm{1}\mathrm{-}\left(\frac{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{2}}}}{{\mathit{\tau }}_{{\mathit{a}}_{\mathrm{1}}}}\right)}\mathrm{·}$(6)A similar expression has been already found by Xie (2014). Since resonance trapping only occurs in cases of convergent migration (i.e., τ_a1>τ_a2), the denominator is always positive and free from singularities.

3.2. Offset and the K-factors

Assuming the constant ratio m₁/m₂, expression (6) shows that the value of the offset is a linear function of m₂ indicating that, for fixed values of the migration timescales, we would expect larger offsets for more massive planetary pairs. Since this does not seem to be the case, then the observed values of Δ_{(p + 1) /p}, particularly for close-in systems, must be dominated by the characteristics of the disk-planet interactions.

To test our model defined by expressions (3)–(6), we can perform simple tests with an N-body integrator including an ad hoc external acceleration (e.g., Cresswell & Nelson 2008): ${\begin{array}{c}{¨}\\ {r}\end{array}}_{\mathit{i}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\tau }}_{{\mathit{a}}_{\mathit{i}}}}[\frac{{{v}}_{\mathit{i}}}{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathrm{𝒦}}_{\mathit{i}}\frac{\mathrm{(}{{v}}_{\mathit{i}}\mathrm{·}{{r}}_{\mathit{i}}\mathrm{)}{{r}}_{\mathit{i}}}{{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}}]\mathit{,}$(7)where r_i is the radial vector and v_i the velocity vector of the planet. Setting the values of τ_ai at certain prefixed values, we can then vary the values of and analyze its effects on the resonance offset.

We performed three N-body simulations of planetary migration using this simple recipe. We used the same masses and initial orbit of the outer planet as used to construct the dynamical maps. The initial value of a₁ was chosen outside the 2:1 MMR and all angles equal to zero. The orbital decay timescales were τ_a1 = 10⁵ and τ_a2 = 7 × 10⁴ yr to guarantee convergent and adiabatic migration. The values of the eccentricity damping timescales τ_ei were chosen to give equal K-factors for both planets (i.e., ). Each run considered a different value of the K-factor: , and .

The initial conditions were integrated for 10⁵ yr. We assumed that the disk began dispersing at T_disp = 7 × 10⁴ yr and its surface density reached zero at T_disp = 8 × 10⁴ yr. This process was modeled by reducing the magnitude of the dissipative force smoothly with a hyperbolic tangent function, and its timescale was sufficiently slow to allow the equilibrium solution of the system to adapt adiabatically. The last 10% of the outputs were used to calculate averaged values of the mean motions and eccentricities. The white symbols in the left panel of Fig. 2 show the final results.

The run with ended with e₁ ≃ 0.03 and P₂/P₁ ≃ 2.004, that is, very close to exact resonance and reminiscent of the individual pairs of HR 8799. The second run, with yielded a lower equilibrium eccentricity (e₁ ≃ 0.01) and a slightly larger offset (P₂/P₁ ≃ 2.01), similar to that of several of the resonant RV planets (e.g., HD 82943). Finally, the simulation with showed the lowest final eccentricity and largest offset; i.e., P₂/P₁ ≃ 2.03, which is on the order of that observed for GJ876. All these values are in perfect agreement with the estimations deduced from expressions (3) and (6).

K-factors on the order of 10⁴ are already much higher than predicted by linear models of Type-I disk-planet migration (e.g., Papaloizou & Larwood 2000; Tanaka & Ward 2004; Cresswell & Nelson 2008). From these works, in the low eccentricity limit, we may write ${\mathit{\tau }}_{{\mathit{a}}_{\mathit{i}}}\mathrm{=}{\mathit{Q}}_{\mathit{a}}\frac{{\mathit{t}}_{{\mathrm{wave}}_{\mathit{i}}}}{{\mathit{H}}_{{\mathit{r}}_{\mathit{i}}}^{\mathrm{2}}}\mathrm{;}{\mathit{\tau }}_{{\mathit{e}}_{\mathit{i}}}\mathrm{=}{\mathit{Q}}_{\mathit{e}}\frac{{\mathit{t}}_{{\mathrm{wave}}_{\mathit{i}}}}{\mathrm{0.780}}\mathit{,}$(8)where H_ri is the disk aspect-ratio in the position of each planet. The value Q_e is a constant and Q_a = Q_a(α) is a function of the slope α of the surface density profile. The linear model by Tanaka et al. (2002) found ${\mathit{Q}}_{\mathit{a}}^{-1}\mathrm{\simeq }\mathrm{2.7}\mathrm{+}\mathrm{1.1}\mathit{\alpha }$, while numerical fits with hydro-simulations usually lead to slightly different functional forms (e.g., D’Angelo & Lubow 2010). Finally, ${\mathit{t}}_{{\mathrm{wave}}_{\mathit{i}}}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{0}}}{{\mathit{m}}_{\mathit{i}}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathrm{\Sigma }\mathrm{\left(}{\mathit{a}}_{\mathit{i}}\mathrm{\right)}{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}}\frac{{\mathit{H}}_{{\mathit{r}}_{\mathit{i}}}^{\mathrm{4}}}{\mathrm{\Omega }\mathrm{\left(}{\mathit{a}}_{\mathit{i}}\mathrm{\right)}}\mathit{,}$(9)is the typical timescale for planetary migration. Here $\mathrm{\Omega }\mathrm{\left(}{\mathit{a}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{𝒢}\mathrm{(}{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathit{i}}\mathrm{)}\mathit{/}{\mathit{a}}_{\mathit{i}}^{\mathrm{3}}}$ is the (circular) Keplerian angular velocity and the gravitational constant.

We assume a laminar disk whose surface density Σ and the aspect-ratio H_r have a power-law dependence with the distance from the star, $\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\mathit{r}}^{\mathrm{-}\mathit{\alpha }}\mathrm{;}{\mathit{H}}_{\mathit{r}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{H}}_{\mathrm{0}}{\mathit{r}}^{\mathit{f}}\mathit{,}$(10)where Σ₀ and H₀ are the values at r = 1 AU. The exponent f is usually referred to as the disk-flare and is zero for flat disks where the aspect-ratio is constant.

Finally, the quantity Q_e in expression (8) is an artificial ad hoc constant parameter introduced by Cresswell & Nelson (2006) to fit the results of hydrodynamical simulations. They found that the best results were obtained with Q_e ~ 0.1, a value that we also adopt here.

From Eqs. (8)–(10) we can now obtain a simple explicit form for the -factor as ${\mathrm{𝒦}}_{\mathit{i}}\mathrm{\equiv }\frac{{\mathit{\tau }}_{{\mathit{a}}_{\mathit{i}}}}{{\mathit{\tau }}_{{\mathit{e}}_{\mathit{i}}}}\mathrm{=}\frac{\mathrm{0.78}{\mathit{Q}}_{\mathit{a}}}{{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{Q}}_{\mathit{e}}}{\mathit{r}}_{\mathit{i}}^{\mathrm{-}\mathrm{2}\mathit{f}}\mathit{.}$(11)Assuming flat disks (f = 0) and typical values of H₀ ~ 0.05, we obtain factors on the order of for any distance from the star, which are the values usually employed in N-body simulations of planetary migration and resonance capture (e.g., Quillen et al. 2013; Silburt & Rein 2015). A flared disk (f> 0) not only leads to a dependence of on the semimajor axis, but also increases its value for close-in planets. Adopting f = 0.25 we obtain similar values as before for a = 1 AU but significantly higher () for a = 0.05 AU.

Although this road seems promising, the increase in the K-factor is not sufficient. A resonance offset of Δ_{2 / 1} ≃ 0.05, such as those observed in close-in systems, would require , which at least seems improbable. Xiang-Gruess & Papaloizou (2015) proposed that planetary traps close to the central star could halt planetary migration while not affecting the eccentricity damping timescale, leading to an artificial increase of the corresponding magnitude of . Even if such an idea seems possible, it would not explain the observed trend of decreasing offset as function of the semimajor axis, nor would it be expected to work for orbital periods on the order of 10² days.

3.3. Goldreich and Schlichting prescription

Goldreich & Schlichting (2014) proposed to modify expressions (8) to account for the contribution of eccentricity damping to changes in the semimajor axis associated with (partial) conservation of the angular momentum. Thus, the effective characteristic timescale for orbital migration should actually be given by$\frac{\mathrm{1}}{{\mathit{\tau }}_{{{\mathit{a}}_{\mathrm{eff}}}_{\mathit{i}}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\tau }}_{{\mathit{a}}_{\mathit{i}}}}\mathrm{+}\mathrm{2}\mathit{\beta }\frac{{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}{{\mathit{\tau }}_{{\mathit{e}}_{\mathit{i}}}}\mathit{,}$(12)where τ_ai and τ_ei maintain the same form as Eqs. (8) and β is a factor that quantifies the fraction of the orbital angular momentum preserved during the migration. The value β = 1 if the orbital angular momentum is completely preserved (e.g., tidal effects for synchronous orbits), while β< 1 in other cases. Goldreich & Schlichting (2014) refer to estimations by Tanaka & Ward (2004) that suggest β ≃ 0.3 for Type-I disk-planet interactions.

The modified migration timescale, of course, also changes the value of the K-factor, leading to a new “effective” form, ${\mathrm{𝒦}}_{{\mathrm{eff}}_{\mathit{i}}}\mathrm{\equiv }\frac{{\mathit{\tau }}_{{\mathit{a}}_{{\mathrm{eff}}_{\mathit{i}}}}}{{\mathit{\tau }}_{{\mathit{e}}_{\mathit{i}}}}\mathrm{=}\frac{{\mathrm{𝒦}}_{\mathit{i}}}{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{\beta }{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}{\mathrm{𝒦}}_{\mathit{i}}}\mathit{,}$(13)where is given by Eq. (11). Since the equilibrium eccentricities within a resonance usually are on the order of , at least at the order of magnitude level, we expect that the new value of the K-factor should not be much different from its original. Moreover, since the denominator in (13) is always larger than unity, it should be verified that . Consequently, angular momentum considerations in the migration prescription do not change the picture significantly and cannot account for K-factors necessary to explain the observed resonance offsets.

The change in the effective migration timescale also affects Eq. (6) for the resonance offset. After some algebraic manipulations, we found a new version given by

$\begin{array}{ccc}& & {\mathrm{\Delta }}_{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathit{p}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}}{\mathit{D}}{\left({\mathit{C}}_{\mathrm{1}}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\\ & & \end{array}$(14)

where the expressions for τ_ai and are those defined in (8) and (11).

Fig. 3 Orbital period ratio (p + 1) /p + Δ(p + 1) /p as a function of the period of the inner planet, as predicted by expression (14) for flared disks with H0 = 0.05, α = 1, and f = 0.25. Top panels: represent the 2:1 MMR while bottom plots represent the 3:2 resonance. Each graph assumes a different value for the mass (in units of m⊕) of the inner planet, whose value is indicated in the top right-hand corner. Colored curves correspond to different mass ratios: m2/m1 = 1.6 (green), m2/m1 = 1.7 (blue), m2/m1 = 1.8 (purple), and m2/m1 = 2 (orange). Open black circles denote the distribution of confirmed planets.

4. Differential migration and resonance offset

Even if planetary migration occurs with limited values of the K-factors, a significant value for the resonance offset may still be obtained if the denominator in expression (14) is sufficiently small. As mentioned previously, this factor measures the differential migration timescale between both planets. In order for resonance capture to be possible, the outer planet must migrate faster than its inner companion (i.e., τ_a2<τ_a1) leading to positive values of (1−τ_a2/τ_a1) and, consequently, real values of Δ_{(p + 1) /p}.

From the expressions for τ_ai, the condition for convergent migration implies that the mass ratio m₂/m₁ between the planets must satisfy the condition $\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}\mathrm{\ge }{\left(\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}\right)}_{\min }\mathrm{\equiv }(\frac{{\mathit{P}}_{\mathrm{2}}}{{\mathit{P}}_{\mathrm{1}}}{)}^{\frac{\mathrm{2}}{\mathrm{3}}\mathrm{(}\mathrm{2}\mathit{f}\mathrm{+}\mathit{\alpha }\mathrm{-}\mathrm{0.5}\mathrm{)}}\mathrm{·}$(15)For given initial conditions, this equation defines the limiting mass ratio leading to capture in different commensurabilities. In the general case of flared disks, the possibility of trapping in a given resonance not only depends on the mass ratio itself but also on the initial orbital separation between the planets.

Complying with condition (15), the next question is how does the value of the offset vary as a function of the mass ratio and orbital distance. Particularly, we wish to analyze whether there exists any set of parameters that appear consistent with the observational distribution of resonant and near-resonant systems.

Figure 3 shows an application of expression (14) for H₀ = 0.05, α = 1 and f = 0.25, in the form of the predicted equilibrium values of the orbital period ratio P₂/P₁ = (p + 1) /p + Δ_{(p + 1) /p} as a function of P₁. Results for the 2:1 MMR are shown in the top panels, while those for the 3:2 are presented in the bottom plots. Observational data is shown in open black circles. Each panel corresponds to a different value of m₁, while each color curve corresponds to a different value of the mass ratio m₂/m₁.

In all cases the offset increases with smaller semimajor axes, although its magnitude is a strong function of both m₁ and m₂/m₁. As expected from Eq. (15), mass ratios close to (m₂/m₁)_min show high values of Δ_{(p + 1) /p}, although (also expected) its value increases linearly with m₁. For the 2:1 MMR, and for these values of the disk parameters, the observed distribution of resonance offsets amongst exoplanetary pairs does not seem compatible with inner masses much lower than m₁ ~ 10m_⊕. Although this limiting case is marginally adequate, a much larger diversity in values is obtained assuming m₁ = 20m_⊕, where practically all of the observed offsets are encompassed by the model, even if the mass ratios are not restricted to values very close to unity.

The bottom plots, corresponding to the 3:2 MMR, show a different story. In all cases the offsets increase with smaller distances to the star, but show very little dependence with both m₁ and the mass ratio. The best results are obtained assuming m₁ ~ 10m_⊕ (at least for these disk parameters). The lower dispersion in Δ_{3 / 2}, with respect to Δ_{2 / 1}, is also in good agreement with the observations.

These results were obtained from our analytical model leading to expression (14) and still require numerical confirmation. The upper left-hand plot of Fig. 4 show four N-body simulations, where planetary migration was modeled according to Eq. (7), incorporating the Goldreich-Schlichting prescription through expressions (12) and (13). We assumed β = 0.3, H₀ = 0.05, and a surface density of Σ₀ = 400 gr/cm². The factors Q_a and Q_e were chosen as described at the end of Sect. 3.2. In all cases we simulated capture into the 2:1 MMR with m₁ = 20m_⊕ and different values for the mass ratio m₂/m₁. As predicted, the offset effectively increases closer to the star, showing a very good agreement with the observations.

The other three plots show different aspects of the simulation with m₂/m₁ = 1.6, where the value for the mass ratio leads to the largest offset. The top right-hand graph shows the evolution of both critical angles as function of P₂/P₁, where we have neglected the data points prior to resonance lock. The main resonant angle θ₁ exhibits a low-amplitude libration around zero, even for offsets Δ_{2 / 1} ≃ 0.1, which are values that are usually associated with nonresonant motion. The auxiliary resonant angle θ₂ also librates throughout the evolution of the system, although with increasing amplitude. Thus, even in this extreme numerical example, the two-planet system remains locked in an ACR-type solution regardless of its proximity to (or separation from) exact resonance.

The two bottom panels of Fig. 4 show complementary information: the eccentricity of the inner planet as function of P₂/P₁ (left) and the relation between both planetary eccentricities (right). Thin black continuous lines indicate the expected functional form as given by expressions (3) obtained for ACR solutions with the Sessin resonant Hamiltonian. The agreement with the N-body run is excellent, indicating that the increase in resonance offset follows closely the loci of low-amplitude ACR solutions and the evolution of the system continues to be dominated by the commensurability even for high values of Δ_{2 / 1}.

We performed a large series of N-body simulations covering different values of the planetary masses and flare index, observing the same dynamical evolution as described in Fig. 4. Thus, contrary to expectations, large offsets may still be linked to resonant capture, even if the resulting ACR solution corresponds to a kinematic (and not dynamical) libration domain (e.g., Henrard & Lamaître 1983).

Fig. 4 Top left panel: result of four N-body simulations leading to capture in the 2:1 MMR. In all cases m1 = 20m⊕. Color lines indicate different mass ratios: m2/m1 = 2 (orange), m2/m1 = 1.8 (purple), m2/m1 = 1.7 (blue), and m2/m1 = 1.6 (green). Remaining panels: different dynamical characteristics of the orbital evolution within the resonance of the run with m2/m1 = 1.6. The black continuous lines in the bottom panels represent the analytical predictions given by expressions (3).

5. Application to individual systems

The mechanism proposed in this paper depends heavily on the individual masses and mass ratio of the planets located in the vicinity of the 2:1 and 3:2 MMR. The increase in the resonance offset close to the star is most prominent when the mass ratio is close to unity and m₁ at least on the order of 10m_⊕. More importantly, in order for convergent migration to take place, the outer planet must be more massive than its inner companion by a factor strongly dependent on the flare index assumed for the disk.

It is not easy, however, to evaluate whether observed exoplanets in (near)-resonant configurations validate this scenario. Most systems where Type-I migration is applicable have been detected by transit, making accurate measurements of the masses a rare commodity.

Fig. 5 Top: distribution of radii between confirmed planets in the vicinity of the 2:1 (blue) and 3:2 (red) MMR. Dashed line corresponds to equal-size bodies (R1 = R2). Middle: radius as a function of the orbital period. Light gray circles indicate the values of planets outside both resonances. Bottom: planetary mass mi vs. radius Ri for all planets (bodies larger than Jupiter shown in gray). Dashed orange line corresponds to the so-called R2.06-law by Lissauer et al. (2011), while the fit proposed by Weiss & Marcy (2014) is shown as a continuous orange line. Our own fit, which is calculated for all planets smaller than Jupiter, is presented as a thick green line. Values within 1σdis dispersion are contained within the light green region.

In order to have at least some statistical validation, we analyzed the known exoplanetary population of confirmed planets, as published in the exoplanet.eu web site in mid-December 2016. Results are shown in Fig. 5. In the top panel we plotted the distribution of radii (R₂ versus R₁) for all planetary pairs in the vicinity of the 2:1 and 3:2 resonances. We only considered those bodies for which the uncertainties in radius were smaller than 40% and resonance offsets in the interval were (0,0.1 ]. Most of the planetary pairs have radii below 4 R_⊕, although a few belong to the giant planet domain up to and beyond on Jupiter radius R_Jup. Although in most cases R₂/R₁> 1, a few are characterized by bigger inner planets. However, without estimations of the planetary density, it is not possible to say whether the masses also comply with this condition and would not be possible candidates for resonance capture through inward planetary migration.

The middle plot shows the distribution of radii as function of the orbital period. The color code is the same as above, while all bodies outside the 2:1 and 3:2 MMRs are indicated by light gray circles. While the population of both resonances show decreasing radii close to the central star, this is particularly evident for the 2:1 commensurability, where R_i ~ 2 R_⊕ for P_i< 10 days, in contrast to larger distances where larger planets abound. The decrease in R_i closer to the star, where the observed resonance offset is larger, appears contrary to expectations from our model (14), which indicates that Δ_{(p + q) /p} ∝ m_i/m₀. As we see later on, this problem can be partially overcome considering tidal evolution during the lifespan of the system.

Finally, since our model requires estimations of the masses, the bottom panel shows the mass-radius distribution of all known planets (resonant or otherwise) whose values have errors below 40%. Two different empirical fits are usually employed for this relation. The ${\mathit{R}}_{\mathit{i}}^{\mathrm{2.06}}$ power law from Lissauer et al. (2011) is shown as a broken orange curve and was deduced from the values of Earth and Saturn. As can be seen from the figure, it shows a fair fit for bodies with R_i> 5 R_⊕, but underestimates the mass of exoplanets with smaller radii. Conversely, the fit presented by Weiss & Marcy (2014) shows a better agreement for smaller bodies, but still falls short of the average values even in this domain. Many small planets have been confirmed since Weiss & Marcy (2014) and a similar analysis performed today would lead to a different set of numerical values.

To obtain an updated empirical relation, we then searched for two best fits, one for R<R_crit and one for larger bodies, where the value of R_crit was chosen close to 4 R_⊕ but with some slight freedom to reduce the overall variance. The result is shown as two thick green lines. Since we also need to model the dispersion around the mean values, the light green shaded areas corresponds to the regions within 1σ_m. Explicitly, we obtained $\log (\frac{{\mathit{m}}_{\mathit{i}}}{{\mathit{m}}_{\mathrm{\oplus }}})\mathrm{\simeq }\{\begin{array}{c}\\ \mathrm{0.60}\mathrm{+}\mathrm{0.92}\hspace{0.17em}\log (\frac{{\mathit{R}}_{\mathit{i}}}{{\mathit{R}}_{\mathrm{\oplus }}})\mathrm{if}\mathit{R}\mathrm{\le }\mathrm{3.8}{\mathit{R}}_{\mathrm{\oplus }}\\ \mathrm{-}\mathrm{0.13}\mathrm{+}\mathrm{2.18}\hspace{0.17em}\log (\frac{{\mathit{R}}_{\mathit{i}}}{{\mathit{R}}_{\mathrm{\oplus }}})\mathrm{if}\mathit{R}\mathit{>}\mathrm{3.8}{\mathit{R}}_{\mathrm{\oplus }}\end{array}$(16)plus a normal distribution centered around this value with dispersion σ_m ≃ 0.42. Both slopes are very similar to those proposed by Weiss & Marcy (2014) and Lissauer et al. (2011). The value of R_crit was fine-tuned to obtain a smooth transition between both segments. Equation (16) henceforth allows us to perform Monte Carlo simulations assigning masses to different observed systems with given values of R_i.

To apply our model statistically we still require knowledge of the disk parameters. Probable values are also weakly constrained. Flare index are expected to be strongly dependent on heating and cooling efficiencies and radiative properties of the disk, and hydro-simulations have registered values on the order of f ~ 0.3 (e.g., Bitsch et al. 2013). Similarly, although classical estimates for the minimum mass solar nebula (MMSN) point towards values of α ~ 3 / 2, there is no obvious reason to believe such a value should be universal. The observed distribution of exoplanetary systems shows a wide diversity of power-law index (e.g., Raymond & Cossou 2014), although the present location may not be indicative of their formation sites if planetary migration was wide spread. Additionally, the scale height H₀ is believed to be H₀ ∈ [ 0.01,0.1 ] and most studies assume values H₀ ≃ 0.05. However, since the resonance offset is a strong function of H₀ (in fact, Δ_{(p + q),p} ∝ 1 /H₀, see Eq. (14)), we kept its value as a free parameter.

Fig. 6 Blue columns: observed orbital period ratio P2/P1 for 36 exoplanetary pairs in the vicinity of the 2:1 resonance. Gray histograms result of 400 Monte Carlo simulations as obtained from Eq. (3) with random assigned values for masses and disk parameters (see text for details). Orange histograms previous values evolved through tidal effects for the duration of the lifespans of the systems.

Fig. 7 Same as previous figure but for systems in the vicinity of the 3:2 resonance.

Our Monte Carlo simulations were performed as follows. For each exoplanetary system in the vicinity of the 2:1 and 3:2 MMRs, we assigned each planet a physical radius R_i in accordance to the estimated values and errors given in the data set. The masses were then generated from expression (16) including a normal distribution with variance σ_m. In each run the disk properties (f,α,H₀) were also chosen randomly assuming normal distributions with mean (f_med,α_med,H_0med) and variances (σ_f,σ_α,σ_H0). The resulting parameters were then fed into Eq. (14) to obtain a value of the resonance offset. If the set of parameters did not lead to convergent migration and a stable resonance configuration (as described by the conditions deduced in Deck & Batygin 2015), we repeated the process. We ran a total of 400 successful cases were run for each system and plotted the distribution of P₂/P₁ as histograms.

Figure 6 shows results for 36 systems in the vicinity of the 2:1 resonance. Blue columns indicate the observed value of the orbital period ratio, while the gray histograms correspond to our Monte Carlo simulations with f = 0.25 ± 0.05, α = 1.0 ± 0.5 and H₀ = 0.03 ± 0.02. The mean values were chosen close to those estimated from Fig. 3, except for the disk scale height taken was slightly lower. The variances were chosen so as to include the values usually present in the literature (e.g., H₀ = 0.05 and α = 3 / 2 for the MMSN). Although these values may be considered arbitrary, they are consistent with current knowledge of protoplanetary disks.

Although the results look promising, some of the larger offsets, such as those observed for Kepler-126, Kepler-154, and Kepler-55, appear as low-probability events. However, we must also allow for resonance divergence through tidal evolution. We adopted the analytical prescription described in Papaloizou (2011), which includes tidal distortions in both the planets and central star. However, this introduces additional free parameters into our model. Without any information on their internal structure or chemical composition, it is not possible to give reliable values for the planetary tidal coefficient Q_p. Even within our own solar system these parameters are still under discussion (e.g., Lainey 2016), although some order of magnitude is generally agreed upon. Since we wanted to avoid fixing a single value for Q_p, we introduced the following empirical formula: $\log \mathrm{\left(}{\mathit{Q}}_{\mathrm{p}}\mathrm{\right)}\mathrm{=}\max \mathrm{(}\mathrm{2}\mathit{,}\mathrm{6}\mathrm{-}\mathrm{0.7}{\mathit{\rho }}_{\mathrm{p}}\mathrm{)}\mathit{,}$(17)where \begin{lxirformule}$\rho_{\rm p}$\end{lxirformule} is the volumetric density of the planet in units of gr/cm³. Even though this relation is not based on any physical evidence, it yields qualitatively correct results in the case of the Earth (Q_p ~ 10²) and Jupiter (Q_p ~ 10⁵); it also gives a smooth transition between both extremes.

Following Benítez-Llambay et al. (2011) we adopted a stellar tidal coefficient Q_∗ = 10⁶, while the stellar rotation was taken as P_∗ = 15 ± 5 days. However, stellar tides proved to be negligible in all cases and the results were found to be virtually insensitive to both Q_∗ and P_∗. Finally, few systems have estimated ages T_age, so we assumed random values following an uniform distribution in the interval T_age ∈ [ 2.5,7.5 ] Gyr. All planets were initially placed in spin-orbit equilibrium (pseudo-synchronous rotation) regardless of their distance from the star.

New final values for the offsets, incorporating tidal evolution, are shown as orange histograms in Fig. 6. When tidal effects are negligible, the orange values appear to be superimposed on the original gray histograms (e.g., HIP-41378, KE-24). In other systems, such as Kepler-256, Kepler-327, or Kepler-55, tidal evolution causes a significant increase in offset leading to final values that are much closer to those currently observed. Other systems, however, continue to be difficult to explain by this mechanism (e.g., Kepler-126, Kepler-154, Kepler-332, and Kepler-363). Most of these are characterized by very large offsets on the order of Δ_{2 / 1} ≃ 0.08, and could in fact constitute part of the background secular population, which appear overlaid to the (near)-resonant systems.

Three other problematic systems are worth mentioning. Kepler-450 has a ratio R₂/R₁ ≃ 3 leading in most Monte Carlo runs to m₂/m₁ ≃ 10. As seen from Eq. (14), for such systems our model inexorably predicts very small offsets, much smaller than the observed value of Δ_{2 / 1} ≃ 0.04. An opposite result perceived for Kepler-9, where our model always seems to predict larger offsets than observed. This is not surprising since the bodies in this system have assigned masses in the giant-planet range and should have undergone Type-II migration. The inability of our model to predict compatible values is thus expected. Another complicated system is Kepler-84, which has a ~48% error in the determination of the stellar radius; this error can lead to a huge mistake in the process of measuring planetary radii.

In general terms, however, our model appears to work well. For most real systems it predicts offsets similar to the observed values, even if the agreement does not occur for the peak of the simulated distribution. We do not believe this is an important issue since we do not know the actual disk parameters that housed each individual planets. Consequently, the fact that we are able to model the diversity and approximate magnitude of Δ_{2 / 1} for so many real cases is, in our opinion, a very positive result. It is important to mention that, as proposed by Silburt & Rein (2015), tidal effects cannot by themselves explain the (near)-resonant configurations. However, they are an important additional mechanism, at least in some planets.

Figure 7 shows results for 36 systems in the vicinity of the 3:2 mean-motion resonance. As before, our model underestimates the offset for a few systems (e.g., Kepler-102, Kepler-131, Kepler-363, and Kepler-393), but again these could be background objects whose position in the vicinity of the resonant domain is merely coincidental. The vast majority of the cases, however, are well reproduced by the model and in several examples the observed value closely matches the peak of the histograms.

As an additional test, we analyzed the sensitivity of these results with respect to the model parameters. We found that the distribution of the offsets was only weakly dependent on the disk flare and surface density profile, although they did show a stronger relation with both H₀ and Q_e. Higher values of both these quantities led to lower Δ_{(p + q) /p}, so more precise estimations in the future could help further evaluate the merits of this approach.

Finally, although in this paper we concentrated on the (near)-resonant population of exoplanets around the 2:1 and 3:2 MMRs, it is important to keep in mind that most of the close-in systems lie in secular configurations. While explaining the origin of such systems lies beyond the scope of the present work, it is possible that at least some of them could be the result of failed resonance captures or the outcome of collisional or planet-scattering events. Not all mass ratios or disk parameters necessarily lead to stable resonance trapping (e.g., Deck & Batygin 2015), and the consequent chaotic reshuffling of the orbits could explain at least some of our (relatively few) negative results.

6. Conclusions

We have presented a very simple model for the resonance offset Δ_{(p + 1) /p} generated by Type-I migration in a flared disk and applicable to first-order resonances. We described planetary migration using analytical prescriptions with correcting factors that were fine-tuned with hydrodynamical simulations of two-planet resonant systems. We also included the post-formation orbital evolution from tidal effects, both in the star and planetary bodies.

We found that the observed distribution of exoplanetary systems in the vicinity of the 2:1 resonance presents compelling evidence of an increase in the value of Δ_{2 / 1} closer to the central star. Planetary pairs near the 3:2 resonance show much weaker dependence and much smaller offsets for all distances. Our model adequately reproduces both these properties and the general distribution of the (near)-resonant Kepler systems is consistent with smooth migration in a thin laminar flared disk with significant flare (perhaps up to f ~ 0.25) and very thin disks (H₀ ~ 0.03).

Validation of this model proved difficult since most observed systems lack reliable estimations of the planetary masses. However, in most cases Monte Carlo simulations with an updated empirical mass-ratio relation and random values for the disk parameters, led to distributions of the resonance offsets similar to the observed values. Although some cases were not adequately explained, indicating that our model is still insufficient, we were able to reproduce values of Δ_{(p + 1) /p} that were at least qualitatively similar and, in some cases, very close to the detected magnitudes.

Finally we found that large offsets, even on the order of Δ_{2 / 1} ≃ 0.1, are still associated with libration of one or both critical arguments. Without a more sophisticated resonance model we cannot say whether this libration is dynamical (i.e., inside

the separatrix) or kinematic (i.e., outside the main resonance domain); however, it is possible that several of these (near)-resonant systems may be resonant after all.

Acknowledgments

This work has been supported by research grants from ANCyT, CONICET, and Secyt-UNC. The authors are grateful to an anonymous referee for inspiring questions that helped improve the work. We also wish to thank IATE and CCAD (Universidad Nacional de Córdoba) for extensive use of their computational facilities.

References

All Figures

Fig. 1 Distribution of orbital period ratios of adjacent pair of planets as a function of the orbital period of the inner planet in the vicinity of the 2:1 MMR (top) and 3:2 MMR (bottom). Data obtained from exoplanet.eu. The red circles identify confirmed planets detected by transits or TTV, while the blue circles correspond to recently validated Kepler planets proposed by Morton et al. (2016). Planets discovered by any other method are depicted with open circles. In the top panel, open green circles denote the GJ876 system, while the dashed black circle groups the resonant (or near-resonant) members of HD 73526, HD 82943, HD 155358, and HD 92788.

In the text

Fig. 2 Dynamical maps of max(Δe) for the 2:1 (left) and 3:2 MMR (right) for fictitious two-planet systems with m1 = 0.05mJup and m2 = 0.10mJup orbiting a central star with mass m0 = 1m⊙. The orbit of the outer planet was initially circular with a2 = 1 AU, and all angular variables where chosen equal to zero. Total integration time was 103 yr. The green lines indicate the location of the zero-amplitude ACR-type librational solutions, estimated from the simple analytical model (see Eq. (14)). The white dots in the left panel are the result of three N-body simulations of resonance trapping in the 2:1 resonance. The broad black curve in the right plot corresponds to the Marchal-Bozis stability limit. See text for details.

In the text

Fig. 3 Orbital period ratio (p + 1) /p + Δ(p + 1) /p as a function of the period of the inner planet, as predicted by expression (14) for flared disks with H0 = 0.05, α = 1, and f = 0.25. Top panels: represent the 2:1 MMR while bottom plots represent the 3:2 resonance. Each graph assumes a different value for the mass (in units of m⊕) of the inner planet, whose value is indicated in the top right-hand corner. Colored curves correspond to different mass ratios: m2/m1 = 1.6 (green), m2/m1 = 1.7 (blue), m2/m1 = 1.8 (purple), and m2/m1 = 2 (orange). Open black circles denote the distribution of confirmed planets.

In the text

Fig. 4 Top left panel: result of four N-body simulations leading to capture in the 2:1 MMR. In all cases m1 = 20m⊕. Color lines indicate different mass ratios: m2/m1 = 2 (orange), m2/m1 = 1.8 (purple), m2/m1 = 1.7 (blue), and m2/m1 = 1.6 (green). Remaining panels: different dynamical characteristics of the orbital evolution within the resonance of the run with m2/m1 = 1.6. The black continuous lines in the bottom panels represent the analytical predictions given by expressions (3).

In the text

Fig. 5 Top: distribution of radii between confirmed planets in the vicinity of the 2:1 (blue) and 3:2 (red) MMR. Dashed line corresponds to equal-size bodies (R1 = R2). Middle: radius as a function of the orbital period. Light gray circles indicate the values of planets outside both resonances. Bottom: planetary mass mi vs. radius Ri for all planets (bodies larger than Jupiter shown in gray). Dashed orange line corresponds to the so-called R2.06-law by Lissauer et al. (2011), while the fit proposed by Weiss & Marcy (2014) is shown as a continuous orange line. Our own fit, which is calculated for all planets smaller than Jupiter, is presented as a thick green line. Values within 1σdis dispersion are contained within the light green region.

In the text

	Fig. 6 Blue columns: observed orbital period ratio P2/P1 for 36 exoplanetary pairs in the vicinity of the 2:1 resonance. Gray histograms result of 400 Monte Carlo simulations as obtained from Eq. (3) with random assigned values for masses and disk parameters (see text for details). Orange histograms previous values evolved through tidal effects for the duration of the lifespans of the systems.
In the text

	Fig. 7 Same as previous figure but for systems in the vicinity of the 3:2 resonance.
In the text