EFFECTS OF TURBULENCE, ECCENTRICITY DAMPING, AND MIGRATION RATE ON THE CAPTURE OF PLANETS INTO MEAN MOTION RESONANCE

In this section, we consider the direct numerical integration of migrating planetary systems, i.e., we integrate the full set of 18 phase space variables for the three-body problem consisting of two migrating planets orbiting a central star. For most of our simulations, the planets are started in the same plane so that the dynamics is only two dimensional; however, we have also run cases that explore all three spatial dimensions. The integrations are carried out using a Burlisch–Stoer (B–S) integration scheme (e.g., Press et al. 1992). In addition to gravity, we include forcing terms that represent inward migration, eccentricity damping, and turbulence. All three of these additional effects arise to the forces exerted on the planet(s) by a circumstellar disk. In this context, however, we do not model the disk directly, but rather include forcing terms to model its behavior, as described below.

To account for planet migration, we assume that the semimajor axis of the outer planet decreases with time according to the ansatz

$$\begin{equation} {1 \over a} {da \over dt} = - {1 \over { \tau _{a} }} \,, \end{equation} \tag{ 1 }$$

where τ_a is the migration timescale. Further, we assume that only the outer planet experiences torques from the circumstellar disk.

Small planets, those with masses smaller than that of Saturn, cannot clear gaps in their circumstellar disks and tend to migrate inward quickly in a process known as Type I migration (e.g., Ward 1997a, 1997b). A number of studies have shown that the Type I migration rate depends on the disk thermal properties and on local gradients of the gas density (e.g., Baruteau & Masset 2008; Paardekooper & Papaloizou 2008; Masset & Casoli 2009; Paardekooper et al. 2010). As a result, for some disks, Type I migration can be much slower (sometimes even directed outward) and a wide range of migration rates is possible. Larger bodies clear gaps and migrate more slowly. Estimates of the migration timescale for planets with a ∼ 1 AU typically fall in the range 10⁴–10⁶ yr (e.g., see Goldreich & Tremaine 1980; Papaloizou & Larwood 2000).

Planets are thus expected to experience a range of migration rates, depending on both planet masses and disk properties. Since we want to isolate the effects of the migration rate on entry into resonance, we adopt a purely parametric approach. We thus consider a wide range of migration rates, where the migration timescale varies over the range τ_a = 10³–10⁶ yr. Note that the shorter timescales are included here to study the physics of resonance capture (at these fast rates) and are not generally expected in most circumstellar disks. On the other hand, fast migration could occur for planets with mass ∼10 M_⊕ migrating within circumstellar disks that have sufficiently small aspect ratios (H/r < 0.03) and large masses (Masset & Papaloizou 2003).

In addition to inward migration, circumstellar disks also tend to damp orbital eccentricity e of the migrating planet. This damping is generally found in numerical simulations of the process (e.g., Lee & Peale 2002; Kley et al. 2004) and can be parameterized such that

$$\begin{equation} {1 \over e} {de \over dt} = - {1 \over { \tau _{e} }} = K \left({1 \over a} {da \over dt} \right) \quad {\rm so} \, \, {\rm that} \quad { \tau _{e} }= { \tau _{a} }/ K \,, \end{equation} \tag{ 2 }$$

where τ_e is the eccentricity damping timescale. Some analytic calculations suggest that eccentricity can be excited through the action of disk torques (Goldreich & Sari 2003; Ogilvie & Lubow 2003), although multiple planet systems would be compromised if this were always the case (Moorhead & Adams 2005). Additional calculations show that disks generally lead to both eccentricity damping and excitation, depending on the disk properties, gap shapes, and other variables (e.g., Moorhead & Adams 2008). The value of the damping parameter K can also be inferred from hydrodynamical simulations, which predict a range of values. Studies of resonant systems (Kley et al. 2004) advocate K values of order unity. In isothermal disk models, however, K ∼ 10 for typical cases (e.g., Cresswell & Nelson 2008). More recent work indicates that in radiative disk models, the eccentricity damping parameter K can be as large as 100 (Bitsch & Kley 2010).

In spite of the aforementioned uncertainties, this study focuses on the case of pure damping, adopts fixed values of K for a given simulation, and considers its effects on the dynamics of mean motion resonances. We expect that the inclusion of damping will act to enhance the survival of mean motion resonances (Lecoanet et al. 2009). Using the ansatz of Equation (2), this study considers a wide range of the damping parameter K such that 0 ⩽ K ⩽ 100, where we consider the cases K = 1 and K = 10 as our "standard" values.

Turbulence is included by applying discrete velocity perturbations at regular time intervals; for the sake of definiteness, the forcing intervals are chosen to be twice the orbital period of the outer planet (four times the period of the inner planet for systems with the 2:1 period ratio). Both components of velocity in the plane of the orbit are perturbed randomly, but the vertical component of velocity is not changed. The amplitude of the velocity perturbations thus represents one of the variables that characterize the system. These amplitudes are chosen to be consistent with the expected torques, as described below.

The torques due to turbulent fluctuations have been studied previously using MHD simulations (e.g., Nelson & Papaloizou 2004; Laughlin et al. 2004; Nelson 2005; Oishi et al. 2007), and these results can be used to estimate the range of amplitudes. The torque exerted on a planet by a circumstellar disk can be expressed as a fraction of the benchmark torque T_D = 2πGΣrm_P, where Σ is the surface density of the disk, r is the radial coordinate, and m_P is the planet mass (Johnson et al. 2006). The scale T_D thus serves as a maximum torque in this problem. The amplitude of the expected angular momentum fluctuations is then given by ΔJ = f_TΓ_RT_Dt_T, where f_T is the fraction of torque scale T_D realized by the disk, Γ_R is a reduction factor due to planets creating gaps in the disk, and t_T is the time required for the disk to produce an independent realization of the turbulence. Previous work suggests that f_T ∼ 0.05 (Nelson 2005), Γ_R ∼ 0.1 (Adams et al. 2008), and t_T is comparable to the orbit time of the outer planet (Laughlin et al. 2004; Nelson 2005). Including all of these factors, we expect that [(ΔJ)/J] ∼ 10⁻⁴ under typical conditions (a disk mass of ∼0.05 M_☉ with well-developed MRI turbulence such that α ∼ 10⁻³). Under some circumstances, the equatorial plane of the disk is not sufficiently ionized to support MRI turbulence and a dead zone develops; in this case, the fraction f_T would be dramatically decreased. Given the uncertainties in turbulent behavior, and the wide range of possible disk conditions, the fluctuation amplitude could vary by an order of magnitude (perhaps more) in either direction. As a result, we consider turbulent fluctuation amplitudes in the range 0 ⩽ [(ΔJ)/J] ⩽ 10⁻³.

For a given realization of the migration scenario, we need to determine whether or not the system resides in mean motion resonance. First, we determine the ratio of the orbital periods of the two planets. It is straightforward to determine when systems have nearly integer period ratios and this condition can be used as a proxy for being in a mean motion resonance. However, this condition is necessary but not sufficient, so we must also monitor the relevant resonance angles (Murray & Dermott 1999; hereafter MD99). For first-order resonances, these angles have the form

$$\begin{equation} \theta _1 = (j + 1) \lambda _2 - j \lambda _1 - \varpi _1 \,, \end{equation} \tag{ 3 }$$

$$\begin{equation} \theta _2 = (j + 1) \lambda _2 - j \lambda _1 - \varpi _2 \,, \end{equation} \tag{ 4 }$$

$$\begin{equation} \theta _3 = \varpi _1 - \varpi _2 \,. \end{equation} \tag{ 5 }$$

For second-order resonances, these angles take the form

$$\begin{equation} \theta _1 = (j + 2) \lambda _2 - j \lambda _1 - 2 \varpi _1 \,, \end{equation} \tag{ 6 }$$

$$\begin{equation} \theta _2 = (j + 2) \lambda _2 - j \lambda _1 - \varpi _1 - \varpi _2 \,, \end{equation} \tag{ 7 }$$

$$\begin{equation} \theta _3 = (j + 2) \lambda _2 - j \lambda _1 - 2 \varpi _2 \,, \end{equation} \tag{ 8 }$$

$$\begin{equation} \theta _4 = \varpi _1 - \varpi _2 \,. \end{equation} \tag{ 9 }$$

In order to monitor the angles and determine if the system is in a resonant state, we must choose the appropriate time windows. Note that the resonance angle θ₀ = ϖ₁ − ϖ₂ oscillates on a much longer timescale than the other angles (where θ₀ = θ₃ (θ₄) for first (second) order resonances). As a result, the angle θ₀ is measured over a time period corresponding to 1500 orbits of the outer planet, whereas the other angles (which oscillate faster) are monitored over a time window of 300 orbits of the outer planet. These timescales are chosen to be (roughly) several times the expected libration periods of the angles (and the expected libration periods can be calculated from the restricted three-body problem—see MD99). Each angle is considered to be in libration if its value stays bounded within 120° of the effective stability point for the time periods given above. In this context, the effective stability point is determined by the mean value of the angle over the given time window for monitoring; note that these systems are highly interactive (e.g., due to turbulent forcing) so that the stability points are not fixed. Notice also that the value of 120° was chosen arbitrarily. If any of the angles θ_i obtain a value greater than 120°, measured from either side of the effective stability point, then that angle is considered not to be in resonance. The code continues to monitor all of the relevant angles for the duration of the time when the periods have a well-defined ratio (of small integers). As a result, each resonance angle could go in and out of libration many times.

Given the formulation described above, we numerically study the entry of planets into mean motion resonance and the subsequent survival of the resonant configurations. These results depend on a number of parameters, including the migration rate, the eccentricity damping rate, the level of turbulence, and the planetary masses. As indicated by the semi-analytic models (Section 3, Quillen 2006), we expect the survival of mean motion resonance to be compromised with sufficiently fast migration rates. The introduction of turbulence can act to further reduce the ability of systems to stay in resonance (Adams et al. 2008), whereas eccentricity damping generally acts in the opposite direction by helping to maintain resonance (Lecoanet et al. 2009). The results also depend on the masses of the planets. As the masses increase, the systems become more highly interactive and mean motion resonance is harder to maintain.

For the first set of simulations, we begin with a standard two-planet system consisting of a Jovian mass planet m₁ = 1 M_J and a "super-Earth" with the mass m₂ = 10 M_⊕. The properties of this system are close to the restricted three-body problem and hence the resonance is expected to be described reasonably well by the pendulum model of MD99. The star is taken to have a mass M_* = 1.0 M_☉. The Jovian planet acts as the inner planet and begins with an eccentricity e₁ = 0.05 and a period of P₁ = 1000 days (so that a₁ = 1.96 AU). The smaller planet starts with an eccentricity e₁ = 0.10 and a semimajor axis of a₂ = 1.8 a₁, equivalent to a period ratio of P₂/P₁ = 2.4, which places the system comfortably outside the 2:1 mean motion resonance. Both of the planets are placed in the same orbital plane. As the outer planet migrates inward, it can (in principle) enter into the 2:1 resonance; if the migrating planet passes through the 2:1 resonance, it can then (potentially) enter into resonant states with smaller period ratios.

In this parametric study, we allow migration of the outer planet, given by Equation (1), to continue throughout the simulations. Since the inner Jovian planet is expected to open a gap in the circumstellar disk however, the migration rate could be altered, where the variations depend on the gap structure. Although not considered herein, some disks with gaps can even halt migration altogether and produce planet traps (Masset et al. 2006). In addition, since the (smaller) outer planet often acquires substantial eccentricity, it will move in and out of the gap over the course of its orbit. This effect leads to time-dependent migration torques that vary on the orbital timescale; the time variations tend to average out over the libration timescale of the resonances, but the migration rate could be altered slightly.

Because these systems are highly chaotic, different realizations of the problem lead to different outcomes. For each set of parameters, we perform an ensemble of (at least) 1000 effectively equivalent simulations, where the simulations differ only by the relative position of the two planets in their orbits and by the relative angle between the orientation of the two orbits (i.e., the arguments of periastron ϖ₂ − ϖ₁). The length of the numerical integration t_T is set by the migration timescale τ_a, such that t_T = τ_a for the slowest migration rate (τ_a = 10⁶ yr) and t_T ≈ 10τ_a for the fastest migration rate (τ_a = 10³ yr). The overall integration times are thus shorter for the faster migration rates, but remain long enough to encompass many libration timescales for the relevant resonance angles.

2.2.1. Time Evolution and Resonance Criteria

For this set of system parameters, Figure 1 illustrates the basic time evolution for an ensemble of planetary systems. Here, the fractions of the systems that reside in mean motion resonance are plotted as a function of time for a moderately short migration timescale τ_a = 2 × 10⁴ yr. The simulations shown in the top panels include eccentricity damping with parameter K = 1; the two panels on the left include turbulence with the standard level of fluctuations (ΔJ)/J ∼ 2 × 10⁻⁴. The various curves in each panel correspond to resonances with period ratios of 2:1 (blue), 5:3 (red), and 3:2 (green). The systems are considered to be in resonance if the period ratios are near the relevant integer values and if any of the resonance angles are librating (see Section 2.1 and Equations (3)–(9)). Note that this migration rate was used because slower migration rates lead to few systems in the 5:3 and 3:2 resonances. Since all of the systems start outside of resonance, the fractions start at zero and increase with time as migration pushes the planets together. The 2:1 resonance is encountered first, so that corresponding fraction grows first. As the systems evolve, resonance is often compromised, so that the fractions reach a peak value and then decrease. After some of the systems leave the 2:1 state, they become locked into the 5:3 resonance, and then sometimes the 3:2 resonance. As a result, the peak fraction occurs later for resonances that are further inward, and the peak is lower for the weaker resonances (as expected). When systems begin to leave resonance, some of them decay by losing a planet through ejection, accretion onto the star, or collision with the other planet (the probabilities of these end states are quantified in Section 2.3). This effect is illustrated in Figure 1 by the black curves, which show the fraction of systems that retain both planets as a function of time.

Zoom In Zoom Out Reset image size

We note that planetary systems are often said to "be in resonance" according to different criteria. This trend is illustrated in Figure 2 for the case of the 2:1 mean motion resonance. As shown in Figure 1, the fraction of systems that reside in resonance is a function of time for a given migration rate. The peak value of this time-dependent curve can be used as one measure of the fraction of systems that are in resonance. However, systems enter and leave resonance at different times, so that the total fraction of systems that enter resonance will be larger than the maximum fraction that reside in resonance at a given time (the peak of this curve). A necessary (but not sufficient) condition for a system to be bound into resonance is for the ratio of the periods to be near 2:1. In this context, the period ratio is "near" 2:1 if |P₂/P₁ − 2| ⩽ 0.01, where we discuss this issue more quantitatively below. As outlined above, we first invoke the constraint that P₂/P₁ ≈ 2. This fraction is shown as the dashed blue curves marked by squares in Figure 2. The four panels show the effects of including eccentricity damping (with K = 1, panels in top row) and turbulent forcing (with (ΔJ)/J ∼ 2 × 10⁻⁴, panels on left side). Next, we note that each of the resonance angles can be either librating or circulating. For those that are librating, the range of angles (the libration width) is highly variable. For this paper, we use the requirement that the resonance angles are confined to be within 120° of the effective stability point (as defined above). With this specification, the corresponding fractions of systems in resonance are shown as the cyan curves (θ₁ angle), the red curves (θ₂ angle), and the green curves (θ₃ angle). The solid curves in each panel show the fraction of systems for which any of the resonance angles are librating at the end of the migration epoch. Note that only a relatively small fraction of the systems maintain resonance for the entire migration epoch. In addition, the inclusion of eccentricity damping (top panels) is crucial for the survival of resonant states.

Zoom In Zoom Out Reset image size

For completeness, we note that the curves shown in Figure 2 have slightly different meanings for the different resonance angles. In order for any one of the angles to be considered in resonance, it must librate over (approximately) three libration periods. However, these periods are not the same for the three angles. In particular, the libration period for θ₃ is much longer than the other two. In addition, for this class of systems, the orbit of the outer (lighter) planet varies much more than that of the inner Jovian planet. As a result, the argument of periastron of the outer planet can circulate on a long timescale, but the resonance angle θ₂ can still be considered (according to the criteria used here) to be librating.

In order to understand how the period ratios vary, we monitor the period ratio for systems that are found in the 2:1 mean motion resonance. Monitoring is triggered by the condition that P₂/P₁ < 2.05; however, once triggered, this bound is relaxed and the period ratio for systems in resonance can take any value as long as the angles are librating (see above). We find that systems typically exhibit both a slight offset from exact commensurability and variations about this offset. The offset is typically less than ∼1% and the standard deviation is ∼2%. We note that offsets and variations of this magnitude are expected, given the size of the terms in the disturbing function, and hence the size of the non-Keplerian velocities due to resonance. Both the offset value and amount of variation depend on the levels of damping and turbulence, and on the duration of resonance. In the absence of turbulence, we find an offset such that P₂/P₁ ∼ 2.008. For systems that do not include eccentricity damping, the variation of the period ratio σ ∼ 0.04, but decreases for systems that stay in resonance over long times (10 times the migration timescale τ_a). For cases with eccentricity damping parameter K = 10, the resonances are longer lasting, and we find σ ∼ 0.015. For systems that include turbulent forcing, the period ratio P₂/P₁ ∼ 2.007 with σ ∼ 0.07 for short-lived resonances, but decreases to P₂/P₁ ∼ 2.002 with σ ∼ 0.015 for longer lived resonances.

2.2.2. Resonance Survival

Figure 3 shows the effects of both turbulence and eccentricity damping on the survival of resonances as a function of migration rate. In this case, we define resonance using the requirement that the planets have nearly integer period ratios and the libration width is less than 120° for any of the resonant angles. The lower right panel shows the survival of resonances as a function of migration rate with no eccentricity damping and no turbulence. This case is thus analogous to the model equations derived in Section 3 below. As expected, systems tend to leave resonance if the migration rate becomes too large. Here, systems leave the 2:1 resonance (blue curve) when the migration rate exceeds roughly 2 × 10⁻⁴ yr⁻¹ (τ_a = 5000 yr). After leaving the 2:1 resonances, systems can become locked into the 5:3 resonance (red curve) and/or the 3:2 resonance (green curve). Note that the curve for 2:1 resonances shows a broad maximum near the migration rate of 10⁻⁵ yr⁻¹ (τ_a = 0.1 Myr), with decreasing probability toward both slower and faster rates. The decrease with increasing migration rate is expected. The decrease toward slower migration rates occurs because some of the systems are locked into higher order resonances, which include the 7:3 and 9:4 mean motion resonances (these fractions are not shown in the figure). With the starting period ratio of 2.4, the systems must pass through these states to reach the 2:1 resonance; with extremely slow migration rates, these weak resonances can (sometimes) survive and thus reduce the probability of the systems entering the 2:1 resonance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The effects of including turbulent fluctuations are shown by the analogous curves in the lower left panel. Turbulence only has a chance to act on long timescales, so that the simulations with long migration times (low migration rates) are affected the most. More specifically, for migration rates slower than about 10⁻⁵ yr⁻¹ (τ_a = 0.1 Myr), turbulence has time to act, and the probability of maintaining a resonant configuration is lower, as shown on the left-hand side of the plot. We note that with the inclusion of turbulence, the weak higher order resonances (7:3 and 9:4) generally do not survive (unlike the case of no turbulence in the lower right panel).

The effects of including eccentricity damping is shown by the top right panel, where we have taken K = 1 (so that eccentricity is damped on the same timescale that migration takes place; see Equation (2)). The inclusion of this damping effect acts to preserve resonance—note that all of the survival fractions are higher when ${\dot{e}} \ne 0$ than in the absence of damping. This effect is especially important for the long-term survival of the resonant states (compare the solid curves in the top panels with those in the bottom panels), especially for the case of the 2:1 resonance. The survival probabilities of the (weaker) 5:3 and 3:2 resonances are also enhanced by the inclusion of eccentricity damping, but the absolute values of these probabilities remain low. Keep in mind that these results correspond to K = 1; larger eccentricity damping rates lead to more dramatic consequences (see below).

When both turbulent forcing and eccentricity damping are included, we obtain the results shown in the upper left panel of Figure 3. In this case, the effects of turbulence dominate at low migration rates, so that fewer resonant systems survive. At high migration rates, however, turbulence does not have sufficient time to act and the effects of eccentricity damping lead to a net gain in the survival fractions. For migration rates faster than about 3 × 10⁻⁵ yr⁻¹ (τ_a ≈ 0.033 Myr), eccentricity damping dominates over the effects of turbulence, so that more resonant systems survive.

For comparison, we consider the survival of resonant systems for the case with eccentricity damping parameter K = 10. These results are shown in Figure 4, where all of the system parameters are the same as in Figure 3 except for the larger rate of eccentricity damping. As expected (e.g., Lecoanet et al. 2009), the simulations with K = 10 result in a larger survival rate than the corresponding cases with K = 1 (compare the upper right panels of Figures 4 and 3). For slower migration rates, where the dominant outcome is the 2:1 resonance (blue curves), the survival rate increases only modestly, from P_b ∼ 0.6 to P_b ∼ 0.8 with increasing values of K (for the case with no turbulence). For higher migration rates, the 3:2 resonance is most common state, and the survival rate increases substantially for the K = 10 case (compared to K = 1 systems). For the simulations that include turbulent fluctuations, however, the differences in survival fractions for the 2:1 resonance between the K = 10 and K = 1 cases are minimal (compare the upper left panels of Figures 3 and 4). In the absence of turbulence, the increase in resonance survival (for larger K) arises most strongly at slow migration rates; however, the regime of slow migration is where turbulence has enough time to act, and hence to compromise resonant states. For 3:2 resonances, which arise primarily at fast migration rates where turbulence does not have enough time to act, the increased eccentricity damping rate leads to substantially larger survival fractions.

Before leaving this section, we note that the simulations shown thus far all start with the planets confined to the orbital plane. In order to see how nonzero inclination angles affect the results, we have carried out a series of test simulations where the planets are given non-zero displacements in the vertical dimension in their starting states (so that these simulations are fully three dimensional). The results of these test simulations indicate that the third dimension is unimportant as long as the initial departures from the plane are not too large. More specifically, the starting vertical coordinates z₀ are uniformly sampled within the range [ − H, H], where H is the scale height of the disk. These test simulations use a variety of scale heights, with H/r = 0, 0.05, 0.10, and 0.20; we also use a 10 M_⊕ outer planet, eccentricity damping parameter K = 10, and our standard level of turbulence. For this choice of parameters, the results are virtually unchanged.

2.2.3. Effects of Eccentricity Damping and Turbulence

Next, we consider the case of eccentricity damping acting alone. Figure 5 shows the survival probabilities as a function of migration rate for four values of the eccentricity damping timescale, where the parameter K = 0.1, 1, 10, and 100. Taken together, the four panels of Figure 5 show that eccentricity damping acts to increase the fraction of systems that remain in mean motion resonance. The effect is most pronounced for the 2:1 resonance and for slow migration rates. In the regime of slow migration, a significant fraction of the systems leave the 2:1 resonance, presumably through the excitation of eccentricity via planet–planet interactions (Adams & Laughlin 2003; Moorhead & Adams 2005; Chatterjee et al. 2008; Ford & Rasio 2008; Matsumura et al. 2010). The inclusion of eccentricity damping counteracts this excitation and allows more systems to remain in resonance.

Zoom In Zoom Out Reset image size

For sufficiently large eccentricity damping rates (characterized by K = 100), essentially all systems remain in the 2:1 resonance until the migration rate exceeds a well-defined value, found numerically to be $|{\dot{a}}|/a \sim 3 \times 10^{-5}$ yr⁻¹ or τ_a ≈ 0.033 Myr (shown in the lower right panel of Figure 5). These results show that the loss of resonant states for the other cases (5:3 and 3:2) also occurs at well-defined values of the migration rate. In addition, as a rough approximation, the migration rates at which these three resonances are compromised are found to be evenly spaced logarithmically (by factors of ∼3). This behavior can be understood in a qualitative manner through simple physical considerations (see below) and through model equations (Section 3). Finally, we note that the 5:3 resonance, which is second order and hence weak, is sparsely populated; as a result, many of the systems with migration timescales ∼3000 yr are not found in any resonance.

The basic clock that determines the dynamics of these planetary systems is set by the libration timescale of the resonance. For the simplest model of the resonance, that resulting from the circular restricted three-body problem, the frequency for external resonances is given by

$$\begin{eqnarray} \omega _0^2 &=& -3 j_1^2 {\cal C}_R n e^{|j_3|} \qquad {\rm where} \nonumber\\ {\cal C}_R &=& \left({m_P \over M_\ast } \right) n [ f_d (\alpha) + \alpha ^{-1} f_i (\alpha) ] \,, \end{eqnarray} \tag{ 10 }$$

where n is the mean motion of the outer planet, α = a₁/a₂, and the functions f_d(α) and f_i(α) are given by the Laplace coefficients (see Section 8.5 of MD99). For odd order resonances, the function f_d < 0, so that the corresponding frequencies are real. For even order resonances, f_d > 0, but the equilibrium angle is shifted by π, and the frequencies are again real (see MD99 for further discussion). Notice also that f_i is nonzero only for the 2:1 resonance. The integers j₁ and j₃ depend on the type of resonance. Although both integers are negative for the cases of interest, the libration timescale only depends on the absolute value. More specifically, the integer pair (|j₁|, |j₃|) takes on the values (1,1), (3,2), and (2,1) for the 2:1, 5:3, and 3:2 resonances, respectively. For the values e = 0.10 and μ = m_P/M_* = 10⁻³, as used in the numerical simulations, we find that ω²₀/(3μen²)≈ 2.0, 3.0, and 8.1 for the three resonances. For these parameter values, the square of the frequencies is spaced by factors of ∼2. This simple analytic result is thus in qualitative, but not quantitative, agreement with the numerical results. One should keep in mind that the orbital elements that enter into these formulae (e.g., e and n) vary over the course of the simulations, so that comparisons are complicated.

If the migration rate ${\dot{a}}/a$ were the only relevant variable, then one would expect that capture into resonance would be compromised at a fixed value of the dimensionless parameter aω₀/$|{\dot{a}}| = \omega _0 { \tau _{a} }$. The case that most closely meets this expectation is that of migration with no eccentricity damping and no turbulent forcing (shown in the lower right panel of Figure 3). For this class of simulations, systems tend to enter the 3:2 resonance states for higher migration rates than for the 2:1 resonances, where this trend is predicted (qualitatively) by the simple theory outlined above. However, the fraction of systems in 5:3 resonance is not larger than the fraction in 2:1 resonance at large migration rates, in spite of the 5:3 having a shorter libration period. The 5:3 resonance is generally weaker, in the sense of being easier to disrupt, than the first-order resonances and does not survive for large migration rates. In terms of survival of the resonances, shown by the solid curves, the fraction in 2:1 is generally larger for all migration rates due to its greater stability.

For the case of migration with large eccentricity damping rates (see the lower right panel of Figure 5), the probability of resonance survival shows the expected qualitative behavior: each of the resonances dominates (has the largest fraction) for a well-defined range of migration rates. The 2:1 resonance is by far the most important for migration timescales longer than about 10⁴ yr. For shorter timescales, there is a narrow window of migration rates where the fraction of systems in 5:3 resonance shows a peak, and then the 3:2 resonance dominates for faster migration rates. Although the ordering of these results is consistent with theoretical expectations, the maximum migration rates are spaced at larger intervals than the factors of $\sqrt{2}$ suggested by the above analysis. Here, the large eccentricity damping rates significantly change the dynamics and hence the numerical values. Nonetheless, the qualitative trend holds up.

As the levels of turbulence increase, systems have greater difficulty maintaining mean motion resonance. This trend is quantified by the simulations depicted in Figures 6 and 7. These numerical experiments are carried out using the standard case of a Jovian planet on the inside and an inward migrating "super-Earth" with mass m₂ = 10 M_⊕. The eccentricity damping rate is set at the standard values of K = 1 (Figure 6) and K = 10 (Figure 7). As the amplitude of the turbulent fluctuations increases (from upper left to lower right in both Figures), the general trend is for the fraction of systems in resonance to decline significantly. The 2:1 mean motion resonance, which is the strongest and the first to be encountered, is compromised for a sufficiently rapid migration rate. As the level of turbulence increases, the migration rate at which systems leave the 2:1 resonance becomes lower (the curves shift to the right in the figures). We also note that the destructive action of turbulence is more pronounced for the solid curves, i.e., for the fraction of systems that remain in resonance at the end of the migration time. Finally, as expected, we find that more resonant systems survive for larger rates of eccentricity damping (compare Figures 6 and 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With the initial conditions used herein, where the planets are started outside the 2:1 resonance, the faster 5:3 and 3:2 resonances are not affected as severely by the presence of turbulence. These other resonances only arise when the migration rate is rapid, so that the migration timescale is short and turbulence has little time to act. For the 5:3 and 3:2 resonances, the probability curves shown in Figure 6 decrease slowly with increasing turbulent amplitude. As expected, the largest effect arises for the largest turbulent amplitude [(ΔJ)/J]_k = 10⁻³, where the probability of remaining in any of the resonant states is extremely low at the end of the migration epoch; the fraction of systems not bound into resonance is marked by the solid cyan curve, which is close to unity for all migration rates. Note that the 3:2 resonance lasts the longest in the face of increasing turbulence. This apparent resilience arises because the 3:2 cases are only present for fast migration rates, the regime where turbulence has less time to act (it is not due to the increased durability of the resonance).

2.2.4. Equal Mass Planets

Next we consider the case of two equal mass planets, with m₁ = m₂ = M_J. The results for survival of mean motion resonance are shown in Figure 8. The panels on the left include the effects of turbulent forcing; the panels on the top include the effects of eccentricity damping, where the parameter K = 1 so that the eccentricity damping timescale is the same as the migration timescale. These results for two Jovian planets are significantly different than those shown in Figure 3 for the case of a lower mass outer planet. One important effect of higher planetary masses is to increase the levels of planet–planet interactions in the systems. This effect, in turn, leads to greater libration widths for systems that stay in resonance and a lower probability of remaining in a resonant state. As a result, the probability of the system residing in either the 5:3 or the 3:2 resonance is significantly lower than in the case of a less interactive system (compare Figures 3 and 8). On the other hand, the fraction of systems that remain in the 2:1 resonance is larger for the more interactive (two Jupiter) systems.

Zoom In Zoom Out Reset image size

During the course of the numerical integrations, the planetary systems can end their evolution in a variety of ways. In many cases, the systems remain bound together, even though mean motion resonance is often compromised as described above. In many other cases, however, planets can be lost through scattering encounters, through collisions with each other, or via accretion onto the central star. This section outlines the probabilities for each of these possible end states of these dynamical systems.

For eccentricity damping parameter K = 1, Figure 9 shows the likelihood of the planetary systems ending their evolution in various possible end states for the standard case with inner planet mass m₁ = 1 M_J and outer planet mass m₂ = 10 M_⊕. These probabilities are shown as a function of migration rate for ensembles of simulations with migration only, migration and eccentricity damping, migration and turbulence, and for simulations including all three effects. The evolution of these systems produces a wide variety of outcomes, including survival of both planets for the entire evolutionary time (shown by the black curves), ejection of a planet (green curves), accretion by the central star (red curves), and collisions between the planets (blue curves). As illustrated by the four panels in the figure, the corresponding probabilities depend sensitively on the migration rates, eccentricity damping rates, and the levels of turbulence.

Zoom In Zoom Out Reset image size

Figure 9 shows several trends. In general, the probability for both planets to survive tends to decrease with increasing values of the migration rate. This trend is expected because slow migration rates allow the systems to adjust as they evolve; these cases with slow migration systematically exhibit less overall action than cases with higher migration rates. One important exception to this trend arises for the case of slow migration rates, the inclusion of turbulent fluctuations, and no eccentricity damping (see the lower left panel in Figure 9). In this regime, migration timescales are long enough that turbulence has time to act, which leads to loss of mean motion resonance (see the previous section), a greater possibility of orbit crossing, and subsequent planetary ejection. For this class of systems, the outer planet has substantially less mass than the inner planet and is far more susceptible to being lost. The outer planet is removed through ejection, accretion onto the central star, and through collisions with the inner Jovian planet. Note that the first two of these channels dominate the third.

Figure 9 shows another trend: as the migration rate increases, the probability of losing a planet through ejection decreases, whereas the probability of losing a planet through accretion onto the star increases. One important physical property that determines the relative number of accretion events versus ejections is the location of the planet(s) in the gravitational potential well of the star at the end of the migration epoch (when the planets are likely to suffer close encounters). The depth of the stellar potential well at a = 1 AU is approximately (30 km s⁻¹)², whereas the depth of the potential well at the surface of Jupiter is (43 km s⁻¹)². These scales are thus comparable. For fast migration rates, the outer planet is able to push the inner planet somewhat farther inward, deeper into the stellar potential well, and hence the probability of ejection decreases.

Figure 10 shows the analogous plots for the channels of planetary loss for systems initially containing two Jovian planets (m₁ = m₂ = 1M_J). The trends are roughly similar to the case with lower mass outer planets: planetary survival decreases with increasing migration rate. Turbulence leads to planetary loss in the regime of slow migration and no eccentricity damping, where the regime of slow migration corresponds to migration timescales longer than about τ_a = 3 × 10⁴ yr. And, as the migration rate increases, there is a shift from loss of planets through ejection to loss of planets through accretion onto the central star. However, for these systems with two Jovian planets, ejections, collisions, and accretion events are on a more equal footing. One clear difference from the case of low-mass outer planets is that planet–planet collisions are more common (compare the blue curves in Figures 9 and 10). The other significant difference is that the inner planet is more often lost during accretion events, rather than the outer planet (shown by the dotted curves in Figure 10).

Zoom In Zoom Out Reset image size

For solar systems with sufficiently large values of the eccentricity damping parameter K, most of the planets survive over the relatively short timescales considered in this paper. For example, for cases with K = 10, most systems remain intact and neither eject nor accrete a planet. However, as shown by the comparison of Figures 3 and 4, the fraction of systems that remain in mean motion resonance for K = 10 is only moderately increased over the values obtained for K = 1. The solar systems that are not in resonance will often eject or accrete planets on longer timescales, even in the absence of additional migration (e.g., Holman & Wiegert 1999; David et al. 2003). This issue should be addressed with additional, longer term numerical integrations, but is beyond the scope of the present work.

As a starting point, we consider a test particle of mass m₀ orbiting in the same plane as a larger planet with mass m_P, which is orbiting a star of mass M_*. The masses thus obey the ordering

$$\begin{equation} m_0 \ll m_P \ll M_\ast \,. \end{equation} \tag{ 11 }$$

The orbital elements of the test particle are as follows: λ is the mean longitude, M is the mean anomaly, a is the semimajor axis, ϖ is the longitude of pericenter, and e is the orbital eccentricity. The analogous variables for the planet have the same symbols but are denoted with the subscript "P" (see below). The Poincaré coordinates (MD99) can be written as

$$\begin{equation} \lambda = M + \varpi \qquad {\rm and} \qquad { {\widetilde{\gamma }} }= - \varpi \,, \end{equation} \tag{ 12 }$$

with momentum variables of the form

$$\begin{equation} \ell = (G M_\ast a)^{1/2} \quad {\rm and} \quad \Gamma = (G M_\ast a)^{1/2} [ 1 - (1 - e^2)^{1/2} ]. \end{equation} \tag{ 13 }$$

The Hamiltonian can be written in the form

$$\begin{equation} H = - {(G M_\ast)^2 \over 2 \ell ^2} - { {\cal R} }\,, \end{equation} \tag{ 14 }$$

where ${ {\cal R} }$ is the disturbing function due to the gravitational interaction between the test particle and the planet.

Specializing to the case of a 2:1 mean motion resonance where the planet is the inner body, we perform a canonical transformation using the generating function

$$\begin{equation} F_2 = I (2 \lambda - \lambda _P) \,, \end{equation} \tag{ 15 }$$

which leads to the new variables

$$\begin{equation} I = - \ell /2 \qquad {\rm and} \qquad \psi = \lambda _P - 2 \lambda \,. \end{equation} \tag{ 16 }$$

The new Hamiltonian for the unperturbed problem, without the disturbing function, has the form

$$\begin{equation} H_{0;\rm{new}} = - {(GM_\ast)^2 \over 8 I^2} - I n_P \,, \end{equation} \tag{ 17 }$$

where n_P is the mean motion of the planet.

Next, we express all quantities in dimensionless form and expand around the resonance. Here, distances are measured in units of the semimajor axis a, time is measured in units of (a²/GM_*)^1/2, and mass is measured in units of M_*. If we define

$$\begin{equation} \delta \equiv I - I_0 \qquad {\rm and} \qquad {1 \over 4 I_0^3} = n_P (t_0) \,, \end{equation} \tag{ 18 }$$

the new Hamiltonian now reads

$$\begin{equation} K_{0;\rm{new}} = {\mbox{constant}} - (n_P - 1) \delta - {3 \delta ^2 \over 8 I_0^4} \,. \end{equation} \tag{ 19 }$$

We must now include the relevant terms from the disturbing function, which provides an expansion in orders of eccentricity (of both the test mass and the planet). Here we keep only the leading order term (see MD99; Quillen 2006) and write the Hamiltonian (from Equation (19)) in the form

$$\begin{eqnarray} K (\delta, \psi, \Gamma, { {\widetilde{\gamma }} }) &=& - 6 \alpha ^2 \delta ^2 - (n_P - 1) \delta - 2 \mu f_2 \alpha ^{1/2} \Gamma \nonumber\\ &&+ A \Gamma ^{1/2} \cos (\psi - \varpi) \,, \end{eqnarray} \tag{ 20 }$$

where A is the expansion coefficient in the disturbing function and where we have used the fact that I₀ = α^−1/2/2 for these units and choice of resonance.

Following Quillen (2006), we perform another canonical transformation using the generating function

$$\begin{equation} F_2 = J_1 (\psi - \varpi) + J_2 \psi \,, \end{equation} \tag{ 21 }$$

which leads to the new variables

$$\begin{equation} J_1 + J_2 = \delta, \quad \phi = \psi - \varpi, \quad J_1 = \Gamma, \quad {\rm and} \quad \theta = \psi, \end{equation} \tag{ 22 }$$

and hence the new Hamiltonian

$$\begin{eqnarray} H &=& - 6 \alpha ^2 \left(\Gamma ^2 + J_2^2 \right) - [ 12 \alpha ^2 J_2 + (n_P - 1) + 2 \mu f_2 \alpha ^{1/2} ] \Gamma \nonumber\\ &&- (n_P - 1) J_2 + A \Gamma ^{1/2} \cos \phi \,. \end{eqnarray} \tag{ 23 }$$

Since J₂ is conserved and constant terms can be dropped, the Hamiltonian can be simplified to the form

$$\begin{equation} H = 6 \alpha ^2 \Gamma ^2 + [ 12 \alpha ^2 J_2 + (n_P - 1) + 2 \mu f_2 \alpha ^{1/2} ] \Gamma - A \Gamma ^{1/2} \cos \phi \,. \end{equation} \tag{ 24 }$$

Next, we rescale the momentum variable Γ according to the transformation

$$\begin{equation} \Gamma \rightarrow \left[ {6 \alpha ^2 \over A} \right]^{2/3} \Gamma \,, \end{equation} \tag{ 25 }$$

and rescale the time variable so that the Hamiltonian H is given by

$$\begin{equation} H = \Gamma ^2 + b \Gamma - \Gamma ^{1/2} \cos \phi \,. \end{equation} \tag{ 26 }$$

The parameter b is thus given by

$$\begin{equation} b = [ 12 \alpha ^2 J_2 + (n_P - 1) + 2 \mu f_2 \alpha ^{1/2} ] 6^{-1/3} (\alpha A)^{-2/3} \,. \end{equation} \tag{ 27 }$$

The first and third terms in square brackets are generally small compared to unity. The central term vanishes on resonance, by definition, but can be of order unity when the system is far from resonance. As a result, the parameter b provides a measure of how far the system resides from a resonant condition. For this paper, we let the parameter b evolve linearly with time so that the systems approach resonance (b = 0) at a well-defined rate.

Using the Hamiltonian with the form given by Equation (26), the equations of motion become

$$\begin{equation} {d \Gamma \over dt} = - \Gamma ^{1/2} \sin \phi \, \end{equation} \tag{ 28 }$$

and

$$\begin{equation} {d \phi \over dt} = 2 \Gamma + b - {1 \over 2 \Gamma ^{1/2}} \cos \phi \,. \end{equation} \tag{ 29 }$$

It is useful to define the reduced momentum variable p ≡ Γ^1/2 so that the equations of motion simplify to the forms

$$\begin{equation} 2 {dp \over dt} = - \sin \phi \, \end{equation} \tag{ 30 }$$

and

$$\begin{equation} {d\phi \over dt} = 2 p^2 + b - {1 \over 2p} \cos \phi \,. \end{equation} \tag{ 31 }$$

Although this ansatz simplifies the equations of motion, note that the variables (ϕ, p) are no longer canonical. We also note that this change of variables is convenient for calculating curves in phase space to analyze the dynamics (this exercise is carried out in the Appendix).

Using the model equations derived above, we can study the entry into mean motion resonance as a function of the normalized migration rate db/dt. Here, the initial conditions are given by the starting momentum Γ₀ and the starting value of the angle ϕ. We choose fixed values of the momentum variable Γ₀ and then study the probability of entering into resonance as a function of migration rate db/dt. Since these systems often display extreme sensitivity to their starting conditions, we must perform many realizations of the numerical integrations for each pair (Γ₀, db/dt), where each realization uses a different value of the starting angle ϕ. For the sake of definiteness, we start the systems with b = b₀ = 10 (well outside of resonance) and let the resonance parameter evolve according to the relation b(t) = b₀ − (db/dt)t. The systems thus pass through resonance at time t = b₀/|db/dt|.

One example integration is shown in Figure 11, which plots the quantity sin ϕ (top panel) and the momentum variable Γ (bottom panel) as a function of time for a system that becomes locked into mean motion resonance. In this case, the libration width of the system steadily decreases with time until it reaches a steady state near time t = 100 (in dimensionless units). In this case, db/dt = 0.1, so that t ∼ 100 corresponds to the time when the system passes through resonance (as expected). The momentum variable Γ stays small until the system enters resonance and then grows steadily (see also the discussion of Quillen 2006).

Zoom In Zoom Out Reset image size

As found previously (Quillen 2006), the probability of entering and surviving in resonance decreases with increasing migration rate. This trend is illustrated in Figure 12, which shows the probability of achieving a resonant state versus the migration rate db/dt. The three curves shown in the figure use different starting values of the momentum variable Γ₀ = 0.1, 1, and 3. Recall that Γ is related to the orbital eccentricity of the migrating planet (Equation (13)). Previous work shows that small starting momentum generally leads to resonance capture, whereas larger values generally do not (Quillen 2006); the value Γ₀ = 1 corresponds to the transition region. For each value of the rate db/dt, we have performed an ensemble of 1000 integrations, each with a different starting value of the angular variable ϕ. The probability of capture decreases with increasing db/dt, but the curves show a great deal of additional structure. The probability of achieving resonance decreases near db/dt = 1 and approaches zero for somewhat larger values db/dt ∼ 3 − 5.

Zoom In Zoom Out Reset image size

Another clear trend is that increasing the initial value of the momentum variable Γ₀ acts to decrease the probability of entering into resonance. In other words, larger eccentricities tend to compromise the chances of attaining resonance. This finding is consistent with the full numerical integrations of the previous section, where eccentricity damping was found to allow for more resonant states (see Figure 5).

The leading order trend illustrated by Figure 12 is that resonant capture is more difficult with fast migration. This result, obtained from the model equations of this section, is thus consistent with the results of the numerical simulations of Section 2 (see Figures 2–8). We can understand this effect through a simple analysis: in the limit of large db/dt = γ, which we consider to be a constant, the equation of motion for ϕ simplifies to the form

$$\begin{equation} {d \phi \over dt} = - \gamma t \quad \Rightarrow \quad \phi = - {1 \over 2} \gamma t^2 \,, \end{equation} \tag{ 32 }$$

where we have used the same sign convention as before. The momentum variable is then given by the remaining equation of motion, which can be written in the integral form

$$\begin{eqnarray} p - p_0 &=& {1 \over 2} \int \sin \left({1 \over 2} \gamma t^2 \right) dt = {1 \over (2 \gamma)^{1/2} } \int \sin u^2 du \nonumber\\ &=& {1 \over 2} \left({\pi \over \gamma } \right)^{1/2} S \left[ \left({\gamma \over \pi } \right)^{1/2} t \right] \,, \end{eqnarray} \tag{ 33 }$$

where S(z) is the Fresnel integral (e.g., Abramowitz & Stegun 1965). In the limit z → ∞, S(z) → 1/2, so the expression on the right-hand side of Equation (33)) approaches a constant value (π/16γ)^1/2. As a result, the momentum variable p approaches a constant, and hence does not grow, so the system does not enter resonance. For a given starting value of the momentum variable Γ₀, the critical value of the migration rate γ can be estimated by

$$\begin{equation} \gamma _c = {\pi \over 16 \Gamma _0} \qquad {\rm or} \qquad \gamma _c \, \Gamma _0 \approx 1/5 \,. \end{equation} \tag{ 34 }$$

For comparison, in the set of simulations shown in Figure 12 with Γ₀ = 1, the probability of survival in resonance P_b falls below unity when the migration rate becomes greater than γ = db/dt ≈ 0.2; P_b falls below 1/2 for γ>1 and goes to zero for larger values.

Another trend present in Figure 12 is that small variations in the migration rate can significantly change the probability of resonant capture, especially for larger starting values of the momentum variable. The curves shown in the Figure display a great deal of variation with db/dt; if the curves were plotted with finer resolution in db/dt, the plot would show even greater variation (and would not show resolved oscillations). This sensitivity to the migration rate can be illustrated further by plotting the time evolution of two nearly identical systems, as shown in Figure 13. In this case, two systems are started with Γ₀, the same angle ϕ₀, and two different migration rates db/dt = 0.300 (solid curve) and db/dt = 0.301 (dashed curve). The evolution of the two systems is nearly identical until about halfway through the total time interval, when the second system becomes locked into mean motion resonance (indicated by the growing values of Γ), whereas the first system continues to circulate with its momentum variable exhibiting a decreasing amplitude.

Zoom In Zoom Out Reset image size

This effect can be (roughly) understood as follows. Suppose we consider circulating solutions such that ϕ ≈ ωt. The equation of motion for the angle ϕ then implies that

$$\begin{equation} \omega \sim 2 p^2 + b - {1 \over 2p} \cos \phi \,. \end{equation} \tag{ 35 }$$

The corresponding solution for the reduced momentum variable p then becomes

$$\begin{equation} p(t) \sim A + {1 \over 2 \omega } \cos \omega t \,, \end{equation} \tag{ 36 }$$

where A is a constant. In this context, the parameter b starts at a positive value (outside resonance) and then decreases. The relation (35) indicates that ω must decrease with time, so that the amplitude of the oscillations of momentum increases with time as the frequency decreases. Near the point where ω → 0, however, the oscillation amplitudes are large and the frequency is small. The system must then match onto one of the possible solutions for late times when b is large and negative. One solution corresponds to ω → b (see relation (35)); in this case, Equation (36) indicates that the momentum variable will oscillate with increasing frequency and decreasing amplitude (shown by the solid curve in Figure 13). Although the momentum variable oscillates, the resonance angle circulates for this case. A second solution exists for sufficiently large p; in this case, the equation of motion (31) for the variable ϕ takes the approximate form

$$\begin{equation} {d \phi \over dt} \approx 2 p^2 + b \,. \end{equation} \tag{ 37 }$$

This equation can be combined with the momentum Equation (30) to obtain the result

$$\begin{equation} {d^2 \phi \over dt^2} + 2 p \sin \phi + {db\over dt} = 0 \,, \end{equation} \tag{ 38 }$$

which is a type of pendulum equation, and hence allows for librating solutions for the angle ϕ(t). This class of solution is depicted by the dashed curve in Figure 13.

In the limit of large p ≫ 1, we can ignore the cosine term in the denominator of Equation (A1) and find the approximate solution

$$\begin{equation} {4 \over 3} \left(p^3 - p_0^3 \right) + 2 b (p - p_0) = \cos \phi - 1 \,. \end{equation} \tag{ A4 }$$

This result can be rewritten in the form

$$\begin{eqnarray} - \sin ^2 (\phi /2) &=& (p - p_0) \left\lbrace {2 \over 3} \left[ p^2 + p p_0 + p_0^2 \right] + b \right\rbrace \nonumber\\ &\approx & (p - p_0) \left\lbrace 2 p_0^2 + b \right\rbrace \,. \end{eqnarray} \tag{ A5 }$$

Note that in the limit of large p, |p − p₀| ≪ p, and the two expressions in the above equation are the same to leading order in the parameter |p − p₀|/p.

In the limit of small p ≪ 1, we can ignore the p² term in Equation (A1) and find the solution

$$\begin{equation} p \cos \phi = b \left(p^2 - p_0^2 \right) + p_0 \,. \end{equation} \tag{ A6 }$$

For sufficiently small p, this expression reduces to the simpler form

$$\begin{equation} p \cos \phi = p_0 \,. \end{equation} \tag{ A7 }$$

The solution for the phase space curves, given by Equation (A3), can allow for multiple roots. We first note that the parameters b, E, and cos ϕ can all be both positive and negative. As a result, the number of roots for p will vary.

Case I: b > 0, E > 0. In this case, only one root for p exists for all values of the angle ϕ (or cos ϕ). For small values of the energy E, the solutions for p get small for negative values of cos ϕ.

Case II: b < 0, E > 0. For cos ϕ>0, only one solution for p exists. For cos ϕ < 0, however, there can be multiple roots provided that |b| is large enough and the energy E is small enough. The conditions required for multiple roots are that

$$\begin{equation} |b| > {8 \over 27} \cos ^2\phi \qquad {\rm and} \qquad 4 |b| E < \cos ^2 \phi \,. \end{equation} \tag{ A8 }$$

In the regime where the "extra" roots arise, p ≪ 1, and the solutions reduce to the approximate form

$$\begin{equation} p = {\cos \phi \pm [ \cos ^2 \phi - 4 |b| E ]^{1/2} \over 2 |b|} \,. \end{equation} \tag{ A9 }$$

Case III: b > 0, E < 0. There are no solutions for cos ϕ < 0. For the case cos ϕ>0, there are either two solutions for small E, or no solutions.

Case IV: b < 0, E < 0. For cos ϕ>0, there are two solutions for small E and no solutions if E is too large and negative. For the case cos ϕ < 0, the two solutions disappear if |b| is too small, where the critical value (for the limiting case cos ϕ = −1) is given by

$$\begin{equation} |b|_c = {3 \over 4^{1/3}} \approx 1.587 \,. \end{equation} \tag{ A10 }$$

For cases of interest, the parameter b becomes large and negative. For bound states, the energy also becomes large and negative. In this regime, the phase curves become almost independent of angle ϕ, with

$$\begin{equation} p^2 \approx {|b| \pm [ b^2 + 4 E ]^{1/2} \over 2} \,. \end{equation} \tag{ A11 }$$

Figure 14 shows one sample phase plot for the case where b = 0, which corresponds to systems that are passing through the resonant condition. For this case, as the energy variable E decreases from positive to negative values, the phase curves change their shape: for positive E, solutions for the momentum variable p exist for all values of the angle ϕ; for negative E, solutions for p exist for a limited range of angles. These isolated regions, which become narrower as the energy E grows more negative, correspond to oscillatory solutions in ϕ (such as that shown in Figure 11).

Zoom In Zoom Out Reset image size