4. Disturbing Function

4.1. Introduction

The full Hamiltonian for a system of \(N\) planets can be written in terms of canonical heliocentric variables as

\[H = -\sum_{i=1}^N\frac{G^2(M_*+m_i)^2\mu_i^3}{2\Lambda_i^2} + \sum_{i=1}^N\sum_{j=1}^i \left( -\frac{Gm_im_j}{|\pmb{r}_i-\pmb{r}_j|} + \frac{\tilde {\bf r}_i \cdot \tilde {\bf r}_j}{m_0} \right)\]

The second term, representing the gravitational interaction potential between the planets in the system, makes the problem non-integrable. Pairwise interactions between two planets, \(i\) and \(j\), can be expressed in terms of the so-called ‘disturbing function’, \({\cal R}^{(i,j)}\), defined so that

\[-\frac{Gm_im_j}{|\pmb{r}_i-\pmb{r}_j|} + \frac{\tilde {\bf r}_i \cdot \tilde {\bf r}_j}{m_0} = -\frac{Gm_im_j}{a_j}{\cal R}^{(i,j)}~.\]

The disturbing function, \({\cal R}^{(i,j)}\), can be expanded in a cosine series in the planets’ angular orbital elements. Defining \(\pmb{\theta}_{i,j} = (\lambda_j,\lambda_i,\varpi_i,\varpi_j,\Omega_i,\Omega_j)\), we can write this cosine series as

\[{\cal R}^{(i,j)} = \sum_{\bf k}c_{\pmb{k}}(\alpha,e_i,e_j,I_i,I_j)\cos(\pmb{k}\cdot \pmb{\theta}_{i,j})\]

where \(\alpha = a_i/a_j\). Rotation and reflection symmetries of the planets’ gravitational interactions dictate that \(c_{\bf k}\ne 0\) only if \(\sum_{l=1}^{6}k_l = 0\) and \(k_5+k_6=2n\) where \(n\) is an integer.

4.1.1. Expansion in orbital elements

In classical disturbing function expansions, the cosine amplitudes \(c_{\pmb{k}}\) are further expanded as Taylor series in powers of the planets’ eccentricities \(e\) and \(s = \sin(I_i/2)\) as

\[c_{\bf k}(\alpha,e_i,e_j,I_i,I_j) = e_i^{|k_3|}e_j^{|k_4|}s_i^{|k_5|}s_j^{|k_6|} \sum_{\nu_1,\nu_2,\nu_3,\nu_4=0}^\infty \tilde{C}_{\bf k}^{(\nu_1,\nu_2,\nu_3,\nu_4)}(\alpha) s_i^{2\nu_1}s_j^{2\nu_2}e_i^{2\nu_3}e_j^{2\nu_4}\]

celmech offers the capability to compute the individual disturbing function coefficients through the function df_coefficient_Ctilde().

4.1.2. Expansion in canonical variables

Since celmech formulates equations of motion in terms of canonical variables rather than orbital elements, it is often more convenient to formulate the Taylor series expansion of the disturbing function as

\[\begin{split}c_{\bf k}(\alpha,e_i,e_j,I_i,I_j)\exp[i\pmb{k}\cdot\pmb{\theta}_{i,j}] = \exp[i(k_1\lambda_j + k_2\lambda_i)] \bar{X_i}^{-k_3} \bar{X_j}^{-k_4} \bar{Y_i}^{-k_5} \bar{Y_j}^{-k_6} \\ \times\sum_{\nu_1,\nu_2,\nu_3,\nu_4=0}^\infty \hat{C}_{\bf k}^{(\nu_1,\nu_2,\nu_3,\nu_4)}(\alpha_{ij}) |Y_i|^{2\nu_1} |Y_j|^{2\nu_2} |X_i|^{2\nu_3} |X_j|^{2\nu_4},\end{split}\]

where

\[\begin{split}\begin{align} X_i &=& \sqrt{\frac{2}{\Lambda_i}}(\kappa_i - \sqrt{-1}\eta_i)\approx e_ie^{\sqrt{-1}\varpi_i} + \mathcal{O}(e_i^3)\\ Y_i &=& \sqrt{\frac{1}{2\Lambda_i}}(\sigma_i - \sqrt{-1}\rho_i)\approx s_ie^{\sqrt{-1}\Omega_i}+ \mathcal{O}(s_ie_i^2) \end{align}\end{split}\]

and bars denote complex conjugates. Writing the expansion above, we’ve assumed \(k_3,k_4,k_5,k_6 \le 0\). If \(k_3>0\) then one makes the replacement \(\bar{X_i}^{-k_3}\rightarrow {X_i}^{k_3}\) and a similar replacements are made if \(k_4,k_5,\) or \(k_6>0\). Note that \(\hat{C}_{\bf k}^{(0,0,0,0)} = \tilde{C}_{\bf k}^{(0,0,0,0)}\). The relationship between the coefficients \(\hat{C}_{\bf k}^{\pmb{\nu}}\) and \(\tilde{C}_{\bf k}^{\pmb{\nu}}\) when \(\pmb{\nu}\ne0\) is more complicated.

4.1.3. Complete expansion

Often, the explicit dependence of the disturbing function on the variables \(\Lambda_i\) is ignored. This allows the resulting equiations of motion to be expressed as polynomials with fixed coefficients involving the variables \(\eta_i,\kappa_i,\rho_i,\sigma_i\) and various trignometric functions of the \(\lambda_i\) variables. At this level of approximation, Hamilton’s equations for \(\lambda_i\) become \(\dot{\lambda_i} = \frac{\partial}{\partial \Lambda_i}H\approx \frac{\partial}{\partial \Lambda_i}H_\mathrm{Kep}\). We can construct a more accurate model while retaining the polynomial nature of the equations of motion by Taylor-expanding the disturbing function coefficients \({C}_{\bf k}^{\pmb{\nu}}(\alpha_{i,j})\) with respect to the Poincare momenta \(\Lambda_i\) about some reference values \(\Lambda_{i,0}\). Defining \(\delta_i = (\Lambda_i-\Lambda_{i,0}) / \Lambda_{i,0}\), we write

\[\hat{C}_{\bf k}^{\pmb{\nu}}(\alpha_{ij}) = \sum_{l_1,l_2 = 0}^\infty {C}_{\bf k}^{\pmb{\nu},(l_1,l_2)}(\alpha_{ij,0})\delta_i^{l_1}\delta_j^{l_2}\]

The complete expansion of the interaction term between planets \(i\) and \(j\) can therefore be written as

\[\begin{split}\begin{multline} -\frac{Gm_im_j}{|\mathbf{r}_i - \mathbf{r}_j|} + \frac{\tilde{\mathbf{r}}_i \cdot \tilde{\mathbf{r}}_j}{m_*} = \\ -\frac{Gm_im_j}{a_{j,0}}\sum_{\mathbf{k}}\sum_{\nu_1,...,\nu_4 = 0}^\infty \sum_{l_1,l_2=0}^\infty C_{\bf{k}}^{\pmb{\nu},\pmb{l}}(\alpha_{ij,0}) |Y_i|^{|k_5|+2\nu_1} |Y_j|^{|k_6|+2\nu_2} |X_i|^{|k_3|+2\nu_3} |X_j|^{|k_4|+2\nu_4} \delta_i^{l_1}\delta_j^{l_j} \\ \times \cos(k_1\lambda_j+k_2\lambda_i+k_3\varpi_i+k_4\varpi_j+k_5\Omega_i+k_6\Omega_j)~. \end{multline}\end{split}\]

Individual coefficients, \(C_{\bf{k}}^{\pmb{\nu},\pmb{l}}\), can be obtained as dictionaries with the function df_coefficient_C(). This is the basic function that the various methods for adding terms to a PoincareHamiltonian are built upon.

4.2. Calculating Coefficients

The functions df_coefficient_Ctilde() and df_coefficient_C() calculate the disturbing function coefficients \(\tilde{C}_{\bf k}^{\pmb{\nu}}\) and \({C}_{\bf k}^{\pmb{\nu},(l_1,l_2)}\), respectively. These functions return results as dictionaries representing sums of Laplace coefficients,

\[b_{s}^{(j)}(\alpha) = \frac{1}{\pi}\int_{0}^\pi \frac{\cos(j\theta)}{(1 + \alpha^2 - 2\alpha\cos\theta)^{s}},\]

and their derivatives with respect to \(\alpha\). Specifically, the dictionary keys represent Laplace coefficients and values represent their coefficients. For example, the disturbing function coefficient

\[\tilde{C}_{(5,-2,0,-1,-2,0)}^{(0,0,0,0)}(\alpha) = \frac{10}{4}\alpha b_{3/2}^{(3)} + \frac{1}{4}\alpha^2\frac{d}{d\alpha}b_{3/2}^{(3)}\]

is obtained with df_coefficient_Ctilde() using

k = (5,-2,0,-1,-2,0)
nu = (0,0,0,0)
Ccoeff = df_coefficient_Ctilde(*k,*nu)
print(Ccoeff)

which dispays:

{(1, (1.5, 3, 0)): 2.5, (2, (1.5, 3, 1)): 0.25, ('indirect', 1): 0}

The value of the disturbing function coefficient as a function of \(\alpha\) is given by the sum of dictionary key-value pairs \(\{(p_i,(s_i,j_i,n_i)): A_i\}\) and \(\{(\mathrm{'indirect'},p_\mathrm{ind}):A_\mathrm{ind}\}\) according to

\[\sum_{i}A_i \alpha^p\frac{d}{d\alpha}b^{(j)}_{s}(\alpha) + A_\mathrm{ind}\alpha^{-p_\mathrm{ind}/2}~.\]

These dictionaries can be evaluated at a specific semi-major axis ratio, \(\alpha\), using the function evaluate_df_coefficient_dict(). Note that individual Laplace coefficient values can also be evaluated using the laplace_b() function.

4.3. Generating Arguments

In many applications, a collection of disturbing function terms with related arguments are desired. For example, a user might want all of the arguments corresponding to a particular mean motion resonance or all of the secular terms up to a given order in some planet pair’s eccentricities. The disturbing_function module provides methods to generate lists of disturbing function arguments (i.e., \((\pmb{k},\pmb{\nu})\) combinations) related to one another.

The function df_arguments_dictionary() supplies a comprehensive list of all possible disturbing function cosine arguments appearing up to a given order, \(N_\mathrm{max}\), in planets’ eccentricities and inclinations. The functions list_secular_terms() and list_resonance_terms() can provide lists of \((\pmb{k},\pmb{\nu})\) pairs for secular or resonant disturbing function terms in a user-specified range of orders.

4.4. API

celmech.disturbing_function.df_arguments_dictionary(Nmax)[source]

Get the arguments appearing in the disturbing function up to order Nmax. Arguments are returned as a nested dictionary. The outer level keys are orders. The inner level keys are \(\Delta k = |k_1-k_2|\), denoting the order of MMR for which the argument appears. The values of the inner level dictionaries are lists of tuples containing \((k_3,k_4,k_5,k_6)\). Specifically, the dictionary is returned in the form:

{
0:{0:[(0,0,0,0)],
1:{1:[(-1,0,0,0),(0,-1,0,0)]},
…
Nmax:{
0:[(k3,k4,k5,k6)],
…
Nmax:[(k3,k4,k5,k6),…]
}
}
Parameters

Nmax (int) – Maximum order of dictionary arguments.

Returns

arguments – Nested defaultdicts containing the cosine arguments appearing in the disurbing function.

Return type

defaultdict

celmech.disturbing_function.df_coefficient_C(k1, k2, k3, k4, k5, k6, nu1, nu2, nu3, nu4, l1=0, l2=0, include_indirect_terms=True)[source]

Get the coefficient of the disturbing function term

\[|Y_i|^{|k_5|+2\nu_1} |Y_j|^{|k_6|+2\nu_2} |X_i|^{|k_3|+2\nu_3} |X_j|^{|k_4|+2\nu_4} \delta_i^{l_1} \delta_j^{l_2} \times \cos[ k_1 \lambda_j + k_2 \lambda_i + k_3 \varpi_i + k_4 \varpi_j + k_5 \Omega_i + k_6 \Omega_j ]~,\]

where the indices \(i\) and \(j\) corresponds to the inner and outer planets, respectively. The result it returned as a dictionary of Laplace coefficient arguemnts and their numerical coefficents.

4. Arguments:

k1int

Coefficient of outer planet’s mean longitude in cosine argument

k2int

Coefficient of inner planet’s mean longitude in cosine argument

k3int

Coefficient of inner planet’s mean longitude in cosine argument

k4int

Coefficient of outer planet’s mean longitude in cosine argument

k5int

Coefficient of inner planet’s longitude of ascending node in cosine argument

k6int

Coefficient of outer planet’s longitude of ascending node in cosine argument

nu1int

Select specific term where the exponent of \(|Y_i|\) is \(|k_5|+2\nu_1\)

nu2int

Select specific term where the exponent of \(|Y_j|\) is \(|k_6|+2\nu_2\)

nu3int

Select specific term where the exponent of \(|X_i|\) is \(|k_3|+2\nu_3\)

nu4int

Select specific term where the exponent of \(|X_j|\) is \(|k_4|+2\nu_4\)

l1int, optional

Select specifc term where exponent of \(\delta_1 = (\Lambda_1 - \Lambda_{1,0})/\Lambda_{1,0})\) is l1. Default value is l1 = 0

l2int, optional

Select specifc term where exponent of \(\delta_2 = (\Lambda_2 - \Lambda_{2,0})/\Lambda_{2,0})\) is l2. Default value is l2 = 0

include_indirect_termsboole, optional

whether the contribution of indirect terms should be accounted for when computing the coefficients. default is true.

returns

The coefficient is given by the sum over laplace coefficients contained in the dictionary entries:

\[\sum C \times \alpha^p \frac{d^{n}}{d\alpha^{n}} b_{s}^{j}(\alpha)\]

where the dictionary entries are in the form { (p,(s,j,n)) : C }

rtype

dictionary

celmech.disturbing_function.df_coefficient_Ctilde(k1, k2, k3, k4, k5, k6, nu1, nu2, nu3, nu4, include_indirect=True)[source]

Get the coefficient of the disturbing function term:

\[s_i^{|k_5|+2\nu_1} s_j^{|k_6|+2\nu_2} e_i^{|k_3|+2\nu_4} e_j^{|k_4|+2\nu_3} \times \cos[ k_1\lambda_j + k_2\lambda_i + k_3 \varpi_i + k_4 \varpi_j + k_5 \Omega_i + k_6 \Omega_j ]\]

where the indices \(i\) and \(j\) corresponds to the inner and outer planets, respectively. The result is returned as a dictionary of Laplace coefficient arguemnts and their numerical coefficents.

4. Arguments:

k1int

Coefficient of outer planet’s mean longitude in cosine argument

k2int

Coefficient of inner planet’s mean longitude in cosine argument

k3int

Coefficient of inner planet’s mean longitude in cosine argument

k4int

Coefficient of outer planet’s mean longitude in cosine argument

k5int

Coefficient of inner planet’s longitude of ascending node in cosine argument

k6int

Coefficient of outer planet’s longitude of ascending node in cosine argument

nu1int

Select specific term where the exponent of \(s_i\) is \(|k_5|+2\nu_1\)

nu2int

Select specific term where the exponent of \(s_j\) is \(|k_6|+2\nu_2\)

nu3int

Select specific term where the exponent of \(e_i\) is \(|k_3|+2\nu_3\)

nu4int

Select specific term where the exponent of \(e_j\) is \(|k_4|+2\nu_4\)

include_indirectbooole, optional

Whether to include the indirect contribution to the disturibing function coefficient.

returns

The coefficient is given by the sum over laplace coefficients contained in the dictionary entries:

..math::

sum A times alpha^p frac{d^{n}}{dalpha^{n}} b_{s}^{j}(alpha)

where the dictionary entries are in the form { (p,(s,j,n)) : A }

rtype

dictionary

celmech.disturbing_function.evaluate_df_coefficient_dict(coeff_dict, alpha)[source]

Evaluate a dictionary representing a sum of Laplace coefficient terms like those returned by df_coefficient_C and df_coefficient_Ctilde evaluated at a specific value of semi-major axis ratio, alpha.

Parameters

  • coeff_dict (dictionary) –

    Dictionary with entries {(p,(s,j,n)) : coeff} representing a sum of Laplace coefficients:

    \[\mathrm{coeff} \times \alpha^p \frac{ d^n b_s^{(j)}(\alpha)} { d\alpha^n}\]

  • alpha (float) – Value of semi-major axis ratio, \(\alpha=a_i/a_j\), appearing as an argument of Laplace coefficients.

Returns

The sum of Laplace coefficeint terms represented by dictionary entries.

Return type

float

celmech.disturbing_function.laplace_b(s, j, n, alpha)[source]

Calculates \(n\)-th derivative with respect to \(\alpha\) of Laplace coefficient \(b_s^j(\alpha)\). Uses recursion and scipy special functions.

Parameters

  • s (float) – half-integer parameter of Laplace coefficient.

  • j (int) – integer parameter of Laplace coefficient.

  • n (int) – Specify the \(n\).

  • alpha (float) – Semi-major axis ratio, \(\alpha=a_{in}/a_{out}\).

Returns

The value of the coefficient.

Return type

float

celmech.disturbing_function.list_resonance_terms(p, q, min_order=None, max_order=None, eccentricities=True, inclinations=True)[source]

Generate the list of disturbing function terms for a p:p-q resonance with eccentricity/inclination order between min_order and max_order

Parameters

  • p (int) – Determines resonance

  • q (int) – Order of the resonance

  • min_order (int, optional) – Minimum order in eccentricities and inclinations to include. Defaults to order of the resonance q (leading order)

  • max_order (int) – Maximum order in eccentricities and inclinations to include. Defaults to order of the resonance q (leading order)

  • eccentricities (bool, optional) – By default, includes all eccentricity terms. Can set to False to exclude any eccentricity terms (e.g., fixed circular orbits).

  • inclinations (bool, optional) – By default, includes all inclination terms. Can set to False to exclude any inclination terms (e.g., co-planar systems).

Returns

term – A list of disturbing function terms. Each entry in the list is of the form (k_vec, nu_vec) (see PoincareHamiltonian.add_cosine_term in poincare.py)

Return type

list

celmech.disturbing_function.list_secular_terms(min_order, max_order, eccentricities=True, inclinations=True)[source]

Generate the list of secular disturbing function terms with eccentricity/inclination order between min_order and max_order.

Parameters

  • min_order (int) – Minimum order in eccentricities and inclinations to include.

  • max_order (int) – Maximum order in eccentricities and inclinations to include.

  • eccentricities (bool, optional) – By default, includes all eccentricity terms. Can set to False to exclude any eccentricity terms (e.g., fixed circular orbits).

  • inclinations (bool, optional) – By default, includes all inclination terms. Can set to False to exclude any inclination terms (e.g., co-planar systems).

Returns

term – A list of disturbing function terms. Each entry in the list is of the form (k_vec, nu_vec) (see PoincareHamiltonian.add_cosine_term in poincare.py)

Return type

list