# Impact of convection and resistivity on angular momentum transport in dwarf novae

N. Scepi<fnref n="1" />, G. Lesur<fnref n="1" />, G. Dubus<fnref n="1" />, and M. Flock<fnref n="2" />

Received 5 September 2017 / Accepted 16 October 2017

## ABSTRACT

The eruptive cycles of dwarf novae are thought to be due to a thermal-viscous instability in the accretion disk surrounding the white dwarf. This model has long been known to imply enhanced angular momentum transport in the accretion disk during outburst. This is measured by the stress to pressure ratio $\alpha$, with $\alpha \approx 0.1$ required in outburst compared to $\alpha \approx 0.01$ in quiescence. Such an enhancement in $\alpha$ has recently been observed in simulations of turbulent transport driven by the magneto-rotational instability (MRI) when convection is present, without requiring a net magnetic flux. We independently recover this result by carrying out PLUTO magnetohydrodynamic (MHD) simulations of vertically stratified, radiative, shearing boxes with the thermodynamics and opacities appropriate to dwarf novae. The results are robust against the choice of vertical boundary conditions. The thermal equilibrium solutions found by the simulations trace the well-known S-curve in the density-temperature plane that constitutes the core of the disk thermal-viscous instability model. We confirm that the high values of $\alpha \approx 0.1$ occur near the tip of the hot branch of the S-curve, where convection is active. However, we also present thermally stable simulations at lower temperatures that have standard values of $\alpha \approx 0.03$ despite the presence of vigorous convection. We find no simple relationship between $\alpha$ and the strength of the convection, as measured by the ratio of convective to radiative flux. The cold branch is only very weakly ionized so, in the second part of this work, we studied the impact of non-ideal MHD effects on transport. Ohmic dissipation is the dominant effect in the conditions of quiescent dwarf novae. We include resistivity in the simulations and find that the MRI-driven transport is quenched ($\alpha \approx 0$) below the critical density at which the magnetic Reynolds number $R_m \le 10^4$. This is problematic because the X-ray emission observed in quiescent systems requires ongoing accretion onto the white dwarf. We verify that these X-rays cannot self-sustain MRI-driven turbulence by photo-ionizing the disk and discuss possible solutions to the issue of accretion in quiescence.

**Key words.** accretion, accretion disks – convection – magnetohydrodynamics (MHD) – turbulence – stars: dwarf novae

## 1. Introduction

A fundamental, yet challenging issue in accretion theory is the transport of angular momentum. Indeed, for matter to accrete it needs to transfer its angular momentum outward. Historically, the transport of angular momentum has been parametrized by the dimensionless parameter $\alpha$, the ratio of the fluid stress (responsible for the transport) to the local thermal pressure (Shakura & Sunyaev 1973). Thus, this “$\alpha$-prescription” reduces the physics to setting a phenomenological value for $\alpha$.

Dwarf novae (DNe) provide the best observational constrains on $\alpha$ (King et al. 2007). DNe are binary systems where matter is transferred by Roche lobe overflow from a solar-type star to a white dwarf. Their light curves show periodic outbursts during which the luminosity typically rises by several magnitudes (Warner 2003). According to the disk instability model (DIM; see Lasota 2001, for a review), these outbursts are caused by a thermal-viscous instability due to the steep temperature dependence of the opacity when hydrogen ionizes around 7000 K. During quiescence, the disk fills up as mass is transferred from the companion, gradually heating up until this critical temperature is reached at some radius. This triggers the propagation of heat fronts through the disk, bringing it into a hot state with a high accretion rate. The disk then empties until the temperature falls below the critical temperature and a cooling front returns the disk to quiescence. The timescales involved are set by

the ability of the disk to transport angular momentum, providing an observational handle on $\alpha$. Outburst-decay timescales imply $\alpha \sim 0.1$ (Kotko & Lasota 2012) whereas recurrence timescales imply $\alpha \sim 0.01$ in quiescence (Cannizzo et al. 1988, 2012). Although it has long been known that the DIM requires transport to be more efficient in outburst than in quiescence (Smak 1984), the physical reason for this change in $\alpha$ has remained elusive.

It is now widely accepted that angular momentum transport in disks is due to turbulence driven by the development of the magneto-rotational instability (MRI; Balbus & Hawley 1991). The properties of MRI-induced transport have been extensively studied over the past 25 years using local, shearing box simulations. Isothermal, stratified local simulations with zero net magnetic flux show a universal value of $\alpha \sim 0.03$ (Hawley et al. 1996; Simon et al. 2012), comparable to the value required for quiescent DNe. Higher values of $\alpha \sim 0.1$ require an external net vertical magnetic flux (Hawley et al. 1995), whose origin in DNe is unclear. Isothermal simulations are convenient idealisations but do not account for turbulent heating or radiative losses, and thus cannot be used to investigate the thermal equilibrium states of DNe. Latter & Papaloizou (2012) carried out more realistic simulations by using an analytic approximate local cooling law in a non-stratified shearing box, and showed numerically that the disk is indeed thermally unstable in the conditions of DNe. Their zero net flux simulations give $\alpha \sim 0.01$ in both cold and

<footnote n="1">¹ Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), 38000 Grenoble, France
e-mail: nicolas.scepi@gmail.com</footnote>

<footnote n="2">² Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California 91109, USA</footnote>
hot states. Hirose et al. (2014) performed the first simulations including radiative transfer, vertical stratification, and the realistic thermodynamics appropriate to DNe. They found that convection increases $\alpha$ to 0.1 in the hot state near the critical temperature, in the absence of net magnetic flux, providing a tantalizing solution to the change in $\alpha$ in DNe (Coleman et al. 2016).

Another explanation for the difference in transport efficiency between hot and cold states was proposed by Gammie & Menou (1998). In the quiescent state of DNe, the plasma is expected to be largely neutral and thus the magnetic field decouples from the disk. The conductivity drops as electrons become scarce and increasingly collide with neutrals. With the electron fraction being a strong function of the temperature, they pointed out that the MRI may not be able to grow or, at least, sustain fully developed turbulence in a quiescent DNe disk.

In light of these results, we have carried out numerical simulations (Sect. 2) to assess the impact of convection (Sect. 3) and resistivity (Sect. 4) on the transport of angular momentum in DNe in the hot, outburst- and cold quiescent states (respectively). In particular, we present the first local magnetohydrodynamics (MHD) shearing box simulations with realistic thermodynamics and radiative transfer that include Ohmic diffusion and conclude on the role that X-rays from the white dwarf boundary layer may play in the ionization balance of quiescent disks.

## 2. Methods

We adopt the local shearing-box approximation (Hawley et al. 1995) to simulate a vertically-stratified patch of accretion disk situated at a distance $R_0 = 1.315 \times 10^{10}$ cm from a $0.6 M_{\odot}$ white dwarf, giving an angular velocity $\Omega(R_0) = 5.931 \times 10^{-3}$ s$^{-1}$. These values are identical to those chosen by Hirose et al. (2014) to facilitate comparisons. The simulations include radiative transport in the flux diffusion approximation and thermodynamic quantities appropriate to the temperature and density regime sampled by DNe.

### 2.1. Basic equations

Curvature terms are not taken into account in the shearing box approximation and the differential Keplerian velocity is modeled as a linearized shear flow $v_y^0 = -(3/2) \Omega x$, where the $x$, $y$, and $z$ directions correspond to the radial, azimuthal, and vertical directions, respectively. The set of equations in the co-rotating frame is:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (1) $$

$$ \rho \frac{\partial \mathbf{v}}{\partial t} + (\rho \mathbf{v} \cdot \nabla) \mathbf{v} = -\nabla \left( P + \frac{B^2}{8\pi} \right) + \left( \frac{\mathbf{B}}{4\pi} \cdot \nabla \right) \mathbf{B} + \rho(-2\Omega\hat{\mathbf{z}} \times \mathbf{v}) + 3\Omega^2 x \hat{\mathbf{x}} - \Omega^2 z \hat{\mathbf{z}}, \qquad (2) $$

$$ \frac{\partial E}{\partial t} + \nabla \cdot [(E + P_t)\mathbf{v} - (\mathbf{v} \cdot \mathbf{B})\mathbf{B}] = -\rho\mathbf{v} \cdot \nabla\Phi - \kappa_P \rho c (a_R T^4 - E_R), \qquad (3) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{B} - \frac{4\pi}{c} \eta \mathbf{J}). \qquad (4) $$

The last three terms of Eq. (2) represent, respectively, the Coriolis force, the tidal force and the vertical component of the gravitational force; $\hat{\mathbf{x}}$ and $\hat{\mathbf{z}}$ are the units vectors in the $x$ and $z$ directions.

Radiative transfer is treated separately from the MHD step using an implicit time-stepping following Flock et al. (2013). In

this step, we solve the coupled matter-radiation equations in the flux-limited diffusion approximation

$$ \frac{\partial E_R}{\partial t} - \nabla \frac{c\lambda(R)}{\kappa_R \rho} \nabla E_R = \kappa_P \rho c (a_R T^4 - E_R), \qquad (5) $$

$$ \frac{\partial \epsilon}{\partial t} = -\kappa_P \rho c (a_R T^4 - E_R), \qquad (6) $$

where $\rho$ is the density, $\mathbf{v}$ the fluid velocity vector, $P$ the thermal pressure, $\mathbf{B}$ the magnetic field vector, $\Phi = (-3x^2 + z^2) \times \Omega^2(R_0)/2$ is the gravitational potential in the co-rotating frame, $\eta$ the Ohmic resistivity, $\mathbf{J} = (c/4\pi)\nabla \times \mathbf{B}$ the current density vector, $E = \rho\epsilon + 0.5\rho v^2 + \mathbf{B}^2/8\pi$ the total energy, $\epsilon$ the internal energy, $P_t = P + \mathbf{B}^2/8\pi$ the thermal pressure plus the magnetic pressure, $c$ the speed of light, $a_R = (4\sigma/c)$ the radiation constant with $\sigma$ the Stefan-Boltzmann constant, $T$ the temperature, $E_R$ the radiation density energy, $\kappa_P$ the Planck opacity and $\kappa_R$ the Rosseland opacity. The radiative energy flux, in the flux diffusion approximation, is $F_{\text{rad}} = (c\lambda(R)/\kappa_R(T)\rho)\nabla E_R$. The flux limiter is defined as $\lambda(R) \equiv (2+R)/(6+3R+R^2)$ with $R \equiv |\nabla E|/(\kappa_R \rho E)$ (Turner & Stone 2001). We do not take into account radiation pressure as it is negligible for the temperatures reached by DNe; Hirose et al. (2014) found that it contributes at most 6% of the gas+radiation pressure.

To close our set of equations we use the following equation of state (EOS) and internal energy function:

$$ P = \frac{\rho}{\mu(\rho, T)} k_B T, \qquad (7) $$

$$ \epsilon = \epsilon(\rho, T), \qquad (8) $$

where $\mu$ is the mean molecular weight. To compute $\mu$ and $\epsilon$, we use pre-calculated tables and interpolate linearly between values of the tables. Tables are computed from the Saha equations assuming ionization equilibrium for the solar composition of Grevesse & Sauval (1998) with a hydrogen abundance $X = 0.7$ and a metallicity $Z = 0.02$. The adiabatic index $\Gamma$, the entropy and the thermal capacity $C_v$ are similarly computed.

We use the opacity tables of Ferguson et al. (2005), which cover the low-temperature region from $2.7 < \log(T) < 4.5$, and OPAL (Iglesias & Rogers 1996), which cover the high-temperature region from $3.75 < \log(T) < 8.7$. We use a linear interpolation to connect the two and extend the resulting table, where necessary, using a zero-gradient extrapolation. Our tables of opacities and thermodynamic quantities are fully consistent with those used by Hirose et al. (2014).

### 2.2. Boundary conditions

We use shear-periodic conditions in the $x$-direction, periodic conditions in the $y$-direction, and either periodic or modified outflow conditions for the $z$-direction to test their influence on the results (Hirose et al. 2014, used only outflow conditions). Outflow conditions usually assume a zero-gradient extrapolation to the ghost cell if material is going outside of the box and prevent matter from outside from coming into the simulation box. Following Brandenburg et al. (1995), we also impose the magnetic field to point in the $z$-direction at the interface with the ghost cells. This circumvents a spurious increase in the magnetic field intensity that we observed, but did not further investigate, when using pure outflow conditions.

The mass in the shearing box may decrease due to the outflow boundaries. To avoid this, we normalize the total mass to the initial mass at each time step by multiplying the density by a
corrective factor. For the run 442O, the mass flux is of the order of ~0.1 g cm⁻² per orbital period. The time-averaged advective flux at the vertical boundary [⟨εv_z⟩](±z_max) remains negligible in outflow conditions (the averaging procedure is described in Sect. 2.5).

## 2.3. Numerical method

We solve the MHD equations on a 3D Cartesian grid with the conservative, Godunov-type code PLUTO (Mignone 2009). We choose a second order Runge-Kutta time integration method. Constrained transport ensures that $\nabla \cdot \mathbf{B} = 0$. The Riemann solver is HLLD except where the pressure difference between two adjacent cell exceeds five times the local pressure, in which case we allow for shock flattening by using the more diffusive solver HLL. To solve the radiative transfer equations, we follow the same implicit scheme as Flock et al. (2013) except that we use the bi-conjugate gradient solver KSPIBCGS (Yang & Brent 2002) as implemented in PETSC (Balay et al. 2016), which we find to provide stabler, faster convergence than BiCGSTAB for our application.

In order to avoid very small time steps, we use floors of 10⁻⁶ and 5 × 10⁻² of the initial mid-plane values for the density and the temperature, respectively. We find that floors are activated over ≲2% of the box and, thus, have no significant impact on the thermodynamical equilibrium of the simulations.

## 2.4. Initial conditions

The initial flow is Keplerian with a weak zero net flux magnetic field of the form $B_z \propto \sin(2\pi x)$. We usually start with an isothermal layer and assume hydrostatic equilibrium to fix the initial vertical density profile. We let the MRI develop and then trigger radiative transfer after 48 local orbits (time is normalized by the quantity 1/Ω₀ thus one orbital period is equivalent to 2π in code units). This time is sufficient for the MRI to saturate and for the disk to reach a quasi-steady state. The isothermal temperature $T_{c0}$ is set to the mid-plane temperature found by using the vertical structure code of Hameury et al. (1998) for given surface density Σ₀ and effective temperature $T_{eff}$. This code solves the vertical structure equations assuming an α-prescription for the angular momentum transport and associated heating rate, radiative transfer in the diffusion approximation, and convection described by mixing length theory with a mixing coefficient $\alpha_{ml} = 1.5$ (based on models of the Sun). We then activate radiative transfer and let the disk equilibrate and reach a quasi-steady state (if there is one).

Some runs are initialized from another simulation by changing the surface density manually instead of starting from an isothermal layer. This has proven useful close to the hydrogen ionization regime, where the disk easily undergoes critical heating/cooling. In these cases, starting from an isothermal state can cause the disk to miss the equilibrium state. Starting from a nearby quasi steady-state allows for a smoother transition to capture the thermal equilibrium.

## 2.5. Runs and diagnostics

Table 1 lists the runs that we performed. We adopt the notation of Hirose et al. (2014) for labeling the runs, with the addition of a final O to indicate simulations with vertical outflow boundary conditions or a final P for a periodic vertical boundary. Σ₀ is the initial surface density and H ≡ c_s(T_{c0})/Ω is the pressure

scale-height (with c_s the sound speed). The horizontal extent of the box is of ±6H for the hot branch and ±3H for the cold and middle branch. L_x, L_y and L_z follow the ratio 1:2:4 (see Table 1).

The horizontal average of a quantity f is defined as:

$$ \langle f \rangle_{x,y} (z,t) = \frac{\iint f(x, y, z, t) \, dx \, dy}{\iint dx \, dy}, \qquad (9) $$

where the integration in x and y is done over the whole simulation box. We use ⟨⟩ brackets to denote spatial averages and [] brackets to indicate an average over a time t_{avg} performed once the simulation has reached a quasi-steady state. The averaging timescale t_{avg} is indicated for each run in the last column of Table 1. We typically average over a hundred orbits except for some cases where the variability occurs on long timescales (e.g., 401O, see Sect. 3.2.3). [Σ], [T_{mid}], [T_{eff}], and [τ_{tot}] are thus the time-averaged values of the surface density, midplane temperature, effective temperature and total optical depth, respectively. σ_{T_{mid}}, σ_{T_{eff}} and σ_{α} are the standard deviations of the fluctuations with time of these quantities.

We compute τ_{tot}, Σ, T_{eff} as

$$ \tau_{\text{tot}} = \int \langle \rho \rangle_{x,y} \langle \kappa_R \rangle_{x,y} \, dz, \qquad (10) $$

$$ \Sigma = \int \langle \rho \rangle_{x,y} \, dz, \qquad (11) $$

and

$$ T_{\text{eff}} = \left( \frac{1}{\sigma_B} (F_{\text{rad } z}^{+} - F_{\text{rad } z}^{-}) \right)^{1/4}, \qquad (12) $$

where $F_{rad z}^{+/-}$ is the radiative flux in the vertical direction on the upper/lower boundary of our simulation box. We also define ⟨Q⁻⟩_{x,y} the total cooling rate as

$$ \langle Q^{-} \rangle_{x,y} \equiv \frac{d}{dz} \langle F_{\text{rad } z} \rangle_{x,y} + \frac{d}{dz} \langle \epsilon v_z \rangle_{x,y}, \qquad (13) $$

$v_z$ the vertical velocity and the quantity $\epsilon v_z$ represents the advective flux of internal energy $F_{adv}$. The integration over z extends over the full box.

We define α, the ratio of stress to pressure, as

$$ \alpha \equiv \frac{[\langle W_{xy} \rangle_{x,y,z}]}{[\langle P \rangle_{x,y,z}]}, \qquad (14) $$

and the instantaneous stress to pressure ratio as $\tilde{\alpha} = \langle W_{xy} \rangle_{x,y,z} / \langle P \rangle_{x,y,z}$. The turbulent stress $W_{xy}$ is $\rho(v_x v_y - v_{Ax} v_{Ay})$ where $v_A$ is the Alfvén velocity.

## 2.6. Numerical convergence

We take 32 points in the x direction, 128 points in the y direction and 256 points in the z direction in all our simulations. This number of points must be sufficient to resolve the most unstable wavelength of the MRI, which is of order H (as the root mean square turbulent magnetic field is of the order of 0.1–0.2 for all simulations), while still allowing for a significant vertical extension. Hence, we use different values of L_x, L_y and L_z on the cold branch and on the hot branch. We note n the number of points per final time-averaged scale height h ≡ c_s([T_{mid}])/Ω. On the high Σ part of the hot branch n is of the order of 20 (Table 1), enough to resolve the MRI. But at lower Σ some of our simulations are
<table><caption>**Table 1.** Initial parameters and results for our ideal MHD simulations.</caption><thead><tr><th>Run</th><th>Restart</th><th>Σ<sub>0</sub></th><th>T<sub>c0</sub></th><th>H/R<sub>0</sub></th><th>[T<sub>mid</sub>] ± σT<sub>mid</sub></th><th>T<sub>midmax</sub></th><th>T<sub>midmin</sub></th><th>[T<sub>eff</sub>] ± σT<sub>eff</sub></th><th>α ± σα</th><th>[τ<sub>tot</sub>]</th><th>n</th><th>L<sub>x</sub>/H</th><th>L<sub>y</sub>/H</th><th>L<sub>z</sub>/H</th><th>t<sub>avg</sub></th></tr></thead><tbody><tr><td colspan="16" style="text-align:center; font-weight:bold;">Hot branch</td></tr><tr><td>432P</td><td>–</td><td>4000</td><td>23 1194</td><td>5.62E-02</td><td>21 6384 ± 7862</td><td>23 8779</td><td>20 5844</td><td>27 729 ± 3249</td><td>0.066 ± 0.014</td><td>39 294</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>16</td></tr><tr><td>430P</td><td>–</td><td>2900</td><td>18 4762</td><td>5.03E-02</td><td>16 8546 ± 9001</td><td>18 8999</td><td>15 1736</td><td>22 584 ± 2218</td><td>0.047 ± 0.014</td><td>34 026</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>429P</td><td>–</td><td>1030</td><td>10 4765</td><td>3.78E-02</td><td>98 395 ± 4544</td><td>11 0410</td><td>91 551</td><td>14 793 ± 1617</td><td>0.042 ± 0.015</td><td>20 450</td><td>20</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>431P</td><td>–</td><td>750</td><td>91 431</td><td>3.54E-02</td><td>86 437 ± 3314</td><td>95 825</td><td>79 883</td><td>13 268 ± 1070</td><td>0.041 ± 0.010</td><td>18 443</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>80</td></tr><tr><td>431O</td><td>–</td><td>750</td><td>91 431</td><td>3.54E-02</td><td>87 500 ± 3097</td><td>96 656</td><td>83 558</td><td>13 128 ± 899</td><td>0.042 ± 0.008</td><td>17 520</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>439P</td><td>–</td><td>540</td><td>79 831</td><td>3.30E-02</td><td>72 093 ± 3953</td><td>83 749</td><td>67 827</td><td>10 718 ± 1296</td><td>0.031 ± 0.014</td><td>22 519</td><td>20</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>439O</td><td>–</td><td>540</td><td>79 831</td><td>3.30E-02</td><td>74 220 ± 3082</td><td>81 313</td><td>68 523</td><td>10 920 ± 969</td><td>0.033 ± 0.010</td><td>19 578</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>468P</td><td>–</td><td>386</td><td>70 561</td><td>3.11E-02</td><td>66 985 ± 1882</td><td>72 396</td><td>62 772</td><td>10 493 ± 791</td><td>0.042 ± 0.009</td><td>16 169</td><td>22</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>468O</td><td>–</td><td>386</td><td>70 561</td><td>3.11E-02</td><td>66 784 ± 3186</td><td>76 966</td><td>62 428</td><td>10 281 ± 1046</td><td>0.042 ± 0.014</td><td>16 755</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>470P</td><td>–</td><td>275</td><td>63 720</td><td>2.95E-02</td><td>56 203 ± 2773</td><td>62 911</td><td>52 180</td><td>8842 ± 1003</td><td>0.036 ± 0.013</td><td>18 288</td><td>20</td><td>1.5</td><td>6</td><td>12</td><td>48</td></tr><tr><td>437P</td><td>–</td><td>174</td><td>54 545</td><td>2.73E-02</td><td>44 630 ± 4281</td><td>52 848</td><td>36 837</td><td>8117 ± 633</td><td>0.046 ± 0.017</td><td>15 323</td><td>19</td><td>1.5</td><td>6</td><td>12</td><td>60</td></tr><tr><td>437O</td><td>–</td><td>174</td><td>54 545</td><td>2.73E-02</td><td>47 414 ± 2788</td><td>52 911</td><td>40 886</td><td>8245 ± 499</td><td>0.052 ± 0.011</td><td>12 375</td><td>20</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>441P</td><td>437P</td><td>134</td><td>–</td><td>2.73E-02</td><td>42 882 ± 3757</td><td>51 820</td><td>34 140</td><td>8161 ± 633</td><td>0.070 ± 0.021</td><td>10 344</td><td>19</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>441O</td><td>437O</td><td>134</td><td>–</td><td>2.73E-02</td><td>35 605 ± 7841</td><td>51 140</td><td>18 170</td><td>7795 ± 511</td><td>0.061 ± 0.022</td><td>19 797</td><td>17</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>446P</td><td>437P</td><td>127</td><td>–</td><td>2.73E-02</td><td>30 515 ± 9137</td><td>46 504</td><td>15 721</td><td>7563 ± 421</td><td>0.066 ± 0.027</td><td>25 335</td><td>16</td><td>1.5</td><td>6</td><td>12</td><td>92</td></tr><tr><td>446O</td><td>437O</td><td>127</td><td>–</td><td>2.73E-02</td><td>34 881 ± 5947</td><td>42 975</td><td>17 104</td><td>7567 ± 330</td><td>0.061 ± 0.021</td><td>16 285</td><td>17</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>442P</td><td>437P</td><td>113</td><td>–</td><td>2.73E-02</td><td>35 004 ± 6857</td><td>45 903</td><td>17 805</td><td>7745 ± 419</td><td>0.079 ± 0.031</td><td>12 796</td><td>17</td><td>1.5</td><td>6</td><td>12</td><td>83</td></tr><tr><td>442O</td><td>–</td><td>113</td><td>34 500</td><td>2.17E-02</td><td>32 724 ± 4878</td><td>43 105</td><td>20 887</td><td>7292 ± 323</td><td>0.075 ± 0.022</td><td>15 318</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>452P</td><td>437P</td><td>100</td><td>–</td><td>2.73E-02</td><td>28 719 ± 9205</td><td>41 063</td><td>14 888</td><td>7302 ± 435</td><td>0.092 ± 0.024</td><td>18 108</td><td>15</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>452O</td><td>442O</td><td>102</td><td>–</td><td>2.17E-02</td><td>33 245 ± 5735</td><td>44 029</td><td>19 624</td><td>7390 ± 348</td><td>0.087 ± 0.031</td><td>11 514</td><td>21</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>453P</td><td>437P</td><td>90</td><td>–</td><td>2.73E-02</td><td>35 292 ± 3918</td><td>42 119</td><td>23 611</td><td>7332 ± 627</td><td>0.098 ± 0.039</td><td>7357</td><td>17</td><td>1.5</td><td>6</td><td>12</td><td>84</td></tr><tr><td>453O</td><td>442O</td><td>90</td><td>–</td><td>2.17E-02</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td colspan="16" style="text-align:center; font-weight:bold;">Middle branch</td></tr><tr><td>438P</td><td>–</td><td>275</td><td>14 000</td><td>1.38E-02</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>1.5</td><td>6</td><td>12</td><td>48</td></tr><tr><td>428O</td><td>–</td><td>229</td><td>13 000</td><td>1.33E-02</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>1.5</td><td>6</td><td>12</td><td>95</td></tr><tr><td>403O</td><td>–</td><td>220</td><td>13 000</td><td>1.33E-02</td><td>7893 ± 1261</td><td>9956</td><td>5362</td><td>4369 ± 162</td><td>0.027 ± 0.010</td><td>2050</td><td>33</td><td>0.75</td><td>3</td><td>6</td><td>191</td></tr><tr><td>402O</td><td>–</td><td>210</td><td>13 000</td><td>1.33E-02</td><td>8085 ± 1370</td><td>10 515</td><td>5016</td><td>4423 ± 141</td><td>0.029 ± 0.011</td><td>2560</td><td>35</td><td>0.75</td><td>3</td><td>6</td><td>318</td></tr><tr><td>401O</td><td>–</td><td>200</td><td>12 000</td><td>1.28E-02</td><td>9431 ± 1224</td><td>11 964</td><td>6348</td><td>4514 ± 187</td><td>0.036 ± 0.012</td><td>8458</td><td>38</td><td>0.75</td><td>3</td><td>6</td><td>350</td></tr><tr><td>404O</td><td>–</td><td>150</td><td>9000</td><td>1.11E-02</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>R</td><td>0.75</td><td>3</td><td>6</td><td>95</td></tr><tr><td colspan="16" style="text-align:center; font-weight:bold;">Cold branch</td></tr><tr><td>435P</td><td>–</td><td>191</td><td>2645</td><td>6.01E-03</td><td>3783 ± 103</td><td>3983</td><td>3559</td><td>3206 ± 153</td><td>0.036 ± 0.009</td><td>16</td><td>15</td><td>0.75</td><td>3</td><td>6</td><td>95</td></tr><tr><td>435O</td><td>–</td><td>191</td><td>4000</td><td>7.39E-03</td><td>3785 ± 127</td><td>4055</td><td>3593</td><td>3334 ± 213</td><td>0.036 ± 0.011</td><td>15</td><td>42</td><td>0.75</td><td>3</td><td>6</td><td>95</td></tr><tr><td>462P</td><td>–</td><td>174</td><td>2785</td><td>6.17E-03</td><td>3102 ± 262</td><td>3471</td><td>2643</td><td>2855 ± 179</td><td>0.042 ± 0.014</td><td>8</td><td>45</td><td>0.75</td><td>3</td><td>6</td><td>95</td></tr><tr><td>465P</td><td>462P</td><td>116</td><td>–</td><td>6.17E-03</td><td>2095 ± 115</td><td>2346</td><td>1922</td><td>2332 ± 115</td><td>0.035 ± 0.008</td><td>2</td><td>39</td><td>0.75</td><td>3</td><td>6</td><td>95</td></tr><tr><td>434P</td><td>–</td><td>93</td><td>1976</td><td>5.20E-03</td><td>1897 ± 49</td><td>2028</td><td>1814</td><td>2115 ± 109</td><td>0.037 ± 0.010</td><td>1</td><td>42</td><td>0.75</td><td>3</td><td>6</td><td>95</td></tr><tr><td>434O</td><td>–</td><td>93</td><td>1976</td><td>5.20E-03</td><td>1861 ± 33</td><td>1966</td><td>1814</td><td>2022 ± 108</td><td>0.034 ± 0.010</td><td>1</td><td>42</td><td>0.75</td><td>3</td><td>6</td><td>95</td></tr><tr><td>476P</td><td>–</td><td>45</td><td>1828</td><td>4.99E-03</td><td>1715 ± 47</td><td>1840</td><td>1604</td><td>1733 ± 116</td><td>0.037 ± 0.011</td><td>4</td><td>41</td><td>0.75</td><td>3</td><td>6</td><td>56</td></tr></tbody></table>

**Notes.** $\Sigma_0$ is the initial surface density in $g \text{ cm}^{-2}$, $T_{c0}$ the initial midplane temperature in K, and $H/R_0$ where $H$ is the corresponding scale height. Brackets [ ] are for quantities averaged over $t_{\text{avg}}$ (given in local orbits) with $\sigma$ the associated standard deviation. $L_x/H$, $L_y/H$ and $L_z/H$ are the number of initial scale heights in the box in the x, y, z directions. n is the number of points per time-averaged scale height. The second column indicates the simulation from which the run was initialized, otherwise we started from an isothermal layer at $T_{c0}$. P denotes a run with periodic vertical boundary condition; O denotes outflow conditions. R signals simulations with runaway heating or cooling.

restarted from previous runs with higher temperatures, decreasing $n$ to 15 in run 452P. We checked those cases with additional runs with a more appropriate box size and resolution and found that the results are quantitatively the same. For instance, we find a difference of ~5% in $T_{\text{mid}}$ and ~2% in $\alpha$ between the run 442O started with an isothermal layer (with $n = 21$) or initialized from run 437O (with $n = 13$).

High Alfvénic velocities are allowed near the boundaries when the box height is large, about $12H$, increasing the numerical diffusion due to the HLLD solver and providing a numerical source of heating. On the hot branch, the contribution of this numerical coronal heating with respect to the total numerical heating is negligible: we performed the simulation 437O with $L_x/H = 1.5$, $L_y/H = 6.0$ and $L_z/H = 8.0$ and found a difference on $[T_{\text{mid}}]$ of only 2% for the hot branch. The effect is much stronger on the middle and cold branch, heating being weaker, leading us to take smaller box sizes. Consequently, $n$ is approximately two times larger in the cold branch than on the hot branch.

According to Ryan et al. (2017), isothermal stratified simulations are not fully resolved since $\alpha$ decreases as $n^{-1/3}$ up to the highest achievable resolutions. Therefore, a caveat when comparing runs on the cold and hot branch is that the two-fold

increase in resolution may lower the value of $\alpha$ by ~20% in the cold branch simulations. We ran the simulation 434O with half the resolution on each direction and found $\alpha = 0.019$. This difference is of the same order as the estimated error on $\alpha$. Hence we cannot conclude about consistency with Ryan et al. (2017). A proper convergence study would be useful to know if this result holds for our simulations including realistic thermodynamics and radiative transfer.

Finally, Simon et al. (2012) showed that the box size $L_x$ and $L_y$ must be $\ge 2H$ to have converged diagnostics of the MRI. Again, these results were obtained using an isothermal equation of state and we do not know for sure if they hold here. According to this criterion, our box size is appropriate on the hot branch but may be too small in $x$ on the cold branch (Table 1). In particular, convective eddies in our simulations could be limited by small box sizes. We doubled the size in $x$ and $y$ of the run 442O to check for this effect and found no notable difference.

## 3. Ideal MHD runs

### 3.1. Thermal equilibrium curve

Figure 1 shows the thermal equilibria reached by the simulations listed in Table 1. The top panel shows the time-averaged
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Thermal equilibria in the [$T_{\text{mid}}$] vs. [$\Sigma$] (top) or [$T_{\text{eff}}$] vs. [$\Sigma$] (bottom) plane. Squares and circles are for periodic and outflow runs, respectively. Triangular dots represent runs with runaway cooling (triangle facing down) or heating (triangle facing up). Error bars represent the standard deviation of the temperature fluctuations. The symbols are color-coded to the value of $\alpha$. The color-coded curves are vertical thermal equilibria using an $\alpha$-prescription. The dashed blue line indicates where the magnetic Reynolds number $R_m = 10^4$ based on an isothermal model (see Sect. 4.1).</figcaption></figure>

midplane temperature [$T_{\text{mid}}$] $\equiv [\langle T(z=0) \rangle_{x,y}]$ versus the surface density [$\Sigma$] (top panel). The bottom panel shows [$T_{\text{eff}}$] versus [$\Sigma$] (bottom panel). The intensity of $\alpha$ is indicated for each run by the color of the dot. The colored curves in Fig. 1 are thermal equilibria curves obtained by using the Hameury et al. (1998) code, which assumes an $\alpha$-prescription and allows for thermal convection (Sect. 2.4). Whereas the latter requires choosing a value for $\alpha$, the simulations do not have this degree-of-freedom since the value of $\alpha$ is set by the turbulent angular momentum transport generated by the MRI.

The simulations trace an equilibrium thermal curve in the shape of an S (an S-curve) with a hot stable branch and a cold

stable branch, providing independent confirmation of the results of Hirose et al. (2014). The S-curve from the simulations is comparable to that predicted by the vertical structure calculations using an $\alpha$-prescription. These also predict a third stable branch at intermediate temperatures ($T_{\text{mid}} \approx 10^4$ K). Indeed, Fig. 1 shows also a set of three stable runs in an intermediate regime of temperature. A middle-branch in radiative MRI simulations was also reported by Hirose (2015) in the context of protoplanetary disks. The middle branch is not as extended as the other branches and may be seen as a prolongation of the cold branch to higher temperatures. In fact, the middle branch does not overlap with the cold branch, so there is at most two equilibrium states for
a given Σ, not three. However, the middle branch runs have a higher opacity than the cold branch and are convectively unstable (Sect. 3.2).

The hot and cold/middle branches are terminated at, respectively, low and high Σ by unstable runs displaying runaway cooling or heating (noted R in Table 1 and indicated by arrows in Fig. 1). We have followed these runaways to the stable equilibrium on the opposite branch using appropriate box sizes. Periodic simulations seem to be slightly more stable around points of critical heating/cooling as 453O is unstable whereas 453P is not. Restarting simulations from previous runs by changing the surface density allows us to follow the hot branch to lower surface densities than in Hirose et al. (2014) where 442O is their last stable run.

We have good agreement between our simulations and the α-prescription calculation concerning the equilibrium temperatures at given α and the critical points for runaway heating/cooling. The equilibrium values follow the calculated curves corresponding to the value of α found in the simulation. Somewhat unexpectedly, the agreement extends to the middle branch even though its location in the (Σ, T) plane depends on the chosen value of the convection mixing length parameter (Cannizzo & Wheeler 1984). Here, we took α<sub>ml</sub> = 1.5 as in the solar convection zone (Sect. 2.4).

## 3.2. Convection and transport of angular momentum

### 3.2.1. Enhancement of α

We see an enhancement of α in the low Σ part of the hot branch, albeit with slightly lower maximum values of α ≈ 0.098 than the α ≈ 0.121 found by Hirose et al. (2014). On the cold and middle branch, we find values of α ranging from 0.029 to 0.042 typical of a zero net vertical flux stratified MRI simulation (Simon et al. 2012). The values of α are comparable on the hot branch for 174 g cm⁻² ≤ Σ ≤ 1030 g cm⁻². For higher Σ the value of α increases up to 0.066; we did not further investigate this trend as this part of the S-curve implies very high accretion rates that are unlikely to be relevant to dwarf novae. The value of α also increases as we approach the tip of the hot branch, going to a maximum of 0.098 for the lowest-density run on the hot branch (453P). This enhancement is clearly seen when plotting α against temperature (Fig. 2, see also Fig. 1 in Coleman et al. 2016). The enhancement is accompanied by stronger fluctuations as we near the unstable tip of the hot branch (see Sect. 3.2.3).

### 3.2.2. Convective transport of heat

Hirose et al. (2014) attributed the enhanced α to convection. We also observe that the enhanced α runs are convectively unstable and further note that the value of α does not depend on the chosen vertical boundary conditions. We measure the stability of a run against convection from the Brunt-Väisälä frequency N, defined by

$$ \frac{N^2}{\Omega^2} = \frac{d \ln([\langle P \rangle]^{1/\langle \Gamma \rangle} [\langle \rho \rangle])}{d \ln z}, \quad (15) $$

where Γ ≡ (∂ ln P/∂ ln ρ)s is the adiabatic index and s is the specific entropy. N² has the sign of the entropy gradient along z, thus a negative Brunt-Väisälä frequency denotes a convectively unstable vertical profile. The Brunt-Väisälä frequency gives a linear stability criterion for the convective stability only in regions where the total pressure is dominated by thermal pressure and does not give direct information about the convection flux.

<figure><img src="image_2.png" /><figcaption>Fig. 2. α as a function of [T<sub>mid</sub>] (top panel) or [T<sub>eff</sub>] (bottom panel). Error bars represent the standard deviations of the fluctuations around the mean values. Blue, green, and red colors are for the cold, middle and hot branches, respectively. Squares and circles indicate periodic and outflow runs. The shaded area corresponds to the thermally unstable region.</figcaption></figure>

Figure 3 shows the radiative, advective, and total fluxes (Eq. (13)) as well as the Brunt-Väisälä frequency as a function of height for several illustrative cases. The top panel shows a highly convective episode (averaging over 16 local orbits) of the run 442P (hot branch). The advective flux is the main source of heat transport in the regions where the disk is convectively unstable (N² < 0) as opposed to the upper regions where the radiative flux dominates. Averaging the same run on a longer time (80 local orbits) including convective and non-convective episodes shows the contribution of the advective flux decreases (middle panel). We also see evidence for downward transport of heat between the midplane and H/R₀ ~ 2.3 × 10⁻². This downward transport is also present in the top panel but is less extended in height. Convection is expected to drive heat upward, where the entropy is lower, hence this downward transport is not due to convective motions. In fact, run 434F on the cold branch is convectively stable yet also shows this downward advective transport (bottom panel). We attribute this transport to vertical mixing by the MRI-driven turbulence. Because motions are approximately adiabatic in this region, this mixing tends to flatten the entropy profile, resulting in a downward transport of heat. Thus, the advective
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Time-averaged vertical profiles of the radiative, advective, and total flux for different runs. **Top panel:** run 442P during a highly convective episode. **Middle panel:** same run but averaging over convective and non convective episodes. **Bottom panel:** run 434P where there is no convection. The vertical line indicates the height above which magnetic pressure becomes larger than thermal pressure. The dashed line shows $N^2/\Omega^2$ where N is the Brunt-Väisälä frequency.</figcaption></figure>

flux includes contributions from both convective motions (when present) and turbulent mixing.

Besides the region of enhanced $\alpha$ near the unstable tip of the hot branch, we find convection also plays a major role in the middle branch runs, notably run 441O, with the fluxes behaving as in the top panel of Fig. 3.

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Evolution of $\tilde{\alpha}$ as a function of time for runs located on different parts of the S-curve: (top panel) 434P is a cold branch run, 468P a non-convective hot branch run, and (bottom panel) 442 a convective hot branch run. For the latter we show the simulation run with outflow (442O) and periodic boundary conditions (442P). We also plot, for run 442O, the contribution of the Maxwell stress to $\tilde{\alpha}$ as a dashed line.</figcaption></figure>

### 3.2.3. Convective/radiative cycles

In the following, we make use of the convective fraction $f_{conv}$ of the flux, defined as

$$f_{conv} = \frac{1}{2} \left( \max \left( 0, \max \left( \left( \frac{F_{adv}(z>0)}{F_{rad z}} \right) \right) \right) + \max \left( 0, \max \left( \left( \frac{F_{adv}(z<0)}{F_{rad z}} \right) \right) \right) \right). \quad (16)$$

We also use an instantaneous convective fraction $\tilde{f}_{conv}$ obtained by smoothing over a time window with a width of 1.6 local orbits. We note that $\tilde{f}_{conv}$ can have values >1 due to the fluctuations remaining in the vertical profile.

Figure 4 shows the evolution of $\tilde{\alpha}$ with time for runs 434P, 468P, 442P and 442O. The first two are convectively stable and show only weak fluctuations in $\tilde{\alpha}$. However, run 442 is located near the tip of the hot branch and is convectively unstable. This run shows $\tilde{\alpha}$ has large fluctuations, going through cycles with bursts of strong angular momentum transport. We see from Fig. 4 that the main contribution to $\tilde{\alpha}$ is the Maxwell part. The maximum and the recurrence of these bursts in $\tilde{\alpha}$ increase towards the unstable tip of the S-curve. For comparison, there is one cycle every ~20 orbits in run 437O ($\Sigma = 174 \text{ g cm}^{-2}$) and one every ~10 orbits in run 442O ($\Sigma = 113 \text{ g cm}^{-2}$).
<figure><img src="image_5.png" /><figcaption>Fig. 5. Autocorrelation function $f_{auto}$ as a function of time for the run 4420. $\tau_0$ is the parabolic decay time of the autocorrelation function (with the parabolic fit shown as a dashed line). $\tau_1$ is the location of the first peak of $f_{auto}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6. Period $\tau_1$ of the cycles in $\alpha$ as a function of the convective fraction for convectively unstable simulations on the hot branch.</figcaption></figure>

To be more quantitative about the frequency of the cycles of $\tilde{\alpha}$, we compute its autocorrelation function

$$f_{\text{auto}}(\tau) = \frac{\int (\tilde{\alpha}(t) - \alpha)(\tilde{\alpha}(t + \tau) - \alpha)dt}{\int (\tilde{\alpha}(t) - \alpha)^2dt} \quad (17)$$

Figure 5 shows this function for run 4420. We extract two characteristic times. The first is the parabolic decay time $\tau_0$ of $f_{auto}$ ($\tau_0 \approx 2$ local orbits in Fig. 5); it corresponds to the time on which small fluctuations of $\tilde{\alpha}$ remain correlated. The second, $\tau_1$, corresponds to the first peak of $f_{auto}$ ($\tau_1 \approx 9$ local orbits in Fig. 5); this corresponds roughly to the period of the convective cycles.

We do not see any clear evolution of $\tau_0$ as the convective fraction increases on the hot branch. However, $\tau_1$ clearly decreases with $f_{conv}$ (Fig. 6). On the middle branch, the cycles are much longer for the same $f_{conv}$. For instance, run 4010 has $f_{conv} = 0.14$ and $\tau_1 \approx 111$ local orbits, requiring a longer averaging timescale $t_{avg}$ than for the other runs.

The stronger variability near the unstable tip of the hot branch is also manifest as increased temperature fluctuations (error bars in Figs. 1 and 2). These fluctuations in temperature are anti-correlated with $\tilde{f}_{conv}$ on the hot branch (Fig. 7, top panel). However, the fluctuations are in phase at the highest temperatures on the middle branch (Fig. 7, bottom panel). The behavior on the hot and middle branch is very different although the fractional amplitudes ([$\max$ - min]/[$\max$ + min]) of the fluctuations

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. $\tilde{f}_{conv}$ vs. $T_{mid}$ for two representative convectively unstable simulations: hot branch run 4420 (top) and middle branch run 4010 (bottom).</figcaption></figure>

in $\tilde{f}_{conv}$ and $T_{mid}$ are the same for both. This difference can be explained by noting that the temperature gradient, assuming only radiative transport, is $\propto K/T^3$, and that convection requires this gradient to be greater than the adiabatic gradient of the gas. On the hot branch, the opacity decreases when the temperature rises, flattening the temperature profile and quenching convection. Upward fluctuations in temperature lead the disk from a convective regime to a radiative regime and vice versa. The opposite correlation is expected on the middle branch since the opacity rises steeply with $T$, triggering a convective regime when fluctuations raise the temperature.

This difference in behavior between the convective parts of the hot branch and the middle branch also extends to $\alpha$. Run 4420 on the hot branch shows a trend for the instantaneous $\tilde{\alpha}$ to increase with the instantaneous convective fraction $\tilde{f}_{conv}$ (Fig. 8, green markers), but this trend is absent in run 4010 on the middle branch (Fig. 8, blue markers). Furthermore, the maximum value of $\tilde{\alpha}$ is lower than on the hot branch run even though $\tilde{f}_{conv}$ reaches higher values. The runs also show a difference in duty cycle. The hot branch run spends roughly equal time in the low and high $\tilde{f}_{conv}$ parts of the diagram. In contrast, the middle branch run spends most of its time in the low $\tilde{f}_{conv}$ part of the diagram, with only temporary excursions to the higher values. This will naturally lead to a smaller time-averaged $\alpha$ in the middle branch. Figure 9 shows no trend in the convective runs between the time-averaged values of $\alpha$ and $f_{conv}$. We conclude that there is no simple relationship between the enhancement of $\alpha$ and the relative strength of the convection measured by $f_{conv}$. Finding the physical origin for this enhancement is likely to require a
<figure><img src="image_9.png" /><figcaption>Fig. 8. $\tilde{\alpha}$ vs. $\tilde{f}_{conv}$ for the same runs as Fig. 7.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. $\alpha$ vs. convective fraction $f_{conv}$ for all convectively unstable runs (same color and symbol coding as Fig. 2).</figcaption></figure>

detailed study of the impact of vertical convection on the MRI dynamo, which is beyond the scope of this paper.

### 3.3. Impact of the vertical boundary conditions

We tested outflow and periodic vertical boundary conditions (Sect. 2.4) because of their potential impact on the interplay between the MRI dynamo and convection. In periodic runs, the box-averaged $\langle B_y \rangle_{x,y,z} = 0$ by construction whereas it fluctuates over time with outflow boundary conditions (Fig. 10). The flips in sign are related to the dynamo of stratified MRI turbulence and are still poorly understood (Brandenburg et al. 1995). We found no quantitative difference between periodic and outflow runs so these flips have no impact on $\alpha$ or on the thermal equilibrium values.

The similar behavior between outflow and periodic runs extends to the “butterfly” diagram that emerges from maps of the time evolution of $\langle B_y \rangle_{x,y}$, the vertical profile of the azimuthal field. Non-convective MRI simulations are known to form “butterfly wings” as the magnetic field flips sign at each quasi-periodic pulsation and propagates outward (e.g. see Fig. 12). Coleman et al. (2017) reported that the quasi-periodic pulsation does not always lead to a sign flip of $\langle B_y \rangle_{x,y}$ when convection is active. We confirm that this is also the case in our convective outflow simulations (Fig. 11, top panel) and in our periodic boundary simulations (bottom panel), even though the magnetic field must adopt a symmetric configuration in the periodic simulations

<figure><img src="image_11.png" /><figcaption>Fig. 10. Average azimuthal magnetic field $\langle B_y \rangle_{x,y,z}$ as a function of time for run 446 (outflow O or periodic P boundary conditions).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Vertical profile of azimuthal magnetic field $\langle B_y \rangle_{x,y}$ (in code units) as a function of time with outflow (top) and periodic (bottom) boundary conditions for run 446 (hot branch).</figcaption></figure>

to ensure that $\langle B_y \rangle_{x,y,z} = 0$. We also see this behavior in the middle branch during convective episodes.

### 4. Resistive MHD runs

Ideal MHD may not apply to the cold branch due to the very low ionization fractions. The development of the MRI can be suppressed as the resistivity due to electron-neutral collisions increase. This is quantified by the magnetic Reynolds number

$$R_m \equiv \frac{c_s h}{\eta}, \qquad (18)$$

with $\eta$ being the resistivity, $c_s$ the sound speed and $h = c_s/\Omega$ the local pressure scale-height. For $R_m < 10^4$, it is expected that diffusion of the magnetic field becomes too important for the disk to sustain MHD turbulence (Hawley et al. 1996; Fleming et al. 2000). Gammie & Menou (1998) pointed out that $R_m$ is of the order of $10^3$ in a quiescent DNe disk.
<figure><img src="image_13.png" /><figcaption>Fig. 12. Time evolution of $\langle R_m \rangle$, $\langle B_y \rangle$ and $\alpha$ for run 435PR ($\Sigma = 191 \text{ g cm}^{-2}$). The vertical black line marks the time at which the simulation is restarted from 435P with resistivity turned on.</figcaption></figure>

We investigate whether or not our cold branch runs can maintain turbulence when resistivity is included. The resistivity is computed as in Blaes & Balbus (1994):

$$ \eta = 230 \left( \frac{n_n}{n_e} \right) T^{1/2} \text{ cm}^2 \text{ s}^{-1}, \qquad (19) $$

with the assumption that we are dominated by electron-neutral collision. The values of $\eta$ are pre-computed with the equation of state (Saha equilibrium) and saved in a table. The time step $\Delta t$ associated with the diffusive term is of the order of $\Delta x^2/\eta$ where $\Delta x$ is the cell minimum length. To avoid dramatically small time steps, we impose a minimum value of $R_m = 50$ in our runs. This floor is largely under the limit for MRI turbulence to be suppressed and does not affect our results. We find $\Delta t \sim 10^{-3}$ to be of the order of the time step of the ideal MHD run. The resistive runs are restarted from the periodic runs in ideal MHD and labeled PR to differentiate them.

<caption>Fig. 13. Same as Fig. 12 for run 462PR ($\Sigma = 174 \text{ g cm}^{-2}$).</caption>

## 4.1. Turbulence decay due to resistivity

We find that the resistivity has a critical impact by suppressing turbulence on the cold branch below some critical density located between $\Sigma = 174 \text{ g cm}^{-2}$ (run 462P) and $\Sigma = 191 \text{ g cm}^{-2}$ (run 435). Figures 12, 13 show the temporal evolution of $\tilde{\alpha}$, $\langle R_m \rangle_{x,y}$, and $\langle B_y \rangle_{x,y}$ for those two runs before and after resistivity is switched on.

Turning on resistivity has no influence on run 435PR (Fig. 12). $R_m$ stays well above the critical value of $10^4$ in the midplane. Some regions above the photosphere have a lower $R_m$, which is without consequences as most of the heating due to MRI happens in the densest part of the disk.

Run 462P has a lower $\Sigma$ and temperature ($[T_{\text{mid}}] \approx 3100 \text{ K}$ compared to 3750 K for 435P). Here, $R_m$ starts to drop as soon as resistivity is switched on (Fig. 13). The turbulence is suppressed once $R_m$ decreases below $\approx 5000$ and the transport of
angular momentum ceases. $\langle B_y \rangle_{x,y}$ diffuses as expected: first near the photosphere and then, as $R_m$ decreases, in the whole box.

## 4.2. Influence of ambipolar diffusion and Hall effect

Let us first compare the order of magnitude of Ohmic (O), ambipolar (A) and the Hall effect (H). From Balbus & Terquem (2001), we have

$$ \frac{\mathrm{O}}{\mathrm{H}} = \left( \frac{\rho}{3 \times 10^{-6} \mathrm{g} \mathrm{cm}^{-3}} \right)^{1/2} \left( \frac{c_s}{v_A} \right), \qquad (20) $$

$$ \frac{\mathrm{O}}{\mathrm{A}} = \left( \frac{\rho}{10^{-8} \mathrm{g} \mathrm{cm}^{-3}} \right) \left( \frac{T}{10^3 \mathrm{K}} \right)^{-1/2} \left( \frac{c_s}{v_A} \right)^2. \qquad (21) $$

We take 462PR as a fiducial run (as it is the first one to be altered by Ohmic diffusion) to compare the importance of each non-ideal effect. The typical values deduced from this run for the density, temperature, ionization fraction and magnetisation $v_A/c_s$ are $5 \times 10^{-7} \mathrm{g} \mathrm{cm}^{-3}$, 3000 K, $10^{-5.5}$ and 0.1, respectively. Therefore, in this run, ambipolar diffusion is expected to have a negligible impact while the Hall effect is only an order of magnitude weaker than Ohmic diffusion. We note that these conclusions also hold as the disc cools down when MRI is suppressed.

However, comparing the amplitude of each non-ideal effect might not be enough to conclude on the dynamics since they affect the MRI in different ways. Let us therefore analyze the impact of ambipolar diffusion and of the Hall effect individually.

Ambipolar diffusion is characterized by the ambipolar Elsasser number (Balbus & Terquem 2001):

$$ \Lambda_A \equiv \frac{\gamma_i \rho_i}{\Omega}, \qquad (22) $$

where $\rho_i = n_i m_i$ is the ion density and $\gamma_i = \langle \sigma v \rangle_{ni} / (m_n + m_i) = 2.7 \times 10^{13} (41.33 m_u / m_n + m_i) \mathrm{cm}^3 \mathrm{g}^{-1} \mathrm{s}^{-1}$ is the ion-neutral drag coefficient. Following Balbus & Terquem (2001), we assume a momentum rate coefficient for ion-neutral collision $\langle \sigma v \rangle_{ni} = 1.9 \times 10^{-9} \mathrm{cm}^3 \mathrm{s}^{-1}$ (Draine 2010), $m_n = 2.33 m_u$ and $m_i = 39 m_u$ (the neutral and ion masses respectively). We obtain $\Lambda_A \sim 10^3$ for our typical values. Hawley & Stone (1998) and Bai & Stone (2011) found that angular momentum transport is not significantly impacted when $\Lambda_A > 100$. Their result is based on numerical simulations with a net magnetic flux and should still hold when there is no net flux: a net magnetic flux implies currents outside of our domain of simulation that are not affected by resistivity and help sustain MRI turbulence. We conclude that the influence of ambipolar diffusion would have on run 462PR is negligible compared to Ohmic diffusion.

The Hall effect is characterized by the Hall Lundquist number

$$ \mathcal{L}_H \equiv \sqrt{\frac{4\pi n_e e H}{\rho c}}. \qquad (23) $$

The impact of the Hall effect on the MRI saturation level is twofold. For strong Hall effect ($\mathcal{L}_H \lesssim 5$), turbulent transport is essentially suppressed and the flow relaxes into a self-organized state (Kunz & Lesur 2013). For weaker Hall effect ($\mathcal{L}_H \gtrsim 10$) the system produces sustained turbulent transport provided that the magnetic Reynolds number is sufficiently large, similarly to the pure Ohmic case. This critical Reynolds number for sustained turbulence is reduced to $10^3$ for $\mathcal{L}_H \approx 30$ while it is identical to the pure Ohmic case for $\mathcal{L}_H \gtrsim 60$ (Sano & Stone 2002), indicating that the Hall effect becomes negligible above this threshold.

Our typical values imply $\mathcal{L}_H \sim 10^3$, so the Hall effect does not affect the MRI close to the region where Ohmic diffusion becomes problematic. As the discs cools down in the MRI stable region, $\mathcal{L}_H$ decreases, making the Hall effect stronger and potentially reviving the MRI. However, the ratio $\mathcal{L}_H/R_m$ does not depend on the temperature, meaning that at the temperature where $\mathcal{L}_H \lesssim 50$, one also expects $R_m \lesssim 50$, well below the critical $R_m$ found in Hall-MRI simulation (Sano & Stone 2002). We conclude that the Hall effect does not affect the critical temperature (or $\Sigma$) below which DNe are MRI-stable, nor does it revive the MRI at lower temperatures.

## 4.3. Influence of X-ray irradiation

The ionization fraction $x_e$ is critical to determine whether or not the coupling between fluid and magnetic field is strong enough for MRI turbulence to be active. In the resistive runs described above (Sect. 4.1) $x_e \approx 4 \times 10^{-6}$ in the MRI-stable run 462PR (at midplane) and $x_e \approx 6 \times 10^{-6}$ in the MRI-unstable run 435PR. Therefore, a very small ionization fraction of order $5 \times 10^{-6}$ can be sufficient to maintain turbulence. These ionization fractions are thermal and do not take into account external sources of ionizations. In particular, quiescent DNe show hard X-ray emission originating from the accretion boundary layer onto the white dwarf that may provide sufficient ionization to maintain the MRI active in the quiescent disk.

The quiescent X-ray spectra are bremsstrahlung with temperatures $kT_{XR}$ ranging from 1 to 10 keV and luminosities $L_{XR}$ ranging from $10^{28}$ to $10^{32} \mathrm{erg} \mathrm{s}^{-1}$ (Byckling et al. 2010), properties akin to those of a protostar. The typical surface densities under consideration are similar to those of a protostellar disk at 1 AU but the DNe disk is $\sim 10^3$ closer to the X-ray source, so the impinging ionizing flux is much higher.

We calculate the expected ionization fraction $x_e$ by following studies of protoplanetary disks. We first compute the ionization rate $\zeta$ as in Glassgold et al. (1997a,b), assuming a photon flux

$$ f_0 = C \frac{L_{XR}}{4\pi R_0^2 k T_{XR}}, \qquad (24) $$

where C is a factor taking into account the irradiation geometry and the albedo A of the disk. The X-ray emission region is about the size of the white dwarf (Mukai 2017), which is much greater than the disk height. In this case, we may expect (King 1997)

$$ C \approx \frac{2}{3\pi} \left( \frac{R_{WD}}{R_0} \right) (1 - A) \approx 10^{-2} (1 - A), \qquad (25) $$

with a white dwarf radius $R_{WD} = 10^9 \mathrm{cm}$.

We obtain $x_e$ by solving Eq. (11) from Fromang et al. (2002), which relates $x_e$ to the fraction of metals $x_M = n_M/n_n$, the dissociative recombination rate for molecular ions, the radiative recombination rate for metal atoms and the rate of charge exchange between them. We further assume that the disk is passive, that is, has no internal heating, and the only source of heat is X-ray irradiation. Hence, in steady-state, the disk is isothermal at the temperature set by $\sigma T_{irr}^4 = kT_{XR}f_0$. The vertical structure is in hydrostatic equilibrium with $\Sigma = 100 \mathrm{g} \mathrm{cm}^{-2}$ to compare to our results with the first MRI-stable run, 462PR.

We find that X-ray ionization is negligible, in agreement with the order-of-magnitude estimate of Gammie & Menou (1998). For example, with $CL_{XR} = 10^{30} \mathrm{erg} \mathrm{s}^{-1}$ and $kT_{XR} = 10 \mathrm{keV}$, we find $\zeta = 3.4 \times 10^{-9} \mathrm{s}^{-1}$ at the disk midplane. The isothermal atmosphere has $T = 1700 \mathrm{K}$ and a midplane density
$n = 3.7 \times 10^{17} \text{ cm}^{-3}$, giving $x_e = 4.1 \times 10^{-10}$ when no metals are present ($x_M = 0$) and $x_e = 1.1 \times 10^{-7}$ when $x_M = 6.86 \times 10^{-6}$ (solar abundance). The fraction $x_e$ decreases for lower values of $CL_{XR}$ or $kT_{XR}$. Hence, X-ray ionization is unable to provide the required ionization fraction except in the upper layer. However, only a small amount of material is active in this layer: taking the values assumed above, $x_e > 5 \times 10^{-6}$ for $z > 1.5 \times 10^8 \text{ cm}$ (about $2H$) and the column density of the active layer represents less than 3% of $\Sigma$.

## 5. Conclusions

We have performed stratified, radiative, ideal and resistive MHD, shearing box simulations in conditions appropriate to DNe. We find that the thermal equilibrium solutions found by the simulations trace the well-known S-curve derived from $\alpha$-prescription models, including a middle branch that extends the cold branch to higher temperatures and is characterized by vigorous convection. We confirm that $\alpha$ increases to $\approx 0.1$ near the unstable tip of the hot branch as reported by Hirose et al. (2014). This increase is thus robust against the choice of numerical code and, as we investigated, against the choice of outflow or periodic vertical boundary conditions. Although convection plays a major role in transporting heat in the runs with an enhanced $\alpha$, we find no clear relationship between $\alpha$ and $f_{conv}$, the average fraction of the flux carried by convection. Notably, $\alpha$ is not enhanced on the middle branch although it is strongly convectively unstable. Detailed studies of the interplay between convection and the MRI dynamo will be required to understand the exact mechanism behind the enhancement of $\alpha$ towards the tip of the hot branch.

Ohmic dissipation can prevent MRI-driven turbulent transport in the cold, quiescent state of DNe (Gammie & Menou 1998). We verified that Ohmic dissipation is the dominant non-ideal effect and carried out resistive MHD simulations to test its influence on $\alpha$. We find that the region of the cold branch with $\Sigma < 191 \text{ g cm}^{-1}$ does not maintain the MRI-driven turbulence. Fleming et al. (2000) found in isothermal simulations that the presence of a net vertical magnetic flux, providing external support to the MRI, sustains turbulence down to $R_m = 100$. Net flux may help push the MRI-inactive region to slightly lower $\Sigma$ but the dynamical consequences might be difficult to evaluate without a global model to study the evolution of the net flux. A net magnetic field is also expected to impact substantially the saturated value of $\alpha$ when the corresponding $\beta \le 10^5$ (Bai & Stone 2013; Salvesen et al. 2016), that is, $B \ge 20 \text{ G}$ on the hot branch at a $\Sigma$ of $174 \text{ g cm}^{-2}$ and $B \ge 4 \text{ G}$ on the cold branch at a $\Sigma$ of $93 \text{ g cm}^{-2}$. Using a dipole approximation, this corresponds to a surface magnetic field of $\sim 10^{-5} \text{ G}$ for the white dwarf, smaller than the typical field of intermediate polars (in which the disk is truncated by the white dwarf magnetosphere), or $10^{-3} \text{ G}$ for the companion, somewhat on the high side of the measured magnetic fields in low-mass stars (see Fig. 3 of Donati & Landstreet 2009). Hence, in principle, the binary components might be able to provide enough net magnetic flux to impact transport even though, in practice, how this is achieved is unknown.

A quiescent disk will not accrete if it enters the MRI-inactive region. This is problematic because the X-ray emission observed in quiescent DNe requires ongoing accretion onto the white dwarf. The X-ray emission region is constrained to the immediate vicinity of the white dwarf by eclipses, as expected from the accretion boundary layer (see Mukai 2017, and references therein). We verified that this emission is unable to self-sustain MRI-driven accretion by photo-ionizing the disk.

If MRI remains inactive, angular momentum transport may be due to an entirely different mechanism in quiescence. Spiral shocks are often invoked as a means to accrete. They have been proposed to explain the observed anti-correlation between outburst recurrence time and binary mass ratio in some subtypes of dwarf novae (Cannizzo et al. 2012; Menou 2000). However, other parameters, such as the mass-transfer rate from the secondary, also depend on the mass ratio and therefore impact accretion, rendering the interpretation of this anti-correlation difficult. Besides, global simulations studying the propagation of spiral shocks in hydrodynamics (Savonije et al. 1994; Lesur et al. 2015; Arzamasskiy & Rafikov 2017) and MHD (Ju et al. 2017) showed that when the ratio $H/R$ is of the order of $10^{-2}$, as is the case in the cold state (especially in the outer region which is the region of interest as spiral shocks propagate inward), spiral shock dissipate rapidly and do not lead to any accretion in the inner parts of the disk.

The problem may turn out to be a red herring if the disk actually never samples the regime where MRI is inactive. Our results show that accretion is possible on the high-$\Sigma$ part of the cold branch, and this may be enough for the DNe cycle to operate. DIM models of DNe outburst cycles with an $\alpha$-prescription show disk rings follow complex tracks in the $(\Sigma, T)$ plane when a cooling front propagates through. The ring does not simply cool on a thermal timescale at constant $\Sigma$ from the tip of the hot branch to the cold branch of the S-curve. Radial heat fluxes cause it to follow a sideways trajectory to a higher $\Sigma$ on the cold branch (see e.g., Fig. 9 in Menou et al. 1999), and this would ensure it never samples the MRI-inactive part of the cold branch. These radial fluxes become important only because the cooling-front size is of order of the vertical height of the disk, a situation that strains the assumption in these calculations that the radial and vertical directions are decoupled. Consequently, there are some uncertainties on how they are taken into account by the DIM and their influence on the light curves (Cannizzo et al. 2010); uncertainties that only a global model of DNe outbursts under MRI-driven transport could resolve.

More prominently, such a global model would also address whether realistic light curves of DNe outbursts can be obtained at all with MRI transport. In this regard, the enhancement of $\alpha$ up to the values required by the DIM is both a major advance; a set back. A set back because the high values of $\alpha$ are limited to the tip of the hot branch, whereas DIM models assumed $\alpha$ to be increased over the whole hot branch. The difference has important consequences on the light curves predicted by the DIM when $\alpha(T)$ is taken from the MRI simulations: Coleman et al. (2016) obtained very jagged light curves that, in our opinion, compare less favorably to the observations than those of the initial DIM (see their Fig. 9). Yet, this is also undisputedly a major advance because this at last provides a physical explanation of the phenomenological change in $\alpha$ between hot and cold state, first noticed over 30 yr ago; and because, after a certain lull, interest in dwarf novae, the source of much of the basic understanding of accretion physics, has been re-awakened.

**Acknowledgements.** We thank the referee for valuable comments and suggestions. N.S. acknowledges financial support from the pole PAGE of the Université Grenoble Alpes. G.L., G.D. and M.F. are grateful to the participants of the KITP 2017 program on Confronting MHD Theories of Accretion Disks with Observations for the many useful conversations pertaining to this work (and thus this research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915). This work was granted access to the HPC resources of IDRIS under the allocation A0020402231 made by GENCI (Grand Equipement National de Calcul Intensif). Some of the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the
Rhône-Alpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.

## References

Arzamasskiy, L., & Rafikov, R. R. 2017, ApJ, submitted, ArXiv e-prints [arXiv:1710.01304]

Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144

Bai, X.-N., & Stone, J. M. 2013, ApJ, 767, 30

Balay, S., Abhyankar, S., Adams, M. F., et al. 2016, PETSc Web page, http://www.mcs.anl.gov/petsc

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235

Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163

Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741

Byckling, K., Mukai, K., Thorstensen, J., & Osborne, J. 2010, MNRAS, 408, 2298

Cannizzo, J. K., & Wheeler, J. C. 1984, ApJS, 55, 367

Cannizzo, J. K., Shafter, A. W., & Wheeler, J. C. 1988, ApJ, 333, 227

Cannizzo, J. K., Still, M. D., Howell, S. B., Wood, M. A., & Smale, A. P. 2010, ApJ, 725, 1393

Cannizzo, J. K., Smale, A. P., Wood, M. A., Still, M. D., & Howell, S. B. 2012, ApJ, 747, 117

Coleman, M. S. B., Kotko, I., Blaes, O., Lasota, J.-P., & Hirose, S. 2016, MNRAS, 462, 3710

Coleman, M. S., Yerger, E., Blaes, O., Salvesen, G., & Hirose, S. 2017, MNRAS, 467, 2625

Donati, J.-F., & Landstreet, J. 2009, ARA&A, 47, 333

Draine, B. T. 2010, Physics of the interstellar and intergalactic medium (Princeton University Press)

Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585

Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464

Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43

Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18

Gammie, C. F., & Menou, K. 1998, ApJ, 492, L75

Glassgold, A. E., Najita, J., & Igea, J. 1997a, ApJ, 480, 344

Glassgold, A. E., Najita, J., & Igea, J. 1997b, ApJ, 485, 920

Grevesse, N., & Sauval, A. 1998, Space Sci. Rev., 85, 161

Hameury, J.-M., Menou, K., Dubus, G., Lasota, J.-P., & Hure, J.-M. 1998, MNRAS, 298, 1048

Hawley, J. F., & Stone, J. M. 1998, ApJ, 501, 758

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690

Hirose, S. 2015, MNRAS, 448, 3105

Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., & Sano, T. 2014, ApJ, 787, 1

Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943

Ju, W., Stone, J. M., & Zhu, Z. 2017, ApJ, 841, 29

King, A. R. 1997, MNRAS, 288, L16

King, A., Pringle, J., & Livio, M. 2007, MNRAS, 376, 1740

Kotko, I., & Lasota, J.-P. 2012, A&A, 545, A115

Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295

Lasota, J.-P. 2001, New Astron. Rev., 45, 449

Latter, H. N., & Papaloizou, J. C. 2012, MNRAS, 426, 1107

Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9

Menou, K. 2000, Science, 288, 2022

Menou, K., Hameury, J.-M., & Stehle, R. 1999, MNRAS, 305, 79

Mignone, A. 2009, Mem. Soc. Astron. Ital. Suppl., 13, 67

Mukai, K. 2017, PASP, 129, 062001

Ryan, B. R., Gammie, C. F., Fromang, S., & Kestener, P. 2017, ApJ, 840, 6

Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016, MNRAS, 457, 857

Sano, T., & Stone, J. M. 2002, ApJ, 577, 534

Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685

Smak, J. 1984, Acta Astron., 34, 161

Turner, N., & Stone, J. 2001, ApJS, 135, 95

Warner, B. 2003, Cataclysmic variable stars (Cambridge University Press), 28

Yang, L. T., & Brent, R. P. 2002, Proc. of the Fifth International Conference on Algorithms and Architectures for Parallel Processing, IEEE