{'fraction_non_alphanumeric': 0.03201431696162259, 'fraction_numerical': 0.01586796579836946, 'mean_word_length': 4.720737117506541, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear hydrodynamical evolution of eccentric Keplerian discs in two dimensions: validation of secular theory

We perform global two-dimensional hydrodynamical simulations of Keplerian discs with free eccentricity over thousands of orbital periods. Our aim is to determine the validity of secular theory in describing the evolution of eccentric discs, and to explore their nonlinear evolution for moderate eccentricities. Linear secular theory is found to correctly predict the structure and precession rates of discs with small eccentricities. However, discs with larger eccentricities (and eccentricity gradients) are observed to precess faster (retrograde relative to the orbital motion), at a rate that depends on their eccentricities (and eccentricity gradients). We derive analytically a nonlinear secular theory for eccentric gas discs, which explains this result as a modification of the pressure forces whenever eccentric orbits in a disc nearly intersect. This effect could be particularly important for highly eccentric discs produced in tidal disruption events, or for narrow gaseous rings; it might also play a role in causing some of the variability in superhump binary systems. In two dimensions, the eccentricity of a moderately eccentric disc is long-lived and persists throughout the duration of our simulations. Eccentric modes are however weakly damped by their interaction with non-axisymmetric spiral density waves (driven by the Papaloizou-Pringle instability, which occurs in our idealised setup with solid walls), as well as numerical diffusion.


INTRODUCTION
Eccentric gas discs are thought to arise in a number of astrophysical contexts. The orbital evolution of a newly born planet due to its tidal interaction with the protoplanetary disc is intricately coupled with the evolution of eccentric motions within the disc. But the importance of this interaction in producing the eccentricities of observed exoplanets is currently unclear (Papaloizou et al. 2001;Papaloizou 2002;Goldreich & Sari 2003;Kley & Dirksen 2006;D'Angelo et al. 2006;Bitsch et al. 2013). Eccentric gas discs are also thought to explain the superhump phenomenon in SU UMatype binary stars (Whitehurst 1988;Lubow 1991;Goodchild & Ogilvie 2006;Smith et al. 2007) and the spectral variability of rapidly rotating Be stars (Okazaki 1991;Papaloizou et al. 1992;Ogilvie 2008). In addition, eccentric gas discs are formed by the tidal disruption of stars or gaseous planets (Guillochon et al. 2011;Liu et al. 2013;, and such a process might be responsible for the gas cloud G2 near Sgr A* in the Galactic Centre Coughlin & Nixon 2015;Pfuhl et al. 2015).  spite this wide range of applications, the nonlinear dynamics of eccentric discs remains poorly understood. In the Solar system, several of the rings of Saturn and Uranus are observed to be elliptical, and a theory for narrow eccentric rings (with small eccentricities but arbitrary eccentricity gradients) composed of weakly collisional particles has been developed (Borderies et al. 1983;Chiang & Goldreich 2000). For small eccentricities, eccentric discs can be described as slowly precessing one-armed density waves (or shocks) (e.g. Okazaki 1991;Papaloizou et al. 1992;Lee & Goodman 1999;Goodchild & Ogilvie 2006;Ogilvie 2008;Saini et al. 2009). A particular solution for a global uniformly eccentric (and apsidally aligned) Keplerian disc has also been studied by Statler (2001).
A general secular theory for three-dimensional eccentric discs with arbitrary eccentricities and eccentricity gradients has been developed (Ogilvie 2001), and we aim to explore (and test using numerical simulations) some of its implications in this work. Secular theory describes dynamical processes that occur on timescales that are much longer than the orbital timescale. In the more familiar context of a system of multiple planets around a star, the Keplerian orbit of each planet is subject to perturbing forces that are pro-portional to the planet-to-star mass ratio; secular theory accounts for the long-term behaviour of eccentricity and inclination resulting from the orbit-averaged perturbing forces. For a thin gaseous disc around a star, however, pressure provides the most important departure from Keplerian orbital motion, being of second order in the aspect ratio H/R of the disc. Other subdominant forces, such as those due to the self-gravity of the disc (Tremaine 2001) and viscous or turbulent stresses (Ogilvie 2001), can also be important. In this situation, secular theory describes the long-term evolution of eccentricity due to gas pressure (and other forces) in thin discs, but neglects contributions that are of fourth order or higher in H/R.
In recent work, we developed a local model of an eccentric disc , which is similar to the shearing box that is commonly used to study the dynamics of circular astrophysical discs (Goldreich & Lynden-Bell 1965;Hawley et al. 1995). We used this model to analyse the vertical oscillatory flows that are driven by the variation in the vertical gravity around a non-circular orbit , and subsequently studied the local linear stability of these discs and their vertical flows . We found that eccentric discs are generically unstable (in three dimensions), being subject to a small-scale instability in which inertial waves are driven by a parametric resonance. This instability is expected to lead to wave activity or turbulence, and to damping of the disc eccentricity and eccentricity gradients (Papaloizou 2005b), unless these are maintained by external forcing (or additional instabilities e.g. viscous overstability). Nonlinear simulations are required to determine the efficiency of this damping process so that its importance in real discs can be quantified (e.g. continuing the work begun by Papaloizou 2005b, by taking into account the vertical structure of the disc).
In this paper we study the nonlinear hydrodynamical evolution of eccentric discs in two dimensions, deferring three-dimensional simulations to future work. Our aim is to determine the validity of linear secular theory (Ogilvie 2001;Papaloizou 2002; in describing the structure and precession rates of eccentric discs, as well as to determine their two-dimensional nonlinear dynamics. In the process, we will derive analytically, and compare with simulations, a two-dimensional nonlinear secular theory for isothermal eccentric discs that is valid for any eccentricity and eccentricity gradient, for a thin untwisted disc with non-intersecting orbits -this is a particular case of the general theory of Ogilvie (2001) that is amenable to analytical study. In two dimensions, eccentric discs do not exhibit the parametric instability or vertical oscillatory flows. However, a fundamental study of the two-dimensional nonlinear dynamics of eccentric discs has not yet been undertaken. We believe this to be worthwhile, in spite of the clear importance of three-dimensional effects in reality, because the two-dimensional dynamics that we will study are likely to also play a role in three dimensions. In particular, the force due to gas pressure, which tends to cause retrograde precession of these discs, can be captured in two dimensions. These simulations also provide a necessary benchmark to allow comparison with future three-dimensional simulations.
We outline our intentionally simplified numerical setup designed to study the dynamics of eccentric discs in § 2. We derive analytically a nonlinear secular theory for untwisted eccentric discs in two-dimensions in Appendix A. This theory is explored, and its predictions for the shapes and precession rates of eccentric discs are compared with simulations, in § 3. The long-term nonlinear evolution of eccentric discs is presented in § 4, where we also analyse the background instability that arises due to our idealised setup with rigid walls in the absence of any free eccentricity. We then finish with our conclusions.

Background disc
Our model consists of an eccentric (nearly) Keplerian disc in a two-dimensional cylindrical domain as an initial condition. We study its resulting nonlinear evolution using the PLUTO code (Mignone et al. 2007), adopting cylindrical polar coordinates (R, φ) with corresponding velocity components uR and u φ . Our governing equations are To allow us to more easily understand the outcome of our simulations, we consider a power-law disc which is isothermal (where the pressure p = c 2 s Σ, and Σ is the total surface density; with sound speed cs = const), with background disc surface density, for the circular case, where σ will be varied. (We expect the qualitative results of this paper to carry over to discs with different thermodynamic behaviour, and those with a radially varying sound speed.) The axisymmetric basic state of the disc has angular velocity Ω = Ω0R − 3 2 (which does not strictly apply when σ = 0 due to radial pressure gradients), where Ω0 is the Keplerian angular velocity at R = Ri, and u φ = RΩ. Our disc occupies the full azimuthal extent 1 φ ∈ [0, 2π), with radial extent R ∈ [Ri, Ro], on which the radial boundary conditions are "reflecting" conditions 2 . These boundary conditions were chosen to confine the eccentric disc for a well-defined study. While these are not appropriate for all applications in which eccentric discs are thought to arise, they may reasonably approximate the boundaries of a disc that has been tidally truncated (e.g. Papaloizou 2005b). For example, our model might be relevant to narrow rings that may be produced by disc-planet interactions that have formed multiple gaps. A reflecting edge may also be appropriate to describe the inner edge of the disc, if there is a sharp drop in the density there, or where discs match onto the surfaces of central objects e.g. white dwarfs in Cataclysmic Variable systems. However, the aim of this paper is not to focus on any particular application, but to explore the fundamental dynamics of eccentric discs that might have more general applicability. We choose units such that Ω0 = Ri = Σ0 = 1, and vary the parameters cs (to mimic a disc with aspect ratio H/R ∈ [0.025, 0.1]), σ and Ro, in addition to the eccentricity of the disc, which we will now describe.

Eccentric mode
We initialise an eccentric disc using nonlinear secular theory (Appendix A). To do this we modify the velocity components and surface density of the background disc so that they describe an eccentric mode, to make the streamlines elliptical. The eccentric mode is a global slowly precessing density wave that satisfies the boundary conditions.
We define the complex eccentricity E(λ) = e(λ)e iω(λ) , where e and ω are the eccentricity and longitude of pericentre. In secular theory, orbits are labelled using their semilatus rectum λ, related to cylindrical radius by .
In linear secular theory for a 2D isothermal disc, where E ≡ ∂ λ E, and in this case R and λ are equivalent (this can be transformed into the Schrödinger equation). Eq. A10 is the equivalent equation in the nonlinear secular theory of Appendix A (in which care must be taken to take into account the difference between R and λ, and the resulting equation is a type of nonlinear Schrödinger equation). We seek solutions that precess at the rate ωp, such that E ∝ e iωpt . Together with the boundary conditions appropriate for fixed circular boundaries, we obtain an eigenvalue problem for the eigenvalue ωp and eigenfunction E. We solve this problem using a shooting method (Press et al. 1992), which works equally well for the solution of Eq. A10 (which is a nonlinear equation for E and E ). This method requires us to provide an initial guess for ωp and E (λi) to select the appropriate solution. To do this, we use a combination of trial and error and the results of an independently coded solution to the linear problem using a Chebyshev collocation method. In each case, we select the fundamental eccentric mode, which has a single maximum in the eccentricity. This is the slowest precessing mode with the longest radial wavelength, which is chosen because it can be simulated more accurately than any mode of shorter wavelength. For the linear theory, we normalise the amplitude of the resulting mode so that its maximum eccentricity is A. In the nonlinear secular theory, its amplitude is determined by E (λi), which we choose so that the maximum eccentricity is A (to within an accuracy of approximately 10 −4 ). A nonlinear eccentric mode cannot be exactly represented as a single Fourier mode with an azimuthal wavenumber m = 1. However, we can obtain uR, u φ , Σ for the eccentric Keplerian disc as follows. First, we obtain the velocity components of the eccentric orbital motion at each grid point (R, φ) using E(λ): after we have converted points in λ, φ to R, φ (Murray & Dermott 1999). We obtain the surface density by considering mass conservation , by noting that M/P = JΩΣ (see Appendix A for definitions). We assume that the mass distribution in the disc is a power law in λ of the form M(λ) = 2πλΣ b (λ), so that which is the natural extension of our cylindrical disc model. We use Eqs. 7-9 as our initial conditions in PLUTO, noting that there are errors O(c 2 s ) because we have made the secular (thin-disc) approximation (Ogilvie 2001). The eccentric disc is exactly represented in nonlinear secular theory for any eccentricity and eccentricity gradient. Note that we neglect the radial pressure gradient which would slightly modify the angular velocity profile when σ = 0.
We use a second-order Runge-Kutta time-stepping algorithm, with linear interpolation within grid-cells. We adopted the dimensionally-unsplit HLLC solver (which was chosen because preliminary investigation found it to be more robust than the Roe solver, and also much less diffusive than the HLL solver). The eccentric mode is input into PLUTO after appropriate interpolation to the grid used in the code.

PREDICTIONS AND VALIDATION OF NONLINEAR SECULAR THEORY
In this section, we will illustrate the predictions of secular theory and compare them with the results of hydrodynamical simulations. We will focus only on aspects that are directly relevant for such a comparison, deferring a more detailed discussion of the simulations to § 4.

Modification of pressure forces in nonlinear theory
Linear secular theory (Eq. 5, together with boundary conditions) predicts the structure and precession rate of eccentric modes as a function of the disc properties. For the case of an untwisted disc with Ro = λo = 2 and σ = 0, we have plotted the fundamental eccentric mode E(λ) in the top panel of Fig. 1. The amplitude is arbitrary, so this represents the shape of the disc for any amplitude according to linear theory. However, neighbouring orbits (for an untwisted disc) intersect if which is observed to occur for this linear mode at R = Ri if A 0.143 (for Ro = 2; similarly this occurs when A 0.18 for R0 = 3, and A 0.24 for R0 = 10). This suggests that secular theory will no longer apply to these discs, since orbital intersections will rapidly lead to eccentricity damping via shocks. However, this prediction is based on linear theory.
Nonlinear secular theory (Eq. A10, together with the boundary conditions) predicts the shape of the eccentric mode to depend on its amplitude A, as we illustrate in the middle panel of Fig. 1 for several amplitudes as indicated in the caption. Inspection of the functional form of the terms on the right hand side of Eq. A10 informs us why this is the case: these terms increasingly differ from those in linear theory as the orbits approach an intersection, i.e. the shape depends on the amplitude because pressure forces act to minimise intersections. This also causes the precession of the mode to differ from the predictions of linear theory, as we will demonstrate in § 3.2.
With our adopted boundary conditions, nonlinear theory predicts a maximum attainable amplitude for the eccentricity, corresponding to the upper curve in the middle panel of Fig. 1, which occurs when A ≈ 0.332. We have plotted the orbital geometry for the eccentric disc with A ≈ 0.332, with orbits spaced equally in λ, in the bottom panel of Fig. 8. This shows clearly the compression and near-intersection of orbits near the inner and outer boundaries for this mode, in addition to the fact that it can no longer be described purely as an m = 1 mode. The maximum value can be explained simply by considering the (untwisted) eccentric mode which is marginally intersecting, so that e − λe = ±1, that also satisfies the boundary conditions that e(λi) = e(λo) = 0. The solution is the piecewise linear profile where The maximum eccentricity of this mode is which explains our observation of a maximum eccentricity for the fundamental eccentric mode in a domain with a given size (in numerical calculations we observe emax(λo = 2) ≈ 0.332, emax(λo = 3) ≈ 0.496 and emax(λo = 5) ≈ 0.66).
This solution is what we must obtain from geometrical considerations simply because we force the eccentricity to vanish at the boundaries. However, this result depends on the boundary conditions, since this limit would no longer apply if e (λo) = 0, for example, and we would not obtain a maximum eccentricity smaller than one.

Precession rates: validation of theory
According to linear theory, the eccentric mode precesses slowly in a retrograde sense due to gas pressure. To see that this must be the case, we can multiply Eq. 5 by E * , seek solutions E ∝ e iωpt and integrate over λ to obtain: with our chosen boundary conditions (Goodchild & Ogilvie 2006). This quantity is a negative real number if ∂ λ Σ 0, i.e. σ 0, indicating retrograde precession 3 . correspond with the time interval between output files used for this analysis. P pred p is the prediction from nonlinear secular theory (Eq. A10) and P NS,Lin,pred p is the (non-secular) prediction from solution of the two-dimensional linearised isothermal hydrodynamic equations in Appendix B (we have entered the latter in the smallest A entry in the table for each case, but it should be remembered that these results are based on a linear calculation). Simulations were run for at least one full precession period.
The precession rate depends on σ, Ro and cs, but is independent of the amplitude of the eccentric mode in linear theory. However, we have shown that nonlinear secular theory predicts the structure of the eccentric mode to be modified if the orbits are close to intersecting, so we might expect ωp to depend on its amplitude also. We demonstrate that this is indeed the case in the black solid lines in Fig. 2, where we have plotted the precession period Pp = 2π ωp against the amplitude of the eccentric mode for several calculations with different Ro ∈ [2, 3, 5] and various A with σ = 0 and assuming a thin disc with cs = 0.05. Gas pressure causes the mode to precess faster, in a retrograde sense, as the eccentricity and eccentricity gradients are increased. We obtain the linear prediction in the limit A → 0, with a departure that is initially ∝ A 2 , before steepening for large A. The precession period continues to shorten until A → emax, when Pp → 0, so that the secular approximation is no longer valid. Nonlinear secular theory therefore predicts its own breakdown for large amplitudes, which occurs when the orbits in the disc come close to intersecting.
In Fig. 2 we have also plotted the initial precession period computed from hydrodynamical simulations (symbols according to the legend) that were run for at least one full precession period. A table of simulations used for this comparison, and those in the rest of this section, is given in Table 1. We calculate the precession period by computing the m = 1 component of the radial velocity (this quantity is used to represent the power in both m = ±1), from which the mean longitude of pericentre can be computed from which is approximately valid for the eccentric mode even when A is not small. Since this mode precesses retrogradely, ω decreases cyclically from π to −π, and we calculate the initial precession period by eye by checking at what time ω ≈ π after t = 0. The data used for this analysis is output at each time unit in the simulation, giving errors of ±1 in the determination of P obs , which we further verified by visual inspection of uR(R, φ, t) at these time snapshots. Simulations with A 0.15 for Ro = 2 were not plotted, since the eccentric mode underwent non-negligible damping during a single precession period, as we will discuss further in § 4. Fig. 2 shows that our results are in good agreement with the nonlinear secular theory for each Ro considered.
To further validate the predictions of secular theory, we have listed the observed and predicted values of Pp as various disc (and simulation) parameters are varied in Table 1. The variation of Pp with σ and cs is reasonably captured in each case. However, the departure from the predictions of secular theory increases as we increase cs: for . Precession period Pp versus peak eccentric mode amplitude A showing a comparison between nonlinear secular theory (black solid lines) and the results of hydrodynamical simulations in a domain Ro = 2, 3 and 5 (red crosses, light blue circles and blue squares, respectively), cs = 0.05 and σ = 0. This demonstrates that the eccentric mode precesses faster as its amplitude is increased. We observe good agreement between simulations and nonlinear secular theory as A and Ro are varied (with a departure of approximately 3%). Note that secular theory predicts its own breakdown for large amplitudes, as Pp → 0. cs = [0.025, 0.05, 0.075, 0.1], the fractional error correspondingly takes the approximate values [1, 3.2, 7.5, 10]%, indicating a roughly c 2 s dependence for this quantity, as expected from secular theory (Ogilvie 2001;. To confirm that this departure indeed results from the partial inaccuracy of the secular approximation, we have also computed the precession frequency of an m = 1 linear eccentric mode by solving the two-dimensional linearised isothermal (non-secular) hydrodynamic equations as outlined in Appendix B, for the cases listed in the table. This accurately matches the observed precession rates of our smallest A simulations in each case. For thin discs with cs 0.1, the precession period predicted by secular theory matches the observed values to within a few percent, typically tending to predict slightly slower precession.
We have also computed the precession rate for a set of calculations with A = 0.05 as the resolution is varied, as indicated in Table 1. This indicates convergence is attained for the precession period even at low resolutions such as 100×128, but a departure appears for even lower resolutions. This suggests that simulations with typically-adopted resolutions would be expected to capture the precession rates of global eccentric modes, but may not accurately capture the precession of shorter wavelength modes.

Summary
We have shown that the retrograde precession of eccentric discs as a result of their gas pressure is enhanced for sufficiently large eccentricities and eccentricity gradients. This dependence of the precession rate on the eccentricity and its gradient is a new result 4 (though there was some numerical evidence of this in Papaloizou 2005b), and occurs because pressure forces are enhanced when the orbits come close to intersecting. While the numerical results would change if we were to adopt different boundary conditions, this general behaviour is likely to be robust since it arises from the functional form of the stress integrals (see Appendix A).
We have thoroughly validated the predictions of nonlinear secular theory regarding the precession rates of eccentric discs as various parameters are varied, in the regime where ωp Ω. We find that the departure from secular theory scales as O(c 2 s ), typically being a few percent or smaller for thin discs if ωp Ω (which we have shown to break down at large amplitudes). We now turn to a more detailed analysis of our simulation results, focussing on the long-term evolution of the eccentricity in two dimensions.

Background instability: nonlinear evolution of the "Papaloizou-Pringle" instability
Our basic setup has rigid walls to confine an eccentric mode and permit a well-defined study. However, it is well known that a circular supersonic rotating shear flow in a container with one or more rigid boundaries is unstable to the development of non-axisymmetric instabilities (Papaloizou & Pringle 1984Goldreich & Narayan 1985;Goldreich et al. 1986;Papaloizou & Pringle 1987;Kato 1987;Narayan et al. 1987). These instabilities are driven either by Kelvin-Helmholtz-type mechanisms (which require a potential vorticity minimum, which is not the case in our problem), or by over-reflection of spiral density waves from the corotation region. While most of the original work on this problem focused on thick discs (referred to as "accretion tori") where this instability excites low azimuthal wavenumber (m = O(1)) modes on a local orbital timescale, slower growing instabilities with high azimuthal wavenumbers (m 1) have also been found to occur in thin Keplerian discs (Papaloizou & Pringle 1987;Narayan et al. 1987;Hanawa 1987). The nonlinear evolution of these instabilities has been studied by Godon (1998) (and earlier for thick "accretion tori" by Hawley 1987), where they were observed to drive subsonic wave activity in hydrodynamic discs.
The basic mechanism of instability in a thin nearly Keplerian disc can be explained as follows. An outgoing spiral density wave with frequency ωm and azimuthal wavenumber m = 0 with a radial angular momentum flux I that approaches its corotation region (the region where the wave frequency satisfies ωm ≈ mΩ) is primarily reflected from its inner Lindblad resonance (where ωm −mΩ = −κ, and κ ≈ Ω is the epicyclic frequency), with the angular momentum flux in the reflected wave being R. However, this wave is also partially transmitted through the corotation region (in which it is evanescent), requiring an outgoing wave to be launched at the outer Lindblad resonance (where ωm − mΩ = κ) with transmitted angular momentum flux T = I − R. In a Keplerian disc, waves propagating inside corotation (ωm mΩ) carry a negative angular momentum flux (so that I < 0), whereas those outside corotation (ωm mΩ) carry a positive flux (so that T > 0). This means that R = I −T < I, so that the reflected wave is amplified. If there is an inner impermeable boundary, the reflected wave will be completely reflected once more so that it will re-enter the corotation region, enabling further amplification at the expense of the Keplerian flow. This mechanism of instability, which is due to the "over-reflection" of waves from the corotation region, also occurs if there is a rigid outer boundary, but does not occur if both boundaries perfectly transmit wave energy, and it can be eliminated with sufficient viscosity.
While our primary aim is to study the evolution of eccentric discs, this instability (which is an artifact of our model) will drive subsonic wave activity that will interact with the eccentric mode. Because of this, we briefly describe the properties of this instability for a circular disc in this section, before discussing the evolution of eccentric discs in the next. In Figs. 3 and 4, we present calculations that show the energy in the lowest ten azimuthal wavenumbers wherê which is defined so that it represents the power in the solution with a given value of |m|, so that only m > 0 need to be considered for this quantity. We compare three different numerical resolutions NR × N φ : 200 × 256 (labelled as L), 400 × 512 (labelled as M) and 800 × 1024 (labelled as H). Each simulation has Ro = 2, cs = 0.05, and is started with random white noise perturbations to each component of the velocity on the grid-scale 5 with amplitude 10 −5 cs.
Since this instability preferentially excites high azimuthal wavenumbers and there is no explicit viscosity in our simulations, we expect some dependence on resolution, though we will show that the results for M and H appear to have converged. Fig. 3 illustrates that simulation L exhibits instability with the excitation of m ∼ 5 propagating waves by t ∼ 2000, but this initially saturates at low amplitude. At later times, m ∼ 3 standing waves develop throughout the domain. We illustrate the velocity field during these stages in Fig. 5 for this simulation, in the early linear growth phase (top) and during the later phase (bottom). The RMS radial velocity driven by the instability, normalised by the sound speed, is plotted in Fig. 6 for all three simulations, along with the efficiency of angular momentum transport 5 No attempt is made to ensure that these are identical for each resolution adopted.   During simulation L , the maximum RMS velocities attained are very subsonic, with u ∼ 10 −3 −10 −2 cs, and the angular momentum transport is weak, with α 10 −5 . Simulations M and H reach larger turbulent velocities of u ∼ 0.03cs and somewhat more efficient angular momentum transport with α ∼ 10 −4 (Fig. 6) -note that α = 10 −4 implies tvisc = R 2 i αcsH = 1 αc 2 s ≈ 4 × 10 6 . The initial instability in these simulations preferentially excites m ∼ 9 components, with m ∼ 6 dominating at later times, which differs from simulation L. However, the instability appears to be well captured with the resolution adopted for simulation M, since simulation H does not differ significantly. In addition, results do not appear to be strongly dependent on the domain size, with simulations in a domain with Ro = 3 (where resolutions L and M correspond with 300×256 and 600×512 grid-points, respectively) giving similar spectra at late times (compare Fig. 3 with Fig. 4) and broadly comparable turbulent velocities and angular momentum transport (Fig. 6) -note that the initial turbulent stages are different, with m ∼ 6 modes excited initially that remain dominant at late times. (We have also checked that the outcome of the in- propagating spiral waves are produced at inner boundary) and at a later stage (bottom; global radial standing waves are shown) of the instability in our lowest resolution simulation, illustrating the non-axisymmetric waves generated by the instability. stability does not depend significantly on σ, at least where σ 1.5.) The instability excites spiral density waves that saturate with subsonic velocities and lead to weak but nonnegligible angular momentum transport. This instability provides background wave activity which complicates the analysis of eccentric modes in the next section. Depending on the particular waves excited by this instability, their nonlinear interaction with the predominantly m = 1 eccentric mode could either drain energy from, or transfer energy into, this component. Hence, we would not expect secular theory to be valid in the presence of a strong non-axisymmetric component, given that it neglects wave-wave couplings.
Given that this instability is minimised for the resolution adopted for simulation L, we primarily focus on simulations using this resolution in the next section, where we turn to analyse the long-term evolution of eccentric discs.

Long-term nonlinear evolution of free eccentric modes
As in § 3.2, we initialise the flow with an eccentric mode with peak eccentricity amplitude A and study its nonlinear evolution, focussing on its behaviour over many (50 − 70) precession periods. We begin by analysing the temporal evolution of the eccentric mode energy for various A in Fig. 7. We plot E1 (as defined in Eq. 17), which approximately represents the energy in the eccentric mode even for moderate A, when m = 1 components are also present (see Fig. 8). The validity of sec- ular theory is indicated by the lack of evolution in E1 during the initial stages. The eccentric mode persists throughout our simulations, but experiences two different kinds of damping as the simulation progresses. Firstly, numerical diffusion acts to damp the eccentric mode, which is particularly pronounced for cases with higher A. Since the streamlines are then highly concentrated (see e.g. the bottom panel of Fig. 1), and there are strong eccentricity gradients, these modes experience appreciable damping by numerical diffusion on the grid-scale. However, this is only the most important damping mechanism during the initial stages when A 0.125. During later stages (i.e. after t 1000), the dominant damping mechanism is the interaction of the m = 1 eccentric mode with other non-axisymmetric waves that are driven by the instability discussed in § 4.1. These waves are excited in the absence of an eccentric mode, but interact with the eccentric mode through nonlinear interactions to generate additional non-axisymmetric components. In most cases this interaction acts to damp the eccentric mode, but it can also transfer energy into the eccentric mode in some cases e.g. see A = 0.005 when t ∼ 2500. In Fig. 8, we have plotted the temporal evolution of Em for m ∈ [1, 10] to illustrate the growth of m = 1 components of the solution. The top left panel repeats our results when A = 0 (Fig. 3) to provide a point of reference. The remaining panels show the same quantity for several different values of A. When A = 0.005, background instability occurs at approximately the same time as the A = 0 case, preferentially exciting similar modes (together with those with azimuthal wavenumbers that differ by 1), and the m = 1 components of the solution saturate with similar energies (∼ 10 −7 ). The m = 1 component does not undergo significant damping until t ∼ 2500, after the background instability has set in. This supports our interpretation that this damping of the eccentric mode is due to its interaction with non-axisymmetric waves driven by the background instability, and not due to a separate instability of the eccentric mode itself.
In Fig. 9, we have plotted the same quantity but for a set of simulations that have double the radial and azimuthal resolution (starting with simulation M of § 4.1). The eccentric mode evolves similarly in these cases as in the lower resolution simulations in Fig. 8. However, the background instability is excited earlier in the simulations with small A, which leads to slightly different long-term quantitative evolution for the eccentric mode amplitude (in fact there is slightly more efficient damping in some of these simulations compared with those with the lower resolution, demonstrating that the amplitude decay is not primarily due to numerical diffusion).
As A is increased, this instability occurs sooner in the simulation, presumably because there is more energy in m = 1 components of the solution at t = 0. However, with the exception of the m = 2 component (which arises pri- The instability that occurs in the absence of an eccentric mode also occurs when A = 0, which excites nonaxisymmetric waves that subsequently interact with the eccentric mode, generally leading to damping. However, the eccentric mode itself does not appear to be subject to separate instabilities in two dimensions. marily to describe the eccentric mode itself when A 0.05), other non-axisymmetric components typically saturate with similar power to the A = 0 case. In all simulations, there is an amplitude-and time-dependent damping (or transient growth, as briefly observed at t ≈ 2000 in the A = 0.005 simulation in Fig. 8) of the eccentric mode due to its interaction with waves driven by the background instability. However, our interpretation, which is guided by Figs. 8 and  9, is that the eccentric mode does not appear to be subject to additional instabilities due to the presence of a nonzero eccentricity (at least in two dimensions), only to the finite-A modification of the instability that occurs when A = 0 due to our adoption of rigid walls.  Figure 10. Magnitude of the radial velocity in the m = 1 component of the solution (|û 1 |) as a function of R at the beginning of the simulation (black solid lines; taken from nonlinear secular theory) and at t = 20000 (red dashed lines; representing the end of our simulations). This illustrates the damping of the eccentric mode, and also that its shape remains similar to that predicted by secular theory.  Figure 11. Longitude of pericentre ω as a function of R for simulations at times t = [0,30,60,90,120,150,180,210,230,260,290]. This illustrates that the eccentric disc remains approximately untwisted during one precession period, though slight twists can develop.
In Fig. 10, we plot |û1(R)| at t = 0 (solid black lines) and t = 20000 (red dashed lines; marking the end of each simulation), which shows the damping of the eccentric mode. In addition, this demonstrates that the eccentric mode remains coherent, persisting throughout the duration of these simulations in a form that is similar to the initial conditions. The similarity in the mode shape at these two times supports the validity of secular theory in describing the shapes of these modes. Our initial disc model is untwisted, but we do not constrain it to remain so. In Fig. 11, we show ω(R) for A = 0.01 and A = 0.1, for times t ∈ [0, 300] (slightly more than Pp) spaced in intervals of 30 time units. At t = 0, the disc is untwisted with ω = π, but slight twists develop during the course of the simulation, preferentially near to the outer boundary. However, the disc remains approximately untwisted as it evolves throughout these simulations with a net twist that is small even when A = 0.1 (smaller than 0.5 rads). It is possible that numerical damping, which is likely to act in a similar way to a shear viscosity (Ogilvie 2001), could be responsible for this twist in the outer regions (where the grid cells are largest), thereby preventing the disc from remaining entirely untwisted (also, if the initial conditions are not exact nonlinear solutions, some twist would be expected to develop). Alternatively, since e → 0 as R → 2, the phase of the complex eccentricity becomes undefined at this location, so we might expect the twist at this location to be arbitrary.
Finally, we analyse the decay rate of the eccentric mode, which we plot in Fig. 12 for several simulations, showing the effects of varying the resolution. This picture is not clear-cut because the background instability is stronger for the higher resolution case. This illustrates that eccentric modes in two dimensions persist for a very long time, even in the presence of amplitude-dependent wave-wave interactions and numerical damping.

Summary
In this section we have presented the results of nonlinear hydrodynamical simulations designed to study the longterm evolution of eccentric discs in two dimensions. Eccen- tric discs remain coherent and do not appear to be subject to instabilities as a result of their free eccentricity in two dimensions, and would presumably persist forever if there was no background instability. They are, however, damped (typically) by their nonlinear interaction with nonaxisymmetric spiral density waves -these waves are driven by a background ("Papaloizou-Pringle") instability in our setup with rigid walls 6 . In real discs there are various mechanisms that could produce such an incoherent ensemble of non-axisymmetric waves e.g. additional sources of turbulence such as magneto-rotational instability (e.g. Heinemann & Papaloizou 2009a,b), gravitational instability (e.g. Papaloizou & Savonije 1991;Laughlin & Rozyczka 1996), convection (e.g. Mamatsashvili & Rice 2011, or possibly the tidal interaction between multiple proto-planets and the disc. Previous simulations of eccentric modes by Papaloizou (2005b) used linear secular theory to construct the initial conditions. This enhances the damping of the eccentric mode due to the generation of shocks at the inner boundary of the domain, since we have observed the linear mode to cause orbital intersections if A is sufficiently large. On the other hand, the modes that we have input using nonlinear secular theory do not experience such strong shock-induced damping, as indicated in Fig. 7. These modes do experience amplitude-dependent damping by numerical diffusion, which can be particularly strong for large A cases, where the orbits are closer to intersecting. But in our simulations the most important damping mechanism is the interaction of the eccentric mode with other non-axisymmetric waves.

CONCLUSIONS
In this paper we have studied fundamental aspects of eccentric gas discs that should have general applicability. We have derived analytically a nonlinear secular theory for isothermal two-dimensional untwisted eccentric Keplerian discs (valid for arbitrary eccentricities and eccentricity gradients for which neighbouring orbits do not intersect; Appendix A), and verified its predictions with idealised hydrodynamical simulations using the PLUTO code. Linear secular theory is found to accurately describe the structures and precession rates of eccentric discs with small eccentricities (and eccentricity gradients), which precess in a retrograde sense due to gas pressure. Discs with larger eccentricities (and eccentricity gradients) are observed to precess at a faster rate, which we have explained as a modification of the pressure forces, and resulting disc structure to prevent orbital intersections, that occurs when the orbits in a disc nearly intersect.
The nonlinear modification of the pressure forces might be particularly important for the highly eccentric discs produced in tidal disruption events (Guillochon et al. 2011;Liu et al. 2013;Coughlin & Nixon 2015), or for narrow eccentric gas rings (cf. Borderies et al. 1983, which was applied to planetary rings). Another potential application is to the period excess in superhump binary systems, which is thought to be explained due to the precession of an eccentric gaseous disc (Whitehurst 1988;Lubow 1991;Murray 2000;Goodchild & Ogilvie 2006;Smith et al. 2007). However, there is observational evidence of temporal variability of the period excess (Mason et al. 2008;Nakata et al. 2014), which may in part be caused by evolution of the pressure forces due to the varying eccentricity and eccentricity gradients in the disc. It would be worthwhile to apply nonlinear secular theory to this problem (e.g. extending Goodchild & Ogilvie 2006) in order to explore this possibility further.
Eccentric discs do not appear to exhibit hydrodynamic instability as a result of their free eccentricity in two dimensions, and we would expect them to essentially persist forever in our simulations except for their interaction with non-axisymmetric spiral density waves excited by a background (Papaloizou-Pringle) instability (which arises in our setup with rigid walls), in addition to numerical viscosity. Presumably the eccentricity would be damped on the disc (turbulent) viscous timescale, but our simulations are not run for this duration (the viscous timescale is very long because the background instability transports angular momentum only weakly). Eccentricity can also be affected by nonadiabatic thermal processes, not considered in this paper. Indeed, (Statler 2001) has argued that a uniform eccentricity is preferred because the surface density of the disc is then constant around each orbit. However, even in this case the fluid undergoes periodic compression because of the variation of vertical gravity around the elliptical orbit (which is not captured in two dimensions).
Our observation that gas discs with free eccentricities can remain eccentric for thousands of orbital periods (Fig. 7) is potentially important regarding the interpretation of observed large-scale asymmetries in gas discs. In particular, the presence of such asymmetries is often used to infer the presence of perturbing planets (e.g. Regály et al. 2014;Hashimoto et al. 2015;Pinilla et al. 2015). However, since we have shown that eccentric modes can be long-lived features, there are alternative mechanisms in addition to perturbing planets that can induce free eccentricity, e.g. gravitational instabilities in the gas disc during the pre-main sequence phase (Adams et al. 1989).
However, it should be noted that any conclusions regarding the longevity of disc eccentricities are based on twodimensional simulations with only a weak background instability. In three dimensions eccentric discs are subject to a parametric instability that excites small-scale inertial waves (Papaloizou 2005a;. This is expected to lead to enhanced damping of the disc eccentricity and eccentricity gradients (Papaloizou 2005b), but further work adopting a realistic vertical disc structure is required to determine its nonlinear outcome. The simulations presented in this paper provide a starting point to allow us to undertake such calculations.

ACKNOWLEDGEMENTS
AJB is supported by the Leverhulme Trust and Isaac Newton Trust through the award of an Early Career Fellowship. The early stages of this research were supported by STFC through grants ST/J001570/1 and ST/L000636/1. We would like to thank the referee for constructive comments that led us to strengthen our conclusions.

APPENDIX A: NONLINEAR SECULAR THEORY OF ECCENTRIC DISCS
In this section we derive analytically a nonlinear secular theory for eccentric discs in two dimensions utilising the local model of , for the case of an untwisted isothermal disc in which λω = 0, i.e. the orbits are aligned (for which nonlinear pressure effects are likely to be minimised). While the formalism can be extended to the case of a twisted disc including viscosity, and to discs with different thermodynamic behaviour (Ogilvie 2001), we do not believe it to be possible to derive the corresponding nonlinear theory analytically in terms of algebraic functions (particularly the stress integrals) except for this special case. We denote cos θ by c and sin θ by s to simplify the expressions below, where the θ = φ − ω(λ), is the true anomaly. At each λ, the reference orbit is Keplerian, so and the angular velocity of the orbital motion is The orbital period for a fluid element with this orbit is and its specific angular momentum is = √ GM λ. The evolution of the mass, angular momentum and eccentricity of the eccentric disc are governed by the following stress-integrals, which we evaluate for a general untwisted eccentric disc with non-intersecting orbits: JR 2 T λφ e iφ dφ dz = ic 2 s e e 2 − 1 λM (λe ) 2 e 2 − eλe − 1 + (e 2 − 1) ((e − λe ) 2 − 1) , where T ij = −pg ij and J = RR λ . In each case the integrals are over the full extent of the disc in φ and z. We have defined R λ ≡ ∂ λ R and R φ ≡ ∂ φ R, so that the relevant components of the metric tensor are The one-dimensional mass density is The mass and angular momentum do not evolve since the first integral above vanishes, which results from us neglecting shear viscosity. The complex eccentricity evolves according to M∂tE = 2e iφ ∂ λ (JR 2 T λφ ) − iλ e iφ ∂ λ (JR λ T λλ ) −i e iφ JR 2 T φφ − JR 2 λ e iφ T λφ dφ dz, (A10) which is a type of nonlinear Schrödinger differential equation (that is second order in λ). For general initial conditions, the time-evolution will produce a twist in the disc, but it is possible to seek modal solutions that remain untwisted. In this case, inputting the stress integrals above leads to an equation for the evolution of the complex eccentricity for any eccentricity and eccentricity gradient for an untwisted disc, as long as the orbits do not intersect. The final form of this equation after substituting the stress integrals is too complicated to be worth writing down in closed form here. It can, however, be computed in Mathematica and exported to text format for input in our Matlab scripts that are used to solve Eq. A10.
We seek modes E ∝ e iωpt that precess slowly at the rate ωp, so that Eq. A10 is converted to an eigenvalue problem of the form ωpe = f (λ, e, e , e ), subject to appropriate boundary conditions, where f is in general a nonlinear function of its arguments. The behaviour of the stress integrals for large e, e is straightforward to understand. As the orbits become closer to intersecting ((e − λe ) 2 → 1) or e → 1 or (λe ) 2 → 1, strong pressure forces arise, which modify the eccentricity profile so that the orbits do not intersect. An interesting result of these forces is that they modify the global precession rate of the mode throughout the disc so that it precesses more rapidly in a retrograde sense as the orbits become closer to intersecting, even if this near intersection only occurs in a small region of the disc. We note that the formalism of Borderies et al. (1983) for narrow eccentric rings (with small eccentricities but arbitrary eccentricity gradients) applied to an isothermal gas also predicts a precession rate (their Eq.13) that depends on the eccentricity gradient, and which diverges in the limit of intersecting orbits. This is consistent with our results. The 2D isothermal model is a simplification of a realistic 3D eccentric disc. The isothermal case (with adiabatic index of 1) is the most compressible 2D model that we can consider, where the pressure increase for a given streamline convergence is minimised (over less compressible models with adiabatic indices larger than 1). However, in 3D the disc can also expand vertically, which may somewhat alleviate this pressure increase. However the qualitative behaviour of the 2D isothermal model is likely to carry over to more realistic models.

APPENDIX B: LINEAR NON-SECULAR THEORY IN 2D
As listed in Table 1, we solve the eigenvalue problem for the full (non-secular) isothermal hydrodynamical equations in two dimensions to compute the precession frequency of the fundamental m = 1 linear eccentric mode, to allow comparison with simulations and secular theory. We assume linear perturbations u R , u φ , Σ ∝ e i(φ−ωpt) , which satisfy where we have taken into account the radial pressure gradient in full. This is solved using a Chebyshev collocation method (adopting 101 points in radius, which we find to be more than sufficient to compute the fundamental mode), which converts the problem to a standard eigenvalue problem that can be solved (for which we use a QZ method) to compute the precession frequency ωp.