Hamiltonian hydrodynamics of eccentric discs

TL;DR

This paper develops a Hamiltonian-based theoretical framework for the ideal hydrodynamics of eccentric astrophysical discs, capturing their slow nonlinear evolution and eccentric modes through a variational principle.

Contribution

It introduces a Hamiltonian formalism for eccentric disc hydrodynamics, extending from 2D to 3D cases, and identifies nonlinear eccentric modes as energy stationary solutions.

Findings

Hamiltonian formalism derived from variational principle.

Identification of nonlinear eccentric modes and their properties.

Numerical examples of eccentric modes up to limiting amplitudes.

Abstract

We show that the ideal hydrodynamics of an eccentric astrophysical disc can be derived from a variational principle. The nonlinear secular theory describes the slow evolution of a continuous set of nested elliptical orbits as a result of the pressure in a thin disc. In the artificial but widely considered case of a 2D disc, the hydrodynamic Hamiltonian is just the orbit-averaged internal energy of the disc, which can be determined from its eccentricity distribution using the geometry of the elliptical orbits. In the realistic case of a 3D disc, the Hamiltonian needs to be modified to take into account the dynamical vertical structure of the disc. The simplest solutions of the theory are uniformly precessing nonlinear eccentric modes, which make the energy stationary subject to the angular momentum being fixed. We present numerical examples of nonlinear eccentric modes up to their limiting amplitudes. Although it lacks dissipation, which is important in many astrophysical contexts, this formalism allows a simpler theoretical approach to the nonlinear dynamics of eccentric discs than that derived from stress integrals, and also connects better with established methods of celestial mechanics for cases in which the disc interacts gravitationally with one or more orbital companion.

Equations

$$x=r\cos \varphi ,y=r\sin \varphi .$$

$$r=\frac{a(1-e^2)}{1+e\cos f},$$

$$f=\varphi -\varpi$$

$$r=a(1-e\cos E),$$

$$\cos f=\frac{\cos E-e}{1-e\cos E},\sin f=\frac{\sqrt{1-e^2}\sin E}{1-e\cos E}$$

$$M=E-e\sin E,$$

$$\Lambda =\sqrt{G{M}_{1}a},\lambda =M+\varpi ,$$

$$\Gamma =\Lambda \left(1-\sqrt{1-e^2}\right),\gamma =-\varpi .$$

$$J=\frac{\partial (x,y)}{\partial (\Lambda ,\lambda )}=J^{\circ }(a)j(E),$$

$$J^{\circ }=\frac{2}{n}$$

$$j$$

$$+\frac{\sqrt{\Gamma (2\Lambda -\Gamma )}{\gamma }_{\Lambda }\sin E}{2}$$

$$=\frac{1-e(e+a{e}_a)}{\sqrt{1-e^2}}-\frac{a{e}_a\cos E}{\sqrt{1-e^2}}-ae{\varpi }_a\sin E$$

$$\frac{\mathrm{d}}{\mathrm{d}a}(a\mp ae)=0,$$

$$\frac{\partial }{\partial \lambda }(J\Sigma n)=0.$$

$$\int \Sigma \mathrm{d}A=\iint \Sigma J\mathrm{d}\Lambda \mathrm{d}\lambda =2\pi \int J\Sigma \mathrm{d}\Lambda ,$$

$$2\pi J\Sigma =M_{\Lambda },$$

$$H=\int H_{\Lambda }\mathrm{d}\Lambda ,$$

$$M_{\Lambda }\frac{\partial \Gamma }{\partial t}=-\frac{\delta H}{\delta \gamma },$$

$$M_{\Lambda }\frac{\partial \gamma }{\partial t}=\frac{\delta H}{\delta \Gamma },$$

$$\frac{\delta H}{\delta \Gamma }=\frac{\partial H_{\Lambda }}{\partial \Gamma }-\frac{\partial }{\partial \Lambda }\left(\frac{\partial H_{\Lambda }}{\partial {\Gamma }_{\Lambda }}\right),$$

$$\frac{\delta H}{\delta \gamma }=\frac{\partial H_{\Lambda }}{\partial \gamma }-\frac{\partial }{\partial \Lambda }\left(\frac{\partial H_{\Lambda }}{\partial {\gamma }_{\Lambda }}\right),$$

$$M_{\Lambda }\frac{\partial \Gamma }{\partial t}=\frac{\partial }{\partial \Lambda }\left(\frac{\partial H_{\Lambda }}{\partial {\gamma }_{\Lambda }}\right),$$

$$\frac{\partial }{\partial t}(M_{\Lambda }\Gamma )+\frac{\partial }{\partial \Lambda }\left(-\frac{\partial H_{\Lambda }}{\partial {\gamma }_{\Lambda }}\right)=0.$$

$$C=\int M_{\Lambda }\Gamma \mathrm{d}\Lambda =\int \Gamma \mathrm{d}M$$

$$\frac{\partial H_{\Lambda }}{\partial t}=\frac{\partial H_{\Lambda }}{\partial \Gamma }\frac{\partial \Gamma }{\partial t}+\frac{\partial H_{\Lambda }}{\partial \gamma }\frac{\partial \gamma }{\partial t}+\frac{\partial H_{\Lambda }}{\partial {\Gamma }_{\Lambda }}\frac{\partial {\Gamma }_{\Lambda }}{\partial t}+\frac{\partial H_{\Lambda }}{\partial {\gamma }_{\Lambda }}\frac{\partial {\gamma }_{\Lambda }}{\partial t}.$$

$$\frac{\partial H_{\Lambda }}{\partial t}+\frac{\partial }{\partial \Lambda }\left(-\frac{\partial H_{\Lambda }}{\partial {\Gamma }_{\Lambda }}\frac{\partial \Gamma }{\partial t}-\frac{\partial H_{\Lambda }}{\partial {\gamma }_{\Lambda }}\frac{\partial \gamma }{\partial t}\right)$$

$$=0.$$

$$M_a\frac{\partial e}{\partial t}=\frac{\sqrt{1-e^2}}{n{a}^{2}e}\frac{\delta H}{\delta \varpi },$$

$$M_a\frac{\partial \varpi }{\partial t}=-\frac{\sqrt{1-e^2}}{n{a}^{2}e}\frac{\delta H}{\delta e},$$

$$\frac{\delta H}{\delta e}=\frac{\partial H_a}{\partial e}-\frac{\partial }{\partial a}\left(\frac{\partial H_a}{\partial e_a}\right),$$

$$\frac{\delta H}{\delta \varpi }=\frac{\partial H_a}{\partial \varpi }-\frac{\partial }{\partial a}\left(\frac{\partial H_a}{\partial {\varpi }_a}\right).$$

$$C=\int M_{a}n{a}^2\left(1-\sqrt{1-e^2}\right)\mathrm{d}a,$$

$$H_{\Lambda }^{(2D)}=M_{\Lambda }⟨\epsilon ⟩,or{H}_a^{(2D)}=M_a⟨\epsilon ⟩,$$

Full text

††pubyear: 2018††pagerange: Hamiltonian hydrodynamics of eccentric discs–C

1 Introduction

Eccentric gaseous discs, in which the dominant motion of the fluid consists of elliptical Keplerian orbits around a central mass, have numerous applications in astrophysics. Circumbinary discs around binary black holes in galactic nuclei, or around young binary stars in protoplanetary environments, have been found in several numerical simulations to become eccentric, even if the binary orbit is circular (MacFadyen & Milosavljević, 2008; Miranda, Muñoz & Lai, 2017; Thun, Kley & Picogna, 2017). Eccentricity leads to modulation of the accretion rate and affects the properties of the disc as well as the evolution of the binary (Artymowicz & Lubow, 1996). Similar physics allows the growth of eccentricity of planets orbiting in gaps or cavities in protoplanetary discs (Papaloizou, Nelson & Masset, 2001; Goldreich & Sari, 2003; Kley & Dirksen, 2006; D’Angelo, Lubow & Bate, 2006; Rice, Armitage & Hogg, 2008; Duffell & Chiang, 2015; Teyssandier & Ogilvie, 2016, 2017; Rosotti et al., 2017; Ragusa et al., 2018). The mechanism is related to those already investigated for circumstellar discs in ‘superhump’ binary stars (Whitehurst, 1988; Lubow, 1991a, b; Smith et al., 2007; Kley, Papaloizou & Ogilvie, 2008) and for planetary rings (Goldreich & Tremaine, 1981; Borderies, Goldreich & Tremaine, 1983). Circumstellar discs with companions on sufficiently inclined orbits can undergo Kozai–Lidov cycles in which the disc becomes highly eccentric (Martin et al., 2014; Fu, Lubow & Martin, 2015a, b; Zanazzi & Lai, 2017; Lubow & Ogilvie, 2017). Eccentric discs are also increasingly discussed in connection with tidal disruption events in which a star on a nearly parabolic orbit passes too close to the black hole at the centre of a galaxy (e.g. Guillochon, Manukian & Ramirez-Ruiz, 2014; Bonnerot et al., 2016; Cao et al., 2018), or a planetesimal is disrupted around a white dwarf (Manser et al., 2016; Cauley et al., 2018; Miranda & Rafikov, 2018).

In many of these applications the eccentricities and their gradients are not small enough to be well described by the linear theories discussed by Papaloizou (2002), Goodchild & Ogilvie (2006), Ogilvie (2008) and Teyssandier & Ogilvie (2016). A nonlinear theory of eccentric discs was developed by Ogilvie (2001) and rederived from an alternative viewpoint by Ogilvie & Barker (2014). Evolutionary equations in one spatial dimension were obtained for the shape and mass distribution of the disc; however, these equations are complicated because they contain nonlinear functions of the eccentricity, eccentricity gradient and twist that must be (pre-)computed numerically by solving a set of nonlinear ordinary differential equations around each orbit and evaluating a number of integrals of stress components combined with geometrical factors. Perhaps as a result of this complexity, the nonlinear aspects of the theory have not been investigated widely in the literature, although inviscid nonlinear modes have been computed and tested in 2D simulations by Barker & Ogilvie (2016).

The purpose of this paper is to present an alternative, and in some ways much simpler, approach to the nonlinear dynamics of eccentric discs. We focus on the special case in which dissipation can be neglected; it is then possible to apply a Hamiltonian formulation, which has attractive mathematical properties and assists in the physical interpretation of the solutions. As in the case of the classical dynamics of systems of particles or rigid bodies, Hamiltonian methods are useful mainly because they clarify the mathematical structure of the problem and provide new insights.

It is well known from the linear theory that eccentricity propagates through gaseous discs predominantly by means of pressure, in the form of a dispersive wave. Viscous forces or other dissipative effects generally contribute a slow diffusion and/or damping of eccentricity. The ideal fluid theories that we discuss in this paper describe the nonlinear version of the propagation of eccentricity by means of pressure, and neglect the weaker effects of viscous forces and other dissipative effects. Self-gravity can also make an important contribution to the propagation of eccentricity, e.g. in planetary rings; it is relatively easy to include within the Hamiltonian framework but we do not consider it in this paper.

This paper is structured as follows. We first discuss the geometry of eccentric discs (Section 2) by defining a set of canonical variables and analysing the related orbital coordinate system. We then derive the evolutionary equations for an eccentric disc from a Hamiltonian (Section 3) and point out some of their mathematical properties. The known linear theory is verified in Section 4. Section 5 discusses special modal solutions of the nonlinear equations and gives numerical examples. Conclusions are given in Section 6 and some mathematical derivations appear in the appendices.

2 Geometry of eccentric discs

Let $(x,y)$ and $(r,\phi)$ be Cartesian and polar coordinates in a plane, related by

[TABLE]

The polar equation for an elliptical Keplerian orbit of semimajor axis $a$, eccentricity $e$ and longitude of periapsis $\varpi$ (‘curly pi’) is

[TABLE]

where

[TABLE]

is the true anomaly. Equivalently,

[TABLE]

where $E$ is the eccentric anomaly, which satisfies

[TABLE]

as well as Kepler’s equation

[TABLE]

where (in this Section only) $M=n(t-\tau)$ is the mean anomaly, $n=(GM_{1}/a^{3})^{1/2}$ is the mean motion, $M_{1}$ is the central mass and $\tau$ is the time of periapsis passage.

For some purposes, such as bringing out the Hamiltonian structure of the dynamics, it is convenient to use canonical variables such as the modified Delaunay variables

[TABLE]

These are action–angle variables, with one pair $(\Lambda,\lambda)$ describing the orbital motion and the other pair $(\Gamma,\gamma)$ describing the eccentricity. More specifically, $\Lambda$ identifies the orbit and determines its period and energy, while $\lambda$ (the mean longitude) is an orbital phase that increases linearly in time by $2\pi$ per orbital period. Then $\Gamma$ is a positive-definite measure of the eccentricity of the orbit, being proportional to $e^{2}$ when $e^{2}\ll 1$, while $\gamma$ determines the orientation of the ellipse.

A planar eccentric disc involves a continuous set of nested elliptical orbits. We can describe the shape of the disc by considering $e$ and $\varpi$ to be functions of $a$. The derivatives of these functions are written as $e_{a}=\mathrm{d}e/\mathrm{d}a$ and $\varpi_{a}=\mathrm{d}\varpi/\mathrm{d}a$. Equivalently, when using modified Delaunay variables, we can consider $\Gamma$ and $\gamma$ to be functions of $\Lambda$, with derivatives $\Gamma_{\Lambda}=\mathrm{d}\Gamma/\mathrm{d}\Lambda$ and $\gamma_{\Lambda}=\mathrm{d}\gamma/\mathrm{d}\Lambda$. Thus $e_{a}$ (or $\Gamma_{\Lambda}-\Gamma/\Lambda$) is a measure of the eccentricity gradient, while $\varpi_{a}$ (or $\gamma_{\Lambda}$) is a measure of the twist.

We can use $(\Lambda,\lambda)$ as a canonical orbital coordinate system covering the disc. The first coordinate $\Lambda$ labels the elliptical orbits, while around each orbit $\lambda$ ranges from [math] to $2\pi$. The Cartesian coordinates $(x,y)$ can be deduced from $(\Lambda,\lambda)$ as follows. The value of $\Lambda$ determines $a$ and hence $e$ and $\varpi$ through the functions $e(a)$ and $\varpi(a)$. The value of $\lambda$ determines the mean anomaly $M=\lambda-\varpi$. Kepler’s equation can be solved to find the eccentric anomaly $E$. We can then find $f$ from equations (5). Thus we obtain $(r,\phi)$ and $(x,y)$.

In previous work (Ogilvie, 2001; Ogilvie & Barker, 2014) we have instead used $(\lambda,\phi)$ as an orbital coordinate system, where (in this paragraph only) $\lambda=a(1-e^{2})$ is the semilatus rectum. The geometry of eccentric discs is slightly easier to analyse when $e$ and $\varpi$ are considered to be functions of $\lambda$ rather than $a$. Also $\lambda$ is directly related to the orbital angular momentum, which is often more important than the orbital energy in the context of accretion discs. However, for the non-dissipative discs described by Hamiltonian hydrodynamics, it is more convenient to label the orbits using $a$, which is directly related to the orbital period and energy, and is a material invariant in the Hamiltonian theory. This approach also works better when combined with secular gravitational interactions, which similarly leave $a$ invariant.

By differentiating the above relations, we find that the Jacobian of the canonical orbital coordinate system $(\Lambda,\lambda)$ is

[TABLE]

where

[TABLE]

is the Jacobian for a circular disc and

[TABLE]

is dimensionless and purely related to the elliptical geometry of the disc. It can be verified that $J$ remains positive throughout the orbit provided that $(e+ae_{a})^{2}+(ae\varpi_{a})^{2}<1$. This condition was given by Statler (2001) and is exactly equivalent to the criterion $(e-\lambda e_{\lambda})^{2}+(\lambda e\varpi_{\lambda})^{2}<1$ given by Ogilvie (2001) in terms of the semilatus rectum. It is the condition for the eccentricity, eccentricity gradient and twist to be sufficiently small that neighbouring orbits do not intersect. Of course $e^{2}<1$ is also required for the orbits to be closed.

In the special case of an untwisted eccentric disc, for which $\varpi_{a}=0$, the condition for marginal orbital intersection, $e+ae_{a}=\pm 1$, is equivalent to

[TABLE]

i.e. vanishing derivative of either periapsis or apoapsis distance with $a$. The orbits intersect when either the periapsis or apoapsis distance decreases with $a$.

The orbital coordinate system is not orthogonal, but fortunately the Hamiltonian approach does not require any tensor calculus. The contravariant components of the orbital velocity are $u^{\Lambda}=0$ and $u^{\lambda}=n=(GM_{1})^{2}/\Lambda^{3}$, which depends only on $\Lambda$.

Assuming that the surface density $\Sigma$ is stationary in an inertial frame on the orbital timescale, the equation of mass conservation implies

[TABLE]

Therefore $J\Sigma$ is independent of $\lambda$. Since the mass of the disc is given by the integral

[TABLE]

we see that

[TABLE]

where $M(\Lambda)$ is the mass of the disc contained within the orbit labelled by $\Lambda$, and $M_{\Lambda}=\mathrm{d}M/\mathrm{d}\Lambda$ is the one-dimensional mass density with respect to $\Lambda$. In the Hamiltonian theory $M(\Lambda)$ does not depend on time because $\Lambda$ is a material invariant.

3 Hamiltonian evolution

The secular theory of celestial mechanics describes the slow evolution of the orbits of a number of bodies around a central mass under the assumption that the mutual gravitational interactions of the bodies are relatively weak and do not involve mean-motion resonances. The secular Hamiltonian is just the mutual gravitational energy of the interacting bodies, averaged over their orbital motion. Since the Hamiltonian does not depend on the orbital longitude of any of the bodies, each orbit preserves its semimajor axis (and energy). Angular momentum is exchanged between the orbits, resulting in precession and oscillations of eccentricity (as well as inclination, if a non-planar system is considered).

The equivalent situation for a non-self-gravitating, ideal fluid disc is one in which forces due to pressure gradients are small compared to the gravity of the central mass. This condition is satisfied in a thin disc, provided that the lengthscale associated with the eccentricity gradient or twist is long compared to the thickness of the disc. The disc can then be treated as a continuous set of nested elliptical orbits that evolve slowly due to pressure. We will see that the Hamiltonian describing this evolution is just equal to the internal energy of the system, in the artificial case of a 2D disc, while in a 3D disc there is a modification due to the dynamical vertical structure. Again, the orbital energy, and therefore the semimajor axis, is a material invariant; this is true because no energy is dissipated, nor is any energy transferred to neighbouring orbits, because the energy flux density associated with a fluid pressure is always parallel to the velocity. Therefore the secular hydrodynamics of an eccentric disc amounts to determining how the shape of the disc evolves in time. We can obtain equations for the time-derivatives of the functions $e(a,t)$ and $\varpi(a,t)$, describing this evolution. In fact we will do this first in the equivalent canonical variables, obtaining evolutionary equations for $\Gamma(\Lambda,t)$ and $\gamma(\Lambda,t)$.

If the secular Hamiltonian of a fluid disc is written as

[TABLE]

with the Hamiltonian density $H_{\Lambda}$ being expressible a function of $(\Gamma,\Gamma_{\Lambda},\gamma,\gamma_{\Lambda},\Lambda)$, then the canonical form of the evolutionary equations is

[TABLE]

where the $\delta$ notation represents a functional derivative, i.e.

[TABLE]

as appears in the Euler–Lagrange equation of variational calculus. Indeed the Hamiltonian $H$ is a functional of the two functions $\Gamma(\Lambda)$ and $\gamma(\Lambda)$ that describe the shape of the eccentric disc at any instant of time.

The rotational invariance of the system means that $H_{\Lambda}$ does not depend explicitly on $\gamma$, as we shall see in detail below. Therefore equation (17) simplifies to

[TABLE]

which can be written in the conservative form

[TABLE]

With suitable boundary conditions, such that the flux $-\partial H_{\Lambda}/\partial\gamma_{\Lambda}$ vanishes at the boundaries of the disc, the total angular-momentum deficit (AMD)

[TABLE]

is conserved. The AMD is a positive-definite measure of the eccentricity of a system, being proportional to $e^{2}$ when $e^{2}\ll 1$, and is widely used in modern celestial mechanics (Laskar, 1997). It is the difference between the angular momentum of the set of elliptical orbits and that of a set of circular orbits with the same semimajor axes. The AMD is conserved because the angular momentum is conserved and $a$ is a material invariant. The transport of angular momentum, or of AMD, is related to the twisting of the eccentric disc, because the flux $-\partial H_{\Lambda}/\partial\gamma_{\Lambda}$ is found to vanish in the untwisted case $\varpi_{a}=0$.

The total Hamiltonian $H$ is also conserved. Since the Hamiltonian density does not depend explicitly on time, it evolves according to

[TABLE]

The notation here deserves some comment. On the left-hand side of this equation we are considering $H_{\Lambda}$ as a function of $(\Lambda,t)$, describing its spatial distribution with the disc and its evolution in time. On the right-hand side we are considering $H_{\Lambda}$ as a function of the geometrical–dynamical parameters $(\Gamma,\Gamma_{\Lambda},\gamma,\gamma_{\Lambda},\Lambda)$, each of which is itself a function of $(\Lambda,t)$. This equation can be written in the conservative form

[TABLE]

Therefore both $C$ and $H$ are conserved if $\partial H_{\Lambda}/\partial\Gamma_{\Lambda}$ and $\partial H_{\Lambda}/\partial\gamma_{\Lambda}$ vanish at the boundaries. These two conservation laws are instances of Noether’s Theorem relating conservation laws to continuous symmetries.

The equivalent evolutionary equations for the more familiar but non-canonical variables $e(a,t)$ and $\varpi(a,t)$ are

[TABLE]

where $M_{a}=\mathrm{d}M/\mathrm{d}a$ is the one-dimensional mass density with respect to $a$. Note that $M_{a}/M_{\Lambda}=\mathrm{d}\Lambda/\mathrm{d}a=na/2$. In this case the relevant functional derivatives are

[TABLE]

From the fact that $H_{a}$ does not depend explicitly on $\varpi$, it is similarly possible to deduce the conservation of the AMD in the form

[TABLE]

which is equivalent to equation (23).

We present a detailed argument in Appendix A for why the secular hydrodynamics a non-self-gravitating, ideal fluid disc can be derived from Hamilton’s equations in the above form. In the artificial case of a 2D disc, the Hamiltonian density is

[TABLE]

where $\varepsilon$ is the specific internal energy and the angle brackets denote an orbital average. In Appendix B we show that this can be written as

[TABLE]

where, for a perfect gas of adiabatic index111There is an unfortunate notational clash between the adiabatic index and one of the modified Delaunay variables. The context should make it clear which meaning $\gamma$ has in any equation. $\gamma$,

[TABLE]

is the dimensionless, ‘geometric’ part of the Hamiltonian density, depending only on $e$, $ae_{a}$, $ae\varpi_{a}$ and $\gamma$, while

[TABLE]

is a fixed function of $\Lambda$ (or $a$), $P^{\circ}$ being the (2D, or vertically integrated) pressure of a circular disc with the same distributions of mass and entropy. In the isothermal limit $\gamma=1$ we have instead

[TABLE]

Explicit expressions for $F^{\mathrm{(2D)}}$ are derived in Appendix C; these involve only elementary functions if $\gamma=1$ or $2$ and involve Legendre functions otherwise.

In the realistic case of a 3D disc, the Hamiltonian is modified because of the dynamical vertical structure of the disc. The Hamiltonian density is instead

[TABLE]

where $\bar{\varepsilon}$ is the (mass-weighted) vertically averaged specific internal energy and $\langle\bar{\varepsilon}\rangle$ is its orbital average. In this case

[TABLE]

with

[TABLE]

where $h(E)$ describes the variation of the dimensionless vertical scaleheight around the orbit and must in general be obtained as the solution of the second-order ordinary differential equation (ODE) (112). Again, the geometric part $F^{\mathrm{(3D)}}$ depends only on $e$, $ae_{a}$, $ae\varpi_{a}$ and $\gamma$. In the isothermal limit $\gamma=1$ we have instead

[TABLE]

In Fig. 1 the geometric parts of the 2D and 3D Hamiltonians are compared in the case of untwisted discs. The four panels show how much the dynamics of an eccentric disc is affected by its thermodynamic behaviour and, in particular, by the dynamical vertical structure of a 3D disc.

4 Verification of linear theory

When $e$, $ae_{a}$ and $ae\varpi_{a}$ are much less than unity, the geometric part of the Hamiltonian density for a 2D disc can be approximated as the quadratic function

[TABLE]

The same expression holds in the isothermal case $\gamma=1$ except that the unimportant constant term is missing. Using this approximation for the Hamiltonian in equations (26) and (27) leads, after some algebra, to the linearized evolutionary equation

[TABLE]

where $\mathcal{E}=e\,\mathrm{e}^{\mathrm{i}\varpi}$ is the complex eccentricity, while $\Sigma^{\circ}=M_{a}/2\pi a$ is the surface density of a circular disc with the same mass distribution. This equation is exactly equivalent to equation (21) of Goodchild & Ogilvie (2006) for a 2D adiabatic disc.

In the same linear regime for a 3D disc, an expansion of the ODE (112) leads to the approximations

[TABLE]

[TABLE]

This is significantly different from the 2D version and leads instead to the linearized evolutionary equation

[TABLE]

which is exactly equivalent to equation (2) of Teyssandier & Ogilvie (2016) for a 3D adiabatic disc. The last term, in particular, is a prograde contribution to the precession of the disc due to pressure.

The difference between the 2D and 3D theories is most clearly seen in the isothermal case $\gamma=1$, for which

[TABLE]

Furthermore, in this case the eccentricity gradient and twist do not affect the vertical oscillation, so $F^{\mathrm{(3D)}}$ separates cleanly into horizontal and vertical parts, where the vertical part depends only on $e$ and has the expansion

[TABLE]

It is a decreasing function of $e$ and so contributes to prograde precession, according to equation (27).

5 Nonlinear eccentric modes

An eccentric mode is a special solution in which the eccentricity distribution $e(a)$ is independent of time and the disc is untwisted and precesses uniformly at angular frequency $\omega$, such that

[TABLE]

The Hamiltonians we are interested in have the property that $\partial H_{a}/\partial\varpi_{a}=0$ in an untwisted disc with $\varpi_{a}=0$. Therefore equation (26) is automatically satisfied in an eccentric mode, while equation (27) reduces to

[TABLE]

This is a nonlinear ODE for $e(a)$ in which the precession rate $\omega$ appears as an eigenvalue. In evaluating the right-hand side of this equation, the twist $\varpi_{a}$ can be set to zero.

The modal equation (48) can also be written in the form

[TABLE]

where

[TABLE]

is the angular momentum. Note that $L+C=\int na^{2}\,\mathrm{d}M=L^{\circ}$ is a constant property of the disc, being the total angular momentum of a circular disc with the same mass distribution.

The solutions of the modal equation generally involve oscillations of the function $e(a)$, similar to waves on a string. It is conventional in celestial mechanics to regard $e$ as a non-negative quantity and $\varpi$ as being defined modulo $2\pi$ wherever $e\neq 0$. In the case of an untwisted eccentric mode, however, where $e(a)$ may have zeros on particular orbits, it is more convenient to allow $e$ to take both positive and negative values such that $e_{a}$ and $\varpi$ vary continuously. This alternative viewpoint is possible because of the invariance of the equations (and the geometry) when $(e,\varpi)$ are replaced with $(-e,\varpi+\pi)$.

Equation (49) is the Euler–Lagrange equation for the nonlinear variational problem in which we seek stationary values of $H$ subject to the constraint $L=\mathrm{constant}$ (or, equivalently, $C=\mathrm{constant}$). The boundary conditions should be such that either the eccentricity vanishes (a rigid, circular boundary) or $\partial H_{a}/\partial e_{a}=0$ (a free, vacuum boundary). The precession frequency $\omega$ appears here as a Lagrange multiplier. This equation implies that, when the shape of the disc undergoes an infinitesimal change from a modal solution, the infinitesimal changes in $H$ and $L$ are related by $\delta H=\omega\,\delta L$. We expect there to exist multiple branches of modal solutions, along each of which the shape of the mode varies continuously, as do the values of $H$, $L$ and $\omega$. Some or all of these branches will reach the limit of a circular disc in which $H\to H^{\circ}$ and $L\to L^{\circ}$, while $\omega$ tends to a constant value which is the precession rate of the mode in linear theory. Along each branch we expect $H$ and $L$ to be related by $\mathrm{d}H/\mathrm{d}L=\omega$.

When the Hamiltonian density for an untwisted disc is written as $H_{a}=H_{a}^{\circ}F$, where $H_{a}^{\circ}$ depends only on $a$, and $F$ is a dimensionless function of $e$ and $f=e+ae_{a}$, the modal equation can be written explicitly as

[TABLE]

We illustrate this theory by considering the case of a 2D disc. We first solve the simple model problem studied by Barker & Ogilvie (2016) in which a 2D isothermal disc is contained within rigid circular boundaries with a radius ratio of $2$. The mass distribution is such that $M_{a}\propto a$, which corresponds to a uniform surface density in the limit of a circular disc. We solve the second-order nonlinear ODE (48) as an eigenvalue problem for the mode frequency $\omega$, using a shooting method. The boundary conditions are $e=0$ at $a=r_{\mathrm{in}}$ and at $a=r_{\mathrm{out}}=2r_{\mathrm{in}}$. We obtain a sequence of modes with increasing numbers of nodes in the eigenfunction $e(a)$ (Fig. 2). The amplitude of the mode can be varied by choosing the value (positive, without loss of generality) of $ae_{a}$ at the inner boundary. As this value approaches its maximum possible value of $1$, the modes become highly nonlinear, with $f=e+ae_{a}$ being close to $\pm 1$ for most values of $a$, and with sharp transitions in between. This means that the modal structures consist of alternating regions in which either the periapses or the apoapses are tightly bunched (Fig. 3).

The variation of the precession rate $\omega$ with the mode amplitude is shown in Fig. 4, along with the variation of the Hamiltonian with the angular-momentum deficit. The precession is retrograde and increases with both amplitude and mode number. The relation $\mathrm{d}H/\mathrm{d}C=-\omega$ (i.e. $\mathrm{d}H/\mathrm{d}L=+\omega$) is verified numerically for each branch. The results for the lowest-order mode agree with those of Barker & Ogilvie (2016), which were obtained by solving a nonlinear ODE for the eccentricity as a function of the semilatus rectum.

When $e+ae_{a}=\pm 1$ we have $ae=\pm a+\mathrm{constant}$. The limiting form of the lowest-order mode is given by

[TABLE]

where $\bar{r}=(r_{\mathrm{in}}+r_{\mathrm{out}})/2$. It has a finite AMD (equal to $0.02318\,Ml_{\mathrm{in}}$, where $l_{\mathrm{in}}=\sqrt{GM_{1}r_{\mathrm{in}}}$, in the case we have considered), which is the largest AMD that the system can support, and an unbounded Hamiltonian because of the logarithmic divergence (equation 136) as $f\to 1$. The limiting forms of the higher-order modes, which were not considered by Barker & Ogilvie (2016), have similar piecewise linear structures, but the positions of the break points cannot be determined in a straightforward way.

We next consider a more realistic problem involving an extended 2D disc with $a_{\mathrm{out}}/a_{\mathrm{in}}=100$ and free boundaries. We also set $\gamma=2$ for this illustrative calculation, and in fact consider a polytrope of index $1$. (As noted above, the geometric Hamiltonian of a 2D disc can be obtained in terms of elementary functions in the cases $\gamma=1$ and $\gamma=2$, which are expected to bracket the range of behaviour expected in a real system.) We take $M_{a}/a\propto\varepsilon^{\circ}\propto T(a)/a$, where $T(a)=\tanh[(a-a_{\mathrm{in}})/w_{\mathrm{in}}]\tanh[(a_{\mathrm{out}}-a)/w_{\mathrm{out}}]$ is a tapering function that is close to $1$ in most of the disc but declines smoothly to [math] at the inner and outer boundaries over length-scales $w_{\mathrm{in}}=0.01a_{\mathrm{in}}$ and $w_{\mathrm{out}}=0.01a_{\mathrm{out}}$. Apart from the taper, this model corresponds in the circular limit to a surface density $\propto r^{-1}$ and a sound speed $\propto r^{-1/2}$. The free boundary conditions correspond to regularity conditions at $a=a_{\mathrm{in}}$ and $a=a_{\mathrm{out}}$, which are singular points of equation (51). For this model they take the form

[TABLE]

where $\varepsilon^{\circ}_{\mathrm{in}}$ is the value of $\varepsilon^{\circ}$ at $a=a_{\mathrm{in}}$ neglecting the taper; this leads to a nonlinear relation between $e$ and $f$ at the boundaries.

The shapes of the first four modes are shown in Fig. 5, where the amplitudes have been chosen such a maximum eccentricity of $0.5$ is attained in each case. Again, the modes form a sequence in which the number of nodes in $e(a)$ increases by one. These modes are less confined than those of the previous problem; with free boundary conditions, $e$ does not need to vanish at the boundaries, and the broad extension of the disc means that the eccentricity gradient can avoid the highly nonlinear limit $f\to 1$ even though $e$ is as large as $0.5$. The precession rate again increases with amplitude and mode number and satisfies the relation $\mathrm{d}H/\mathrm{d}C=-\omega$ (Fig. 6). The numerical values of $-\omega$ in units of $\varepsilon_{\mathrm{in}}/l_{\mathrm{in}}$ are much smaller than the values of $-\omega$ in units of $c_{\mathrm{s}}/l_{\mathrm{in}}$ found in the previous problem (Fig. 4). This is because (i) the eccentricity gradients are more moderate, which reduces the precession rate; (ii) the internal energy in most of the disc is significantly less than $\varepsilon_{\mathrm{in}}$ because of the assumed temperature gradient; and (iii) the majority of the angular momentum is located at $a\gg a_{\mathrm{in}}$.

6 Conclusion

In this paper we have developed a new framework for studying the nonlinear dynamics of eccentric astrophysical discs. In the case of an ideal, non-self-gravitating fluid, a thin eccentric disc consists of a continuous set of nested elliptic orbits that evolve slowly as a result of pressure. The semimajor axis of each orbit is conserved, as is the total angular-momentum deficit of the disc, while the orbits undergo precession and exchanges of eccentricity. The nonlinear dynamics of this system can be derived from a convenient variational principle and expressed in a concise Hamiltonian form, using either canonical (e.g. modified Delaunay) variables or the more familiar Keplerian orbital elements. The Hamiltonian of a 2D disc is just the orbit-averaged internal energy of the disc, which can be determined from its eccentricity distribution using the geometry of the elliptical orbits. In the realistic case of a 3D disc, the Hamiltonian needs to be modified to take into account the dynamical vertical structure of the disc, but is still equal to a multiple of a suitably averaged internal energy.

The simplest solutions of the theory are uniformly precessing nonlinear eccentric modes, which make the Hamiltonian stationary subject to the angular momentum deficit being fixed. We have presented numerical examples of branches of modal solutions in illustrative situations up to their maximum possible amplitudes.

In many astronomical situations of interest, eccentric discs are expected to interact gravitationally with orbital companions; examples include circumbinary discs around young binary stars or binary black holes in galactic nuclei, circumstellar discs in cataclysmic variables and young binary stars, and protoplanetary discs with gap-opening planets. These interactions generally involve both secular effects, which contribute to precession and exchanges of eccentricity between the disc and the companion, and resonant effects, which contribute to the growth or damping of eccentricity. Both have been considered in a linear regime by Goodchild & Ogilvie (2006) and Teyssandier & Ogilvie (2016). It is likely that the Hamiltonian formalism developed here will be useful in understanding the nonlinear dynamics found in long-term simulations by Rosotti et al. (2017) and Ragusa et al. (2018), as it is a description of an ideal fluid that connects more readily with the established methods of celestial mechanics. It may also be helpful in understanding instabilities of eccentric discs that involve mode couplings (Ferreira & Ogilvie, 2008; Kato, 2008; Barker & Ogilvie, 2014). It would be valuable in future to explore further the nonlinear dynamics of 3D eccentric discs using this approach.

The theory developed here is intended to be easier to understand and to use than the nonlinear description of eccentric discs worked out by Ogilvie (2001). The older work considers the eccentricity as a function of semilatus rectum (or angular momentum) rather than semimajor axis (or energy), and requires a number of stress integrals to be computed. The new theory is all derived from a scalar Hamiltonian but does not include dissipative effects. It would be valuable to try to incorporate the effects of viscous dissipation, heating and cooling while preserving as much as possible of the unifying structure of this theory.

Acknowledgements

This research was supported by STFC through grants ST/L000636/1 and ST/P000673/1.

Appendix A Derivation of the Hamiltonian form of the equations

In the Hamiltonian theory of dynamics, any canonical transformation from phase-space variables $({\mn@boldsymbol{q}},{\mn@boldsymbol{p}})$ to $({\mn@boldsymbol{Q}},{\mn@boldsymbol{P}})$ satisfies the ‘direct conditions’

[TABLE]

We apply these relations with $({\mn@boldsymbol{q}},{\mn@boldsymbol{p}})$ being the Cartesian phase-space coordinates $(x,y,v_{x},v_{y})$ and $({\mn@boldsymbol{Q}},{\mn@boldsymbol{P}})$ being the action–angle variables $(\lambda,\gamma,\Lambda,\Gamma)$. We use the special symbol $\eth$ (rather than $\partial$) to denote partial derivatives of the four-dimensional phase-space transformation $(x,y,v_{x},v_{y})\mapsto(\lambda,\gamma,\Lambda,\Gamma)$ in order to distinguish them from those of the two-dimensional real-space transformation $(x,y)\mapsto(\Lambda,\lambda)$ associated with the orbital coordinate system. Thus, in particular,

[TABLE]

For a planar system, in the presence of a perturbing force per unit mass ${\mn@boldsymbol{f}}$, the equation of motion of a fluid element or test particle is

[TABLE]

where the dots represent the Newtonian force due to the central mass. Using the chain rule and the direct conditions, we deduce that the canonical variables that would be constant in a Keplerian orbit evolve according to

[TABLE]

(We will not require $\dot{\lambda}$, which in any case is dominated by the Keplerian mean motion $n$.)

The partial derivatives of the 2D transformation are related to those of the 4D transformation by

[TABLE]

The Jacobian matrix of the 2D transformation is

[TABLE]

with determinant

[TABLE]

and inverse

[TABLE]

For a perturbing force of the form

[TABLE]

where $P$ is a vertically integrated pressure (or other isotropic stress), we deduce that

[TABLE]

Using the chain rule on the partial derivatives of $P$ gives

[TABLE]

Multiplying by $J$ and using the inverse Jacobian matrix gives

[TABLE]

Now we notice from equation (62) that

[TABLE]

Furthermore, we have the differential identities

[TABLE]

which do not rely on any special properties of the functions $x(\lambda,\gamma,\Lambda,\Gamma)$ and $y(\lambda,\gamma,\Lambda,\Gamma)$. We now carry out an orbital time-average of the equations, denoted by angle brackets:

[TABLE]

This averaging operation, which also corresponds to a mass-weighted spatial average around the orbit, allows us to integrate by parts with respect to the periodic variable $\lambda$, with the results

[TABLE]

A.1 2D analysis

In the artificial but widely considered case of a 2D ideal fluid, variations of the vertically integrated pressure $P$ and the specific internal energy $\varepsilon$ around each orbit are adiabatic and related by $\mathrm{d}\varepsilon=-P\,\mathrm{d}A$, where $A=1/\Sigma=(2\pi/M_{\Lambda})J$ is the specific area. The evolutionary equations, multiplied by $2\pi$, then reduce to

[TABLE]

The first equation implies that $\Lambda$, and therefore $a$, is a material invariant in this theory, which is true because there is no dissipation of energy, nor does the pressure transfer any energy between neighbouring orbits. The three equations are in the canonical Hamiltonian form

[TABLE]

where the Hamiltonian $H$ is given by

[TABLE]

with Hamiltonian density

[TABLE]

involving the orbit-averaged specific internal energy $\langle\varepsilon\rangle$.

A.2 3D analysis

When the third ($z$) dimension perpendicular to the plane of the disc is considered, it can be seen that eccentric discs are never in hydrostatic equilibrium in this vertical direction. A non-zero eccentricity leads to an oscillatory variation of the vertical gravitational acceleration around the orbit, while a non-zero eccentricity gradient or twist leads to an oscillatory horizontal compression of the fluid. In either case the disc undergoes a forced vertical ‘breathing mode’ (Ogilvie, 2001; Ogilvie & Barker, 2014). The differential equation describing the dependence of the vertical scaleheight222There is an unfortunate notational clash between the scaleheight and the Hamiltonian. The context should make it clear which meaning $H$ has in any equation. $H$ around the orbit, for an ideal fluid disc that is stationary in an inertial frame on the orbital timescale, can be written as

[TABLE]

where

[TABLE]

describes the vertical gravity in a disc (assumed here to be thin), while $P$ and $\Sigma$ are the vertically integrated pressure and density. (Equation 84 is equivalent to equations 141–142 of Ogilvie & Barker (2014).) The mass-weighted vertically averaged specific internal energy of a perfect gas is

[TABLE]

Adiabatic flow means that $(P/H)\propto(\Sigma/H)^{\gamma}$ around the orbit, so $\mathrm{d}\bar{\varepsilon}=-(P/H)\,\mathrm{d}(AH)$ when $\Lambda$ is held constant.

The associated energy equation is obtained by multiplying equation (84) by $\partial H/\partial\lambda$:

[TABLE]

The three terms in the square brackets are the kinetic, gravitational and internal energies per unit mass, after vertical averaging. The two source terms on the right-hand side are due to variable vertical gravity and horizontal orbital compression, respectively. An orbital average gives

[TABLE]

On the other hand, multiplying equation (84) by $H$ and averaging (with integration by parts) gives the virial relation

[TABLE]

Now consider how the scaleheight responds to variations in the geometry of the disc. Note that, of the two agents driving the vertical oscillation, $\Psi$ depends on $(\Gamma,\gamma)$ as well as $(\Lambda,\lambda)$, while $A$ also depends on $\Gamma_{\Lambda}$ and $\gamma_{\Lambda}$. Consider, for example, variations of $\Gamma$. We differentiate equation (84) with respect to $\Gamma$ to obtain

[TABLE]

We then multiply by $H$ and average (with integration by parts):

[TABLE]

Using the property

[TABLE]

we can rearrange to obtain

[TABLE]

where

[TABLE]

Using the virial relation, this quantity can be rewritten either as

[TABLE]

which is a multiple of the average internal energy, or as

[TABLE]

which is the average of (potential–kinetic) energy. Indeed, this is (minus) the averaged Lagrangian that appears in the theory of nonlinear dispersive waves (Whitham, 1965).

It can be shown similarly that

[TABLE]

When the additional horizontal gravitational force (with mass-weighted vertical averaging) due to $\Psi$ is added to the evolutionary equations, we have

[TABLE]

Using equations (88), (93) and (97), we see that $\langle\dot{\Lambda}\rangle$ evaluates to zero, and that the equations are of Hamiltonian form, with Hamiltonian density given by

[TABLE]

Appendix B Dimensionless forms of the Hamiltonian density

B.1 2D analysis

For a perfect gas of adiabatic index $\gamma$, adiabatic flow around each orbit means that $P=K\Sigma^{\gamma}$, where $K(a)$ is related to the entropy distribution, which does not evolve. The specific internal energy is then

[TABLE]

Recall that equation (15) relates the surface density $\Sigma$ to the mass distribution $M_{\Lambda}$ and the Jacobian $J=J^{\circ}j$. If the disc had the same distributions of mass and entropy but were circular, it would have $\Sigma=\Sigma^{\circ}=M_{\Lambda}/2\pi J^{\circ}$ and $P=P^{\circ}=K(\Sigma^{\circ})^{\gamma}$. Thus $\Sigma=\Sigma^{\circ}j^{-1}$ and

[TABLE]

where

[TABLE]

is a given function of $a$ and

[TABLE]

is the dimensionless, geometric part of the Hamiltonian density, which depends on $e$, $ae_{a}$, $ae\varpi_{a}$ and $\gamma$.

In the case of a (globally) isothermal disc we instead have

[TABLE]

where $c_{\mathrm{s}}=\mathrm{constant}$ is the isothermal sound speed. Then

[TABLE]

where

[TABLE]

and

[TABLE]

This expression agrees with the limit $\gamma\to 1$ of equation (105) after the constant term $1/(\gamma-1)$ (which, like the terms depending only on $a$ in equation 107, does not affect the dynamics) is removed from it. This follows from the result

[TABLE]

B.2 3D analysis

For a 3D disc, adiabatic flow means instead that $(P/H)=K(\Sigma/H)^{\gamma}$, where $K$ is a fixed function of $a$. If we write the scaleheight as $H=H^{\circ}(a)h(E)$, where $H^{\circ}=(P^{\circ}/\Sigma^{\circ}n^{2})^{1/2}$ is the (hydrostatic) vertical scaleheight of a circular disc with the same mass and entropy distributions, then we have $\Sigma=\Sigma^{\circ}j^{-1}$ and $P=P^{\circ}j^{-\gamma}h^{-(\gamma-1)}$. Equation (84) reduces to a dimensionless ODE for the dimensionless scaleheight $h$:

[TABLE]

or, in terms of derivatives with respect to the eccentric anomaly,

[TABLE]

We then have the Hamiltonian density

[TABLE]

with

[TABLE]

and

[TABLE]

In the isothermal case we have instead

[TABLE]

[TABLE]

[TABLE]

Appendix C Analaytical expressions for the Hamiltonian of a 2D disc

The quantity $j$ in square brackets in equation (105) or (109) can be written as

[TABLE]

where $\theta=E-\alpha$, with the amplitude $q$ and phase $\alpha$ satisfying

[TABLE]

which implies

[TABLE]

To avoid orbital intersection we require $|q|<1$. In the untwisted case $\varpi_{a}=0$ it is natural to take $\alpha=0$ and regard $q$ as a signed quantity. $q$ is a measure of the degree of areal compression and the level of nonlinearity in an eccentric disc; it generalizes the nonlinearity parameter used in the theory of planetary rings in the regime $e\ll 1$.

For general $\gamma$ the integral (105) can be evaluated in terms of Legendre functions. The special case $\gamma=2$ and the limiting isothermal case $\gamma=1$ involve only elementary functions.

Let

[TABLE]

for $|q|<1$. This is an even function of $q$ with $I_{p}(0)=1$. For $p>0$ or $p<-1$ it is a monotonically increasing function of $q$ for $0<q<1$. For $p>1/2$ it diverges proportionally to $(1-q)^{(1/2)-p}$ as $q\to 1$. For $-1<p<0$ it is monotonically decreasing. For $p=0$ or $p=-1$ it is constant and equal to $1$. Special cases include

[TABLE]

In general, $I_{p}(q)$ can be written in terms of the Legendre function:

[TABLE]

A series expansion, convergent for all $|q|<1$, is

[TABLE]

where

[TABLE]

which behaves as

[TABLE]

Thus

[TABLE]

To express the integral (105) in terms of $I_{p}(q)$, we use equation (119) and change variables from $E$ to $\theta$. Note that

[TABLE]

Thus, for $\gamma>1$,

[TABLE]

with

[TABLE]

When $\gamma=2$ we have a result in terms of elementary functions:

[TABLE]

In the untwisted case this reduces to

[TABLE]

where $f=e+ae_{a}$.

For the isothermal case $\gamma=1$ we use instead

[TABLE]

leading to

[TABLE]

In the untwisted case this reduces to

[TABLE]

The apparent singularities at $e=f$ in the untwisted expressions can be removed by using the trigonometric parametrization $e=\sin 2\alpha$, $f=\sin 2\beta$, $\sqrt{1-e^{2}}=\cos 2\alpha$, $\sqrt{1-f^{2}}=\cos 2\beta$. This leads to

[TABLE]

for $\gamma=2$ and

[TABLE]

for $\gamma=1$.