# Linear and nonlinear eccentric mode evolution in unstratified MHD discs

Elliot M. Lynch,1 Janosz W. Dewberry,2
1Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230, Saint-Genis,-Laval, France
2 Canadian Institute for Theoretical Astrophysics, 60 St. George Street,
Toronto, ON M5S 3H8, Canada,
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this paper we develop a framework for studying unstratified, magnetised
eccentric discs and compute uniformly precessing eccentric modes in a
cylindrical annulus which provide convenient initial conditions for numerical
simulations. The presence of a magnetic field in an eccentric disc can be
described by an effective gas with a modified equation of state. At magnetic
field strengths relevant to the magneto-rotational instability the magnetic
field has negligible influence on the evolution of the eccentric disc, however
the eccentric disc can significantly enhance the magnetic field strength over
that in the a circular disc. We verify the suitability of these eccentric disc
solutions by carrying out 2D simulations in RAMSES. Our simulated modes (in
2D) follow a similar evolution to the purely hydrodynamical modes, matching
theoretical expectations, provided they are adequately resolved. Such
solutions will provide equilibrium states for studies of the eccentric
magneto-rotational instability and magnetised parametric instability in
unstratified discs and are useful for exploring the response of disc
turbulence on top of a fluid flow varying on the orbital timescale.

###### keywords:

accretion, accretion discs – MHD – magnetic fields – celestial mechanics

††pubyear: 2023††pagerange: Linear and nonlinear eccentric mode evolution in
unstratified MHD discs–B.2

##  1 Introduction

Eccentric gaseous discs, where the gas orbits on Keplerian ellipses, are found
in a variety of astrophysical contexts. To date there have been many
theoretical and numerical studies considering unmagnetised eccentric discs
(Ogilvie, 2001; Ogilvie & Barker, 2014; Barker & Ogilvie, 2016; Wienkers &
Ogilvie, 2018; Ogilvie & Lynch, 2019; Pierens et al., 2020; Dewberry et al.,
2020a). Recently several studies have considered the behaviour of magnetic
fields in eccentric discs. The effect of magnetic stresses on eccentric discs
was considered by Ogilvie (2001) who developed a turbulent stress model based
on the ideal induction equation in an orbital coordinate system. Ogilvie &
Barker (2014) include an magnetohydrodynamic (MHD) form of their eccentric
shearing box model. This was used to study the linear phase of the magneto-
rotational instability (MRI) by Chan et al. (2018), while Lynch & Ogilvie
(2021) used the formalism to study the effect of a coherent magnetic field on
the disc vertical structure. Global simulation of MRI in eccentric discs were
performed by Dewberry et al. (2020b) who found that sufficiently nonlinear
eccentric waves can shut off the MRI. Oyang et al. (2021) compared the
excitation of eccentricity in MRI turbulent discs to viscous, hydrodynamical
discs and found that the latter were excited to larger eccentricities. Finally
Chan et al. (2022) performed a global simulation of the MRI in an elliptical
annulus of large ($e=0.5$) constant eccentricity motivated by the highly
eccentric discs found in tidal disruption events.

One challenge for (hydro or MHD) simulations of eccentric discs is the strong
differential precession due to pressure forces which arises for arbitrary
eccentricity profiles. Notably this occurs for a uniformly eccentric ring
which, naively, might be considered the simplest eccentricity profile to
simulate. This strong differential precession was a problem encountered by
Chan et al. (2022) who utilised an elliptical coordinate system to model a
disc of uniform eccentricity, which became significantly misaligned from the
simulation grid after only 15 outer disc orbits. Strong differential
precession quickly generates large pressure gradients as a result of orbital
compression, which can be very difficult to resolve numerically, leading to
artificial damping of the disc eccentricity. The strength of this differential
precession will depend on the magnetic field strength and configuration. The
rapid evolution of the disc orbits can potentially lead to transient phenomena
that are primarily a consequence of the choice of initial conditions. This
makes it difficult to disentangle the effect of the MRI and parametric
instability from the evolution of the non-steady initial conditions. It would
thus be beneficial to study how the MRI develops on top of a steady, or slowly
evolving, eccentric background.

A solution to this problem can be found in the existence of eccentric modes.
These are untwisted eccentric discs with a time independent eccentricity
profile which undergo uniform (i.e. rigid body) precession as a result of
pressure gradients and other non-Keplerian forces. Eccentric modes are thus a
particularly suitable setting for numerical simulations. They are also well
motivated physically as they often provide a good approximation to the relaxed
state of many eccentric discs when excitation and damping processes are
considered (Kley et al., 2008; Miranda et al., 2017; Teyssandier & Ogilvie,
2016, 2017; Ragusa et al., 2017).

In this paper we extend the Hamiltonian eccentric disc theory of Ogilvie &
Lynch (2019) to allow for the inclusion of a large scale, structured, magnetic
field in an unstratified disc. We calculate modal (uniformly precessing)
solutions for these eccentric MHD discs. Such solutions are not intended as a
realistic model of a magnetised eccentric disc, owing to the neglect of
important 3D effects (Ogilvie, 2001, 2008; Ogilvie & Barker, 2014; Teyssandier
& Ogilvie, 2016; Ogilvie & Lynch, 2019) and unrealistic global field structure
(see Ogilvie, 1997, for the 3D field structure in a circular disc), however
these solutions are intended to provide a convenient setting for numerical
simulations of the eccentric MRI. To this end we run 2D MHD simulations in the
code RAMSES (Teyssier, 2002; Fromang et al., 2006; Faure et al., 2014) using
our calculated eccentric mode as an initial condition to test their
suitability for numerical calculations. By using an eccentric mode as our
initial condition we aim to avoid strong differential precession, due to
pressure, destroying the eccentric disc as seen in Chan et al. (2022).

This paper is structured as follows. In Section 2 we give an overview of
eccentric disc geometry and orbital coordinate systems. In Section 3 we extent
the Hamiltonian formalism of Ogilvie & Lynch (2019) to unstratified ideal MHD
discs, which we use to derive linear theory in Section 4 and compute nonlinear
eccentric modes in Section 5. We compare the eccentric disc theory against 2D
MHD simulations in Section 6. Finally, we present our conclusions in Section 8
and a derivation of the magnetic vector potential is given in the appendix to
aid with future numerical work.

##  2 Eccentric Disc Geometry

The geometry of an eccentric disc consists of a set of non-intersecting,
confocal Keplerian ellipses where the dominant fluid motion consist of the
Keplerian motion. These Keplerian orbits slowly evolve due to the effects of
pressure gradients and, in MHD discs, magnetic fields. To describe both the
geometry, and dynamics, of an eccentric discs it is often convenient to make
use of an orbital coordinate system. This is a coordinate system, based on the
orbital elements of celestial mechanics, that describes a point in the mid-
plane of the disc by an orbit labelling coordinate, specifying the orbit the
point lies on, and a coordinate denoting where along that orbit the point lies
Typically such an orbital coordinate system will define a time dependant map
from some circular reference disc onto the physical eccentric disc, with the
dynamics of the eccentric disc being described by the slow evolution of this
orbital coordinate system.

We formulate this orbital coordinate system in terms of a Lagrangian map
between the reference and physical variables
$\boldsymbol{a}\mapsto\boldsymbol{x}$, where $\boldsymbol{a}$ are the orbital
coordinates associated with a fluid element and $\boldsymbol{x}$ is the fluid
element position vector in Cartesian coordinates. This Lagrangian map can be
thought of as mapping some reference circular state into the physical
eccentric disc, similar to Ogilvie (2018) and Ogilvie & Lynch (2019). We
denote the Jacobian associated with this Lagrangian map by
$J_{ij}=\frac{\partial x_{i}}{\partial a^{j}}$ and introduce the notation

$J_{3}=\det(J_{ij})=\frac{J}{J^{\circ}}\frac{H}{H^{\circ}},$ (1)

where $J$ is the Jacobian determinant of the 2D transform and $H$ is the disc
scale height and we adopt the convention that the superscript $\cdot^{\circ}$
denotes a quantity in the reference circular disc. Note that this Jacobian is
for an orbital coordinate system using the stretched vertical coordinate
$\tilde{z}$. In Ogilvie (2018) this is approximated by its value at the
midplane, which is valid when the disc is sufficiently thin. For the
unstratified discs considered here the horizontal and vertical parts of the
transform are separable so this approximation is unnecessary.

We can define an orbital coordinate system, following Ogilvie & Lynch (2019),
where the orbits are labelled by the semimajor axis, a, and the position
around the orbit are labelled by the eccentric anomaly $E$. The shape of each
orbit is controlled by the orbits eccentricity, $e$, and longitude of
pericentre $\varpi$. The $(a,E)$ orbital coordinate system are related to the
cylindrical radius through

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

and to the azimuthal angle, $\phi$, through the true anomaly $f=\phi-\varpi$
which satisfies,

$\cos f=\frac{\cos E-e}{1-e\cos E},\quad\sin f=\frac{\sqrt{1-e^{2}}\sin
E}{1-e\cos E}.$ (3)

We can extend this coordinate system to 3D by taking the disc midplane as a
reference plane and labelling points by their height above/below the midplane,
$z$. One can also introduce a stretched vertical coordinate $\tilde{z}=z/H$,
where $H$ is some characteristic vertical lengthscale such as the disc
thickness or scale height. In unstratified models this is typically taken to
be the sonic length $H=c_{s}\Omega^{-1}$, with $c_{s}$ the sound speed.

It will also often be useful to make use of the mean anomaly $M$, which is
related to the eccentric anomaly through

$M=E-e\sin E=n(t-\tau),$ (4)

where $n=(GM_{1}/a^{3})^{1/2}$ is the mean motion, $M_{1}$ is the mass of the
central object and $\tau$ is the time of pericentre passage. This allows us to
define an $(a,M)$ orbital coordinate system where the position around the
orbit is now denoted by the mean anomaly. This can be extended to 3D in the
same way as the $(a,E)$ coordinates.

The Jacobian determinant of the $(a,M)$ orbital coordinate system can be
expressed as $J=J^{\circ}j$, where $J^{\circ}=a$ and we have introduced the
dimensionless Jacobian determinant from Ogilvie & Lynch (2019),

$\displaystyle\begin{split}j=(J/J^{\circ})&=\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}-\frac{ae_{a}\cos
E}{\sqrt{1-e^{2}}}-ae\varpi_{a}\sin E\\\
&=\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}[1-q\cos(E-E_{0})]\end{split}$ (5)

which is related to the elliptical geometry of the disc and we have introduced
the notation that a subscript $a$ denotes a partial derivative with respect to
the semimajor axis. As in Ogilvie & Lynch (2019) we have introduced the
orbital intersection parameter, q, and $E_{0}$, the eccentric anomaly at which
the maximum orbital compression occurs. These are related to the orbital
elements and their derivatives through

$q\cos E_{0}=\frac{ae_{a}}{1-e(e+ae_{a})},\quad q\sin
E_{0}=\frac{\sqrt{1-e^{2}}ae\varpi_{a}}{1-e(e+ae_{a})}.$ (6)

##  3 Derivation of the Eccentric Disc Hamiltonian

To derive the equation governing the evolution of the eccentric orbits we
shall start from the Lagrangian formulation of ideal MHD. After performing a
vertical integration we exploit a scale separation which occur in “thin” discs
where the Lagrangian can be separated into an $O(1)$ contribution from the
Keplerian terms and $O(\epsilon^{2})$ contributions from the internal and
magnetic energies. In an unstratified disc $\epsilon$ should be thought of as
a characteristic measure of the reciprocal Mach number (or reciprocal Alfvén
number for strongly magnetised discs), rather than the aspect ratio used in
thin disc theory.

The Lagrangian for ideal MHD is (e.g. Ogilvie, 2016)

$L=\iiint\rho^{\circ}\left[\frac{1}{2}u^{2}-\Phi(\boldsymbol{x})-\varepsilon(\boldsymbol{a},J_{3})-\frac{J_{3}B^{2}(\boldsymbol{a},J_{3},J_{ij})}{2\mu_{0}\rho_{0}}\right]d^{3}\boldsymbol{a},$
(7)

where $u$ is the fluid velocity, $\rho^{\circ}$ is the density in the
reference disc, $B^{i}$ is the disc magnetic field, $\varepsilon$ is the
specific internal energy and $\boldsymbol{a}$ are the 3D Lagrangian
coordinates. Ogilvie (2018) developed a fairly general thin disc model based
on affine transforms of fluid elements which provides a convenient setting for
formulating eccentric disc models. The affine transform for a coplanar,
unwarped, disc is

$\boldsymbol{x}=\bar{\boldsymbol{x}}+H\tilde{z}\hat{\boldsymbol{e}}_{z},$ (8)

where $\bar{\boldsymbol{x}}$ are the coordinates in the disc midplane and
$\tilde{z}$ is a stretched vertical coordinate. For an unstratified disc,
where $\Phi(\boldsymbol{x})$ is independent of the vertical coordinate, with
periodic boundary conditions in the vertical direction, we can take
$H=H^{\circ}$ and treat all quantities as independent of $\tilde{z}$.
Therefore we have $J_{3}=J/J^{\circ}$, $J_{zz}=1$, $J_{zx}=J_{zy}=0$ and we
can vertically integrate the Lagrangian (Equation 7) to obtain,

$L=\iint\Sigma^{\circ}\left[\frac{1}{2}u^{2}-\Phi(\bar{\boldsymbol{x}})-\bar{\varepsilon}(\bar{\boldsymbol{a}},J)-\frac{(J/J^{\circ})B^{2}H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\right]d^{2}\bar{\boldsymbol{a}},$
(9)

where $\Sigma^{\circ}=\rho^{\circ}H^{\circ}$ is the surface density of the
reference disc and $\bar{\boldsymbol{a}}$ are the 2D Lagrangian coordinates.
Here overbars denote a 2D (midplane) coordinate system.

To obtain the disc magnetic field we can look for periodic solutions to the
induction equation in an eccentric shearing box (Ogilvie & Barker, 2014).
These periodic magnetic field solutions were derived by Lynch & Ogilvie
(2021), and consists of a vertical field with a quasi-toroidal (orbit
following) field,

$B^{i}=B_{t0}(a,z/H)\frac{J^{\circ}}{J}\frac{H^{\circ}}{H}\frac{v_{\rm
orbital}^{i}}{n}+B_{z0}(a)\frac{J^{\circ}}{J}\hat{e}_{z}^{i}.$ (10)

Here $B_{t0}$ and $B_{z0}$ are the toroidal and vertical magnetic fields in
the reference circular disc (we have chosen not to use the superscript ∘ to
make later expressions less cumbersome). $\boldsymbol{v}_{\rm orbital}$ is the
Keplerian velocity vector. The model of Lynch & Ogilvie (2021) was based on
the eccentric shearing box and Equation 10 was setup to satisfy the
(approximate) solenoidal condition of the local model (Equation C17 of Ogilvie
& Barker (2014)). One can show that it also satisfies the exact solenoidal
condition by adopting the $(a,M,\tilde{z})$ coordinate system and taking the
divergence,

$\nabla_{i}B^{i}=\frac{1}{J_{3}}\frac{\partial}{\partial
M}\left(J_{3}B_{t0}(a,z/H)\frac{J^{\circ}}{J}\frac{H^{\circ}}{H}\right)+\frac{\partial}{\partial\tilde{z}}\left(B_{z0}(a)\frac{J^{\circ}}{J}\right)=0,$
(11)

as a result of the independence of $J_{3}B^{M}$ (first term) and $B^{z}$
(second term) on $M$ and $\tilde{z}$ respectively, and we have used $v_{\rm
orbital}^{i}=n\hat{e}_{M}^{i}$ in this coordinate system.

In cylindrical, unstratified, geometry $H/H^{\circ}=1$ and
$B_{t0}(a,z/H)=B_{t0}(a)$. Making use of a dimensionless Jacobian determinant
$j=J/J^{\circ}$ we obtain the following for the magnetic pressure due to the
magnetic field given by Equation 10,

$\frac{B^{2}}{2\mu_{0}}=a^{2}\frac{B_{t0}^{2}(a)}{2\mu_{0}}j^{-2}\frac{1+e\cos
E}{1-e\cos E}+\frac{B_{z0}^{2}}{2\mu_{0}}(a)j^{-2},$ (12)

where we have made use of $v_{\rm orbital}^{2}=a^{2}n^{2}\frac{1+e\cos
E}{1-e\cos E}$.

Substituting Equation 12 for the magnetic pressure into the Lagrangian we
arrive at

$\displaystyle\begin{split}L&=\iint\Sigma^{\circ}\Biggl{[}\frac{1}{2}u^{2}-\Phi(\bar{\boldsymbol{x}})-\bar{\varepsilon}(\bar{\boldsymbol{a}},J)\\\
&-a^{2}\frac{1+e\cos E}{1-e\cos
E}\frac{(J/J^{\circ})^{-1}B^{2}_{t0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}-\frac{(J/J^{\circ})^{-1}B^{2}_{z0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\Biggr{]}d^{2}\bar{\boldsymbol{a}}.\end{split}$
(13)

Expanding Equation 13 into the Keplerian Lagrangian $L_{K}$ and a
perturbation,

$\displaystyle\begin{split}L&=L_{k}-\iint\Sigma^{\circ}\Biggl{[}V(\bar{\boldsymbol{x}})+\bar{\varepsilon}(\bar{\boldsymbol{a}},J)\\\
&+a^{2}\frac{1+e\cos E}{1-e\cos
E}\frac{(J/J^{\circ})^{-1}B^{2}_{t0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}+\frac{(J/J^{\circ})^{-1}B^{2}_{z0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\Biggr{]}d^{2}\bar{\boldsymbol{a}},\end{split}$
(14)

where $V(\boldsymbol{x})=\Phi(\boldsymbol{x})-\Phi_{K}(\boldsymbol{x})$. At
leading order we have,

$\frac{\delta L}{\delta\boldsymbol{x}}\approx\frac{\delta
L_{K}}{\delta\boldsymbol{x}}=0,$ (15)

which is Keplerian orbital motion in the plane. Performing the Whitham/orbit
average (Whitham, 1965) of Equation 14, using the fact that the Lagrangian is
nearly integrable,

$\displaystyle\begin{split}L&=\int
m_{a}\left[na^{2}\dot{M}+na^{2}\sqrt{1-e^{2}}\dot{\varpi}+\frac{GM_{1}}{2a}\right]\,da\\\
&-\int m_{a}\Biggl{[}\langle
V(\bar{\boldsymbol{x}})\rangle+\langle\bar{\varepsilon}(\bar{\boldsymbol{a}},J)\rangle\\\
&+a^{2}\frac{B^{2}_{t0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\left\langle\frac{1+e\cos
E}{1-e\cos
E}j^{-1}\right\rangle+\frac{B^{2}_{z0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\langle
j^{-1}\rangle\Biggr{]}\,da,\end{split}$ (16)

where $\langle\cdot\rangle=\frac{1}{2\pi}\int\cdot\,dM$ denotes an orbit
average and we have introduced the mass per unit semimajor axis, $m_{a}$,
which is related to the surface density in the reference circular disc by

$m_{a}=2\pi a\Sigma^{\circ}.$ (17)

The dynamics of $M$ is dominated by orbital motion. Therefore we separate out
the terms in the Lagrangian describing the orbital motion from those
describing the slow evolution of the disc orbits,

$\displaystyle\begin{split}L&=\int
m_{a}na^{2}\sqrt{1-e^{2}}\dot{\varpi}\,da-\int m_{a}\Biggl{[}\langle
V(\bar{\boldsymbol{x}})\rangle+\langle\bar{\varepsilon}(\bar{\boldsymbol{a}},J)\rangle\\\
&+a^{2}\frac{B^{2}_{t0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\left\langle\frac{1+e\cos
E}{1-e\cos
E}j^{-1}\right\rangle+\frac{B^{2}_{z0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\langle
j^{-1}\rangle\Biggr{]}\,da.\end{split}$ (18)

The associated Hamiltonian is obtained via a Legendre transform,

$\displaystyle\begin{split}\mathcal{H}&=\int
m_{a}na^{2}\sqrt{1-e^{2}}\dot{\varpi}\,da-L\\\ &=\int m_{a}\Biggl{[}\langle
V(\bar{\boldsymbol{x}})\rangle+\langle\bar{\varepsilon}(\bar{\boldsymbol{a}},J)\rangle\\\
&+a^{2}\frac{B^{2}_{t0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\left\langle\frac{1+e\cos
E}{1-e\cos
E}j^{-1}\right\rangle+\frac{B^{2}_{z0}(\bar{\boldsymbol{a}})H^{\circ}}{2\mu_{0}\Sigma^{\circ}}\langle
j^{-1}\rangle\Biggr{]}\,da.\end{split}$ (19)

Henceforth we shall only consider Keplerian potentials so that
$V(\bar{\boldsymbol{x}})=0$. For a perfect gas we can write (Ogilvie & Lynch,
2019),

$m_{a}\langle\bar{\varepsilon}(\bar{\boldsymbol{a}},J)\rangle=H_{a}^{\circ}F^{(\gamma)},$
(20)

where we have introduced the geometric part of the Hamiltonian,

$F^{(p)}=\frac{1}{p-1}\langle j^{-(p-1)}\rangle,$ (21)

along with the circular Hamiltonian density,

$H_{a}^{\circ}=2\pi aP_{g}^{\circ},$ (22)

where $P_{g}^{\circ}$ is the vertically integrated gas pressure in the
reference disc. Introducing the Alfvén velocity in the reference disc:
$v^{i}_{a0}=B^{i}_{0}/\sqrt{\mu_{0}\rho^{\circ}}$, we parameterise the
magnetic field strength in terms of a dimensionless toroidal $(V_{t})$ and
vertical $(V_{z})$ Alfvén velocities, where $V_{t}=av^{E}_{a0}/c_{s}$ and
$V_{z}=v^{z}_{a0}/c_{s}$.

From Ogilvie & Lynch (2019) we have the following expression for $F^{(2)}$,

$F^{(2)}(e,q,E_{0})=\langle
j^{-1}\rangle=\frac{\sqrt{1-e^{2}}}{1-e(e+ae_{a})}\frac{q-e(1-\sqrt{1-q^{2}})\cos
E_{0}}{q\sqrt{1-q^{2}}}.$ (23)

Similarly we obtain can obtain an expression for $\left\langle\frac{1+e\cos
E}{1-e\cos E}j^{-1}\right\rangle$,

$\displaystyle\begin{split}\left\langle\frac{1+e\cos E}{1-e\cos
E}j^{-1}\right\rangle&=\frac{\sqrt{1-e^{2}}}{1-e(e+ae_{a})}\frac{q+e(1-\sqrt{1-q^{2}})\cos
E_{0}}{q\sqrt{1-q^{2}}}\\\ &=F^{(2)}(e,q,E_{0}+\pi).\end{split}$ (24)

Thus the vertical magnetic field acts like a $\gamma=2$ perfect gas, while the
quasi-toroidal magnetic field acts like a $\gamma=2$ gas with an anti-phased
orbital compression.

We can then write the Hamiltonian as

$\displaystyle\begin{split}\mathcal{H}=\int
H_{a}^{\circ}\Biggl{(}F^{(\gamma)}(e,q,E_{0})&+\frac{\gamma}{2}V_{t}^{2}F^{(2)}(e,q,E_{0}+\pi)\\\
&+\frac{\gamma}{2}V_{z}^{2}F^{(2)}(e,q,E_{0})\Biggr{)}\,da,\end{split}$ (25)

where the factor of $\gamma$ in the magnetic terms appear as a result of
factoring out $H^{\circ}_{a}$. This Hamiltonian can be split into
contributions from the gas internal energy, and the energy in the quasi-
toroidal and vertical magnetic fields, with $\mathcal{H}=\mathcal{H}_{\rm
gas}+\mathcal{H}_{\rm tor}+\mathcal{H}_{\rm vert}$, corresponding to the first
second and third term in the brackets of Equation 25. As in gas discs each of
these terms are a product of the, geometry independent, Hamiltonian density in
the reference disc and a geometric part which encapsulates the dependence on
the orbital geometry.

It is convenient to reformulate this Hamiltonian as a Hamiltonian for a single
effective gas. We can do this by introducing a new geometric part of the
Hamiltonian,

$\displaystyle\begin{split}&F_{V_{t},V_{z}}(e,q,E_{0})=\frac{1}{1+\frac{\gamma}{2}(V_{t}^{2}+V_{z}^{2})}F^{(\gamma)}(e,q,E_{0})\\\
&+\frac{\frac{\gamma}{2}V_{t}^{2}}{1+\frac{\gamma}{2}(V_{t}^{2}+V_{z}^{2})}F^{(2)}(e,q,E_{0}+\pi)+\frac{\frac{\gamma}{2}V_{z}^{2}}{1+\frac{\gamma}{2}(V_{t}^{2}+V_{z}^{2})}F^{(2)}(e,q,E_{0}),\end{split}$
(26)

which is a weighted sum of two adiabatic gas $F^{(p)}$ with different ratio of
specific heats and a third nonadiabatic gas term for the toroidal field. In
this case the Hamiltonian can be written as,

$\mathcal{H}=\int\tilde{H}_{a}^{\circ}F_{V_{t},V_{z}}(e,q,E_{0})da,$ (27)

where we have introduced

$\tilde{H}_{a}^{\circ}=2\pi
a(P^{\circ}_{g}+P_{m}^{\circ})=\left(1+\frac{\gamma}{2}V_{t}^{2}+\frac{\gamma}{2}V_{z}^{2}\right)H_{a}^{\circ}=\frac{1+\beta^{\circ}}{\beta^{\circ}}H_{a}^{\circ},$
(28)

which is the Hamiltonian density in the reference circular disc and
$\beta^{\circ}$, the plasma-$\beta$ in the reference circular disc. Unlike the
simpler perfect gas case $F_{V_{t},V_{z}}$ is no longer only a function of the
geometry for a given ratio of specific heat, it also depends on the “partial
pressure” of the constitutive effective gasses on a given orbit.

Hamilton’s equations in the noncannoical $e(a,t)$, $\varpi(a,t)$ are (Ogilvie
& Lynch, 2019)

$\displaystyle m_{a}\frac{\partial e}{\partial t}$
$\displaystyle=\frac{\sqrt{1-e^{2}}}{na^{2}e}\frac{\delta\mathcal{H}}{\delta\varpi},$
(29) $\displaystyle m_{a}\frac{\partial\varpi}{\partial t}$
$\displaystyle=-\frac{\sqrt{1-e^{2}}}{na^{2}e}\frac{\delta\mathcal{H}}{\delta
e}.$ (30)

The ideal MHD eccentric disc Hamiltonian preserves the symmetries of the
unmagnetised eccentric disc Hamiltonian of Ogilvie & Lynch (2019); i.e. time
translation and global rotation, with the Hamiltonian only depending on
$\varpi$ through it derivative $\varpi_{a}$. As such, following Ogilvie &
Lynch (2019), one can show that the total Hamiltonian, $\mathcal{H}$, which in
an unstratified disc corresponds to the sum of the magnetic and internal
energies, and the angular momentum deficit (AMD), a positive definite measure
of the total eccentricity commonly used in celestial mechanics,

$\mathcal{C}=\int m_{a}na^{2}\left(1-\sqrt{1-e^{2}}\right)da,$ (31)

are conserved.

The simplest solutions to the eccentric disc equations are the eccentric
modes, which are solutions where $e$ is independent of time and the disc is
untwisted and uniformly precessing at an angular frequency $\omega$. These
solutions are a particularly convenient setting for numerical simulations as
they avoid strong differential precession seen in generic eccentricity
profiles. The eccentric mode equation is obtained from Equation 30,

$m_{a}\omega=-\frac{\sqrt{1-e^{2}}}{na^{2}e}\frac{\delta\mathcal{H}}{\delta
e}.$ (32)

Equation 29 is automatically satisfied as the disc is untwisted (so
$\frac{\delta\mathcal{H}}{\delta\varpi}=0$) and $e$ is independent of time.

As shown in Ogilvie & Lynch (2019), Equation 32 can be written as

$\omega\frac{\delta L}{\delta e}=\frac{\delta\mathcal{H}}{\delta e},$ (33)

where

$L=\int m_{a}na^{2}\sqrt{1-e^{2}}\,da$ (34)

is the angular momentum. Equation 33 can be interpreted as a variational
problem which makes $\mathcal{H}$ (here corresponding to the total disc
internal + magnetic energy) stationary at a fixed angular momentum
$L=\mathrm{const}$.

For constant $V_{t},V_{z}$ the eccentric mode equations are explicitly
(Ogilvie & Lynch, 2019):

$\displaystyle\begin{split}-\frac{\omega
m_{a}}{\tilde{H}_{a}^{\circ}}\frac{na^{2}e}{\sqrt{1-e^{2}}}&=\frac{\partial
F_{V_{t},V_{z}}}{\partial e}-ae_{a}\frac{\partial^{2}F_{V_{t},V_{z}}}{\partial
e\partial f}\\\ &-a(2e_{a}+ae_{aa})\frac{\partial^{2}F_{V_{t},V_{z}}}{\partial
f^{2}}-\frac{d\ln(\tilde{H}^{\circ}_{a})}{d\ln a}\frac{\partial
F_{V_{t},V_{z}}}{\partial f},\end{split}$ (35)

where we have introduced $f=e+ae_{a}$. If $V_{t}$ or $V_{z}$ depends on the
semimajor axis then the equation for an eccentric mode becomes,

$\displaystyle\begin{split}-\frac{\omega
m_{a}}{\tilde{H}_{a}^{\circ}}\frac{na^{2}e}{\sqrt{1-e^{2}}}&=\frac{\partial
F_{V_{t},V_{z}}}{\partial e}-ae_{a}\frac{\partial^{2}F_{V_{t},V_{z}}}{\partial
e\partial f}\\\ &-a(2e_{a}+ae_{aa})\frac{\partial^{2}F_{V_{t},V_{z}}}{\partial
f^{2}}-\frac{d\ln(H^{\circ}_{a})}{d\ln a}\frac{\partial
F_{V_{t},V_{z}}}{\partial f}\\\
&-\frac{\gamma}{2}\frac{\beta^{\circ}}{1+\beta^{\circ}}\left(\frac{dV_{t}^{2}}{d\ln
a}\frac{\partial F^{(2)}}{\partial
f}\Biggl{|}_{E_{0}=\pi}+\,\frac{dV_{z}^{2}}{d\ln a}\frac{\partial
F^{(2)}}{\partial f}\Biggl{|}_{E_{0}=0}\right).\end{split}$ (36)

Note that we have the perfect gas circular Hamiltonian density
($H_{a}^{\circ}$) on the forth term on the right hand side. We also have
$F^{(2)}|_{E_{0}=0}=F^{(2)}(e,q(e,f),0)$ and
$F^{(2)}|_{E_{0}=\pi}=F^{(2)}(e,q(e,f),\pi)$.

For untwisted discs $F_{V_{t},V_{z}}(e,f)$ has an apparent singularity when
$e=f$, where the eccentricity gradients vanish. Following Ogilvie & Lynch
(2019) this apparent singularity can be removed using the trigonometric
parametrisation $e=\sin 2\alpha$, $f=\sin 2\beta$. Expressions for
$F^{(1)}(e,q,0)$ and $F^{(2)}(e,q,0)$, in terms of this parametrisation are
given in Appendix C of Ogilvie & Lynch (2019). For including quasi-toroidal
fields we will also need

$F^{(2)}(e,q,\pi)=\cos(2\alpha)\cos(\alpha-\beta)\sec(2\beta)\sec(\alpha+\beta).$
(37)

##  4 Linear Theory

When $e$, $ae_{a}$ and $ae\varpi_{a}$ are much less than unity, the geometric
part of the Hamiltonian density in a 2D disc can be approximated as (Ogilvie &
Lynch, 2019)

$F^{(\rm
2D)}\approx\frac{1}{2}e(e+ae_{a})+\frac{1}{4}\gamma\left[(ae_{a})^{2}+(ae\varpi_{a})^{2}\right],$
(38)

where we have dropped an unimportant constant term, that has no influence on
the dynamics, so that we can use this expression for isothermal discs. In
addition to Equation 38 we require the linear limit of
$F^{(2)}(e,q,E_{0}+\pi)$, to include the quasi-toroidal field, this can be
obtained in a similar way and is

$F^{(2)}(e,q,E_{0}+\pi)\approx\frac{1}{2}e(e+ae_{a})+\frac{1}{2}\left[(ae_{a})^{2}+(ae\omega_{a})^{2}\right]+aee_{a},$
(39)

where the first two terms arise from the adiabatic variation of the magnetic
pressure, similar to the vertical field, while the last term arises from the
magnetic tension and the non-adiabatic variation of the magnetic pressure.
Combining these we arrive at an expression for $F_{V_{t},V_{z}}$,

$\displaystyle\begin{split}F_{V_{t},V_{z}}(e,q,E_{0})&\approx\frac{1}{2}e(e+ae_{a})+\frac{\tilde{\gamma}}{4}\left[(ae_{a})^{2}+(ae\omega_{a})^{2}\right]\\\
&+\frac{\frac{\gamma}{2}V_{t}^{2}}{1+\frac{\gamma}{2}(V_{t}^{2}+V_{z}^{2})}aee_{a},\end{split}$
(40)

where we have introduced a modified ratio of specific heats

$\tilde{\gamma}=\gamma\frac{1+V_{t}^{2}+V_{z}^{2}}{1+\frac{\gamma}{2}(V_{t}^{2}+V_{z}^{2})}.$
(41)

To connect with the existing work on linear eccentric disc theory it is useful
to rewrite Equation 40 in terms of the complex eccentricity
$\mathcal{E}=e\exp(i\varpi)$, the geometric part of the Hamiltonian in the
linear limit is then

$\displaystyle\begin{split}F_{V_{t},V_{z}}(e,q,E_{0})&\approx\frac{1}{2}[|\mathcal{E}|^{2}+\mathrm{Re}(a\mathcal{E}\mathcal{E}^{*}_{a})]+\frac{\tilde{\gamma}}{4}|a\mathcal{E}_{a}|^{2}\\\
&+\frac{\frac{\gamma}{2}V_{t}^{2}}{1+\frac{\gamma}{2}(V_{t}^{2}+V_{z}^{2})}\mathrm{Re}(a\mathcal{E}\mathcal{E}^{*}_{a}).\end{split}$
(42)

The, non-canonical, Hamilton’s equations for the complex eccentricity
are111This form of Hamilton’s equations for the eccentric disc theory was
originally suggested by Prof. Gordon Ogilvie in an earlier draft of Ogilvie &
Lynch (2019) as a way of connecting Hamiltonian eccentric disc theory with the
Shrödinger equation.

$m_{a}\dot{\mathcal{E}}=-\frac{2i\sqrt{1-e^{2}}}{na^{2}}\frac{\delta\mathcal{H}}{\delta\mathcal{E}^{*}},$
(43)

with the functional derivative of $\mathcal{E}^{*}$ being related to the
functional derivative of $e$ and $\varpi$ through,

$\frac{\delta}{\delta\mathcal{E}^{*}}=\frac{1}{2}\mathrm{e}^{i\varpi}\left(\frac{\delta}{\delta
e}+\frac{i}{e}\frac{\delta}{\delta\varpi}\right).$ (44)

Substituting in the Linear form of the Hamiltonian into Equation 43 and
performing the functional derivative we obtain a linear equation for the
evolution of the complex eccentricity in an unstratified ideal MHD disc,

$\displaystyle\begin{split}2\Sigma^{\circ}na^{3}\frac{\partial\mathcal{E}}{\partial
t}&=\frac{\partial}{\partial
a}\left(i\tilde{\gamma}P^{\circ}a^{3}\frac{\partial\mathcal{E}}{\partial
a}\right)+ia^{2}\frac{dP^{\circ}}{da}\mathcal{E}+i\partial_{a}\left[a^{2}\frac{(aB_{0}^{\phi})^{2}}{\mu_{0}}\right]\mathcal{E}\end{split}$
(45)

where we have used $M_{a}=2\pi a\Sigma^{\circ}$ and
$\tilde{H}_{a}^{\circ}=2\pi aP^{\circ}$. The first two terms on the right-hand
side correspond to an adiabatic gas with an effective ratio of specific heats,
$\tilde{\gamma}$ set by the plasma-$\beta$. The final term arises from the
non-adiabatic change to the magnetic pressure from the stretching of the
magnetic field lines222This can be shown by deriving the linear, magnetised,
eccentric disc equations following a similar procedure to Goodchild & Ogilvie
(2006), a task that is significantly more involved than taking the linear
limit of the Hamiltonian theory.. Terms arising due to the magnetic tension
cancel with additional non-adiabatic terms along with the modification to the
background rotation profile as a result of the magnetic tension in the
unperturbed disc.

Specialising to an eccentric mode in a disc with a purely vertical field
($V_{t}=0$) the equation simplifies to

$2\Sigma^{\circ}na^{3}\omega\mathcal{E}=\frac{\partial}{\partial
a}\left(\tilde{\gamma}P^{\circ}a^{3}\frac{\partial\mathcal{E}}{\partial
a}\right)+\frac{dP^{\circ}}{da}a^{2}\mathcal{E}.$ (46)

When the gas pressure, $P_{g}$, and dimensionless Alfvén velocity, $V_{z}$,
are constant one can rewrite the above equation as

$2\Sigma^{\circ}na^{3}\tilde{\omega}\mathcal{E}=\frac{\partial}{\partial
a}\left(\gamma P_{g}^{\circ}a^{3}\frac{\partial\mathcal{E}}{\partial
a}\right),$ (47)

where we have introduced a rescaled precession frequency
$\tilde{\omega}=\omega/(1+V_{z}^{2})$. Therefore, under these restrictions on
the pressure and magnetic field, the eccentric mode in the magnetised and
unmagnetised discs are identical and differ only by their precession
frequency.

##  5 Nonlinear Modes in an Isothermal Disc

Our primary motivation for extending the unstratified eccentric disc theory to
include magnetic fields is to provide initial conditions for simulations of
the eccentric MRI. As such we focus on calculating the nonlinear eccentric
modes found in a disc contained between two, circular, rigid walls as done to
setup the hydrodynamical simulations of Barker & Ogilvie (2016), rather than
the more realistic free boundaries considered in Ogilvie & Lynch (2019). As a
model of a realistic MHD disc, however, the unstratified model derived in the
previous section has major limitations; namely it fails to account for the
dynamical vertical structure of the disc which is known to be important to
correctly describe the dynamics of eccentric discs Ogilvie (2001, 2008);
Ogilvie & Barker (2014); Ogilvie & Lynch (2019) and can significantly increase
the field strength of the quasi-toroidal field, particularly for more
nonlinear eccentric discs (Lynch & Ogilvie, 2021). There is also the issue of
how the disc interacts with the external magnetic fields.

We, thus, consider a simple MHD generalisation of the mode computed in Barker
& Ogilvie (2016). This consists of a globally isothermal disc with a constant
reference surface density, $\Sigma^{\circ}$, with the Hamiltonian density in
the reference circular disc being $H^{\circ}_{a}=2\pi
ac_{s}^{2}\Sigma^{\circ}$, where $c_{s}$ is a constant sound speed.

For the modes that will be simulated in Section 6 we impose a purely vertical
field with,

$V_{t}=0,\quad V_{z}=\frac{nl_{\rm MRI}}{2\pi c_{s}\sqrt{16/15}}W(a).$ (48)

For comparison, in this section, we also compute a quasi-toroidal case with

$V_{t}=\frac{nl_{\rm MRI}}{2\pi c_{s}\sqrt{16/15}}W(a),\quad V_{z}=0.$ (49)

In both cases $l_{\rm MRI}$ is a constant lengthscale, which in the vertical
field case corresponds to the lengthscale of the fastest growing MRI mode in
the reference circular disc. $W(a)$ describes the taper on the inner and outer
disc boundaries, for which we use

$W(a)=\frac{1}{2}\left[1+\tanh\left(\frac{a-a_{\rm in}}{w_{\rm
transition}}\right)\right]\left[1-\tanh\left(\frac{a-a_{\rm out}}{w_{\rm
transition}}\right)\right],$ (50)

when we wish to include a taper. The disc is contained within two rigid
circular walls located at $a_{\rm min}$ and $a_{\rm max}$, such that $a_{\rm
min}\leq a_{\rm in}\leq a_{\rm out}\leq a_{\rm max}$. We therefore have
boundary conditions $e(a_{\rm min})=e(a_{\rm max})=0$. The precession
frequency of the mode, $\omega$, is an eigenvalue of the problem. One can
solve Equation 36for the eccentric mode by specifying $e_{a}$ on the inner
boundary and employing a shooting method to obtain $\omega$.

To provide initial conditions for our simulations we solve for eccentric modes
with $l_{\rm MRI}=c_{s}/n(a_{\rm min})$, $a_{\rm min}=1$, $a_{\rm in}=1.5$,
$a_{\rm out}=4.5$, $a_{\rm max}=5$ and $c_{s}=0.05$, with purely vertical
fields. These modes are shown in Figure 1. The simulated modes with different
choices of parameters (as discussed in Section 6) have functionally
indistinguishable eccentricity profiles.

Figure 1: Eccentricity profiles of the eccentric modes used in simulations
with $\max(e)=0.2$, $\max(e)=0.35$ and $\max(e)=0.5$ respectively. For both
modes $l_{\rm MRI}=c_{s}/n(a_{\rm min})$, $a_{\rm min}=1$, $a_{\rm in}=1.5$ ,
$a_{\rm out}=4.5$, $a_{\rm max}=5$ and $w_{\rm transition}=0.05$. The
eccentric modes in the unmagnetised and MHD discs with a purely vertical field
are almost indistinguishable. The dotted line is the “limiting slope"
solution.

The MHD eccentric modes depicted in Figure 1 are nearly indistinguishable from
the unmagnetised case. The same is true of eccentric modes computed using a
similar strength quasi-toroidal field. This is perhaps not surprising given
the magnetic field strength in these discs is set by the requirement that the
circular reference disc is MRI unstable. This means the plasma-$\beta$ in the
reference disc never drops below $\beta=84$. In the eccentric disc the range
of plasma-$\beta$ attained is greater as a result of the lateral orbital
compression which occurs in the presence of eccentricity gradients.

Figures 2 and 3 show the minimum and maximum plasma-$\beta$ on an orbit, in
the absence of a taper, for the vertical and quasi-toroidal field models
respectively. In the absence of a taper the orbital compression results in
regions of high magnetic field strengths in the inner disc, particularly for
the mode with $\max[e]=0.5$ . This effect is lessened when the taper is
included as the magnetic field strength drops to zero close to the boundary
where the effects of orbital compression are greatest. Despite attaining
plasma-$\beta$ as low as $\beta\sim 4$ for the $\max[e]=0.5$ mode, without
taper, differs only slightly from the unmagnetised case (having a slightly
lower eccentricity gradient on the inner boundary). This is, in part, a
geometric effect where the shape of highly nonlinear eccentric mode is
dictated by the requirement that $|q|<1$ to avoid an orbital intersection. As
discussed in Barker & Ogilvie (2016); Ogilvie & Lynch (2019) this results in a
limiting slope solution given by

$ae=\begin{cases}a-a_{\rm min},&a_{\rm min}<a<\bar{a}\\\ -a+a_{\rm
max},&\bar{a}<a<a_{\rm max}\end{cases}$ (51)

where $\bar{a}=(a_{\rm min}+a_{\rm max})/2$. In principle higher order
limiting slope solutions might depend on the magnetic field strength as it is
not obvious how the mode selects nodes for the higher order modes.

Figure 2: Minimum (dashed lines) and maximum (solid lines) plasma beta on
orbits for the eccentric modes with a purely vertical field and no taper. The
dotted line is the plasma beta in the circular disc. Orbital compression can
lead to significant local enhancement of the magnetic fields, while at the
eccentricity maxima the plasma-$\beta$ is constant around the orbit. The
simulations do not reach as extreme plasma-$\beta$s as a result of the taper.
Figure 3: Same as Figure 2 but for a quasi-toroidal field. There is a greater
variation of the magnetic field around the orbit in the inner disc. Unlike the
vertical field case the stretching of the field lines around an eccentric
orbit means the magnetic field always varies around the orbit, even at the
eccentricity maxima where the surface density is constant.

The large orbital compression responsible for the regions of greatly enhanced
magnetic fields in the modes calculated above are primarily a consequence of
the adoption of rigid circular boundaries. As shown in Ogilvie & Lynch (2019)
adoption of more realistic free boundaries conditions result in more moderate
eccentricity gradients for a given $\max[e]$. In Appendix A we solve for the
MHD eccentric modes with free boundaries and a taper in both the magnetic
field and the surface density. The variation of the magnetic field strength
and the orbit are much reduced compared with the equivalent rigid boundary
eccentric mode due to the smaller eccentricity gradients. This confirms that
the strong enhancement of the magnetic fields seen in the modes computed for
Figures 2 and 3 are primarily a consequence of the rigid wall boundaries.

While the strong enhancement of the magnetic field in the fundamental (zero
node) modes considered thus far are primarily a consequence of the choice of
boundary conditions, higher order modes (i.e. with multiple nodes) will attain
larger eccentricity gradients for a given $\max[e]$. Figure 4 shows the
maximum magnetic field enhancement for eccentric modes with $\max[e]=0.1$.
These modes are discrete due to the boundary conditions. Increasing mode
number results in an increased $\max[q]$, resulting in an increasing magnetic
field enhancement due to the greater lateral compression. This magnetic field
enhancement approximately follows $1/(1-q_{\rm max})$ for small $\max[e]$.
This magnetic field enhancement may have important consequence for the
eccentric MRI if the higher field strengths are able to stabilise the MRI.

Figure 4: Mode number dependence of the magnetic field enhancement for modes
with $\max[e]=0.1$. Dashed line shows the function $1/(1-q_{\rm max})$, which
is the approximate magnetic field enhancement for disc with low eccentricity.

##  6 Nonlinear simulations

We have run nonlinear MHD simulations to demonstrate the integrity of the
eccentric disc solutions calculated in Section 5, specifically the vertical
field case including a taper. We limit our focus to purely 2D simulations in
this paper in order to isolate the eccentric modes’ role as MHD equilibria,
which are of course unstable in three dimensions. In a companion paper, we
explore the growth and turbulent saturation of the instabilities of these
equilibria in fully 3D simulations.

Table 1: Details of the nonlinear simulations considered in this paper. All simulations were run on a uniform cylindrical mesh with $r_{1}/r_{0}=5$, and $c_{s}=0.05a_{\text{min}}n(a_{\text{min}}).$ The magnetized runs include purely vertical fields as described by (48), with $l_{\text{MRI}}=c_{s}/n(a_{\text{min}})$, $w_{\text{transition}}=0.1a_{\text{min}},$ $a_{\text{in}}=2a_{\text{min}}$, and $a_{\text{out}}=4a_{\text{min}}$ (except for m35t, which uses $w_{\text{transition}}=0.05a_{\text{min}}$, $a_{\text{in}}=1.5a_{\text{min}}$, and $a_{\text{out}}=4.5a_{\text{min}}$). The final three columns compare precession frequencies predicted by the eccentric mode calculations with frequencies and decay rates measured in the simulations. Simulation | $\max[e]$ |  $B_{z}$? | $N_{R}\times N_{\phi}$ | Predicted $\omega$ | Simulation $Re[\omega]$ | Simulation $Im[\omega]$
---|---|---|---|---|---|---
h2 | $0.20$ | No | $480\times 800$ | -0.004235 | -0.0044 | $-4.8\times 10^{-5}$
h35 | $0.35$ | No | $480\times 800$ | -0.005012 | -0.0054 | $-1.0\times 10^{-4}$
h5 | $0.50$ | No | $480\times 800$ | -0.007459 | -0.0045 | $-5.5\times 10^{-3}$
m2 | $0.20$ | Yes | $480\times 800$ | -0.004244 | -0.0044 | $-4.8\times 10^{-5}$
m35 | $0.35$ | Yes | $480\times 800$ | -0.005042 | -0.0054 | $-1.0\times 10^{-4}$
m35t | $0.35$ | Yes | $480\times 800$ | -0.005047 | -0.0055 | $-9.9\times 10^{-5}$
m35l | $0.35$ | Yes | $240\times 400$ | -0.005042 | -0.0033 | $-3.8\times 10^{-4}$
m35h | $0.35$ | Yes | $960\times 800$ | -0.005042 | -0.0053 | $-6.2\times 10^{-5}$
m5 | $0.50$ | Yes | $480\times 800$ | -0.007582 | -0.0045 | $-5.5\times 10^{-3}$
m5h | $0.50$ | Yes | $960\times 800$ | -0.007582 | -0.0046 | $-5.6\times 10^{-3}$

###  6.1 Setup

We use a uniform grid version of the code RAMSES (Teyssier, 2002; Fromang et
al., 2006; Faure et al., 2014) 333Available at
https://sourcesup.renater.fr/projects/dumses/ , which employs a high-order
Godunov method to solve the magnetohydrodynamic equations under the
cylindrical approximation (i.e., without vertical gravity). Taken with a
purely isothermal equation of state, these are

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho{\boldsymbol{u}})=0,$ (52)
$\frac{\partial(\rho{\boldsymbol{u}})}{\partial
t}+\nabla\cdot(\rho{\boldsymbol{u}\boldsymbol{u}}-{\boldsymbol{B}\boldsymbol{B}})+\nabla\left(P+\frac{\boldsymbol{B}\cdot\boldsymbol{B}}{2}\right)=-\rho\nabla\Phi,$
(53) $\frac{\partial{\boldsymbol{B}}}{\partial
t}+\nabla\cdot({\boldsymbol{u}\boldsymbol{B}}-{\boldsymbol{B}\boldsymbol{u}})=0,$
(54)

where $\rho$ is the gas density, $P=c_{s}^{2}\rho$ is the pressure, $c_{s}$ is
the sound speed, ${\bf u}$ is the fluid velocity field, $\Phi=-GM_{1}/R$ is
the (Newtonian and cylindrical) gravitational potential of a central mass $M$,
and ${\bf B}$ is the magnetic field. We initialise purely two-dimensional
simulations with the surface densities, radial and azimuthal velocities, and
purely vertical magnetic fields (from Equations 60-63) corresponding to
eccentricity profiles like those shown in Fig. 1. Table 1 lists relevant
properties. We do not simulate the quasi-toroidal field case in this paper.
This case is complicated by the difficulty in ensuring the solenoidal
condition is satisfied when switching from the orbital to the polar grid, this
is best done by use of a vector potential (derived in Appendix B.2) which is
not implemented in the version of RAMSES we are using.

We impose quasi-rigid wall boundary conditions at both the inner and outer
radial boundaries $r_{0}=a_{\text{min}}$ and $r_{1}=a_{\text{max}}$, fixing
$u_{R}=0$ and setting $u_{\phi}$ by the Keplerian angular velocity of the
circular reference disk. Nonzero eccentricity gradients at the boundaries
imply nonzero surface density variations with $\phi$. We therefore use a zero-
gradient boundary condition for the density, setting its value in the ghost
cells to the value of the last cell in the active domain. We lastly set
$B_{z}$ to zero in the ghost cells ($B_{R}$ and $B_{\phi}$ remain identically
zero throughout, and so the magnetic field remains trivially solenoidal in
these 2D simulations).

To track the evolution of eccentricity in our simulations, we compute the
semimajor axis of each grid cell from (e.g., Miranda et al., 2017)

$a(R,\phi)=\left(\frac{2}{R}-\frac{u^{2}}{GM_{1}}\right)^{-1},$ (55)

and the eccentricity vector from

${\boldsymbol{e}}=[e\cos\varpi,e\sin\varpi]=(GM_{1})^{-1}\left[u^{2}{\boldsymbol{R}}-({\boldsymbol{u}\cdot
R}){\boldsymbol{u}}\right]-\hat{\boldsymbol{R}}.$ (56)

Binning the eccentricity values in every cell by semimajor axis, we average
within each bin to produce one-dimensional eccentricity profiles
$\tilde{e}=\tilde{e}(a)$ at each timestep. These we use in turn to compute the
integrated angular momentum deficit $\mathcal{C}$. We additionally consider
the time-evolution of the total Hamiltonian

$\mathcal{H}=\int_{a_{\text{min}}}^{a_{\text{max}}}\\!\int_{0}^{2\pi}(c_{s}^{2}\Sigma\ln\Sigma+P_{m})\,R\,d\phi\,dR,$
(57)

which is conserved in the ideal secular theory.

###  6.2 Results

Figure 5: Spacetime diagrams showing radial profiles of radial velocity at a
fixed $\phi$ vs. time for both hydrodynamic (names starting with $h$) and MHD
($m$) simulations. The periodic changes in the sign of $v_{r}$ over
$t\Omega_{0}\sim 1000$ illustrate the retrograde precession of the eccentric
distortions. Figure 6: Same as Fig. 5, but for vertical magnetic field in the
MHD simulations.

For a given eccentricity profile and resolution, our hydrodynamic and MHD
simulations show remarkably similar evolution. Table 1 compares the precession
frequencies we observe in simulations against the eigenvalues computed in
generating our initial conditions. We measure precession rates in the
simulations by fitting lines to binned arguments of pericentre
$\tilde{\omega}$ that have been averaged over the interior of the disc. The
measured precession frequencies agree reasonably well with the predicted
eigenvalues, except in simulations that we identify as under-resolved (namely
m35l, h5, m5, m5h). For the simulations, we also estimate eccentricity decay
rates (listed as imaginary parts of the frequencies) by fitting slopes to the
natural logarithm of the AMD as a function of time, and assuming that the AMD
decays as $e^{2}.$

The spacetime diagrams in Fig. 5 show radial profiles of radial velocity as a
function of time for most of the hydrodynamic and MHD simulations listed in
Table 1. Sliced at a fixed $\phi=0,$ these spacetime diagrams illustrate the
coherent precession of untwisted eccentric distortions for maximum
eccentricities of $0.2$ and $0.35.$ The modes initialised with $\max[e]=0.5$
involve very strong eccentricity gradients (and hence density variations) near
the inner boundary, coming closer to the limiting eccentric mode shape for our
radial extent $r_{1}/r_{0}=5.$ Their interaction with the inner boundary leads
to shock formation that is visible in the bottom three panels of Fig. 5.
Although the eccentric distortions in these simulations continue to precess,
they clearly take on a different character from the initial conditions.
Papaloizou (2005) observed such shocks in simulations initialised with linear
eccentric modes prescribed a finite amplitude, and Barker & Ogilvie (2016)
excluded them by considering only smaller values of $\max[e].$

For a given value of the maximum eccentricity in the simulation domain, these
spacetime diagrams show very little difference between the hydrodynamic and
MHD simulations; the vertical magnetic field simply causes slightly more rapid
precession. The spacetime diagrams in Fig. 6 illustrate the corresponding
evolution of the vertical field with time. The simulations m35 and m35t have
similar values of $\max[e],$ but two different widths of “envelope” (see
Equation 50) for the vertical magnetic field ($w_{\rm transition}=0.1$ and
$0.05,$ respectively). The different distributions of vertical magnetic flux
do little to alter the characteristics or behaviour of the eccentric mode.

Figure 7: Time-evolution of angular momentum deficit (Equation 31) for all of
the simulations listed in Table 1. The AMD of the magnetised and unmagnetised
discs are nearly indistinguishable (e.g. see the curves for h2 and m2; h35 and
m35; and h5 and m5). Figure 8: Time-evolution of the conserved Hamiltonian
(Equation 57) for all of the simulations listed in Table 1. As with Figure 7
the evolution of the total-Hamiltonian for the magnetised and unmagnetised
disc are nearly indistinguishable. Figure 9: Plot showing simulations’ time
evolution in $\mathcal{H}$-AMD phase space. All ten simulations from Figs. 7
and 8 are represented by points with colours that run from dark to light as
time progresses. The blue and orange curves respectively illustrate the
relationship between $\mathcal{H}$ and AMD for ideal untwisted eccentric
modes, and eccentric distortions with twists set everywhere equal to the
values required for orbital intersection. For better visibility, the y-scale
transitions from log to linear at $\mathcal{H}=0.06.$

Fig. 7 provides a quantitative measure of eccentricity decay, plotting the per
cent change in integrated AMD versus time. For lower eccentricities
($\max[e]=0.2,0.35$) we attribute eccentricity decay both to numerical
diffusion, and to weak damping by the initial growth of the Papaloizou-Pringle
instability (see Barker & Ogilvie, 2016). The level of decay over the
simulation runtimes (of $400T_{0},$ where $T_{0}=2\pi/\Omega_{0}$ and
$\Omega_{0}=\sqrt{GM/a_{\text{in}}^{3}}$ are the orbital period and angular
velocity at the inner boundary) is consistent with the hydrodynamic results
reported by Barker & Ogilvie (2016).

The shock formation in the simulations with larger $\max[e]=0.5$ leads to much
stronger eccentricity damping initially (until $t\Omega_{0}\simeq 400$), and
shallower, “bursty” decay at later times. We attribute this stochastic
evolution to periodic interaction between the strongly modified distortion and
the inner boundary. Fig. 7 quantitatively demonstrates the similarity of the
hydrodynamic and MHD results for a given eccentricity profile and resolution,
regardless of vertical flux distribution (compare m35 and m35t).

The curves in Fig. 8 show the evolution of the total Hamiltonian (Equation
57), which is clearly not conserved in our simulations. Adiabatic damping,
i.e. damping which is slow relative to the precession timescale, should lead
to a slow evolution of the eccentricity along the family of ideal eccentric
modes, towards modes of lower amplitude. As in the unmagnetised case (Ogilvie
& Lynch, 2019), equation 33 implies that an infinitesimal change in the total
Hamiltonian is related to an infinitesimal change in the AMD by
$d\mathcal{H}=-\omega d\mathcal{C}$. The modes in our simulations have
retrograde precession $(\omega<0)$ meaning a decreases in AMD should lead to a
decreasing $\mathcal{H}$, when damping is slow enough. Thus slowly damped
modes should follow the blue curve in Fig. 9 which shows the $\mathcal{H}$-AMD
phase space. However, Fig. 8-9 demonstrates secular growth for the simulations
with $\max[e]=0.2,0.35$ (except for the low-resolution simulation m35l, which
shows similar decay to the simulations with $\max[e]=0.5$).

One potential explanation for this growth is that non-adiabatic damping in our
simulations shifts the initialised eccentric profiles away from the family of
ideal eccentric modes that minimise $\mathcal{H}$ for a given AMD. In
particular when damping is strong enough the eccentric disc will develop a
twist as a result of the disc transporting AMD to compensate for spatial
variations of the damping rate (Ferreira & Ogilvie, 2009). Fig. 9 shows that
the resolved simulations evolve from the untwisted "modal" $\mathcal{H}$-AMD
relation (blue curve) to the "maximally twisted" $\mathcal{H}$-AMD relation
(Orange Curve). The latter is obtained by taking a given eccentric mode and
twisting it until it reaches an orbital intersection everywhere. This is
consistent with the disc gradually twisting, over the course of the
simulation, causing a growth in the total Hamiltonian. This is supported by
the simulations orbital elements, computed using 56, and from looking at the
residual in the radial velocity when the radial velocity of the untwisted
eccentricity profile is subtracted; both of which show the disc becoming
increasingly twisted with time. This twisting of the disc occurs over 1000s of
orbits and is thus much milder than that seen for non-modal initial profiles
(e.g. the const-$e$ profiles studied by Chan et al., 2022) which become highly
twisted over 10s of orbital periods.

Fig. 10 plots profiles of binned eccentricity $\tilde{e}(a)$ at the beginning
(solid lines) and end (dashed lines) of our simulations. For $\max[e]=0.2$ and
$0.35,$ the plot shows the decay of eccentric modes that retain roughly the
same profile in eccentricity, except in m35l (which has half the radial and
azimuthal resolution). Although the profiles for the simulations with
$\max[e]=0.5$ deviate qualitatively, h5, m5 and m5h remain strongly distorted
and relatively untwisted by the end of the simulations. The panels in Fig. 11
show snapshots of radial velocity (top) and vertical magnetic field (bottom)
at the end of the simulations m2 (left), m35 (middle), and m5 (right).

The differences with increasing resolution illustrated by Figs. 7-10 (compare
m35,m35l,m35h, and m5,m5h) indicate that care should be taken in resolving
disk distortions with strong eccentricity gradients. We do not claim to have
completely resolved the eccentric modes’ precession in any of our simulations;
m35l, m35, and m35h demonstrate a clear reduction in AMD decay with increasing
resolution. However, this decay is slow compared with the dynamical timescales
of interest for magnetorotational and parametric instabilities. Further, m35
and m35h exhibit qualitatively similar if not quantitatively identical
evolution.

Figure 10: Binned eccentricity profiles at the beginning (solid) and end
(dashed) of our simulations. The modes with $\max[e]=0.2,0.35$ decay to
qualitatively similar eccentricity profiles (except in the low-resolution run
m35l). On the other hand, shocks in the simulations initialised with eccentric
distortions with $\max[e]=0.5$ lead to strongly modified eccentricity
profiles. Figure 11: Snapshots showing end-state radial velocities (top) and
vertical magnetic fields (bottom) in the simulations m2, m35, and m5. Steep
density gradients in m5 lead to the development of shocks that are likely
unresolved.

##  7 Discussion

At the magnetic field strength relevant to the MRI the magnetic field has
negligible influence on the eccentric modes, which are almost
indistinguishable from their unmagnetised counterparts. In 2D (i.e.
specifically suppressing the MRI and parametric instability) the evolution of
the magnetised and unmagnetised eccentric modes in RAMSES are qualitatively
the same. There are some differences seen between the magnetised and
unmagnetised simulations with $\max[e]=0.5$, however these simulations are not
adequately resolved.

Although the magnetic field has little effect on the eccentric mode, the
presence of an eccentric mode can have a strong influence on the magnetic
field: with lateral compression by the orbital motion, the presence of
eccentricity gradients can enhance the magnetic field strength in regions of
the disc. Similarly the magnetic field strength is reduced in regions of the
disc where the orbital velocity diverges. Despite the strong magnetic field
enhancements, the magnetic field configurations we setup are stable in our 2D
simulations and their slow evolution is consistent with that expected due to
the evolution of the eccentricity profile. In the simulated modes the
enhancement of the magnetic field is primarily a result of the imposition of
circular rigid wall boundaries. However, this effect is potentially very
important in short wavelength/tightly-wound eccentric discs such as those
expected in the inner regions of black hole discs as simulated by Dewberry et
al. (2020b).

One issue that we have encountered is the difficulty of both resolving and
converging the eccentric modes in numerical simulations. This is important for
studies of the eccentric MRI, as having a high enough resolution to resolve
the MRI (e.g. as measured by MRI quality factors) may not be sufficient to
ensure the simulation is well resolved. One also needs adequate horizontal
resolution to resolve the the eccentric mode as well. This is particularly
important if one is interested in analysing the effects of the MRI on the
eccentric disc as the strong damping of the eccentricity by the grid may
overwhelm the effects of the MHD turbulence. Given the relatively strong
damping seen in our 2D simulations assessing the influence of the MRI on the
eccentric disc will prove challenging unless MRI is very efficient at damping
(or in principle exciting) eccentricity.

In this paper we have limited our focus to the 2.5D cylindrical disc setup.
This setup has a number of advantages numerically (easier to implement the
vertical boundary and to achieve adequate vertical resolution), however it
does not give a good approximation to a physical 3D disc. As in hydrodynamic
eccentric discs the variation of vertical gravity and pressure around an
eccentric orbit leads to a dynamically varying scale height around an orbit.
This causes prograde precession of the eccentric disc (Ogilvie, 2001, 2008;
Ogilvie & Barker, 2014; Ogilvie & Lynch, 2019). More importantly the vertical
compression induced by the scale height oscillation can greatly enhance the
quasi-toroidal magnetic fields in nonlinearly eccentric discs (Lynch &
Ogilvie, 2021). Additionally the periodic solution to the induction equation
within the disc needs to match onto the current free external field. This can
be constructed in a circular disc via matched asymptotics (Ogilvie, 1997).
However for non-axisymmetric discs the set of external field solutions (which
can be described using cylindrical harmonics) are generically incompatible
with the field configuration within an eccentric disc (excepting the purely
quasi-toroidal case where no magnetic flux leaves the disc). The internal and
external fields could be connected by a force free transition layer in the
upper disc atmosphere. However such a field configuration is likely unstable
even in the absence of the MRI.

The full 3D problem is important, however, and deserves further attention. A
simpler initial approach might be to simulate a 3D MHD disc while exciting
eccentricity at the outer boundary (Similar to Dewberry et al., 2020b) and
observe the magnetic field response.

##  8 Conclusion

In this paper we have extended the Hamiltonian eccentric disc theory of
Ogilvie & Lynch (2019) to include a magnetic field in an unstratified,
cylindrical geometry. We have solved for the uniformly precessing eccentric
mode solutions of our model and shown that, for magnetic field strengths
relevant to the onset of MRI, the resulting eccentricity profile, and
precession rate, is nearly identical to the unmagnetised case. While such
eccentric modes are of limited utility in describing realistic 3D eccentric
discs due to several important physical effects not being present in the
unstratified geometry, they provide a useful setting for the study of the
eccentric MRI and magnetised parametric instability to further our
understanding of how disc turbulence operates in eccentric discs. More
broadly, this will help inform our understanding of how disc turbulence
operates in flows that vary on the orbital timescale. To this end we confirm
the suitability of our eccentric mode solutions for numerical applications by
using them as initial conditions for 2D MHD simulations in RAMSES. In 2D
simulations we obtain long lived uniformly precessing eccentric flows that
agree closely with the analytical predictions. These flows will provide the
background state for 3D, unstratified, simulations studying the stability of
these eccentric discs to both the MRI and parametric instability which will be
presented in a future publication.

## Acknowledgements

The authors would like to thank Guillaume Laibe and Enrico Ragusa for many
helpful comments on the draft of this manuscript and the anonymous reviewer
for comments and suggestions, which improved the clarity of the paper.

E. Lynch would like to thank the European Research Council (ERC). This
research was supported by the ERC through the CoG project PODCAST No 864965.
This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No 823823.

J. Dewberry gratefully acknowledges support from the Natural Sciences and
Engineering Research Council of Canada (NSERC), [funding reference #CITA
490888-16].

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Barker & Ogilvie (2016) Barker A. J., Ogilvie G. I., 2016, MNRAS, 458, 3739
  * Chan et al. (2018) Chan C.-H., Krolik J. H., Piran T., 2018, ApJ, 856, 12
  * Chan et al. (2022) Chan C.-H., Piran T., Krolik J. H., 2022, arXiv e-prints, p. arXiv:2201.03728
  * Dewberry et al. (2020a) Dewberry J. W., Latter H. N., Ogilvie G. I., Fromang S., 2020a, MNRAS, 497, 435
  * Dewberry et al. (2020b) Dewberry J. W., Latter H. N., Ogilvie G. I., Fromang S., 2020b, MNRAS, 497, 451
  * Faure et al. (2014) Faure J., Fromang S., Latter H., 2014, A&A, 564, A22
  * Ferreira & Ogilvie (2009) Ferreira B. T., Ogilvie G. I., 2009, MNRAS, 392, 428
  * Fromang et al. (2006) Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371
  * Goodchild & Ogilvie (2006) Goodchild S., Ogilvie G., 2006, MNRAS, 368, 1123
  * Kley et al. (2008) Kley W., Papaloizou J. C. B., Ogilvie G. I., 2008, A&A, 487, 671
  * Lynch & Ogilvie (2021) Lynch E. M., Ogilvie G. I., 2021, MNRAS, 501, 5500
  * Miranda et al. (2017) Miranda R., Muñoz D. J., Lai D., 2017, MNRAS, 466, 1170
  * Ogilvie (1997) Ogilvie G. I., 1997, MNRAS, 288, 63
  * Ogilvie (2001) Ogilvie G. I., 2001, MNRAS, 325, 231
  * Ogilvie (2008) Ogilvie G. I., 2008, MNRAS, 388, 1372
  * Ogilvie (2016) Ogilvie G. I., 2016, Journal of Plasma Physics, 82, 205820301
  * Ogilvie (2018) Ogilvie G. I., 2018, MNRAS, 477, 1744
  * Ogilvie & Barker (2014) Ogilvie G. I., Barker A. J., 2014, MNRAS, 445, 2621
  * Ogilvie & Lynch (2019) Ogilvie G. I., Lynch E. M., 2019, MNRAS, 483, 4453
  * Oyang et al. (2021) Oyang B., Jiang Y.-F., Blaes O., 2021, MNRAS, 505, 1
  * Papaloizou (2005) Papaloizou J. C. B., 2005, A&A, 432, 757
  * Pierens et al. (2020) Pierens A., McNally C. P., Nelson R. P., 2020, arXiv e-prints, p. arXiv:2005.14693
  * Ragusa et al. (2017) Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS, 464, 1449
  * Teyssandier & Ogilvie (2016) Teyssandier J., Ogilvie G. I., 2016, MNRAS, 458, 3221
  * Teyssandier & Ogilvie (2017) Teyssandier J., Ogilvie G. I., 2017, MNRAS, 467, 4577
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Whitham (1965) Whitham G. B., 1965, Journal of Fluid Mechanics, 22, 273
  * Wienkers & Ogilvie (2018) Wienkers A. F., Ogilvie G. I., 2018, MNRAS, 477, 4838
  * Zanazzi & Ogilvie (2020) Zanazzi J. J., Ogilvie G. I., 2020, arXiv e-prints, p. arXiv:2009.06636

##  Appendix A Modes with free boundaries

As discussed in Ogilvie & Lynch (2019) the large eccentricity gradients seen
in our simulated modes, which are responsible for the strong variation of the
density around the orbit, are primarily a consequence of imposing rigid
circular boundaries. More realistic free boundaries, appropriate for an
eccentric disc of finite extent, tend to result in smaller eccentricity
gradients for a given value of $\max[e]$, reducing the level of magnetic field
concentration by the eccentric mode.

Following the same procedure laid out in Zanazzi & Ogilvie (2020), one obtains
the free boundary conditions by introducing a taper $T(a)$444Note $T(a)\neq
W(a)$ the taper in the magnetic field as the latter is designed to isolate the
magnetic field from the boundaries, which drops to zero on the disc boundary,
into the disc mass and internal energy by taking $m_{a}\rightarrow m_{a}T(a)$
and $\langle\bar{\varepsilon}\rangle\rightarrow\langle\bar{\varepsilon}\rangle
T(a)$. In an isothermal disc the latter implies the sound speed also drops to
zero on the boundary with $c_{s}^{2}\propto T(a)$. This is in fact a
requirement for the disc to truncate without the forces due to pressure
gradients exceeding those due to gravity and dominating the dynamics of the
fluid in the outer disc. Taking the lengthscale of the taper to zero we
require

$\frac{\partial F_{V_{t},V_{z}}}{\partial f}\Biggl{|}_{a_{\rm min},a_{\rm
max}}=0,$ (58)

in order that the precession frequency remain finite.

In the presence of a taper in the magnetic field (e.g. the magnetic fields
given by Equations 48-49) then Equation 12 simplifies to

$\frac{\partial F^{(\gamma)}}{\partial f}\Biggl{|}_{a_{\rm min},a_{\rm
max}}=0,$ (59)

which is matches the free boundary condition in the unmagnetised disc.

We now consider vertical field setup for the modes computed in Section 5, but
impose Equation 59 for the free boundaries. This results in the modes depicted
in Figure 12 (left). Such mode have non-zero eccentricity on both boundaries
and are thus a challenge to simulate numerically. They do, however, posses the
shallower eccentricity gradients and monotonically decreasing profiles
characteristic of the fundamental mode in more realistic setups. As shown by
Figure 12 (right) the variation of the plasma-$\beta$ around the orbit is much
milder than that seen in the rigid wall case and suggests the fundamental mode
in eccentric disc does not strongly enhance the magnetic field over the
circular value. As discussed in Section 5, this is not the case for higher
order modes which can support larger eccentricity gradients, and thus magnetic
field enhancements, independently of which boundary condition is adopted.

Figure 12: (Left) Example eccentric modes with free boundaries. (Right) Same
as Figure 2 but for the modes with free boundaries. The equivalent modes with
quasi-toroidal fields look broadly similar. There is much less enhancement of
the magnetic field as modes with realistic boundaries are less steep than
those with rigid boundaries.

##  Appendix B Input variables for numerical simulations

###  B.1 Input Variables for the RAMSES simulations

In the $a,\phi$ grid the input variables required for simulations, for a given
eccentricity profile, $e(a)$, are

$\displaystyle\rho$
$\displaystyle=\rho^{\circ}\left(\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}-\frac{ae_{a}}{\sqrt{1-e^{2}}}\frac{\cos\phi+e}{1+e\cos\phi}\right)^{-1}$
(60) $\displaystyle u_{R}$ $\displaystyle=\frac{an}{\sqrt{1-e^{2}}}e\sin\phi$
(61) $\displaystyle u_{\phi}$
$\displaystyle=\frac{an}{\sqrt{1-e^{2}}}(1+e\cos\phi)$ (62) $\displaystyle
B_{z}$ $\displaystyle=\frac{nl_{\rm
MRI}}{2\pi\sqrt{16/15}}W(a)\left(\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}-\frac{ae_{a}}{\sqrt{1-e^{2}}}\frac{\cos\phi+e}{1+e\cos\phi}\right)^{-1}$
(63)

These are then interpolated onto the cylindrical grid by performing a 1D
interpolation from $a$ to $R$ for each azimuthal slice in the ($a$, $\phi$)
grid and copying the resulting 2D fields in the vertical direction (if
present).

###  B.2 Vector Potential

For numerical implementations it is often useful to specify the magnetic field
using a vector potential in order to ensure that the magnetic field is
divergence free. For example transforming from the orbital coordinate system
to cylindrical polars can induce a non-zero velocity divergence if
interpolation is done on the B-field. Instead interpolating the A-field
ensures the resulting magnetic field, in cylindrical polars, obeys the
solenoidal condition.

The vector potential obeys

$\dot{A}_{i}-\epsilon_{ijk}u^{j}\epsilon^{kab}\partial_{a}A_{b}=\partial_{i}f$
(64)

where an overdot indicates a partial derivative with respect to time,
$\epsilon_{ijk}$ indicates the permutation symbol and $f$ is the electrostatic
potential. We are interested in magnetic fields which are steady on the
orbital timescale, for such magnetic fields Equation 64 simplifies to

$2u^{j}\partial_{[j}A_{i]}=\partial_{i}f,$ (65)

where the square brackets around the indices denote that the indices are being
anti-symmetrised, i.e. $T_{[ij]}=\frac{1}{2}(T_{ij}-T_{ji})$. As we did when
computing the $B$-field, we neglect the contribution from the slow rotation of
the field due to disc precession.

In the $(a,M,\tilde{z})$ orbital coordinate system, with stretched vertical
coordinate $\tilde{z}=z/H$, the fluid velocity simplifies to
$u^{i}=n\hat{e}^{i}_{M}$, where $n$ is the mean motion. The steady $B$-field
in a eccentric disc, in the $(a,M,\tilde{z})$ orbital coordinate system is

$B^{i}=B_{t0}(a,\tilde{z})j^{-1}h^{-1}\hat{e}_{\tilde{M}}^{i}+B_{z0}(a)j^{-1}H^{-1}\hat{e}_{\tilde{z}}^{i}.$
(66)

Using the relationship between the $A$ and $B$ fields we find that the $A$
field must obey the following expressions in order to yield the steady
magnetic field solution in an eccentric disc,

$\displaystyle 0$ $\displaystyle=\partial_{[M}A_{\tilde{z}]},$ (67)
$\displaystyle\frac{B_{t0}(a,\tilde{z})}{jh}$
$\displaystyle=-\frac{2}{JH}\partial_{[a}A_{\tilde{z}]},$ (68)
$\displaystyle\frac{B_{z0}(a)}{jH}$
$\displaystyle=\frac{2}{JH}\partial_{[a}A_{M]}.$ (69)

Substituting the expression for the velocity and Equation 67 into Equation 65
we find that the electrostatic potential is a function of semimajor axis only,
$f=f(a)$. This means the $M$ and $\tilde{z}$ components of Equation 65 are
satisfied. To obtain agreement with the $B$ field solution we set $A_{a}=0$
and obtain the following for $A_{M}$ and $A_{\tilde{z}}$,

$\displaystyle A_{M}$ $\displaystyle=\int J^{\circ}B_{z0}(a)\,da,$ (70)
$\displaystyle A_{\tilde{z}}$ $\displaystyle=-\int
H^{\circ}J^{\circ}B_{t0}(a,\tilde{z})\,da.$ (71)

In order to satisfy the $a$ component of Equation 65 we require the
electrostatic potential satisfy

$f(a)=-\int nJ^{\circ}B_{z0}\,da.$ (72)

We now transform to the $(a,E,z)$ orbital coordinate system, which is more
useful for transforming $\boldsymbol{A}$ into other coordinate systems. To do
this we assume that we have a thin disc, allowing us to evaluate the
coordinate transform at the midplane similar to Ogilvie (2018). This works
equally well for the unstratified discs considered in the rest of this paper.
The vector potential in this coordinate system is

$\displaystyle A_{a}$ $\displaystyle=-e_{a}\sin E\int J^{\circ}B_{z0}(a)\,da,$
(73) $\displaystyle A_{E}$ $\displaystyle=(1-e\cos E)\int
J^{\circ}B_{z0}(a)\,da,$ (74) $\displaystyle A_{z}$
$\displaystyle=-\frac{1}{H}\int H^{\circ}J^{\circ}B_{t0}(a,\tilde{z})\,da.$
(75)

The general expression for $\boldsymbol{A}$ in Cartesian and cylindrical
coordinates are fairly complicated, but can be obtained from the above
expressions as follows

$\displaystyle A_{x}$ $\displaystyle=-\frac{1}{J}\left(e_{a}\sin
E\frac{\partial y}{\partial E}+(1-e\cos E)\frac{\partial y}{\partial
a}\right)\int J^{\circ}B_{z0}(a)\,da,$ (76) $\displaystyle A_{y}$
$\displaystyle=\frac{1}{J}\left(e_{a}\sin E\frac{\partial x}{\partial
E}+(1-e\cos E)\frac{\partial x}{\partial a}\right)\int
J^{\circ}B_{z0}(a)\,da,$ (77) $\displaystyle A_{z}$
$\displaystyle=-\frac{1}{H}\int H^{\circ}J^{\circ}B_{t0}(a,\tilde{z})\,da,$
(78)

while the expression in cylindrical polars can be obtained in the usual manor
from $A_{R}=\cos\phi A_{x}+\sin\phi A_{y}$ and $A_{\phi}=-r\sin\phi
A_{x}+r\cos\phi A_{y}$. To evaluate these we require the following expressions
for the partial derivatives of the Cartesian coordinates,

$\displaystyle\frac{\partial x}{\partial a}$ $\displaystyle=(\cos
E-e)\cos\varpi-\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}\sin E\sin\varpi-
ae_{a}\cos\varpi-a(\cos E-e)\varpi_{a}\sin\varpi-a\varpi_{a}\sqrt{1-e^{2}}\sin
E\cos\varpi,$ (79) $\displaystyle\frac{\partial x}{\partial E}$
$\displaystyle=-a\sin E\cos\varpi-a\sqrt{1-e^{2}}\cos E\sin\varpi,$ (80)
$\displaystyle\frac{\partial y}{\partial a}$
$\displaystyle=\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}\sin E\cos\varpi+(\cos
E-e)\sin\varpi-ae_{a}\sin\varpi-a\varpi_{a}\sqrt{1-e^{2}}\sin
E\sin\varpi+a\varpi_{a}(\cos E-e)\cos\varpi,$ (81)
$\displaystyle\frac{\partial y}{\partial E}$
$\displaystyle=a\sqrt{1-e^{2}}\cos E\cos\varpi-a\sin E\sin\varpi.$ (82)

For an untwisted disc where the pericentre direction is aligned with the
$x$-axis these expressions simplify significantly and we arrive at the
following expressions for $\boldsymbol{A}$ in Cartesian coordinates

$\displaystyle A_{x}$ $\displaystyle=-\frac{\sin
E}{J}\left(ae_{a}\sqrt{1-e^{2}}\cos
E+\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}(1-e\cos E)\right)\int
J^{\circ}B_{z0}(a)\,da,$ (83) $\displaystyle A_{y}$
$\displaystyle=\frac{1}{J}\left(-ae_{a}\sin^{2}E+(\cos E-e-ae_{a})(1-e\cos
E)\right)\int J^{\circ}B_{z0}(a)\,da,$ (84) $\displaystyle A_{z}$
$\displaystyle=-\frac{1}{H}\int H^{\circ}J^{\circ}B_{t0}(a,\tilde{z})\,da,$
(85)

while in cylindrical polars $\boldsymbol{A}$ is given by

$\displaystyle A_{R}$ $\displaystyle=\frac{\sin E}{J}\left(-\cos\phi
ae_{a}\sqrt{1-e^{2}}\cos E-\cos\phi\frac{ae_{a}}{\sqrt{1-e^{2}}}(1-e\cos
E)-\sin\phi ae_{a}\sin E-ae_{a}\sqrt{1-e^{2}}\right)\int
J^{\circ}B_{z0}(a)\,da,$ (86)
$\displaystyle\begin{split}A_{\phi}&=\frac{r}{J}\left(\sin E\sin\phi
ae_{a}\sqrt{1-e^{2}}\cos E+\sin\phi\sin
E\frac{1-e(e+ae_{a})}{\sqrt{1-e^{2}}}(1-e\cos
E)-ae_{a}\cos\phi\sin^{2}E+\cos\phi(\cos E-e-ae_{a})(1-e\cos E)\right)\\\
&\times\int J^{\circ}B_{z0}(a)\,da,\end{split}$ (87) $\displaystyle A_{z}$
$\displaystyle=-\frac{1}{H}\int H^{\circ}J^{\circ}B_{t0}(a,\tilde{z})\,da.$
(88)