Janosz W Dewberry
email: dewberry@cantab.net

Henrik N Latter

Gordon I Ogilvie

Sebastien Fromang

HFQPOs and discoseismic mode excitation in eccentric, relativistic discs. II. Magnetohydrodynamic simulations

Keywords:

Trapped inertial oscillations (r-modes) provide a promising explanation for high-frequency quasi-periodic oscillations (HFQPOs) observed in the emission from black hole X-ray binary systems. An eccentricity (or warp) can excite r-modes to large amplitudes, but concurrently the oscillations are likely damped by magnetohydrodynamic (MHD) turbulence driven by the magnetorotational instability (MRI). We force eccentricity in global, unstratified, zero-netflux MHD simulations of relativistic accretion discs, and find that a sufficiently strong disc distortion generates trapped inertial waves despite this damping. In our simulations, eccentricities above ∼ 0.03 in the inner disc excite trapped waves. In addition to the competition between r-mode damping and driving, we observe that larger amplitude eccentric structures modify and in some cases suppress MRI turbulence. Given the variety of distortions (warps as well as eccentricities) capable of amplifying r-modes, the robustness of trapped inertial wave excitation in the face of MRI turbulence in our simulations provides support for a discoseismic explanation for HFQPOs.

regions of the black hole accretion disc, where r-modes are trapped (Ferreira & Ogilvie 2009). The excitation mechanism can be described as a three-mode coupling in which a trapped axisymmetric r-mode grows in amplitude by interacting with the eccentricity or warp, and a second, non-axisymmetric inertial wave. In Dewberry et al. (2020) (hereafter Paper I), we numerically demonstrated, with hydrodynamic simulations of eccentric discs, this amplification and its subsequent saturation. Paper I verifies and expands on previous semi-analytical calculations (Ferreira & Ogilvie 2008;Oktariani et al. 2010).

The simulations described in Paper I demonstrate that in a hydrodynamic approximation, the saturation of r-mode growth can lead to the non-linear excitation of higher-frequency, vertically unstructured inertial-acoustic waves (sometimes called p-modes, but here referred to as f-modes to differentiate them from purely acoustic oscillations). The presence of such a secondary coupling supports the idea that the non-linear interaction of multiple discoseismic oscillations might be responsible for the appearance of multiple HFQPOs in some sources (e.g., Ortega-Rodríguez et al. 2014). Additionally, the simulations in Paper I suggest that r-modes can redistribute angular momentum locally, reshaping their own trapping regions and blurring their own frequencies over time.

However, magnetic fields and magnetohydrodynamic (MHD) turbulence complicate the picture. Not only do large-scale, coherent fields alter the geometry of the r-mode trapping region (Fu & Lai 2009;Dewberry et al. 2018Dewberry et al. , 2019)), but previous numerical investigations of MHD-turbulent relativistic discs failed to uncover significant r-mode signatures (Arras et al. 2006;Reynolds & Miller 2009). These numerical investigations did not include an explicit forcing mechanism for the oscillations, though, and therefore could only establish that turbulence driven by the magnetorotational instability (MRI) does not actively excite them (in fact, r-modes are more likely to be damped by turbulence ;Latter & Ogilvie 2006).

In this paper we seek to answer the question of whether sufficient excitation by disc eccentricity might overcome r-mode damping by turbulent fluctuations, generalizing the hydrodynamic simulations described in Paper I to include magnetic fields and MHD turbulence. We employ a vertically unstratified, cylindrical framework, and use the outer radial boundary condition described in Paper I to force eccentricity into discs with pre-saturated MRI turbulence. When streamlines become sufficiently non-circular, we find that r-mode signatures appear in the power spectral density (PSD), with frequencies enhanced by an additional restoring force from magnetic tension. These signatures reveal for the first time the excitation of trapped inertial waves in the presence of turbulence driven by the MRI, and support their viability as an explanation for HFQ-POs. Inertial-acoustic wave excitation again accompanies trapped inertial wave excitation in our simulations. While the non-linear mode coupling seen in Paper I is less apparent, the additional excitation of f-modes over a range of frequencies suggests that both non-linear interactions and mode mixing are at play.

Although we focus on intermediate regimes in which the oscillations coexist with turbulent fluctuations, their excitation is further enhanced by a weakening of the MRI in our eccentric discs: we find that strongly forced eccentricities suppress the Maxwell stress and evacuate magnetic fields. We speculate that the MRI suppression is related to the intense periodic compressions and rarefactions that fluid blobs experience due to the steepening of the eccentric wave near the ISCO (such steepening is predicted by non-linear analysis; see Ogilvie & Lynch 2019;Lynch & Ogilvie 2019). Eccentricity gradients and twists generically produce variations in density (e.g., Ogilvie 2001), although the compression may be less significant far from the black hole. The weakening of the MHD turbulence in our simulations sets the basis for exciting new work exploring the effects of strong disc distortion on the MRI in other environments, not least in black hole accretion discs formed in tidal disruption events (TDEs).

We describe our numerical setup in Section 2, present our simulation results in Section 3 and provide discussion and conclusions in Sections 4 and 5. For theoretical background on discoseismology and the r-mode excitation mechanism, we refer readers to Kato (2001), Dewberry et al. (2019), and Paper I in this series.

NUMERICAL METHODS

In this section we describe our numerical methods, simulation setup, and diagnostics.

Equations and code

As in Paper I, we present simulations run with a version of the code RAMSES (Teyssier 2006;Fromang et al. 2006) that solves the equations of ideal MHD on a uniform cylindrical grid (e.g., Faure et al. 2014),1 under the cylindrical approximation (no vertical gravity; see, e.g., Armitage 1998;Hawley 2001;Sorathia et al. 2012). The continuity, momentum and induction equations solved by RAMSES are given by ∂ρ ∂t

+ ∇ • (ρu) = 0, (1)
∂(ρu) ∂t + ∇ • (ρuu -BB) + ∇ P + B • B 2 = -ρ∇Φ, (2)
∂B ∂t + ∇ • (uB -Bu) = 0, (3)
where ρ, u, P, B and Φ are the mass density, fluid velocity, gas pressure, magnetic field and gravitational potential, respectively. We supplement Equations ( 1)-( 3) with an isothermal equation of state P = c 2 s ρ, for c s the purely constant isothermal sound speed. The sound speed serves as a direct proxy for temperature and disc thickness in our idealised model.

Once again we approximate relativistic effects by utilizing a Paczynski-Wiita gravitational potential. With the exclusion of vertical gravity, this potential is given in cylindrical polars (r, φ, z) by

Φ = -GM r -2r g . (4)
With the gravitational constant G, central black hole mass M and speed of light c set to one, the gravitational radius r g = GM/c 2 and frequency ω g = c 3 /(GM) define code units for space and time. Unless otherwise stated, velocities are then given in units of c.

Particle orbits with angular velocity Ω PW determined by the force balance rΩ 2 PW = ∂ r Φ deviate from the simple Keplerian rotation law Ω K ∝ r -3/2 . In particular, our choice of Φ obliges the square of the horizontal epicyclic frequency κ 2 = 2Ω[2Ω + r∂ r Ω] to fall below zero in the inner disc, defining the innermost stable circular orbit (ISCO) at a radius r ISCO = 6r g . This non-monotonic variation in turn produces the maximum in κ that defines a trapping region for r-modes (see fig. 1 in Paper I). Finally, the orbital period at the ISCO is given by T orb ∼ 61.56ω -1 g .

Initial conditions

We begin with simulations of circular discs, initializing pure rotation on cylinders with u φ = rΩ PW . The saturated states established in these circular disc simulations provide points of comparison, as well as initial conditions for runs in which we impose non-circular streamlines. As in Paper I, we generate disc eccentricity solely with our choice of outer radial boundary condition (see Section 2.3). We place the inner boundary at a smaller radius than the ISCO, within the plunging region; the radial domain is [r 0 , r 1 ] = [4r g , 18r g ] in most of our runs. In all cases we include the full azimuthal range. Unless otherwise stated, we place our vertical boundaries at z = ±H, where H = c s /Ω PW is the isothermal scale height at 8r g (a radial location close to the expected r-mode trapping region). For the value c s = 0.03c taken in all of our simulations, H ∼ 0.5r g , implying an aspect ratio of H/r ∼ 0.064 at r = 8r g . We have taken a larger sound speed than those used in Paper I (where c s = 0.01 -0.02c), in order to allow greater angular resolution. Inertial wave trapping worsens with increasing c s (and hence disc thickness), but large azimuthal resolutions are necessary to resolve MRI turbulence (e.g., Sorathia et al. 2012). As noted in Paper I, such a thin disc may not provide an appropriate description for the SPL emission state, during which large accretion rates make the flow geometry uncertain (see, e.g., Laor & Netzer 1989;Esin et al. 1997;Nayakshin et al. 2000). In any case, a complete model will need to include dissipation in a hot corona, which is excluded by our focus on the mid-plane.

Initially, we set the density to a constant ρ 0 outside r ISCO , and to a floor value within. In order to avoid depletion of the disc due to turbulent transport, we follow Faure et al. (2014) in adding a source term to the continuity equation. Outside 7r g , this source term introduces mass at a rate

∂ρ ∂t = - (ρ -ρ 0 ) τ , ( 5
)
where τ is taken as ten times the local orbital timescale.

With regard to magnetic fields, we set B r = B φ = 0 and initialize purely vertical fields with radial profiles similar to the zeronet-flux configurations considered by Sorathia et al. (2012). Specifically, we set

B z = B 0 S Ω sin[2π(r -r ISCO )], (6)
where S = 1 if 6 < r/r g < 16 and 0 otherwise, and B 0 scales the magnetic field such that the wavelength of the fastest growing MRI mode (λ MRI = 2π √ 16/15V A,z /Ω, where V A,z = B z / √ ρ) remains ≤ H throughout the simulation domain. This configuration produces a plasma beta (ratio of gas to magnetic pressure) with a minimum of β ∼ 80 close to the ISCO.

On top of this MRI-unstable background, we impose whitenoise velocity perturbations |δu i | ≤ 10 -6 c s . We then evolve the MRI-turbulent, circular disc to a saturated, quasi-steady state before 'turning on,' at t = 100T orb , an outer boundary condition that generates the eccentricity. We finally run both the forced and unforced simulations for an additional 100T orb .

Boundary conditions

We implement periodic azimuthal boundary conditions, and radial boundary conditions similar to those implemented in Paper I. A 'diode' outflow condition at the inner boundary sets the radial mass flux in the ghost cells to the value in the innermost active cell if that value is negative (inflowing), or to zero otherwise. Meanwhile, we calculate the azimuthal velocity perturbation to the background orbital motion at the last active cell, and add this to an extrapolation of rΩ PW in each ghost cell. Density, vertical mass flux and magnetic fields are simply matched to their values in the innermost cells of the active domain. We find that altering this inner radial boundary condition for the magnetic field has little effect on results.

At the outer radial boundary, we set u r = u z = B φ = B z = 0, calculate B r to satisfy the solenoidal condition, and match u = rΩ PW φ. In simulations of circular discs we simply set ρ = ρ 0 at r 1 . On the other hand, to produce non-circular streamlines in a given simulation we impose a non-axisymmetric, precessing density profile, ρ E (φ, t), in the outer radial ghost cells. As described in section 3.4 of Paper I, a non-axisymmetric density profile enforced in the ghost cells produces a pressure gradient in the outer disc, in turn forcing the inward propagation of eccentricity. The forcing profile ρ E (φ, t) is chosen to be an eigenmode of the non-linear (Newtonian) secular theory, evaluated at the outer boundary. 2We calculate these eigenmodes (and their associated precession frequencies) semi-analytically, following the procedure described by Barker & Ogilvie (2016), and assume a ratio of outer to inner boundaries r 1 /r ISCO = 3. The fundamental non-linear solutions can be uniquely identified by the maximum eccentricity in the eigenmode, denoted by A f . However, we stress that the actual eccentricity produced in the simulation domain is significantly lower than A f , first because we only force the non-axisymmetric density profile in the ghost cells, and the resulting perturbations in the active domain are smaller amplitude. Additionally, the simulations (unlike the semi-analytical secular theory) are relativistic, and permit continual propagation through the inner boundary. Lastly, damping by dissipation, and more importantly non-linear steepening and potential shocking of the eccentric wave near the ISCO reduce the eccentricity amplitude (see Section 3.2).

We implement periodic vertical boundary conditions, which are appropriate for the cylindrical model. In Paper I we found that these lead to the formation of laminar 'elevator flows' (steady columns of constant u z ) as r-modes' growth saturates. Enforcing a rigid lid at z = ±H halts the formation of these (presumed) numerical artifacts, and in hydrodynamic simulations has little effect on the dynamics other than to force the excited inertial oscillations to form as global standing modes (rather than vertically travelling waves). We find that fluctuations associated with MRI turbulence disrupt the elevator flows in all but the most strongly eccentric discs, and so the need for an artificially rigid vertical boundary is less compelling.

In the absence of vertical gravity, periodic vertical boundary conditions usually preclude the formation of coherent standing waves, which generally need to be inserted by hand (see Appendix A1 in Paper I). We do not fine-tune our initial conditions, and so do not recover standing mode oscillations. But our simulation setup is still sufficient to demonstrate the excitation of radially confined inertial waves. We note that periodic vertical boundary conditions are only formally appropriate under the cylindrical approximation, and can be problematic in stratified simulations (e.g., Salvesen et al. 2016). Given the apparent connection between HFQPOs and coronal plasma, future investigations including vertical structure should certainly consider using boundary conditions that allow the outflow of mass and magnetic flux.

Diagnostics

We use several of the same diagnostics used in Paper I, first of all denoting volume, azimuthal and vertical averages for a given quantity X by

X V = V XdV V dV , (7)
X φ = 1 2π 2π 0 Xdφ, (8)
X z = 1 2H H -H Xdz, (9)
where the volume V covers the radial range [6r g , 16r g ]. Angular momentum transport in MRI-turbulent discs is typically quantified by considering the Maxwell and Reynolds stresses, and so we define

α R = ρu r δu φ V P V , (10)
α M = - B r B φ V P V , (11)
where δu φ = u φ -rΩ PW . We similarly write the volume-averaged plasma beta parameter as

β = P V P M V , (12)
where P M = |B • B|/2 is the magnetic pressure. The quantity ρu 2 r V / P V provides a good quantitative diagnostic for the strength of disc eccentricity, since eccentric deformations dominate the radial kinetic energy (Papaloizou 2005). The complex eccentricity (defined as E = e exp [i ], where e and are the eccentricity and longitude of pericentre) is another measure but is non-trivial to calculate for the twisted distortions in our simulations, because of the non-linearity of the flow and the presence of turbulent fluctuations. For simplicity, we instead use the estimate

Ẽ = u r z cos φdφ π u φ φ,z , (13)
calculated after shifting (globally) in φ to account for the retrograde precession of the distorted disc. This average tracks the imaginary part of E, but we stress that it is correct only to first order in eccentricity, and loses accuracy very close to the ISCO where the eccentric waves steepen and orbits come close to intersecting.

In the hydrodynamic simulations presented in Paper I, we found an appropriate diagnostic for r-mode excitation in ρu 2 z V , as well as the fraction of vertical kinetic energy contained within an annular domain encompassing the r-mode trapping region (r ∈ [7r g , 9r g ]). This diagnostic paints a murkier picture in MRIturbulent simulations, because the turbulence itself contributes to the vertical kinetic energy. We therefore follow Reynolds & Miller (2009) in relying most heavily on timing analysis to search for oscillation excitation. Specifically, we consider the power spectral density (PSD), defined as P(ω) ∝ |F ( f )| 2 , where F ( f ) is the Fourier transform of the signal f (t). P(ω) quantifies the power in a given oscillation frequency ω. We primarily take f to be a suitable average of the radial mass flux, and use angular frequencies in units of ω g throughout.

Hydrodynamic test simulations

In most of the simulations described in this paper, we use boundary conditions that force eccentricity more strongly than in Paper I, where quite weak eccentric structures were sufficient to excite rmodes. The inward-propagating eccentric waves in our MHD simulations, while producing only modest eccentricities (e 0.1), are nonetheless sufficiently non-linear to steepen and circularise near the ISCO (Lynch & Ogilvie 2019). In addition to this natural, non-linear decrease in eccentricity relative to linear predictions, we expect the eccentric waves to suffer turbulent damping due to the MRI. An interesting question is which of the two effects dominates.

To determine the degree to which the distortions' propagation is affected by MHD turbulence vis-a-vis non-linear, conservative circularisation, we have run 2D, purely hydrodynamic simulations using the same sound speed, resolution and forcing boundary conditions as in 3D. These simulations are described in Appendix A and also discussed in Section 3.2.1. In summary, the hydrodynamic 2D simulations exhibit nearly identical disc deformations to their MHD counterparts, suggesting that the turbulent damping of eccentricity does not play a significant role, and that hydrodynamic nonlinear effects control the eccentric wave evolution near the ISCO.

SIMULATION RESULTS

In this section we present the results of our simulations. Table 1 lists the (i) simulation label, (ii) resolution, (iii) forcing amplitude for eccentricity, (iv) Maxwell stress, (v) Reynolds stress, (vi) plasma beta, (vii) volume-averaged radial kinetic energy density (normalised by the kinetic energy associated with sonic motion), and maxima in Ẽ both throughout the entire domain (viii) and near the trapping region (ix). As described in Section 2.2, we use the saturated states achieved by t = 100T orb in the circular runs (zn, znHR, znL, zn2h) as the initial condition for simulations in which we force eccentricity from the outer boundary. We then run both the resulting non-circular discs, as well as the original circular discs, for an additional 100T orb . The values listed in Table 1 have been averaged over the last 50T orb of each simulation.

The resolution used in our runs provides 16 vertical grid cells per H, and a grid cell aspect ratio of dr:rdφ:dz ∼2:2:1 at r = 8r g . The run zn is our fiducial simulation of a circular disc, which we use to produce the initial condition for most of our eccentric disc simulations. In znHR we have doubled both radial and azimuthal resolution (due to numerical expense, this simulation has only been run for 100T orb ). Meanwhile in znL we have extended the radial domain to 30r g , and in zn2h and zn2hAf37 we have doubled the vertical extent (to z ∈ [-2H, 2H]).

Circular disc simulations

Although our primary focus is on eccentric discs, in this section we pause to describe the results of our circular disc simulations. These circular runs are important to consider, first because the quasi-steady states which they achieve provide the initial conditions for our eccentric disc simulations. Second, they provide an opportunity to check the performance of our code, which has been used to explore the linear (Latter et al. 2015) and non-linear (Faure et al. 2014) evolution of the global MRI in Newtonian, but not relativistic contexts. Lastly, our circular disc simulations provide an important point of comparison necessary for identifying r-modes in our timing analysis of the eccentric runs (Section 3.2).

Table 1. Table listing (i) simulation label, (ii) resolution, (iii) maximum eccentricity of the non-linear eigenmode used to generate the density profile enforced at the outer boundary, (iv) Maxwell stress α M = -B r B φ V / P V , (v) Reynolds stress α R = ρu r δu φ V / P V , (vi) plasma beta, (vii) volume-averaged radial kinetic energy density, (viii) maximum in Ẽ (plotted in Fig. 5), and (ix) maximum in Ẽ between 7r g and 9r g . All values have been averaged over the last 50T orb .

Label

N r × N φ × N z A f α M α R β ρu 2 r V / P V max[ Ẽ] max D [ Ẽ]
u 2 z V / P V Figure 1.
Volume-averaged quantities for our fiducial simulations of circular discs. Top left: Maxwell stress. Top right: Reynolds stress. Middle left: Mass fraction. Middle right: inverse plasma beta. Bottom left: radial kinetic energy. Bottom right: vertical kinetic energy. These history plots illustrate the saturated, quasi-steady states after t = 100T orb that provide both the initial conditions and a point of comparison for our eccentric disc simulations.

Volume-averaged quantities

The plots in Fig. 1 show volume-averaged quantities over time for our simulations of circular discs. Fig. 1 (top left) illustrates the saturation of the Maxwell stress at values α M ∼ 0.02 -0.03. These values are consistent with analogous simulations run by Sorathia et al. (2012) with a Newtonian potential. We note that we have also run simulations initialized with a net-flux vertical field, but, following an initial transient these simulations produce nearly identical saturated states to the zero-net-flux runs. We attribute this outcome to advection of the net-vertical flux through the ISCO and inner boundary, since in test simulations with a Newtonian poten-tial we have reproduced the elevated Maxwell stresses observed by Sorathia et al. (2012) in their net-flux runs. Future simulations might consider the effects of a more sustainable source of poloidal magnetic field on disc variability.

Fig. 1 (top right) shows somewhat larger Reynolds stresses than are typically observed in simulations of the MRI, indicating saturated ratios α M /α R ∼ 2 -3 (rather than the usual 4 -6). As described in Section 3.1.2, the additional contribution to the Reynolds stress may be partially due to inertial-acoustic f-modes (i.e., density waves) originating from the inner disc. Meanwhile, Fig. 1 (middle left) shows an equilibration of the total mass (normalized by its initial value), balancing mass accretion through the inner boundary against mass input from the source term in the continuity equation. Fig. 1 (middle right) shows saturation of the inverse plasma beta at β -1 ∼ 0.05, while the bottom panels in Fig. 1 show similar routes to saturation for the volume-averaged radial and vertical kinetic energy densities.

Our circular disc simulations produce magnetic tilt angles (Guan et al. 2009) of θ B = arcsin(α M β)/2 ≈ 11 -12 • , suggesting marginally resolved turbulence (Sorathia et al. 2012). As indicated by the values in Table 1 and Fig. 1, znHR and znL (our simulations with larger resolution and radial domain, resp.) exhibit very similar saturated states to zn, although in znHR the vertical kinetic energy saturates at a marginally higher value. In zn2h, doubling the vertical extent of the disc produces larger stresses, matched by a decrease in β. This is unsurprising in light of unstratified shearing box simulations that, with the exclusion of vertical gravity, show a dependence of the turbulent stress on the vertical extent of the domain (Shi et al. 2016). The cylindrical approximation is invalid on scales larger than a scale height, however, and including a very large vertical extent would eliminate the unstratified model's advantage as an inexpensive numerical framework. In any case, our circular disc simulations fulfill our primary requirement of producing quasi-steady states characterized by sustained, generic MRI turbulence. Further, we observe very similar dynamics in zn2hAf37 (an eccentric disc simulation initialized from zn2h) and znAf37, its counterpart with the vertical range z ∈ [-H, H].

Timing analysis for zn

In order to gain a qualitative understanding of the oscillations (or lack thereof) in the circular disc simulation zn, we consider space-time diagrams in time and radius from the last 100T orb . Those shown in Fig. 2 illustrate radial profiles of radial mass flux. We consider mass flux rather than velocity in order to exclude fluctuations in the low-density region within the ISCO. To isolate modes with different character, we provide spacetime diagrams showing mid-plane profiles of ρu r φ (top), φ = 0 profiles of ρu r z (middle) and profiles of ρu r φ,z (bottom). These averages work as dynamical filters, allowing us to home in on the oscillations of interest. Specifically, taking an azimuthal average isolates axisymmetric modes with azimuthal wavenumber m = 0, which are of the most interest to the r-mode model for HFQPOs. Taking a vertical average, on the other hand, shifts focus to inertial-acoustic f-modes, since they have no vertical structure (i.e., vertical wavenumber k z = 0). Finally, φz averages isolate 'fundamental' axisymmetric inertialacoustic waves with k z = m = 0.

We gain more quantitative insight from power spectral densities P(r, ω) computed from the spacetime data in Fig. 2, which we show via the heatmaps in Fig. 3. We normalize P(r, ω) arbitrarily, but saturate the colour-plots according to the largest peak at frequencies ω > 0.005ω g . The solid and dashed lines plot radial profiles for κ PW and Ω PW calculated for particle orbits in a Paczynski-Wiita potential. These can differ from characteristic orbital and epicyclic frequencies calculated from the actual fluid flow, due to background pressure gradients, magnetic stresses, or a dynamic redistribution of angular momentum by oscillations in the disc (see Paper I). We find, however, that in our MHD-turbulent simulations, restoring forces from magnetic tension do more to modify inertial waves' frequencies and trapping than higher order modifications to κ.

Like Reynolds & Miller (2009), we do not observe explicit signatures of global trapped inertial waves in our simulations of circular discs. Qualitatively, the spacetime diagram for ρu r φ (top panel) indicates some variability in the inner disc, but no dominant axisymmetric oscillations. The broad-band power concentrated beneath the profile for κ (at all radii) in the PSD shown in Fig. 3 (left) suggests that the fluctuations apparent in Fig. 2 (top) are inertial in nature, since local analyses predict that axisymmetric inertial waves are evanescent except where ω 2 < κ 2 . The fact that this power extends to frequencies greater than κ can be explained by an enhancement of small-scale inertial wave frequencies by Alfvénic restoring forces (see, e.g., Fu & Lai 2009;Paris & Ogilvie 2018;Dewberry et al. 2018). Figs. 2 (top) and 3 (left) do not exclude the presence of trapped r-modes, but the plots support the conclusions of Reynolds & Miller (2009) that the MRI does not preferentially excite them over the spectrum of small-scale inertial or, rather, 'Alfvénic-epicyclic' waves.

Meanwhile, the spacetime diagrams in Fig. 2 for ρu r z (middle) and ρu r φ,z (bottom) illustrate the propagation of many inertial-acoustic waves. As described in Paper I, both axisymmetric and non-axisymmetric f-modes may be excited in the inner regions of black hole accretion discs. They are prone to viscous overstability, for instance, in viscous models of relativistic discs possessing a transonic inflow (e.g., Chan 2009; Miranda et al. 2015). In discs with a reflecting inner boundary, on the other hand, nonaxisymmetric f-modes are subject to excitation via a transmission of wave energy at the corotation radius (where the pattern speed matches the orbital angular velocity), due to the profile of vortensity in a relativistic disc (Lai & Tsang 2009;Fu & Lai 2011, 2013).

The PSD for ρu r z shown in Fig. 3 (middle) reveals a spectrum of f-modes across a broad range of m, this time visible as horizontal streaks. The PSD for ρu r φ,z (Fig. 3, right) indicates that one of these can be identified as an axisymmetric f-mode with ω max[κ], located at r ∼ 8r g (close to where κ achieves this maximum). A detailed analysis of the non-axisymmetric f-modes indicated by Fig. 3 (middle) is beyond the scope of this work, but PSDs calculated by taking spatial Fourier transforms in azimuth (rather than time) show that waves with azimuthal wavenumber m = 2 -4 dominate. This is not inconsistent with the excitation of non-axisymmetric f-modes by the corotational instability considered by Lai & Tsang (2009), which produces the largest growth rates for m = 2, 3. But in Paper I we found (as did Miranda et al. 2015) that placing the inner boundary within the ISCO and imposing an outflow boundary condition greatly reduces the impact of this instability, which relies on reflection at the inner boundary. Further, Fu & Lai (2011) found that while corotationally excited f-modes' propagation remains relatively unaltered by large-scale magnetic fields, their growth rates may be reduced. These considerations suggest that the MRI turbulence itself is responsible for exciting the inertial-acoustic waves (see, e.g., Heinemann & Papaloizou 2009a,b, 2012).

Eccentric disc simulations

In this section we present the main results of the paper, describing our eccentric disc simulations. As detailed in Sections 2.2 and 2.3, we initialize these runs with the quasi-steady states achieved by zn and zn2h at 100T orb , and then turn on the eccentricity-forcing outer radial boundary condition. We discuss the inward propagation of eccentricity caused by this boundary condition in Section 3.2.1, and its dynamical consequences in Section 3.2.2.

Eccentricity propagation

In each forced simulation, the precessing, non-axisymmetric density profile enforced at the outer radial boundary induces a travelling eccentric wave that propagates inward at roughly c s /2. Reaching the ISCO within ∼ 10 -20T orb , these waves saturate as steady, slowly precessing distortions. Because the eccentric waves continually propagate through the ISCO and out of the simulation domain (rather than reflecting off an inner boundary), the disc deformations in all simulations manifest as twisted eccentric modes with a radially varying phase.

The volume-averaged quantities plotted in Fig. 4 illustrate the rapid saturation of the eccentric distortions in simulations znAf10, znAf30, znAf35, znAf37 and znAf40. The top panel shows total disc mass (normalized by the total mass at t = 0 initialized in zn), which increases marginally, due to the eccentricity-forcing outer boundary condition, before reaching a new equilibrium. states with larger forcing amplitudes giving rise to larger eccentricities. As shown in the bottom panel, the twisted, travelling eccentric waves enhance the Reynolds stress α R , because they involve a net angular momentum flux. Fig. 5 provides a clearer picture of the structure of the eccentric distortions in simulations znAf10-40, plotting radial profiles of Ẽ, a linear approximation to the complex eccentricity (see Equation 13). We calculate these profiles at a given time t after shifting globally in azimuth by δφ = |ω P |t, where ω P < 0 is the precession frequency enforced at r 1 . This produces profiles of Ẽ that are in phase despite the different precession frequencies associated with different forcing amplitudes for eccentricity. The profiles in Fig. 5 have further been averaged over the final 50T orb of each simulation.

Notably, the plot shows that the eccentric distortions decrease in amplitude toward the inner disc. This decrease disagrees with linear theory, which predicts an increase of eccentricity near the ISCO (see Ferreira & Ogilvie 2009, or the more weakly forced simulations in paper I). However, a recent nonlinear analysis (Lynch & Ogilvie 2019) shows that the steepening of the wave and the occurrence of near intersections between neighbouring streamlines may cause the eccentricity to decrease inwards, even if the angular-momentum flux is conserved. Furthermore, steepening of the wave can greatly enhance dissipation, causing an attenuation of the angular-momentum flux and therefore the wave amplitude.

One might also ask whether turbulent fluctuations in our simulations are also damping the eccentric waves' inward propagation. Figure 5. First-order estimate (see Equation 13) of the imaginary part of the complex eccentricity, Ẽ ∼ -Im[E] = -e sin , averaged over 50T orb . Larger forcing amplitudes A f at the outer boundary produce larger eccentricities in the simulation domain. This linear approximation illustrates the basic radial structure of the eccentric distortions in our simulations, including a decrease in eccentricity near the ISCO due to non-linear wave steepening. We note that as a linear approximation, Ẽ does not accurately capture this non-linear steepening of the eccentric waves, which is clearer in Fig. 6 (top).

Indeed, Ferreira & Ogilvie (2009) found that turbulent damping, parametrised by a bulk viscosity, can hinder the propagation of warps and eccentricities. However, we have run 2D, purely hydrodynamic simulations with the same resolution, boundary conditions and forcing amplitudes as znAf10-40. These hydrodynamic runs show very similar distortions to their MHD counterparts (see Appendix A), demonstrating that the decrease in eccentricity near the ISCO is due primarily to non-linear steepening of the eccentric waves, rather than to turbulent damping. The colour-plots in Fig. 6 (top) illustrate this steepening, showing snapshots of u r z (overlaid by streamlines calculated from the vertically averaged flow) taken 50T orb after initializing our eccentric forcing in simulations (from left to right) znAf10, znAf30, znAf37 and znAf40. The colour-plots on the bottom of Fig. 6 show mid-plane snapshots of Maxwell stress taken at the same time. The streamlines depict the distortion caused by our outer boundary condition, since circular orbits appear as vertical lines on unfolded polar plots (although we note that equivalent radial deviations at larger radii imply smaller eccentricities). In the snapshots from the more strongly forced simulations, sharp changes in u r and very closely spaced streamlines illustrate the inherently non-linear wave steepening not captured by the linear approximation Ẽ plotted in Fig. 5.

Dynamical consequences of eccentricity

Broadly, we identify three dynamical regimes for the simulations with forced eccentricity, depending on the degree of distortion. For forcing amplitudes A f 0.35 (which produce eccentricities of e ∼ 0.007 -0.03 near the trapping region), we do not observe clear signatures of trapped r-modes. However, eccentricities of e 0.03 produced by larger forcing amplitudes do appear to excite trapped inertial waves, which leave distinct timing features in the PSD. Additionally, in this regime the eccentricity begins to significantly in- terfere with the MRI, reducing the Maxwell stress (as illustrated by the sequence of snapshots in Fig. 6, bottom). In the third regime, even larger forcing amplitudes (A f 0.4) produce eccentric distortions that suppress the MHD turbulence within r 10r g . Our simulations in this regime yield similar results to those obtained in hydrodynamic runs utilizing periodic vertical boundary conditions (see Paper I). In the following subsections, we discuss each of the regimes of behavior depicted in Fig. 6 in turn.

Regime I: low-eccentricity

We first discuss simulations znAf10, znAf30 and znAf35, in which forcing boundary conditions induce eccentricities of e 0.03 near the expected r-mode trapping region. As illustrated by Fig. 6 (top left), the velocity perturbations associated with the eccentric distortion in znAf10 possess amplitudes comparable to those of the turbulent fluctuations. Consequently, the eccentric distortion appears only with a time-average (see the darkest curve in Fig. 5). While the distortions in znAf30 and znAf35 are more visible, none of the three runs exhibit signatures of trapped inertial wave excitation.

This absence is indicated both by a lack of distinct variability in timing analyses like that shown in Fig. 3, and by Fig. 7, which plots the fraction of vertical kinetic energy contained within the annular domain r ∈ [7r g , 9r g ]. In Paper I, this fraction increased rapidly in simulations initialized with non-negligible eccentricity, as r-modes trapped near ∼ 8r g grew large enough in amplitude to dominate the total vertical kinetic energy. In contrast, znAf10, znAf30 and znAf35 show little to no change from the saturated value of ∼ 0.25 produced by zn. Averaged over the final 50T orb , the three simulations respectively produce kinetic energy fractions of ∼ 0.27, 0.25, and 0.24.

Given that weaker or comparable eccentricities readily excited r-modes in the hydrodynamic simulations in Paper I, this absence indicates that turbulence driven by the MRI damps the oscillations. It is possible that trapped inertial waves excited by the eccentric structures in these simulations simply do not reach large enough amplitudes to stand out above the noise, but unlikely that such oscillations would coexist with turbulent fluctuations without some interference. Turbulent damping of r-modes would align with the understanding that most contributions from dissipation due to turbulence in an accretion disc will go toward damping oscillations (Latter & Ogilvie 2006).

Regime II: trapped wave excitation

We now turn to the intermediate regime in which the eccentric distortions excite trapped waves but do not completely suppress MRI turbulence. Simulation znAf37 provides the fiducial example, with a curve in Fig. 7 that shows (highly variable) enhancements in vertical kinetic energy within the trapping region r ∈ [7r g , 9r g ].

Like Fig. 2, Fig. 8 shows spacetime diagrams calculated from different averages of the radial mass flux, this time from simulation znAf37. Since the slow precession of the eccentric disc dominates the vertical average ρu r z , we focus on radial profiles of ρu r φ , (top), cos(πz/H)ρu r φ,z (middle), and ρu r φ,z (bottom). Azimuthal averages filter out the disc eccentricity, while the factor of cos(πz/H) removes vertically homogeneous features (i.e., fmodes). This factor also has the disadvantage of removing inertial waves with higher vertical wavenumbers, but still serves to isolate the axisymmetric, vertically structured oscillations of interest. Fig. 9 shows PSDs calculated from the spacetime data pictured in Fig. 8, sampled between t = 120 -200T orb .

The spacetime diagrams in Fig. 8 show enhanced variability over those in Fig. 2. Specifically, the oscillations at ∼ 8r g circumscribed by a white box in Fig. 8 (middle) produce features near the maximum of the epicyclic frequency in the PSDs shown in Fig. 9 (left and middle). In addition to κ PW (white line), the plots in The spacetime diagram and PSD for ρu r φ,z shown in Fig. 8 (bottom) and Fig. 9 (right) again illustrate the radial propagation of axisymmetric inertial-acoustic waves (f-modes), but with greater amplitudes and power than in the analogous timing data for the circular disc simulation zn (compare against Fig. 2, bottom and Fig. 3, right). In Paper I we found that a non-linear self-coupling of trapped inertial oscillations can produce axisymmetric f-modes at twice the r-modes' frequency. The axisymmetric f-modes excited in znAf37 do not produce near-integer ratios as clearly as in Paper I, though, the amplified inertial-acoustic waves' power extending over a range of frequencies both comparable to and greater than those of the (magnetically enhanced) trapped inertial waves. variability in both PSDs. The variability in cos(πz/H)ρu r φ,z illustrates the preferential excitation of (magnetically enhanced) trapped inertial waves. The signal is clearly only quasi-periodic, showing power enhancement over a range of frequencies including max[κ PW ] (solid grey line) and max[κ PW + ω Az ] (dashed line), rather than the isolated peaks exhibited by our more controlled hydrodynamic simulations (cf., fig. 9, right in Paper I). This may be due in part to the dynamical nature of the background flow: while idealised, coherent magnetic fields produce smooth frequency changes in linear r-modes, we do not expect turbulent field fluctuations to provide a time-independent restoring force (Dewberry et al. 2018). From a numerical standpoint, without finely tuned initial conditions or the rigid lid employed in Paper I, we also do not expect the formation of coherent standing wave oscillations in our unstratified, cylindrical framework.

The integrated PSDs for ρu r φ,z (Fig. 10, right) further indicate that significantly more inertial-acoustic waves are excited in the eccentric disc simulation znAf37 than in zn. Enhanced f-mode variability at frequencies greater than those of the trapped inertial waves suggests non-linear coupling, while amplified peaks falling in the inertial range instead indicate mode mixing.3 However, the f-mode spectrum covers too broad a frequency range to identify a frequency commensurability with any confidence.

Lastly, the enhanced Maxwell stress seen in zn2h motivates discussion of trapped inertial wave excitation in simulations with a larger vertical extent, although the cylindrical approximation is formally invalid on vertical scales H. Simulation zn2hAf37 was initialized from the saturated state of zn2h and similarly covers the vertical extent z ∈ [-2H, 2H]. PSDs from the simulation show similar enhancements in variability to those exhibited by znAf37. The Maxwell stress in zn2hAf37 decreases at the same rate as in znAf37 (see below), but, starting from a higher initial value, remains larger at the end of run-time (see Table 1).

Regime III: MRI suppression

In this section we describe the weakening and even suppression of the MRI turbulence by increasingly strong disc distortions, best exhibited by the simulation znAf40. Qualitatively, the sequence of colour-plots in Fig. 6 (top) illustrate a laminarization of the flow, while the bottom colour-plots show a corresponding suppression of the Maxwell stress. Quantitatively, the eccentric distortions' impact on the MRI can be tracked most visibly via the Maxwell stress α M = -B r B φ V / P V , which we plot in Fig. 11. The figure clearly shows a decrease in α M that becomes more severe with increasing eccentricity forcing amplitude. While simulations znAf10-37 maintain modest but finite values of α M , in znAf40 the disc distortion nearly completely suppresses the MHD turbulence by the end of the run-time.

In the absence of MRI turbulence in the inner disc, znAf40 forms the laminar 'elevator flows' described in Paper I. These elevator flows, which are an artifact of periodic vertical boundary conditions in the cylindrical model, halt the growth of trapped inertial waves. They can be eliminated by an artificially rigid vertical boundary condition, which we have chosen to avoid in order to mitigate any spurious numerical effects on the MHD turbulence. Alternatively, the inclusion of vertical gravity and density stratification in future work should naturally eliminate the elevator flows. Their appearance in fact highlights the similarity of oscillations excited in znAf40 to those simulated in Paper I, and suggests that in some cases the suppression of the MRI by disc distortion might allow hydrodynamic instabilities (in particular r-mode excitation) to proceed unimpeded.

DISCUSSION

In this section we discuss the observational relevance of the conditions for the trapped inertial wave excitation demonstrated by our simulations, potential explanations for the weakening and suppression of MRI turbulence in our more strongly eccentric runs, and directions for future work.

Trapped r-mode excitation and HFQPOs

The simulations described in Section 3.2 demonstrate that a sufficiently strong disc eccentricity can enhance variability in the inner regions of relativistic, MHD-turbulent accretion discs. Our timing analysis associates this enhanced variability with Alfvénicepicyclic inertial waves, excited and radially confined near the maximum in κ. This demonstration of trapped inertial wave excitation in the presence of MRI turbulence suggests that magnetic fields may not be as destructive to r-modes as predicted by Fu & Lai (2009) and Reynolds & Miller (2009). Our simulations therefore provide support for the discoseismic explanation for HFQPOs, but with important caveats.

For one thing, we find that in the presence of the MRI turbulence, trapped inertial wave excitation requires a larger eccentricity (e 0.03) near the ISCO than in the laminar, hydrodynamic simulations presented in Paper I of this series, suggesting that MRI turbulence actively damps the trapped modes. This value of eccentricity is likely a lower limit, since the zero-net-flux vertical magnetic fields initialized in our simulations produce only moderate levels of turbulence (α ∼ 0.03 -0.04).

As detailed in Section 2.3 of Paper I, eccentric distortions excited at a 3:1 resonance with the orbit of the X-ray binary (e.g., Lubow 1991) can very reasonably be expected to propagate all the way to the inner disc, especially in the presence of elevated accretion rates associated with the SPL state (see fig. 5 1996;Neil et al. 2007;Zurita et al. 2008;Kosenkov & Veledina 2018) indicate that black hole accretion discs in many of the systems are likely to be eccentric, although actual eccentricities are difficult to estimate, and to our knowledge there have been no direct measurements.

On the other hand, decreases in eccentricity within ∼ 20r g due to non-linear effects, predicted by Lynch & Ogilvie (2019) and exhibited by both our MHD and hydrodynamic (see Appendix A) simulations, may place a limit on eccentricities near the r-mode trapping region. In discs sustaining more vigorous turbulence than that produced by our zero-net-flux, non-radiative simulations, this upper limit on eccentricity may not be enough to overcome turbulent r-mode damping. On yet another hand, the non-linear steepening of eccentric waves could provide a more favorable environment for the distortions' associated hydrodynamic instabilities (like rmode excitation), by first weakening MRI turbulence.

Regardless, we reiterate that the non-linear coupling that excites trapped inertial waves in our simulations is not unique to eccentric distortions; warps provide an alternative route to excitation (Ferreira & Ogilvie 2008). Indeed, early simulations of 'tilted' discs exhibited enhanced variability (Henisey et al. 2009), although in those particular runs the variability may not have been due to standing mode oscillations (Henisey et al. 2012). Practical considerations motivated our focus on the timing properties of eccentric discs, since simulating warped discs requires the full inclusion of vertical gravity, at much greater numerical expense (e.g., Liska et al. 2019). But we expect the fundamental mechanism of r-mode excitation via a non-linear coupling with a non-axisymmetric distortion to translate from eccentric to warped accretion discs.

MRI modification

A complete analysis of the effects of disc distortions on MRI stability lies beyond the scope of this project. However, the distinct radial localization of the MRI-suppression in our simulations is suggestive and deserves some commentary. The Maxwell stress is quenched first near the ISCO, where the eccentricity amplitude is lower but the quasi-stationary eccentric waves steepen and nearly shock (see Fig. 6,top); this points to the strong compression associated with the eccentric waves' non-linear steepening as a possible cause of the MRI's suppression.

The steepening of the eccentric deformations in our simulations can be described formally by enhancements of the eccentricity gradient (change in eccentricity with radial coordinate) and twist (change in longitude of pericentre with radial coordinate). The variation of eccentricity with radius in our simulations constitutes an important difference between our setup and that of Chan et al. (2018), who found persistent instability in their local analysis of uniformly eccentric discs with net vertical magnetic flux.

Although the eccentric distortions in our simulations are essentially stationary from a global standpoint, their nonaxisymmetry means that orbiting fluid elements will suffer intense and periodic flow perturbations every orbit. Nearby fluid elements tethered by magnetic fields will therefore be sucked in and out of a weak shock on the MRI timescale. It is reasonable to expect that this interference would disrupt the classical instability mechanism of the MRI involving angular momentum transport via field lines tethering and pulling apart nearby fluid blobs. Indeed, the local analysis of Chan et al. (2018) suggests that even a uniform distribution of eccentricity (without significant compression) can decrease linear MRI growth rates.

The non-linear steepening of eccentric waves may additionally alter the turbulent dissipation of the magnetic field; studies of resistive plasmas in a variety of contexts indicate that compression can enhance magnetic reconnection (e.g., Hesse et al. 2011;Uzdensky & McKinney 2011). There may also be a parallel with the suppression of zero-net-flux MRI by gravitoturbulence in vertically stratified shearing boxes (Riols & Latter 2018), since the gravitational instability proceeds via spiral density waves with similarly strong density variations.

While the impact of distortions on the MRI may be less significant far away from a black hole, where eccentric modes have much longer wavelengths, compression due to strong eccentricity gradients may still affect MHD turbulence in a wider variety of astrophysical environments. The role played by eccentricities in MHD accretion discs should therefore be explored further in more controlled (i.e., non-relativistic) numerical simulations. Eccentric (Ogilvie & Barker 2014;Wienkers & Ogilvie 2018) and/or warped (Torkelsson et al. 2000;Ogilvie & Latter 2013;Paardekooper & Ogilvie 2019) shearing boxes likely offer the best approach to exploring the effects of disc distortion on the MRI.

Future work

Along with further investigation into the effects of disc distortion on MRI turbulence, several other avenues should be explored so as to better evaluate the viability of the discoseismic model for HFQ-POs. While our simulations demonstrate that sufficient eccentricity can excite trapped inertial waves despite damping by MRI turbulence, we have excluded relevant physics in an effort to simplify a messy problem. For one thing, in the absence of vertical gravity or the rigid vertical boundary condition used in Paper I, the trapped inertial waves excited in our MHD simulations do not form global standing modes in z. Future work would likely require the full inclusion of density stratification to realize the formation of bona fide global modes.

Along with vertical gravity, future investigations of the disco-seismic model for HFQPOs in X-ray binaries should include more detailed radiative and thermal physics relevant to the elevated states of accretion in which HFQPOs are observed. The connection between disc oscillations and hot coronal plasma (which produces the emission on which HFQPOs are imprinted) should also be an essential component of any comprehensive model. Coherent magnetic fields might provide such a connection by facilitating a transfer of variability from oscillations in the disc to time-dependent reconnection and emission in the corona, but exploring this possibility will require more sophisticated simulations (perhaps drawing inspiration from studies of wave-interactions in the solar corona; e.g., Potter et al. 2019;Morton et al. 2019;Riedl et al. 2019).

CONCLUSIONS

This paper describes the results of global, magnetohydrodynamic simulations of relativistic, eccentric accretion discs. These simulations show that trapped inertial waves (r-modes) can grow in the presence of MHD turbulence if excited by sufficiently strong disc distortion. These oscillations are important as they may drive high-frequency, quasi-periodic oscillations seen in the emission from black hole X-ray binaries. The trapped inertial wave growth is accompanied by outwardly propagating, axisymmetric inertialacoustic waves (f-modes). Additionally, our simulations reveal a disruption of the magnetorotational instability in the inner regions of sufficiently distorted discs, likely owing to strong, nonaxisymmetric compression associated with steepening and circularisation of the eccentric waves. Like Reynolds & Miller (2009), we do not observe signatures of trapped r-modes in simulations of circular discs. However, in the presence of disc eccentricity, which we force into the domain with a non-axisymmetric density profile imposed at the outer boundary, we find three dynamical regimes. Simulations with low-amplitude forcing for eccentricity do not exhibit signatures of trapped inertial waves, likely due to r-mode damping by turbulent fluctuations. In runs in which larger forcing amplitudes produce sufficient eccentricity close to the ISCO, on the other hand, we do observe signatures of trapped inertial wave excitation in the power spectral density, along with a weakening of the Maxwell stresses associated with MRI turbulence. For eccentricities large enough that the MRI is completely suppressed and the inner disc is essentially laminar, we observe dynamics similar to the hydrodynamic simulations presented in Dewberry et al. (2020) (Paper I).

Our simulations are simplified (isothermal, unstratified, zeronet-flux), and validation of the discoseismic model for HFQPOs will require the inclusion of more sophisticated physics. Still, this work demonstrates that sufficient disc distortion may, in principle, excite trapped inertial waves in black hole accretion discs despite damping by MHD turbulence.

APPENDIX A: HYDRODYNAMIC TEST SIMULATIONS

In this appendix we describe 2D, purely hydrodynamic simulations (summarised in Table A1) that provide a point of direct comparison for the eccentric distortions driven in our MHD simulations. These runs are motivated primarily by an inward decrease of eccentricity amplitude not seen in the much more weakly forced deformations in Paper I. In exhibiting very similar distortions to the MHD simulations described throughout this paper, the simulations presented in this appendix strongly indicate that the inward decrease in eccentricity is caused primarily by steepening and circularisation of the eccentric waves due to non-linear, purely hydrodynamic effects, rather than turbulent damping. Simulations hAf10-40 use exactly the same sound speed, horizontal resolution, and radial boundary conditions as znAf10-40. Meanwhile, hlAf40 has been run on the much larger radial grid r/r g ∈ [4, 54], with a resolution of N r × N φ = 744 × 800. To generate an eccentricity-forcing boundary condition for this simulation, we take the ratio of outer to inner boundaries as r 1 /r ISCO = 9 in calculating an eigenmode of the non-linear theory considered by Barker & Ogilvie (2016) (see Section 2.3, or Paper I). Note that this calculation produces a weaker forcing than that used in hAf40, which has r 1 /r ISCO = 3, since a larger radial domain implies that a shallower eccentricity gradient is required to produce the same maximum eccentricity in the two (Newtonian) calculations.

The hydrodynamic simulations hAf10-40 produce very similar distortions to their MHD counterparts znAf10-40. The radial profiles of Ẽ plotted in Fig. A1 are understandably smoother than, but still remarkably similar to, those shown in Fig. 5. Fig. A2 further illustrates the similarity, displaying very similar colour plots to those shown in Fig. 6 (top). Comparing Table A1 with Table 1, the hydrodynamic simulations show values of Ẽ that are smaller by 10% than the corresponding MHD simulations in the outer disc (potentially due to an absence of turbulent velocity fluctuations), and larger by 5% near the trapping region. This indicates that the MHD turbulence damps the eccentricity amplitudes by no more than ∼ 10 -15%. The saturated values of α R and ρu 2 r V / P V listed in Table A1 are also similar to the difference between the circular and eccentric MHD simulations zn and znAf10-40, although these quantities are not strictly comparable because of the MRI weakening that takes place in the latter simulations. Lastly, Fig. A3 describes the eccentric distortion forced in the larger domain run hlAf40, showing an unfolded polar snapshot of radial velocity overlaid by streamlines (top left), the density profile enforced in the outer ghost cells (top right), and a radial plot of Ẽ (bottom). The figure demonstrates first that the much weaker forcing at large radii still gives rise to the values of Ẽ ∼ 0.05-0.010 observed for hAf40 and znAf40. Secondly, the plots in Fig. A3 show that the wave-steepening observed in our other simulations is not unique to the shorter radial domain, illustrating an eccentric distortion with a much larger wavelength at larger radii. Importantly, this rapid increase in wavelength with radius implies that turbulent eccentricity damping will be even weaker throughout most of the disc. Eccentric distortions should therefore face little difficulty in propagating to the ISCO, especially in the presence of the elevated mass accretion rates associated with the SPL state.

Figure 2 .

 2 Figure 2. Spacetime diagrams illustrating azimuthally averaged (top), vertically averaged (middle) and both azimuthally and vertically averaged (bottom) radial mass flux over the last 100T orb of runtime for simulation zn. These spacetime diagrams indicate periodic behavior over a range of radii and frequencies, but provide no evidence of trapped wave excitation.

Figure 3 .

 3 Figure 3. Power spectral densities calculated from the spacetime data shown in Fig. 2. The solid and dashed white lines show the profiles of κ PW and Ω PW associated with a Paczynski-Wiita potential. The left-hand PSD for ρu r φ reveals broad-band inertial power, the middle PSD for ρu r z shows a spectrum of inertial-acoustic f-modes propagating at all radii, and the right-hand PSD for ρu r φ,z suggests the localization of an axisymmetric f-mode close to the maximum in κ PW .

Figure 4 .

 4 Figure 4. History plots showing mass fraction (top), volume-averaged radial kinetic energy density (middle), and Reynolds stress (bottom) in simulations znAf10-40 (in which we force eccentric distortion). The change and rapid saturation in these quantities from the initial values provided by the circular run zn illustrate the transition to quasi-steady, twisted eccentric discs, with larger distortions produced by larger forcing amplitudes.

Figure 6 .

 6 Figure 6. Unfolded polar snapshots showing vertically averaged radial velocity and overplotted streamlines (top), and mid-plane Maxwell stress normalized by P 0 = c 2 s ρ 0 (bottom) for simulations (from left to right) znAf10, znAf30, znAf37, and znAf40. The snapshots are taken from simulations with larger and larger forcing amplitudes for eccentricity (from left to right), 50T orb after turning on the outer boundary condition that drives the deformation. The top plots illustrate the eccentric distortion of the disc, while the bottom plots suggest a weakening of the MRI turbulence with increasing eccentricity.

Figure 7 .

 7 Figure 7. The fraction of vertical kinetic energy contained within the annular domain D defined by r ∈ [7r g , 9r g ]. Simulations znAf10, znAf30, and znAf35 show little change from the saturated value taken from the circular run zn. On the other hand, znAf37 and znAf40 show larger deviations, hinting at enhanced variability due to trapped wave excitation.

Figure 8 .

 8 Figure 8. Spacetime diagrams showing radial profiles of ρu r φ , (top) cos(πz/H)ρu r φ,z (middle) and ρu r φ,z (bottom) over time in simulation znAf37. The white box in the middle plot highlights variability not apparent in the circular disc simulation zn, which we associate with (magnetically enhanced) trapped inertial waves excited by the eccentricity in the disc.

  Fig. 9 include radial profiles of κ PW + ω Az , where ω Az = π ṼAz /H is an estimate of the Alfvén frequency for modes with k z = π/H, calculated using a time-average of ṼAz = ( |B z | 2 φ,z / ρ φ,z ) 1/2 from the last 100T orb of zn. Local analyses and semi-analytic calculations predict that in the presence of a strong toroidal magnetic field, rmode frequency enhancements due to the additional restoring force of magnetic tension will be ω Az (see section 3.2 in Dewberry et al. 2019). Fig. 9 (middle) validates this prediction, showing the preferential excitation of magnetically altered inertial oscillations close in both frequency and radius to max[κ PW + ω Az ].

Figure 9 .

 9 Figure 9. Power spectral densities calculated from the spacetime data shown in Fig. 8, over the interval t = 120 -200T orb . As in Fig. 3, the solid and dashed white lines show the profiles of κ PW and Ω PW associated with the Paczynski-Wiita potential, while the dotted line shows κ PW + ω Az for a mode with k z = π/H. The peaks in power near max[κ] apparent in the PSDs for ρu r φ (left) and for cos(πz/H)ρu r φ,z (middle) correspond to the trapped oscillations highlighted in the box in Fig. 8 (middle). The PSD calculated from spacetime data for ρu r φ,z (right) shows a more complicated spectrum of inertial-acoustic waves than was excited in the hydrodynamic simulations of Paper I.

Figure 11 .

 11 Figure 11. Plot showing the evolution of the Maxwell stress in simulations znAf10-40. The decreases in α M over time illustrate the weakening of MHD turbulence by eccentric disc distortion.

  in Ferreira & Ogilvie 2009, or fig. 2 in Paper I). Observations of superhumps in BHBs (O'Donoghue & Charles

Figure A1 .

 A1 Figure A1. Same as Fig. 5, but for the 2D, hydrodynamic simulations hAf10-40. The plot illustrates similar profiles of Ẽ to the MHD runs znAf10-40.

Figure A2 .Figure A3 .

 A2A3 Figure A2. Same as Fig. 6 (top), but for the 2D, hydrodynamic simulations hAf10-40. The snapshots show very similar distortions to those in znAf10-40.

Table A1 .

 A1 Table listing (i) simulation label, (ii) eccentricity forcing amplitude, (ii) saturated Reynolds stress, (iv) radial KE, (v) max | Ẽ| within 18r g , and (vi) max | Ẽ| between 6 -10r g . The values have been averaged over 30 -50T orb for hAf10-40, and 100 -150T orb for hlAf40.

	Label	A f	α R	ρu 2 r V / P V	max | Ẽ| max D | Ẽ|
	hAf10	0.10	0.0002	0.0089	0.011	0.010
	hAf30	0.30	0.0070	0.1298	0.043	0.027
	hAf35	0.35	0.0123	0.2217	0.057	0.031
	hAf37	0.375 0.0164	0.2933	0.065	0.033
	hAf40	0.40	0.0222	0.3944	0.075	0.034
	hlAf40	0.40	0.0107	0.1643	0.088	0.037

Freely available at https://sourcesup.renater.fr/projects/dumses/ MNRAS 000,1-15 (2020)

A non-linear theory is necessitated by the amplitude of forcing used in our simulations. Static sinusoidal profiles for density (associated with an m = 1 linear mode) also induce eccentricity, but we have found that density profiles and precession frequencies from more precise non-linear calculations produce a much cleaner quasi-steady state.MNRAS 000,1-15 (2020)

MNRAS 000,1-15 (2020)

Mixing between vertically structured r-modes and vertically unstructured f-modes can occur due to a breaking of symmetry with respect to the mid-plane by mean magnetic fields with both vertical and toroidal components(Dewberry et al. 2019).MNRAS 000,1-15 (2020)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000,1-15 (2020)

ACKNOWLEDGEMENTS

The authors thank the anonymous reviewer for comments and suggestions, which improved the quality of the paper. J. Dewberry thanks Adrian Barker, Roman Rafikov, Omer Blaes and Dong Lai for very helpful discussions. This work was funded by the Cambridge Commonwealth, European and International Trust, the Vassar College De Golier Trust, the Cambridge Philosophical Society, and the Tsung-Dao Lee Institute.

DATA AVAILABILITY

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