Free Access
 Issue A&A Volume 602, June 2017 A101 11 Planets and planetary systems https://doi.org/10.1051/0004-6361/201629642 22 June 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 m1 and m2 orbiting a star m0, orbital periods P1<P2, and in the vicinity of a first-order (p + 1) /p MMR. We can then define the so-called resonance offset as (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. 1Distribution 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. Open with DEXTER

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)/Pi for adjacent planets in multiplanetary systems, as function of the orbital period Pi 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 104 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 P2/P1 ∈ [ 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 P1 ≤ 15 days and a second subset with P1 ≥ 15 days. These authors attributed this dichotomy to tidal effects, which should be important for orbital periods below (say) P1 ≤ 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 P1 ~ 1 day, down to Δ2 / 1 ~ 0.01 for orbital periods P1 ~ 103 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 P1 ~ 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 P2/P1 ∈ [ 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 P1 ~ 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. 2Dynamical 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. Open with DEXTER

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 (P2/P1,e1) plane. The masses of the planets were chosen equal to m1 = 0.05mJup and m2 = 0.10mJup, and the central star assumed of solar mass. We specifically chose a mass ratio m2/m1> 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 a2 = 1 AU and all angles equal to zero.

The color code corresponds to the maximum values of | e1(t)−e1(t = 0) | (denoted here as maxe)) attained during a 103 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 maxe) 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 (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 e1, along with a relationship between the eccentricities of both planets. These give (3)(see Lee 2004), where the coefficient Ci depends solely on the semimajor axes ratio a1/a2. For the MMRs under consideration, they take the following values: C1 ≃ 1.5,C2 ≃ 0.29, for the 2:1 resonance; and C1 ≃ 1.2,C2 ≃ 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 ei 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 e1 ~ 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 (4)where the new coefficients (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 (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 m1/m2, expression (6) shows that the value of the offset is a linear function of m2 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): (7)where ri is the radial vector and vi 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 a1 was chosen outside the 2:1 MMR and all angles equal to zero. The orbital decay timescales were τa1 = 105 and τa2 = 7 × 104 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 105 yr. We assumed that the disk began dispersing at Tdisp = 7 × 104 yr and its surface density reached zero at Tdisp = 8 × 104 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 e1 ≃ 0.03 and P2/P1 ≃ 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 (e1 ≃ 0.01) and a slightly larger offset (P2/P1 ≃ 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., P2/P1 ≃ 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 104 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 (8)where Hri is the disk aspect-ratio in the position of each planet. The value Qe is a constant and Qa = Qa(α) is a function of the slope α of the surface density profile. The linear model by Tanaka et al. (2002) found , while numerical fits with hydro-simulations usually lead to slightly different functional forms (e.g., D’Angelo & Lubow 2010). Finally, (9)is the typical timescale for planetary migration. Here is the (circular) Keplerian angular velocity and the gravitational constant.

We assume a laminar disk whose surface density Σ and the aspect-ratio Hr have a power-law dependence with the distance from the star, (10)where Σ0 and H0 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 Qe 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 Qe ~ 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 (11)Assuming flat disks (f = 0) and typical values of H0 ~ 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 102 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(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, (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

(14)

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

 Fig. 3Orbital 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. Open with DEXTER

## 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 m2/m1 between the planets must satisfy the condition (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 H0 = 0.05, α = 1 and f = 0.25, in the form of the predicted equilibrium values of the orbital period ratio P2/P1 = (p + 1) /p + Δ(p + 1) /p as a function of P1. 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 m1, while each color curve corresponds to a different value of the mass ratio m2/m1.

In all cases the offset increases with smaller semimajor axes, although its magnitude is a strong function of both m1 and m2/m1. As expected from Eq. (15), mass ratios close to (m2/m1)min show high values of Δ(p + 1) /p, although (also expected) its value increases linearly with m1. 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 m1 ~ 10m. Although this limiting case is marginally adequate, a much larger diversity in values is obtained assuming m1 = 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 m1 and the mass ratio. The best results are obtained assuming m1 ~ 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, H0 = 0.05, and a surface density of Σ0 = 400 gr/cm2. The factors Qa and Qe were chosen as described at the end of Sect. 3.2. In all cases we simulated capture into the 2:1 MMR with m1 = 20m and different values for the mass ratio m2/m1. 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 m2/m1 = 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 P2/P1, where we have neglected the data points prior to resonance lock. The main resonant angle θ1 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 θ2 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 P2/P1 (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. 4Top 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). Open with DEXTER

## 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 m1 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. 5Top: 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. Open with DEXTER

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 (R2 versus R1) 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 RJup. Although in most cases R2/R1> 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 Ri ~ 2 R for Pi< 10 days, in contrast to larger distances where larger planets abound. The decrease in Ri closer to the star, where the observed resonance offset is larger, appears contrary to expectations from our model (14), which indicates that Δ(p + q) /pmi/m0. 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 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 Ri> 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<Rcrit and one for larger bodies, where the value of Rcrit 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 (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 Rcrit 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 Ri.

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 H0 is believed to be H0 ∈ [ 0.01,0.1 ] and most studies assume values H0 ≃ 0.05. However, since the resonance offset is a strong function of H0 (in fact, Δ(p + q),p ∝ 1 /H0, see Eq. (14)), we kept its value as a free parameter.

 Fig. 6Blue 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. Open with DEXTER

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

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 Ri 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,α,H0) were also chosen randomly assuming normal distributions with mean (fmed,αmed,H0med) 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 P2/P1 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 H0 = 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., H0 = 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 Qp. 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 Qp, we introduced the following empirical formula: (17)where \begin{lxirformule}$\rho_{\rm p}$\end{lxirformule} is the volumetric density of the planet in units of gr/cm3. Even though this relation is not based on any physical evidence, it yields qualitatively correct results in the case of the Earth (Qp ~ 102) and Jupiter (Qp ~ 105); it also gives a smooth transition between both extremes.

Following Benítez-Llambay et al. (2011) we adopted a stellar tidal coefficient Q = 106, 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 Tage, so we assumed random values following an uniform distribution in the interval Tage ∈ [ 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 R2/R1 ≃ 3 leading in most Monte Carlo runs to m2/m1 ≃ 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 H0 and Qe. 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 (H0 ~ 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.

## All Figures

 Fig. 1Distribution 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. Open with DEXTER In the text
 Fig. 2Dynamical 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. Open with DEXTER In the text
 Fig. 3Orbital 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. Open with DEXTER In the text
 Fig. 4Top 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). Open with DEXTER In the text
 Fig. 5Top: 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. Open with DEXTER In the text
 Fig. 6Blue 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. Open with DEXTER In the text
 Fig. 7Same as previous figure but for systems in the vicinity of the 3:2 resonance. Open with DEXTER In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.