Modeling Radial Velocity Data of Resonant Planets to Infer Migration Histories

Our first goal is to develop a method to evaluate whether a planetary system's RV observations can be explained by a pair of planets in an ACR configuration. In order to do so, we adopt a Bayesian model comparison framework (e.g., Gregory 2005). Given a set of observational data D and a model ${{ \mathcal M }}_{i}$ for how that data depend on a set of model parameters θ, Bayes theorem states

$$\begin{eqnarray*}&&p(\theta | D,{{ \mathcal M }}_{i})=\displaystyle \frac{p(D| \theta ,{{ \mathcal M }}_{i})p(\theta | {{ \mathcal M }}_{i})}{p(D| {{ \mathcal M }}_{i})}.\end{eqnarray*}$$

In standard Bayesian parlance, $p(\theta | D,{{ \mathcal M }}_{i})$ is referred to as the posterior, $p(D| \theta ,{{ \mathcal M }}_{i})$ as the likelihood, $p(\theta | {{ \mathcal M }}_{i})$ as the prior, and $p(D| { \mathcal M })$ as the evidence. Parameter inference methods such as MCMC compute or approximate the posterior distribution, $p(\theta | D,{{ \mathcal M }}_{i})$, which represents one's beliefs about the underlying system parameters, conditioned on the observation. In model comparison problems, one is instead interested in determining evidences, $p(D| {{ \mathcal M }}_{i})$, for two or more competing models. Two models ${{ \mathcal M }}_{i}$ and ${{ \mathcal M }}_{j}$ can be compared by computing their Bayes factor, $p(D| {{ \mathcal M }}_{i})/p(D| {{ \mathcal M }}_{j})$, the ratio of the two models' posterior probabilities in light of the data, under the assumption that the two theories are equally likely to be true a priori.

We sample posteriors and compute evidence for two classes of models for each planet pair's RV data. The first model class is a conventional N-body model that computes the RV signal by integrating the equations of motion governing the planetary system. The N-body model's parameters include the planets' masses m_i and orbital parameters at a reference epoch, chosen to be the median time of the RV observations. We parameterize the planets' orbits in terms of their orbital periods P_i, mean anomalies M_i, eccentricities e_i, and arguments of periapse ω_i. Our model's prior probabilities are uniform in the planets' orbital parameters and log uniform in the planets' masses. We adopt the usual strategy of sampling in variables $\sqrt{{e}_{i}}\cos {\omega }_{i}$ and $\sqrt{{e}_{i}}\sin {\omega }_{i}$ in our MCMC fitting to improve convergence behavior. The nested sampling algorithm that we use requires variables be normalizable, so, in addition to restricting ${e}_{i}\in [0,1)$ and restricting angular variables to lie between 0 and 2π, we restrict the range of our priors in the variables P_i and m_i: priors on P_i are uniform within ±25% of values determined by a maximum likelihood fit and priors on m_i are log uniform within a factor of 1/5 and 5 times the maximum likelihood fit values.⁵ In addition to the planets' orbital parameters, the model includes an RV zero point, γ_k, and instrumental jitter term, σ_k, for each set of RV observations. We adopt uniform priors for the zero points, γ_k, between the minimum and maximum values of each instrument's reported RV measurements. We adopt log-uniform priors for the jitter terms, σ_k, ranging between 10⁻³ and 30 times each instrument's median reported uncertainty. Parameters in the N-body model that lead to close encounters within 3 planet Hill radii or less during the integration are automatically rejected. For simplicity, we fix the host star mass to the reported best-fit value from the literature and assume the planets' orbits are coplanar with an inclination of 90°, i.e., the system is observed edge on. We will refer to this model as the "unrestricted" model.

The second model class presumes that the RV data are the product of a pair of planets residing in an ACR configuration. This model is similar to the RV model developed by Baluev (2008). Restricting planets to an ACR configuration reduces the degrees of freedom of the model. We parameterize an ACR model for a given j:j − k resonance in terms of the inner planet's RV semiamplitude, K₁, period P₁, time of transit ${T}_{c,1}$, and argument of periapse ω₁, along with m₂/m₁, the mass ratio of the outer planet relative to the inner, and a parameter, s, described below that sets the normalized AMD of the planet pair. Our ACR models also include the same instrumental parameters, γ_k and σ_k, used in the unrestricted model, and we adopt the same priors for these parameters in the ACR model as in the unrestricted model. The parameter $s\in [0,1]$ sets the AMD of a planet pair to be ${ \mathcal D }={s}^{2}\ {{ \mathcal D }}_{\max }$, where ${{ \mathcal D }}_{\max }$ is the AMD at the point where all the ACR tracks for different planet-mass ratios intersect at a global maximum of the averaged disturbing function (see Michtchenko et al. 2006).⁶ Because, at small to moderate eccentricity values, ACR configurations obey ${\beta }_{1}\sqrt{\alpha }{e}_{1}\propto {\beta }_{2}{e}_{2}$ (e.g., Deck & Batygin 2015), this definition of s implies that e_i ∝ s for small to moderate eccentricities. In our ACR model MCMC fits, we sample in the variables $\sqrt{s}\cos {\omega }_{1}$ and $\sqrt{s}\sin {\omega }_{1}$ rather than s and ω₁. Consequently, our priors are uniform in s and ω₁. We adopt uniform priors in P₁ and ${T}_{c,1}$. The orbital period of the outer planet's RV signal is set to ${P}_{2}=\tfrac{j}{j-k}{P}_{1}$, and the semiamplitude of the outer planet is given by

$$\begin{eqnarray*}&&{K}_{2}=\displaystyle \frac{{m}_{2}}{{m}_{1}}\,{\left(\displaystyle \frac{j-k}{j}\right)}^{1/3}\displaystyle \frac{\sqrt{1-{e}_{1}^{2}}}{\sqrt{1-{e}_{2}^{2}}}{K}_{1}.\end{eqnarray*}$$

The argument of periapse of the outer planet is set to ω₂ = ω₁ + π. The ACR model includes, in addition to the continuous parameters just described, a discrete degree of freedom that fixes the relative orbital phases of the planet pair. Specifically, with the resonant angle σ₁ fixed at its equilibrium value σ₁ = 0 (when k is odd) or σ₁ = π/k (when k is even), the orbital phase of the outer planet is related to the orbital phase of the inner planet by

$$\begin{eqnarray}&&{M}_{2}=\left(1-\displaystyle \frac{k}{j}\right){M}_{1}+\displaystyle \frac{1-j+k}{j}\pi +n\displaystyle \frac{2\pi }{j},\end{eqnarray} \tag{ 5 }$$

for $n=0,1,\ldots j-1$. We handle this discrete degree of freedom by running multiple nested sampling fits, each with n fixed to one of the j possible values. Then, assuming that any of the j possible values are equally likely a priori, the evidence of the ACR model is computed according to

$$\begin{eqnarray}&&p(D| {{ \mathcal M }}_{\mathrm{ACR}})=\displaystyle \frac{1}{j}\displaystyle \sum _{m=0}^{j-1}p(D| {{ \mathcal M }}_{\mathrm{ACR}},n=m),\end{eqnarray} \tag{ 6 }$$

where $p(D| {{ \mathcal M }}_{\mathrm{ACR}},n=m)$ is the evidence computed with n in Equation (5) fixed to m. Likelihoods for both models are computed in the standard way from a χ² value that includes jitter terms added in quadrature to the reported measurement uncertainties. We use the radvel (Fulton et al. 2018) code, modified to include our N-body and ACR models, to forward-model RVs, evaluate model likelihoods, and run MCMC simulations. We use the nested sampling code dynesty (Speagle 2020) to estimate the Bayesian evidence of each model. Specifically, we use the dynesty.NestedSampler class with default settings. For the ACR model, RV signals are synthesized as the sum of two planets' Keplerian RV signals rather than by means of N-body simulation. We discuss the validity of this approximation in Appendix B, where we demonstrate that it has negligible impact over the timescales spanned by the observations we fit.

In this section we fit and compare our N-body and ACR models to two systems hosting potentially resonant giant-planet pairs.

3.2.1. HD 45364

The first system we fit is HD 45364, a system of two planets in a 3:2 MMR orbiting a K0 star of mass ${M}_{* }=0.82\,{M}_{\odot }$ originally discovered by Correia et al. (2009). The best-fitting solution reported by Correia et al. (2009) displays large resonant angle libration amplitudes and secular eccentricity variations and is therefore not in an ACR configuration. A number of subsequent studies have further explored the dynamics and possible formation histories of the HD 45364 system. Rein et al. (2010) perform a suite of hydrodynamical simulations of planet–disk interactions in order to explore potential migration histories leading to the present-day orbital configuration of HD 45364 system. Their simulations produce systems with orbital parameters significantly different from the Correia et al. (2009) best-fit solution but with RV signals that fit the observations with approximately equal statistical significance. Correa-Otto et al. (2013) show that the dynamical configurations similar to those found by both Correia et al. (2009) and Rein et al. (2010) can be obtained through simple parameterized migration simulations with exponential migration and eccentricity-damping forces. As expected, the planets in the simulations of Correa-Otto et al. (2013) are driven into ACR configurations. While Correa-Otto et al. (2013) note that the ACR configurations' lack of libration amplitude is at odds with the orbital solutions reported by Correia et al. (2009) and Rein et al. (2010), they do not attempt to determine whether the RVs of HD 45364 can be adequately fit with planets in a zero libration amplitude configuration. Delisle et al. (2015) also analyze the HD 45364 system in order to constrain properties of the system's migration during formation, similar to Section 4.2 below. However, their constraints are derived under the assumption that the system possesses the nonzero libration amplitude measured by Correia et al. (2009), and they do not explore the possibility that the system may in fact reside in an ACR configuration.

Figure 2 summarizes our fit results for HD 45364. The ACR model is formally preferred by our nested sampling simulations, with a log Bayes factor of ∼7. The rightmost panel of Figure 2 shows that the posterior distribution of the N-body model planet eccentricities are concentrated quite close to the ACR values. We integrate 1500 initial conditions randomly selected from the unrestricted model MCMC posterior sample for 100 orbits of planet b in order to determine the resonant angle libration amplitudes and average period ratios plotted in the middle panels of Figure 2. Resonant angle libration amplitudes were recorded by determining the maximum deviations of the angles ${\theta }_{b}=3{\lambda }_{c}-2{\lambda }_{b}-{\varpi }_{b}$ and ${\theta }_{c}=3{\lambda }_{c}-2{\lambda }_{b}-{\varpi }_{c}$ from their respective equilibrium values of 0 and π over the course of the integration. Cumulative distributions of libration amplitudes are plotted in the center-right panel of Figure 2. While the majority of initial conditions show resonant libration for both angles with amplitudes $\lesssim \pi /2$, neither cumulative distribution reaches 1 due to a small fraction of simulations in which one or both resonant angles circulate. The moderate to large libration amplitudes inferred when fitting the unrestricted model may initially appear to be at odds with the preference for the zero libration amplitude ACR model inferred from the Bayesian evidence calculations. This apparent contradiction suggests that the inference of large libration amplitudes in the unrestricted model is principally driven by the larger phase-space volume occupied by these configurations and not by any significant improvement in fit derived from large libration amplitude solutions. This example serves to illustrate that the influence of priors must be carefully considered when deriving conclusions about the dynamics of planetary systems from RV fits.

Zoom In Zoom Out Reset image size

3.2.2. HD 33844

The second system we fit, HD 33844, consists of two giant planets orbiting in or near a 5:3 MMR and was originally reported by Wittenmyer et al. (2016). We adopt the stellar mass of ${M}_{* }=1.84\,{M}_{\odot }$ for HD 33844 based on the best-fit value reported by Stassun et al. (2017). Combining RV fits with longer-term N-body simulations to rule out RV solutions leading to dynamical instabilities on timescales less than 100 Myr, Wittenmyer et al. (2016) find that many of the stable orbital configurations in the neighborhood of their best-fit solution reside in the 5:3 MMR.⁷

Figure 3 summarize our fit results for HD 33844. Our inferred planet masses and orbital parameters are in good agreement with values reported by Wittenmyer et al. (2016). The unrestricted and ACR models provide similar quality fits to the data and, due to its lower dimensionality, the ACR model is formally preferred with a log Bayes factor of ∼6 based on the results of our nested sampling simulations.

Zoom In Zoom Out Reset image size

In order to more clearly interpret the resonant state of RV solutions inferred with the unrestricted model, we turn to the analytic resonance model of Hadden (2019). While the full dynamics of a 5:3 MMR are captured by a two-degree-of-freedom system with two resonant angles, σ₁ and σ₂ defined in Section 2, Hadden (2019) shows that, to an excellent approximation, the dynamics of the MMR only depend on a single resonant angle ${\theta }_{\mathrm{res}}=5{\lambda }_{c}-3{\lambda }_{b}-2z$, where $z\,\approx \arg ({e}_{c}{e}^{i{\varpi }_{c}}-{e}_{b}{e}^{i{\varpi }_{b}})$. To determine the resonant behavior of the RV solutions derived from the unrestricted model, we again integrate 1500 initial conditions randomly selected from the MCMC posterior sample for 100 orbits of planet b and record angle libration amplitudes as the maximum deviation of the angle θ_res from its equilibrium value, π. (ACR solutions have periodic variations in θ_res on the synodic timescale about the equilibrium of the averaged system with an amplitude of approximately ∼0.1π.) Approximately 50% of these initial conditions show librating behavior, and the cumulative distribution of their libration amplitudes are recorded in the middle-right panel of Figure 3. The rightmost panel shows that the N-body model prefers a larger eccentricity for planet b when the system is not restricted to lie inside the ACR.

So far, we have not considered the long-term stability of planets' orbital configurations in our fits. The high planet masses and close spacings of the HD 45364 and HD 33844 systems place them near the boundary of dynamical instability. This dynamical fragility can be leveraged to derive tighter constraints on the planets' masses and orbital properties by rejecting any orbital configurations that rapidly go unstable (rejection sampling based on long-term stability is common practice in the literature: e.g., Wittenmyer et al. 2016). Because we are interested in evaluating whether systems' RVs can be explained by planets in an ACR configuration rather than deriving precise orbital constraints, we will not apply this rejection sampling procedure here. Instead, we explore the dynamical stability in the vicinity of ACR configurations by means of the dynamical maps shown in Figure 4. The maps in Figure 4 indicate dynamical stability results for a grid of short N-body integrations. In each integration, planets are initialized with masses equal to the median planet masses derived from our MCMC fits and orbital elements corresponding to an ACR equilibrium, except that we vary the planets' period ratio away from the equilibrium value. We explore a range of planet eccentricities by running integrations over a grid of normalized AMD values up to the value ${{ \mathcal D }}_{\max }$ described in Section 3.1. Each integration is initialized with λ₂ = λ₁, and all osculating orbital elements are set to their ACR values, except for the planets' period ratios, which we vary along the x-axes of the maps. The planets' osculating elements in the ACR configuration are computed using the transformation from mean to osculating elements described in Appendix A.2. The simulations are run using the WHFast integrator (Rein & Tamayo 2015) with time steps set to 1/30 of the periapse passage timescale of the inner planet, ${T}_{\mathrm{peri}}=2\pi /{\dot{f}}_{\max }$, where ${\dot{f}}_{\max }$ is the time derivative of the true anomaly at pericenter (Wisdom 2015). Regular and chaotic trajectories are identified based on the MEGNO chaos indicator (Cincotta et al. 2003).

Zoom In Zoom Out Reset image size

From Figure 4, we see that, in the case of HD 45364, the 3:2 ACR provides a large island of stability in a phase space that is otherwise unstable except for a small pocket of stability just wide of the resonance at low eccentricity. This is in agreement with the original analysis of Correia et al. (2009). The structure of the 5:3 ACR in the HD 33844 system, on the other hand, is somewhat more complicated. The ACR shows two separate islands of stability separated by a thin chaotic region near ${ \mathcal D }\approx 0.1$ along with a small pocket of stability wide of the resonance at low eccentricity.

Irrespective of whether or not the HD 45364 and HD 33844 systems are truly in ACR configurations, we can conclude that they most likely reside in resonance as the local phase space outside of these resonances is almost entirely chaotic and unstable. The stability maps shown in Figure 4 highlight how knowledge of the ACR structure allows much more efficient identification of a stable parameter space than a brute-force grid search of the high-dimensional parameter space of planets' orbital elements when fitting RVs. It is well known that MMRs can provide phase protection from close encounters that would otherwise destabilize strongly interacting planets. Because ACR configurations are elliptic equilibria of the averaged dynamical system describing the planets' resonant motion, they should be the most dynamically stable resonant configurations possible: the degree of nonlinearity in the neighborhood of an elliptic equilibrium, and thus the potential for chaos and dynamical instability, increases with distance from the equilibrium (Giorgilli et al. 1989). Therefore, vicinities of ACR configurations provide good starting points for any search for dynamically stable orbital configurations when modeling strongly interacting planet pairs. In their analysis of the HD 82943 system, Baluev & Beaugé (2014) previously noted the utility of enforcing ACR configurations to ensure dynamical stability when fitting RV data. Conversely, if no stable ACR configuration exists for a given set of planet masses, then stable orbital configurations likely do not exist in the nearby phase space.

In this section we explore how the present-day dynamical configurations can be used to glean information about planet–disk interactions that presumably established these configurations during the planet formation process. Because theoretical understanding of planet migration is incomplete, we rely on our simple parameterized model of eccentricity damping and migration in Equation (2). While the migration and eccentricity-damping rates experienced by planets may vary over the course of the formation process in ways that are not captured by our simple parameterized damping model, the model is still useful for capturing qualitative aspects of planet–disk interactions during the formation process. The quantitative constraints derived with our simple migration model provide useful metrics against which predictions of more rigorous theoretical and numerical studies of planet migration can be compared.

If it is assumed that planet pairs were subject to migration and eccentricity-damping forces for a time period sufficiently long to reach dynamical equilibrium, a given ACR equilibrium configuration corresponds directly to a ratio of migration timescale to eccentricity-damping timescale at the time the equilibrium was established. This is illustrated in Figure 5. The figure shows the posterior distributions of each planet pair's eccentricities from our MCMC fits with ACR models. In different colored lines, we plot the loci of equilibria along these ACR tracks for different values of τ_α/τ_e, where we define ${\tau }_{e}^{-1}\equiv {\tau }_{e,1}^{-1}+{\tau }_{e,2}^{-1}$ and assume that ${m}_{2}{\tau }_{e,2}={m}_{1}{\tau }_{e,2}$ so that the strength of eccentricity damping experienced by a planet is proportional to its mass. This simple prescription approximately matches the expected scaling of eccentricity damping experienced due to type I migration (Kley & Nelson 2012) or dynamical friction from interactions with planetesimals (Ida 1990). For each system, Figure 5 also shows the distribution of τ_α/τ_e values inferred by mapping our posterior samples of mass ratios and eccentricities to values of τ_α/τ_e. For HD 45364 we infer ${\tau }_{\alpha }/{\tau }_{e}={8}_{-3}^{+8}$, and for HD 33844 we find ${\tau }_{\alpha }/{\tau }_{e}={20}_{-12}^{+65}$.

Zoom In Zoom Out Reset image size

Of course, the planet pairs could have failed to reach an equilibrium between migration and eccentricity-damping forces and instead have been stranded at their current resonant configurations if the gas disk dissipated before any equilibrium was reached. In this case, the posterior distributions of τ_α/τ_e derived from our ACR model fits would serve as lower limits on the timescale of eccentricity damping relative to the migration rate: had the eccentricity damping been stronger, the systems would have reached equilibria at lower eccentricities while ${ \mathcal D }$ increased monotonically after capture into resonance. Far from an equilibrium between migration and eccentricity-damping forces, the eccentricity damping will play a negligible role in the evolution of ${ \mathcal D }$ and

$$\begin{eqnarray}&&\dot{{ \mathcal D }}\approx \displaystyle \frac{{\beta }_{1}{\beta }_{2}\sqrt{\alpha }}{{\beta }_{1}\sqrt{\alpha }(1+s)+{\beta }_{2}s}\displaystyle \frac{1}{2{\tau }_{\alpha }}\,.\end{eqnarray} \tag{ 7 }$$

Inserting the appropriate parameters in Equation (7), HD 45364 would reach the median posterior value of ${ \mathcal D }$ increasing from 0 in a time of ∼0.15τ_α. Similarly, HD 33844 would reach its median value of ${ \mathcal D }$ in ∼0.03τ_α. The times the two systems would have required to reach their current MMRs by large-scale migration from a period ratio of approximately ∼2:1 are ∼0.2τ_α for HD 45364 and ∼0.1τ_α for HD 33844. This scenario therefore requires that planets' migration timescales be comparable to or longer than the lifetime of the protoplanetary disk while migration rates predicted by conventional type II migration theory are usually at least one or two orders of magnitude shorter than disk lifetimes (Kley & Nelson 2012).

One of the puzzling features of the HD 45364 and HD 33844 systems is that they both lie interior to the 2:1 MMR. In standard core accretion theories of planet formation, massive planets like those in these systems are expected to form farther apart than the 2:1 MMR (Armitage 2007 and references therein). If the resonances in these systems were established through convergent migration, then, during formation, the planets should have first encountered and been captured into the 2:1 MMR. One possibility is that their migration was sufficiently rapid to avoid capture in the 2:1 MMR. This is the scenario considered, for example, in the simulations of Rein et al. (2010) to explain the orbital configuration of HD 45364. Another possibility is that capture in the 2:1 MMR was rendered temporary in both systems because the ACR configuration was an unstable spiral in the presence of dissipative forces. This mechanism presents a compelling alternative explanation because, as Rein et al. (2010) note, the migration rates required to avoid capture in the 2:1 MMR are significantly faster than the rates of migration predicted by conventional type II migration theory. We explore this possibility further in Figures 5 and 6.

Zoom In Zoom Out Reset image size

Figure 5 illustrates the stability of ACR equilibria as a function of τ_α/τ_e for the current resonance occupied by each system as well as the 2:1 MMR (plus the 5:3 MMR for HD 45364 system), which the planets would have traversed on their way to their current MMR configuration. The curves are computed as follows: first, we determine the equilibrium configuration for a given τ_α/τ_e using the resonance equations of motion, taking p = 1 in Equation (2) describing the dissipative forces experienced by the planets. Next, we compute the Jacobian of the equations of motion at this equilibrium. The stability is then determined based on the real parts of the eigenvalues of the Jacobian as described in Section 2. The plot shows the maximum real eigenvalue component in units of ${\tau }_{\alpha }^{-1}$, which we have taken to be ${\tau }_{\alpha }=3\times {10}^{5}{P}_{2}$. (We confirm that, for sufficiently large τ_α, growth rates simply scale with ${\tau }_{\alpha }^{-1}$.) For HD 45364, the 2:1 ACR equilibrium is unstable for $6\lt {\tau }_{\alpha }/{\tau }_{e}\lt 211$. This range encompasses 65% of the posterior points plotted in the top right panel of Figure 5. For HD 33844, instability occurs at the 2:1 MMR for $3.6\lt {\tau }_{\alpha }/{\tau }_{e}\lt 67$, encompassing 81% of the posterior.

In Figure 6, we perform similar calculations, but we now relax our assumption of ${\tau }_{e,1}/{\tau }_{e,2}={m}_{2}/{m}_{1}$ and instead explore a range of relative eccentricity-damping strengths. We examine credible regions of τ_a/τ_e inferred from our ACR model fit posterior as a function of ${\tau }_{e,1}/{\tau }_{e,2}$. The figure also shows the regions for which different combinations of ${\tau }_{e,1}/{\tau }_{e,2}$ and ${\tau }_{\alpha }/{\tau }_{e}$ drive the system to an ACR configuration that is made unstable by dissipative forces.

While the elliptic point of the ACR configuration may be rendered unstable by dissipative forces, the ultimate outcome of the instability can lead to either passage through the resonance or a limit cycle within the resonance. In order for the planets of the HD 45364 and HD 33844 systems to have escaped permanent capture in a 2:1 MMR, instabilities must have led to escape rather than a limit cycle. Using an integrable approximation for the dynamics of a first-order MMR, Deck & Batygin (2015) derive an approximate criterion to determining whether instability is expected to result in an escape or a limit cycle. Their criterion can be stated in terms of a critical value of ${ \mathcal D }={{ \mathcal D }}_{\mathrm{bif}}$, where the Hamiltonian $\bar{H}({\boldsymbol{z}},{ \mathcal D })$ exhibits a bifurcation and goes from having one elliptic fixed point to having two elliptic fixed points and one unstable fixed point. If ${ \mathcal D }\gt {{ \mathcal D }}_{\mathrm{bif}}$, then the unstable librations grow until the system eventually escapes resonance while, for ${ \mathcal D }\lt {{ \mathcal D }}_{\mathrm{bif}}$, the librations saturate at a finite value when the system reaches a limit cycle. To compute ${{ \mathcal D }}_{\mathrm{bif}}$ for the 2:1 MMR, we use an integrable approximation for the resonant dynamics similar to the one employed by Deck & Batygin (2015) implemented in the CELMECH⁸ package. Damping parameter combinations expected to lead to a limit cycle in the 2:1 MMR, where the equilibrium configuration is unstable but ${ \mathcal D }\lt {{ \mathcal D }}_{\mathrm{bif}}$, are indicated by pink shading in Figure 6.

Figures 7 and 8 show example encounters with the 2:1 MMR for HD 45364 and HD 33844, respectively, at fixed ${\tau }_{\alpha }/{\tau }_{e}$ and for a range of ${\tau }_{e,1}/{\tau }_{e,2}$ values that are indicated by the colored points in Figure 6. The simulations in Figures 7 and 8 use the symplectic WHFast integrator, instead of the IAS15 integrator, with time steps set to 1/100 of the initial orbital period of the inner planet. We choose the planets' individual migration timescales such that ${\tau }_{\alpha }=3\times {10}^{5}{P}_{2}$ and $\tfrac{{m}_{1}}{{a}_{\mathrm{1,0}}{\tau }_{m,1}}+\tfrac{{m}_{2}}{{a}_{\mathrm{2,0}}{\tau }_{m,2}}=0$. The latter condition ensures that the time derivative of the systems' energy is approximately zero and minimizes any inward drift of the system to ensure that dynamics remain well resolved by the simulations' fixed time steps. Trajectories are colored according to the predicted stability of their 2:1 ACR equilibrium: pink trajectories are stable while green trajectories are unstable. In both Figures 7 and 8, the planets are initialized with a period ratio of 2.2 on circular orbits and migrate with ${\tau }_{\alpha }=3\times {10}^{5}{P}_{2}$. For HD 45364 shown in Figure 7, we set ${\tau }_{e}=1/({\tau }_{e,1}^{-1}+{\tau }_{e,2}^{-1})=8.13{\tau }_{\alpha }$ while for HD 33844 shown in Figure 8, we set ${\tau }_{e}=1/({\tau }_{e,1}^{-1}+{\tau }_{e,2}^{-1})=31{\tau }_{\alpha }$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For HD 45364, shown in Figure 7, the unstable simulations grow in libration amplitude and eventually escape from the resonance. The unstable simulations for HD 33844, shown in Figure 8, exhibit more complicated behavior. Upon reaching the unstable ACR equilibrium, they initially grow in libration amplitude. Three of the simulations go on to escape from the resonance, just as in the unstable HD 45364 simulations. The other unstable simulations, however, reach a limit cycle. This behavior is consistent with the predictions of Figure 6, where many of the simulations lie in the pink shaded region of parameter space, leading to an unstable 2:1 MMR but with ${ \mathcal D }\lt {{ \mathcal D }}_{\mathrm{bif}}$.

For HD 33844, Figure 6 shows that there is a significant region of parameter space which is both consistent with our ACR model fit and for which the 2:1 ACR configuration is unstable but the 5:3 ACR configuration remains stable. HD 33844 b and c could, therefore, have readily evaded permanent capture in the 2:1 MMR and continued migrating toward the 5:3 MMR, where they were ultimately permanently captured. There is also a region of parameter space for HD 45364 that is consistent with our ACR model fit and for which the 2:1 ACR configuration is unstable but the 3:2 ACR configuration remains stable. HD 45364 b and c could, therefore, also have escaped the 2:1 MMR and continued migrating toward the 3:2 MMR. However, in order to migrate into the 3:2 MMR, the system would also have to have traversed the second-order 5:3 MMR. There is a narrow region of parameter space in Figure 6 that leads to unstable ACR equilibria in both the 2:1 and 5:3 MMRs but a stable 3:2 ACR equilibrium. However, such parameters are not generic, and a scenario in which the planet pair escaped both the 2:1 and 5:3 MMRs through instability but reached a stable equilibrium in the 3:2 MMR appears to suffer from a fine-tuning problem. Alternatively, HD 45364 b and c could have reached the 3:2 MMR if their migration rate was sufficiently fast to avoid capture into the 5:3 MMR. This would require a migration rate sufficiently fast to avoid the 5:3 MMR but not too fast so as to avoid the 3:2 MMR. This is possible provided ${\left(\displaystyle \frac{{m}_{1}+{m}_{2}}{{M}_{* }}\right)}^{-2}\lesssim {\tau }_{\alpha }/P\lesssim {\left(\displaystyle \frac{{m}_{1}+{m}_{2}}{{M}_{* }}\right)}^{-4/3}$ (Xu & Lai 2017).

How do our inferences about the migration histories of HD 45364 and HD 33844 compare to theoretical predictions? For planets experiencing type I migration, the relative strength of eccentricity and migration rates, τ_α/τ_e, are expected to be ≳100 (e.g., Baruteau et al. 2014), significantly larger than the values we infer in Figures 5 and 6. This is especially true for HD 45364, where our constraints are the most precise. We note that Delisle et al. (2015) examined HD 45364 and found that τ_e,1/τ_e,2 ≳ 6 and τ_e/τ_α ∼ 10 based on the reported orbital solution of Correia et al. (2009), in agreement with our results. However, the planets in both systems are sufficiently massive that they are expected to have open gaps in their protoplanetary disk, so the predictions of type I migration are probably not applicable (e.g., Kley & Nelson 2012).

Hydrodynamical simulations of the migration of massive, resonant planet pairs in a gas disk have been studied by a number of authors (Bryden et al. 2000; Kley 2000; Snellgrove et al. 2001; Papaloizou 2003; Kley et al. 2004; Sándor et al. 2007; Crida et al. 2008; Rein et al. 2010; André & Papaloizou 2016). Kley et al. (2004) infer an effective ${\tau }_{\alpha }/{\tau }_{e}$ of approximately ∼1–10 from their simulations of two giant planets in a disk by comparing the results of their hydrodynamic simulations with simple damped N-body simulations. This result is in good agreement with the values inferred from our ACR fits for both systems. Crida et al. (2008) and Morbidelli & Crida (2007) note the importance of an inner disk to provide eccentricity damping for the innermost planet. Figure 6 shows that both systems would have required the eccentricity damping experienced by the inner planet to be comparable in magnitude to the damping experienced by the outer planet if they were to have escaped the 2:1 MMR due to instability and yet reached a stable equilibrium in their present-day resonances.

While we have sketched plausible migration histories for the HD 45364 and HD 33844 systems that can explain how they might evade capture in a 2:1 MMR and reach their current configurations, we cannot rule out other scenarios. For example, the planets may have had significantly different masses at the time they encountered the 2:1 MMR and only have grown to their current masses after becoming trapped in the resonances they currently occupy. Another possibility is that the disk conditions, and thus migration and eccentricity-damping rates, varied significantly between the planets' encounters with the 2:1 MMR and their subsequent encounters with more closely spaced resonances. Additionally, our treatment of eccentricity damping and migration is likely oversimplified. In Equation (2), we parameterized the dependence of the planets' semimajor axis evolution under disk forces by including the parameter "p," which describes the degree to which eccentricity-damping forces conserve angular momentum (p = 1 corresponds to eccentricity damping at constant angular momentum). As Xu et al. (2018) note, for type I migration, this treatment is only valid when eccentricities are less than the local disk aspect ratio, $e\ll H/r$, where H is the disk scale height and r the orbital radius. Xu et al. (2018) find that using a more realistic treatment of migration and eccentricity-damping forces in the type I regime generally leads to larger equilibrium eccentricities and more stable capture when compared to the approximate treatment we employ in Equation (2). For giant planets that open a common gap (the regime likely occupied by the planets we have considered here), the dependence of migration and eccentricity-damping forces on planet eccentricities is less well studied, though hydrodynamical simulations by Crida et al. (2008) suggest a potentially complicated dependence. Ultimately, a comparative study of a larger sample of resonant systems both in and interior to the 2:1 MMR, in conjunction with theoretical investigations of planet–disk interactions in systems of resonant gap-opening planet pairs, could shed further light on why some systems become captured in the 2:1 MMR while others manage to be captured into closer resonances.

In Section 3.2, we showed that the RV signals of both HD 45364 and HD 33844 were consistent with the hypothesis that the systems reside in ACR configurations. In each case, our nested sampling simulations yielded model evidences that prefer ACR configurations for both systems. However, in both cases, the preference for the ACR models is somewhat marginal. Here we briefly explore the prospects for follow-up observations to further bolster the evidence for ACR configurations in these systems or detect deviations from the predictions of the ACR models. In Figure 9, we show the predicted RV signals of HD 45364 and HD 33844 from 2020 April 15 to 2025 April 15. The figures compare the predicted RV signals projected using both the unrestricted and ACR models. In order to identify the most opportune times for follow-up observations, we seek to identify times where the models' posterior predictive distributions of RV values are most distinct. This is achieved by comparing the distribution of the predicted RV of a system at a given time using a two-sample Kolmogorov–Smirnov (KS) test (Kolmogorov 1933; Smirnov 1948). In the bottom panels of Figure 9, we plot the negative logarithm of the p values returned by these KS tests as a function of time. Small p values indicate that the distributions predicted by the two are highly unlikely to be derived from the same underlying distribution and provide a heuristic measure of the future epochs at which additional observations could most effectively distinguish between the two models.

Zoom In Zoom Out Reset image size

In this appendix, we derive the equations of motion governing the dynamics of a j:j − k resonance between a pair of coplanar planets of mass m₁ and m₂ orbiting a star of mass m₀. We begin by considering the conservative evolution of the system in Jacobi coordinates using modified canonical Delauney variables (e.g., Morbidelli 2002). The canonical momenta are given by

$$\begin{eqnarray*}\begin{array}{rcl}{{\rm{\Lambda }}}_{i} & = & {\tilde{m}}_{i}\sqrt{G{\tilde{M}}_{i}{a}_{i}}\\ {{\rm{\Gamma }}}_{i} & = & {\tilde{m}}_{i}\sqrt{G{\tilde{M}}_{i}{a}_{i}}(1-\sqrt{1-{e}_{i}^{2}}),\end{array}\end{eqnarray*}$$

where G is the gravitational constant, ${\tilde{m}}_{i}=\tfrac{{\eta }_{i-1}}{{\eta }_{i}}{m}_{i}$, ${\tilde{M}}_{i}\,=\tfrac{{\eta }_{i}}{{\eta }_{i-1}}{m}_{0}$, and ${\eta }_{i}={\sum }_{j=0}^{i}{m}_{i}$. The angle variables conjugate to Λ_i and Γ_i are, respectively, the planets' mean longitudes, λ_i, and ${\gamma }_{i}=-{\varpi }_{i}$, where ϖ_i are the planets' longitudes of periapse. The Hamiltonian is given by

$$\begin{eqnarray}&&\begin{array}{l}H({\lambda }_{i},{\gamma }_{i};{{\rm{\Lambda }}}_{i},{{\rm{\Gamma }}}_{i})=-\displaystyle \sum _{i=1}^{2}\displaystyle \frac{{G}^{2}{\tilde{M}}_{i}^{2}{\tilde{m}}_{i}^{3}}{2{{\rm{\Lambda }}}_{i}^{2}}\\ \quad -{{Gm}}_{1}{m}_{2}\,\left(\displaystyle \frac{1}{| {r}_{1}-{r}_{2}| }-\displaystyle \frac{{r}_{1}\cdot {r}_{2}}{| {r}_{2}{| }^{3}}\right),\end{array}\end{eqnarray} \tag{ A1 }$$

where the position vectors r_i are functions of the canonical variables, and we omit terms that are second order and higher in the planet−star mass ratios.

To derive equations of motion governing the dynamics of a j:j − k resonance, we follow Michtchenko et al. (2006) and transform to new canonical angle variables,

$$\begin{eqnarray}&&\left(\begin{array}{l}{\sigma }_{1}\\ {\sigma }_{2}\\ Q\\ l\end{array}\right)=\left(\begin{array}{llll}-s & 1+s & 1 & 0\\ -s & 1+s & 0 & 1\\ -1/k & 1/k & 0 & 0\\ \displaystyle \frac{1}{2} & \displaystyle \frac{1}{2} & 0 & 0\end{array}\right)\cdot \left(\begin{array}{l}{\lambda }_{1}\\ {\lambda }_{2}\\ {\gamma }_{1}\\ {\gamma }_{2}\end{array}\right),\end{eqnarray} \tag{ A2 }$$

where s = (j − k)/k, along with new conjugate momenta defined implicitly in terms of the old momentum variables as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{i} & = & {I}_{i}\\ {{\rm{\Lambda }}}_{1} & = & \displaystyle \frac{L}{2}-P/k-s({I}_{1}+{I}_{2})\\ {{\rm{\Lambda }}}_{2} & = & \displaystyle \frac{L}{2}+P/k+(1+s)({I}_{1}+{I}_{2}).\end{array}\end{eqnarray} \tag{ A3 }$$

The full Hamiltonian is independent of the new angle $l=({\lambda }_{1}+{\lambda }_{2})/2$, and its corresponding canonical momentum,

$$\begin{eqnarray*}&&L={{\rm{\Lambda }}}_{1}-{{\rm{\Gamma }}}_{1}+{{\rm{\Lambda }}}_{2}-{{\rm{\Gamma }}}_{2}\,,\end{eqnarray*}$$

is the total angular momentum, which is conserved due to the rotational symmetry of the system.

To derive equations of motion governing the resonant dynamics, we perform a near-identity canonical transformation $({\sigma }_{i},{I}_{i},Q,P)\to ({\bar{\sigma }}_{i},{\bar{I}}_{i},\bar{Q},\bar{P})$ such that, to leading order in the planet−star mass ratio, the new Hamiltonian is given by

$$\begin{eqnarray}&&\bar{H}({\bar{\sigma }}_{i},{\bar{I}}_{i};L,\bar{P})=\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{HdQ},\end{eqnarray} \tag{ A4 }$$

i.e., the new Hamiltonian is the old Hamiltonian averaged over the "fast" variable Q. We construct this transformation explicitly below in Appendix A.2 so that we may transform the osculating elements of a pair of planets to the variables of our Hamiltonian model and vice versa.

The averaging procedure yields

$$\begin{eqnarray*}&&\bar{P}/k=\displaystyle \frac{{\bar{{\rm{\Lambda }}}}_{2}-{\bar{{\rm{\Lambda }}}}_{1}}{2}-\left(s+\displaystyle \frac{1}{2}\right)({\bar{{\rm{\Gamma }}}}_{1}+{\bar{{\rm{\Gamma }}}}_{2})\end{eqnarray*}$$

as a constant of motion. In order to gain better physical intuition for the meaning of the conserved quantity, $\bar{P}$, it is instructive to express its value in terms of the value of the planets' eccentricities and distance to nominal resonance, $\tfrac{j-k}{j}\tfrac{{P}_{2}}{{P}_{1}}-1\equiv {\rm{\Delta }}$. This can be accomplished as follows: first, define the reference semimajor axes ${a}_{\mathrm{2,0}}$ and ${a}_{\mathrm{1,0}}\,={\left(\displaystyle \frac{{\tilde{M}}_{1}}{{\tilde{M}}_{2}}\right)}^{1/3}{\left(\displaystyle \frac{j-k}{j}\right)}^{2/3}{a}_{\mathrm{2,0}}$ corresponding to the nominal location of resonance, such that $L={{\rm{\Lambda }}}_{\mathrm{1,0}}+{{\rm{\Lambda }}}_{\mathrm{2,0}}$, where ${{\rm{\Lambda }}}_{i,0}=\tilde{{m}_{i}}\sqrt{G{\tilde{M}}_{i}{a}_{i,0}}$ for i = 1, 2. Conservation of $\bar{P}$ and L implies that ${\bar{{\rm{\Lambda }}}}_{2}-{{\rm{\Lambda }}}_{\mathrm{2,0}}=-\tfrac{s+1}{s}({\bar{{\rm{\Lambda }}}}_{1}-{{\rm{\Lambda }}}_{\mathrm{1,0}})$. Defining $\delta {{\rm{\Lambda }}}_{i}\,={{\rm{\Lambda }}}_{i}-{{\rm{\Lambda }}}_{i,0}$, we can write

$$\begin{eqnarray*}&&{\rm{\Delta }}\approx 3\left(\displaystyle \frac{\delta {{\rm{\Lambda }}}_{2}}{{{\rm{\Lambda }}}_{\mathrm{2,0}}}-\displaystyle \frac{\delta {{\rm{\Lambda }}}_{1}}{{{\rm{\Lambda }}}_{\mathrm{1,0}}}\right)=3\delta {{\rm{\Lambda }}}_{2}\left(\displaystyle \frac{(s+1){{\rm{\Lambda }}}_{\mathrm{1,0}}+s{{\rm{\Lambda }}}_{\mathrm{2,0}}}{(s+1){{\rm{\Lambda }}}_{\mathrm{2,0}}{{\rm{\Lambda }}}_{\mathrm{1,0}}}\right).\end{eqnarray*}$$

We define the "normalized AMD" as

$$\begin{eqnarray}&&{ \mathcal D }=\displaystyle \frac{\tfrac{1}{2}({{\rm{\Lambda }}}_{\mathrm{2,0}}-{{\rm{\Lambda }}}_{\mathrm{1,0}})-\bar{P}/k}{(s+1/2)({\tilde{m}}_{1}+{\tilde{m}}_{2})\sqrt{{{Gm}}_{0}{a}_{\mathrm{2,0}}}}.\end{eqnarray} \tag{ A5 }$$

Rewriting Equation (A5) in terms of orbital elements gives

$$\begin{eqnarray*}\begin{array}{rcl}{ \mathcal D } & \approx & {\beta }_{1}\sqrt{{\alpha }_{0}}\left(1-\sqrt{1-{e}_{1}^{2}}\right)\\ & & +{\beta }_{2}\left(1-\sqrt{1-{e}_{2}^{2}}\right)\\ & & -\displaystyle \frac{k{\beta }_{2}{\beta }_{1}\sqrt{{\alpha }_{0}}{\rm{\Delta }}}{3\left({\beta }_{1}\sqrt{{\alpha }_{0}}(s+1)+{\beta }_{2}s\right)},\end{array}\end{eqnarray*}$$

where ${\beta }_{i}=\tfrac{{\tilde{m}}_{i}}{{\tilde{m}}_{1}+{\tilde{m}}_{2}}\sqrt{\tfrac{{\tilde{M}}_{i}}{{m}_{0}}}$ and ${\alpha }_{0}={a}_{\mathrm{1,0}}/{a}_{\mathrm{2,0}}$. For small eccentricities and Δ ≈ 0, we have that ${ \mathcal D }\approx {\beta }_{1}\sqrt{{\alpha }_{i}}{e}_{1}^{2}/2+{\beta }_{2}{e}_{2}^{2}/2$.

In order to implement the equations of motion given in Equation (3) governing resonant dynamics, we evaluate the integral in Equation (A4) and its derivative with respect to canonical variables by numerical quadrature using a Gauss–Legendre quadrature rule. We use the exoplanet (Foreman-Mackey et al. 2019) package's Kepler solver algorithm to compute the planets' ${{\boldsymbol{r}}}_{i}$ as functions of canonical variables to evaluate the integrand of Equation (A4). Derivatives of $\bar{H}$ are then computed using the Theano (Theano Development Team 2016) package's automatic differentiation capabilities.

The orbital elements appearing in the canonical variables of the equations of motion governing the planets' resonant dynamics represent "mean" elements. These mean elements differ from the planets' osculating elements, which exhibit additional high-frequency variations on the timescale of the planets' synodic period. This is because the Hamiltonian system governing the resonant dynamics is derived from a phase-space reduction of the full three-body problem via a near-identity canonical transformation that eliminates the degree of freedom associated with the canonical coordinate $Q=({\lambda }_{2}-{\lambda }_{1})/k$. The ACR equilibrium configurations of the resonant Hamiltonian correspond to periodic orbits of the full three-body problem that are 2π periodic in Q. In order to properly initialize N-body simulations with planets residing in an ACR configuration, the periodic variations of the osculating orbital elements must be accounted for; ignoring these corrections can result in large libration amplitudes when running N-body simulations. We do this by constructing the canonical transformation from the variables of our resonance model to the canonical coordinates of the unaveraged three-body problem.

In order to construct our canonical transformation, we begin by writing the full Hamiltonian of the planar three-body problem as

$$\begin{eqnarray}&&\begin{array}{l}H({\sigma }_{i},{I}_{i},Q,P;L)={H}_{0}({I}_{i},P;L)+\epsilon {\bar{H}}_{1}({\sigma }_{i},{I}_{i},P;L)\\ \quad +\,\epsilon {H}_{1,\mathrm{osc}.}({\sigma }_{i},{I}_{i},Q,P;L),\end{array}\end{eqnarray} \tag{ A6 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{H}_{0} & = & -\displaystyle \sum _{i=1}^{2}\displaystyle \frac{{G}^{2}{\tilde{M}}_{i}^{2}{\tilde{m}}_{i}^{3}}{2{{\rm{\Lambda }}}_{i}^{2}}\\ \epsilon {\bar{H}}_{1} & = & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{2\pi }{\displaystyle \int }_{0}^{2\pi }\left(\displaystyle \frac{1}{| {r}_{1}-{r}_{2}| }-\displaystyle \frac{{r}_{1}\cdot {r}_{2}}{{r}_{2}^{3}}\right){dQ}\\ \epsilon {H}_{1,\mathrm{osc}.} & = & {{Gm}}_{1}{m}_{2}\left(\displaystyle \frac{1}{| {r}_{1}-{r}_{2}| }-\displaystyle \frac{{r}_{1}\cdot {r}_{2}}{{r}_{2}^{3}}\right)-{\bar{H}}_{1}\end{array}\end{eqnarray} \tag{ A7 }$$

and serves as a bookkeeping parameter that we will set to unity after deriving the transformation to first order in . We focus on constructing the transformation for phase-space points in the vicinity of an ACR. We adopt the vector notation ${\boldsymbol{z}}=({\sigma }_{1},{\sigma }_{2},{I}_{1},{I}_{2})$ and define an equilibrium point ${\bar{{\boldsymbol{z}}}}_{\mathrm{eq}}({P}_{0})$ as a point that satisfies

$$\begin{eqnarray*}&&{{\rm{\nabla }}}_{{\boldsymbol{z}}}\left[{H}_{0}({\bar{{\boldsymbol{z}}}}_{\mathrm{eq}},{P}_{0};L)+\epsilon {\bar{H}}_{1}({\bar{{\boldsymbol{z}}}}_{\mathrm{eq}},{P}_{0};L)\right]=0.\end{eqnarray*}$$

We next transform to canonical variables $\delta {\boldsymbol{z}}={\boldsymbol{z}}-{\bar{{\boldsymbol{z}}}}_{\mathrm{eq}}$ and $\delta P=P-{P}_{0}$ centered on the equilibrium configuration and approximate the Hamiltonian by expanding in $\delta {\boldsymbol{z}}$ as

$$\begin{eqnarray}\begin{array}{rcl}H & = & \displaystyle \frac{1}{2}\delta {{\boldsymbol{z}}}^{T}\cdot {{\rm{\nabla }}}_{{\boldsymbol{z}}}^{2}\left[{H}_{0}+\epsilon {\bar{H}}_{1}\right]\cdot \delta {\boldsymbol{z}}+\left({\omega }_{\mathrm{syn}}+\delta {\boldsymbol{z}}\cdot {{\rm{\nabla }}}_{{\boldsymbol{z}}}{\omega }_{\mathrm{syn}}\right)\\ & & \times \delta P/k+\epsilon \left[{H}_{1,\mathrm{osc}.}+\delta {{\boldsymbol{z}}}^{T}\cdot {{\rm{\nabla }}}_{{\boldsymbol{z}}}{H}_{1,\mathrm{osc}.}\right],\end{array}\end{eqnarray} \tag{ A8 }$$

where ${\omega }_{\mathrm{syn}}\equiv k\displaystyle \frac{\partial {H}_{0}}{\partial P}$ is the synodic frequency n₂ − n₁, and we have dropped constant terms from Equation (A8). We next perform a symplectic linear transformation,

$$\begin{eqnarray}&&{\boldsymbol{w}}={S}^{-1}\cdot \delta {\boldsymbol{z}},\end{eqnarray} \tag{ A9 }$$

to new canonical variables ${\boldsymbol{w}}=({w}_{1},{w}_{2},{\tilde{w}}_{1},{\tilde{w}}_{2})$, where w_i are the new coordinates and ${\tilde{w}}_{i}$ are the conjugate momenta. The matrix S satisfies

$$\begin{eqnarray}\begin{array}{rcl}{S}^{-1}\cdot \left({\rm{\Omega }}\cdot {{\rm{\nabla }}}_{{\boldsymbol{z}}}^{2}\left[{H}_{0}+\epsilon {\bar{H}}_{1}\right]\right)\cdot S & = & D\\ {S}^{T}\cdot {\rm{\Omega }}\cdot S & = & {\rm{\Omega }},\end{array}\end{eqnarray} \tag{ A10 }$$

where D is a diagonal matrix of the form $D\,=\mathrm{diag}(i{\omega }_{1},i{\omega }_{2},-i{\omega }_{1},-i{\omega }_{2})$. We choose S so that the components of w satisfy ${\tilde{w}}_{i}=-{{iw}}_{i}^{* }$. After the canonical change of variables, Hamiltonian (A8) transforms to

$$\begin{eqnarray}&&\begin{array}{l}H({\boldsymbol{w}},\delta P,Q)=i\displaystyle \sum _{i=1}^{2}{\omega }_{i}{w}_{i}{\tilde{w}}_{i}+\left({\omega }_{\mathrm{syn}}+{\boldsymbol{w}}\cdot {{\rm{\nabla }}}_{{\boldsymbol{w}}}{\omega }_{\mathrm{syn}}\right)\delta P/k\\ \,+\,\epsilon \left[{H}_{1,\mathrm{osc}.}({\bar{{\boldsymbol{z}}}}_{\mathrm{eq}},Q)+{\boldsymbol{w}}\cdot {{\rm{\nabla }}}_{{\boldsymbol{w}}}{H}_{1,\mathrm{osc}.}(Q)\right],\end{array}\end{eqnarray} \tag{ A11 }$$

where ${{\rm{\nabla }}}_{{\boldsymbol{w}}}={S}^{T}\cdot {{\rm{\nabla }}}_{{\boldsymbol{z}}}$. This transformation puts the Hamiltonian in the form of a perturbed pair of harmonic oscillators with frequencies ω_i and allows us to apply canonical perturbation theory to construct the transformation from osculating to mean variables below. To construct the transformation from δz to w, we evaluate the matrix ${\rm{\Omega }}\cdot {{\rm{\nabla }}}_{{\boldsymbol{z}}}^{2}\,\left[{H}_{0}+\epsilon {\bar{H}}_{1}\right]$ by means of automatic differentiation and then numerically determine its eigenvectors in order to calculate the matrices S and D.

We now construct a canonical transformation $({\boldsymbol{w}},\delta P,Q)\,\to (\bar{{\boldsymbol{w}}},\delta \bar{P},\bar{Q})$ by means of a canonical Lie transformation (e.g., Morbidelli 2002) in order to eliminate the dependence of the Hamiltonian on Q to first order in . The transformation is generated by the function ${\chi }_{1}(\bar{{\boldsymbol{w}}},\delta \bar{P},\bar{Q})$ so that $({\boldsymbol{w}},\delta P,Q)=\exp \,[\epsilon {{ \mathcal L }}_{{\chi }_{1}}](\bar{{\boldsymbol{w}}},\delta \bar{P},\bar{Q})$, where ${{ \mathcal L }}_{{\chi }_{1}}=\left\{\cdot ,{\chi }_{1}\right\}$ is the Lie derivative with respect to χ₁, i.e., the canonical Poisson bracket of a given function of phase-space coordinates and the function χ₁. We will solve for χ₁ so that Hamiltonian (A11) is transformed to

$$\begin{eqnarray}&&\begin{array}{l}\bar{H}(\bar{{\boldsymbol{w}}},\delta \bar{P},\bar{Q})\equiv \exp \,[\epsilon {{ \mathcal L }}_{{\chi }_{1}}]H(\bar{{\boldsymbol{w}}},\delta \bar{P},\bar{Q})=i\displaystyle \sum _{i=1}^{2}{\omega }_{i}{\bar{w}}_{i}{\bar{\tilde{w}}}_{i}\\ \,+\,\left({\omega }_{\mathrm{syn}}+\bar{{\boldsymbol{w}}}\cdot {{\rm{\nabla }}}_{{\boldsymbol{w}}}{\omega }_{\mathrm{syn}}\right)\delta \bar{P}/k+{ \mathcal O }({\epsilon }^{2}).\end{array}\end{eqnarray} \tag{ A12 }$$

Expanding Equation (A12) and collecting terms that are first order in , we obtain the following equation for χ₁:⁹

$$\begin{eqnarray}&&\begin{array}{l}i\displaystyle \sum _{i=1}^{2}{\omega }_{i}\,\left({\bar{\tilde{w}}}_{i}\displaystyle \frac{\partial {\chi }_{1}}{\partial {\bar{\tilde{w}}}_{i}}-{\bar{w}}_{i}\displaystyle \frac{\partial {\chi }_{1}}{\partial {\bar{w}}_{i}}\right)+\displaystyle \frac{1}{k}\left({\omega }_{\mathrm{syn}}+\bar{{\boldsymbol{w}}}\cdot {{\rm{\nabla }}}_{{\boldsymbol{w}}}{\omega }_{\mathrm{syn}}\right)\displaystyle \frac{\partial {\chi }_{1}}{\partial \bar{Q}}\\ +\,{H}_{1,\mathrm{osc}.}+\bar{{\boldsymbol{w}}}\cdot {{\rm{\nabla }}}_{\bar{{\boldsymbol{w}}}}{H}_{1,\mathrm{osc}.}=0.\end{array}\end{eqnarray} \tag{ A13 }$$

Inserting the ansatz solution

$$\begin{eqnarray}&&{\chi }_{1}(\bar{Q},{\bar{w}}_{i},{\bar{\tilde{w}}}_{i})=G(\bar{Q})+\displaystyle \sum _{i}{\bar{w}}_{i}{F}_{i}(\bar{Q})+{\bar{\tilde{w}}}_{i}{\tilde{F}}_{i}(\bar{Q})\end{eqnarray} \tag{ A14 }$$

for χ₁ into Equation (A13), we have

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{\mathrm{syn}}}{k}\displaystyle \frac{{dG}(Q)}{{dQ}}={H}_{1,\mathrm{osc}.}({\bar{{\boldsymbol{z}}}}_{\mathrm{eq}},Q)\end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{\mathrm{syn}}}{k}\displaystyle \frac{{{dF}}_{i}(Q)}{{dQ}}=\displaystyle \frac{\partial {H}_{1,\mathrm{osc}.}(Q)}{\partial {w}_{i}}-\displaystyle \frac{1}{k}\displaystyle \frac{\partial {\omega }_{\mathrm{syn}}}{\partial {w}_{i}}\displaystyle \frac{{dG}(Q)}{{dQ}}+i{\omega }_{i}{F}_{i}(Q)\end{eqnarray} \tag{ A16 }$$

and ${\tilde{F}}_{i}=-{{iF}}_{i}^{* }$. To solve for F_i, we Fourier decompose F_i and ${H}_{1,\mathrm{osc}.}$ as ${F}_{i}={\sum }_{l}{\hat{F}}_{i,l}\exp \,[{ilQ}]$ and

$$\begin{eqnarray*}&&\left[\displaystyle \frac{\partial }{\partial {\bar{w}}_{i}}-\displaystyle \frac{\partial \mathrm{ln}{\omega }_{\mathrm{syn}}}{\partial {\bar{w}}_{i}}\right]{H}_{1,\mathrm{osc}.}({\bar{{\boldsymbol{z}}}}_{\mathrm{eq}},Q)=\displaystyle \sum _{l}{\hat{X}}_{i,l}\exp \,[{ilQ}].\end{eqnarray*}$$

We use a fast Fourier transform algorithm to compute the amplitudes ${\hat{X}}_{i,l}$ and then compute solutions for F₁(Q) and F₂(Q) as Fourier sums with amplitudes

$$\begin{eqnarray}&&{\hat{F}}_{i,l}=\displaystyle \frac{-{ik}}{l{\omega }_{\mathrm{syn}}-k{\omega }_{i}}{\hat{X}}_{i,l}\,.\end{eqnarray} \tag{ A17 }$$

We convert the canonical variables of the ACR equilibrium solutions of the resonance Hamiltonian to the canonical variables of the unaveraged three-body problem which are then, in turn, used to compute the osculating orbital elements of planets in the ACR configuration. To first order in , this is achieved by evaluating the function $\bar{x}\mapsto x=\bar{x}+\{\bar{x},{\chi }_{1}\}$ at the phase-space point $(\bar{{\boldsymbol{w}}},\delta \bar{P},\bar{Q})=(0,0,Q)$ for any dynamical variable or function of dynamical variables, x, of interest. For a given ACR equilibrium solution of the resonance Hamiltonian, ${\bar{{\boldsymbol{z}}}}_{\mathrm{eq}}({P}_{0})$, Equations (A9), (A14), and (A15), give the solutions

$$\begin{eqnarray}&&P(Q)={P}_{0}+\{\delta \bar{P},{\chi }_{1}\}{}_{\bar{{\boldsymbol{w}}}=0}={P}_{0}+\displaystyle \frac{k}{{\omega }_{\mathrm{syn}}}{H}_{1,\mathrm{osc}.}({\bar{{\boldsymbol{z}}}}_{\mathrm{eq}},Q)\end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray}&&{\boldsymbol{z}}(Q)={\bar{{\boldsymbol{z}}}}_{\mathrm{eq}}+\left\{\delta \bar{{\boldsymbol{z}}},{\chi }_{1}\right\}={\bar{{\boldsymbol{z}}}}_{\mathrm{eq}}+S\cdot {\rm{\Omega }}\cdot \left(\begin{array}{l}{F}_{1}(Q)\\ {F}_{2}(Q)\\ -{{iF}}_{1}^{* }(Q)\\ -{{iF}}_{2}^{* }(Q)\end{array}\right)\end{eqnarray} \tag{ A19 }$$

for osculating canonical variables, which can then be used to compute the planets' osculating orbital elements. Figure A1 shows an example of osculating elements computed as a function of Q using the transformations given in Equations (A18) and (A19) for a pair of Jupiter-mass planets in a 3:2 ACR. Numerical routines for reproducing Figure A1 as well as computing the transformation from mean to osculating variables for different resonances, planet masses, and normalized AMD values are available online at github.com/shadden/ResonantPlanetPairsRVModeling.

Zoom In Zoom Out Reset image size