Nonlinear evolution of the magnetorotational instability in eccentric disks

The magnetorotational instability (MRI) has been extensively studied in circular magnetized disks, and its ability to drive accretion has been demonstrated in a multitude of scenarios. There are reasons to expect eccentric magnetized disks to also exist, but the behavior of the MRI in these disks remains largely uncharted territory. Here we present the first simulations that follow the nonlinear development of the MRI in eccentric disks. We find that the MRI in eccentric disks resembles circular disks in two ways, in the overall level of saturation and in the dependence of the detailed saturated state on magnetic topology. However, in contrast with circular disks, the Maxwell stress in eccentric disks can be negative in some disk sectors, even though the integrated stress is always positive. The angular momentum flux raises the eccentricity of the inner parts of the disk and diminishes the same of the outer parts. Because material accreting onto a black hole from an eccentric orbit possesses more energy than material tracing the innermost stable circular orbit, the radiative efficiency of eccentric disks may be significantly lower than circular disks. This may resolve the"inverse energy problem"seen in many tidal disruption events.


INTRODUCTION
Eccentric gaseous disks arise in a surprisingly wide variety of astrophysical contexts.A number of mechanisms can explain the existence of eccentric disks, the most commonly invoked one being external perturbation.In eccentric binaries, secular gravitational interaction endows forced and free eccentricities upon circumbinary and circumobject disks (e.g., Murray & Dermott 1999); in circular binaries, tidal forces couple to circumobject disks through the 3 : 1 mean motion resonance and allow free eccentricity to grow exponentially (Lubow 1991).Another possibility is that disks become more eccentric over time.Viscous overstability (Kato 1978), which amplifies small-scale eccentric perturbations in isolated disks (e.g., Lyubarskĳ et al. 1994;Ogilvie 2001), is often cited in this connection.A third option is for disks to be born eccentric.Outgassing from planetesimals can create eccentric disks (Trevascus et al. 2021), and so can the tidal disruption of stars (Shiokawa et al. 2015;Piran et al. 2015;Svirski et al. 2017) and molecular clouds (e.g., Bonnell & Rice 2008) by supermassive black holes.On the phenomenological side, eccentric disks are sometimes invoked to explain asymmetric lines in white dwarfs (e.g., Gänsicke et al. 2006), as well as asymmetric broad emission lines in active galactic nuclei (e.g., Eracleous et al. 1995;Tucker et al. 2021) and tidal disruption events (TDEs) (e.g., Guillochon et al. 2014;Liu et al. 2017).
Because ideal magnetohydrodynamics (MHD) is supported even by low levels of ionization (Blaes & Balbus 1994;Gammie 1996), we expect magnetic fields to play a role in many of the eccentric disks enumerated above.The presence of magnetic fields changes the way disks evolve because of the magnetorotational instability (MRI) (Balbus & Hawley 1991;Hawley & Balbus 1991).Simply put, in a disk whose inner parts rotate faster than the outer parts, differential rotation can latch onto horizontal bits of the magnetic field, stretch them out, and amplify them.The gas connected to one of these bits on the inside is pulled back by magnetic tension, loses angular momentum, migrates inward, and picks up orbital speed.In the meantime, the gas on the outside is dragged forward, gains angular momentum, drifts outward, and slows down.The rising velocity difference across the horizontal magnetic field in turn enhances its stretching, precipitating an instability.
Analytic calculations for circular disks show that a perturbation can grow by orders of magnitude per orbit in the linear stage (Balbus & Hawley 1991), making the MRI among the most vigorous MHD instabilities.The initially exponential amplification eventually enters the nonlinear stage and breaks down into MHD turbulence.Orbital shear enforces a correlation between the radial and azimuthal components of the turbulent velocity, and between the same components of the magnetic field.Turbulent stresses transport angular momentum outward; gas robbed of angular momentum sinks to smaller radii, and the disk accretes.
The saturation process is amenable only to numerical investigation.Previous simulations of circular disks, numbering in the hundreds, are divided into shearing-box simulations, which consider a small neighborhood of the disk as representative of the whole (e.g., Hawley & Balbus 1992;Hawley et al. 1995Hawley et al. , 1996;;Brandenburg et al. 1995;Stone et al. 1996), and global simulations, which take in the entire disk (e.g., Hawley & Balbus 1991;Armitage 1998;Matsumoto 1999;Hawley 2000;De Villiers & Hawley 2003).This large body of work converged on the consensus that, irrespective of the circumstances simulated, MHD turbulence in circular disks saturates within 10 to 20 orbits, and the stresses at saturation correspond to a Shakura & Sunyaev (1973) alpha parameter between 0.01 and 1.For circular disks around black holes, the alpha parameter may change substantially near their inner edges at the innermost stable circular orbit (ISCO).In disks with gas-dominated pressure, the alpha parameter can increase very rapidly as matter approaches and crosses the ISCO (Noble et al. 2010); alternatively, in super-Eddington radiation-dominated disks, it may exhibit a sharp peak at a radius a short distance outside the ISCO (Jiang et al. 2014).
There is no reason to expect the MRI and the associated MHD stresses to be absent from eccentric disks, though their character and the exact manner in which they approach saturation may differ from circular disks.Very little heed, however, has hitherto been paid to any aspect of the role of magnetic fields in eccentric disks.We were the first to establish analytically that eccentric disks are susceptible to the MRI (Chan et al. 2018).Compared to the circular MRI, the growth rate of the eccentric MRI is smaller at the order-unity level and the range of unstable wavelengths is wider.That work, however, is incomplete because it examined only linear stability.It remains an open question whether the robust growth of MHD stresses in the linear stage would, as the MRI turns nonlinear, translate to saturated stress levels significant enough to affect disk evolution.
As of writing, only a couple of simulations have looked at how eccentricity interacts nonlinearly with the MRI.Dewberry et al. (2020a) set up a circular disk in the potential of Paczyńsky & Wiita (1980) and excited eccentric waves from the outer edge; they saw that the MRI is active at large radii where eccentricities are higher, but suppressed near the ISCO where eccentricities are lower and differential apsidal precession is stronger.Oyang et al. (2021) fed the disk around one member of a binary through Roche-lobe overflow; they found that the  -component of the magnetic stress aids tidal gravity in growing eccentricity but all other stress components oppose it.Neither work achieved eccentricities ≳ 0.1 because doing so requires overcoming two challenges.
The first problem stems from the fact that existing Newtonian MHD codes are capable of handling only Cartesian, cylindrical, and spherical coordinate systems.If a moderately eccentric disk were simulated in one of these coordinate systems, the streamlines would be oblique to the grid, and so would the magnetic fields dragged out by orbital shear.Numerical artifacts reflecting grid symmetry would creep in; at the same time, excessive numerical dissipation would prevent the disk from maintaining its shape over a long time.These drawbacks limited previous simulations of isolated eccentric disks to eccentricities small enough to be implementable as an  = 1 perturbation to the initial velocity (e.g., Papaloizou 2005a), but alternative approaches do exist (Barker & Ogilvie 2016;Dewberry et al. 2020b,a).Here we demonstrate that Newtonian simulations can be performed in arbitrary coordinate systems with generalrelativistic magnetohydrodynamics (GRMHD) codes, provided that the metric is judiciously chosen.Employing a coordinate system molded to the shape of moderately eccentric disks significantly suppresses the numerical errors arising from ordinary cylindrical coordinates.
The second difficulty is with the simulation setup.One may think that the setting of localized perturbations in Chan et al. (2018) lends itself naturally to shearing-box simulations (e.g., Ogilvie & Barker 2014;Wienkers & Ogilvie 2018).Drawing inspiration from circular disks, one may imagine shearing boxes in eccentric disks to have edges running along curves of constant semilatus rectum and constant azimuth.However, the very notion of an eccentric shearing box is suspect.Circular shearing boxes assume that the disk, being homogeneous, is equivalent to a tiling of the shearing box; this justifies periodic azimuthal and shift-periodic radial boundary conditions.The assumption breaks down for eccentric disks.Chan et al. (2018) showed that a perturbation grows differently at different positions along the orbit, depending on the local orbital shear.This means the conditions at the leading edge of an eccentric shearing box are different from the trailing edge, and they also vary along each of the other two edges, so boundary conditions that directly copy one edge to the other would be inappropriate.We can avoid these questions about the eccentric shearing box by performing global simulations instead.It is worth noting that the first simulations of the circular MRI were also global (Hawley & Balbus 1991).
We recount our simulation setup in Section 2. The results from the simulations are presented in Section 3 (see movies) and discussed in Section 4. Our concluding remarks are gathered in Section 5.

METHODS
We outline our simulation strategy in Section 2.1.We continue with the details of the simulation setup in the subsequent subsections and in the Appendix; readers uninterested in the technicalities may skip to the results in Section 3.

Overall strategy
Our goal is to simulate the nonlinear evolution of the eccentric MRI in a purely Newtonian setting.The only reason we turn to a GRMHD code is because numerical issues force us to tailor the coordinate system to the eccentric disk shape, but existing Newtonian codes lack the facility to deal with bespoke coordinate systems.The GRMHD code we use for this purpose is Athena++ (White et al. 2016;Stone et al. 2020).
There can be drawbacks to solving Newtonian problems with a GRMHD code, principally the large truncation error potentially created by the smallness of the typical kinetic and internal energies compared to the rest energy.This error can, however, be mitigated by careful design of the simulation setup, as described in later subsections.With our setup, the truncation error is ≲ 10 −7 for the vast majority of cells.
A major difference of GRMHD codes compared with Newtonian codes is that gravity enters not as explicit momentum and energy source terms, but through the metric.Our choice of the metric in Sections 2.3 and 2.4 ensures that orbits are closed ellipses that do not apsidally precess, allowing our simulations to closely approximate Newtonian behavior.
We simulate a suite of disks.The initial hydrodynamic configuration of the disks follows the common prescription in Section 2.5: Gas is placed within a limited radial range, so that the inner and outer edges of the disk are well-separated from the inner and outer boundaries, respectively, of the simulation domain.All disks are tracked for 15 orbits so turbulence may have enough time to reach saturation.The individual disks comprising the suite are designed with contrast in mind.They are classified along two orthogonal dimensions: circular versus eccentric, unmagnetized versus magnetized.
Comparison between circular and eccentric disks gives us an idea whether the saturation level of the MRI depends on eccentricity.The eccentric disks have a moderate eccentricity of 0.5, so that the character of the MRI specific to eccentric disks can reveal itself without being overwhelmed by hydrodynamic effects.
Comparison between unmagnetized and magnetized disks helps us disentangle MHD effects from hydrodynamic effects.Magnetized disks are seeded with an initial magnetic field as described in Section 2.6.Two magnetic topologies are considered, vertical-and dipolar-field, because topology can influence the saturated MHD turbulence (e.g., Hawley et al. 1995Hawley et al. , 1996;;Sano et al. 2004;Bai & Stone 2013).

Equations
We employ natural units, which means the length and velocity units are the gravitational radius and the speed of light, respectively.The sign convention is (−, +, +, +), Greek indices range over {0, 1, 2, 3}, Latin indices range over {1, 2, 3}, and Einstein summation is implied.
The equations of GRMHD are Here  is the coordinate time,  is the comoving mass density,   is the velocity,  is the determinant of the metric   , and Γ   is the Christoffel symbol of the second kind.From the Hodge dual of the electromagnetic tensor *   we obtain the magnetic field   = *  0 and the projected magnetic field   =   *    .Lastly, the stress-energy tensor is (4) where  and  are the gas pressure and adiabatic index, respectively.

Eccentric coordinate system
We solve the GRMHD equations in an eccentric coordinate system (Ogilvie 2001).It is similar to the cylindrical coordinate system, except that circular coordinate surfaces are replaced by axially aligned elliptical ones, chosen such that their cross sections along the midplane match the initial disk streamlines.Adapting the coordinate system to the disk reduces numerical artifacts and dissipation in our simulations.The eccentricities and orientations of the coordinate surfaces can in principle vary from one elliptical cylinder to the next, but here we specialize to the case in which both are spatially uniform.
Let (, , , ) be cylindrical coordinates.We work with gravity weak enough to be well-described by a quasi-Newtonian potential Φ(, ); the nonzero components of the metric in this limit are Let (, log , , ) be eccentric coordinates specialized for use in our simulations; they are related to cylindrical coordinates by Here  is the eccentricity of our eccentric coordinates, set to 0 for circular disks and 0.5 for eccentric disks.Coordinate surfaces of constant  are elliptical cylinders of semilatus rectum , or equivalently, semimajor axis  = /(1 −  2 ).We opt for log  instead of  in order to generate a logarithmic grid, but for ease of understanding we continue to label that coordinate by  and express results in terms of .We reuse  and  without risk of ambiguity because these two coordinates do not participate in the coordinate transformation from cylindrical to eccentric.The metric and connection in both coordinates are provided in Appendix A.

Gravitational potential and orbits
We ignore vertical gravity in these first simulations of the eccentric MRI, so the potential depends only on .In this sense, our simulations resemble earlier simulations of unstratified circular disks (e.g., Hawley & Balbus 1991;Hawley et al. 1995Hawley et al. , 1996;;Sano et al. 2004;Fromang & Papaloizou 2007;Fromang et al. 2007;Lesur & Longaretti 2007;Simon et al. 2009;Guan et al. 2009;Bodo et al. 2011).However, instead of the point-mass gravitational potential Φ(, ) = −1/, we adopt because, as proven in Appendix B, this potential admits closed eccentric orbits at all distances.The modification matters because general-relativistic apsidal precession, albeit small at large distances, accumulates over the tens of orbits during which the MRI saturates.
The properties of orbits in our potential are also derived in Appendix B; here we repeat the parts that support our exposition.For an orbit along a coordinate curve of semilatus rectum  or semimajor axis  = /(1 −  2 ), the orbital period is The specific energy and angular momentum conserved with respect to our metric are The specific energy includes the rest energy;  = 1 corresponds to marginally bound material, and  → 0 as material becomes more bound.The nonzero components of the orbital velocity are so the physical velocity at pericenter is Consider a collection of such orbits nested within one another, all following coordinate curves.The velocity field thus generated has finite divergence: Consequently, if the motion of a gas without pressure is described by this velocity field, the density along a streamline of constant  cannot be uniform, but must instead vary with  in proportion to So as to guarantee an initial condition that is a genuine hydrodynamical steady state, we take this modulation into account in Section 2.5 even though the modulation amplitude is tiny for our disk parameters.

Hydrodynamic initial condition
Because we ignore vertical gravity, our initial disk is translationally symmetric along the -direction.It can be described as an elliptical annular cylinder, each shell of which orbits the coordinate axis along a coordinate surface of constant  with a velocity as given by Equations ( 15) and ( 16).The scale of the disk is characterized by its fiducial orbit, whose semilatus rectum is the geometric mean of the semilatera recta of its inner and outer edges.The disk should be large enough that orbital velocities are non-relativistic, and small enough that severe truncation errors do not arise from the evolution of the total energy as a result of the smallness of the kinetic energy with respect to the rest energy.In light of the fact that the linear growth rate of the MRI is inversely proportional to the orbital period (Chan et al. 2018), we additionally require that runs of different  have the same semimajor axis and thus orbital period.We settle on a semilatus rectum of  * = 200(1 −  2 ) for the fiducial orbit, and we report time in units of the orbital period at this orbit.
The initial density profile is with  * the density at the pericenter of the fiducial orbit, (, ) = ( * , 0).To give the disk edges that are not too sharp and to build in a numerical vacuum, the density is modulated spatially as where is a smoothed top-hat function and is the fractional (log )-position within the simulation domain.
The modulation is governed by three parameters:  v = 0.01 is the vacuum-to-disk density ratio, and  gb = 0.5 and  gt = 0.1 are the fractional (log )-extents of the disk body and the disk-vacuum transition, respectively.
Because the most unstable mode of the circular MRI is incompressible, the thermodynamic properties of the gas are expected to have little bearing on the growth of the eccentric MRI.Even so, we would like the pressure to be small enough initially that the configuration above is close to equilibrium, and the internal energy to be large enough at all times that catastrophic truncation errors are avoided.We balance these competing desires by setting the Mach number at the fiducial pericenter to  * = 30.Additionally, we adopt an adiabatic index of  = 1 + 10 −5 , corresponding to a nearly isothermal gas, so that the internal energy is larger at fixed pressure.The initial pressure is therefore at the fiducial pericenter and everywhere else.This initial condition is not strictly hydrostatic; transient outgoing waves are launched from the inner edge as the disk seeks force balance.Moreover, our use of a soft equation of state means that density perturbations are stronger and pressure gradients have a lesser effect on disk evolution than if an adiabatic equation of state were used.

Magnetic initial condition
The magnetized runs are initialized with the two kinds of magnetic topologies illustrated in Figure 1.The vertical-field topology refers to a magnetic field with one nonzero component where  mb = 0.4 and  mt = 0.1 are the fractional (log )-extents of the magnetized disk body and the transition from the magnetized disk body to the unmagnetized disk edges, respectively.The net vertical magnetic flux persists throughout the simulation.Thanks to its simplicity and its ability to generate the fastest-growing instabilities, the vertical-field topology was considered in the first studies of the circular MRI (Balbus & Hawley 1991;Hawley & Balbus 1991) and in our analytic study of the eccentric MRI (Chan et al. 2018).The dipolar-field topology is derived from a magnetic potential with one nonzero component where is the fractional -position within the simulation domain and  mb = 0.4.The inclusion of the metric determinant makes the magnetic-field strength more uniform over azimuth.The dipolarfield topology sees frequent application in global simulations of the circular MRI (e.g., Hawley 2000).
For both topologies, the initial plasma beta, defined as the ratio of the initial volume integral of gas pressure to that of magnetic pressure, is 100.The pressure due to the magnetic field is subtracted from the gas to preserve the total.The gas pressure in magnetized regions is perturbed at the 0.01 level to seed the MRI.Periodic boundary conditions apply to the -direction.They are also employed in the -direction, in accordance with our neglect of vertical gravity.Outflow boundary conditions are used in the -direction: We copy all quantities to the ghost zone, zero the -component of the velocity if it points into the simulation domain, and zero the -and -components of the magnetic field always.This last step reduces unphysical influences from the boundaries.
To prevent numerical issues, we require that the pressure in the simulation domain always satisfy ( /) 1/2 ≥ 2 × 10 −4 .In addition, whenever the recovery of primitive variables fails, the primitive variables from the previous time step are carried forward.
As the simulation progresses, different parts of the disk may evolve differently in eccentricity and orientation, so the disk could eventually become misaligned with the grid, resulting in greater numerical dissipation.The disk also occupies a wider range of semilatera recta due to pressure gradients and outward angular momentum transport, bringing the now differently shaped disk into contact with the boundaries.The inner boundary poses a lower limit on the pericenter, but this restriction is arguably physical because there are indeed radii an accretion flow cannot return from.If the central object is a star, the disk cannot extend inside the star or its magnetosphere.If the central object is a black hole, Figure 2 informs us that material with sufficiently low angular momentum and high energy can plunge directly into the black hole; low angular momentum and high energy are, in Schwarzschild spacetime as a function of spherical radius  and specific angular momentum .The specific energy  must satisfy 1 2 ( 2 − 1) >  eff , and marginally bound trajectories have  = 0. Material accretes by losing angular momentum; thus, its trajectory is described by potentials of decreasing .In circular disks, trajectories have the lowest energy allowed by the stable potential well; such trajectories evolve along the sequence of dots downward until they arrive at the smallest radius that supports circular orbits, the ISCO.In eccentric disks, trajectories have energies above the bottom of the potential well, which allows them to make radial excursions as suggested by the double-headed arrow.The smallest radius a bound trajectory can reach without falling in is the marginally bound orbit; such a trajectory has  = 0 and  2 = 16.Because trajectories energetic enough to overcome the centrifugal barrier have energies greater than an ISCO orbit, material plunging into the black hole on these trajectories has less energy available for radiation compared to material accreting on circular trajectories.
of course, the hallmarks of an eccentric orbit.By contrast, the interaction with the outer boundary is unphysical and can lead to numerical artifacts; therefore, we restrict our attention to the first 15 orbits at  = 200, before the disk starts running into the outer boundary and numerical artifacts appear.Steady state appears to obtain at this time for the plasma beta and the alpha parameter, despite the continual evolution in disk shape.

Overview
Figure 3 tells us how much the disks have changed by the end of the simulations.The circular unmagnetized disk after 15 orbits is almost indistinguishable from the initial disk.The eccentric unmagnetized disk changes shape slightly because pressure gradients, which are nonzero along the inner and outer disk edges, cause differential apsidal precession (e.g., Ogilvie 2001); the density striation is a result of this adjustment process.
The most conspicuous contrast between unmagnetized and magnetized disks, of whatever eccentricity, is that unmagnetized disks remain smooth while magnetized disks develop large density fluctuations.This, of course, is due to the MRI creating MHD turbulence in the magnetized disks.The fluctuations are larger in vertical-field disks than in dipolar-field disks.
Comparison can also be made between circular and eccentric disks regardless of magnetic topology.Unlike the circular disks, which stay circular despite the MHD turbulence, the inner parts of eccentric disks grow more eccentric and precess prograde.
As the inner parts of eccentric disks shrink, their pericenters move inside the inner boundary, and their material accretes across the boundary while retaining its eccentricity.This loss of material from the most eccentric orbits is the reason why there is a sparsely populated region between the inner parts and the inner boundary, visible in the late-time eccentric disks in Figure 3.
There are also differences among the eccentric disks, concerning chiefly eccentricity evolution and to a lesser degree precession.The unmagnetized and vertical-field disks in eccentricity gradient and a much higher degree of precession.
The eccentricity gradient can be quantified by computing the instantaneous eccentricity ē, which is the eccentricity of the orbit material would follow given its instantaneous velocity   if only gravitational forces act; this orbit is also known as the osculating orbit.We calculate ē from   using Equations ( 15), ( 16), and (B9).Figure 4 contains plots of the mass-weighted vertical average of ē: The instantaneous eccentricity deviates little from its initial value in the unmagnetized and vertical-field disks, but it develops a clear gradient in the dipolar-field disk, with the eccentricity higher than its initial value in the inner parts and lower in the outer parts.The implications of the eccentricity gradient will be discussed in Section 4.2.

Plasma beta and the alpha parameter
The top half of Figure 5 portrays the mass-weighted vertical average of the plasma beta at the end of the simulations, defined as The plasma beta is ∼ 10 in vertical-field disks and ∼ 100 in dipolar-field disks; the variation within a disk is about one order of magnitude.Comparable levels of plasma beta are witnessed in circular and eccentric disks with the same magnetic topology, suggesting that the MRI is unimpeded by eccentricity, in contrast to the findings of Dewberry et al. (2020a).The stronger magnetic fields produced by the vertical-field topology accord with simulations of circular disks in the literature (e.g., Hawley et al. 1995Hawley et al. , 1996;;Sano et al. 2004;Bai & Stone 2013).30) and (37), respectively.The panels are the outcomes of imposing various magnetic topologies on an initial disk and evolving it for 15 Regions with negligible levels of magnetic field have exceedingly large values of plasma beta.The boundaries of the simulation domain are traced by thin ellipses.
We can also examine the plasma beta along a one-dimensional profile running from the inside of the disk to the outside at different times during the simulations.Considering that the disk evolves in eccentricity and orientation, it makes little sense to look at profiles over semilatus rectum; instead, we construct profiles over cylindrical radius.The mass-weighted average plasma beta at cylindrical radius  is given by where the hypersurface element is Temporal and spatial averaging smooths out turbulent fluctuations.Temporal averaging is performed over two intervals, each lasting one-third of the simulation duration; comparison between the intervals gives us an idea how close the MRI is to saturation.Spatial averaging is limited to the cylindrical shell of radius  picked out by the delta-function.The results are plotted in the top half of Figure 6, and the legend lists the intervals of temporal averaging.The relatively small difference between the two intervals at radii 100 ≲  ≲ 200 suggests that steady state is achieved to some degree at those radii, despite the relatively short simulation duration.We also reach similar conclusions as we did with Figure 5: The plasma beta is quite uniform over the disk, it is not significantly modified by the introduction of eccentricity, but it is one to two orders of magnitude lower in vertical-field disks than in dipolar-field disks.The alpha parameter is conventionally taken to be the sum of Reynolds and Maxwell stresses divided by the gas pressure.However, it is difficult to determine the mean flow and departures from it in a disk whose inner and outer parts evolve differently in eccentricity and orientation.We therefore consider only the Maxwell, not Reynolds, contribution to the alpha parameter, working under the assumption that the Maxwell stress dominates the total, as is uniformly the case for circular disks (e.g., Hawley et al. 1995).The Maxwell stress is defined in terms of the projected magnetic field   , similar to the stress-energy tensor.To make the stress more physically interpretable, we measure   in a local orthonormal basis, whose basis vectors are those of cylindrical coordinates but with lengths normalized to unity: We then define the Maxwell stress to be − φ  R. A factor of  is attached to   in Equation ( 33) because our eccentric coordinates use log , not .Orthonormality guarantees that The bottom half of Figure 5 depicts the mass-weighted vertical average of the Maxwell-only alpha parameter at the end of the simulations: The alpha parameter is positive almost everywhere in circular disks, as required for outward angular momentum transport.By contrast, the alpha parameter is uniformly positive in the lower half of the eccentric disks where material falls to pericenter, but it is consistently negative in certain sectors of the upper half where material flies out to apocenter.In the eccentric vertical-field disk, positive sectors occupy significantly more area than negative sectors.The magnitude of the alpha parameter, whether positive or negative, varies from ∼ 0.2 to 5; however, if we construct areaweighted histograms of the magnitude separately for positive and negative sectors, we find that the positive histogram is shifted by a factor of ∼ 2 toward larger values relative to its negative counterpart.In the eccentric dipolar-field disk, the total area of positive sectors is only slightly larger than that of negative sectors.In addition, the magnitude of the alpha parameter has a narrower distribution, ∼ 0.3 to 3; the positive histogram is again displaced by a factor of ∼ 2 compared to the negative one.We shall speculate about why the alpha switches sign in Section 4.1.
The net effect of angular momentum transport is revealed by integrating the Maxwell stress over an orbit.When doing so by eye on the basis of Figure 5, it is important to take into account the relative areas and alpha-parameter ranges of the positive and negative sectors.More quantitatively, we construct in Figure 6 mass-weighted radial profiles of the alpha parameter at different times using the prescription with  the same hypersurface element from Equation (32).The alpha parameter is comparable in circular and eccentric disks with the same magnetic topology.In vertical-field disks, the alpha parameter is ∼ 0.5 to 1; in dipolar-field disks, it is still positive, but only ∼ 0.05 to 0.15.Stronger stress for vertical than dipolar magnetic field agrees with previous simulations of circular disks (e.g., Hawley et al. 1995Hawley et al. , 1996;;Sano et al. 2004;Bai & Stone 2013).

Specific angular momentum and binding energy
To investigate the effect of internal stresses, we examine how the specific angular momentum squared  2 and binding energy  b = 1 −  evolve.The mass-weighted vertical averages of the two quantities at the end of the simulations are Here  φ =   is the velocity measured in the local orthonormal cylindrical basis, defined analogously to  φ in Equation ( 34), and the covariant velocity component   = − is a conserved quantity of our time-independent metric.Figure 7 displays the specific angular momentum squared and binding energy, normalized to their initial values at the inner edge.Our focus is on the eccentric disks.Because of pressure gradients, even the inner edge of the unmagnetized disk experiences a reduction in specific angular momentum squared by ∼ 25% from its initial value.The decrease is greater in magnetized disks where the Maxwell stress also contributes: The dipolar-field disk reports a drop by ∼ 33%, and the verticalfield disk, which has a smaller plasma beta and larger alpha parameter than the dipolar-field disk, records a suppression by ∼ 50%.It should be noted that the range of specific angular momentum is constrained by the geometry of the simulation domain: Once material loses enough angular momentum that its pericenter recedes inside the inner boundary, it leaves the simulation domain.
The specific binding energy at the inner edge changes by rather less.It is almost identical to its initial value in the unmagnetized disk, ∼ 30% higher in the vertical-field disk, and ∼ 10% lower in the dipolar-field disk.The smaller fractional changes suggest that the torques transporting angular momentum outward, regardless of whether they are hydrodynamic or magnetic in nature, occur preferentially near apocenter where the attendant work done is smaller.

Quality factors
We close this section by examining how well our magnetized disks resolve the MRI.The figures of merit for circular disks are the quality factors, defined as the ratio of the characteristic wavelength of the MRI to cell sizes in different directions (Hawley et al. 2011).We generalize their mass-weighted vertical averages to eccentric disks as  39) and ( 40), respectively.The top row shows two initial disks with different eccentricities.The panels under each top-row panel are the outcomes of imposing various magnetic topologies on an initial disk and evolving it for 15 orbits.In all panels, the inner edge is arbitrarily defined to be where the density is 0.1  * , with  * the fiducial density from Equation ( 20), and regions less dense than that are left blank.Furthermore, the specific angular momentum squared or binding energy in each panel is normalized by its value at the inner edge of the corresponding initial disk, which is why the inner edge appears yellow in the top row.The boundaries of the simulation domain are traced by thin ellipses.The fractional changes in the specific angular momentum squared and binding energy together determine the change in eccentricity seen in Figure 4.
where  is the orbital period from Equation ( 12), and Δ and Δ are the cell sizes in the -and -directions, respectively.In keeping with work on circular disks, we base our quality factors on the physical, cylindrical components  μ of the projected magnetic field   , given by Equations ( 34) and (35).
Figure 8 showcases the quality factors at the end of the simulations.For ease of comparison with the criteria that indicate adequate resolution for circular disks, to wit, ⟨  ⟩ ≳ 15 and ⟨  ⟩ ≳ 20 (Hawley et al. 2013), the color scales of the figure are centered on these values.In terms of ⟨  ⟩, the vertical-field disks are extremely well-resolved everywhere, but the same is true for the dipolar-field disks only for a limited range of semilatera recta.In terms of ⟨  ⟩, all disks are marginally under-resolved.Interestingly, ⟨  ⟩ is noticeably higher in the pericentric half of the disks; this is because the strong orbital shear near pericenter amplifies  φ , not  ẑ .

Stresses in circular and eccentric disks
Our earlier investigation into the linear stage of the MRI found that the MRI grows in both circular and eccentric disks; the  41) and ( 42), respectively.The panels are the outcomes of imposing various magnetic topologies on an initial disk and evolving it for 15 orbits.The boundaries of the simulation domain are traced by thin ellipses.
growth rate in eccentric disks is about half that in circular disks (Chan et al. 2018).The present simulations suggest that the nonlinear stage of the MRI is also not so different in circular and eccentric disks, in that it saturates to comparable levels of plasma beta and alpha parameter in the two kinds of disks.Thus, the MRI functions in much the same capacity in eccentric disks as in circular disks, mediating outward angular momentum transport.That being said, the MRI in eccentric disks exhibits two intriguing features not witnessed in circular disks.The first one is the sign flip in the Maxwell-only alpha parameter.In circular disks, the alpha parameter is positive almost everywhere; by contrast, in eccentric disks, the alpha parameter can be consistently negative in some sectors of the disk even though the azimuthal integral of the Maxwell stress remains positive.
To understand the basic principle underlying this surprising behavior, we break temporarily from the general-relativistic treatment used in the rest of the article and work in the Newtonian limit.We denote the Newtonian velocity and magnetic field by ṽ and B, respectively, and their components measured against the local orthonormal cylindrical basis by ṽ and B .We conflate B with  î elsewhere in the text, so the Maxwell stress is simply − B B .Our starting point is the induction equation: To track the time evolution of the magnetic field at a point comoving with the flow, we consider We set ∇ • ṽ = 0 because the MRI is largely incompressible, and we take / = 0 because we ignore vertical gravity in our simulations.Consequently, In circular disks, the time-and azimuth-averaged  ṽ  / is negative, the averaged (1/) ( ṽ /) is zero, and the averaged ṽ is also negative, albeit with a very small magnitude.Because − B B > 0 nearly everywhere, all terms on the right-hand side of Equation ( 45) are therefore on average zero or positive, so − B B grows over time until limited by dissipation.
In eccentric disks, the time-averaged  ṽ  / is also negative everywhere.However, the averaged (1/) ( ṽ /) is no longer zero: It changes from negative to positive shortly before pericenter and back after pericenter.The averaged ṽ also changes sign, from negative to positive at pericenter and back at apocenter.The second and third terms therefore have azimuthdependent signs.When orbits are significantly eccentric, the averaged (1/) ( ṽ /) is comparable in magnitude to the averaged  ṽ  /, and the averaged ṽ is comparable to the averaged ṽ  , so all three terms can be important.
Consider material starting from apocenter with positive Maxwell stress, that is, − B B > 0. Initially all three terms work together to make − B B more positive.As the material swings toward pericenter and its velocity becomes more azimuthal than radial, (1/) ( ṽ /) turns positive while | ṽ | drops, flipping the sign of the second term and reducing the magnitude of the third term.If | B | ≲ | B |, which is marginally satisfied at these azimuths, the second term can dominate, pulling − B B toward negative values.For some distance past pericenter, the second term continues to be negative; at the same time, ṽ > 0 and  ṽ / > 0, so the third term now tends to reduce the magnitude of − B B at an increasing rate regardless of its sign.The second and third terms combined make − B B < 0 at some point near pericenter.Further beyond pericenter, the second term changes back to positive and the third term decreases in magnitude.By the time the material returns to apocenter, the first and second terms have restored − B B to positive, transporting angular momentum outward.
Both magnetic topologies share these qualitative features, but they differ in the quantitative details.The dipolar-field topology yields a very close alignment between the local directions of velocity and magnetic field, whereas the vertical-field topology produces merely a correlation between the two directions that has much more scatter.In addition, as we can see in Figure 5, the sectors of the dipolar-field disk with − B B < 0 coincide with the region where ṽ > 0, whereas the vertical-field disk sectors with − B B < 0 take up only a small part of the ṽ > 0 region, and their boundary is much more irregular.We speculate that both of these contrasts are due to the larger amplitude of turbulence the MRI drives in the presence of a vertical magnetic field, an effect well-documented in circular disks (e.g., Hawley et al. 1995Hawley et al. , 1996;;Sano et al. 2004;Bai & Stone 2013).Because the third term of Equation ( 45) is sensitive to correlations between fluctuations in velocity and magnetic field, large-amplitude turbulence may change the stress evolution.This supposition is corroborated by the fact that the vertical root-mean-square of   , a good proxy for the magnitude of turbulent fluctuations in unstratified simulations, in the vertical-field disk is ∼ 4 times that in the dipolar-field disk.
The other curiosity is that, whereas MHD stresses in circular disks transport angular momentum and energy at such rates as to keep shrinking orbits circular, the two rates can be independent of each other in eccentric disks.Circular disks have  b  2 ≈ 1 2 at all radii, according to Equation (B10); by contrast,  b and  2 in eccentric disks vary by amounts that depend on location, and the mismatch dictates how the local eccentricity evolves.The fractional changes of  b and  2 estimated in Section 3.3 imply that the unmagnetized disk has the least eccentric inner edge and the dipolar-field disk has the most eccentric, which agrees with Figure 4.One possible explanation of the ranking relies on the fact that stresses near pericenter tend to lower the eccentricity of the inner edge and stresses near apocenter tend to raise it (Svirski et al. 2017).Stresses in the unmagnetized disk, being purely hydrodynamical, should be concentrated near pericenter where pressure gradients are the steepest; conversely, the Maxwell stress in the vertical-and dipolar-field disks are relatively evenly distributed between pericenter and apocenter, as made apparent by the bottom half of Figure 5.Our speculation in the previous paragraph may bear on the ordering of the two magnetized disks: Stronger turbulence in the vertical-field disk disrupts coherent angular momentum transport, so the inner edge is less eccentric than in the dipolar-field disk.

Dynamics and energy dissipation at the inner edge
As demonstrated in Section 3.3, MHD stresses in eccentric disks are typically more effective at moving angular momentum than energy.Because angular momentum constrains the size of the inner edge, and because the binding energy at the inner edge determines the amount of energy available for radiation, the different transport rates of these two quantities could have observable effects on eccentric disks.
To explore these effects, let us first consider how disks behave around a star.For a circular disk, the inner edge is located at the larger of the stellar or Alfvén radius.The total energy a given amount of material dissipates in the disk is its binding energy on a circular orbit at that radius.Additional dissipation happens in the boundary layer at the inner edge as the material comes into corotation with the star or its magnetosphere.If the star-regulated rotation speed in this boundary layer is small compared to the orbital speed, the additional dissipation is equal to the dissipation that took place in the disk itself.
For an eccentric disk, the stellar or Alfvén radius likewise defines the inner edge, but in this case, it is through a match to the pericenter of the inner edge.Because the corresponding semimajor axis is larger than the stellar or Alfvén radius, material dissipates less energy in the disk proper.The total energy dissipated, however, is exactly the same as in the circular case because the eccentric orbit of the material ultimately transforms into a circular one at the stellar or Alfvén radius.What makes eccentric disks different is that dissipation in the boundary layer accounts for a larger fraction of the total.Energy can be dissipated there in a variety of ways, depending on how the inner edge interacts with the star or its magnetosphere.In one extreme, the material immediately comes into corotation.The fast-moving material dissipates large amounts of energy in a small region, potentially producing hard radiation.In the other extreme, the material loses a tiny part of its kinetic energy each time it grazes past the star or its magnetosphere, and the material migrates inward gradually on orbits of decreasing semimajor axes and eccentricities.Dissipation in this case could be more spatially distributed, happening both in the grinding encounters and along the circularizing orbits; if so, the eccentric disk may resemble a circular one in appearance.
We now turn to disks around black holes.As shown in Figure 2, material on eccentric orbits can reach the event horizon even if it has more angular momentum than that of an ISCO orbit provided that it has enough energy to overcome the centrifugal barrier inside the ISCO radius; equivalently, it can do so if its eccentricity is high enough.Moreover, when the angular momentum is at least that of an ISCO orbit, the peak of the centrifugal barrier is always above the energy of an ISCO orbit.Consequently, the amount of energy available for such material to radiate is always smaller than the binding energy at the ISCO, the conventional estimate for radiative efficiency in circular disks (Novikov & Thorne 1973;Page & Thorne 1974).The energy available to eccentric material is even smaller when we take into account extra energy extraction by the Maxwell stress in circular disks (Thorne 1974;Krolik 1999;Gammie 1999;Krolik & Hawley 2002;Noble et al. 2009;Avara et al. 2016;Schnittman et al. 2016;Kinch et al. 2021).Unlike eccentric disks around stars, the energy retained by the plunging material cannot be recovered for radiation in some boundary layer.
The diminution of the radiative efficiency of eccentric disks around black holes could be particularly relevant to TDEs, in which the energy radiated is considerably below that expected from the accretion of a reasonable fraction of a stellar mass of material onto a black hole through a conventional circular disk (Piran et al. 2015).Svirski et al. (2017) suggested that this "inverse energy problem" could be resolved by internal stresses that transport angular momentum preferentially and cause the most bound debris to plunge; our results provide quantitative support to that qualitative argument.

Effects of vertical gravity
Our simulations of the nonlinear development of the MRI in eccentric disks ignored vertical gravity.We did so for two reasons.The first is to isolate the effects of orbital-plane motions, the chief driver of the MRI, when they are azimuthally modulated in eccentric disks.The second is to produce a baseline for comparison with future simulations that do include vertical gravity.
Because vertical gravity is a form of tidal gravity, it is much weaker than radial gravity as long as the disk is thin.For this reason, the modifications it introduces are likely minor and unable to quench the MRI.Here we examine what these modifications may be.
One effect of vertical gravity that is discernible in circular disks carries over to eccentric disks: It regulates the buoyant eviction of accumulated magnetic flux that may determine the saturation level of the MRI (e.g., Tout & Pringle 1992;Davis et al. 2010;Beckwith et al. 2011;Hirose & Turner 2011;Begelman et al. 2015;Hogg & Reynolds 2018).
The variation of vertical gravity around an eccentric orbit causes the disk to expand in height as it travels from pericenter to apocenter and to collapse on the way back.This "breathing" can have large amplitudes even for mildly eccentric disks (Ogilvie & Barker 2014;Lynch & Ogilvie 2021).For the thickest and most eccentric disks, breathing can even become non-adiabatic: Extreme compression can create shocks that eject material vertically (Ryu et al. 2021).
Vertical gravity can also influence the eccentric MRI more subtly, without shocks.Vertical motion opens up more avenues of energy exchange between the orbit and the magnetic field.Compression and expansion can alter the Alfvén speed, hence the critical wavelength for stability against the MRI, and it can also modulate the wavelengths of advected perturbations.All these variations are periodic, raising the possibility of instability through parametric resonance (Papaloizou 2005b).
It bears repeating that the changes due to vertical gravity should be small for thin, moderately eccentric disks, so we expect our results here to remain generally valid.

CONCLUSIONS
In circular disks, the MRI stirs up correlated MHD turbulence, turbulent stresses transport angular momentum outward, and the disk accretes.Our simulations demonstrate that much the same process operates in eccentric disks, and with comparable efficiency.Like circular disks, the quantitative level of Maxwell stress achieved in eccentric disks depends on the magnetic topology, but eccentric and circular disks with the same magnetic topology reach comparable levels of plasma beta and alpha parameter.
Although the mass-weighted, disk-averaged Maxwell stress in an eccentric disk produces an outward angular momentum flux similar in magnitude to that in a circular disk with the same magnetic topology, the Maxwell stress in an eccentric disk can cause inward angular momentum transport in certain disk sectors.This behavior is seen only in eccentric disks likely because the radial velocity and its azimuthal gradient are both nonzero and have signs that vary over azimuth.
By removing proportionately more angular momentum than energy from the inner parts of the disk, MHD stresses can promote accretion of highly eccentric material.This material has higher energy than a circular orbit with the same angular momentum, so the radiative efficiency of a disk around a black hole can be suppressed relative to a circular disk.These findings corroborate earlier suggestions about how the power output from eccentric disks may explain why the total observed radiated energy in many TDEs is one to two orders of magnitude lower than expected from circular accretion of a stellar mass of material (Svirski et al. 2017).

Figure 1 .
Figure 1.Poloidal slices showing with arrows the initial magnetic field for our two magnetic topologies.Background colors plot the initial density on the same scale as Figure 3.The ordinate is more stretched than the abscissa, and the vectors are arbitrarily scaled.

Figure 2 .
Figure 2. Effective potential  eff = −1/ + 1 2 ( 2 / 2 ) (1−2/ ) in Schwarzschild spacetime as a function of spherical radius  and specific angular momentum .The specific energy  must satisfy 12 ( 2 − 1) >  eff , and marginally bound trajectories have  = 0. Material accretes by losing angular momentum; thus, its trajectory is described by potentials of decreasing .In circular disks, trajectories have the lowest energy allowed by the stable potential well; such trajectories evolve along the sequence of dots downward until they arrive at the smallest radius that supports circular orbits, the ISCO.In eccentric disks, trajectories have energies above the bottom of the potential well, which allows them to make radial excursions as suggested by the double-headed arrow.The smallest radius a bound trajectory can reach without falling in is the marginally bound orbit; such a trajectory has  = 0 and  2 = 16.Because trajectories energetic enough to overcome the centrifugal barrier have energies greater than an ISCO orbit, material plunging into the black hole on these trajectories has less energy available for radiation compared to material accreting on circular trajectories.

Figure 3 .
Figure 3. Midplane slices of density, in units of the fiducial density  * from Equation (20).The top row shows two initial disks with different eccentricities.The panels under each top-row panel are the outcomes of imposing various magnetic topologies on an initial disk and evolving it for 15 orbits.The boundaries of the simulation domain are traced by thin ellipses in order to better distinguish low-density regions from regions not covered by the simulations.Short lines from the origin indicate the approximate orientations of the inner parts of the eccentric disks.

Figure 4 .
Figure 4. Mass-weighted vertical average of the instantaneous eccentricity, as defined in Equation (29).The top-left panel shows the eccentric initial disk.The other panels are the outcomes of imposing various magnetic topologies on the initial disk and evolving it for 15 orbits.The boundaries of the simulation domain are traced by thin ellipses.In keeping with Figure 7, regions with density below 0.1  * are left blank, where  * is the fiducial density from Equation (20).

Figure 5 .
Figure 5. Mass-weighted vertical averages of the plasma beta in the top half and Maxwell-only alpha parameter in the bottom half, as defined in Equations (30) and (37), respectively.The panels are the outcomes of imposing various magnetic topologies on an initial disk and evolving it for 15 Regions with negligible levels of magnetic field have exceedingly large values of plasma beta.The boundaries of the simulation domain are traced by thin ellipses.

Figure 6 .
Figure 6.Mass-weighted radial profiles of the plasma beta in the top half and Maxwell-only alpha parameter in the bottom half, as defined in Equations (31) and (38), respectively.The profiles are time-averaged over the two intervals in the legend.The two rows in the bottom half have different vertical scales.The magnetic field strength varies weakly with eccentricity, but is much stronger in vertical-field disks than in dipolar-field disks.

Figure 7 .
Figure 7. Mass-weighted vertical averages of the specific angular momentum squared in the left half and specific binding energy in the right half, as defined in Equations (39) and (40), respectively.The top row shows two initial disks with different eccentricities.The panels under each top-row panel are the outcomes of imposing various magnetic topologies on an initial disk and evolving it for 15 orbits.In all panels, the inner edge is arbitrarily defined to be where the density is 0.1  * , with  * the fiducial density from Equation (20), and regions less dense than that are left blank.Furthermore, the specific angular momentum squared or binding energy in each panel is normalized by its value at the inner edge of the corresponding initial disk, which is why the inner edge appears yellow in the top row.The boundaries of the simulation domain are traced by thin ellipses.The fractional changes in the specific angular momentum squared and binding energy together determine the change in eccentricity seen in Figure4.

Figure 8 .
Figure 8. Mass-weighted vertical averages of the vertical quality factor in the top half and azimuthal quality factor in the bottom half, as defined in Equations (41) and (42), respectively.The panels are the outcomes of imposing various magnetic topologies on an initial disk and evolving it for 15 orbits.The boundaries of the simulation domain are traced by thin ellipses.