8. Numerical Resonance Models

8.1. Introduction

Approximate equations of motion governing the dynamics of a mean motion resonance can be generated by including terms from a distrubing function expansion using the PoincareHamiltonian class.

However, celmech also allows users to derive equations governing mean motion resonance dynamics without resorting to a truncated Taylor expansions in eccentricities and inclinations. celmech provides two such classes, PlanarResonanceEquations and SpatialResonanceEquations, are described below.

8.2. Planar Resonance Equations

The PlanarResonanceEquations class provides equations of motion governing a \(j\mathrm{:}j-k\) mean-motion resonance between two coplanar planets. In order to do so, the class first constructs a numerical Hamiltonian governing the system by means of numerical integration.

8.2.1. Introduction

Starting from the full 3-body Hamiltonian governing a pair of planar planets,

\[\begin{split}\begin{align} {\cal H}({\pmb \lambda},{\pmb \gamma};\pmb{\Lambda},\pmb{\Gamma}) &= H_\mathrm{Kep}(\pmb{\Lambda}) + H_\mathrm{int}({\pmb \lambda},{\pmb \gamma};\pmb{\Lambda},\pmb{\Gamma})\\ H_\mathrm{Kep}(\Lambda_i) &= -\sum_{i=1}^2\frac{G^2(M_*+m_i)^2\mu_i^3}{2\Lambda_i^2} \\ H_\mathrm{int}(\lambda_i,\gamma_i;\Lambda_i,\Gamma_i) &= -\frac{Gm_1m_2}{|{\bf r}_1-{\bf r}_2|} + \frac{\tilde {\bf r}_1 \cdot \tilde {\bf r}_2}{M_*}~, \end{align}\end{split}\]

the first step in constructing this Hamiltonian is to perform a canonical transformation to new canonical coordinates

\[\begin{split}\begin{align} \begin{pmatrix}\sigma_1\\ \sigma_2 \\ Q \\ l \end{pmatrix} = \begin{pmatrix} -s & 1+s & 1 & 0 \\ -s & 1+s & 0 & 1 \\ -1/k & 1/k & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0 & 0 \end{pmatrix} \cdot \begin{pmatrix}\lambda_1 \\ \lambda_2 \\ \gamma_1 \\ \gamma_2 \end{pmatrix} \end{align}\end{split}\]

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

\[\begin{split}\begin{eqnarray} \Gamma_i &=& I_i \nonumber\\ \Lambda_1 &=& \frac{L}{2} - P/k - s (I_1 + I_2) \nonumber \\ \Lambda_2 &=& \frac{L}{2} + P/k + (1+s) (I_1 + I_2)~.\label{eq:equations_of_motion:momenta} \end{eqnarray}\end{split}\]

Conservation of angular momentum implies that the quantity

\[L = \Lambda_1 + \Lambda_2 + \Gamma_1 + \Gamma_2~,\]

conjugate to the angle \(l\), is conserved

Next, we separate \(H_\mathrm{int}\) into two pieces,

\[\begin{split}H_\mathrm{int} &= \bar{H}_\mathrm{int}(\sigma_i,I_i,P; L) + H_\mathrm{int,osc.}(\sigma_i,I_i,Q,P; L) \\ \bar{H}_\mathrm{int}(\sigma_i,I_i,P; L) &= \frac{1}{2\pi}\int_{0}^{2\pi} H_\mathrm{int} dQ\\ H_\mathrm{int,osc.} &= H_\mathrm{int} - \bar{H}_\mathrm{int}\end{split}\]

The PlanarResonanceEquations class then models the dynamics of the resonant Hamiltonian,

\[{\cal H}_\mathrm{res}(\sigma_i,I_i;P , L) = H_\mathrm{Kep}(I_i ; P, L) + \bar{H}_\mathrm{int}(\sigma_i,I_i;P , L)\]

where \(P\) is now a conserved quantity because \({\cal H}_\mathrm{res}\) is independent of the coordinate \(Q\).

8.2.2. Units and Dynamical Variables

As written, \({\cal H}_\mathrm{res}(\sigma_i,I_i;P , L)\) represents a two degree-of-freedom system that depends on two parameters \(P \mathrm{ and } L\) (in addition to the planets’ and star’s masses). The PlanarResonanceEquations class reduces the number of free parameters by selecting appropriate units for action variables and time. In particular, the units are chosen by first defining the reference semi-major axes, \(a_{i,0}\), so that

\[\begin{split}L = \tilde{m}_1 \sqrt{G\tilde{M}_1a_{1,0}} + \tilde{m}_2 \sqrt{G\tilde{M}_2a_{2,0}}\\ a_{1,0}/a_{2,0} = \left(\frac{j-k}{j}\right)^{2/3}\left(\frac{\tilde{M}_1}{\tilde{M}_2}\right)^{1/3}\end{split}\]

Next, the action variables are rescaled by

\[I_i \rightarrow \frac{1}{(\tilde{m}_1 + \tilde{m}_2)\sqrt{GM_*a_{2,0}}}I_i\]

and time is measured in units such that \(\sqrt{\frac{GM_*}{a_{2,0}^{3}}}=1\) (i.e., the orbital period at \(a=a_{2,0}\) is equal to \(2\pi\)). Finally, rather than specifying the conserved quantity \(P\), the equations of motion used by PlanarResonanceEquations are formulated in terms of

\[{\cal D} = \frac{\Lambda_{2,0}-\Lambda_{1,0} - \Lambda_{2} + \Lambda_{1}}{s+1/2} + \Gamma_1 + \Gamma_2\]

where \(\Lambda_{i,0} = \tilde{m}_i \sqrt{G\tilde{M}_ia_{i,0}}\). In terms of orbital elements,

\[{\cal 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) - \frac{k\beta_2\beta_1 \sqrt{\alpha_0} \Delta }{3 \left(\beta_1\sqrt{\alpha_0} (s+1)+\beta_2 s\right)}\]

where \(\Delta = \frac{j-k}{j}\frac{P_2}{P_1} - 1\), \(\beta_i = \frac{\tilde{m}_i}{\tilde{m}_1+\tilde{m}_2}\sqrt{\frac{\tilde{M}_i}{M_*}}\) and \(\alpha_0 = \frac{a_{1,0}}{a_{2,0}}\).

Finally, PlanarResonanceEquations uses canonical coordinate \(y_i = \sqrt{2I_i}\sin\sigma_i\) and conjugate momenta \(x_i = \sqrt{2I_i}\cos\sigma_i\) instead of \((\sigma_i,I_i)\) in order to formulate the equations of motion.

Various methods of the PlanarResonanceEquations class take dynamical variables as a vector in the standardized form:

\[\pmb{z} = (y_1,y_2,x_1,x_2,{\cal D})\]

when computing equations of motion and other quantities.

8.2.3. Dissipative Forces

The PlanarResonanceEquations class includes capabilities to model dissipative forces inducing migration and eccentricity damping in addition to the conservative dynamics of resonant interactions. Specifically, migration forces of the form

\[\begin{split}\begin{align} \frac{d\ln a_i}{dt} &= -\frac{1}{\tau_{m,i}} - 2p\frac{e_i^2}{\tau_{e,i}}\\ \frac{d\ln e_i}{dt} &= -\frac{1}{\tau_{e,i}} \end{align}\end{split}\]

These forces are specified by setting the attributes K1, K2, p, and tau_alpha.

8.3. Spatial Resonance Equations

The SpatialResonanceEquations class provides equations of motion governing a mean-motion resonance between a planet pair on mutually inclined orbits. To study the dynamics of a \(j:j-k\) MMR, we will perform a canonical transformation from the usual Poincare canonical variables to new canonical variables as in the planar case above. However, before proceeding with this transformation, we follow Malige (2002) and use the invariance of the system’s angular momentum vector direction to reduce the number of degrees of freedom by one. This reduction is accomplished by first adopting a coordinate system in which the system’s angular momentum vector lies along the \(z\)-axis. Next, complex canonical variables \(y_i=\sqrt{Q_i}e^{-\mathrm{i} q_i}\) are introduced and the reduced Hamiltonian is written by replacing the complex canonical variables \(y_i\) with the expressions

\[\begin{split}\begin{eqnarray} y_1 &=& \mathrm{i} \frac{y}{\sqrt{2}} \sqrt{1 + \frac{\Lambda_2 - \Gamma_2 - \Lambda_1 + \Gamma_1} {\Lambda_1 + \Lambda_2 - \Gamma_2 - \Gamma_1 - y\bar{y}} } \\ y_2 &=& -\mathrm{i} \frac{y}{\sqrt{2}} \sqrt{1 + \frac{\Lambda_2 - \Gamma_2 - \Lambda_1 + \Gamma_1} {\Lambda_1 + \Lambda_2 - \Gamma_2 - \Gamma_1 - y\bar{y}} }~, \end{eqnarray}\end{split}\]

which introduce the new complex canonical variable \(y\). Finally, the new canonical variable \(y\) can be written as \(y=\sqrt{Q}e^{-\mathrm{i}q}\), introducing the canonical coordinate-momentum pair \((q,Q)\). In terms of orbital elements, the new variables are

\[Q=2\mu_1\sqrt{M_1a_1(1-e_1^2)}\sin^2(I_1/2) +2\mu_2\sqrt{M_2a_2(1-e_2^2)}\sin^2(I_2/2)\]

and \(q=-\frac{1}{2}(\Omega_1+\Omega_2)\).

Now we introduce the new canonical angle variables

\[\begin{split}\begin{align} \begin{pmatrix}\sigma_1\\ \sigma_2 \\ \phi \\ \psi \\ l \end{pmatrix} = \begin{pmatrix} -s & 1+s & 1 & 0 & 0 \\ -s & 1+s & 0 & 1 & 0 \\ -s & 1+s & 0 & 0 & 1 \\ -1/k & 1/k & 0 & 0 & 0 \\ -s & 1+s & 0 & 0 & 0 \end{pmatrix} \cdot \begin{pmatrix}\lambda_1 \\ \lambda_2 \\ \gamma_1 \\ \gamma_2 \\ q \end{pmatrix} \label{eq:equations_of_motion:new_angles} \end{align}\end{split}\]

where \(s = (j-k)/{k}\), along with their conjugate action variables

\[\begin{split}\begin{align} \begin{pmatrix}I_1 \\ I_2 \\ \Phi \\ \Psi \\ L \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ -j & k-j & 0 & 0 & 0 \\ 1 & 1 & -1 & -1 & -1 \end{pmatrix} \cdot \begin{pmatrix}\Lambda_1 \\ \Lambda_2 \\ \Gamma_1 \\ \Gamma_2 \\ Q \end{pmatrix}~. \label{eq:equations_of_motion:new_actions} \end{align}\end{split}\]

After this canonical transformation, the Hamiltonian is independent of the angle variables \(l\) so the quantity

\[\begin{eqnarray} L &=& \Lambda_1+\Lambda_2-Q-\Gamma_1-\Gamma_2 \end{eqnarray}\]

is conserved, reflecting the conservation of the system’s total angular momentum. The phase space is further reduced by numerically averaging over the fast angle \(\psi\), similar to the PlanarResonanceEquations class. We also eliminate the explicit dependence of our equations of motion on the value of the total angular momentum by introducing a simultaneous re-scaling of the Hamiltonian and the action variables that preserves the form of Hamilton’s equations. In order to do so, let us first define the reference semi-major axes \(a_{1,0}\) and \(a_{2,0}\) such that

\[\begin{equation*} L=\mu_1\sqrt{GM_1a_{1,0}} + \mu_2\sqrt{GM_2a_{2,0}} \end{equation*}\]

and

\[\begin{equation*} a_{1,0}=\left(\frac{M_1}{M_2}\right)^{1/3}\left(\frac{j-k}{j}\right)^{2/3}a_{2,0}~. \end{equation*}\]

In other words, \(a_{1,0}\) is defined as the nominal location of an exact resonance with a planet situated at \(a_{2,0}\). We then re-scale the Hamiltonian and the action variable, dividing each by a factor of \((\mu_1 + \mu_2)\sqrt{GM_*a_{2,0}}\) so that the new Hamiltonian and canonical action variables are given by

\[\begin{eqnarray} \{\mathcal{H},{I}_1,{I}_2,{\Phi},{\Psi}\} \rightarrow \frac{1}{(\mu_1+\mu_2)\sqrt{GM_*a_{2,0}}} \times \{\mathcal{H},I_1,I_2,\Phi,\Psi\}~. \end{eqnarray}\]

and the re-scaled angular momentum, \({L}= \beta_2 + \beta_1\sqrt{a_{1,0}/a_{2,0}}\), is a constant. In terms of the newly re-scaled variables, the averaged Hamiltonian is then

\[\mathcal{H}(\sigma_1,\sigma_2,\phi,{I}_1,{I}_2,{\Phi};{\Psi}) = H_\mathrm{kep}({I}_1,{I}_2,{\Phi};{\Psi}) + \epsilon \bar{H}_\mathrm{int}(\sigma_1,\sigma_2,\phi,{I}_1,{I}_2,{\Phi};{\Psi})\]

where

\[\begin{split}\begin{eqnarray} H_\mathrm{kep} &=& -n_{2,0}\sum_{i=1}^2\frac{\beta_i^3(1+m_i/M_*)^{1/2}}{2{\Lambda}_i^2}\\ \bar{H}_\mathrm{int}&=&\frac{n_{2,0}}{2\pi}\int_{-\pi}^{\pi}\left(\frac{a_{2,0}}{GM_*}\pmb{v}_1\cdot\pmb{v}_2-\frac{a_{2,0}}{|\pmb{r}_2-\pmb{r}_1|} \right)d\psi \end{eqnarray}\end{split}\]

with \(\beta_i=\frac{\mu_i}{\mu_1+\mu_2}\sqrt{1+m_i/M_*}\), \(n_{2,0} = \sqrt{GM_*a_{2,0}^{-3}}\), \(\epsilon = \frac{m_1m_2}{M_*(\mu_1+\mu_2)}\), and

\[\begin{split}\begin{align*} {\Lambda}_1 &= -s{L}-{\Psi}/k - s({I}_1 + {I}_2 + {\Phi})= \beta_1\sqrt{a_{1}/a_{2,0}}\\ {\Lambda}_2 &= (s+1){L} + {\Psi}/k + (s+1)({I}_1 + {I}_2 + {\Phi})= \beta_2\sqrt{a_{2}/a_{2,0}}~. \end{align*}\end{split}\]

In order to avoid coordinate singularities that occur when \(I_i=0\) or \(\Phi=0\), the SpatialResonanceEquations class formulates the Hamiltonian in terms of the canonical variables

\[\begin{split}\begin{eqnarray} (y_i,x_i) &=& (\sqrt{2I_i}\sin\sigma_i,\sqrt{2I_i}\cos\sigma_i)\nonumber\\ (y_\mathrm{inc},x_\mathrm{inc}) &=& (\sqrt{2\Phi}\sin\phi,\sqrt{2\Phi}\cos\phi)~. \end{eqnarray}\end{split}\]

Because the averaged Hamiltonian is independent of the angle variable \(\psi\), the action variable \(\Psi\) is a constant of motion. Rather than parameterizing the Hamiltonian in terms of \(\Psi\), we introduce a variable \({\cal D}\) that has a more straightforward interpretation as the “average” angular momentum deficit of the resonant planet pair. In order to define \({\cal D}\), let us first introduce \(\rho_\mathrm{res}\) as the value of \({\Lambda}_1/{\Lambda}_2\) when the planet pair’s period ratio is equal to the nominal resonant period ratio. Explicitly, \(\rho_\mathrm{res} = \frac{\beta_1}{\beta_2}\left(\frac{j-k}{j}\right)^{1/3}\left(\frac{M_1}{M_2}\right)^{1/6}\). Then, we define \({\cal D}\), along with an additional constant \(\Lambda_{2,\mathrm{res}}\), such that

\[\begin{split}\begin{equation} \begin{pmatrix} -\Psi/k\\ {L} \end{pmatrix} = \begin{pmatrix} s + (1+s) \rho_\mathrm{res} & 0 \\ (1+\rho_\mathrm{res}) & -1 \end{pmatrix} \cdot \begin{pmatrix} \Lambda_{2,\text{res}} \\ {\cal D} \end{pmatrix}~. \end{equation}\end{split}\]

Inverting this equation and solving for \(\mathcal{D}\) yields

\[\begin{equation} {\cal D} = \frac{1}{2}\left(x_1^2+y_1^2+x_2^2+y_2^2+x_\mathrm{inc}^2+y_\mathrm{inc}^2\right) + (\Lambda_{2,\mathrm{res}}-{\Lambda}_2) + (\rho_\mathrm{res}\Lambda_{2,\mathrm{res}}-{\Lambda}_1) \end{equation}\]

We express \(\mathcal{D}\) in terms of the planet pair’s period ratio by using

\[\Delta \equiv \frac{j-k}{j}\frac{P_2}{P_1} - 1 \approx 3\left(\frac{\delta\hat{\Lambda}_{2}}{\Lambda_{2,\mathrm{res}}}-\frac{\delta\hat{\Lambda}_{1}}{\rho_\mathrm{res}\Lambda_{2,\mathrm{res}}}\right)=3\frac{\delta\hat{\Lambda}_{2}}{\Lambda_{2,\mathrm{res}}}\left(1+\frac{s}{(1+s)\rho_\mathrm{res}}\right)\]

allowing us to write

\[\begin{eqnarray} {\cal D} = \frac{1}{2}\left(x_1^2+y_1^2+x_2^2+y_2^2+x_\mathrm{inc}^2+y_\mathrm{inc}^2\right) - \frac{1}{3} \frac{\rho_\mathrm{res}\Lambda_{2,\mathrm{res}}}{s + \rho_\mathrm{res}(1+s)} \Delta \end{eqnarray}\]

In other words, \(\mathcal{D}\) sets the value of the pair’s angular momentum deficit when the planets’ period ratio is at the exact resonant value (\(\Delta=0\)).

8.4. API

class celmech.numerical_resonance_models.PlanarResonanceEquations(j, k, n_quad_pts=40, m1=1e-05, m2=1e-05, K1=100, K2=100, tau_alpha=100000.0, p=1)

A class for the model describing the dynamics of a pair of planar planets in/near a mean motion resonance.

Includes the effects of dissipation.

j

Together with k specifies j:j-k resonance

Type

int

k

Order of resonance.

Type

int

alpha

Semi-major axis ratio \(a_1/a_2\)

Type

float

eps

Mass parameter m1*mu2 / (mu1+mu2)

Type

float

m1

Inner planet mass

Type

float

m2

Outer planet mass

Type

float

tau_alpha

Migration timescale defined by \(\tau_\alpha^{-1} =\tau_{m,2}^{-1} - \tau_{m,1}^{-1}\)

Type

float

K1

Ratio of inner planet eccentricity damping timescale to migration timescale, tau_alpha.

Type

float

K2

Ratio of outer planet eccentricity damping timescale to migration timescale, tau_alpha.

Type

float

p

Sets coupling between eccentricity damping and semi-major axis damping so that

\[\frac{d}{dt}a_i = -a_i/\tau_{m,i} - 2pe_i^2/\tau_{e,i}\]

A value of p=1 corresponds to eccentricity damping at constant angular momentum.

Type

float

timescales

Dictionary containing the migration and eccentricity damping time-scales of the inner and outer planets based on the damping parameters of the resonance model.

Type

dict

H(z)

Calculate the value of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The value of the Hamiltonian evaluated at z.

Return type

float

H_flow(z)

Calculate flow induced by the Hamiltonian

\[\dot{z} = \Omega \cdot \nabla_{z}H(z)\]

Parameters

z (ndarray) – Dynamical variables

Returns

Flow vector

Return type

ndarray

H_flow_jac(z)

Calculate the Jacobian of the flow induced by the Hamiltonian

\[\Omega \cdot \Delta H(z)\]

Parameters

z (ndarray) – Dynamical variables

Returns

Jacobian matrix

Return type

ndarray

H_kep(z)

Calculate the value of the Keplerian component of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The value of the Keplerian part of the Hamiltonian evaluated at z.

Return type

float

H_kep_flow(z)

Calculate the flow vector generated by the Keplerian component of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The flow vector of the Keplerian component of the Hamiltonian

Return type

ndarray

H_pert(z)

Calculate the value of the perturbation component of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The value of the perturbation part of the Hamiltonian evaluated at z.

Return type

float

H_pert_flow(z)

Calculate the flow vector generated by the perturbation component of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The flow vector generated by the perturbation part of the Hamiltonian evaluated at z.

Return type

ndarray

H_pert_full(psi, z)

Calculate the value of the perturbation component of the full, unaveraged Hamiltonian.

Parameters

• psi (float) –

  Angle variable

  \[(\lambda_2-\lambda_1)/k\]

• z (ndarray) – Dynamical variables

Returns

The value of the perturbation part of the Hamiltonian evaluated at \((\psi,z)\).

Return type

float

H_pert_osc(psi, z)

Calculate the oscillating part of value of the perturbation component of the unaveraged Hamiltonian.

Parameters

• psi (float) –

  Angle variable

  \[(\lambda_2-\lambda_1)/k\]

• z (ndarray) – Dynamical variables

Returns

The value of the perturbation part of the Hamiltonian evaluated at \((\psi,z)\).

Return type

float

dyvars_from_rebound_simulation(sim, iIn=1, iOut=2, osculating_correction=False)

Convert dynamical variables

to a Rebound simulation.

Parameters

• sim (rebound.Simulation) – Simulation object to

• iIn (int, optional) – Integer index of inner planet to use for computing dynamical variables. Default is 1

• iOut (int, optional) – Integer index of outer planet to use for computing dynamical variables. Default is 2

• osculating_correction (boole, optional) – If True, correct the orbital elements of the rebound simulation by transforming to the mean variables of the resonance model.

Returns

z – Dynamical variables

Return type

ndarray

dyvars_to_Poincare_actions(dyvars)

Convert a set of dynamical variables to Poincare actions

Parameters

dyvars (ndarray) – Dynamical variables to convert

Returns

actions – Poincare actions stored as dictionary entries.

Return type

dict

dyvars_to_orbital_elements(z)

Convert dynamical variables

\[z = (y_1,y_2,x_1,x_2,{\cal D})\]

to orbital elements

\[(a_1,e_1,\theta_1,a_2,e_2,\theta_2)\]

dyvars_to_rebound_simulation(z, Q=0, pomega1=0, osculating_correction=True, include_dissipation=False, **kwargs)

Convert dynamical variables

\[z = (y_1,y_2,x_1,x_2,{\cal D}\]

to a Rebound simulation.

Parameters

• z (ndarray) – Dynamical variables

• pomega1 (float, optional) – Planet 1’s longitude of periapse. Default is 0

• Q (float, optional) – Angle variable Q = (lambda2-lambda1) / k. Default is Q=0. Resonant periodic orbits are 2pi-periodic in Q.

• osculating_correction (bool, optional) – If True, apply correction from the ‘averaged’ variables employed by the resonance model to osculating variables used in the un-averaged three-body problem. Correction is done to first order in planet masses and assumes that the system is near an equilibrium configuration of the resonance model. Default is True.

• include_dissipation (bool, optional) – Include dissipative effects through reboundx’s external forces. Default is False

Keyword Arguments

inc (float, default = 0.0) – Inclination of planets’ orbits.

Returns

Returns a tuple. The first item of the tuple is a rebound simulation. The second item is a reboundx.Extras object if ‘include_dissipation’ is True, otherwise the second item is None.

Return type

tuple

find_equilibrium(guess, dissipation=False, tolerance=1e-09, max_iter=10)

Use Newton’s method to locate an equilibrium solution of the equations of motion.

By default, an equilibrium of the dissipation-free equations is sought. In this case, the AMD value of the equilibrium solution will be equal to the value of the initially supplied guess of dynamical variables. If dissipative terms are included then the equilibrium will depend on the parameters K1 and K2 as well as tau_alpha.

Parameters

• guess (ndarray) – Initial guess for the dynamical variables at the equilibrium.

• dissipation (bool, optional) – Whether dissipative terms are considered in the equations of motion. Default is False.

• tolerance (float, optional) – Tolerance for root-finding such that solution satisfies \(|f(z)|\) < tolerance Default value is 1E-9.

• max_iter (int, optional) – Maximum number of Newton’s method iterations. Default is 10. If maximum is reached, result will be returned with a warning.

• include_dissipation (bool, optional) – Include dissipative terms in the equations of motion. Default is False

Returns

zeq – Equilibrium value of dynamical variables.

Return type

ndarray

flow(z)

Calculate the flow vector of the equations of motion

\[\dot{z} = \Omega \cdot \nabla_{z}H(z) + f_{dis}(z)\]

Parameters

z (ndarray) – Dynamical variables

Returns

Flow vector

Return type

ndarray

flow_jac(z)

Calculate the Jacobian of the equations of motion given by

\[\Omega \cdot \Delta H(z) + \nabla f_{dis}(z)\]

Parameters

z (ndarray) – Dynamical variables

Returns

Jacobian matrix

Return type

ndarray

grad_H_pert_full(psi, z)

Calculate the gradient of the perturbation component of the full, unaveraged Hamiltonian with respect to variables \(z\).

Parameters

• psi (float) –

  Angle variable

  \[(\lambda_2-\lambda_1)/k\]

• z (ndarray) – Dynamical variables

Returns

grad – The gradient the perturbation part of the Hamiltonian with respect to \(z\) evaluated at \((\psi,z)\).

Return type

ndarray

grad_omega_syn(z)

Calculate the gradient of the synodic frequency,

\[\nabla_{\pmb{z}} \omega_{syn}\]

as a function of dynamical variables \(\pmb{z}\)

Parameters

z (ndarray) – Dynamical variables

Returns

Synodic frequency

Return type

float

integrate_initial_conditions(dyvars0, times, dissipation=False)

Integrate initial conditions and calculate dynamical variables and orbital elements at specified times.

Integration is done with the scipy.integrate.odeint method.

Parameters

• dyvars0 (ndarray) – Initial conditions in the form of an array of dynamical variables.

• times (ndarray) – Times at which to calculate output.

Returns

soln – Dictionary containing both the dynamical variables and the orbital elements of the integrated solution. The dictionary keys are - ‘times’: Times at which solutin is computed. - ‘dynamical_variables’: ndarray of dynamical at each of the times in ‘times’. - ‘orbital_elements’: dict containing arrays of the various orbital elements.

Return type

dict

mean_to_osculating_dyvars(Q, z, N=256)

Apply perturbation theory to transfrom from the phase space coordiantes of the averaged model to osculating phase space coordintates of the full phase space. Assumes that the transformation is being applied at the fixed point of the averaged model. Details are described in the appendix of Hadden & Payne (2020)

8. Arguemnts

Qfloat or ndarry

Value(s) of angle \(Q=(lambda2-lambda1)/k\) at which to apply transformation. Equilibrium points of the averaged model correspond to orbits that are periodic in Q in the full phase space.

zndarray

Dynamical variables of the averaged model, math::z = (y_1,y_2,x_1,x_2,{cal D})

Nint, optional

Number of Q points to evaluate functions at when performing Fourier transformation. Default is N=256

returns

zosc – The osculating values of the dynamical varaibles for the input Q values. The dimension of the returned variables is set by the dimension of the input ‘Q’. If Q is a float, then z is an array of length 5. If Q is an array of length M then zosc is an (M,5) array.

rtype

ndarray, (5,) or (M,5)

omega_syn(z)

Calculate the synodic frequency,

\[\omega_{syn} = n_2 - n_1\]

as a function of dynamical variables \(z\)

Parameters

z (ndarray) – Dynamical variables

Returns

Synodic frequency

Return type

float

orbital_elements_to_dyvars(orbels)

Convert orbital elements

\[(a_1,e_1,\theta_1,a_2,e_2,\theta_2)\]

to dynamical variables

\[z = (y_1,y_2,x_1,x_2,{\cal D})\]

class celmech.numerical_resonance_models.SpatialResonanceEquations(j, k, n_quad_pts=40, m1=1e-05, m2=1e-05)

A class for equations describing the dynamics of a pair of planets in/near a mean motion resonance.

j

Together with k specifies j:j-k resonance

Type

int

k

Order of resonance.

Type

int

alpha

Semi-major axis ratio a_1/a_2

Type

float

eps

Mass parameter m1*m2 / (mu1+mu2)

Type

float

m1

Inner planet mass

Type

float

m2

Outer planet mass

Type

float

H(z)

Calculate the value of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The value of the Hamiltonian evaluated at z.

Return type

float

H_flow(z)

Calculate flow induced by the Hamiltonian

\[\dot{z} = \Omega \cdot \nabla_{z}H(z)\]

Parameters

z (ndarray) – Dynamical variables

Returns

Flow vector

Return type

ndarray

H_flow_jac(z)

Calculate the Jacobian of the flow induced by the Hamiltonian

\[\Omega \cdot \Delta H(z)\]

Parameters

z (ndarray) – Dynamical variables

Returns

Jacobian matrix

Return type

ndarray

H_kep(z)

Calculate the value of the Keplerian component of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The value of the Keplerian part of the Hamiltonian evaluated at z.

Return type

float

H_pert(z)

Calculate the value of the perturbation component of the Hamiltonian.

Parameters

z (ndarray) – Dynamical variables

Returns

The value of the perturbation part of the Hamiltonian evaluated at z.

Return type

float

dyvars_from_rebound_simulation(sim, iIn=1, iOut=2, osculating_correction=False, full_output=False)

Convert rebound simulation to dynamical variables.

Parameters

• sim (rebound.Simulation) – Simulation object

• iIn (int, optional) – Index of inner planet in simulation particles list. Default is 1.

• iOut (int, optional) – Index of outer planet in simulation particles list. Default is 2.

• osculating_correction (boole, optional) – If True, correct the orbital elements of the rebound simulation by transforming to the mean variables of the resonance model. Default is False.

• full_output (boole, optional) – Return a dictionary containing the re-scaling factors from time and distance units of the input simulation to the time and distance units used by SpatialResonanceEquations.

Returns

• z (ndarray) – Dynamical variables for resonance equations of motion.

• scales (dict, optional) – Returned if ‘full_output’ is True. ‘scales’ contains conversion factors from simulation to resonance equation units such that:

  distance_sim = scales[‘distance’] * distance_eqs time_sim = scales[‘time’] * time_eqs

dyvars_to_Poincare_actions(z)

Convert dynamical variables .. math:

z = (y1,y2,yinc,x1,x2,xinc,amd)

to canonical Poincare action variables .. math:

(L1,L2,Gamma1,Gamma2,Q1,Q2)

Parameters

z (ndarray) – Dynamical variables

Returns

Dictionary of canonical actions with strings as keys.

Return type

dict

dyvars_to_orbital_elements(z)

Convert dynamical variables .. math:

z = (y_1,y_2,y_inc,x_1,x_2,y_inc,{cal D})

to orbital elements .. math:

(a_1,e_1,inc_1,theta_1,a_2,e_2,inc_2,theta_2)

dyvars_to_rebound_simulation(z, psi=0, Omega=0, osculating_correction=True, include_dissipation=False, **kwargs)

Convert dynamical variables to a Rebound simulation.

Parameters

• z (ndarray) – Dynamical variables

• Omega (float, optional) – Value of cyclic coordinate Omega

• psi (float, optional) – Angle variable psi = (lambda1-lambda2) / k. Default is psi=0. Resonant periodic orbits are 2pi-periodic in Q.

• include_dissipation (bool, optional) – Include dissipative effects through reboundx’s external forces. Default is False

Keyword Arguments

• Mstar (float) – Default=1.0 Stellar mass.

• inc (float, default = 0.0) – Inclination of planets’ orbits.

• period1 (float) – Default = 1.0 Orbital period of inner planet

• units (tuple) – Default = (‘AU’,’Msun’,’days’) Units of rebound simulation.

Returns

Returns a tuple. The first item of the tuple is a rebound simulation. The second item is a reboundx.Extras object if ‘include_dissipation’ is True, otherwise the second item is None.

Return type

tuple

integrate_initial_conditions(dyvars0, times)

Integrate initial conditions and calculate dynamical variables and orbital elements at specified times.

Integration is done with the scipy.integrate.odeint method.

Parameters

• dyvars0 (ndarray) – Initial conditions in the form of an array of dynamical variables.

• times (ndarray) – Times at which to calculate output.

mean_to_osculating_dyvars(psi, z, N=256)

Apply perturbation theory to transfrom from the phase space coordiantes of the averaged model to osculating phase space coordintates of the full phase space.

Assumes that the transformation is being applied at the fixed point of the averaged model.

8. Arguemnts

Qfloat or ndarry

Value(s) of angle (lambda2-lambda1)/k at which to apply transformation. Equilibrium points of the averaged model correspond to orbits that are periodic in Q in the full phase space.

zndarray

Dynamical variables of the averaged model, \(z = (y_1,y_2,y_I,x_1,x_2,x_I,\mathcal{D})\)

Nint, optional

Number of Q points to evaluate functions at when performing Fourier transformation. Default is N=256

returns

zosc – The osculating values of the dynamical varaibles for the input Q values. The dimension of the returned variables is set by the dimension of the input ‘Q’. If Q is a float, then z is an array of length 5. If Q is an array of length M then zosc is an (M,5) array.

rtype

ndarray, (7,) or (M,7)