celmech: A Python package for celestial mechanics

We present celmech, an open-source Python package designed to facilitate a wide variety of celestial mechanics calculations. The package allows users to formulate and integrate equations of motion incorporating user-specified terms from the classical disturbing function expansion of the interaction potential between pairs of planets. The code can be applied, for example, to isolate the contribution of particular resonances to a system's dynamical evolution and develop simple analytical models with the minimum number of terms required to capture a particular dynamical phenomenon. Equations and expressions can be easily manipulated by leveraging the extensive symbolic mathematics capabilities of the sympy Python package. The celmech package is designed to interface seamlessly with the popular $N$-body code REBOUND to facilitate comparisons between calculation results and direct $N$-body integrations. The code is extensively documented and numerous example Jupyter notebooks illustrating its use are available online.


INTRODUCTION
Celestial mechanics is the oldest branch of physics and increasingly sophisticated analytical techniques for predicting planetary motions have been developed over the course of the subject's long history. While precise orbital solutions are now typically found through direct numerical integrations of the equations of motion, analytic methods remain central to solving open problems in dynamics. We are often more interested in developing theoretical insight into the qualitative elements (e.g., mean motion resonances, secular evolution) responsible for dynamical phenomena such as chaos, transit timing variations, extreme orbital eccentricities etc. This is especially true when studying exoplanetary systems, where precise mass and orbital determinations are often unavailable, and analytical intuition can guide the construction and interpretation of suites of N -body integrations across a relevant subset of the vast parameter space.
The fact that planetary masses are much smaller than the central stellar mass, and that orbital eccentricities and inclinations are typically small, often makes it possible to model planetary dynamics using only a small number of slowly varying terms in the disturbing function expansion of the gravitational interaction potential between planet pairs (a Fourier series in the angular orbital elements and a power series in eccentricities and inclinations). Nevertheless, this long-standing approach can be tedious when formulating even relatively simple models as it entails identifying appropriate disturbing function terms, finding expressions for their coefficients, and then evaluating them by means of a variety of special functions. Examining any of the various high-order expansion disturbing function expansions published during the 19th century (e.g., Peirce 1849; Le Verrier 1855; Boquet 1889), one quickly appreciates that formulating any sophisticated, high-order analytic theories by hand represents a monumental undertaking.
Since the sixties, such detailed work has been performed using specialized codes dedicated to computing and manipulating disturbing function expansions (see references in Henrard 1989;Brumberg 1995). The state-of-the-art TRIP code (Gastineau & Laskar 2011, a dedicated computer algebra system for celestial mechanics, has proven especially valuable for studies the solar system's long-term dynamics (Hoang et al. 2021(Hoang et al. , 2022Mogavero & Laskar 2022).
However, such analyses have largely been limited to the specialized groups with the expertise to develop these highly technical routines. 1 By contrast, the last decade has clearly demonstrated the impact open-source scientific software can have on the field (Tollerud et al. 2019), both by opening up specialized calculations to a wide base of users, and by providing application programming interfaces (APIs) that enable the interconnection of separate codes in a broad range of new and creative ways. Particularly successful cases in the exoplanet field include, e.g., ttvfast (Deck et al. 2014), EXOFAST (Eastman et al. 2013), batman (Kreidberg 2015), radvel (Fulton et al. 2018), REBOUND (Rein & Liu 2012), REBOUNDx (Tamayo et al. 2020), exostriker (Trifonov 2019), NbodyGradient , and exoplanet (Foreman-Mackey et al. 2021). Ellis & Murray (2000) provided an important early step towards making computer-assisted disturbing function expansions more accessible to the wider research community. They derived a method for calculating coefficients associated with specific Fourier terms in the disturbing function. This offers an efficient alternative to computing complete high-order disturbing function expansions by means of dedicated computer algebra. Mathematica notebooks implementing their algorithm were made available as part of the online material accompanying the popular textbook Murray & Dermott (1999).
In this paper, we present celmech, an open-source Python package, which facilitates the manipulation of disturbing function expansions to any order in the eccentricities and inclinations, and enables a broad range of analytical investigations. celmech is built upon on the disturbing function expansion algorithm derived by Ellis & Murray (2000), though with important modifications described in Section 4 that allow the equations of motion to be formulated in a powerful Hamiltonian framework instead of in orbital elements. We have designed the code's API to easily interface with the REBOUND package, streamlining comparisons between semi-analytic celmech models and direct N -body integrations. Additionally, celmech includes functionality for performing canonical transformations, which are useful for both extracting conserved quantities and reducing the dimensionality of various problems (Sec. 5.3). We also plan to incorporate additional modules in upcoming releases.
The plan of the paper is as follows: we present a simple example application in Section 2 with the goal of introducing readers to the code and illustrating a typical workflow before proceeding through more in-depth discussions of technical details in subsequent sections. Section 3 reviews the coordinate system and Hamiltonian framework that celmech uses to formulate the equations of motion. Section 4 describes celmech's capabilities for performing disturbing function expansions in terms of both classical Keplerian orbital elements (Sec 4.1), as well as canonical variables (Sec 4.2). Section 5 introduces celmech's object-oriented approach to modeling Hamiltonian systems in general (Sec. 5.1) as well as those specific to planetary systems (Sec. 5.2). It also describes celmech's ability to manipulate Hamiltonian models via the application canonical transformations (Sec. 5.3), including Lie series transformation series that arise in Hamiltonian perturbation theory (Sec. 5.4). To aid the reader, we provide a list the main mathematical symbols appearing in the paper along with their definitions in Appendix A.

A SIMPLE EXAMPLE
We begin with a practical example to highlight some of the central machinery in celmech, demonstrate the code's API, and illustrate a possible workflow. We consider a system of three Earth-mass planets around a solar-mass star. The inner two are initialized at the 3:2 mean-motion resonance (MMR) with orbital periods of 1 and 1.5 years respectively, while the third planet's orbital period is 3.2 years and does not form any simple integer ratios with the inner two. The outermost orbit begins with zero eccentricity and inclination to the reference plane. The innermost orbit has both its eccentricity and inclination (in radians) set to 0.02, while the middle orbit has these quantities both set to 0.03. The longitudes of ascending node, and the pericenters of the inner two orbits lie along the positive xaxis, and both planets start at conjunction at their respective apocenters where their orbits are spaced furthest apart (Fig. 1). The accompanying Jupyter notebook shows how celmech conveniently interfaces with the REBOUND Nbody package, which provides extensive functionality for initializing orbital configurations that can then be converted to a celmech PoincareHamiltonian object. The PoincareHamiltonian class, detailed below in Section 5.2, will be used throughout the present example to build and manipulate a Hamiltonian model for the system's dynamics.
2.1. First-Order Approximation 1 A limited version of TRIP is available as a binary at www.imcce.fr/Equipes/ASD/trip. Initial orbital configuration in the reference x − y plane (inclinations to this plane are small 2 • ). Right panel: Comparison of the evolution of the middle planet's orbital eccentricity between a celmech model including the leading-order 3:2 MMR terms in the disturbing function between the inner two planets (orange, Equation (1)) to an N -body integration (blue).
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, which we named 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 the 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):

Mean Variables
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 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 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 Lie series (Morbidelli 2002, Sec. 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 [k(λ 2 − λ 1 )], and can be analytically summed over (Malhotra 1993, Sec. 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 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 1 chi.add_MMR_terms(p=2,q=1,l_max=1, indexIn=1, indexOut=2) 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 Fig. 2), yields excellent agreement with N -body.

Going to Higher Order
In the right panel of Fig. 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 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, though 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 higher order secular terms do not appear until 4th order in eccentricities and inclinations.

Limits
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 (Sec. 5.4) corrects 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.

Applications
Consideration of only the leading order terms is attractive for both analytical understanding and computational cost. As one incorporates more terms of progressively higher order, 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 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 . 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.

COORDINATE SYSTEM
This section details the coordinate system and dynamical variables the celmech code uses for formulating Hamiltonian equations of motion. We begin by reviewing the Hamiltonian formulation of the planetary N -body problem, considering N − 1 planets of mass m i with i = 1, ..., N − 1 orbiting a central star of mass M * . Let u u u * and u u u i be the Cartesian position vectors of the star and the planets, respectively, in the inertial barycentric frame. We start by formulating the Hamiltonian in term of canonical heliocentric coordinates, (r r r i ,r r r i ), where r r r i = u u u i − u u u * is the ith planet's heliocentric position vector andr r r i = m iu u u i is its canonically conjugate momentum (Laskar & Robutel 1995). In these coordinates, the Hamiltonian reads where µ i = miM * M * +mi , and M i = M * + m i . 3 We will perform a canonical transformation from the canonical heliocentric coordinates (r r r i ,r r r i ) to the canonical modified Delaunay variables that constitute action-angle coordinates of the integrable two-body Hamilontians, H Kep,i . In order to effect this transformation, we first define a set of 'canonical heliocentric' Keplerian orbital elements, (a i , e i , I i , λ i , i , Ω i ), which represent, respectively, the planetary orbit's semi-major axis, eccentricity, inclination, mean longitude, longitude of periapse, and longitude of ascending node. Standard orbital elements are computed from bodies' three-dimensional Cartesian positions and velocities as well as a "gravitational parameter" (e.g., Murray & Dermott 1999). For each planet, canonical heliocentric elements are computed from the Cartesian position vector r r r i , the "velocity" vectorr r r i /µ i , and the gravitational parameter GM i . Whiler r r i /µ i does not correspond to any physical velocity vector in the system, this definition of canonical heliocentric elements ensures that the modified Delaunay variables, defined below, form a canonical set of action-angle coordinates for the two-body Hamiltonian, H Kep,i . The celmech module nbody simulation utilities provides the function reb calculate orbits to compute canonical heliocentric orbital elements from a REBOUND simulation and the function reb add from elements to add particles to a REBOUND simulation by specifying these elements (Jupyter Notebook example).
The canonical modified Delaunay coordinates 4 are defined in terms of the canonical heliocentric orbital elements. The action variables are and have units of angular momentum. The angle variables canonically conjugate to the action variables Λ i , Γ i , and Q i are λ i , γ i = − i , and q i = −Ω i , respectively. In order to avoid coordinate system singularities that leave the canonical variables γ i and q i undefined when e i = 0 and I i = 0, celmech formulates equations of motion in terms canonical coordinate-momentum pairs rather than (Γ i , γ i ) and (Q i , q i ). In modified Delaunay variables, the two-body Hamiltonian, Equation (6), simply becomes By contrast, H (i,j) int. does not admit a simple expression in terms of the canonical variables. Rather, the benefit of expressing the Hamiltonian in terms of classical action angle variables is that H int. can be decomposed as a Fourier series in the canonical angle variables, enabling the application of methods of Hamiltonian perturbation theory.

EXPANSION OF THE DISTURBING FUNCTION
The celmech.disturbing function module provides tools for calculating coefficients appearing in the cosine series expansion of the interaction Hamiltonian, H (i,j) int. both in terms of both orbital elements as well as canonical variables. We begin by describing the expansion of the disturbing function in terms of orbital elements.

Expansion in Terms of Orbital Elements
In order to write the cosine series expansion, we will take a j > a i , 3 Strictly speaking, Hamiltonian (6) should also include the term −r r r 0 M * · i =0r r r i wherer r r 0 is the total momentum of the system and is conjugate to the canonical coordinate r r r 0 = u u u * (see, e.g., Equation 27 of Hernandez & Dehnen 2017). This term is necessary to derive the correct canonical equation of motion for r r r 0 but is not needed to get the correct canonical equations for r r r i with i > 0 becauser r r 0 ≡ 0. The timeevolution of r r r 0 can be calculated without explicitly integrating its equation of motion by noting that r r r . 4 Different literature sources use different names for these canonical variables. Murray & Dermott (1999) refer to these variables as "Poincaré" variables whereas Morbidelli (2002) refers to these variables as modified Delaunay variables and reserves the name "Poincaré" variables for the set of canonical variables introduced in Equation (9).
To compute the coefficientsC ν ν ν k (α) appearing in Equation (11), we first represent the interaction Hamiltonian as the sum of two parts, H ind. ) where we define the direct and indirect parts of the disturbing function as 5 and we denote the contributions of R , respectively. Given values of k and ν ν ν, the celmech code implements the algorithm described in Ellis & Murray (2000) to compute [C ν ν ν k (α)] dir. in terms of Laplace coefficients and their derivatives while the calculation of [C ν ν ν k (α)] ind. is detailed in Appendix B. The celmech.disturbing function module's df coefficient Ctilde function can be used to compute the coefficientsC ν ν ν k (α). The function takes two integer tuples, k = (k 1 , k 2 , ..., k 6 ) and nu = (ν 1 , ..., ν 4 ), as arguments and returns a representation of the coefficient as a python dictionary with key-value pairs in the form {(p 1 , (j 1 , s 1 , n 1 )) : p l , n l , and j l are integers, and s l are half-integers.
The function celmech.disturbing function module's evaluate df coefficient dict can be used to compute a numerical value of a coefficient,C ν ν ν k , 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 Jupyter example notebook DisturbingFunctionCoefficients.ipynb.

Expansion in Terms of Canonical Variables
While ourC ν ν ν k 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 auxillary variables: we introduce reference semi-major axes a i,0 about which expansion of coefficients will be centered where we use a Roman 'i' to denote the imaginary unit, i.e., i ≡ √ −1. The complex variables z i = e i e i i and ζ i = s i e iΩi can then be expressed explicitly in terms of the variables X i , Y i , and δ i according to 5 Note that our definition of the indirect disturbing function, R ind. , differs from Murray & Dermott (1999) because we formulate our equations of motion in terms of canonical heliocentric variables. Murray & Dermott (1999) formulate equations in terms of heliocentric position and velocity vectors that do not constitute canonically conjugate variable pairs. Consequently, their formulation of the equations of motion requires different indirect terms to be computed depending on whether a planet is subject to an interior or exterior perturber.
where overlines denote complex conjugates. To express the amplitude of a single disturbing function cosine term, cos(k · θ θ θ i,j ), in terms of the canonical variables introduced in Section 3, we seek an expression for coefficients C ν ν ν,l l l k such that, when k 3 , k 4 , k 5 , and k 6 > 0, the following relationship hold: (16) When k 3 < 0, z k3 i and X k3 i in Equation (16) are replaced withz −k3 i andX −k3 i , respectively, and analogous substitutions for variables (z j , ζ i , ζ j ) and (X j , Y i , Y j ) are made when any of k 4 , k 5 , or k 6 < 0. In Appendix C we work out explicit expressions for the coefficients C ν ν ν,l l l k in terms of terms of theC ν ν ν k coefficients and their derivatives. While these expressions are, in general, complicated, note that C (0,0,0,0),(0,0) k =C (0,0,0,0) k . The disturbing function module's df coefficient C function can be used to obtain disturbing function coefficients C ν ν ν,l l l k .

Generating disturbing function arguments of a given order
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 3 , k 4 , k 5 , k 6 ) that satisfy |k 3 | + |k 4 | + |k 5 | + |k 6 | ≤ N max , along with the condition k 5 + k 6 = 2n for integer n. The latter condition ensures C n n n,l l l k = 0 (see Section 4.1). These terms are further grouped by their value of 6 i=3 k i , which must be equal to −(k 1 +k 2 ) to ensure C n n n,l l l k = 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 userspecified order. In particular, given input p and q specifying the p:p−q along with maximum order, list resonance terms returns a python list of tuples in the form ((k 1 , k 2 , ..., k 6 ), (ν 1 , ..., ν 4 )) that includes all terms that satisfy (q − p)k 1 + pk 2 = 0 |k 3 | + |k 4 | + |k 5 | + |k 6 | + 2(ν 1 + ν 2 + ν 3 + ν 4 ) ≤ N max (17) in addition to the usual rules 6 i=1 k i = 0 and k 5 + k 6 = 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 1 = k 2 = 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.

Oblateness and General Relativity
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 where r = |r|,r = r/r, andẑ is the direction perpendicular to the primary's equatorial plane. In terms of orbital elements,r ·ẑ = sin(I) sin(f + − Ω) = 2s √ 1 − s 2 sin(f + − Ω) 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 < r −3 >= a −3 1 − e 2 −3/2 and < r −3 sin 2 (f + − Ω) >= 1 2 < r −3 > to find that the orbit-averaged value of the quadrupole potential is given by < φ J2 >= −J 2 GM * a R * a 2 1 − e 2 −3/2 1/2 + 3(s 4 − s 2 ) . The corresponding contribution to the Hamiltonian can be written in terms of canonical variables as Nobili & Roxburgh (1986) describe how the apsidal precession caused by general relativity (GR) can be approximated by including a potential term equal to 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.

HAMILTONIAN AND POINCARE MODULES
This section begins in Section 5.1 with a description of how the celmech code represents general Hamiltonian systems. In Section 5.2, we describe the PoincareHamiltonian class, which is capable of representing the Hamiltonian of a planetary system by the inclusion of a collection of disturbing function terms capturing inter-planet gravitational interactions. Section 5.3 describes how canonical transformations, which can be applied to modify any Hamiltonian represented by celmech, are implemented. Section 5.4 describes the implementation of the FirstOrderGeneratingFunction class that can be used to apply techniques from canonical perturbation theory to "average" away fast oscillations in the instantaneous, "osculating" orbital elements and transform to "mean" elements that better track the system's long-term evolution.

Representing Hamiltonian Systems
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. 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 functions flow func and jacobian func that take numerical values of the canonical variables as input and return numerical value of the flow given by Hamilton's equations, and its Jacobian with respect to the canonical variables. These function are used in conjunction with the scipy package's (Virtanen et al. 2020) integrate.ode class 6 to evolve the phase space state of the system.

Constructing Hamiltonians for Planetary Systems
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 q q = (λ 1 , η 1 , ρ 1 , ..., λ N −1 , η N −1 , ρ N −1 ) and momenta p p p = (Λ 1 , κ 1 , σ 1 , ..., Λ N −1 , κ N −1 , σ N −1 ). In addition to storing these canonical variables, a Poincare instance owns a collection of PoincareParticle objects through the Poicare.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 Poincare-Hamiltonian'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 = where H Kep,i is given by Equation (10). The user can then add specific interaction terms between 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 1 , ..., k 6 ) along with inner and outer planet indices, i and j, as input and adds the terms to the Hamiltonian. 7 The combinations of l l 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 l 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 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.

Canonical Transformations
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 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. 8 User's can also use one the following class methods for initializing some frequently-used canonical transformations: • from type2 generating function allows the user to specify a generating function F 2 (q q q, P P P ) producing a canonical transformation that satisfies the equations Provided an expression for the generating function, F 2 , this method will attempt to automatically determine the corresponding transformation rules by applying sympy's solve function to express the new variables (Q Q Q, P P P ) in terms of the old variables (q q q, p p p) and vice versa.
• cartesian to polar produces a canonical transformation that transforms user-specified conjugate coordinatemomentum pairs of old variables (q i , p i ) to new conjugate pairs (Q i , P i ) = (tan −1 (q i /p i ), 1 2 (q 2 i + p 2 i ) . All other conjugate pairs will remain unchanged.
• polar to cartesian produces a canonical transformation that transforms user-specified conjugate pairs, (q, i , p i ), to new conjugate coordinate-momentum pairs (Q i , P i ) = √ 2p i sin q i , √ 2p i cos q i . All other conjugate pairs will remain unchanged.
• from linear angle transformation produces the canonical transformation (q q q, p p p) → (Q Q Q, P P P ) given by where T is a user-supplied, invertible matrix.
• 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 planets' mean longitudes, λ i .
• 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 p i are the old momentum variables, the new Hamiltonian will be H = f H and p i = f p i will replace p i 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, re-scaled variables are instead • composite combines a series of canonical transformations.
While time-dependent canonical transformations are not currently supported, they can always be implemented by re-formulating 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 a canonical transformations are applied to eliminate the explicit dependence of the transformed Hamiltonian on one or more canonical variables. This is because, when the transformed Hamiltonian does not depend on one of the angle variables, its 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 , p i ) → (φ 1 , ..., φ M , Q 1 , ...., Q N −M , I 1 , ..., I M , P 1 , ..., P N −M ) where the transformed Hamiltonian, H (Q i , P i ; I j ), is independent of the coordinate variables φ j . so that the quantities I 1 , ..., 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 , P i ) 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 respresenting a point in the reduced phase space with canonical variables (Q i , P i ) 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 possess a full qp attribute that stores the full set of variables, (φ 1 , ..., φ M , Q 1 , ...., Q N −M , I 1 , ..., I M , P 1 , ..., P N −M ), while the qp attribute 9 only stores the variables of the reduced phase space, (Q i , P i ). The reduced Hamiltonian object's full qp attribute will only store the initial values of the ignorable coordinates, φ j , determined at the time the transformation was applied. It is important to recognize the the actual time evolution of these coordinates, given by d dt φ j = ∂ ∂Ij H (Q i , P i ; I j ), is generally non-trivial and depends on the specific trajectory, (Q i (t), P i (t)), of the system in the reduced phase space. Consequently, it will usually not be possible to re-construct the full phase space state of the system in the old canonical variables (q i , p i ) if only the reduced equations of motion governing the (Q i , P i ) are integrated. Nonetheless, in many applications the values of these ignorable coordinates are unimportant for the aspects the dynamics one wishes to study and 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.

Lie series transformations
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 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) mean motion resonances, where the "fast" oscillations due to the resonant terms are observable as transit timing variations (for a Lie sires approach to TTVs, see Nesvorný & Morbidelli 2008).
A powerful approach for transforming back and forth between instantaneous, or "osculating", canonical variables (which are subject to a multitude of fast oscillations) and "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, non-integrability, 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 Kep = N −1 i=1 H Kep,i , a term H slow contains 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 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, p) → (q , p ) such that the transformed Hamiltonian is H ≡ H • T −1 = H Kep + H slow + O((m/M * ) 2 ). We use the Lie series method to construct such a transformation. The Lie transformation, exp[L χ ], of a function of phase space coordinates, f , is defined in terms of the Lie derivative operator, L χ ≡ [·, χ], as 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(k · θ θ θ i,j ) with sin(k · θ θ θ i,j )/(k 1 n j + k 2 n i ) provided k 1 n j + k 2 n i = 0. The condition that k 1 n j + k 2 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 of PoincareHamiltonian that overwrites the parent class's routines for adding cosine terms so that the trigonometric terms cos(k ·θ θ θ i,j ) are instead replaced with sin(k · θ θ θ 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. These are demonstrated in example notebooks here and here.
In addition to 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 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), though not all these studies deploy a Lie series formalism. To zeroth order in inclinations and eccentricities, the disturbing function, R ind. , given in Equation (12) can be written as where P (ψ; α) = (1 + α 2 − 2α cos(λ j − λ i )) −1/2 . Therefore, a Lie series transformation with generating function whereP (α) = 1 2π π −π P (ψ, α)dψ 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 where K is the complete elliptic integral of the first kind and F (·|m) is an incomplete elliptic integrals of the first kind with modulus m. A term of the form given in Eqauation (27) can be added to the symbolic generating function stored by a FirstOrderGeneratingFunction object with the add zeroth order term method.

SUMMARY
In this paper we presented celmech, an open-source python package for conducting celestial mechanics calculations. The code combines a disturbing function expansion algorithm based on a method originally devised Ellis & Murray (2000) with the symbolic mathematics capabilities provided by the sympy library (Meurer et al. 2017) to enable users to efficiently create and integrate Hamiltonian models for planetary systems. The celmech API is designed so that users can easily compare these Hamiltonian models to direct N -body integrations with the REBOUND code (Rein & Liu 2012).
The code is available through the Python Package Index (PyPI) as well as on GitHub at github.com/shadden/celmech. The code can can be freely redistributed under the GPL-3.0 license. Online documentation can be found at celmech.readthedocs.io. A library of Jupyter notebook examples illustrating various features and applications of the code is available at github.com/shadden/celmech/tree/master/jupyter examples.
The celmech code's modular design will allow ongoing development. Some features currently implemented in celmech, in addition to the modules and capabilities detailed in this paper, include tools for formulating and integrating equations of motion governing MMRs that are derived from numerically averaged disturbing functions (rather than truncated power series expansions in eccentricity and inclination) and tools for formulating and integrating secular equations of motion. We plan to detail additional features in forthcoming papers.
Software: matplotlib (Hunter 2007), NumPy (Oliphant 2006  αi,j = ai/aj Semi-major axis ratio θ θ θi,j = (λj, λi, i, j , Ωi, Ωj) Vector of angle variables C ν ν ν k (α) Disturbing function coefficient of cos(k k k · θ θ θi,j) for expansion in ei and si b ind. = − a j GM * m i m jr r r i ·r r r j , closely following Laskar & Robutel (1995), who present a method for deriving a complete expansion of the indirect disturbing function, best accomplished with the aid of computer algebra. We extend their derivation in order to provide explicit formulas for disturbing function coefficients associated with specific cosine arguments. Writingr r r i = µ i n i a i v v v i , the indirect disturbing function can be written as Below we will treat the disturbing function terms arising from the expressions labeled 'Term 1' and 'Term 2' in Equation (B4) separately. But before doing so, we express Z i explicitly in terms of the mean longitude, λ i , rather than the true anomaly, f i , using Hansen coefficients (e.g., Hughes 1981) where the coefficients, N 0,1 k,σ , appearing in the series expansion of X 0,1 k (e) can be calculated with the aid of recursion relations described in Hughes (1981) and Murray & Dermott (1999). 10 We define the related set of coefficients X 0,1 k (e) = (1 − e 2 ) −1/2 X 0,1 k (e) so that we may write after noting that no k = 0 term appears in the sum because X 0,1 0 (e) = −e. We also define coefficients N 0,1 k,p = p n=0 −1/2 n (−1) n N 0,1 k,p−n , so that we may write X 0,1 k (e) = e |k−1| ∞ p=0 N 0,1 k,p e 2p .