Local and global dynamics of eccentric astrophysical discs

Mon. Not. R. Astron. Soc. 000, 000–000 (0000) Printed 9 May 2018 (MN LATEX style file v2.2) Local and global dynamics of eccentric astrophysical discs Gordon I. Ogilvie and Adrian J. Barker Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA arXiv:1409.6487v1 [astro-ph.EP] 23 Sep 2014 9 May 2018 ABSTRACT 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. Key words: accretion, accretion discs – hydrodynamics – celestial mechanics 1 1.1 INTRODUCTION 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 c 0000 RAS 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. 2 Gordon I. Ogilvie and Adrian J. Barker [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 (Kato & Osaki 2013; Osaki & Kato 2013). 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). 1.2 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). Small eccentricities are governed by linear equations which can be derived through a perturbation analysis of a circular disc. This approach has been taken by, e.g., Kato (1983), Lee & Goodman (1999), Tremaine (2001), Papaloizou (2002), Goodchild & Ogilvie (2006) and Ogilvie (2008). 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 1996, 1998, 2000; Smith 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. 2009, 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. 1.3 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. 2 2.1 LARGE-SCALE GEOMETRY AND DYNAMICS 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 − e2 ), 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 eiω . (Note that the eccentricity, e, and the base of natural logarithms, e, are distinguished typographic 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs 3 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 ω). 2.2 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, r = R(λ, φ) = λ . 1 + e(λ) cos[φ − ω(λ)] (1) The semi-latus rectum λ can be thought of as a quasi-radial coordinate that replaces r, while φ is the usual azimuthal angular coordinate. c 0000 RAS, MNRAS 000, 000–000 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 gλλ = Rλ2 , g λλ = gλφ = gφλ = Rλ Rφ , R2 + Rφ2 , R2 Rλ2 g λφ = g φλ = − gφφ = R2 + Rφ2 , (2) Rφ , R 2 Rλ g φφ = 1 ,(3) R2 4 Gordon I. Ogilvie and Adrian J. Barker 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. 2.3 Hydrodynamic equations We consider an ideal fluid satisfying the equation of motion, (∂t + uj ∇j )ui = −g ij  ∇j Φ +  1 ∇j p , ρ (8) the equation of mass conservation, ∂t ρ + ∇i (ρui ) = 0, (9) and the thermal energy equation, Figure 2. Orbital coordinates (λ, φ) in an example of an eccentric disc. The semi-latus rectum λ is used to label the orbits and has the role of a quasi-radial coordinate. The second coordinate is the usual azimuthal angle φ. where the subscripts on R denote partial derivatives of the function R(λ, φ). The components of the Levi-Civita connection (Christoffel symbols) are Γλλλ = Rλλ , Rλ Γλλφ = Γλφλ = Γλφφ = − (R2 + 2Rφ2 − RRφφ ) , RRλ Γφλλ = 0, Γφλφ = Γφφλ = Rλφ Rφ − , Rλ R Rλ , R (4) (5) Γφφφ = 2Rφ . R (6) 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 Ds = 0, (10) where ui 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 + ui ∂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 ui , which are simply the rates of change of the orbital coordinates of fluid elements. Equation (8) can be written using partial derivatives as   1 Dui + Γijk uj uk = −g ij ∂j Φ + ∂j p , (11) ρ while equation (9) can be written either as ρ Dρ = − ∂i (Jui ) J (12) or in the conservative form ∂t (Jρ) + ∂i (Jρui ) = 0. (13) An alternative to equation (10) is ∂(x, y, z) = (det gij )1/2 = RRλ . J= ∂(λ, φ, z) (7) 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. Dp = − γp ∂i (Jui ), J (14) where γ = (∂ ln p/∂ ln ρ)s is the adiabatic index. The conservative form of the equation for the total energy of the fluid is ∂t (JρEtot ) + ∂i [J(ρEtot + p)ui ] = 0, (15) where Etot = 21 gij ui uj + ε + Φ (16) 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 = gil Γljk +glj Γlik are used, provided that the gravitational potential is independent of t. In detail, equations (11), (12) and (14) are Duλ + Γλλλ (uλ )2 + 2Γλλφ uλ uφ + Γλφφ (uφ )2     1 1 = −g λλ ∂λ Φ + ∂λ p − g λφ ∂φ Φ + ∂φ p , ρ ρ (17) c 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs Duφ + 2Γφλφ uλ uφ + Γφφφ (uφ )2     1 1 = −g λφ ∂λ Φ + ∂λ p − g φφ ∂φ Φ + ∂φ p , ρ ρ Duz = −∂z Φ − 1 ∂z p, ρ (18) (19)  1 1 Dρ = −ρ ∂λ (Juλ ) + ∂φ (Juφ ) + ∂z uz J J   1 1 λ φ z ∂λ (Ju ) + ∂φ (Ju ) + ∂z u , Dp = −γp J J λ φ (22) Orbital motion The orbital motion corresponds to the velocity field ui = ω 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. Γλφφ Ω2 = −g λλ ∂λ Φ0 − g λφ ∂φ Φ0 , (23) Ω∂φ Ω + Γφφφ Ω2 = −g λφ ∂λ Φ0 − g φφ ∂φ Φ0 . (24) Since Φ0 is a function of R only, these equations simplify to dΦ0 , dR ∂φ (R2 Ω) = 0, (25) (26) and are satisfied, as expected, when Φ0 = −GM/R and R2 Ω = ℓ = (GM λ)1/2 , i.e.  1/2 GM Ω= [1 + e cos(φ − ω)]2 . (27) λ3 2.5 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 by 2 r̈ − ℓ GM = − 2 + fr , r3 r ℓ̇ = rfφ , c 0000 RAS, MNRAS 000, 000–000 ṙ = ℓ e sin θ, λ (30) e e−iθ = e cos θ − ie sin θ = (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.) R 2 Ω2 = λ λ , 1 + e cos θ (21) with D = ∂t + u ∂λ + u ∂φ + u ∂z . The alternative, conservative form of equation (20) is 2.4 r= 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 z ∂t (Jρ) + ∂λ (Jρuλ ) + ∂φ (Jρuφ ) + ∂z (Jρuz ) = 0. where ℓ = r2 φ̇ 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 (20) and  5 (28) (29) λ iλṙ −1− r ℓ (31) and so E = e eiω =  iλṙ λ −1− r ℓ  eiφ , (32) which therefore evolves according to ℓĖ = rfφ (eiφ + E) + λ(fφ − ifr ) eiφ . (33) 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 = Rλ f λ + Rφ f φ and fφ = Rf φ , then we can write (using equation A5) ℓ̇ = R2 f φ , 2 φ ℓĖ = 2R f (e (34) iφ λ iφ + E) − iλRλ f e . (35) 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) ∂t (Jρ) + ∂λ [Jρ(λ̇ + uλ )] + ∂φ (Jρuφ ) + ∂z (Jρuz ) = 0, (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 ∂t M + ∂λ F = 0, (37) where ZZ Z M= Jρ dφ dz = JΣ dφ (38) Ris the one-dimensional mass density with respect to λ, Σ = ρ dz being the surface density, and ZZ F= Jρ(λ̇ + uλ ) dφ dz (39) is the quasi-radial mass flux. Note that the mass of the disc is ZZZ Z Jρ dλ dφ dz = M dλ, (40) where the integral is carried out over an appropriate range of λ. We can also write ∂t M + ∂λ (Mv̄ λ ) = 0, λ where v̄ = F/M is the mean quasi-radial velocity. (41) 6 Gordon I. Ogilvie and Adrian J. Barker Given that ℓ is independent of φ and z, the angular momentum equation has the one-dimensional form ZZ ∂t (Mℓ) + ∂λ (F ℓ) = JρR2 f φ dφ dz, (42) Using the relation (A9) we obtain ZZ  ℓM(∂t + v̄ λ ∂λ )E = 2(eiφ + E)∂λ (JR2 T λφ ) −iλ eiφ ∂λ (JRλ T λλ ) − i eiφ JR2 T φφ  JR2 iφ (e + E − λE ′ )T λφ dφ dz. − λ or, equivalently, M(∂t + v̄ λ ∂λ )ℓ = ZZ JρR2 f φ dφ dz, (43) 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 = ∇j T ij , where T ij is a symmetric stress tensor describing the collective effects of the disc (pressure, viscosity, self-gravity, etc.). Then (Ogilvie 2001)   1 R2 JRλ2 λφ λλ λ ∂λ (JRλ T ) + ∂φ T ρf = JRλ JRλ2 R2 − R2 φφ T + ∂z T λz , λRλ (44) 1 1 ρf = ∂λ (JR2 T λφ ) + ∂φ (JR2 T φφ ) + ∂z T φz . (45) JR2 JR2 φ The terms in JρR2 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 ZZ ∂t (Mℓ) + ∂λ (F ℓ) = ∂λ (JR2 T λφ ) dφ dz, (46) which can be written in the conservative form ∂t (Mℓ) + ∂λ (F ℓ + G) = 0, (47) where G=− ZZ JR2 T λφ dφ dz (48) is the internal torque. With the help of the equation of mass conservation (37), and the fact that ℓ depends only on λ, it simplifies to ∂G dℓ + = 0, (49) dλ ∂λ which determines the mass flux F (or the mean quasi-radial velocity v̄ λ = 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 F ℓM(∂t + v̄ λ ∂λ )E ZZ h i = Jρ 2R2 f φ (eiφ + E) − iλRλ f λ eiφ dφ dz. When ρf i = ∇j T ij , the terms involving ∂z again integrate to zero. The terms involving ∂φ are less obvious and an integration by parts is needed to obtain ZZ  2(eiφ + E)∂λ (JR2 T λφ ) ℓM(∂t + v̄ λ ∂λ )E = −iλ eiφ ∂λ (JRλ T λλ ) − i eiφ JR2 T φφ  2  JR2 R iφ +iλ 2λ T λφ ∂φ e dφ dz. R Rλ 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 ZZ JR2 T λφ dφ dz, (53) ZZ ZZ ZZ JR2 T λφ eiφ dφ dz, (54) JRλ T λλ eiφ dφ dz, (55) JR2 T φφ eiφ dφ dz. (56) 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 ∂t M + ∂λ (Mv̄ λ ) = 0, (57) ZZ ZZ dℓ Mv̄ λ = ∂λ (JR2 T λφ ) dφ dz + JρR2 f φ dφ dz,(58) dλ ℓM(∂t + v̄ λ ∂λ )E = ZZ  2(eiφ + E)∂λ (JR2 T λφ ) −iλ eiφ ∂λ (JRλ T λλ ) − i eiφ JR2 T φφ  JR2 iφ − (e + E − λE ′ )T λφ dφ dz λ ZZ h i + Jρ 2R2 f φ (eiφ + E) − iλRλ f λ eiφ dφ dz. (59) 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 (R2 + z 2 )−1/2 is expanded in a Taylor series about the midplane z = 0 of a thin disc, we obtain Φ = Φ0 (R) + Φ2 (R) 21 z 2 + O(z 4 ), (50) (51) (52) (60) where Φ0 = −GM/R is the potential in the midplane, while Φ2 = GM/R3 has a different dependence on R. Its contribution to f λ is fλ = − 1 3 GM z 2 ∂R (Φ2 21 z 2 ) = . Rλ 2 R 4 Rλ Its contribution to ℓM(∂t + v̄ λ ∂λ )E is therefore ZZ   Jρ −iλRλ f λ eiφ dφ dz  Z Z 3i =− Ω2 eiφ Jρz 2 dz dφ, 2 (61) (62) c 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs 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. 3 LOCAL MODEL OF AN ECCENTRIC DISC 3.1 Flow decomposition We have seen in Section 2.4 that the eccentric orbital motion with uφ = Ω(λ, φ) and uλ = uz = 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 : uλ = v λ , uφ = Ω + v φ , uz = v z . (63) The residual parts of the hydrodynamic equations (17)–(21) are then, without approximation, Dv λ + Γλλλ (v λ )2 + 2Γλλφ v λ (Ω + v φ ) + Γλφφ (2Ω + v φ )v φ   1 = −g λλ ∂λ (Φ − Φ0 ) + ∂λ p ρ   1 λφ (64) −g ∂φ (Φ − Φ0 ) + ∂φ p , ρ Dv φ + (v λ ∂λ + v φ ∂φ )Ω + 2Γφλφ v λ (Ω + v φ ) + Γφφφ (2Ω + v φ )v φ   1 = −g λφ ∂λ (Φ − Φ0 ) + ∂λ p ρ   1 φφ −g ∂φ (Φ − Φ0 ) + ∂φ p , (65) ρ Dv z = −∂z (Φ − Φ0 ) − 1 ∂z p, ρ (66)   1 1 Dρ = −ρ ∆ + ∂λ (Jv λ ) + ∂φ (Jv φ ) + ∂z v z , J J   1 1 λ φ z Dp = −γp ∆ + ∂λ (Jv ) + ∂φ (Jv ) + ∂z v , J J (67) (68) where 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  1  λλ g ∂λ p + g λφ ∂φ p ,(71) Dv λ + 2Γλλφ Ωv λ + 2Γλφφ Ωv φ = − ρ Dv φ + (v λ ∂λ + v φ ∂φ )Ω + 2Γφλφ Ωv λ + 2Γφφφ Ωv φ  1  λφ =− g ∂λ p + g φφ ∂φ p , ρ Dv z = −Φ2 z − φ z D = ∂t + v ∂λ + (Ω + v )∂φ + v ∂z (69) is the Lagrangian derivative and 1 ∂φ (JΩ) (70) J 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 ∆= c 0000 RAS, MNRAS 000, 000–000 1 ∂z p, ρ λ (72) (73) Dρ = −ρ(∆ + ∂λ v λ + ∂φ v φ + ∂z v z ), φ (74) z Dp = −γp(∆ + ∂λ v + ∂φ v + ∂z v ), (75) 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. 3.2 Local approximation We now select a reference orbit λ = λ0 with angular velocity Ω(λ0 , φ) (Fig. 3). We consider a reference point that follows this orbit, starting from the pericentre φ = ω0 at t = 0. Let ϕ(t) be the solution of dϕ(t)/dt = Ω(λ0 , ϕ(t)) subject to the initial condition ϕ(0) = ω0 . Then the orbital coordinates of the reference point are (λ, φ, z) = (λ0 , ϕ(t), 0). Let P0 be the period of the reference orbit, such that ϕ(P0 ) = ϕ(0) = ω0 . We examine the neighbourhood of the reference point by letting λ = λ0 + ξ, λ 7 φ = ϕ(t) + η, z = ζ, t = τ, (76) 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, 8 Gordon I. Ogilvie and Adrian J. Barker which makes them periodic functions of τ . We have also dropped the subscript zeros, so that D = ∂τ + v ξ ∂ξ + (Ωλ ξ + Ωφ η + v η )∂η + v ζ ∂ζ . (84) The thermal energy equation can be written as either Ds = 0 (85) or Dp = −γp(∆ + ∂ξ v ξ + ∂η v η + ∂ζ v ζ ). Figure 3. Construction of a local model of an eccentric disc. 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 Ω = Ω0 (τ ) + Ωλ0 (τ )ξ + Ωφ0 (τ )η + O(ǫ2 ), (77) 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 ∂t + Ω0 ∂φ = ∂τ . (78) D = ∂τ + v ξ ∂ξ + (Ωλ0 ξ + Ωφ0 η + v η )∂η + v ζ ∂ζ . (79) Equations (71)–(74) are then rewritten as  1  λλ g ∂ξ p + g λφ ∂η p , (80) Dv ξ + 2Γλλφ Ωv ξ + 2Γλφφ Ωv η = − ρ Dv η + (Ωλ + 2Γφλφ Ω)v ξ + (Ωφ + 2Γφφφ Ω)v η  1  λφ g ∂ξ p + g φφ ∂η p , =− ρ Dv ζ = −Φ2 ζ − 1 ∂ζ p, ρ Dρ = −ρ(∆ + ∂ξ v ξ + ∂η v η + ∂ζ v ζ ), 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 τ satisfying ġij = Ω(gik Γkjφ + gkj Γkiφ ), 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 ∂τ (Jρ)+∂ξ (Jρv ξ )+∂η [Jρ(Ωλ ξ+Ωφ η+v η )]+∂ζ (Jρv ζ ) = 0.(87) The Lagrangian derivative of the specific kinetic energy of the relative motion can be shown to be 1 D( 21 gij v i v j ) = −v i vj ∇i ω j − Φ2 ζv ζ − v i ∂i p. (88) ρ The first source term on the right-hand side involves the covariant derivative ∇i ω j = ∂i ω j + Γjik ω k of the orbital velocity field. Here vi = gij v j are the covariant components of the relative velocity. The conservative form of the energy equation for the relative motion is ∂τ (JρErel ) + ∂ξ [J(ρErel + p)v ξ ] +∂η [JρErel (Ωλ ξ + Ωφ η + v η ) + Jpv η ] +∂ζ [J(ρErel + p)v ζ ] = Jρ(−v i vj ∇i ω j + Φ̇2 12 ζ 2 ) − Jp∆, (81) (82) (83) in which, as explained above, the geometrical coefficients are evaluated at the reference point (λ, φ, z) = (λ0 , ϕ(τ ), 0), (89) where Erel = 12 gij v i v j + ε + 21 Φ2 ζ 2 . At the leading order in ǫ, therefore, (86) (90) 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. 3.3 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 ξ ′ = ξ, η ′ = α(τ )η + β(τ )ξ, ζ ′ = ζ, τ ′ = τ, (91) c 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs 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 ∂ξ = ∂ξ′ ∂τ = ∂τ′ + β∂η′ , ∂η = + (α̇η + α∂η′ , ∂ζ = ∂ζ′ , β̇ξ)∂η′ . β̇ = −αΩλ . (92) (93) Ds = 0 (105) (94) (95) 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 α = λ(1 + e cos θ)−2 , β= 3 2  1+ 2eλe′ 1 − e2  (96) GM λ3 1/2 τ λe′ (2 + e cos θ) sin θ λω ′ − (1 − e2 )(1 + e cos θ)2 (1 + e cos θ)2 +constant. − (97) The relation between the time τ (with origin τ = 0 at pericentre) and θ is given by Kepler’s equation in the form nτ = E − e sin E, (98) where n = (GM/a ) is the mean motion, a = λ/(1 − e2 ) is the semi-major axis and E is the eccentric anomaly, such that sin E = 3 1/2 (1 − e2 )1/2 sin θ . 1 + e cos θ (99) Thus τ can be expressed in terms of θ by inverting this sine function. The explicit expression for the Jacobian of the shearing coordinates is J = 1 + (e − λe′ ) cos θ − λeω ′ sin θ . 1 + e cos θ (100) The equations of the local model now have the form ξ Dv + 2Γλλφ Ωv ξ 2Γλφφ Ωv η + i 1 h λλ ′ = − g (∂ξ + β∂η′ )p + g λφ α∂η′ p , ρ Dv η + (Ωλ + 2Γφλφ Ω)v ξ + (Ωφ + 2Γφφφ Ω)v η i 1h = − g λφ (∂ξ′ + β∂η′ )p + g φφ α∂η′ p , ρ c 0000 RAS, MNRAS 000, 000–000 (103) (104) provided that α and β are chosen such that α̇ = −αΩφ , 1 ′ ∂ζ p, ρ h i Dρ = −ρ ∆ + (∂ξ′ + β∂η′ )v ξ + α∂η′ v η + ∂ζ′ v ζ . 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 D = ∂τ′ + v ξ (∂ξ′ + β∂η′ ) + αv η ∂η′ + v ζ ∂ζ′ Dv ζ = −Φ2 ζ ′ − 9 (101) The thermal energy equation can be written as either or i h Dp = −γp ∆ + (∂ξ′ + β∂η′ )v ξ + α∂η′ v η + ∂ζ′ v ζ . 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 ∂τ′ (J ρ)+∂ξ′ (J ρv ξ )+∂η′ [J ρ(αv η +βv ξ )]+∂ζ′ (J ρv ζ ) = 0.(107) Subject to periodic boundary conditions in ξ ′ and η ′ , and suitable boundary conditions in ζ ′ , the conserved mass in the shearing box is Z ZZZ dM = J ρ dξ ′ dη ′ dζ ′ . (108) R The horizontal momentum components P ξ = v ξ dM R and P η = v η dM of the box satisfy the equations ∂τ′ P ξ + 2Γλλφ ΩP ξ + 2Γλφφ ΩP η = 0, ∂τ′ P η (102) (106) + (Ωλ + 2Γφλφ Ω)P ξ + (Ωφ + 2Γφφφ Ω)P η (109) = 0, (110) which allow an epicyclic oscillation around the reference orbit, but it would be natural to constrain the solutions to 10 Gordon I. Ogilvie and Adrian J. Barker satisfy P ξ = P η = 0, meaning that the reference orbit has been correctly defined. The energy equation in shearing coordinates has the form Φ̇2 12 ζ 2 ) − J p∆. (111) 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, eiφ and Rλ are known functions of time in the local model. Since the local modelRfollows an orbiting reference point, the azimuthal integral . . . dφ can be interpreted as R 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. 3.4 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/R3 )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 Dv ξ + 2∆v ξ − Dv η + 2RΩ η 1 v =− ∂ξ p, J ρJ 2 1 ∂ζ p, ρ Dρ = −ρ(∆ + ∂ξ v ξ + ∂η v η + ∂ζ v ζ ),  D = ∂τ + v ξ ∂ξ + ( 21 − 2J ) J˙ , J  J = J = Rλ = 1 − λe′ cos θ R Ω= x = J ξ, t = τ, (124) so that ∂x = J −1 ∂ξ , ∂t = ∂τ − ∆ξ∂ξ , and the velocity components are related by v x = J v ξ + J˙ ξ, v y = Rv η + 12 (1 − J )Ωξ. ξ (126) η 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. 4 4.1 LAMINAR FLOWS Nonlinear vertical oscillations Dv ζ = −Φ2 ζ − (115) Dρ = −ρ(∆ + ∂ζ v ζ ), 1 ∂ζ p, ρ Dp = −γp(∆ + ∂ζ v ),  (125) (114) (127) (128) (129) ζ (116) (117) (118) and θ = Ωt. These equations can be derived from the standard hydrodynamic equations of the shearing sheet, 1 Dv x − 2Ωv y = − ∂x p, ρ (119) 1 x 1 Ωv = − ∂y p, 2 ρ (120) Dv y + (Note that, since these are Cartesian coordinates, v 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′ 6= 0 then a time-dependent homogeneous transformation of the radial coordinate is involved (cf. Latter & Ogilvie 2009, Appendix A): ζ Ωξ + v η ∂η + v ζ ∂ζ , R ∆=  (123) (113) with λe′ sin θ 1 − λe′ cos θ (122) 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 (112) 1 Ω ξ v = − 2 ∂η p, 2R ρR Dv ζ = −Ω2 ζ − Dρ = −ρ(∂x v x + ∂y v y + ∂z v z ), x +∂ζ [J (ρErel + p)v ζ ] = J ρ(−v vj ∇i ω + (121) D = ∂t + v x ∂x + (− 32 Ωx + v y )∂y + v z ∂z . +∂η′ [J (ρErel + p)(αv η + βv ξ )] j 1 ∂z p, ρ with ∂τ′ (J ρErel ) + ∂ξ′ [J (ρErel + p)v ξ ] i Dv z = −Ω2 z − with D = ∂τ + v ∂ζ . The periodic variation of the vertical gravity coefficient Φ2 around the orbit (in the case E 6= 0) and the orbital velocity divergence ∆ (in the case E ′ 6= 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 Fp′ (x) = −xFρ (x). (130) c 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs p̂˙ ρ̂˙ = = −(∆ + w). ρ̂ γ p̂ We normalize the profiles such that Z ∞ Fρ (x) dx = 1 (131) and Z ∞ (132) −∞ Fp (x) dx = −∞ Z ∞ Fρ (x)x2 dx = 1. −∞ (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,   x2 , (133) Fρ (x) = Fp (x) = (2π)−1/2 exp − 2 the homogeneous structure, 1 Fρ = √ , 2 3 Fp = (134) 3 − x2 √ 4 3 (135) (for x2 < 3 only), and the polytropic structure,  n x2 Fρ (x) = Cn 1 − , 2n + 3  n+1 x2 2n + 3 Cn 1 − , Fp (x) = 2(n + 1) 2n + 3 (136) (137) 2 (for x < 2n + 3 only) where n > 0 (not necessarily an integer) is the polytropic index and Cn = [(2n + 3)π]−1/2 Γ(n + 32 ) Γ(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, v ζ = w(τ )ζ, (138) ρ = ρ̂(τ )Fρ (x), (139) p = p̂(τ )Fp (x), (140) 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 trueR semi-thickness. In each case the surface density is Σ(τ ) = ρ dζ = Rρ̂(τ )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(τ ), ρ̂(τ ) and p̂(τ ) obey Ḣ = w, H ẇ + w2 = −Φ2 + (141) p̂ , ρ̂H 2 c 0000 RAS, MNRAS 000, 000–000 (142) 11 (143) 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 satisfies ρ̂˙ Ḣ J˙ Ω̇ Σ̇ = + = −∆ = − − . Σ ρ̂ H J Ω (144) 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, Z Z Z M = JΣ dφ = JΣΩ dτ = JΣΩ dτ, (145) we have JΣΩ = M/P, where P is the orbital period. Since ρ̂ ∝ (JΩH)−1 and p̂ ∝ ρ̂γ ∝ (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 Ḧ + Φ2 ∝ (JΩ)−(γ−1) H −(γ+1) , H (146) 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 to # " −1/2  p̂ GM 2 (∆ + w) .(147) ẇ + w = −Φ2 + 1 − αb λ3 ρ̂H 2 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.] 4.2 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 and p̂ = p̂0 = constant, such that p̂0 /ρ̂0 H02 = Φ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), 12 Gordon I. Ogilvie and Adrian J. Barker Figure 4. Left: Laminar flows for eccentric discs with eccentricities 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6, no eccentricity gradient and γ = 1. The vertical scaleheight, in units of the hydrostatic value for a circular disc with the same semi-latus rectum, is plotted versus the true anomaly. For small e the variation is approximately sinusoidal and agrees with the linear theory, but for larger e extreme behaviour occurs close to the pericentre where H approaches zero. Right: Behaviour near pericentre for e = 0.6, 0.65, 0.7 and 0.75 (inner to outer lines), with variables scaled by the corresponding minimum value of H, which is Hmin = 2.56 × 10−3 , 9.69 × 10−5 , 4.36 × 10−7 and 3.05 × 10−11 , respectively. we have (from Appendix B)   GM Φ2 = [1 + 3e cos θ + O(ǫ2 )] λ3     i GM h −iφ 2 = 1 + Re 3E e + O(ǫ ) , λ3 ∆ =  GM λ3 =  GM λ3 1/2 1/2 h   i Re iλE ′ e−iφ + O(ǫ2 ) . h   i ρ̂ = ρ̂0 1 + Re ρ̃ e−iφ + O(ǫ2 ) , h   i p̂ = p̂0 1 + Re p̃ e−iφ + O(ǫ2 ) ,  GM λ3 1/2 h H̃ = iw̃ = (148) ρ̃ =   i Re w̃ e−iφ + O(ǫ2 ) . (149) (150) (151) (152) (153) 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 −iH̃ = w̃, (154) −iw̃ = −3E + p̃ − ρ̃ − 2H̃ − αb (iλE ′ + w̃), (155) −iρ̃ = − ip̃ = −(iλE ′ + w̃), γ −3E + (γ − 1 − iαb )λE ′ , γ − iαb p̃ 3E + λE ′ = . γ γ − iαb (157) (158) We can similarly define a dimensionless surface density perturbation [λe′ sin θ − λeω ′ cos θ + O(ǫ2 )] The laminar solution is of the form h   i H = H0 1 + Re H̃ e−iφ + O(ǫ2 ) , w= where we have allowed for a bulk viscosity as described above. The solution is (156) Σ̃ = ρ̃ + H̃ = λE ′ . (159) Note that this dynamical solution differs significantly from the hydrostatic non-solution in which the acceleration (−iw̃) and viscous terms are neglected in equation (155): H̃ = iw̃ = [−3E + (γ − 1)λE ′ ]/(γ + 1) and ρ̃ = p̃/γ = (3E + 2λE ′ )/(γ + 1). In the absence of bulk viscosity, the amplitude with which H oscillates is larger by a factor of 1 + γ −1 in the dynamical solution than in the hydrostatic non-solution. 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 " # −1/2  GM ij T = −p 1 − αb (∆ + w) g ij (160) λ3 and its vertical integral is # "  −1/2 Z GM (∆ + w) g ij . T ij dz = −p̂H 1 − αb λ3 In the linear theory developed above this becomes Z h   i T λλ dz = −p̂0 H0 1 + Re T̃ λλ e−iφ + O(ǫ2 ) (161) (162) c 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs Z Z h   i λT λφ dz = −p̂0 H0 Re T̃ λφ e−iφ + O(ǫ2 ) , (163) h   i λ2 T φφ dz = −p̂0 H0 1 + Re T̃ φφ e−iφ + O(ǫ2 ) , (164) with T̃ λλ = p̃ + H̃ − αb (iλE ′ + w̃) + 2(E + λE ′ ), (165) T̃ λφ (166) T̃ φφ = −iE, ′ = p̃ + H̃ − αb (iλE + w̃) + 2E, (167) 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(ǫ), ZZ JR2 T λφ dφ dz = 0, (168) ZZ JR T ZZ   JRλ T λλ eiφ dφ dz = −πλp̂0 H0 T̃ λλ − 3E − 2λE ′ ,(170) ZZ 2 λφ iφ 2 e dφ dz = −πλ p̂0 H0 T̃ λφ ,   JR2 T φφ eiφ dφ dz = −πλp̂0 H0 T̃ φφ − 4E − λE ′ .(171) To the same level of approximation, ZZ M= Jρ dφ dz = 2πλρ̂0 H0 + O(ǫ2 ), ZZ (169) JρΩ2 z 2 eiφ dφ dz = πλp̂0 H0 (2H̃ + 2E). (172) (173) 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 2Σ(GM λ3 )1/2 ∂t E = −∂λ (2λ2 P T̃ λφ ) h i +iλ∂λ λP (T̃ λλ − 3E − 2λE ′ ) + iλP (T̃ φφ − 4E − λE ′ ) +λP T̃ λφ − 3iλP (H̃ + E). Using the above results this simplifies to h i 2Σ(GM λ3 )1/2 ∂t E = i ∂λ λ2 P (H̃ + 4E + λE ′ ) −iλP (3H̃ + 5E + λE ′ ). (174) (175) As can be seen from the expression (157) for H̃, 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 2Σ(GM λ3 )1/2 ∂t E = i ∂λ [(2 − γ −1 )P λ3 E ′ ] +i(4 − 3γ −1 )λ2 E∂λ P + 3i(1 + γ −1 )λP E, (176) 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 2Σ(GM λ3 )1/2 ∂t E = i ∂λ (γP λ3 E ′ ) + iλ2 E∂λ P. c 0000 RAS, MNRAS 000, 000–000 (177) 13 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   1/2  1 GM αb 2 H , (178) 2 γ 2 + αb2 λ3 which, for αb ≪ γ, is half the (mass-weighted mean) kinematic bulk viscosity. 4.3 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 that H̃ = −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 γ 6= 1. The behaviour is qualitatively different from that in Fig. 4 and the oscillations are of more modest amplitude. 4.4 Comparison with Ogilvie (2001) If heating and cooling are neglected and the stress is an instantaneous bulk viscous stress, then the relevant parts of equations (211)–(217) of Ogilvie (2001) reduce to (1 + e cos θ)2 ∂θ ln f2 = −(γ + 1)f3 − (γ − 1)g1 , (179) (1+e cos θ)2 ∂θ f3 = −f32 −(1+e cos θ)3 +f2 [1−αb (g1 +f3 )].(180) The correspondence with the equations of this section is f2 7→ p̂/ρ̂H 2 , f3 7→ w(GM/λ3 )−1/2 , g1 7→ ∆(GM/λ3 )−1/2 . 14 Gordon I. Ogilvie and Adrian J. Barker Figure 5. Minimum and maximum scaleheight, in units of the hydrostatic value for a circular disc with the same semi-latus rectum, for eccentric discs with no eccentricity gradient and with γ = 1, 1.2, 1.4, 1.6, 1.8 and 2. The outermost curves are for γ = 1 and the innermost ones for γ = 2. 5 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. 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 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. In the companion paper (Barker & Ogilvie 2014) we use the local model to analyse the linear hydrodynamic stability of an eccentric disc. In the absence of viscosity and magnetic fields, eccentric discs are susceptible to a hydrodynamic instability that excites internal (inertial) waves and may induce hydrodynamic turbulence. This is likely to be important for the evolution of the eccentricity, but also potentially for transport processes and mixing. ACKNOWLEDGEMENTS This research was supported by STFC through grants ST/J001570/1 and ST/L000636/1. REFERENCES Armitage P. J., Natarajan P., 2005, ApJ, 634, 921 Barker A. J., Ogilvie G. I., 2014, MNRAS, in press Borderies N., Goldreich P., Tremaine S., 1983, AJ, 88, 1560 D’Angelo G., Lubow S. H., Bate M. R., 2006, ApJ, 652, 1698 Goldreich P., Sari R., 2003, ApJ, 585, 1024 Goldreich P., Tremaine S., 1981, ApJ, 243, 1062 c 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs Goodchild S., Ogilvie G., 2006, MNRAS, 368, 1123 Goodman J., 1993, ApJ, 406, 596 Guillochon J., Ramirez-Ruiz E., Lin D., 2011, ApJ, 732, 74 Guillochon J., Loeb A., MacLeod M., Ramirez-Ruiz E., 2014, ApJ, 786, L12 Gurzadian V. G., Ozernoi L. M., 1979, Nature, 280, 214 Haswell C. A., King A. R., Murray J. R., Charles P. A., 2001, MNRAS, 321, 475 Hessman F. V., Mantel K.-H., Barwig H., Schoembs R., 1992, A&A, 263, 147 Kato S., 1978, MNRAS, 185, 629 Kato S., 1983, PASJ, 35, 249 Kato T., Osaki Y., 2013, PASJ, 65, 97 Kerswell R. R., 2002, Annu. Rev. Fluid Mech., 34, 83 Kley W., Dirksen G., 2006, A&A, 447, 369 Kley W., Papaloizou J. C. B., Ogilvie G. I., 2008, A&A, 487, 671 Latter H. N., Ogilvie G. I., 2006, MNRAS, 372, 1829 Latter H. N., Ogilvie G. I., 2009, Icarus, 202, 565 Lee E., Goodman J., 1999, MNRAS, 308, 984 Lubow S. H., 1991, ApJ, 381, 259 Lubow S. H., 1991, ApJ, 381, 268 Lyubarskij Y. E., Postnov K. A., Prokhorov M. E., 1994, MNRAS, 266, 583 Marzari F., Scholl H., Thébault P., Baruteau C., 2009, A&A, 508, 1493 Marzari F., Baruteau C., Scholl H., Thebault P., 2012, A&A, 539, A98 Mosqueira I., 1996, Icarus, 122, 128 Murray J. R., 1996, MNRAS, 279, 402 Murray J. R., 1998, MNRAS, 297, 323 Murray J. R., 2000, MNRAS, 314, L1 Ogilvie G. I., 2001, MNRAS, 325, 231 Ogilvie G. I., 2008, MNRAS, 388, 1372 Ogilvie G. I., Latter H. N., 2013, MNRAS, 433, 2403 Ogilvie G. I., Latter H. N., 2013, MNRAS, 433, 2420 Ogilvie G. I., Lubow S. H., 2003, ApJ, 587, 398 Osaki Y., Kato T., 2013, PASJ, 65, 50 Paardekooper S.-J., Thébault P., Mellema G., 2008, MNRAS, 386, 973 Papaloizou J. C. B., 2002, A&A, 388, 615 Papaloizou J. C. B., 2005, A&A, 432, 743 Papaloizou J. C. B., 2005, A&A, 432, 757 Peiris H. V., Tremaine S., 2003, ApJ, 599, 237 Rivinius T., Carciofi A. C., Martayan C., 2013, A&ARv, 21, 69 Smak J., 2009, Acta Astron., 59, 89 Smith A. J., Haswell C. A., Murray J. R., Truss M. R., Foulkes S. B., 2007, MNRAS, 378, 785 Syer D., Clarke C. J., 1992, MNRAS, 255, 92 Tremaine S., 1995, AJ, 110, 628 Tremaine S., 2001, AJ, 121, 1776 Warner B., 1995, Cataclysmic Variable Stars. Cambridge Univ. Press, Cambridge Whitehurst R., 1988, MNRAS, 232, 35 and therefore   λ 2 (∂φ + 1) = 1, R (A2) which implies R3 λ and simplifies the expression for R2 + 2Rφ2 − RRφφ = Γλφφ = − R2 . λRλ (A3) (A4) By differentiating equation (A1) once with respect to φ, we find λ (R − iRφ ) = 1 + E e−iφ , (A5) R2 which can also be written as ∂φ (R eiφ ) = iR2 iφ (e + E) λ (A6) or i −iφ (e + E e−2iφ ). (A7) λ Differentiating with respect to λ and interchanging the order of differentiation, we find   i Rλ i h −iφ ′ −2iφ = − ∂φ e + (E − λE ) e . (A8) R2 eiφ λ2 −∂φ (R eiφ )−1 = This in turn implies  2 iφ  iR4 R e = 2 2 (eiφ + E − λE ′ ). ∂φ Rλ λ Rλ (A9) APPENDIX B: EXPRESSIONS FOR THE COEFFICIENTS Let θ = ϕ(t) − ω be the true anomaly on the reference orbit. Then the following quantities, when evaluated at the reference point as described above, may be written explicitly in terms of θ: λ R= , (B1) 1 + e cos θ Rλ = 1 + (e − λe′ ) cos θ − λeω ′ sin θ , (1 + e cos θ)2 (B2) Rφ = λe sin θ , (1 + e cos θ)2 (B3) λ[1 + (e − λe′ ) cos θ − λeω ′ sin θ] , (1 + e cos θ)3  1/2 GM Ω= (1 + e cos θ)2 , λ3 J= λΩλ = APPENDIX A: GEOMETRICAL RELATIONS 15  GM λ3 1/2  − (B4) (B5) 3 (1 + e cos θ)2 2  +2(1 + e cos θ)(λe′ cos θ + λeω ′ sin θ) , (B6) From the polar equation (1) of an ellipse, we have λ = 1 + e cos(φ − ω) R c 0000 RAS, MNRAS 000, 000–000 (A1) Ωφ =  GM λ3 1/2 (−2e sin θ)(1 + e cos θ), (B7) 16 Gordon I. Ogilvie and Adrian J. Barker  GM (1 + e cos θ)3 , (B8) Φ2 = λ3  1/2 (1 + e cos θ)[λe′ sin θ − λeω ′ (cos θ + e)] GM ∆= ,(B9) 3 λ 1 + (e − λe′ ) cos θ − λeω ′ sin θ  gλλ = λ −1 [1 + (e − λe′ ) cos θ − λeω ′ sin θ]2 (1 + e cos θ)4 gλφ [1 + (e − λe′ ) cos θ − λeω ′ sin θ]e sin θ = , (1 + e cos θ)4 (B10) (B11) 2 1 + 2e cos θ + e , (1 + e cos θ)4 λ−2 gφφ = g λλ = (B12) (1 + e cos θ)2 (1 + 2e cos θ + e2 ) , [1 + (e − λe′ ) cos θ − λeω ′ sin θ]2 λg λφ = − (1 + e cos θ)2 e sin θ , 1 + (e − λe′ ) cos θ − λeω ′ sin θ λ2 g φφ = (1 + e cos θ)2 , Γλλφ 1 , 1 + (e − λe′ ) cos θ − λeω ′ sin θ ′ λΓφλφ (B16) (B17) ′ 1 + (e − λe ) cos θ − λeω sin θ , = 1 + e cos θ Γφφφ = (B14) (B15) λe′ sin θ − λeω ′ (cos θ + e) , = (1 + e cos θ)[1 + (e − λe′ ) cos θ − λeω ′ sin θ] λ−1 Γλφφ = − (B13) 2e sin θ . 1 + e cos θ (B18) (B19) 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 ∂τ = Ω∂θ . APPENDIX C: MAGNETOHYDRODYNAMIC EQUATIONS In ideal magnetohydrodynamics (MHD) the equation of motion can be written in the form    1 B2 j i ij (∂t + u ∇j )u = −g ∇j Φ + ∇j p + ρ 2µ0 1 j i (C1) + B ∇j B , µ0 and the induction equation in the form (∂t + uj ∇j )B i = B j ∇j ui − B i ∇j uj , (C2) while the solenoidal condition is i ∇i B = 0. (C3) In the local model using non-shearing coordinates, under similar scaling assumptions to those used in deriving the hydrodynamic equations, the three components of the equation of motion (80)–(82) are therefore modified to Dv ξ + 2Γλλφ Ωv ξ + 2Γλφφ Ωv η      B2 B2 1 λλ λφ + g ∂η p + g ∂ξ p + =− ρ 2µ0 2µ0 1 (B ξ ∂ξ + B η ∂η + B ζ ∂ζ )B ξ , + ρµ0 Dv η + (Ωλ + 2Γφλφ Ω)v ξ + (Ωφ + 2Γφφφ Ω)v η      B2 1 λφ B2 + g φφ ∂η p + =− g ∂ξ p + ρ 2µ0 2µ0 1 ξ η ζ η + (B ∂ξ + B ∂η + B ∂ζ )B , ρµ0 Dv ζ = −Φ2 ζ − + 1 ∂ζ ρ  p+ B2 2µ0 (C4) (C5)  1 (B ξ ∂ξ + B η ∂η + B ζ ∂ζ )B ζ , ρµ0 (C6) while the three components of the induction equation are DB ξ = (B ξ ∂ξ + B η ∂η + B ζ ∂ζ )v ξ −B ξ (∆ + ∂ξ v ξ + ∂η v η + ∂ζ v ζ ), (C7) DB η = Ωλ B ξ + Ωφ B η + (B ξ ∂ξ + B η ∂η + B ζ ∂ζ )v η −B η (∆ + ∂ξ v ξ + ∂η v η + ∂ζ v ζ ), (C8) DB ζ = (B ξ ∂ξ + B η ∂η + B ζ ∂ζ )v ζ −B ζ (∆ + ∂ξ v ξ + ∂η v η + ∂ζ v ζ ), (C9) and the solenoidal condition is ∂ξ B ξ + ∂η B η + ∂ζ B ζ = 0. (C10) Under the transformation to shearing coordinates these equations become Dv ξ + 2Γλλφ Ωv ξ + 2Γλφφ Ωv η      B2 B2 1 λλ ′ g (∂ξ + β∂η′ ) p + + g λφ α∂η′ p + =− ρ 2µ0 2µ0 1 ξ ′ ′ η ′ ζ ′ ξ + [B (∂ξ + β∂η ) + B α∂η + B ∂ζ ]B , (C11) ρµ0 Dv η + (Ωλ + 2Γφλφ Ω)v ξ + (Ωφ + 2Γφφφ Ω)v η      B2 B2 1 λφ ′ ′ φφ ′ g (∂ξ + β∂η ) p + + g α∂η p + =− ρ 2µ0 2µ0 1 + [B ξ (∂ξ′ + β∂η′ ) + B η α∂η′ + B ζ ∂ζ′ ]B η , (C12) ρµ0 Dv ζ = −Φ2 ζ − + 1 ′ ∂ζ ρ  p+ B2 2µ0  1 [B ξ (∂ξ′ + β∂η′ ) + B η α∂η′ + B ζ ∂ζ′ ]B ζ , ρµ0 (C13) DB ξ = [B ξ (∂ξ′ + β∂η′ ) + B η α∂η′ + B ζ ∂ζ′ ]v ξ −B ξ [∆ + (∂ξ′ + β∂η′ )v ξ + α∂η′ v η + ∂ζ′ v ζ ], (C14) DB η = Ωλ B ξ + Ωφ B η + [B ξ (∂ξ′ + β∂η′ ) + B η α∂η′ + B ζ ∂ζ′ ]v η −B η [∆ + (∂ξ′ + β∂η′ )v ξ + α∂η′ v η + ∂ζ′ v ζ ], (C15) c 0000 RAS, MNRAS 000, 000–000 Local and global dynamics of eccentric discs DB ζ = [B ξ (∂ξ′ + β∂η′ ) + B η α∂η′ + B ζ ∂ζ′ ]v ζ −B ζ [∆ + (∂ξ′ + β∂η′ )v ξ + α∂η′ v η + ∂ζ′ v ζ ], (∂ξ′ + β∂η′ )B ξ + α∂η′ B η + ∂ζ′ B ζ = 0. c 0000 RAS, MNRAS 000, 000–000 (C16) (C17) 17