{'fraction_non_alphanumeric': 0.035413190968746525, 'fraction_numerical': 0.012367923479034591, 'mean_word_length': 4.512015693967631, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Local and global dynamics of eccentric astrophysical discs

We formulate a local dynamical model of an eccentric disc in which the dominant motion consists of elliptical Keplerian orbits. The model is a generalization of the well known shearing sheet, and is suitable for both analytical and computational studies of the local dynamics of eccentric discs. It is spatially homogeneous in the horizontal dimensions but has a time-dependent geometry that oscillates at the orbital frequency. We show how certain averages of the stress tensor in the local model determine the large-scale evolution of the shape and mass distribution of the disc. The simplest solutions of the local model are laminar flows consisting of a (generally nonlinear) vertical oscillation of the disc. Eccentric discs lack vertical hydrostatic equilibrium because of the variation of the vertical gravitational acceleration around the eccentric orbit, and in some cases because of the divergence of the orbital velocity field associated with an eccentricity gradient. We discuss the properties of the laminar solutions, showing that they can exhibit extreme compressional behaviour for eccentricities greater than about $0.5$, especially in discs that behave isothermally. We also derive the linear evolutionary equations for an eccentric disc that follow from the laminar flows in the absence of a shear viscosity. In a companion paper we show that these solutions are linearly unstable and we determine the associated growth rates and unstable modes.


Astrophysical motivation
Eccentric discs, in which the dominant motion consists of elliptical Keplerian orbits, occur in a wide variety of astrophysical situations. For example, an eccentric gaseous disc is formed directly when a star (or a giant planet) evolves, through scattering or secular interaction, on to an orbit that closely approaches the galactic centre (or the host star), and is tidally disrupted (e.g. Gurzadian & Ozernoi 1979;Guillochon, Ramirez-Ruiz & Lin 2011); such a process might be responsible for the gas cloud G2 near Sgr A * in the Galactic Centre (Guillochon et al. 2014).
In an eccentric binary star, a circumstellar or circumbinary disc acquires a forced eccentricity from the binary orbit via secular gravitational interaction. The importance of this for planet formation in binary stars has been recognized (Paardekooper, Thébault & Mellema 2008). Even if the binary orbit is circular, certain mean-motion resonances can allow a free eccentricity of the disc to grow, initially exponentially (Lubow 1991a). These effects can of course occur in non-stellar binaries such as binary black holes with accretion discs (Armitage & Natarajan 2005), planet-satellite systems with planetary rings (Goldreich & Tremaine 1981;Borderies, Goldreich & Tremaine 1983) and protoplanetary systems (Kley & Dirksen 2006). Whether planet-disc interactions lead to eccentricity excitation (Goldreich & Sari 2003;Ogilvie & Lubow 2003;D'Angelo, Lubow & Bate 2006) depends on the dynamics of eccentric discs, because of the strong coupling between the planet and the disc.
Even in the absence of an orbiting companion, a disc may become eccentric through an instability of the circular state, such as viscous overstability (Kato 1978;Ogilvie 2001). In contrast to the naive expectation that viscosity tends to circularize a disc, viscous overstability may explain the eccentricity of decretion discs formed around rapidly rotating Be stars (e.g. Rivinius, Carciofi & Martayan 2013).
In eccentric binaries the forced eccentricity of the disc is locked to that of the binary and may not be easily detectable. However, discs with a free eccentricity precess as a result of their pressure and any gravitational influences that cause a departure from Keplerian motion. This is the generally accepted explanation of the superhump phenomenon in the SU UMa class of dwarf novae (e.g. Warner 1995), in which the accretion disc expands sufficiently during superoutbursts to encounter the 3:1 resonance with the binary orbit and becomes eccentric. (Some other systems exhibit steady accretion and permanent superhumps.) The elliptical outer rim of the disc in OY Car was measured by Hessman et al. (1992) through the variation of eclipses of the hot spot.
[For a critical analysis of observational evidence for eccentric discs in SU UMa stars from a particular standpoint, see Smak (2009).] Recently, Kepler has been used to observe superoutbursts and superhumps with much greater accuracy . In dwarf novae the optical emission is modulated at the frequency at which the disc precesses in the frame that rotates with the binary orbit, owing to the interaction between the eccentric mode and the tidal deformation. Related phenomena are also reported in low-mass X-ray binaries, although the radiative mechanisms are different (Haswell et al. 2001).

Theoretical and computational background
Several theoretical and computational approaches have been taken to the study of eccentric discs. One is to try to generalize the classical theory of viscous accretion discs to allow for orbits of arbitrary eccentricity (Syer & Clarke 1992;Lyubarskij, Postnov & Prokhorov 1994;Ogilvie 2001). These analyses, of which the last is by far the most general, aim to derive evolutionary equations for the shape and mass distribution of eccentric discs due to viscous and other internal stresses. Earlier, equations governing the evolution of narrow and slightly eccentric planetary rings were formulated by Borderies, Goldreich & Tremaine (1983). A separate body of theoretical work relates to eccentric collisionless stellar discs in galactic nuclei, notably M31 (Tremaine 1995;Peiris & Tremaine 2003).
The broad conclusion of this work is that eccentricity can propagate through a disc by means of pressure and selfgravity, as a slow one-armed density wave, while viscosity causes it to diffuse (except in cases where it is excited by viscous overstability). Differential apsidal precession due to the rapid rotation of the central object, relativistic effects, self-gravity of the disc or the presence of orbiting companions can be important, as can three-dimensional effects due to the vertical structure and oscillation of the disc.
Numerical simulations of eccentric discs have mainly been carried out using smoothed particle hydrodynamics (SPH), which readily produces eccentric discs in circular binary stars with mass ratios typical of SU UMa stars (Whitehurst 1988;Lubow 1991b;Murray 1996Murray , 1998Murray , 2000Smith et al. 2007). More recently, grid-based simulations have also found the development of eccentric discs in the presence of a planetary (Kley & Dirksen 2006) or stellar (Kley, Papaloizou, & Ogilvie 2008;Marzari et al. 2009Marzari et al. , 2012 companion. Papaloizou (2005a) found that eccentric discs are hydrodynamically unstable in the absence of viscosity. The instability is three-dimensional and takes the form of a parametric resonance of inertial waves, as also occurs in tidally distorted discs (Goodman 1993) and in the classic elliptical instability of flows with non-circular streamlines (Kerswell 2002). A related phenomenon occurs in warped discs (Ogilvie & Latter 2013b). Papaloizou (2005b) carried out numerical simulations of the instability of eccentric discs in the absence of vertical gravity and found that it led to subsonic turbulence.

Plan of this paper
This paper is organized as follows. In Section 2 we describe the geometry of an eccentric disc and recall the properties of the orbital coordinates defined by Ogilvie (2001). We formulate the hydrodynamic equations in this coordinate system and obtain the evolutionary equations for an eccentric disc in terms of orbital averages of force and stress components. In Section 3 we derive a local model of an eccentric disc, which will be useful for analytical and computational studies of instabilities and turbulence in eccentric discs. In Section 4 we consider the simplest hydrodynamic solutions of this local model, which are non-hydrostatic and necessarily involve a vertical oscillation of the disc; we also discuss the evolution of eccentric discs under this laminar dynamics.
In a companion paper (Barker & Ogilvie 2014) we use the local model to analyse the linear hydrodynamic stability of an eccentric disc.

Introduction
In a thin astrophysical disc, the orbital motion is hypersonic and fluid elements follow ballistic trajectories to a first approximation. Around a spherical central mass, these trajectories are Keplerian orbits, which can have eccentricity and inclination. A general Keplerian disc involves smoothly nested orbits of variable eccentricity and inclination: it is both elliptical and warped. The case of a warped disc composed of variably inclined but circular orbits is more familiar and was treated recently by Ogilvie & Latter (2013a). In this paper we consider instead the case of an eccentric disc composed of variably elliptical but coplanar orbits. The general case of a warped and eccentric disc remains for future work.
The dominant motion in an eccentric disc is orbital motion in the form of Keplerian ellipses. Orbital precession due to a small departure of the gravitational potential from that of a point mass, or due to weak relativistic effects, can be treated, along with the collective effects of the disc, as a perturbation of the Keplerian motion. The eccentric disc can therefore be considered, to a first approximation, as a continuum of nested elliptical rings (Fig. 1), whose shape can be regarded as fixed in a non-rotating frame on the timescale of the orbital motion.
Following Ogilvie (2001), we label the orbits using their semi-latus rectum λ = a(1 − e 2 ), where a and e are the semi-major axis and the eccentricity. The semi-latus rectum is directly related to the specific angular momentum = (GM λ) 1/2 , where G is Newton's constant and M is the central mass. In an eccentric disc, the eccentricity e and the longitude of pericentre ω can be regarded as functions of λ, and can be conveniently combined into the complex eccentricity E = e e iω . (Note that the eccentricity, e, and the base of natural logarithms, e, are distinguished typographi- Figure 1. Four examples of eccentric discs. Top left: uniformly eccentric disc with e = 0.5 and ω = 0. Top right: untwisted disc with ω = 0 and e increasing linearly with λ from 0 to 0.5. Bottom left: untwisted disc with ω = 0 and e decreasing linearly with λ from 0.5 to 0. Bottom right: twisted disc with e = 0.5 and ω increasing logarithmically with λ such that λω = 1. cally.) The real and imaginary parts of E are equivalent to the components of the eccentricity vector (e cos ω, e sin ω).

Orbital coordinates
Again following Ogilvie (2001), we make use of orbital coordinates (λ, φ) in the plane of the disc (Fig. 2), instead of polar coordinates (r, φ). The transformation between them is given by the polar equation of an ellipse, . (1) The semi-latus rectum λ can be thought of as a quasi-radial coordinate that replaces r, while φ is the usual azimuthal angular coordinate.
Orbital coordinates are, of course, well adapted to the geometry of an eccentric disc, in which the Keplerian ellipses correspond to the curves λ = constant. The disadvantage of these coordinates is that they are not orthogonal. Since we will need to carry out vector and tensor calculus in these coordinates, we must first define a number of geometrical quantities, quoting results from Ogilvie (2001).
The components of the metric tensor gij and of its inverse g ij are where the subscripts on R denote partial derivatives of the function R(λ, φ). The components of the Levi-Civita connection (Christoffel symbols) are The orbital coordinate system is easily extended to three dimensions by adding the third coordinate z. Apart from gzz = g zz = 1, all metric and connection components involving z vanish. The Jacobian of the coordinate system is In order that J > 0, i.e. that the orbits are closed and nested without intersection, we require |E| < 1 and |E − λE | < 1, where the prime denotes a derivative with respect to λ. Therefore both the eccentricity and the eccentricity gradient must be sufficiently small. Further useful geometrical relations are derived in Appendix A.
Although by convention e ≥ 0, the same orbit is obtained by reversing the sign of e and increasing (or decreasing) ω by π, which leaves the complex eccentricity E unchanged. In situations where E(λ) has a simple zero corresponding to a circular orbit, E is well defined and non-zero on that orbit. However, as conventionally defined, e has a discontinuous gradient there and ω changes abruptly by π. In order to avoid this discontinuity and to construct a local model in such a case, we will allow e to become negative on one side of the circular orbit so that e is continuous there.
In principle the orbital coordinates are time-dependent as well as non-orthogonal. The difficulties that this introduces were treated by Ogilvie (2001). In this paper we mainly circumvent these difficulties by assuming that the orbital geometry is fixed on the orbital timescale. However, the slow time-dependence of the coordinates is taken into account in deriving the evolutionary equations for the eccentric disc in Section 2.5.

Hydrodynamic equations
We consider an ideal fluid satisfying the equation of motion, the equation of mass conservation, and the thermal energy equation, where u i is the velocity, ∇i is the covariant derivative, Φ is the gravitational potential, ρ is the density, p is the pressure, s is the specific entropy and D = ∂t + u i ∂i is the Lagrangian derivative acting on scalar fields. We assume that the disc is of sufficiently low mass that self-gravity may be neglected. We work with the contravariant velocity components u i , which are simply the rates of change of the orbital coordinates of fluid elements. Equation (8) can be written using partial derivatives as while equation (9) can be written either as or in the conservative form An alternative to equation (10) is where γ = (∂ ln p/∂ ln ρ)s is the adiabatic index. The conservative form of the equation for the total energy of the fluid is where is the specific total energy and ε is the specific internal energy. This result follows from equations (11) and (13) when the relations dε = T ds−p d(ρ −1 ) and ∂ k gij = g il Γ l jk +g lj Γ l ik are used, provided that the gravitational potential is independent of t. In detail, equations (11), (12) and (14) are and with D = ∂t + u λ ∂ λ + u φ ∂ φ + u z ∂z. The alternative, conservative form of equation (20) is (In Section 2.5 we will use a modified form of this equation that takes into account the slow time-dependence of the orbital coordinates.)

Orbital motion
The orbital motion corresponds to the velocity field u i = ω i , where ω φ = Ω(λ, φ) is the orbital angular velocity and ω λ = ω z = 0. This velocity field should satisfy the equation of motion in the midplane z = 0, where Φ = Φ0(R), when the pressure is neglected, i.e.
Since Φ0 is a function of R only, these equations simplify to and are satisfied, as expected, when Φ0 = −GM/R and

Evolution of mass, angular momentum and eccentricity
A principal aim of a theory of eccentric discs is to obtain a system of equations that govern the evolution of the shape and mass distribution of the disc. Unlike the case of a warped disc composed of circular orbits, these equations do not follow simply from the conservation of mass and angular momentum. The evolution of the complex eccentricity, or eccentricity vector, is more subtle and is not purely conservative in nature.
We consider first the case of a test particle in a Keplerian orbit in the plane z = 0 and subject to a perturbing force within that plane. Its motion is governed bÿ where = r 2φ is the specific angular momentum and fr and f φ are the (orthogonal) polar components of the perturbing force per unit mass. The osculating orbital elements are defined by equating the instantaneous position and velocity of the particle with those of a Keplerian orbit. Thus where λ = 2 /GM is the semi-latus rectum, e is the eccentricity, θ = φ − ω is the true anomaly and ω is the longitude of pericentre. From the above relations we have and so which therefore evolves according to Although unfamiliar in this form, this equation is equivalent to the Gauss perturbation equations of celestial mechanics in the case of planar motion. If instead we use the contravariant orbital components f λ and f φ , which are related by fr = We now consider a continuous disc. The equation of mass conservation in a three-dimensional conservative form that takes into account the time-dependence of the orbital coordinates is (Ogilvie 2001) (36) whereλ is the rate of change of λ with time in an inertial coordinate system, due to the slow evolution of the orbital geometry. Integrating this equation with respect to φ and z over the full extent of the disc, and assuming that no mass is gained or lost vertically, we obtain the one-dimensional conservative form where is the one-dimensional mass density with respect to λ, Σ = ρ dz being the surface density, and is the quasi-radial mass flux. Note that the mass of the disc is where the integral is carried out over an appropriate range of λ. We can also write wherev λ = F/M is the mean quasi-radial velocity.
Given that is independent of φ and z, the angular momentum equation has the one-dimensional form or, equivalently, which is the continuum analogue of equation (34). Let us consider the case of internal perturbing forces that are due to stress divergences, i.e. ρf i = ∇jT ij , where T ij is a symmetric stress tensor describing the collective effects of the disc (pressure, viscosity, self-gravity, etc.). Then (Ogilvie 2001) The terms in JρR 2 f φ involving ∂ φ and ∂z integrate to zero, assuming suitable boundary conditions in z that ensure that no angular momentum is lost or gained vertically, and it follows that which can be written in the conservative form where is the internal torque. With the help of the equation of mass conservation (37), and the fact that depends only on λ, it simplifies to which determines the mass flux F (or the mean quasi-radial velocityv λ = F/M) instantaneously in terms of the torque distribution. As expected, the stress component T λφ determines the redistribution of angular momentum within the disc and thereby regulates the accretion flow. The eccentricity equation is less obvious because it is not conservative. We can expect the continuum analogue of equation (35) to be When ρf i = ∇jT ij , the terms involving ∂z again integrate to zero. The terms involving ∂ φ are less obvious and an integration by parts is needed to obtain Using the relation (A9) we obtain This result is consistent with equation (167) of Ogilvie (2001), which was derived more formally using the method of multiple time-scales, and can be written in various ways (especially concerning where to place the ∂ λ ).
In order to close the equations governing the shape and mass distribution of the disc, the four stress integrals that are needed are therefore External forces acting on the disc, which would include the force due to any departure of the gravitational potential from that of a point mass, also contribute to the evolution of and E in the obvious way. Thus the governing equations in the presence of internal forces (described by T ij ) and external forces (described by f i ) are An important aspect of the three-dimensional theory of eccentric discs is that the point-mass potential makes a nonzero contribution to f λ . When the potential Φ = −GM (R 2 + z 2 ) −1/2 is expanded in a Taylor series about the midplane z = 0 of a thin disc, we obtain where Φ0 = −GM/R is the potential in the midplane, while Φ2 = GM/R 3 has a different dependence on R. Its contribution to f λ is Its contribution to M(∂t +v λ ∂ λ )E is therefore which is the first term on the right-hand side of equation (167) of Ogilvie (2001). This is a three-dimensional effect due to the weakening of the (cylindrical) radial gravitational force away from the midplane. The cylindrical radial component is the relevant one because the gas away from the midplane is moving in a plane of (approximately) constant z, rather than in an inclined Keplerian orbit. Although the integrals of stresses and forces that appear in the evolutionary equations involve integrals with respect to the azimuthal angle and are in this sense global or largescale quantities, we will see below that these integrals naturally emerge in the form of time-averages in a local model that follows the orbital motion.

Flow decomposition
We have seen in Section 2.4 that the eccentric orbital motion with u φ = Ω(λ, φ) and u λ = u z = 0 satisfies the equation of motion in the midplane z = 0 when the pressure is neglected. We now write the fluid motion as the sum of this orbital motion and a relative velocity v i : The residual parts of the hydrodynamic equations (17)-(21) are then, without approximation, where is the Lagrangian derivative and is the orbital velocity divergence, which vanishes only in the case E = constant.
To simplify the equations in a way appropriate for a local model of a thin disc, we apply the following scaling argument. Let 1 be a characteristic value of the aspect ratio H/r of the disc, and let us consider a system of units in which r and Ω are O( 0 ), so that H and the sound speed are O( 1 ). Therefore p/ρ and Φ − Φ0 are O( 2 ). We are interested in describing local nonlinear fluid dynamical phenomena that take place on a lengthscale comparable to H and on a timescale comparable to Ω −1 . Therefore, when acting on v i , ρ or p, the operator ∂i is O( −1 ), while the operator D is O( 0 ). We assume that the relative velocity components v i are O( 1 ), i.e. comparable in magnitude to the sound speed but much smaller than the orbital velocity. This allows them to be nonlinear and to generate Reynolds stresses comparable to the pressure. The equations simplify at the leading order in to where Φ2 is defined in equation (60). Although several terms involving products of components of v i have been dropped in equations (71) and (72), these equations are still nonlinear because v i appears in the Lagrangian derivative. The 'inertial' terms involving the interaction between the orbital motion and the relative velocity are linear, however. The Jacobian J has dropped out of equations (74) and (75) because it varies on a lengthscale O( 0 ). For similar reasons there are no horizontal derivatives of Φ2 in these equations.
We examine the neighbourhood of the reference point by letting where (ξ, η) are (non-orthogonal) local coordinates in the orbital plane, while ζ and τ are introduced only for notational uniformity. As we are interested in a region with an extent comparable to H in each dimension around the reference point, ξ, η and ζ are small, O( ).
Since we are evaluating equations (71)-(74) at the leading order in , the geometrical coefficients appearing in those equations may be replaced with their values at the reference point (λ, φ) = (λ0, ϕ(τ )), making them functions of τ only, with period P0.
Because of its appearance in the Lagrangian derivative, The local orbital coordinates (ξ, η, ζ) are centred on a reference point that follows a particular Keplerian orbit (red ellipse). The blue boxes represent the azimuthal compression due to the nonuniformity of the orbital motion, discussed in Section 3.3, although they do not represent the orbital shear.
however, the orbital angular velocity Ω(λ, φ) needs to be expanded about the reference point in the form where Ω0(τ ) = Ω(λ0, ϕ(τ )) is the angular velocity of the reference point, while Ω λ0 (τ ) and Ω φ0 (τ ) are the partial derivatives ∂ λ Ω and ∂ φ Ω of Ω(λ, φ) evaluated at that point. Since, at a fixed time, the differentials of the local coordinates are identical to those of the global orbital coordinates, the contravariant velocity components and spatial derivatives are unchanged by the transformation. Thus (v ξ , v η , v ζ ) = (v λ , v φ , v z ) and (∂ ξ , ∂η, ∂ ζ ) = (∂ λ , ∂ φ , ∂z). However, the time-dependence of the transformation means that the time-derivatives are related by At the leading order in , therefore, Equations (71)-(74) are then rewritten as in which, as explained above, the geometrical coefficients are evaluated at the reference point (λ, φ, z) = (λ0, ϕ(τ ), 0), which makes them periodic functions of τ . We have also dropped the subscript zeros, so that The thermal energy equation can be written as either These are the equations of ideal hydrodynamics in the local model of an eccentric disc in non-shearing coordinates. In Appendix B we give expressions for the coefficients appearing in these equations. Since these involve the true anomaly θ explicitly, and it is not straightforward to express θ in terms of the time τ , it may be preferable to regard θ as the timelike variable instead of τ , replacing the derivative ∂τ with Ω∂ θ . The metric coefficients are now functions of τ satisfyinġ gij = Ω(g ik Γ k jφ + g kj Γ k iφ ), where the dot denotes a timederivative. Similarly the Jacobian and the orbital velocity divergence become functions of τ related by ∆ = (J/J)+Ω φ , with Ω φ =Ω/Ω. Thus the conservative form of equation (83) is (87) The Lagrangian derivative of the specific kinetic energy of the relative motion can be shown to be The first source term on the right-hand side involves the covariant derivative ∇iω j = ∂iω j + Γ j ik ω k of the orbital velocity field. Here vi = gijv j are the covariant components of the relative velocity. The conservative form of the energy equation for the relative motion is where Note that the possible sources of energy for the local model are the orbital shear (accessed through Reynolds stresses), the time-dependence of the vertical gravity coefficient Φ2, and the orbital velocity divergence.

Shearing coordinates
The equations of the local model derived so far have an explicit dependence on the horizontal coordinates ξ and η because of their appearance in the Lagrangian derivative (but only for 'non-axisymmetric' solutions that depend on η). We can derive a horizontally homogeneous model by transforming to shearing (and oscillating) local coordinates (ξ , η , ζ , τ ) that are Lagrangian with respect to the local orbital motion and are defined by where α and β remain to be determined. The scale factor α corresponds to a time-dependent stretching of the azimuthal coordinate (cf. Fig. 3), while the term β corresponds to a shearing of the coordinate system. Partial derivatives transform according to The Jacobian of the shearing coordinates is J = J/α, the factor of α coming from the transformation just introduced. We will not transform the velocity components because we are interested in evaluating Reynolds stress coefficients such as T λφ = T ξη . The Lagrangian derivative transforms into the spatially homogeneous form provided that α and β are chosen such thaṫ Since Ω φ =Ω/Ω, the first equation is satisfied by choosing α ∝ Ω −1 . This periodically variable scale factor for the azimuthal direction means that the new coordinate η is related to the mean anomaly or mean longitude, i.e. to the time taken to travel along the orbit. We choose α to have dimensions of length, so that the new variable η does also. The factor β has a secular dependence on time, representing orbital shear, as well as a periodic one. The expressions for α and β are The relation between the time τ (with origin τ = 0 at pericentre) and θ is given by Kepler's equation in the form where n = (GM/a 3 ) 1/2 is the mean motion, a = λ/(1 − e 2 ) is the semi-major axis and E is the eccentric anomaly, such that Thus τ can be expressed in terms of θ by inverting this sine function. The explicit expression for the Jacobian of the shearing coordinates is The equations of the local model now have the form Dv ξ + 2Γ λ λφ Ωv ξ + 2Γ λ φφ Ωv η = − 1 ρ g λλ (∂ ξ + β∂ η )p + g λφ α∂ η p , The thermal energy equation can be written as either In Appendix C we give, for future reference, the extension of these equations to ideal magnetohydrodynamics. Significantly, the coefficients in these equations are independent of ξ and η : the model is homogeneous in the local horizontal spatial coordinates. However, the coefficients do depend on τ . This time-dependence is mostly periodic, with the orbital period, but terms associated with azimuthal derivatives (∂ η ) have in addition a secular dependence on τ through β.
This spatial homogeneity means that periodic boundary conditions can be applied in ξ and η , leading to the model of an eccentric shearing sheet or box. In this case ξ and η are restricted to the ranges [−L ξ /2, L ξ /2] and [−L η /2, L η /2] respectively, identifying a finite patch of fluid. However, because of the secular dependence of β on τ , the shearing coordinates should be remapped from time to time to prevent them from becoming too distorted. This can be done, for example, by reducing the value of the additive constant in equation (97) from time to time so that β remains within a reasonable range of positive and negative values. When β is changed from β1 to β2, the relationship between the non-shearing coordinates (ξ, η) and the shearing coordinates (ξ , η ) is modified; in order for the periodic boundary conditions to be mutually compatible in the old and new coordinates, it is necessary that β1 − β2 be an integer multiple of the aspect ratio L η /L ξ .
The local model inherits three dimensionless parameters from the geometry of the eccentric disc: e (eccentricity), λe (eccentricity gradient) and λeω (twist). This makes it more complicated than the local model of warped discs (Ogilvie & Latter 2013a).
Using the fact that the orbital velocity divergence is ∆ = ∂τ ln(JΩ) = ∂τ ln(J/α) = ∂τ ln J , we can express the conservation of mass in the form (107) Subject to periodic boundary conditions in ξ and η , and suitable boundary conditions in ζ , the conserved mass in the shearing box is The horizontal momentum components P ξ = v ξ dM and P η = v η dM of the box satisfy the equations which allow an epicyclic oscillation around the reference orbit, but it would be natural to constrain the solutions to satisfy P ξ = P η = 0, meaning that the reference orbit has been correctly defined. The energy equation in shearing coordinates has the form In Section 2.5 we saw that four different integrals of components of the stress tensor, (53)-(56), are required in order to close the system of equations governing the global evolution of the shape and mass distribution of an eccentric disc. In the local model the stress components can readily be calculated; for example, the Reynolds-stress component T λφ corresponds to −ρv ξ v η . The geometrical factors J, R, e iφ and R λ are known functions of time in the local model. Since the local model follows an orbiting reference point, the azimuthal integral . . . dφ can be interpreted as a time-integral . . . Ω dτ over a single orbit, where Ω is also a known function of time. In the case of turbulent flows, additional spatial averaging over the box and time-averaging over multiple orbits can be carried out to obtain the relevant stress integrals.

Relation to the standard shearing sheet
If e = 0, the reference orbit is circular and the model can be related to the standard shearing sheet. We then have λ = R and Ω = (GM/R 3 ) 1/2 = constant on the reference orbit, and the metric components and connection coefficients simplify considerably.
The equations of the local model (in non-shearing coordinates) reduce to with and θ = Ωt. These equations can be derived from the standard hydrodynamic equations of the shearing sheet, with (Note that, since these are Cartesian coordinates, v x is the same as vx, etc.) If e = 0 then the relation is straightforward because J = 1 and ∆ = 0. Note that Rv η equates to v y , and (1/R)∂η to ∂y, because η is an angular variable.
If e = 0 then a time-dependent homogeneous transformation of the radial coordinate is involved (cf. Latter & Ogilvie 2009, Appendix A): so that and the velocity components are related by The solution with v ξ = v η = 0 corresponds to a free epicyclic oscillation proportional to ξ, which is the local representation of the eccentricity gradient. Models of this type have been used in the theory of planetary rings (Mosqueira 1996), in which the parameter |λe | < 1 is called q.

Nonlinear vertical oscillations
Laminar flows are the simplest solutions of the local model, being horizontally invariant and having a purely vertical velocity v ζ that, like ρ and p, depends only on ζ and τ . They satisfy the equations The periodic variation of the vertical gravity coefficient Φ2 around the orbit (in the case E = 0) and the orbital velocity divergence ∆ (in the case E = 0) drive a vertical oscillation of the disc, which is nonlinear unless E and λE are both small.
Before solving these equations we consider the simpler problem of hydrostatic equilibrium in a circular disc where the vertical gravitational acceleration is proportional to the distance above the midplane. Let Fρ(x) and Fp(x) be dimensionless functions of a dimensionless vertical coordinate x, which describe the equilibrium profiles of density and pressure. The equation of hydrostatic equilibrium in dimensionless form is We normalize the profiles such that (The latter two integrals are easily shown to be equal by integrating by parts and applying reasonable boundary conditions.) Simple examples are the isothermal structure, the homogeneous structure, (for x 2 < 3 only), and the polytropic structure, Fp(x) = 2n + 3 2(n + 1) (for x 2 < 2n + 3 only) where n > 0 (not necessarily an integer) is the polytropic index and Cn = [(2n + 3)π] −1/2 Γ(n + 3 2 ) Γ(n + 1) is a normalization constant. It can be shown that the polytropic structure approaches the isothermal structure in the limit n → ∞, and approaches the homogeneous structure in the limit n → 0.
Provided that γ = constant, the laminar flow has the form of a homogeneous expansion and contraction of the disc, where x = ζ/H(τ ) is the vertical coordinate ζ scaled by a time-dependent vertical scaleheight H(τ ), and the dimensionless functions Fρ(x) and Fp(x) are as defined above. Although the disc is not in hydrostatic equilibrium, its internal structure can be related to that of an equilibrium disc. In the isothermal case H(τ ) is the Gaussian scaleheight of the disc, while in the homogeneous case or the polytropic case it is a fraction of the true semi-thickness. In each case the surface density is Σ(τ ) = ρ dζ =ρ(τ )H(τ ) and the second vertical moment of the density is ρζ 2 dζ =ρ(τ )H(τ ) 3 , so H(τ ) is the standard deviation of the mass distribution. Equations (127)-(129) are satisfied provided that the functions H(τ ), w(τ ),ρ(τ ) andp(τ ) obeẏ Note that, if γ > 1 + n −1 (or γ > 1 in the case of an isothermal structure), the disc is stably stratified. However, buoyancy forces do not affect the dynamics of the laminar flow because of its horizontal invariance.
The surface density satisfieṡ This is a statement of mass conservation, and implies that JΣΩ = constant (cf. Ogilvie 2001). If the orbital velocity divergence is non-zero because of an eccentricity gradient, then the surface density varies periodically around the orbit.
In terms of the one-dimensional mass density introduced in Section 2.5, we have JΣΩ = M/P, where P is the orbital period. Sinceρ ∝ (JΩH) −1 andp ∝ρ γ ∝ (JΩH) −γ , the last term in equation (142) is ∝ (JΩH) −(γ−1) H −2 . When w is eliminated between equations (141) and (142) we obtain which describes a nonlinear vertical oscillator forced by the periodically varying orbital geometry. We are mainly interested in solutions that have period P in τ (or period 2π in θ), which are stationary on the orbital timescale when viewed in a non-rotating frame of reference. The general solution, however, includes a free oscillation that can only be eliminated by an appropriate choice of initial condition, or by including some dissipation.
In fact it is straightforward to include a bulk viscosity in the description of laminar flows. If the dynamic bulk viscosity is parametrized as µ b = α b p(GM/λ 3 ) −1/2 , as in Ogilvie (2001), and α b is independent of ζ, then equation (142) is modified tȯ We neglect the effects of viscous heating. [Ogilvie (2001) included shear and bulk viscosity (allowing for a non-zero relaxation time), viscous heating and radiative cooling.]

Linear theory for small eccentricity
A linear theory can be developed when both E and λE are small compared to unity, in which case the eccentric disc can be regarded as a small perturbation of a circular disc in which λ is the radial coordinate. In the circular case, the local model is hydrostatic with H = H0 = constant, w = 0,ρ =ρ0 = constant andp =p0 = constant, such that p0/ρ0H 2 0 = Φ2 = GM/λ 3 . In the presence of a small eccentricity and a small eccentricity gradient, such that e and λe are O( ) with 1 being a small parameter (different from that used previously), we have (from Appendix B) The laminar solution is of the form Bearing in mind that ∂τ = Ω∂ θ with Ω = (GM/λ 3 ) 1/2 [1 + O( )], in order to satisfy equations (141)-(143) we require the dimensionless perturbations to satisfy where we have allowed for a bulk viscosity as described above. The solution is We can similarly define a dimensionless surface density perturbatioñ Σ =ρ +H = λE .
We now refer to the global analysis of Section 2.5 and calculate the evolution of eccentricity associated with the laminar solution. The stress tensor of the laminar flow is and its vertical integral is In the linear theory developed above this becomes where the final terms in each case are due to the azimuthal variation of the metric coefficients (Appendix B). The required stress integrals are therefore, correct to O( ), To the same level of approximation, JρΩ 2 z 2 e iφ dφ dz = πλp0H0(2H + 2E).
We can now apply these results to the evolutionary equation (59) for E, evaluating it correct to first order. Note that there is no angular momentum transport or quasi-radial mass flux to this order because of the absence of shear viscosity. In applying our local results to the global disc, we write Σ and P for vertically integrated density and pressure in the circular global disc, and allow for the λ-dependence of these and other quantities. We obtain Using the above results this simplifies to As can be seen from the expression (157) forH, this equation contains γ and α b only in the combination γ − iα b (as was also found in the two-dimensional linear theory of Goodchild & Ogilvie 2006). In the inviscid case α b = 0 it simplifies to which we have derived independently by a three-dimensional linear perturbation analysis of a circular disc. For comparison, the two-dimensional linear theory of Goodchild & Ogilvie (2006) gives instead Differences between the two-and three-dimensional theories were found to be crucial for the dynamics of eccentric discs around Be stars by Ogilvie (2008). Equation (176) is a dispersive wave equation related to the Schrödinger equation, and indicates how eccentricity propagates through a disc by means of pressure. Bulk viscosity can easily be included by replacing γ with γ − iα b ; this gives an eccentricity diffusion coefficient of which, for α b γ, is half the (mass-weighted mean) kinematic bulk viscosity.

Behaviour for larger eccentricity
The linear theory shows that the periodic variation around the orbit of the vertical gravity coefficient Φ2 and the orbital velocity divergence ∆ induce a dynamical vertical oscillation of the disc. For a disc that behaves isothermally (γ = 1) the fractional oscillation amplitude of the scaleheight in linear theory is three times the eccentricity. It is clear, then, that in this case the oscillation will become strongly nonlinear at eccentricities well below unity.
We have computed the periodic solutions of the ordinary differential equations describing the nonlinear vertical oscillations. The left panel of Fig. 4 shows azimuthal profiles of H for eccentricities up to 0.6, by which point an extreme compression has occurred near pericentre. The right panel illustrates an almost universal behaviour near pericentre for larger eccentricities, when the variables are rescaled in terms of the minimum value of H, which is absurdly small in the case e = 0.75. In this regime vertical gravity is unimportant near pericentre; its variation around the orbit does however induce a dynamical collapse that bounces near pericentre because of the large pressure that develops there. Fig. 5 shows how the minimum and maximum values of H, which are obtained at pericentre and apocentre respectively, depend on e for various γ in the absence of an eccentricity gradient. The extreme behaviour occurs at larger e when γ is larger; this can be understood as an extension of the linear result thatH = −3E/γ when E = 0.
Finally, Fig. 6 shows some laminar flows in discs with an eccentricity gradient but with a circular reference orbit. Here the vertical oscillation is driven only by the orbital velocity divergence and occurs only if γ = 1. The behaviour is qualitatively different from that in Fig. 4 and the oscillations are of more modest amplitude.

CONCLUSION
Although Keplerian discs are usually assumed to be circular, the general Keplerian disc is elliptical, in accordance with Kepler's first law. As found by Syer & Clarke (1992) and other authors, viscosity (and other dissipative effects) do not necessarily lead to the circularization of a disc, even though circular orbits have the least energy for a given angular momentum; this is because dissipative forces can tap the reservoir of orbital energy by causing mass redistribution. Eccentricity may result from the initial conditions of the disc (as in the case of a disc formed through tidal disruption of a body on an elliptical orbit), from secular forcing by a companion with an elliptical orbit, from resonant forcing by a companion with a circular or elliptical orbit, or from instability of a circular disc.
In this paper we have revisited the theory of eccentric discs developed by Ogilvie (2001). In particular, we have formulated a local model, which is a generalization of the well known shearing sheet (or box) to the geometry of an eccentric disc. We have discussed the simplest hydrodynamic solutions in the local model, which are necessarily non-hydrostatic and involve a vertical oscillation at the orbital period. These oscillations can become highly nonlinear and exhibit extreme behaviour at eccentricities significantly less than unity, especially if the disc behaves isothermally. It would be valuable to determine, using numerical simulations, whether these extreme solutions are realized in practice. We have also computed the stresses associated with the laminar flows in a linear regime and derived the associated global evolutionary equation for the eccentricity, which differs significantly from a two-dimensional theory that neglects the vertical structure and oscillation of the disc.
A question not addressed in this paper is a possible vertical dependence of the eccentricity. In the absence of Figure 6. Laminar flows for eccentric discs with γ = 5/3 and with eccentricity gradients λe = 0.3, 0.6 and 0.9, but with a circular reference orbit (e = 0) and no twist (ω = 0). For λe = 0.3 the scaleheight varies sinusoidally with a maximum at pericentre. For larger eccentricity gradients the profile develops a local minimum at pericentre and two maxima elsewhere, but the behaviour is not extreme.
viscosity, turbulence and magnetic fields, layers of the disc at different heights are relatively weakly coupled by pressure gradients and can undergo independent epicyclic oscillations to some extent. As discussed by Latter & Ogilvie (2006), this allows the eccentricity to propagate radially with a nontrivial vertical profile. There may therefore be a transition in behaviour when the viscous, turbulent or magnetic stresses are very small.
On the right-hand sides of these equations, quantities such as λ, e, e and ω are to be evaluated on the reference orbit λ = λ0. Note that the combinations of metric and connection components with various powers of λ make these quantities dimensionless. Note also that the orbital velocity divergence ∆ is non-zero when the eccentricity gradient E is non-zero. It is possible, in principle, to write these coefficients in terms of τ rather than θ. These quantities can be related through Kepler's equation and the mean and eccentric anomalies, or by inverting the function ϕ(τ ). However, in practice, it is easier to use θ instead of τ as a time-like variable. Time-derivatives of these quantities at the reference point may be evaluated using the rule ∂τ = Ω∂ θ .