celmech: A Python Package for Celestial Mechanics

The PoincareHamiltonian class has several functions for adding specific terms from the disturbing function between pairs of planets in a system that can be used build up a dynamical model. We expect that the 3:2 MMR between the inner two planets will be important, so we call the following method of our PoincareHamiltonian object, Hp:

• Hp.add_MMR_terms(p=3, q=1, indexIn=1, indexOut=2).

This adds the lowest-order disturbing function terms in eccentricities and inclinations that are associated with a 3:2 MMR between planets 1 and 2 (the star has index 0) to the Hamiltonian represented by the PoincareHamiltonian object Hp. More generally, the arguments p and q are used to specify terms associated with any p:p − q resonance, and the arguments indexIn and indexOut specify the pair of planets for which terms should be added.

We can check the terms added to the disturbing function terms by inspecting the Hp.df attribute, which prints them with their associated normalization in units of energy (this disambiguation is useful in systems with more than one pair of planets):

$$\begin{eqnarray}\begin{array}{rcl} & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(3,-2,0,-1,0,0)}^{(0,0,0,0),(0,0)}({\alpha }_{\mathrm{1,2}}){e}_{2}\cos (3{\lambda }_{2}-2{\lambda }_{1}-{\varpi }_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(3,-2,-1,0,0,0)}^{(0,0,0,0),(0,0)}({\alpha }_{\mathrm{1,2}}){e}_{1}\cos (3{\lambda }_{2}-2{\lambda }_{1}-{\varpi }_{1}),\end{array}\end{eqnarray} \tag{ 1 }$$

with the gravitational constant G, masses m_i , reference semimajor axes a_i,0 (assumed constant), and α_i,j = a_i,0/a_j,0. At first order in eccentricities (e_i ) and inclinations (I_i ), there are only two disturbing function terms associated with the 3:2 MMR with the resonant angles 3λ₂ − 2λ₁ − ϖ₁ and 3λ₂ − 2λ₁ − ϖ₂, where the λ_i are the mean longitudes, and ϖ_i are the longitudes of pericenter. These terms would typically be found by hand, e.g., in the Appendix of Murray & Dermott (1999). In the notation ⁵ of Murray & Dermott (1999), the expression would read

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{MD} & : & -\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}\left[{f}_{27}^{(3)}({\alpha }_{\mathrm{1,2}}){e}_{1}\cos (3{\lambda }_{2}-2{\lambda }_{1}-{\varpi }_{1})\right.\\ & & +\left.{f}_{31}^{(3)}({\alpha }_{\mathrm{1,2}}){e}_{2}\cos \left(3{\lambda }_{2}-2{\lambda }_{1}-{\varpi }_{2}\right)\right].\end{array}\end{eqnarray} \tag{ 2 }$$

Our $\tilde{C}$ coefficients correspond to the f coefficients of Murray & Dermott (1999), where we choose to include some additional granularity in our $\tilde{C}$ indices for specifying higher-order terms. The index convention is explained in the following subsection (Section 2.2).

The PoincareHamiltonian object also allows us to integrate the equations of motion using only the terms one has added to the disturbing function. We show the evolution for the middle planet's orbital eccentricity in Figure 1, comparing an N-body integration using REBOUND to our simple celmech model incorporating only the leading-order resonant terms between the inner two planets.

We see that the oscillation amplitude and frequency is reproduced to within ≈10%, although this causes the two models to quickly lose phase coherence. For many applications, this level of agreement is sufficient, but we can improve the agreement by incorporating more terms at higher order in the eccentricities and inclinations.

Even going to just the next (second) order in the eccentricities and inclinations introduces a large number of additional terms. The disturbing function terms required to account for the 3:2 MMR interactions between the inner planet pair and the secular interactions between all planet pairs at second order in eccentricities and inclinations are given by

$$\begin{eqnarray}\begin{array}{rcl} & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(3,-2,0,-1,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){e}_{2}\cos (3{\lambda }_{2}-2{\lambda }_{1}-{\varpi }_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(3,-2,-1,0,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){e}_{1}\cos (3{\lambda }_{2}-2{\lambda }_{1}-{\varpi }_{1})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(6,-4,0,-2,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){e}_{2}^{2}\cos (6{\lambda }_{2}-4{\lambda }_{1}-2{\varpi }_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(6,-4,-1,-1,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){e}_{1}{e}_{2}\cos (6{\lambda }_{2}-4{\lambda }_{1}-{\varpi }_{1}-{\varpi }_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(6,-4,-2,0,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){e}_{1}^{2}\cos (6{\lambda }_{2}-4{\lambda }_{1}-2{\varpi }_{1})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(6,-4,0,0,0,-2)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){s}_{2}^{2}\cos (6{\lambda }_{2}-4{\lambda }_{1}-2{{\rm{\Omega }}}_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(6,-4,0,0,-1,-1)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){s}_{1}{s}_{2}\cos (6{\lambda }_{2}-4{\lambda }_{1}-{{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(6,-4,0,0,-2,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){s}_{1}^{2}\cos (6{\lambda }_{2}-4{\lambda }_{1}-2{{\rm{\Omega }}}_{1})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,0,0,1)}({\alpha }_{\mathrm{1,2}}){e}_{2}^{2}-\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,0,1,0)}({\alpha }_{\mathrm{1,2}}){e}_{1}^{2}\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,1,0,0)}({\alpha }_{\mathrm{1,2}}){s}_{2}^{2}-\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(1,0,0,0)}({\alpha }_{\mathrm{1,2}}){s}_{1}^{2}\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(0,0,-1,1,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){e}_{1}{e}_{2}\cos (-{\varpi }_{1}+{\varpi }_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{\mathrm{2,0}}}{\tilde{C}}_{(0,0,0,0,-1,1)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,2}}){s}_{1}{s}_{2}\cos (-{{\rm{\Omega }}}_{1}+{{\rm{\Omega }}}_{2})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,0,0,1)}({\alpha }_{\mathrm{1,3}}){e}_{3}^{2}-\displaystyle \frac{{{Gm}}_{1}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,0,1,0)}({\alpha }_{\mathrm{1,3}}){e}_{1}^{2}\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,1,0,0)}({\alpha }_{\mathrm{1,3}}){s}_{3}^{2}-\displaystyle \frac{{{Gm}}_{1}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(1,0,0,0)}({\alpha }_{\mathrm{1,3}}){s}_{1}^{2}\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,-1,1,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,3}}){e}_{1}{e}_{3}\cos (-{\varpi }_{1}+{\varpi }_{3})\\ & - & \displaystyle \frac{{{Gm}}_{1}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,-1,1)}^{(0,0,0,0)}({\alpha }_{\mathrm{1,3}}){s}_{1}{s}_{3}\cos (-{{\rm{\Omega }}}_{1}+{{\rm{\Omega }}}_{3})\\ & - & \displaystyle \frac{{{Gm}}_{2}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,0,0,1)}({\alpha }_{\mathrm{2,3}}){e}_{3}^{2}-\displaystyle \frac{{{Gm}}_{2}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,0,1,0)}({\alpha }_{\mathrm{2,3}}){e}_{2}^{2}\\ & - & \displaystyle \frac{{{Gm}}_{2}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(0,1,0,0)}({\alpha }_{\mathrm{2,3}}){s}_{3}^{2}-\displaystyle \frac{{{Gm}}_{2}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,0,0)}^{(1,0,0,0)}({\alpha }_{\mathrm{2,3}}){s}_{2}^{2}\\ & - & \displaystyle \frac{{{Gm}}_{2}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,-1,1,0,0)}^{(0,0,0,0)}({\alpha }_{\mathrm{2,3}}){e}_{2}{e}_{3}\cos (-{\varpi }_{2}+{\varpi }_{3})\\ & - & \displaystyle \frac{{{Gm}}_{2}{m}_{3}}{{a}_{\mathrm{3,0}}}{\tilde{C}}_{(0,0,0,0,-1,1)}^{(0,0,0,0)}({\alpha }_{\mathrm{2,3}}){s}_{2}{s}_{3}\cos (-{{\rm{\Omega }}}_{2}+{{\rm{\Omega }}}_{3}),\end{array}\end{eqnarray} \tag{ 3 }$$

where the inclinations are expressed through ${s}_{i}=\sin ({I}_{i}/2)$, and the Ω_i are the longitudes of ascending node. We have gone from just 2 terms in Equation (1) to 26 terms in Equation (3). While the coefficients of such an expansion would be extremely tedious to look up by hand, these terms are easy to incorporate with celmech.

As can be verified in Equation (3), the $\tilde{C}$ subscript indices specify the coefficients of the cosine arguments

$$\begin{eqnarray}&&\begin{array}{l}({k}_{1},{k}_{2},{k}_{3},{k}_{4},{k}_{5},{k}_{6})\longleftrightarrow \cos \left({k}_{1}{\lambda }_{\mathrm{out}}+{k}_{2}{\lambda }_{\mathrm{in}}+{k}_{3}{\varpi }_{\mathrm{in}}\right.\\ \,\,+\,\left.{k}_{4}{\varpi }_{\mathrm{out}}+{k}_{5}{{\rm{\Omega }}}_{\mathrm{in}}+{k}_{6}{{\rm{\Omega }}}_{\mathrm{out}}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

where the "in" and "out" subscripts refer to the inner and outer body, respectively, and the ordering of the k coefficients has been chosen to match the ordering conventionally used in the cosine arguments (e.g., Murray & Dermott 1999). The amplitude of each cosine term is a power series in the eccentricities and inclinations. We choose to separate the terms in these expansions, and specify them by the $\tilde{C}$ superscript indices, as detailed in Section 4.2. Most terms in Equation (3) are the leading-order contributions with superscripts (0, 0, 0, 0), although we see examples of higher-order contributions in the secular terms without a cosine argument (all k_i = 0 in Equation (4)).

The most important second-order terms for our setup are the second-order 3:2 MMR terms, i.e., the 6:4 cosine terms above, involving the combination 6λ₂ − 4λ₁. In celmech, we would incorporate MMR terms up to a specified maximum order by adding the max_order keyword to our call above:

• Hp.add_MMR_terms(p=3, q=1, max_order=2, indexIn=1, indexOut=2).

The secular terms (terms above that do not depend on the mean longitudes λ, including those with no cosine factors) also make a small but noticeable impact. For each pair of planet indices (idx1, idx2), the secular terms can be added by

• Hp.add_secular_terms(max_order=2, indexIn=idx1, indexOut=idx2).

Comparing Figure 1 with the left panel of Figure 2, we see that this second-order model (orange curve) does improve the agreement with N-body.

Zoom In Zoom Out Reset image size

While it is tempting to simply continue adding higher-order terms, in this case it would not help. There is a separate issue causing a discrepancy with the N-body integrations, which can be fixed (at negligible computational cost) by converting from osculating to mean variables. Any time one takes only a subset of terms in the disturbing function, one is considering a (slightly) different problem from the full one. This distinction is particularly important, given our rather special initial conditions.

Each time an inner planet overtakes its outer neighbor at a conjunction, it gains orbital energy on approach as it is being pulled forward by its outer neighbor, only to lose that energy approximately symmetrically on the way out. These encounters cause short-period spikes in the orbits' semimajor axes, and the fact that we started our planets at a conjunction, means that we are starting at the peak of one of those spikes. In other words, our initial semimajor axes are, in a sense, spurious, and not reflective of the mean values they have over the vast majority of the orbit. This in turn impacts the period ratio between the inner two planets and their resonant dynamics.

Rather than incorporating even more terms in our disturbing function, the most elegant way of handling these (and other) ignored terms is through a near-identity canonical transformation from "osculating" to "mean" variables using the Lie series method (Morbidelli 2002; our Section 5.4). One can do this in celmech by creating a FirstOrderGeneratingFunction object, and specifying the terms one wants to average out. The terms in the disturbing function causing the semimajor axis spikes discussed above are of the form $\cos \left[k({\lambda }_{2}-{\lambda }_{1})\right]$, and can be analytically summed over (Malhotra 1993; our Section 5.4). These terms are independent of the eccentricities and inclinations, and thus of zeroth order, and can be added in celmech starting from a set of poincare_variables specifying the (osculating) initial conditions via

• chi=FirstOrderGeneratingFunction(poincare_variables)
• chi.add_zeroth_order_term().

The effects of unmodeled, nearby MMRs are much smaller for this problem than the zeroth-order terms above, but still make a noticeable difference, and are added by calls of the form

• chi.add_MMR_terms(p=2,q=1,l_max=1, indexIn=1, indexOut=2)
• chi.add_MMR_terms(p=4,q=1,l_max=1, indexIn=1, indexOut=2).

Once we have incorporated all the terms for which we want to correct, we transform to mean variables with chi.osculating_to_mean(), and the resulting integration (green curve in left panel of Figure 2) yields excellent agreement with N-body.

In the right panel of Figure 2, we can see that, even at second order (orange), there is a significant discrepancy with N-body (blue) in the orbital inclination evolution of the middle planet. In order to approach the N-body evolution, we have to include terms up to third order in the eccentricities and inclinations. This is the order at which the terms directly coupling the inclinations to the quickly varying eccentricities first appear. These are easily obtained by setting max_order=3 in the calls above, although we spare the reader the details of the resulting 56 terms in the disturbing function. The add_secular_terms method accepts a similar max_order keyword argument, but the higher-order secular terms do not appear until fourth order in eccentricities and inclinations.

The degree of agreement between N-body and celmech models will depend on planet masses, with smaller masses generally giving better agreement. In particular, the conversion from osculating to mean variables is not perfect. The Lie series transformation (Section 5.4) corrects the unmodeled terms (e.g., conjunctions, nearby MMRs) to first order in the planetary masses. This generates second-order (and higher) mass terms that should be added to the equations of motion governing the evolution of the mean variables. For our low-mass planets, these corrections are negligible, but at higher mass ratios with the central star, these terms start to be significant (Sansottera & Libert 2019). Adding higher-order terms in the eccentricities and inclinations will help until one reaches the level at which these second-order mass terms become significant.

The consideration of only the leading-order terms is attractive for both analytical understanding and computational cost. As one incorporates more terms of progressively higher orders, celmech can become slower than direct N-body integrations. Nevertheless, there are several applications for high-order models. First, in many instances, it may be possible to deploy integration algorithms that are more efficient than celmech's default general-purpose integrator, and any such algorithms can potentially make use of the lower-level derivative and Jacobian information generated by the PoincareHamiltonian instance (i.e., Hp). This derivative and Jacobian information can also be valuable for locating equilibrium configurations such as periodic orbits via root-finding methods. High-order expansion can also be the starting point for developing accurate integrable approximations of the dynamics through canonical transformations (Hadden 2019). Finally, incremental comparisons of successively higher-order expansions with N-body integrations can also provide a quick guide to the degree of complexity required to model a particular dynamical phenomenon of interest.

In order to write the cosine series expansion, we will take a_j > a_i , define α_i,j = a_i /a_j , ${s}_{i}=\sin ({I}_{i}/2)$, θ _i,j = (λ_j , λ_i , − γ_i , − γ_j , − q_i , − q_j ), and write

$$\begin{eqnarray}\begin{array}{rcl}{H}_{\mathrm{int}}^{(i,j)} & = & -\displaystyle \frac{{{Gm}}_{i}{m}_{j}}{{a}_{j}}\sum _{{\boldsymbol{k}}}\sum _{{\nu }_{1},{\nu }_{2},{\nu }_{3},{\nu }_{4}=0}^{\infty }{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{i,j})\\ & & \times {s}_{i}^{| {k}_{5}| +2{\nu }_{1}}{s}_{j}^{| {k}_{6}| +2{\nu }_{2}}{e}_{i}^{| {k}_{3}| +2{\nu }_{3}}{e}_{j}^{| {k}_{4}| +2{\nu }_{4}}\cos ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})\end{array}\end{eqnarray} \tag{ 11 }$$

where k = (k₁, k₂,...,k₆) and ν = (ν₁,...,ν₄) are integer multiindices. Rotation and reflection symmetries of the planets' gravitational interactions dictate that ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )\ne 0$ only if ${\sum }_{l=1}^{6}{k}_{l}=0$ and k₅ + k₆ = 2n where n is an integer (e.g., Murray & Dermott 1999). For small eccentricities and inclinations, the coefficient of the cosine term $\cos ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})$ is of the order $\sim {s}_{i}^{| {k}_{5}| }{s}_{j}^{| {k}_{6}| }{e}_{i}^{| {k}_{3}| }{e}_{j}^{| {k}_{4}| }$, and the higher-order corrections are of degree 2 and greater in the planets' eccentricities and inclinations.

To compute the coefficients ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )$ appearing in Equation (11), we first represent the interaction Hamiltonian as the sum of two parts, ${H}_{\mathrm{int}}^{(i,j)}=-\tfrac{{{Gm}}_{i}{m}_{j}}{{a}_{j}}({{ \mathcal R }}_{\mathrm{dir}.}^{(i,j)}+{{ \mathcal R }}_{\mathrm{ind}.}^{(i,j)})$ where we define the direct and indirect parts of the disturbing function as ⁸

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{\mathrm{dir}.}^{(i,j)} & = & \frac{{a}_{j}}{| {{\boldsymbol{r}}}_{j}-{{\boldsymbol{r}}}_{i}| },\\ \mathrm{and}\,{{ \mathcal R }}_{\mathrm{ind}.}^{(i,j)} & = & -\frac{{a}_{j}}{{{GM}}_{* }{m}_{i}{m}_{j}}{\tilde{{\boldsymbol{r}}}}_{i}\cdot {\tilde{{\boldsymbol{r}}}}_{j},\end{array}\end{eqnarray} \tag{ 12 }$$

and we denote the contributions of ${{ \mathcal R }}_{\mathrm{dir}.}^{(i,j)}$ and ${{ \mathcal R }}_{\mathrm{ind}.}^{(i,j)}$ to the value of the coefficient ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )$ as ${[{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )]}_{\mathrm{ind}.}$ and ${[{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )]}_{\mathrm{ind}.}$, respectively. Given the values of k and ν , the celmech code implements the algorithm described in Ellis & Murray (2000) to compute ${\left[{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )\right]}_{\mathrm{dir}.}$ in terms of Laplace coefficients and their derivatives while the calculation of ${[{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )]}_{\mathrm{ind}.}$ is detailed in Appendix B. The celmech.disturbing_function module's df_coefficient_Ctilde function can be used to compute the coefficients ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )$. The function takes two integer tuples, k = (k₁, k₂,...,k₆) and nu = (ν₁,...,ν₄), as arguments and returns a python dictionary representation of the coefficient with key-value pairs in the form $\{({p}_{1},({j}_{1},{s}_{1},{n}_{1})):{A}_{1},\ldots ,({p}_{N},({j}_{N},{s}_{N},{n}_{N})):{A}_{N},(^{\prime} {\mathtt{indirect}}^{\prime} ,{p}_{\mathrm{ind}}):{A}_{\mathrm{ind}}\}$ with

$$\begin{eqnarray}&&{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )={A}_{\mathrm{ind}}{\alpha }^{-{p}_{\mathrm{ind}}/2}+\sum _{l=1}^{N}{A}_{l}{\alpha }^{{p}_{l}}\displaystyle \frac{{d}^{{n}_{l}}}{d{\alpha }^{{n}_{l}}}{b}_{{s}_{l}}^{({j}_{l})}(\alpha )\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray*}&&{b}_{s}^{(j)}(\alpha )=\displaystyle \frac{1}{\pi }{\int }_{-\pi }^{\pi }\displaystyle \frac{\cos (j\theta )}{{\left(1+{\alpha }^{2}-2\alpha \cos \theta \right)}^{s}}d\theta \,,\end{eqnarray*}$$

and p_l , n_l , and j_l are integers, and s_l is a half-integer.

The function celmech.disturbing_function module's evaluate_df_coefficient_dict can be used to compute a numerical value of a coefficient, ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}$, for a particular choice of α. The function deriv_df_coefficient computes derivatives of disturbing function coefficients with respect to α. This function takes a dictionary of the form returned by df_coefficient_Ctilde as an argument and returns a new dictionary of the same form representing the derivative of the coefficient with respect to α. Some basic calculations of disturbing function coefficients are demonstrated in the Jupyter example notebook DisturbingFunctionCoefficients.ipynb.

While our ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}$ coefficients for the direct part of the disturbing function match those of Ellis & Murray (2000) derived in terms of orbital elements, formulating Hamiltonian equations of motion requires a disturbing function expansion expressed in terms of canonical variables. To derive an expansion in terms of canonical variables, we first introduce some auxiliary variables: we introduce the reference semimajor axes a_i,0 about which the expansion of coefficients will be centered along with α_ij,0 = a_i,0/a_j,0, ${{\rm{\Lambda }}}_{i,0}={{\rm{\Lambda }}}_{i}\sqrt{{a}_{i,0}/{a}_{i}}$, and δ_i = (Λ_i − Λ_i,0)/Λ_i,0. Next, we define ${X}_{i}=\sqrt{\tfrac{1}{{{\rm{\Lambda }}}_{i,0}}}({\kappa }_{i}-i{\eta }_{i})$ and ${Y}_{i}=\tfrac{1}{2}\sqrt{\tfrac{1}{{{\rm{\Lambda }}}_{i,0}}}({\sigma }_{i}-i{\rho }_{i})$, where we use a Roman " ${\rm{i}}$ "to denote the imaginary unit, i.e., ${\rm{i}}\equiv \sqrt{-1}$. The complex variables ${z}_{i}={e}_{i}{e}^{i{\varpi }_{i}}$ and ${{\zeta }}_{i}={s}_{i}{e}^{{\rm{i}}{{\rm{\Omega }}}_{i}}$ can then be expressed explicitly in terms of the variables X_i , Y_i , and δ_i according to

$$\begin{eqnarray}&&{z}_{i}=\displaystyle \frac{1}{\sqrt{1+{\delta }_{i}}}{X}_{i}\sqrt{1-\displaystyle \frac{1}{4(1+{\delta }_{i})}{X}_{i}{\bar{X}}_{i}},\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{\zeta }_{i}=\displaystyle \frac{1}{\sqrt{1+{\delta }_{i}}}\displaystyle \frac{{Y}_{i}}{\sqrt{1-\tfrac{1}{2(1+{\delta }_{i})}{X}_{i}{\bar{X}}_{i}}}\end{eqnarray} \tag{ 15 }$$

where overlines denote complex conjugates. To express the amplitude of a single disturbing function cosine term, $\cos ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})$, in terms of the canonical variables introduced in Section 3, we seek an expression for coefficients ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},{\boldsymbol{l}}}$ such that, when k₃, k₄, k₅, and k₆ > 0, the following relationship holds:

$$\begin{eqnarray}&&\begin{array}{l}\begin{array}{l}\left(\displaystyle \frac{1}{{a}_{j}}\right){z}_{i}^{{k}_{3}}{z}_{j}^{{k}_{4}}{\zeta }_{i}^{{k}_{5}}{\zeta }_{j}^{{k}_{6}}\sum _{{\boldsymbol{\nu }}=0}^{\infty }{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}}){s}_{i}^{2{\nu }_{1}}{s}_{j}^{2{\nu }_{2}}{e}_{i}^{2{\nu }_{3}}{e}_{j}^{2{\nu }_{4}}\\ \,\,=\,\left(\displaystyle \frac{1}{{a}_{j,0}}\right){X}_{i}^{{k}_{3}}{X}_{j}^{{k}_{4}}{Y}_{i}^{{k}_{5}}{Y}_{j}^{{k}_{6}}\sum _{{\boldsymbol{l}}=0}^{\infty }\sum _{{\boldsymbol{\nu }}=0}^{\infty }{C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},{\boldsymbol{l}}}\\ \,\,\,\,\times \,({\alpha }_{{ij},0})| {Y}_{i}{| }^{2{\nu }_{1}}| {Y}_{j}{| }^{2{\nu }_{2}}| {X}_{i}{| }^{2{\nu }_{3}}| {X}_{j}{| }^{2{\nu }_{4}}{\delta }_{i}^{{l}_{1}}{\delta }_{j}^{{l}_{2}}.\end{array}\end{array}\end{eqnarray} \tag{ 16 }$$

When k₃ < 0, ${z}_{i}^{{k}_{3}}$ and ${X}_{i}^{{k}_{3}}$ in Equation (16) are replaced with ${\bar{z}}_{i}^{-{k}_{3}}$ and ${\bar{X}}_{i}^{-{k}_{3}}$, respectively; and analogous substitutions for variables (z_j , ζ_i , ζ_j ) and (X_j , Y_i , Y_j ) are made when any of k₄, k₅, or k₆ < 0. In Appendix C, we work out explicit expressions for the coefficients ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},{\boldsymbol{l}}}$ in terms of terms of the ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}$ coefficients and their derivatives. While these expressions are, in general, complicated, note that ${C}_{{\boldsymbol{k}}}^{(0,0,0,0),(0,0)}={\tilde{C}}_{{\boldsymbol{k}}}^{(0,0,0,0)}$. The disturbing_function module's df_coefficient_C function can be used to obtain disturbing function coefficients ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},{\boldsymbol{l}}}$.

In many applications, one is interested in gathering a series of particular disturbing function terms, such as those associated with a particular resonance, up to a maximum order, N_max, in eccentricities and inclinations. The disturbing_function module provides the function df_arguments_dictionary that allows the user to enumerate all cosine arguments, k · θ _ij , appearing in the disturbing function expansion up to a given order N_max in inclinations and eccentricities. In order to do so, the function enumerates tuples (k₃, k₄, k₅, k₆) that satisfy $| {k}_{3}| +| {k}_{4}| +| {k}_{5}| +| {k}_{6}| \leqslant {N}_{\max }$, along with the condition k₅ + k₆ = 2n for integer n. The latter condition ensures ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{n}},{\boldsymbol{l}}}\ne 0$ (see Section 4.1). These terms are further grouped by their value of ${\sum }_{i=3}^{6}{k}_{i}$, which must be equal to − (k₁ + k₂) to ensure ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{n}},{\boldsymbol{l}}}\ne 0$. The Jupyter notebook example DisturbingFunctionArguments.ipynb illustrates how the results returned by df_arguments_dictionary are structured.

The disturbing_function module also provides the convenience function list_resonance_terms to create a list of all k and ν combinations associated with a particular MMR that appear in the expansion of the disturbing function up to a user-specified order. In particular, given the input p and q specifying the p:p − q MMR along with a maximum order, list_resonance_terms returns a python list of tuples in the form ((k₁, k₂,...,k₆), (ν₁,...,ν₄)) that includes all terms that satisfy

$$\begin{eqnarray}&&\begin{array}{c}(q-p){k}_{1}+{{pk}}_{2}=0,\,\mathrm{and}\\ \,\,| {k}_{3}| +| {k}_{4}| +| {k}_{5}| +| {k}_{6}| +2({\nu }_{1}+{\nu }_{2}+{\nu }_{3}+{\nu }_{4})\leqslant {N}_{\max }\end{array}\end{eqnarray} \tag{ 17 }$$

in addition to the usual rules of ${\sum }_{i=1}^{6}{k}_{i}=0$ and k₅ + k₆ = 2n for integer n. Similarly, the function list_secular_terms can be used to create a list of all k and ν combinations that appear as secular terms with k₁ = k₂ = 0 in the disturbing function expansion up to a user-specified order. Both the list_resonance_terms and list_secular_terms functions are illustrated in a Jupyter notebook example.

The effects of primary oblateness and general relativistic precession can be important for the dynamical evolution of planetary systems. Here we give expressions for the orbit-averaged potentials of these effects in terms of the canonical variables used by celmech. The gravitational potential due to the quadrupole moment of an oblate central primary at a heliocentric position r is given by

$$\begin{eqnarray}&&{\phi }_{{J}_{2}}({\boldsymbol{r}})={J}_{2}\displaystyle \frac{{{GM}}_{* }{R}_{* }^{2}}{{r}^{3}}{P}_{2}(\hat{{\boldsymbol{r}}}\cdot \hat{{\boldsymbol{z}}})\end{eqnarray} \tag{ 18 }$$

where r = ∣ r ∣, $\hat{{\boldsymbol{r}}}={\boldsymbol{r}}/r$, and $\hat{{\boldsymbol{z}}}$ is the direction perpendicular to the primary's equatorial plane. In terms of orbital elements, $\hat{{\boldsymbol{r}}}\cdot \hat{{\boldsymbol{z}}}=\sin (I)\sin (f\,+\,\varpi -{\rm{\Omega }})$ $=\,2s\sqrt{1-{s}^{2}}\sin (f\,+\,\varpi -{\rm{\Omega }})$ where f is the true anomaly, and the inclination, I, is measured relative to the primary's equatorial plane. Denoting orbit-averaged values by < · > , we use $\lt {r}^{-3}\gt =\,{a}^{-3}{\left(1-{e}^{2}\right)}^{-3/2}$ and $\lt {r}^{-3}{\sin }^{2}(f+\varpi -{\rm{\Omega }})\gt =\,\tfrac{1}{2}\lt {r}^{-3}\gt$ to find that the orbit-averaged value of the quadrupole potential is given by $\lt {\phi }_{{J}_{2}}\gt =\,-{J}_{2}\tfrac{{{GM}}_{* }}{a}{\left(\tfrac{{R}_{* }}{a}\right)}^{2}{\left(1-{e}^{2}\right)}^{-3/2}\left(1/2+3({s}^{4}-{s}^{2})\right)$. The corresponding contribution to the Hamiltonian can be written in terms of canonical variables as

$$\begin{eqnarray}&&\begin{array}{l}{H}_{{J}_{2}}({\rm{\Lambda }},\eta ,\kappa ,\rho ,\sigma )=-\displaystyle \frac{{G}^{4}{M}_{* }^{4}{\mu }^{6}{J}_{2}{R}_{* }^{2}}{{{\rm{\Lambda }}}^{3}{\left({\rm{\Lambda }}-\tfrac{1}{2}({\eta }^{2}+{\kappa }^{2})\right)}^{3}}\\ \,\,\times \,\left[\displaystyle \frac{1}{2}-\displaystyle \frac{3}{4}\left(\displaystyle \frac{{\rho }^{2}+{\sigma }^{2}}{{\rm{\Lambda }}-\tfrac{1}{2}({\eta }^{2}+{\kappa }^{2})}\right)+\displaystyle \frac{3}{16}{\left(\displaystyle \frac{{\rho }^{2}+{\sigma }^{2}}{{\rm{\Lambda }}-\tfrac{1}{2}({\eta }^{2}+{\kappa }^{2})}\right)}^{2}\right].\end{array}\end{eqnarray} \tag{ 19 }$$

Nobili & Roxburgh (1986) describe how the apsidal precession caused by general relativity (GR) can be approximated by including a potential term equal to

$$\begin{eqnarray}\begin{array}{rcl}{H}_{\mathrm{GR}}(\eta ,\kappa ) & = & -{\displaystyle \int }_{-\pi }^{\pi }\displaystyle \frac{3{G}^{2}{M}_{* }^{2}}{{c}^{2}{r}^{2}}d\lambda =-\displaystyle \frac{3{G}^{2}{M}_{* }^{2}}{{c}^{2}{a}^{2}\sqrt{1-{e}^{2}}}\\ & = & -\displaystyle \frac{3{G}^{4}{M}^{4}{\mu }^{4}}{{{\rm{\Lambda }}}^{3}{c}^{2}\left({\rm{\Lambda }}-\tfrac{1}{2}\left({\eta }^{2}+{\kappa }^{2}\right)\right)}.\end{array}\end{eqnarray} \tag{ 20 }$$

Primary oblateness and GR effects can be modeled by adding terms in the form of Equations (19) and (20) to the Hamiltonian represented by a PoincareHamiltonian object, described below in Section 5.2, using the add_orbit_average_J2_terms and add_gr_potential_terms methods.

The celmech.hamiltonian module defines the PhaseSpaceState and Hamiltonian classes, which serve as the basic objects used by celmech to represent Hamiltonian dynamical systems. The Jupyter example notebook SimplePendulum.ipynb provides a basic example that applies these classes to model the motion of a simple pendulum. PhaseSpaceState instances represent points in the phase space of a Hamiltonian dynamical system and store symbolic variables along with numerical values for canonical variables as key-value pairs of an ordered dictionary under the qp attribute. The Hamiltonian class can be used to define a Hamiltonian function of the phase space coordinates and momenta of a particular PhaseSpaceState instance. A Hamiltonian instance is initialized by passing a symbolic expression for the Hamiltonian along with a PhaseSpaceState instance. The PhaseSpaceState will then be stored by the new Hamiltonian object as its state attribute. The symbolic expression for the Hamiltonian should depend on the same canonical variables stored by the PhaseSpaceState instance. The Hamiltonian represented by a Hamiltonian class can also depend on any number of symbolic parameters whose numerical values should be set at initialization by passing a dictionary that pairs parameter symbols with their numerical values. The numerical parameter values are stored as the H_params attribute and can be updated at any point after initializing a Hamiltonian instance to explore how a system's dynamics depends on parameter values.

In addition to storing a symbolic representation of the Hamiltonian function and its associated flow, a Hamiltonian object can be used to integrate the equations of motion and evolve a phase space trajectory in time. In order to do so, a Hamiltonian instance implements the functions flow_func and jacobian_func that take numerical values of the canonical variables as input and return the numerical value of the flow given by Hamilton's equations, and its Jacobian with respect to the canonical variables. These functions are used in conjunction with the scipy package's (Virtanen et al. 2020) integrate.ode class ⁹ to evolve the phase space state of the system.

The celmech.poincare module provides tools to model planetary systems' dynamics by building a Hamiltonian from the disturbing function expansion described in Section 4. To build these Hamiltonian models, the module defines the Poincare class, a special subclass of the PhaseSpaceState class, and PoincareHamiltonian, a special subclass of the Hamiltonian class. The Poincare class represents the phase space state of an N-body planetary system in terms of the canonically conjugate coordinates q = (λ₁, η₁, ρ₁,...,λ_N−1, η_N−1, ρ_N−1) and momenta p = (Λ₁, κ₁, σ₁,...,Λ_N−1, κ_N−1, σ_N−1). In addition to storing these canonical variables, a Poincare instance owns a collection of PoincareParticle objects through the Poincare.particles attribute. The individual members of this collection represent the individual bodies of the N-body system, i.e., the star and planets. Each planet PoincareParticle object records the values of its mass and six canonical variables as well as various quantities, including orbital elements, derived from these values. The API for accessing the mass and orbital-element values of individual PoincareParticle instances is designed to closely mirror that of the REBOUND particle class.

A PoincareHamiltonian is initialized by passing a Poincare object instance, which is then stored as the PoincareHamiltonian's state attribute. Upon initialization, the symbolic Hamiltonian stored by a PoincareHamiltonian object for an N-body planetary system is simply the sum of N − 1 individual Keplerian Hamiltonians. In other words, the Hamiltonian is given by $H={\sum }_{i=1}^{N-1}{H}_{\mathrm{Kep},i}$ where H_Kep,i is given by Equation (10). The user can then add specific interaction terms between the pairs of planets from the disturbing function expansion described in Section 4 using a variety of methods of the PoincareHamiltonian class. The most basic method for adding disturbing function terms is PoincareHamiltonian.add_cosine_term. This method takes the integers k = (k₁,...,k₆) along with the inner and outer planet indices, i and j, as input and adds the terms

$$\begin{eqnarray}\begin{array}{rcl} & - & \displaystyle \frac{{{Gm}}_{i}{m}_{j}}{{a}_{j,0}}\sum _{{\boldsymbol{l}},{\boldsymbol{\nu }}}{C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},{\boldsymbol{l}}}({\alpha }_{i,j})| {Y}_{i}{| }^{| {k}_{5}| +2{\nu }_{1}}| {Y}_{j}{| }^{| {k}_{6}| +2{\nu }_{2}}\\ & & \times | {X}_{i}{| }^{| {k}_{3}| +2{\nu }_{3}}| {X}_{j}{| }^{| {k}_{4}| +2{\nu }_{4}}{\delta }_{i}^{{l}_{i}}{\delta }_{j}^{{l}_{j}}\cos [{\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{{ij}}]\end{array}\end{eqnarray} \tag{ 21 }$$

to the Hamiltonian. ¹⁰ The combinations of l and ν occurring in the sum are controlled by the user through a variety of optional keyword arguments. If no keyword arguments are passed, the default behavior is to include only the lowest-order contribution from the disturbing function, i.e., a single term with l = (0, 0) and ν = (0, 0, 0, 0). The PoincareHamiltonian methods add_MMR_terms and add_secular_terms combine this basic add_cosine_term method with the functions for enumerating the particular disturbing function arguments described in Section 4.3 to add collections of resonant or secular terms to the Hamiltonian up to a user-specified order in inclinations and eccentricities.

The study of Hamiltonian systems is often greatly simplified by canonical transformations from one set of canonical variables to another that leave Hamilton's equations unchanged. This allows for the algorithmic determination of constants of motion that can sometimes be used to greatly reduce the dimensionality of the problem (e.g., Hadden 2019). The celmech CanonicalTransformation class provides a general framework for implementing canonical transformations. Instances of this class allow the user to apply canonical transformations to expressions by applying substitution rules that replace any occurrences of old canonical variables with new ones and vice versa. A CanonicalTransformation instance can also be used to generate new Hamiltonian objects by applying its transformation to an existing Hamiltonian instance's expression for the Hamiltonian, stored as the H attribute. The resulting Hamiltonian object's state attribute will be a new PhaseSpaceState object representing the new canonical variables, and its H attribute will be the transformed Hamiltonian.

Users can create canonical transformation objects by supplying sets of old and new canonical variables along with substitution rules. ¹¹ Users can also use one the following class methods for initializing some frequently used canonical transformations:

• 1.  
  from_type2_generating_function allows the user to specify a generating function F₂( q , P ) producing a canonical transformation that satisfies the equations

  $$\begin{eqnarray}&&{\boldsymbol{Q}}={{\rm{\nabla }}}_{{\boldsymbol{P}}}{F}_{2}\,;\,{\boldsymbol{p}}={{\rm{\nabla }}}_{{\boldsymbol{q}}}{F}_{2}\,.\end{eqnarray} \tag{ 22 }$$

  Provided with an expression for the generating function, F₂, this method will attempt to automatically determine the corresponding transformation rules by applying sympy's solve function to express the new variables ( Q , P ) in terms of the old variables ( q , p ) and vice versa.
• 2.  
  cartesian_to_polar produces a canonical transformation that transforms user-specified conjugate coordinate–momentum pairs of old variables (q_i , _pi ) to new conjugate pairs $({Q}_{i},{P}_{i})=\left(\mathrm{atan}2({q}_{i},{p}_{i}),\tfrac{1}{2}({q}_{i}^{2}+{p}_{i}^{2})\right)$, where atan2 denotes the two-argument arctangent function that returns a value in the range ( − π, π]). All other conjugate pairs will remain unchanged.
• 3.  
  polar_to_cartesian produces a canonical transformation that transforms user-specified conjugate pairs, (q_i , _pi ), to new conjugate coordinate–momentum pairs $({Q}_{i},{P}_{i})\,=\left(\sqrt{2{p}_{i}}\sin {q}_{i},\sqrt{2{p}_{i}}\cos {q}_{i}\right)$. All other conjugate pairs will remain unchanged.
• 4.  
  from_linear_angle_transformation produces the canonical transformation ( q , p ) → ( Q , P ) given by

  $$\begin{eqnarray}&&{\boldsymbol{Q}}=T\cdot {\boldsymbol{q}}\,;\,{\boldsymbol{P}}={\left({T}^{-1}\right)}^{{\rm{T}}}\cdot {\boldsymbol{p}}\end{eqnarray} \tag{ 23 }$$

  where T is a user-supplied, invertible matrix.
• 5.  
  from_poincare_angles_matrix takes a Poincare instance as input, along with an invertible matrix T, and produces a transformation where the new canonical coordinates are linear combinations of the planets' angular orbital elements given by Q = T · (λ₁,...,λ_N , − ϖ₁,..., − ϖ_N, − Ω₁,..., − Ω_N ), and the corresponding new conjugate momenta are given in terms of the canonical modified Delaunay variables described in Section 3 by ${\boldsymbol{P}}={\left({T}^{-1}\right)}^{{\rm{T}}}\,\cdot ({{\rm{\Lambda }}}_{1},\ldots ,{{\rm{\Lambda }}}_{N},{{\rm{\Gamma }}}_{1},\ldots ,{{\rm{\Gamma }}}_{N},{Q}_{1},\ldots ,{Q}_{N})$.
• 6.  
  Lambdas_to_delta_Lambdas takes a Poincare instance as input and implements the canonical transformation replacing the Λ_i variables with δΛ_i = Λ_i − Λ_i,0 as the momenta canonically conjugate to the planets' mean longitudes, λ_i .
• 7.  
  rescale_transformation produces a canonical transformation that simultaneously rescales the Hamiltonian and all canonical momentum variables by a common factor, f. In other words, if H is the old Hamiltonian and _pi are the old momentum variables, the new Hamiltonian will be $H^{\prime} ={fH}$, and ${p}_{i}^{\prime} ={{fp}}_{i}$ will replace _pi as momentum variables conjugate to the old coordinates q_i . The user can optionally specify a subset of canonical variable pairs (y_i , x_i ) as Cartesian pairs, which will be transformed so that the new, rescaled variables are instead $(y{{\prime} }_{i},x{{\prime} }_{i})=(\sqrt{f}{y}_{i},\sqrt{f}{x}_{i})$.
• 8.  
  composite combines a series of canonical transformations.

While time-dependent canonical transformations are not currently supported, they can always be implemented by reformulating the Hamiltonian under consideration in the extended phase space that includes time and the negative of the system's energy as an additional conjugate variable pair.

Often canonical transformations are applied to eliminate the explicit dependence of the transformed Hamiltonian on one or more canonical variable. This is because, when the transformed Hamiltonian does not depend on one of the angle variables, the corresponding conjugate momentum variable is conserved, and Hamilton's equations for the remaining degrees of freedom are simplified. To make the discussion more concrete, suppose we have a canonical transformation of the canonical variables of an N degree of freedom system (q_i , _pi ) → (ϕ₁,...,ϕ_M , Q₁, ....,Q_N−M , I₁,...,I_M , P₁,...,P_N−M ) where the transformed Hamiltonian, $H^{\prime} ({Q}_{i},{P}_{i};{I}_{j})$, is independent of the coordinate variables ϕ_j . so that the quantities I₁,...,I_M are constant. When such an elimination can be accomplished, it is frequently convenient to consider a system's dynamics in the reduced phase space comprised of the remaining N − M degrees of freedom with canonical variables (Q_i , _Pi ) and treat the M conserved quantities, I_j , as parameters of the Hamiltonian. When transforming a Hamiltonian by applying the old_to_new_hamiltonian method of a CanonicalTransformation object, the user has the option to automatically detect whether the resulting transformed Hamiltonian is independent of any of the new canonical variables and, if so, produce a Hamiltonian object with a state attribute representing a point in the reduced phase space with canonical variables (Q_i , _Pi ) and with the quantities I_j stored as parameters under the Hamiltonian object's H_params attribute. To enable the application of the inverse transformation from the new canonical variables back to old ones in this situation, the new Hamiltonian object possesses a full_qp attribute that stores the full set of variables, (ϕ₁,...,ϕ_M , Q₁, ....,Q_N−M , I₁,...,I_M , P₁,...,P_N−M ), while the qp attribute ¹² only stores the variables of the reduced phase space, (Q_i , _Pi ). The reduced Hamiltonian object's full_qp attribute will only store the initial values of the ignorable coordinates, ϕ_j , determined at the time that the transformation was applied. It is important to recognize that the actual time evolution of these coordinates, given by $\tfrac{d}{{dt}}{\phi }_{j}=\tfrac{\partial }{\partial {I}_{j}}H^{\prime} ({Q}_{i},{P}_{i};{I}_{j})$, is generally nontrivial and depends on the specific trajectory, (Q_i (t), _Pi (t)), of the system in the reduced phase space. Consequently, it will usually not be possible to reconstruct the full phase space state of the system in the old canonical variables (q_i , _pi ) if only the reduced equations of motion governing the (Q_i , _Pi ) are integrated. Nonetheless, in many applications, the values of these ignorable coordinates are unimportant for the aspects of the dynamics one wishes to study, and the application of the inverse transformation to the data stored in the full_qp attribute can provide useful information about how the system evolves in the original canonical variables. A basic illustration of these principles is provided in an example Jupyter notebook ¹³ .

Modeling the dynamics of a planetary system with a Hamiltonian that includes only a finite number of terms from a disturbing function can be justified because the infinite number of ignored terms are rapidly oscillating and have negligible average influence on the motion, at least over sufficiently short timescales. Nevertheless, these fast oscillations in the full problem can make it challenging to match the initial conditions in an averaged formulation, and such discrepancies can make comparisons between analytical and N-body models challenging, especially near regions of phase space where the dynamical behavior changes sharply, e.g., near resonance boundaries. Additionally, while typically we are interested in methods for removing these nuisance terms in the disturbing function, there are also cases where these fast oscillations are of primary interest. An important example is for transiting pairs of planets near (but not in) MMR, where the "fast" oscillations due to the resonant terms are observable as TTVs (for a Lie series approach to TTVs, see Nesvorný & Morbidelli 2008).

A powerful approach for transforming back and forth between the instantaneous, or osculating, canonical variables (which are subject to a multitude of fast oscillations) and the mean variables that remove these fast modes is through the Lie series method. This technique introduces a near-identity canonical transformation (ensuring that the mean variables remain canonical) that, by construction, can cancel any set of rapidly oscillating terms in the Hamiltonian to first order in the planet–star mass ratio. This transformation is often referred to as averaging since, at this leading order in the masses, the technique is equivalent to averaging the problem over the rapidly oscillating angles. An advantage of the Lie series method is that it explicitly constructs the resultant terms at higher orders in the masses, which are not available through an averaging approach. A thorough discussion of canonical perturbation theory using the Lie series method, including important issues related to chaos, nonintegrability, and the (non)convergence of these series, is beyond the scope of this paper and can be found in references such as Morbidelli (2002) or Lichtenberg & Lieberman (2013). We will limit our discussion to the construction of canonical transformations to first order in the planet–star mass ratios, which are equivalent to transformations back and forth between osculating and mean variables, and are implemented in celmech.

To introduce Lie-series generating functions, consider splitting the Hamiltonian of a N-body system into three pieces: a Keplerian piece given by ${H}_{\mathrm{Kep}}={\sum }_{i=1}^{N-1}{H}_{\mathrm{Kep},i}$, a term H_slow containing slowly varying disturbing function terms such as those with (near) resonant cosine arguments or secular terms, and a term H_osc that contains rapidly oscillating terms with negligible influence on the dynamics. We seek a canonical transformation that eliminates H_osc to first order in planet masses. In other words, we seek a canonical transformation $T:(q,{\boldsymbol{p}})\to ({\boldsymbol{q}}^{\prime} ,{\boldsymbol{p}}^{\prime} )$ such that the transformed Hamiltonian is $H^{\prime} \equiv H\circ {T}^{-1}={H}_{\mathrm{Kep}}+{H}_{\mathrm{slow}}+{ \mathcal O }({\left(m/{M}_{* }\right)}^{2})$. We use the Lie series method to construct such a transformation. The Lie transformation, $\exp [{{ \mathcal L }}_{\chi }]$, of a function of phase space coordinates, f, is defined in terms of the Lie derivative operator, ${{ \mathcal L }}_{\chi }\equiv [\cdot ,\chi ]$, as

$$\begin{eqnarray}&&\exp [{{ \mathcal L }}_{\chi }]f=\sum _{n=0}^{\infty }\displaystyle \frac{1}{n!}{{ \mathcal L }}_{\chi }^{n}f=f+[f,\chi ]+\displaystyle \frac{1}{2}[[f,\chi ],\chi ]+...\,.\end{eqnarray} \tag{ 24 }$$

Since a Lie transformation evolves f under the flow induced by Hamilton's equations for the "Hamiltonian" χ, it is canonical by construction. The Lie series method seeks to find a generating function χ such that the Lie transformation $\exp [{{ \mathcal L }}_{\chi }]$ produces the desired transformation, T, that eliminates the unwanted terms. If the generating function, χ, is itself of order ${ \mathcal O }(m/{M}_{* })$, then the transformed Hamiltonian is $H^{\prime} \,=\exp [{{ \mathcal L }}_{\chi }]H\approx {H}_{\mathrm{Kep}}+{H}_{\mathrm{slow}}+{H}_{\mathrm{osc}}+[{H}_{\mathrm{Kep}},\chi ]$ after dropping terms ${ \mathcal O }({\left(m/{M}_{* }\right)}^{2})$, and the generating function χ should satisfy [H_Kep, χ] = − H_osc. Noting that the Poisson bracket of H_Kep with all canonical variables except the mean longitudes vanishes and that

$$\begin{eqnarray}&&\left[{H}_{\mathrm{Kep}},\displaystyle \frac{\sin ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})}{{k}_{1}{n}_{j}+{k}_{2}{n}_{i}}\right]=-\cos ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})\end{eqnarray} \tag{ 25 }$$

where n_i = ∂H_Kep/∂Λ_i is the ith planet's mean motion, we immediately see that any disturbing function terms in H_osc involving planets i and j written in the form of Equation (21) will be eliminated to first order in the planet masses by simply including a corresponding term in χ that replaces $\cos ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})$ with $\sin ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})/({k}_{1}{n}_{j}+{k}_{2}{n}_{i})$ provided k₁ n_j + k₂ n_i ≠ 0. The condition that k₁ n_j + k₂ n_i ≠ 0 should always be satisfied since H_osc is presumed to be comprised only of rapidly oscillating terms.

The celmech class FirstOrderGeneratingFunction allows users to construct a generating function χ that eliminates a user-specified collection of disturbing function terms in the transformed Hamiltonian to first order in planet masses following the approach just described above. Terms are added to the generating function in much the same way that terms are added to a PoincareHamiltonian instance. In fact, FirstOrderGeneratingFunction is a subclass of the PoincareHamiltonian class that overwrites the parent class's routines for adding cosine terms so that the trigonometric terms $\cos ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})$ are instead replaced with $\sin ({\boldsymbol{k}}\cdot {{\boldsymbol{\theta }}}_{i,j})/({k}_{1}{n}_{j}+{k}_{2}{n}_{i})$. The class provides a symbolic representations of the generating function under the chi and N_chi attributes, similar to the PoincareHamiltonian class's H and N_H attributes. The class also provides an array of methods for transforming canonical variables and deriving perturbative solutions to the equations of motion. The Jupyter example notebook FirstOrderGeneratingFunction.ipynb illustrates the basic application of the FirstOrderGeneratingFunction to more accurately predict the secular evolution of the mean eccentricities of a pair of planets near a 3:2 MMR. The Jupyter example notebook TransitTimingVariations.ipynb demonstrates the application of Lie series transformations to analytically predicting TTVs.

In addition to the term-wise elimination of individual rapidly oscillating disturbing function terms described above, it is possible to remove all harmonics depending solely on the angular combination ψ = λ_i − λ_j to zeroth order in eccentricities and inclinations with a single transformation. A number of authors have previously exploited the fact that perturbative solutions for planets' mutual perturbation at zeroth order in eccentricities and inclinations can be expressed in a closed form without needing to resort to an expansion in Fourier harmonics (e.g., Malhotra 1993; Agol et al. 2005; Nesvorný & Vokrouhlický 2014; Hadden et al. 2019), although not all these studies deploy a Lie series formalism. To zeroth order in inclinations and eccentricities, the disturbing function, ${{ \mathcal R }}_{\mathrm{dir}.}^{(i,j)}+{{ \mathcal R }}_{\mathrm{ind}.}^{(i,j)}$, given in Equation (12) can be written as

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal R }}_{\mathrm{dir}.}^{(i,j)}+{{ \mathcal R }}_{\mathrm{ind}.}^{(i,j)}\approx P({\lambda }_{j}-{\lambda }_{i};\alpha )\\ \,\,-\,\displaystyle \frac{1}{\sqrt{\alpha }}\cos ({\lambda }_{i}-{\lambda }_{j})+{ \mathcal O }(e,{I}^{2})\end{array}\end{eqnarray} \tag{ 26 }$$

where $P(\psi ;\alpha )={\left(1+{\alpha }^{2}-2\alpha \cos ({\lambda }_{j}-{\lambda }_{i})\right)}^{-1/2}$. Therefore, a Lie series transformation with generating function

$$\begin{eqnarray}\begin{array}{rcl}\chi & = & -\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{2}({n}_{1}-{n}_{2})}\left(-\displaystyle \frac{1}{\sqrt{\alpha }}\sin ({\lambda }_{i}-{\lambda }_{j})\right.\\ & & +\left.{\displaystyle \int }_{0}^{{\lambda }_{i}-{\lambda }_{j}}\left(P(\psi ,\alpha )-\bar{P}(\alpha )\right)d\psi \right)\,,\end{array}\end{eqnarray} \tag{ 27 }$$

where $\bar{P}(\alpha )=\tfrac{1}{2\pi }{\int }_{-\pi }^{\pi }P(\psi ,\alpha )d\psi$ is the mean value of P(ψ, α), will eliminate the disturbing function terms in Equation (26). The integral in Equation (27) can be written explicitly in terms of elliptic integrals as

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{\psi }\left(P(\psi ^{\prime} ,\alpha )-\bar{P}(\alpha )\right)d\psi ^{\prime} \\ \,\,=\,\displaystyle \frac{2}{1-\alpha }F\left(\displaystyle \frac{\psi }{2}\left|-\displaystyle \frac{4\alpha }{{\left(1-\alpha \right)}^{2}}\right.\right)-\displaystyle \frac{2}{\pi }{\mathbb{K}}({\alpha }^{2})\psi \end{array}\end{eqnarray} \tag{ 28 }$$

where ${\mathbb{K}}$ is the complete elliptic integral of the first kind, and F( · ∣m) is an incomplete elliptic integral of the first kind with modulus m. A term of the form given in Equation (27) can be added to the symbolic generating function stored by a FirstOrderGeneratingFunction object with the add_zeroth_order_term method.

Let us define the new variables ${\hat{X}}_{i}\equiv {\left(1+{\delta }_{i}\right)}^{1/2}{X}_{i}$ and ${\hat{Y}}_{i}\equiv {\left(1+{\delta }_{i}\right)}^{1/2}{Y}_{i}$ so that ${z}_{i}={\hat{X}}_{i}{\left(1-\tfrac{1}{4}| {\hat{X}}_{i}{| }^{2}\right)}^{1/2}$ and ${\zeta }_{i}={\hat{Y}}_{i}{\left(1-\tfrac{1}{2}| {\hat{X}}_{i}{| }^{2}\right)}^{-1/2}$ (see Equations (14) and (15)) along with the new disturbing function expansion coefficients ${\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}})$ that satisfy

$$\begin{eqnarray}&&\begin{array}{l}\left(\displaystyle \frac{1}{{a}_{j}}\right){z}_{i}^{{k}_{3}}{z}_{j}^{{k}_{4}}{\zeta }_{i}^{{k}_{5}}{\zeta }_{j}^{{k}_{6}}\sum _{{\boldsymbol{\nu }}=0}^{\infty }{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}})| {\zeta }_{i}{| }^{2{\nu }_{1}}| {\zeta }_{j}{| }^{2{\nu }_{2}}| {z}_{i}{| }^{2{\nu }_{3}}| {z}_{j}{| }^{2{\nu }_{4}}=\\ \,\,\left(\displaystyle \frac{1}{{a}_{j}}\right){\hat{X}}_{i}^{{k}_{3}}{\hat{X}}_{j}^{{k}_{4}}{\hat{Y}}_{i}^{{k}_{5}}{\hat{Y}}_{j}^{{k}_{6}}\sum _{{\boldsymbol{\nu }}=0}^{\infty }{\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}})| {\hat{Y}}_{i}{| }^{2{\nu }_{1}}| {\hat{Y}}_{j}{| }^{2{\nu }_{2}}| {\hat{X}}_{i}{| }^{2{\nu }_{3}}| {\hat{X}}_{j}{| }^{2{\nu }_{4}}\end{array}\end{eqnarray} \tag{ C1 }$$

where we assume here and throughout this appendix that k_m ≥ 0 for m = 3,...,6. ¹⁵ In order to obtain explicit expressions for the coefficients ${\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}})$, we write the monomial ${z}_{i}^{{k}_{3}}{\zeta }_{i}^{{k}_{5}}| {\zeta }_{i}{| }^{2{\nu }_{1}}| {z}_{i}{| }^{2{\nu }_{3}}$ in terms of variables ${\hat{X}}_{i}$ and ${\hat{Y}}_{i}$ as

$$\begin{eqnarray}\begin{array}{rcl}{z}_{i}^{{k}_{3}}{\zeta }_{i}^{{k}_{5}}{s}_{i}^{2{\nu }_{1}}{e}_{i}^{2{\nu }_{3}} & = & {\hat{X}}_{i}^{{k}_{3}}{\hat{Y}}_{i}^{{k}_{5}}| {\hat{Y}}_{i}{| }^{2{\nu }_{1}}| {\hat{X}}_{i}{| }^{2{\nu }_{3}}{\left(1-\displaystyle \frac{1}{4}| {\hat{X}}_{i}{| }^{2}\right)}^{| {k}_{3}| /2+{\nu }_{3}}{\left(1-\displaystyle \frac{1}{2}| {\hat{X}}_{i}{| }^{2}\right)}^{-| {k}_{5}| /2-{\nu }_{1}}\\ & \equiv & {\hat{X}}_{i}^{{k}_{3}}{\hat{Y}}_{i}^{{k}_{5}}| {\hat{Y}}_{i}{| }^{2{\nu }_{1}}| {\hat{X}}_{i}{| }^{2{\nu }_{3}}\sum _{n=0}^{\infty }{{ \mathcal T }}_{| {k}_{3}| /2+{\nu }_{3},| {k}_{5}| /2+{\nu }_{1}}^{(n)}| {\hat{X}}_{i}{| }^{2n}\end{array}\end{eqnarray} \tag{ C2 }$$

where the coefficients

$$\begin{eqnarray*}&&{{ \mathcal T }}_{p,q}^{(n)}={\left(-\displaystyle \frac{1}{2}\right)}^{n}\sum _{l=0}^{n}\left(\displaystyle \genfrac{}{}{0em}{}{p}{l}\right)\left(\displaystyle \genfrac{}{}{0em}{}{-q}{n-l}\right){\left(\displaystyle \frac{1}{2}\right)}^{l}\end{eqnarray*}$$

are obtained by expanding the factors ${\left(1-\tfrac{1}{4}| {\hat{X}}_{i}{| }^{2}\right)}^{| {k}_{3}| /2+{\nu }_{1}}$ and ${\left(1-\tfrac{1}{2}| {\hat{X}}_{i}{| }^{2}\right)}^{-| {k}_{5}| /2-{\nu }_{3}}$ in the first line of Equation (C2) and collecting terms $\propto | {\hat{X}}_{i}{| }^{2n}$. Using Equation (C2), we can write the individual terms appearing in the sum on the left-hand side of Equation (C1) in terms of the variables ${\hat{X}}_{i},{\hat{Y}}_{i},{\hat{X}}_{j}$, and ${\hat{Y}}_{j}$ as

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}}){z}_{i}^{{k}_{3}}{z}_{j}^{{k}_{4}}{\zeta }_{i}^{{k}_{5}}{\zeta }_{j}^{{k}_{6}}{s}_{i}^{2{\nu }_{1}}{s}_{j}^{2{\nu }_{2}}{e}_{i}^{2{\nu }_{3}}{e}_{j}^{2{\nu }_{4}}=\\ \,\,{\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}}){\hat{X}}_{i}^{{k}_{3}}{\hat{X}}_{j}^{{k}_{4}}{\hat{Y}}_{i}^{{k}_{5}}{\hat{Y}}_{j}^{{k}_{6}}| {\hat{Y}}_{i}{| }^{2{\nu }_{1}}| {\hat{Y}}_{j}{| }^{2{\nu }_{2}}\sum _{n,m\,=\,0}^{\infty }{{ \mathcal T }}_{| {k}_{3}| /2+{\nu }_{3},| {k}_{5}| /2+{\nu }_{1}}^{(n)}{{ \mathcal T }}_{| {k}_{4}| /2+{\nu }_{4},| {k}_{6}| /2+{\nu }_{2}}^{(m)}| {\hat{X}}_{i}{| }^{2n}| {\hat{X}}_{j}{| }^{2m}\,.\end{array}\end{eqnarray} \tag{ C3 }$$

Collecting terms in Equation (C3) with N₃ = ν₃ + n and N₄ = ν₄ + n, we arrive at the following expression relating ${\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}})$ to the coefficients ${\tilde{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}})$:

$$\begin{eqnarray}&&{\hat{C}}_{{\boldsymbol{k}}}^{({N}_{1},{N}_{2},{N}_{3},{N}_{4})}({\alpha }_{{ij}})=\sum _{n=0}^{{N}_{3}}\sum _{m=0}^{{N}_{4}}{{ \mathcal T }}_{| {k}_{3}| /2+{\nu }_{3},| {k}_{5}| /2+{\nu }_{1}}^{(n)}{{ \mathcal T }}_{| {k}_{4}| /2+{\nu }_{4},| {k}_{6}| /2+{\nu }_{2}}^{(m)}{\tilde{C}}_{{\boldsymbol{k}}}^{({N}_{1},{N}_{2},n,m)}({\alpha }_{{ij}}).\end{eqnarray} \tag{ C4 }$$

To derive explicit expressions for the coefficients ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},({l}_{1},{l}_{2})}$ in terms of ${\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}$, we write out the δ_i and δ_j dependence of the terms appearing in the sum on the right-hand side of Equation (C4) explicitly as

$$\begin{eqnarray}&&\begin{array}{l}\left(\displaystyle \frac{1}{{a}_{j}}\right){\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij}})\,\,{\hat{X}}_{i}^{{k}_{3}}{\hat{X}}_{j}^{{k}_{4}}{\hat{Y}}_{i}^{{k}_{5}}{\hat{Y}}_{j}^{{k}_{6}}| {\hat{Y}}_{i}{| }^{2{\nu }_{1}}| {\hat{Y}}_{j}{| }^{2{\nu }_{2}}| {\hat{X}}_{i}{| }^{2{\nu }_{3}}| {\hat{X}}_{j}{| }^{2{\nu }_{4}}\\ \,\,=\,\left\{\displaystyle \frac{1}{{a}_{j,0}}{\left(1+{\delta }_{i}\right)}^{-(| {k}_{3}| +| {k}_{5}| )/2-{\nu }_{1}-{\nu }_{3}}{\left(1+{\delta }_{j}\right)}^{2-(| {k}_{4}| +| {k}_{6}| )/2-{\nu }_{2}-{\nu }_{4}}{\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}\left({\alpha }_{{ij},0}{\left(\displaystyle \frac{1+{\delta }_{i}}{1+{\delta }_{j}}\right)}^{2}\right)\right\}\\ \,\,\,\,\times \,{X}_{i}^{{k}_{3}}{X}_{j}^{{k}_{4}}{Y}_{i}^{{k}_{5}}{Y}_{j}^{{k}_{6}}| {Y}_{i}{| }^{2{\nu }_{1}}| {Y}_{j}{| }^{2{\nu }_{2}}| {X}_{i}{| }^{2{\nu }_{3}}| {X}_{j}{| }^{2{\nu }_{4}}\,.\end{array}\end{eqnarray} \tag{ C5 }$$

Accordingly, we may write

$$\begin{eqnarray}&&{C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},({l}_{1},{l}_{2})}({\alpha }_{{ij},0})=\displaystyle \frac{1}{{l}_{1}!{l}_{2}!}\displaystyle \frac{{\partial }^{{l}_{1}+{l}_{2}}}{\partial {\delta }_{i}^{{l}_{1}}\partial {\delta }_{j}^{{l}_{2}}}{\left[{\left(1+{\delta }_{i}\right)}^{-{p}_{i}/2}{\left(1+{\delta }_{j}\right)}^{-{p}_{j}/2}{\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}\left({\alpha }_{0}{\left(\displaystyle \frac{1+{\delta }_{i}}{1+{\delta }_{j}}\right)}^{2}\right)\right]}_{({\delta }_{i},{\delta }_{j})=0}\end{eqnarray} \tag{ C6 }$$

where _pi = ∣k₃∣ + ∣k₅∣ + 2ν₁ + 2ν₃ and p_j = 4 + ∣k₄∣ + ∣k₆∣ + 2ν₂ + 2ν₄. For the case l₁ = l₂ = 0, from Equation (C6), we simply have ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},(0,0)}({\alpha }_{{ij},0})={\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}({\alpha }_{{ij},0})$. More generally, an explicit expression of ${C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},({l}_{1},{l}_{2})}({\alpha }_{{ij},0})$ in terms of ${\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}$ and its derivatives can be derived by applying the product rule for derivatives along with di Bruno's formula (e.g., Johnson 2002) generalizing the chain rule to higher derivatives to Equation (C6). A tedious but straightforward calculation yields

$$\begin{eqnarray*}&&{C}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }},({l}_{1},{l}_{2})}({\alpha }_{0})={\left.\sum _{{m}_{1}=0}^{{l}_{1}}\sum _{{r}_{1}=1}^{{l}_{1}-{m}_{1}}\sum _{{m}_{2}=0}^{{l}_{2}}\sum _{{r}_{2}=1}^{{l}_{2}-{m}_{2}}{{\rm{\Psi }}}_{{l}_{1},{l}_{2},{p}_{1},{p}_{2},{m}_{1},{m}_{2},{r}_{1},{r}_{2}}{\alpha }_{0}^{({r}_{i}+{r}_{j})}\displaystyle \frac{{d}^{({r}_{1}+{r}_{2})}}{d{\alpha }^{({r}_{1}+{r}_{2})}}{\hat{C}}_{{\boldsymbol{k}}}^{{\boldsymbol{\nu }}}(\alpha )\right|}_{\alpha ={\alpha }_{0}}\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{{\rm{\Psi }}}_{{l}_{1},{l}_{2},{p}_{1},{p}_{2},{m}_{1},{m}_{2},{r}_{1},{r}_{2}}=\displaystyle \frac{1}{{l}_{1}!{l}_{2}!}\left(\displaystyle \genfrac{}{}{0em}{}{{l}_{2}}{{m}_{2}}\right)\left(\displaystyle \genfrac{}{}{0em}{}{{l}_{1}}{{m}_{1}}\right){(-2{r}_{1}-{p}_{2}/2)}_{{m}_{2}}{(-{p}_{1}/2)}_{{m}_{1}}{{ \mathcal B }}_{{l}_{1}-{m}_{1},{r}_{1}}^{(+)}{{ \mathcal B }}_{{l}_{2}-{m}_{2},{r}_{2}}^{(-)}\,,\end{eqnarray*}$$

(x)_k ≡ x(x − 1)...(x − k + 1), and

$$\begin{eqnarray*}&&{{ \mathcal B }}_{k,n}^{(\pm )}={B}_{n,k}\left(\left(\displaystyle \genfrac{}{}{0em}{}{\pm 2}{1}\right)1!,\ldots ,\left(\displaystyle \genfrac{}{}{0em}{}{\pm 2}{n-k+2}\right)(n-k+2)!\right)\end{eqnarray*}$$

with B_n,k denoting partial Bell polynomials.