Magnetorotational Instability in Eccentric Disks

Our analysis is based on the framework laid out by Ogilvie (2001) and Ogilvie & Barker (2014). Ogilvie (2001) introduced the orbital coordinate system (λ, ϕ), illustrated in the top half of Figure 1. A constant-λ contour is an ellipse with semilatus rectum λ and one focus at the origin, and ϕ is the azimuth; the ellipses must vary slowly in orientation and eccentricity over λ so they do not intersect (Ogilvie & Barker 2014). The coordinate system can be extended by adding a vertical coordinate z perpendicular to the plane of ellipses. Using standard methods of Riemannian differential geometry, Ogilvie & Barker (2014) wrote down the components of the ideal MHD equations for adiabatic gas in this nonorthogonal coordinate system.

Zoom In Zoom Out Reset image size

Because particles in the midplane orbit a point mass at the origin along ellipses defining the orbital coordinate system, the orbital coordinate system provides a foundation for extending the shearing box to eccentric disks. To do so, Ogilvie & Barker (2014) chose some reference particle (λ₀, θ(t), 0) and defined a nonorthogonal, shearing-box coordinate system $(\xi ,\eta ,\zeta )\equiv (\lambda -{\lambda }_{0},\phi -\theta (t),z)$ such that ξ∕λ₀, η, and ζ∕λ₀ are $\sim \epsilon$, where  ≪ 1 is the disk aspect ratio. The reference particle and the shearing-box coordinate basis, as well as another basis to be defined in Section 2.3, are shown in Figure 1. Ogilvie & Barker (2014) obtained the velocity perturbation in the shearing box by subtracting from the gas velocity in the inertial frame the velocity of particles following coordinate ellipses, assuming that the velocity perturbation is $\sim \epsilon$ times the particle velocity. Finally, they subtracted the contribution due to orbital motion from the time derivative. This procedure gave them their Equations (83), (84), (86), and (C4)–(C10); these equations are MHD equations which are nonlinear in the velocity perturbation and we derive Equation (1) below from them. We are interested in how the perturbation at the same (λ, ϕ) as the reference particle evolves; thus, we can drop the subscript from λ₀ without ambiguity.

Henceforth, we adopt the orbital coordinate system defined by aligned ellipses of constant eccentricity e. We can convince ourselves that the velocity field defined by particles orbiting the point mass along coordinate ellipses is divergence-free by writing it out explicitly, but we can also see intuitively why this is so: at pericenter, orbits are closer together, but particles move faster.

Because a divergence-free velocity field is incompressible, one solution of the MHD equations in the midplane is that density is a function of λ only, pressure is uniform, magnetic field is vertical and uniform, and gas travels along coordinate ellipses. We choose this solution as the unperturbed background. We consider the particular case where the background density is uniform and, as applies near the midplane, vertical gravity can be ignored; the latter condition means we are looking at the perturbation at fixed ζ above the reference particle.

We specialize the MHD equations of Ogilvie & Barker (2014) according to these assumptions. Because we are considering the perturbation at the same (λ, ϕ) as the reference particle, we elide terms proportional to ξ and η in their Equation (84). Because our background is uniform and we ignore vertical stratification, we discard background spatial gradients and vertical gravity from all their equations. Because the background flow is divergence-free, we set its divergence Δ to zero in their Equations (83), (86), and (C7)–(C9). We retain only terms that are first order in perturbed quantities, and we replace (∂_ξ, ∂_η, ∂_ζ) by $({{ik}}_{\xi },\lambda {{ik}}_{\eta },{{ik}}_{\zeta })$.

To arrive at our form of the linearized MHD equations, we choose a magnetic field unit that absorbs a factor of ${(4\pi )}^{-1/2}$. We let ρ, ${c}_{{\rm{s}}}$, and B be the background density, adiabatic sound speed, and vertical magnetic field, respectively, and v^μ be the contravariant components of the velocity perturbation in the shearing-box coordinate basis. We also denote by u the perturbed logarithmic density times ${c}_{{\rm{s}}}$ and by w^μ the magnetic field perturbation divided by ρ^1∕2. The linearized MHD equations are then

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dM}}\left(\begin{array}{c}u\\ {v}^{\xi }\\ \lambda {v}^{\eta }\\ {v}^{\zeta }\\ {w}^{\xi }\\ \lambda {w}^{\eta }\\ {w}^{\zeta }\end{array}\right)\,=\,\left(\begin{array}{ccccccc}0 & -{{if}}_{\xi }^{{\rm{s}}} & -{{if}}_{\eta }^{{\rm{s}}} & -{{if}}_{\zeta }^{{\rm{s}}} & 0 & 0 & 0\\ {{\mathsf{F}}}_{\lambda }^{{\rm{s}}} & 0 & {\mathsf{A}} & 0 & {{if}}_{\zeta }^{{\rm{m}}} & 0 & {{\mathsf{F}}}_{\lambda }^{{\rm{m}}}\\ {{\mathsf{F}}}_{\phi }^{{\rm{s}}} & {\mathsf{B}} & {\mathsf{C}} & 0 & 0 & {{if}}_{\zeta }^{{\rm{m}}} & {{\mathsf{F}}}_{\phi }^{{\rm{m}}}\\ -{{if}}_{\zeta }^{{\rm{s}}} & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & {{if}}_{\zeta }^{{\rm{m}}} & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & {{if}}_{\zeta }^{{\rm{m}}} & 0 & {\mathsf{D}} & {\mathsf{E}} & 0\\ 0 & -{{if}}_{\xi }^{{\rm{m}}} & -{{if}}_{\eta }^{{\rm{m}}} & 0 & 0 & 0 & 0\end{array}\right)\left(\begin{array}{c}u\\ {v}^{\xi }\\ \lambda {v}^{\eta }\\ {v}^{\zeta }\\ {w}^{\xi }\\ \lambda {w}^{\eta }\\ {w}^{\zeta }\end{array}\right),\end{eqnarray} \tag{ 1 }$$

and the solenoidal condition for the magnetic field reads

$$\begin{eqnarray}&&{{if}}_{\xi }^{{\rm{m}}}{w}^{\xi }+{{if}}_{\eta }^{{\rm{m}}}(\lambda {w}^{\eta })+{{if}}_{\zeta }^{{\rm{m}}}{w}^{\zeta }=0.\end{eqnarray} \tag{ 2 }$$

The time variable in Equation (1) is the mean anomaly M of the reference particle measured from the pericenter, related to the mean motion n by $M={nt};$ in other words, M∕(2π) equals time in units of orbital periods. The background Alfvén speed is ${v}_{{\rm{A}}}\equiv B/{\rho }^{1/2};$ from this, we derive the Alfvén parameter ${f}_{\mu }^{{\rm{m}}}\equiv {k}_{\mu }{v}_{{\rm{A}}}{/}n$, which compares the frequencies of MHD waves and mean orbital motion. The acoustic parameter ${f}_{\mu }^{{\rm{s}}}\equiv {k}_{\mu }{c}_{{\rm{s}}}{/}n$ does the same for sound waves. The other matrix elements are

$$\begin{eqnarray}&&{\mathsf{A}}(M)\equiv -2{{\rm{\Gamma }}}_{\phi \phi }^{\lambda }{\rm{\Omega }}{/}\lambda =2{\rm{\Omega }}{\mathsf{G}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\mathsf{B}}(M)\equiv -\lambda ({{\rm{\Omega }}}_{\lambda }+2{{\rm{\Gamma }}}_{\lambda \phi }^{\phi }{\rm{\Omega }})=-\tfrac{1}{2}{\rm{\Omega }},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\mathsf{C}}(M)\equiv -({{\rm{\Omega }}}_{\phi }+2{{\rm{\Gamma }}}_{\phi \phi }^{\phi }{\rm{\Omega }})=-2{\rm{\Omega }}{\mathsf{H}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\mathsf{D}}(M)\equiv \lambda {{\rm{\Omega }}}_{\lambda }=-\tfrac{3}{2}{\rm{\Omega }},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\mathsf{E}}(M)\equiv {{\rm{\Omega }}}_{\phi }=-2{\rm{\Omega }}{\mathsf{H}},\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{{\mathsf{F}}}_{\lambda }^{{\rm{s}},{\rm{m}}}(M)\equiv -({g}^{\lambda \lambda }{{if}}_{\xi }^{{\rm{s}},{\rm{m}}}+\lambda {g}^{\lambda \phi }{{if}}_{\eta }^{{\rm{s}},{\rm{m}}}),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\mathsf{F}}}_{\phi }^{{\rm{s}},{\rm{m}}}(M)\equiv -(\lambda {g}^{\lambda \phi }{{if}}_{\xi }^{{\rm{s}},{\rm{m}}}+{\lambda }^{2}{g}^{\phi \phi }{{if}}_{\eta }^{{\rm{s}},{\rm{m}}}).\end{eqnarray} \tag{ 9 }$$

Here

$$\begin{eqnarray}&&{\rm{\Omega }}(M)\equiv {n}^{-1}(d\theta /{dt})={(1-{e}^{2})}^{-3/2}{(1+e\cos \theta )}^{2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\lambda }(M)\equiv {n}^{-1}{\partial }_{\lambda }(n{\rm{\Omega }})=-\tfrac{3}{2}{\rm{\Omega }}/\lambda ,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\phi }(M)\equiv {\partial }_{\phi }{\rm{\Omega }}=-2{\rm{\Omega }}{\mathsf{H}},\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&{\mathsf{G}}(M)\equiv 1/(1+e\cos \theta ),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\mathsf{H}}(M)\equiv e\sin \theta /(1+e\cos \theta ),\end{eqnarray} \tag{ 14 }$$

with θ(M) being the true anomaly of the reference particle. Expressions for the inverse metric g^μν and the Christoffel symbol of the second kind ${{\rm{\Gamma }}}_{\nu \rho }^{\mu }$ are found in Equations (B13)–(B19) of Ogilvie & Barker (2014). Since the matrix elements are either constant or M-dependent with period 2π regardless of e, choosing M as the time variable means that our results are independent of the semimajor axis of the reference particle. In the circular limit, ${\mathsf{A}}$ is the centrifugal force, ${\mathsf{B}}$ relates to the Coriolis force, ${\mathsf{D}}$ encodes orbital shear, and $({\rm{\Omega }},{\mathsf{G}},{\mathsf{H}})=(1,1,0)$.

Circular MRI originates from the destabilization of slow magnetosonic waves in differentially rotating disks (Balbus & Hawley 1998). Since these waves are virtually incompressible, we, like Balbus & Hawley (1991), are motivated to look first for similarly incompressible perturbations in eccentric MRI. Orbital shear creates nonzero horizontal components of the velocity perturbation, so an incompressible perturbation must have a vertical wavevector, that is, $\propto {e}^{{{ik}}_{\zeta }\zeta }$.

The adoption of a vertical wavevector means k_ξ, k_η, and w^ζ all vanish. The incompressible limit is characterized by ${\rm{\Omega }},{f}_{\mu }^{{\rm{m}}}\ll {f}_{\mu }^{{\rm{s}}};$ when the wavevector is vertical, it is equivalent to setting u and v^ζ to zero. Under these two assumptions, Equation (2) is automatically satisfied, while Equation (1) simplifies significantly to

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dM}}\left(\begin{array}{c}{v}^{\xi }\\ \lambda {v}^{\eta }\\ {w}^{\xi }\\ \lambda {w}^{\eta }\end{array}\right)=\left(\begin{array}{cccc}0 & {\mathsf{A}} & {if} & 0\\ {\mathsf{B}} & {\mathsf{C}} & 0 & {if}\\ {if} & 0 & 0 & 0\\ 0 & {if} & {\mathsf{D}} & {\mathsf{E}}\end{array}\right)\left(\begin{array}{c}{v}^{\xi }\\ \lambda {v}^{\eta }\\ {w}^{\xi }\\ \lambda {w}^{\eta }\end{array}\right),\end{eqnarray} \tag{ 15 }$$

where $k\equiv {k}_{\zeta }$ and $f\equiv {f}_{\zeta }^{{\rm{m}}}$ for brevity, and we assume without loss of generality that f ≥ 0. The perturbation consists of the two-dimensional velocity sector (v^ξ, λv^η) and the two-dimensional magnetic sector $({w}^{\xi },\lambda {w}^{\eta })$ in velocity units. The matrix in Equation (15) splits into

$$\begin{eqnarray}\left(\begin{array}{cccc}0 & {\mathsf{A}} & 0 & 0\\ {\mathsf{B}} & {\mathsf{C}} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & {\mathsf{D}} & {\mathsf{E}}\end{array}\right)+\left(\begin{array}{cccc}0 & 0 & {if} & 0\\ 0 & 0 & 0 & {if}\\ {if} & 0 & 0 & 0\\ 0 & {if} & 0 & 0\end{array}\right).\end{eqnarray} \tag{ 16 }$$

The first term is parameterized only by e and the second term only by f; therefore, the behavior of eccentric MRI can be fully understood by studying Equation (15) for all (e, f). The first term is time-dependent and describes how orbital variation excites oscillation within each sector. The second term is time-independent and describes how the background magnetic field couples the two sectors and creates magnetic oscillation. Both oscillations are themselves stable, but their coupling may give rise to instability: for e = 0, instability takes the form of circular MRI (Balbus & Hawley 1991); for e > 0, instability results from an extension of circular MRI to eccentric disks (Section 3.1) or the parametric interaction between velocity and magnetic sectors (Section 4.2).

We define a cylindrical coordinate system (R, φ, z) confocal with the orbital coordinate system, and we equip each point with the cylindrical coordinate basis alongside the shearing-box coordinate basis, as in Figure 1. The coordinate systems are related by $(R,\varphi )=(\lambda {\mathsf{G}},\phi )$, so contravariant components transform as

$$\begin{eqnarray}&&\left(\begin{array}{c}{v}^{R}\\ {{Rv}}^{\varphi }\end{array}\right)={\mathsf{G}}\left(\begin{array}{cc}1 & {\mathsf{H}}\\ 0 & 1\end{array}\right)\left(\begin{array}{c}{v}^{\xi }\\ \lambda {v}^{\eta }\end{array}\right),\end{eqnarray} \tag{ 17 }$$

and similarly for (w^ξ, λw^η). Useful properties of cylindrical components are exposed when we convert Equation (15) to the cylindrical coordinate basis using Equation (17),

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dM}}\left(\begin{array}{c}{v}^{R}\\ {{Rv}}^{\varphi }\\ {w}^{R}\\ {{Rw}}^{\varphi }\end{array}\right)=\left(\begin{array}{cccc}{V}_{1} & {V}_{2} & {if} & 0\\ {V}_{3} & {V}_{4} & 0 & {if}\\ {if} & 0 & {W}_{1} & {W}_{2}\\ 0 & {if} & {W}_{3} & {W}_{4}\end{array}\right)\left(\begin{array}{c}{v}^{R}\\ {{Rv}}^{\varphi }\\ {w}^{R}\\ {{Rw}}^{\varphi }\end{array}\right),\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}\left(\begin{array}{cc}{V}_{1} & {V}_{2}\\ {V}_{3} & {V}_{4}\end{array}\right)\equiv \tfrac{1}{2}{\rm{\Omega }}\left(\begin{array}{cc}{\mathsf{H}} & (1-{e}^{2}){{\mathsf{G}}}^{2}+3\\ -1 & -{\mathsf{H}}\end{array}\right)\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}\left(\begin{array}{cc}{V}_{1} & {V}_{2}\\ {V}_{3} & {V}_{4}\end{array}\right)\left(\begin{array}{cc}0 & -1\\ 1 & 0\end{array}\right)-\left(\begin{array}{cc}0 & -1\\ 1 & 0\end{array}\right)\left(\begin{array}{cc}{W}_{1} & {W}_{2}\\ {W}_{3} & {W}_{4}\end{array}\right)=\tfrac{1}{2}{(1-{e}^{2})}^{-1/2}.\end{eqnarray} \tag{ 20 }$$

If f = 0, the velocity and magnetic sectors decouple, and (v^R, Rv^φ) evolves independently of $({w}^{R},{{Rw}}^{\varphi })$. If $f\ne 0$, Equation (18) admits solutions of the form

$$\begin{eqnarray}&&({v}^{R},{{Rv}}^{\varphi },{w}^{R},{{Rw}}^{\varphi })=({\psi }_{1},{\psi }_{2},{ig}{\psi }_{2},-{ig}{\psi }_{1}),\end{eqnarray} \tag{ 21 }$$

where ψ₁(M) and ψ₂(M) are complex-valued functions and g is a root of

$$\begin{eqnarray}&&g-{g}^{-1}=\tfrac{1}{2}{(1-{e}^{2})}^{-1/2}{f}^{-1}.\end{eqnarray} \tag{ 22 }$$

The correctness of this solution is evident upon substituting Equations (20)–(22) into Equation (18). Equation (18) can therefore be recast into the equivalent form

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dM}}\left(\begin{array}{c}{v}^{R}\\ {{Rv}}^{\varphi }\end{array}\right)=\left(\begin{array}{cc}{V}_{1} & {V}_{2}-{fg}\\ {V}_{3}+{fg} & {V}_{4}\end{array}\right)\left(\begin{array}{c}{v}^{R}\\ {{Rv}}^{\varphi }\end{array}\right).\end{eqnarray} \tag{ 23 }$$

Since the right-hand side of Equation (22) is positive, its two roots satisfy $g\gt 1$ and −1 < g < 0, respectively. Two roots beget two solution families: positive-g solutions have magnetic components that are stronger than velocity components, while negative-g solutions have the opposite situation. The Maxwell stress is stronger than the Reynolds stress when g > 0, and it is the other way around when g < 0. Increasing e or decreasing f causes $| g|$ to move further away from unity, enhancing the contrast between velocity and magnetic components for both families. Because g can be readily inferred from (e, f) and the solution family, we shall report only the time evolution of $({v}^{R},{{Rv}}^{\varphi })$.

Equations (1), (15), (18), and (23) of eccentric MRI and Equation (34) of the toy model to be introduced in Section 4.1 all have the form

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{x}}}{{dt}}={\boldsymbol{A}}(t){\boldsymbol{x}}(t),\end{eqnarray} \tag{ 24 }$$

where ${\boldsymbol{x}}(t)$ is a vector and ${\boldsymbol{A}}(t)$ is a periodic matrix with period T. For eccentric MRI, T is the orbital period. We cannot derive a dispersion relation from this equation, so we turn to the theory of Floquet (1883).

Consider the complex-valued equation

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{X}}}{{dt}}={\boldsymbol{A}}(t){\boldsymbol{X}}(t),\end{eqnarray} \tag{ 25 }$$

where ${\boldsymbol{X}}(t)$ is a matrix. A matrix-valued function ${\boldsymbol{F}}(t)$ is called a fundamental matrix if ${\boldsymbol{F}}(t)$ is a solution of Equation (25) and $\det {\boldsymbol{F}}(t)\ne 0$ for all t. We can convince ourselves that ${\boldsymbol{F}}(t){\boldsymbol{C}}$, where ${\boldsymbol{C}}$ is a constant matrix, is a fundamental matrix if and only if $\det {\boldsymbol{C}}\ne 0$. In addition, ${\boldsymbol{F}}(t+T)$ is also a fundamental matrix.

The fundamental matrix ${\boldsymbol{G}}(t)$ satisfying ${\boldsymbol{G}}(0)={\boldsymbol{1}}$ is called the principal fundamental matrix. Now ${\boldsymbol{G}}(t+T)$ and ${\boldsymbol{G}}(t){\boldsymbol{G}}(T)$ are both fundamental matrices with the same value at t = 0; thus, ${\boldsymbol{G}}(t+T)={\boldsymbol{G}}(t){\boldsymbol{G}}(T)$ by the uniqueness of the solution. In other words, the monodromy matrix ${\boldsymbol{G}}(T)$ advances ${\boldsymbol{G}}(t)$ by a period.

Complex matrices, barring some exceptions such as nilpotent matrices, are diagonalizable. Hence, we set ${\boldsymbol{G}}(T)\,={{\boldsymbol{EDE}}}^{-1}$, where ${\boldsymbol{E}}\equiv ({{\boldsymbol{e}}}_{1},{{\boldsymbol{e}}}_{2},\,...)$ is the matrix of column eigenvectors and ${\boldsymbol{D}}\equiv \mathrm{diag}({\alpha }_{1},{\alpha }_{2},\,...)$ is the diagonal matrix of eigenvalues, also called Floquet multipliers. Consider the fundamental matrix

$$\begin{eqnarray}&&{\boldsymbol{M}}(t)\equiv {\boldsymbol{G}}(t){\boldsymbol{E}};\end{eqnarray} \tag{ 26 }$$

clearly, ${\boldsymbol{M}}(0)={\boldsymbol{E}}$ and ${\boldsymbol{M}}(T)={\boldsymbol{ED}}$. This means that if ${{\boldsymbol{x}}}_{j}(t)$ is the jth column of ${\boldsymbol{M}}(t)$, then ${{\boldsymbol{x}}}_{j}(t)$ solves Equation (24), ${{\boldsymbol{x}}}_{j}(0)={{\boldsymbol{e}}}_{j}$, and ${{\boldsymbol{x}}}_{j}(T)={\alpha }_{j}{{\boldsymbol{e}}}_{j}$.

Our task in solving Equation (24) therefore reduces to finding ${\boldsymbol{G}}(T)$ by numerically integrating Equation (25) over one period with the identity matrix as the initial condition and then computing the eigenvalues and eigenvectors of ${\boldsymbol{G}}(T)$. Each eigenvector ${{\boldsymbol{e}}}_{j}$ produces one mode ${{\boldsymbol{x}}}_{j}(t)$ of the full solution. Note that ${{\boldsymbol{e}}}_{j}$ and ${{\boldsymbol{x}}}_{j}(t)$ are defined up to proportionality.

If we let ${\boldsymbol{P}}(t)\equiv {\boldsymbol{M}}(t)\exp ((-t/T)\mathrm{ln}{\boldsymbol{D}})$, where exp and ln are matrix exponentiation and logarithm, respectively, then

$$\begin{eqnarray}{\boldsymbol{P}}(t+T) & = & {\boldsymbol{M}}(t+T){{\boldsymbol{D}}}^{-1}\exp ((-t/T)\mathrm{ln}{\boldsymbol{D}})\\ & = & {\boldsymbol{G}}(t){\boldsymbol{G}}(T){{\boldsymbol{ED}}}^{-1}\exp ((-t/T)\mathrm{ln}{\boldsymbol{D}})\\ & = & {\boldsymbol{G}}(t){\boldsymbol{E}}\exp ((-t/T)\mathrm{ln}{\boldsymbol{D}})\\ & = & {\boldsymbol{P}}(t).\end{eqnarray} \tag{ 27 }$$

Thus, we can write a mode as

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{j}(t)\equiv {\alpha }_{j}^{t/T}{{\boldsymbol{p}}}_{j}(t),\end{eqnarray} \tag{ 28 }$$

where ${{\boldsymbol{p}}}_{j}(t)$ is periodic with period T. If ${\boldsymbol{A}}(t)$ in Equations (24) and (25) is constant, then ${\boldsymbol{P}}(t)$ is also constant, so all modes are either exponential or sinusoidal, and their respective growth rates or oscillation frequencies are given by the diagonal of $\mathrm{ln}{\boldsymbol{D}}/T$. Note that ${{\boldsymbol{x}}}_{j}(t)$ is periodic if and only if ${\alpha }_{j}$ is a root of unity or ${{\boldsymbol{p}}}_{j}(t)$ is constant.

A mode ${{\boldsymbol{x}}}_{j}(t)$ is stable if $| {\alpha }_{j}| \leqslant 1$ and unstable if $| {\alpha }_{j}| \gt 1$. The stability of Equation (24) depends only on the mode with the largest $| {\alpha }_{j}| ;$ we call this mode the most unstable mode and let the growth per period of Equation (24) be $\gamma \equiv {\max }_{j}\mathrm{ln}| {\alpha }_{j}|$.

Since any fundamental matrix ${\boldsymbol{F}}(t)$ is a solution of Equation (25), we have $d(\mathrm{ln}\,\det \,{\boldsymbol{F}})/{dt}=\mathrm{tr}{\boldsymbol{A}}({\rm{t}});$ in particular,

$$\begin{eqnarray}&&{\alpha }_{1}{\alpha }_{2}\cdots =\,\det {\boldsymbol{G}}(T)=\exp {\int }_{0}^{T}{dt}\,\mathrm{tr}{\boldsymbol{A}}({\rm{t}}).\end{eqnarray} \tag{ 29 }$$

If the integral vanishes, ${\alpha }_{1}{\alpha }_{2}\cdots =\,1$, so γ ≥ 0.

The special nature of the 2 × 2 matrix ${\boldsymbol{A}}(t)$ in Equations (23) and (34) leads to additional useful properties. Since ${\boldsymbol{A}}(t)$ is real, ${\boldsymbol{G}}(T)$ and its trace, α₁ + α₂, are both real; since ${\boldsymbol{A}}(t)$ is traceless, Equation (29) yields α₁α₂ = 1. Thus, either $| {\alpha }_{1}| =| {\alpha }_{2}| =1$, which gives two stable modes; or α₁ and α₂ are both real, which gives one stable mode and one unstable mode.

Consider the latter case. If ${\alpha }_{1},{\alpha }_{2}\gt 0$, then both modes are sign-preserving in the sense that each component of ${{\boldsymbol{x}}}_{j}(t)$ retains the same sign after a period; conversely, if α₁, α₂ < 0, then the two modes are sign-reversing because each component of ${{\boldsymbol{x}}}_{j}(t)$ flips sign. Moreover, Equation (26) implies ${\boldsymbol{M}}(T)={\boldsymbol{G}}(T){\boldsymbol{E}}$, the jth column of which is ${\alpha }_{j}{{\boldsymbol{e}}}_{j}={\boldsymbol{G}}(T){{\boldsymbol{e}}}_{j}$. If ${\boldsymbol{G}}(T)$ is diagonal, then ${\boldsymbol{E}}$ is the identity matrix, and the components of ${{\boldsymbol{e}}}_{j}$ are real. If ${\boldsymbol{G}}(T)$ is not diagonal, then the real matrix ${\boldsymbol{G}}(T)$ mixes the two components of ${{\boldsymbol{e}}}_{j}$ to give ${{\boldsymbol{e}}}_{j}$ times a real scalar α_j; this can only be so if the components share the same complex phase. We can therefore take ${{\boldsymbol{x}}}_{j}(0)={{\boldsymbol{e}}}_{j}$ to be real without loss of generality; Equations (23) and (34) then compel ${{\boldsymbol{x}}}_{j}(t)$ to be real for all t.

The results in the previous paragraphs have important implications for eccentric MRI. Equations (15) and (18) have four modes each, dividing into pairs of two: a pair from solving Equation (23) with g > 0 and another pair from solving the same equation with g < 0 (Section 2.3). Here we showed that each pair comprises either two stable modes or one stable mode and one unstable mode; thus, eccentric MRI can have at most two unstable modes, one for each sign of g. Furthermore, while we can choose (v^R, Rv^φ) in Equation (23) to be real at all times, Equation (21) simultaneously makes $({w}^{R},{{Rw}}^{\varphi })$ imaginary, so the mode is restricted to the velocity sector. Since the perturbation is $\propto {e}^{{ik}\zeta }$ (Section 2.2), selecting a different complex phase means observing at a different height above the reference particle. If the velocity sector is real and the magnetic sector imaginary at some ζ, then the velocity sector is imaginary and the magnetic sector real at $\zeta +\tfrac{1}{2}j\pi {k}^{-1}$ for any odd integer j.

Equations (15), (18), and (23) are numerically integrated over one orbit, and the growth per orbit γ of eccentric MRI is given by the eigenvalue of the monodromy matrix with the greatest complex magnitude (Section 2.4). The left panel of Figure 2 displays γ as a function of (e, f); unstable regions are where γ > 0. Circular disks correspond to e = 0; in agreement with Balbus & Hawley (1991), we find MRI if 0 < f² < 3, that is, if the background magnetic field is weak and the wavenumber is small, and stability otherwise.

Zoom In Zoom Out Reset image size

The behavior for arbitrary (e, f) is more complicated. Bands and horns are, respectively, unstable regions found by solving Equation (23) with positive and negative values of g that satisfy Equation (22); bands further divide into the classical band and parametric bands.

The classical band is the extension of the unstable region of circular MRI to e > 0. It contains positive-g unstable modes, whose magnetic components are larger than their velocity components (Section 2.3). Growth is fastest at $(e,f)\,=(0,\tfrac{1}{4}\sqrt{15})$, with $\gamma =\tfrac{3}{2}\pi$. The width of the classical band, as measured in the f-direction, and its γ at fixed f both fall by a factor of ≈2 from e = 0 to e = 1. This is because these modes grow with the help of constant orbital shear just as in circular MRI, but when e is large, orbital shear, encapsulated by λΩ_λ and Ω_ϕ, is small during the long time spent near apocenter.

Parametric bands contain positive-g unstable modes not included in the classical band; like classical-band modes, the magnetic components of these modes are larger than their velocity components (Section 2.3). Parametric bands appear as banana-shaped unstable regions above the classical band that are, loosely speaking, elongated in the e-direction and stacked in the f-direction. Overall, γ in parametric bands is about half the largest γ for circular MRI; more precisely, bandwidth and γ both increase with e at fixed f up to a broad maximum, then decrease slowly toward e = 1, while γ decreases slightly with increasing f at fixed e. Parametric bands are spaced at f-steps of $\approx \tfrac{1}{2}$, with adjacent parametric bands separated by a narrow but finite gap for all e, and they alternate between sign-preserving and sign-reversing in the f-direction (Section 2.5). The regular spacing of parametric bands, their clean separation from one another, and their narrowing toward e = 0 all suggest an origin related to parametric resonance (Section 4.2).

Horns contain negative-g unstable modes, whose velocity components are larger than their magnetic components (Section 2.3). The right panel of Figure 2 depicts horns in isolation. Their width and γ both increase with e, attaining noticeable width only at e ≳ 0.8; they are also regularly spaced in f. Because horns have smaller γ than bands, they are mostly buried underneath bands in the left panel, which portrays only the most unstable modes of Equations (15) and (18); however, horns do emerge between bands when bandgaps widen at e ≳ 0.95. At such high e, horns have γ approximately half that of bands, so horn modes can be as important as band modes in stirring MHD turbulence.

The largest γ at any given e > 0 is generally a factor of a few smaller than that at e = 0. However, the largest f that permits instability rises rapidly with e from the circular-limit value of $\sqrt{3}$, so the f-range over which MRI operates is substantially wider at e > 0. Horns also tend to fill in bandgaps at e ≳ 0.95, making more values of f susceptible to MRI at high e.

While there is no exponential growth at band and horn edges because γ = 0 there by definition, growth in general may still occur (Appendix A).

The division of eccentric MRI modes into band and horn modes has physical significance, which is most easily appreciated in the e = 0 and f = 0 limits.

The e = 0 limit reproduces circular MRI; these circular modes follow the dispersion relation (e.g., Balbus & Hawley 1991)

$$\begin{eqnarray}&&{(\omega /n)}^{2}={f}^{2}+\tfrac{1}{2}\mp {\left(4{f}^{2}+\tfrac{1}{4}\right)}^{1/2}.\end{eqnarray} \tag{ 30 }$$

The upper sign yields (ω∕n)² < 0 if and only if 0 < f² < 3, while the lower sign has ${(\omega /n)}^{2}\geqslant 1$ for all f. Since the 0 < f² < 3 part of the f-axis in Figure 2 is covered only by the classical band, we associate the upper and lower signs with bands and horns, respectively. Band modes with ${f}^{2}\geqslant 3$ and all horn modes are stable at e = 0; these stable circular modes are destabilized by orbital variation through parametric resonance at 0 < e ≪ 1, producing parametric bands and horns, respectively (Section 4.2).

The f = 0 limit is trickier. Equation (23) does not apply, so we cannot classify modes as band or horn. Moreover, Equations (15) and (18) have nondiagonalizable monodromy matrices, leaving us with just the three modes given in Appendix A. All three modes are stable because the magnetic and velocity sectors decouple (Section 2.3) and the two sectors are individually stable (Section 2.2). The first two modes have vanishing (w^ξ, λw^η) and periodic (v^ξ, λv^η); hence, we identify them as epicycles, or inertial modes with vertical wavevectors. The third mode has vanishing (v^ξ, λv^η), corresponding to the situation where a gas packet is displaced along the orbit without any change in velocity. We call this neutrally stable mode a sliding mode, and it is analogous to the azimuthal displacements in circular disks discussed by Balbus & Hawley (1991). The magnetic field perturbation of this mode is frozen into the background flow, and $\lambda {w}^{\eta }$ varies periodically in proportion to the orbital speed.

Although inertial and sliding modes are, strictly speaking, neither band nor horn, we can associate them with band and horn modes at $f\gt 0$ by studying how these latter modes behave as $f\to 0$. We find that the two band modes merge to the sliding mode, whereas the two horn modes tend independently toward the two inertial modes. Just as azimuthal displacements in circular disks are readily destabilized by orbital shear in the presence of a weak magnetic field (Balbus & Hawley 1991), the sliding mode is destabilized at 0 < f ≪ 1 to produce the classical band. Horn modes are more closely related to epicycles, which are stable in Keplerian disks (e.g., Rayleigh 1917); consequently, horn modes are not destabilized at 0 < f ≪ 1 and are destabilized to any appreciable extent only at e ≳ 0.8.

When e > 0, the matrices in Equations (1), (15), (18), and (23) are time-dependent; hence, unstable modes do not grow at a steady exponential rate, nor do their components bear a constant ratio. Instead, components vary at different paces in the course of an orbit in such a way that they are all multiplied by a common factor after a complete orbit, as guaranteed by Equation (28). It is therefore instructive to examine in detail how unstable modes evolve within a single orbit.

To accentuate the difference between bands and horns, we pick six (e, f) from where they do not overlap, that is, where precisely one mode is unstable (Section 2.5); our selection is indicated by dots in Figure 2. The six unstable modes include a circular band mode, a classical-band mode, two modes from adjacent parametric bands, and two modes from adjacent horns. To determine the time dependence of each mode, we numerically integrate Equation (23) over an orbit with the pericenter value of the mode as the initial condition; Figure 3 plots the resulting trajectory of two components (v^R, Rv^φ) of the mode. We choose the complex phase of the perturbation such that $({v}^{R},{{Rv}}^{\varphi })$ is always real (Section 2.5); this is done purely for ease of visualization and has no physical significance. Equation (21) takes us from (v^R, Rv^φ) to (w^R, Rw^φ), which in this case is purely imaginary. For sign-preserving modes (Section 2.5), the trajectory of the subsequent orbit traces out the same shape magnified by a factor of e^γ > 1. For sign-reversing modes, the magnification is −e^γ < −1; that is, the trajectory is enlarged and inverted with respect to the origin.

Zoom In Zoom Out Reset image size

For the circular band mode at (e, f) = (0, 1) in the top left panel, the mode grows exponentially, the trajectory is straight, and the same growth rate applies to all components of the perturbation. For the classical-band mode at (e, f) = (0.5, 1) in the top right panel, however, orbital variation bends the trajectory away from a straight line; this is symptomatic of the uneven growth of different components within an orbit and is a generic feature of MRI growth in eccentric disks.

For the parametric-band modes at (e, f) = (0.5, 2.5) and (e, f) = (0.5, 3) in the center panels, the middle part of each trajectory, traversed while the gas travels out to the apocenter and back, loops counterclockwise around the origin; the mode does not grow along the loop, as evidenced by the confinement of (v^R)² + (Rv^φ)² to a finite range. The ends, corresponding to pericenter passage, deviate from the loop; the deviation is outward whenever the mode grows and inward whenever the mode decays. Although most growth takes place near pericenter, ${({v}^{R})}^{2}+{({{Rv}}^{\varphi })}^{2}$ does not necessarily increase monotonically throughout pericenter passage. For the f = 2.5 mode, the trajectory makes $\tfrac{3}{2}$ turns around the origin; for the f = 3 mode, the trajectory goes around twice. As we discuss below, this winding number is always integer or half-integer.

For the horn modes at (e, f) = (0.99, 0.85) and (e, f) = (0.99, 1.35) in the bottom panels, the trajectories are qualitatively the same as those of parametric-band modes, except that the apocentric loop is clockwise and γ is generally smaller.

Taking appropriate limits in Equation (23) yields physical insight about the apocentric loop. Near apocenter, the diagonal elements of the matrix in Equation (23) are small because ${\mathsf{H}}\approx 0$. If, additionally,

$$\begin{eqnarray}&&({V}_{2}-{fg})({V}_{3}+{fg})\lt 0,\end{eqnarray} \tag{ 31 }$$

Equation (23) describes stable oscillation in which (v^R, Rv^φ) loosely traces out an ellipse of horizontal-to-vertical axis ratio ${[-({V}_{2}-{fg})/({V}_{3}+{fg})]}^{1/2};$ note that V₂ > 0 and V₃ < 0. Horn modes always satisfy Equation (31) because their $g\lt 0$ (Section 2.3); for them, magnetic tension and orbital forces drive oscillation together. Band modes are harder to handle because Equation (31) is true only for parts of the orbit; for mathematical expedience, we consider Equation (31) only at apocenter, trusting that if it holds at apocenter, then the continuity of $({V}_{2}-{fg})({V}_{3}+{fg})$ over θ would ensure that it holds over a finite range around apocenter as well. Band modes have g > 0 (Section 2.3), thus $-{V}_{3}=\tfrac{1}{2}{\rm{\Omega }}\,=$ $\tfrac{1}{2}{(1-{e}^{2})}^{-3/2}{(1-e)}^{2}\,\leqslant$ $\tfrac{1}{2}{(1-{e}^{2})}^{-1/2}\,=$ $f(g-{g}^{-1})\lt {fg}$, where we used Equation (22) in the fourth step. It follows that Equation (31) is equivalent to ${fg}\gt {V}_{2}=\tfrac{1}{2}{(1-{e}^{2})}^{-1/2}+\tfrac{3}{2}{\rm{\Omega }};$ in other words, oscillation occurs if magnetic tension beats orbital forces. The solution of the last inequality is ${f}^{2}\gt \tfrac{3}{2}{(1+e)}^{-3}(2-e)$, which includes all parametric bands and part of the classical band, so modes there exhibit apocentric loops. Furthermore, Equation (22) yields ${[\partial ({fg})/\partial f]}_{e}=2g/(1+{g}^{2})\gt 0$; hence, when f is larger, $({V}_{2}-{fg})({V}_{3}+{fg})$ at apocenter is more negative, Equation (31) is satisfied over a larger fraction of the orbit around apocenter, and the trajectory spends more time looping near apocenter and less time growing near pericenter. This is exactly the trend suggested by the three modes at e = 0.5 in Figure 3.

The argument in the previous paragraph explains why the apocentric loop exists and the ends deviate from the apocentric loop. However, it is not very useful near pericenter. For bands, the argument is frustrated by the fact that the signs of ${V}_{2}-{fg}$ and ${V}_{3}+{fg}$ depend on (e, f). For horns, the argument hardly matters because their unstable modes appear only at e ≳ 0.8 (Section 3.1); thus, little time is spent where ${\mathsf{H}}\approx 0$. Instead, we advance another argument applicable to the near-pericenter evolution of modes along the midlines of parametric bands and horns. Because γ reaches a local maximum there, we can reasonably expect that both ends would be growing, which simplifies our considerations; the behavior along the midline is also likely characteristic of the entire parametric band or horn. Growth of midline modes is concentrated near pericenter because orbital shear is necessary to draw out magnetic field perturbations, and orbital shear is the strongest relative to magnetic tension in that part of the orbit. The time an orbit spends near pericenter is $\sim {[{\rm{\Omega }}(\theta =0)]}^{-1}$, and the instantaneous growth rate is roughly the orbital shear, that is, $\sim {\rm{\Omega }}(\theta =0);$ their product is therefore always of order unity, which may explain why γ varies weakly with e in Figure 2.

The fact that unstable modes grow by a real factor every orbit (Section 2.5) means that (v^R, Rv^φ) at all pericenters must lie on a line that includes the origin, as in Figure 3. Sign-preserving modes have integer winding numbers because $({v}^{R},{{Rv}}^{\varphi })$ at successive pericenters are on the same side of the origin; sign-reversing modes have half-integer winding numbers because (v^R, Rv^φ) switches sides every orbit (Section 2.5). Each band or horn has a single winding number; we call twice this number its order q. The classical band has order q = 0. The lowermost parametric band and horn have orders q = 1 and q = 3, respectively; each band or horn above is one order higher. The winding number increases with f because stronger magnetic tension drives faster oscillation around the apocentric loop.

Denote the shearing-box and cylindrical coordinate bases by $({\hat{{\boldsymbol{e}}}}_{\lambda },{\hat{{\boldsymbol{e}}}}_{\phi })$ and $({\hat{{\boldsymbol{e}}}}_{R},{\hat{{\boldsymbol{e}}}}_{\varphi })$, respectively, and their normalized versions by $({\hat{{\boldsymbol{e}}}}_{\hat{\lambda }},{\hat{{\boldsymbol{e}}}}_{\hat{\phi }})$ and $({\hat{{\boldsymbol{e}}}}_{\hat{R}},{\hat{{\boldsymbol{e}}}}_{\hat{\varphi }});$ the normalized coordinate bases are depicted in the top half of Figure 1. Recall that in differential geometry, a coordinate basis is defined to be tangent to the coordinate curves; specifically, ${\hat{{\boldsymbol{e}}}}_{\lambda }$ is tangent to curves of constant ϕ and z, while ${\hat{{\boldsymbol{e}}}}_{R}$ is tangent to curves of constant φ and z. Since these two sets of curves coincide, we have ${\hat{{\boldsymbol{e}}}}_{\lambda }\,\parallel \,{\hat{{\boldsymbol{e}}}}_{R}$ and ${\hat{{\boldsymbol{e}}}}_{\hat{\lambda }}={\hat{{\boldsymbol{e}}}}_{\hat{R}}$.

The $\hat{\varphi }\hat{R}$-component of the Reynolds stress tensor is ${T}_{\hat{R}}^{\hat{\varphi }}={v}^{R}{{Rv}}^{\varphi };$ the associated angular momentum and energy fluxes are ${{RT}}_{\hat{R}}^{\hat{\varphi }}$ and $R{\rm{\Omega }}{T}_{\hat{R}}^{\hat{\varphi }}$, respectively. Figure 4 plots the two fluxes over an orbit for a circular band mode, a classical-band mode, a parametric-band mode, and a horn mode. Fluxes in the subsequent orbit have the same shape, but the overall normalization is e^2γ > 1 times greater (Section 2.4). Since (v^R, Rv^φ) is chosen to be real (Section 2.5), Figure 4 shows fluxes at ζ such that the Reynolds stress is the greatest. The Maxwell stress reaches its maximum at a different ζ (Section 2.5); the Maxwell stress there is g² times the Reynolds stress here, where $| g| \gt 1$ for band modes and $| g| \lt 1$ for horn modes (Section 2.3).

Zoom In Zoom Out Reset image size

The $\hat{\varphi }$-momentum flux in the $\hat{\lambda }$-direction ${T}_{\hat{\lambda }}^{\hat{\varphi }}$ is obtained by performing a coordinate transformation from cylindrical to shearing-box on the lower index of ${T}_{\hat{R}}^{\hat{\varphi }}.$ Lower indices transform covariantly as the basis; thus, ${T}_{\hat{\lambda }}^{\hat{\varphi }}={T}_{\hat{R}}^{\hat{\varphi }}$. Intuitively, ${T}_{\hat{\lambda }}^{\hat{\varphi }}$ is the flux of angular momentum through the elliptical, constant-λ line segment in the bottom half of Figure 1, and ${T}_{\hat{R}}^{\hat{\varphi }}$ is the same through the circular, constant-R line segment. The greater length of the former line segment is made up for by its obliquity to the latter; hence, the two fluxes are the same.

For the circular band mode at (e, f) = (0, 1) in the first row, fluxes grow exponentially, yielding straight lines on semi-logarithmic plots. For the classical-band mode at (e, f) = (0.5, 1) in the second row, orbital variation bends the trajectory of (v^R, Rv^φ) away from exponential growth (Section 3.3), so fluxes do not increase monotonically. Note that not all classical-band modes have ${T}_{\hat{R}}^{\hat{\varphi }}\gt 0$ throughout the orbit. If e and f are both large, v^R switches sign and then switches back pre-pericenter; if e is large but f is small, Rv^φ changes sign in like manner post-pericenter. In either case, ${T}_{\hat{R}}^{\hat{\varphi }}\lt 0$ over a fraction of the orbit.

For the parametric-band mode at (e, f) = (0.5, 3) in the third row and the horn mode at (e, f) = (0.99, 1.35) in the fourth row, the trajectory of (v^R, Rv^φ) makes a roughly elliptical loop around the origin near apocenter (Section 3.3); ${v}^{R}{{Rv}}^{\varphi }$ therefore oscillates almost sinusoidally over that part of the orbit, and ${T}_{\hat{R}}^{\hat{\varphi }}$ is also more or less sinusoidal considering that R varies slowly there. The sign change of ${T}_{\hat{R}}^{\hat{\varphi }}$ means that stresses sometimes move angular momentum and energy outward and sometimes inward; the sinusoidal nature of ${T}_{\hat{R}}^{\hat{\varphi }}$ near apocenter leads to strong cancellation between outward and inward fluxes. Because the modes in Figure 4 have relatively small f, cancellation may not be very conspicuous; at large f, however, where the apocentric loop covers more of the orbit and the winding number is large (Section 3.3), ${T}_{\hat{R}}^{\hat{\varphi }}$ goes through many periods of sinusoidal oscillation near apocenter, so we anticipate close to complete cancellation. Near pericenter, the mode grows, and ${T}_{\hat{R}}^{\hat{\varphi }}$ can be larger post-pericenter than pre-pericenter or the other way around; this asymmetry can create net transport.

Parametric-band and horn modes fail to grow (Section 3.3), and ${T}_{\hat{R}}^{\hat{\varphi }}$ integrates to a vanishing value over the same part of the orbit, namely, the apocentric loop. The concurrence is unsurprising. Orbital shear feeds the perturbation by draining energy from the background flow and converting it to the kinetic and magnetic energy of the perturbation; in doing so, orbital shear establishes a positive correlation between v^R and ${{Rv}}^{\varphi }$, which leads to net outward transport. Along the apocentric loop, orbital shear is too weak compared to magnetic tension to do either.

Net transport changes osculating orbital elements across the disk. When e ≈ 1, a small increase in the argument of pericenter at one edge of the disk and a corresponding decrease at the other leads to differential apsidal precession and apocentric stream crossing, as seen in Figure 5. This phenomenon is distinct from the differential apsidal precession described by Ogilvie (2001), which requires a radial pressure gradient.

Zoom In Zoom Out Reset image size

The physics of eccentric MRI boils down to the interaction between magnetic and orbital forces: the background magnetic field controls the oscillation between velocity and magnetic sectors (Section 2.2), while orbital variation modulates the strength of this oscillation. We can expose this interaction by eliminating w^ξ and λw^η from Equation (15),

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}}{{{dM}}^{2}}\left(\begin{array}{c}{v}^{\xi }\\ \lambda {v}^{\eta }\end{array}\right)-\left(\begin{array}{cc}0 & {\mathsf{A}}\\ {\mathsf{B}}+{\mathsf{D}} & {\mathsf{C}}+{\mathsf{E}}\end{array}\right)\displaystyle \frac{d}{{dM}}\left(\begin{array}{c}{v}^{\xi }\\ \lambda {v}^{\eta }\end{array}\right)\\ \quad +\,\left(\begin{array}{cc}{f}^{2} & -\dot{{\mathsf{A}}}\\ {\mathsf{BE}}-\dot{{\mathsf{B}}} & {\mathsf{J}}\end{array}\right)\left(\begin{array}{c}{v}^{\xi }\\ \lambda {v}^{\eta }\end{array}\right)=0,\end{eqnarray} \tag{ 32 }$$

where ${\mathsf{J}}(M)\equiv {f}^{2}+{\mathsf{AD}}+{\mathsf{CE}}-\dot{{\mathsf{C}}}$, and overdot denotes differentiation with respect to M. This equation reduces to Equations (106) and (107) of Balbus & Hawley (1998) in the circular limit. The equation describes a pair of coupled, damped oscillators; in the circular limit, the natural frequencies of the two oscillators are f and ${({f}^{2}-3)}^{1/2}$, respectively, but eccentric orbital motion causes periodic modulation of the latter natural frequency.

This observation motivates the study of the toy model

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}x}{{{dt}}^{2}}+{\omega }^{2}(1+h\cos t)x=0\end{eqnarray} \tag{ 33 }$$

as a step toward better understanding eccentric MRI. This equation governs an oscillator whose frequency is periodically modulated around ω; thus, h and ω of the toy model are respectively analogous to e and f of eccentric MRI. The equation can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left(\begin{array}{c}x\\ \dot{x}\end{array}\right)=\left(\begin{array}{cc}0 & 1\\ -{\omega }^{2}(1+h\cos t) & 0\end{array}\right)\left(\begin{array}{c}x\\ \dot{x}\end{array}\right),\end{eqnarray} \tag{ 34 }$$

where $\dot{x}\equiv {dx}/{dt};$ analyzing this equation with the method developed in Section 2.4 results in Figure 6.

Zoom In Zoom Out Reset image size

Figure 6 resembles Figure 2 in multifarious ways. Unstable regions are organized into bands separated by finite gaps. Bandwidth rises with h, with the most rapid change around h = 1. Bands are regularly spaced at ω-steps of $\tfrac{1}{2}$ and are alternately sign-preserving and sign-reversing (Section 2.5). Unstable modes have $(x,\dot{x})$ going clockwise around the origin with winding number $\tfrac{1}{2}q;$ here q is the band order, which is one for the lowermost band and one higher for every band above it. The striking similarities between the two figures are evidence that our toy model captures the essential features of the parametric bands of eccentric MRI, and that what we learn about the former can provide guidance in understanding the latter.

The left side of the qth-order band in Figure 6 appears to pinch off to the point $(h,\omega )=(0,\tfrac{1}{2}q)$. To see whether this is true, we determine band edges at h ≪ 1 by solving Equation (33) perturbatively (Appendix B), as was previously done in the context of parametric resonance (e.g., Landau & Lifshitz 1969). The perturbative and numerical results are in excellent agreement for q ≤ 5, confirming that bands do stretch all the way to the ω-axis. The width of the qth-order band, to leading order in h, is just ${2}^{-3q}{q}^{2q-1}{[(q-1)!]}^{-2}{h}^{q}$ (Bell 1957); hence, the figure cannot resolve its extremely thin tip at h ≪ 1 if q ≥ 3. Since bands extend the instability at 0 < h ≪ 1 due to parametric resonance to finite h, they and their counterparts in eccentric MRI deserve the name "parametric bands."

Figure 6 divides into the small-amplitude regime at h ≤ 1 and the large-amplitude regime at h > 1. The small-amplitude regime is inherently stable because $1+h\cos t\geqslant 0$ throughout the period; parametric resonance can be excited only if 2ω closely matches an integer. In contrast, the large-amplitude regime is inherently unstable because $1+h\cos t\lt 0$ over part of the period, which allows exponential growth for a finite amount of time.

The toy model suggests that parametric bands in eccentric MRI should be understood as the result of orbital variation coupling to magnetic oscillation. At small e, weak orbital variation slightly modulates the frequencies of stable circular modes (Section 3.2). The physics of eccentric MRI in this regime is the same as parametric resonance: almost all of the (e, f)-space at e ≪ 1 and f² ≥ 3 is stable, with instability restricted to ranges around discrete values of f that shrink rapidly as f increases. This tight constraint on f is because parametric resonance demands a close frequency match. At large e, strong orbital variation overwhelms stable circular modes to the degree that exponential growth is possible over a part of the orbit near pericenter. This phenomenon can be viewed as an extension of parametric resonance to large e, but it is opposite to parametric resonance in terms of the f-range that is stable; in other words, now almost all of (e, f)-space is unstable, whereas stability requires a close frequency match.

While there is a visual resemblance only between parametric bands in Figure 2 and bands in Figure 6, it is easy to infer that horns are also part of the parametric phenomenon; after all, they derive from the same Equation (23) with merely a different g. Parametric behavior in MRI involves a gradually increasing degree of destabilization of stable circular modes as e increases from zero, so its defining feature must be that each of its unstable regions narrows leftward and collapses to a point on the f-axis; parametric bands and horns do just that, as we now show. Band and horn modes in the circular limit obey Equation (30). Our toy model suggests that parametric resonance occurs for both bands and horns around $\omega /n=\tfrac{1}{2}q$ for some integer q ≥ 1; the f-coordinates satisfying this condition are marked separately for band and horn modes with gray triangles in Figure 2. Each band and horn clearly converges toward its respective gray triangle as e → 0. It turns out that q here equals the q denoting band or horn order in Section 3.3: the destabilization of a stable circular mode executing $\tfrac{1}{2}q$ oscillations per orbit should produce unstable modes with the same winding number $\tfrac{1}{2}q$ and thus order q.

We can estimate the threshold between the small- and large-amplitude regimes as follows. If we keep only terms proportional to ${d}^{2}(\lambda {v}^{\eta })/{{dM}}^{2}$ and λv^η in the second row of Equation (32), we have

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}}{{{dM}}^{2}}(\lambda {v}^{\eta })+{\mathsf{J}}(\lambda {v}^{\eta })\approx 0.\end{eqnarray} \tag{ 35 }$$

This simplification was used by Balbus & Hawley (1998), and its validity is justified by its results. Because ${\mathsf{J}}$ is a periodic function, this equation describes an oscillator with frequency modulated over M. The properties of the oscillator depend on how the average or baseline of ${\mathsf{J}}$ compares with its amplitude. Observe that

$$\begin{eqnarray}{\mathsf{J}} & = & {f}^{2}-\displaystyle \frac{3-2{e}^{2}}{{(1-{e}^{2})}^{3}}-\displaystyle \frac{7e-4{e}^{3}}{{(1-{e}^{2})}^{3}}\cos \theta \\ & & -\displaystyle \frac{5{e}^{2}-2{e}^{4}}{{(1-{e}^{2})}^{3}}{\cos }^{2}\theta -\displaystyle \frac{{e}^{3}}{{(1-{e}^{2})}^{3}}{\cos }^{3}\theta .\end{eqnarray} \tag{ 36 }$$

We take the baseline of ${\mathsf{J}}$ to be the first two terms on the right-hand side. Since the coefficients of all powers of $\cos \theta$ are negative, we approximate the amplitude of ${\mathsf{J}}$ to be the negative of their sum. In the toy model, the transition from the small- to the large-amplitude regime happens when the amplitude of the frequency modulation equals the baseline frequency, that is, when h = 1. In eccentric MRI, we expect the same transition when the amplitude of ${\mathsf{J}}$ equals its baseline, that is, when ${(1-e)}^{3}{f}^{2}=3-2e;$ this is the threshold plotted in Figure 2, up to a factor of unity. The assumption we made to arrive at Equation (35) is of course ad hoc, but the fact that it reproduces the threshold means that the terms discarded from Equation (32), including drag-like terms and mass-like cross-terms, do not enter into the essence of eccentric MRI.

The threshold between the mostly stable and mostly unstable regimes can also be derived using a more physical argument. MRI grows when orbital shear stretches out magnetic field perturbations; larger f makes the background magnetic field stiffer, restricting the region where orbital shear operates to a smaller fraction of the orbit near pericenter. Circular MRI grows if $f={{kv}}_{{\rm{A}}}/n\lt \sqrt{3};$ by analogy, eccentric MRI should grow if ${{kv}}_{{\rm{A}}}/(n{\rm{\Omega }})\lt \sqrt{3}$ at pericenter, that is, if ${(1-e)}^{3}{f}^{2}\lt 3+3e$. This is very similar to the criterion just derived.

We close this section with an insight regarding eccentric MRI that follows from a contrast between our toy model and the full picture. Eccentric MRI has a classical band at 0 < f² ≲ 3, but the toy model does not have a corresponding band at $0\lt \omega \lesssim \tfrac{1}{2}$. This is because the classical band in eccentric MRI is the extension of circular MRI to e > 0 (Section 3.1), but no such extension is possible for the toy model, which is always stable at h = 0. Thus, orbital variation can drive MRI whether or not the time-averaged orbital shear can do so.

We have treated only the linear stage of MRI in eccentric disks, but MHD stresses in real disks depend on how MRI leads to saturated MHD turbulence, both at what rate and to what final amplitude.

On the one hand, eccentric MRI has a γ that is typically a sizable fraction of the maximum γ of circular MRI (Section 3.1); hence, the number of orbits needed for MRI to go from linear to saturated in eccentric disks may be only a few times that in circular disks. Saturation levels may nevertheless be lower due to the slower linear growth. On the other hand, modes with f² ≥ 3 are linearly stable in circular MRI, so energy can reach those small scales only through the nonlinear, relatively slow process of turbulent cascade from larger scales, whereas modes with ${(1-e)}^{3}{f}^{2}\lesssim 3$ are linearly unstable in eccentric MRI (Section 4.2) and grow right from the start. Since saturation requires a steady state to prevail at all scales, the fact that smaller-scale modes grow sooner in eccentric disks may help MRI saturate faster. Saturation levels may likewise be higher. Whether slower growth or a wider range of unstable wavenumbers is more important can be determined only by nonlinear simulations of the saturation process.

We can only speculate on how the saturated stage of MRI differs in eccentric and circular disks. Self-similar turbulence is characterized by two wavenumbers: a smaller wavenumber corresponding to the scale at which turbulence is driven and kinetic energy is injected and a larger wavenumber corresponding to the scale at which microscopic dissipation converts kinetic energy to internal energy. The inertial range refers to the range between these two wavenumbers; the turbulent power spectrum is a power law in this range. In the linear stage, eccentric MRI is unstable up to wavenumbers ${(1-e)}^{-3/2}$ times larger than circular MRI (Section 4.2); in the saturated stage, it is plausible that the driving range reaches similarly large wavenumbers. However, because small-scale dissipation is independent of large-scale motion, the inertial range should always cut off at about the same wavenumber, so we expect the inertial range at e > 0 to be narrower. For fixed mean motion and vertically integrated pressure, the rate of energy injection at e > 0 may be higher, and the power-law index of the inertial range may be different. It is also possible that band modes with dominant magnetic components could interact nonlinearly with horn modes with dominant velocity components (Section 2.3), leading to quantitative changes in MHD turbulence, especially when e ≈ 1.

So far we have considered incompressible eccentric MRI assuming vertical wavevectors and ignoring vertical gravity; we interpret its unstable modes either as stable circular modes destabilized at 0 < e ≪ 1 by orbital variation through parametric resonance (Section 4.2) or as inertial and sliding modes destabilized at f > 0 (Section 3.2). If we allow for nonvertical wavevectors and vertical gravity, then parametric resonance in hydrodynamic disks can also destabilize inertial and gravity modes (Papaloizou 2005). We may therefore expect MHD disks to generally host destabilized and magnetically modified inertial, sliding, and gravity modes.

The height of thin eccentric disks responds to the modulation of vertical gravity along the orbit, and even mildly eccentric disks can thicken dramatically from pericenter to apocenter (Ogilvie & Barker 2014). Vertical oscillation has the same timescale as orbital variation. The two cooperate in hydrodynamic disks to destabilize inertial modes through parametric resonance (Barker & Ogilvie 2014); the same may happen to the three aforementioned modes in MHD disks.

Vertical oscillation also changes the background density ρ, as well as the vertical wavenumber k of a mode advected with the flow. The resulting modulation of ${v}_{{\rm{A}}}=B/{\rho }^{1/2}$ and $f={{kv}}_{{\rm{A}}}/n$ means modes may switch between stable and unstable within an orbit as the disk shuttles between small- and large-amplitude regimes (Section 4.2). Moreover, stable circular modes may parametrically resonate with orbital motion in a different manner because their f is no longer constant.

Lastly, the background flow of eccentric disks may vary in eccentricity and orientation as a function of semilatus rectum. Horizontal compression and expansion of the background flow may alter inertial, sliding, and gravity modes; it also changes ρ and thus f of stable circular modes.

Our work sheds light on the evolution of the bound debris of TDEs around SMBHs. The debris typically has e ≳ 0.99, so eccentric MRI can grow for values of f that are $\sim {(1-e)}^{-3/2}\gtrsim 1000$ times greater than in the circular limit (Section 4.2). With such a broad f-range linearly unstable, saturation of MRI-driven MHD turbulence may take place in only a few orbits (Section 5.1), so angular momentum transport at the rate associated with a saturated state could begin with relatively little delay.

Our linear formalism tells us how fast MRI amplifies MHD perturbations, not the magnitude of MHD stresses at saturation. Improving on the estimates made by Svirski et al. (2017) of whether angular momentum transport or energy dissipation is more efficient requires nonlinear calculations. Nevertheless, we may expect both effects to be weaker near apocenter if the oscillatory behavior of MHD stresses (Section 3.4) carries over from the linear to the saturated stage and high-f modes dominate at saturation.

MHD stresses may also give rise to differential apsidal precession in the saturated stage, as they do in the linear stage (Section 3.4). Such precession spreads the range of apsidal orientation of the debris, perhaps resulting in weak apocentric shocks. In contrast, general relativistic (GR) bulk apsidal precession rotates every debris orbit through an angle inversely proportional to its pericenter distance. For pericenter distances ≳10 times the gravitational radius, the precession angle is small enough that stream crossing occurs near apocenter (Dai et al. 2015; Shiokawa et al. 2015); for smaller pericenter distances, large swings may lead to closer-in stream crossing and strong shocks. It is unclear whether shocks accompanying MHD and GR precession enhance or diminish the eccentricity of the orbits closest to the SMBH (Bonnerot et al. 2017; Svirski et al. 2017).

If the eccentricity of the inner parts of the debris rises due to either angular momentum transport or shocks, they will plunge directly across the innermost stable circular orbit even though they have lost little orbital energy to radiation. Detailed simulations are required to determine under what circumstances plunging is the likely scenario.