© ESO 2018

1. Introduction

The missing link between accretion theory and observations is the mechanism transporting angular momentum. Turbulent transport is one candidate and is often parameterized by the dimensionless parameter α, ratio of the fluid stress to the local thermal pressure (Shakura & Sunyaev 1973). Theory of thin viscous α-disks has been widely applied to various systems such as protoplanetary disks (Sano et al. 2000; Papaloizou & Terquem 1999), soft X-ray transients (SXT; King & Ritter 1998; Tetarenko et al. 2018), dwarf novae (DNe; Smak 1984), or even active galactic nuclei (AGN; Starling et al. 2004; see, however, Hameury et al. 2009). However, in soft X-ray transient sources and dwarf novae, which are the focus of this paper, the best estimators of α are often found because their nature is time dependent (King et al. 2007).

DNe are compact binary systems in which matter is transferred by Roche-lobe overflow from a solar-type star to a white dwarf. DNe have been observed over decades, and their light curves show periodic outbursts with amplitudes of 2–3 magnitudes (Warner 2003). Their eruptive behavior is well explained by the disk instability model (DIM; Lasota 2001; Dubus et al. 2018), in which a thermal instability caused by the ionization of hydrogen is responsible for a hysteresis cycle. In this system, the hot or cold state has a higher or lower accretion rate than the mass accretion rate coming from the companion, which causes the disk to fill up or empty. This leads to an hysteresis cycle in luminosity. The timescales involved are related to the ability for the disk to efficiently extract its angular momentum, providing an observational handle on α. Outburst-decay timescales imply α  ≈  0.1 in the hot state (Kotko & Lasota 2012), and recurrence timescales gives α  ≈  0.01 in the cold state (Cannizzo et al. 1988, 2012). It is unclear why this change in α occurs and if it is related to the properties of turbulence. Nevertheless, we can use these observational constraints as a benchmark to distinguish between different transport mechanisms for DNe.

Consensus appears to be reached that the magneto-rotational instability (MRI; Balbus & Hawley 1991) is the main driver of turbulent motions as it requires only a subthermal magnetic field threading the disk and a radially decreasing angular velocity; these are quasi-omnipresent conditions in astrophysical disks. However, magneto-hydrodynamical processes require coupling between gas and magnetic field, and this can be problematic in cold regions of protoplanetary disks (Gammie 1996), but also for dwarf novae in quiescence (Gammie & Menou 1998; Scepi et al. 2018, hereafter S18).

Given the complexity of MRI saturation, numerical simulations have been extensively used to study the characteristics of the resulting turbulence and transport. The properties of MRI-driven turbulence are known to depend on the magnetic configuration, especially the amplitude of the local net magnetic field (Hawley et al. 1995), which is largely unconstrained in DNe disks. A self-sustaining version of the MRI turbulence also exists in the absence of an externally imposed magnetic field, known as zero net-flux configuration (ZNF), and it is known to be independent of the initial magnetic seed (Hawley et al. 1996). Hence, ZNF simulations provide a minimum level of angular momentum transport by the MRI that does not require a large-scale magnetic field. We add the caveat that the convergence of ZNF simulations is still debated (see Ryan et al. 2017 for more information on the subject).

Isothermal simulations of MRI with ZNF are known to give a universal value of α  ≈  0.01 (Hawley et al. 1996; Simon et al. 2012) that is characteristic of the quiescent state of DNe. Using an approximate function of cooling, Latter & Papaloizou (2012) were able to retrieve thermal equilibrium curves of DNe from ZNF MRI simulations, showing the two expected stable branches, the hot and the cold state, which correspond to the eruptive and the quiescent states, respectively. However, their simulations exhibit a value of α  ≈  0.01 regardless of the branch. By including vertical stratification, a realistic equation of state and radiative transfer in the regime of DNe, Hirose et al. (2014) found for the first time that ZNF MRI can lead to α  ≈  0.1 in the low-density part of the hot branch. The actual interpretation is that there exists a non-linear coupling with convection which enhances α in this regime. This result was confirmed in S18 with a different code and in Coleman et al. (2018) in a different regime of opacities. However, this enhanced α at the tip of the hot branch fails to reproduce the light curves of DNe by inducing reflares, which is inconsistent with observations (Coleman et al. 2016). Moreover, S18 found that in ZNF the cold branch was too resistive to sustain magnetic turbulence, leading to α  ≈  0 in quiescence, as suggested by Gammie & Menou (1998).

Hence, the simple application of ZNF MRI to the case of DNe leads to inconsistencies with the DIM model and the observations. ZNF is a common assumption since the magnetic configuration is unknown in DNe. Moreover, large-scale magnetic fields leads to higher Alfvén speeds in the atmosphere and thus are more computationally demanding. They also require a robust solver when the magnetization is high. It is well known, however, that α depends on the ratio of the thermal to magnetic pressure, β (Hawley et al. 1996), which can lead to α ≳ 0.1 if the disk is threaded by a strong enough large-scale magnetic field. Additionally, it is expected that a large magnetic field will sustain MRI turbulence to lower ionization levels than in ZNF (Fleming et al. 2000). This readmits MRI as a possible candidate on the cold branch. The purpose of this paper is twofold: to study the local thermal equilibrium of a disk in realistic conditions of DNe threaded by a large-scale magnetic field, and to determine whether this field could resolve the discrepancies between MRI transport and observations.

2. Methods

We adopted the local shearing-box approximation (Hawley et al. 1995) to simulate a vertically stratified patch of an accretion disk located at a distance R₀ = 1.315 × 10¹⁰ cm from a 0.6 M_⊙ white dwarf, giving an angular velocity Ω(R₀)=5.931 × 10⁻³ s⁻¹ to facilitate comparisons with Hirose et al. (2014). The simulations include radiative transport in the flux-limited 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 linear shear flow $v_y^0=-(3/2)\Omega x$, where x, y, and z correspond to the radial, azimuthal, and vertical directions, respectively. The set of equations in the corotating frame is

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·(\rho \mathit{v})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \rho \frac{\partial \mathit{v}}{\partial t}+(\rho \mathit{v}·\mathbf{\nabla })\mathit{v}=-\mathbf{\nabla }(P+\frac{B^2}{8\pi })+(\frac{\mathbf{B}}{4\pi }·\mathbf{\nabla })\mathbf{B}+\rho (-2\mathrm{\Omega }\widehat{\mathit{z}}\times \mathit{v}\hfill \\ \hfill & +3{\mathrm{\Omega }}^{2}x\widehat{\mathit{x}}-{\mathrm{\Omega }}^{2}z\widehat{\mathit{z}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial E}{\partial t}+\mathbf{\nabla }·[(E+P_t)\mathit{v}-(\mathit{v}·\mathbf{B})\mathbf{B}]=-\rho \mathit{v}·\mathbf{\nabla }\mathrm{\Phi }-{\kappa }_{\mathrm{P}}\rho c(a_R{T}^4-E_{\mathrm{R}}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathbf{B}}{\partial t}=\mathbf{\nabla }\times (\mathit{v}\times \mathbf{B}-\frac{4\pi }{c}\eta \mathbf{J}).\hfill \end{array}$$(4)

The last three terms of Eq. (1) represent the Coriolis force, the tidal force, and the vertical component of the gravitational force, respectively; $\widehat{\mathit{x}}$ and $\widehat{\mathit{z}}$ are the unit vectors in the x and z directions.

Radiative transfer is treated separately from the MHD step using an implicit time-stepping following the implementation of Flock et al. (2013). In this step, we solve the coupled matter-radiation equations in the flux-limited diffusion approximation

$$\begin{array}{cc}& \frac{\partial E_{\mathrm{R}}}{\partial t}-\mathbf{\nabla }\frac{c\lambda (R)}{{\kappa }_{\mathrm{R}}\rho }\mathbf{\nabla }{E}_{\mathrm{R}}={\kappa }_{\mathrm{P}}\rho c(a_R{T}^4-E_{\mathrm{R}})\hfill \\ \hfill & \frac{\partial ϵ}{\partial t}=-{\kappa }_{\mathrm{P}}\rho c(a_R{T}^4-E_{\mathrm{R}}),\hfill \end{array}$$

where ρ is the density, v the fluid velocity vector, P the thermal pressure, B the magnetic field vector, Φ = Ω²(R₀)( − 3x² + z²)/2 is the gravitational potential in the corotating frame, η the Ohmic resistivity, J = (c/4π)∇ × B the current density vector, E = ϵ + 0.5ρv² + B^2/8π the total energy, ϵ the internal energy, P_t = P + B^2/8π the thermal pressure plus the magnetic pressure, c the speed of light, a_R = (4σ/c) the radiation constant with σ the Stefan–Boltzmann constant, T the temperature, E_R the radiation density energy, κ_P the Planck opacity, and κ_R the Rosseland opacity. The radiative energy flux in the flux-diffusion approximation is F_rad = (cλ(R)/κ_R(T)ρ)∇E_R. The flux limiter is defined as λ(R)≡(2 + R)/(6 + 3R + R²) with R ≡ |∇E|/(κ_RρE) (Turner & Stone 2001). We performed two simulations with the limiter from Minerbo (1978) and found no notable differences. We did not take into account radiation pressure as it is negligible for the temperatures reached by DNe; Hirose et al. (2014) found that it contributes 6% of the gas+radiation pressure at most.

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

$$\begin{array}{cc}& P=\frac{\rho }{\mu (\rho ,T)}{k}_{\mathrm{B}}T\hfill \\ \hfill & ϵ=ϵ(\rho ,T),\hfill \end{array}$$

where μ is the mean molecular weight. To compute μ and ϵ, we used precalculated tables and interpolated linearly between the table values. Tables were 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 Γ, the entropy, and the thermal capacity C_v were computed similarly.

We used 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 covers the high-temperature region from 3.75 <  log(T)< 8.7. We used a linear interpolation to connect the two and extended the resulting table, where necessary, using a zero-gradient extrapolation.

2.2. Boundary conditions

We used shear-periodic boundary conditions in the x-direction and periodic boundary conditions in the y-direction. We showed in S18 that the choice of vertical boundary conditions only slightly change the final thermal equilibrium with zero net-flux. Since periodic boundary conditions are much more computationally consuming than modified outflow boundaries, we chose to use the latter. The modified outflow boundaries (Brandenburg et al. 1995) assume a zero-gradient extrapolation for all hydrodynamic quantities leaving the box and prevent matter from outside from entering the simulation box. The difference with pure outflow conditions is that they impose the magnetic field to be vertical at the z-boundary. We did not investigate other boundary conditions.

The mass in the shearing box may decrease as a result of the outflow boundaries. To avoid this, we normalized the total mass to the initial mass at each time step by multiplying the density by a corrective factor. For simulations with β  ≈  10³, outflows become important. In the simulation 439F with β  ≈  10³, the disk would be empty after seven orbits if no normalization were used. When β  ≈  10², the disk expands considerably because of magnetic pressure and the location of the photosphere is outside of the box. When losses of matter and energy are important, we normalized the pressure by the same factor as for the density to account for the energy loss of the disk through vertical boundaries. When the pressure was normalized, we observed a change of 40% in the midplane temperature for β  ≈  10² compared to only 7% when β  ≈  10³. We ran two simulations with β as low as ≈10² and a normalized pressure, but tried to avoid this configuration otherwise.

2.3. Numerical method

We solved the MHD equations on a 3D Cartesian grid with the conservative Godunov-type code PLUTO (Mignone 2009). We chose a second-order Runge-Kutta time-integration method. Constrained transport ensures that ∇ ⋅ B = 0 up to machine precision. We used the HLLD solver to solve the Riemann problem. The HLLD solver can be used with a general EOS where Γ is not a constant. However, when we computed the maximum and minimum propagation speed waves, we used the Davis estimate (Davis 1988) and assumed that Γ is constant on each side of the propagating waves. The Riemann solver is HLLD except where the pressure difference between two adjacent cells exceeds five times the local pressure, the local ratio of thermal pressure over magnetic pressure is lower than 10⁻³ or the Alfvénic speed is greater than 30, in which case we used the more diffusive solver HLL. In our simulations each cell spent less than 10% of its time in HLL. To solve the radiative transfer equations, we followed the same implicit scheme as Flock et al. (2013), except that we used the biconjugate 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 used floors of 10⁻⁶ and 5 × 10⁻² of the initial mid-plane values for the density and temperature, respectively.

2.4. Initial conditions

The initial flow was Keplerian with a starting constant vertical magnetic field B_z 0 of variable intensity. We chose the intensity of the magnetic fields to be 8 G, 2 G, or 0.8 G. This is consistent with a white dwarf with a radius of 10⁹ cm and a dipolar surface magnetic field of 20 000, 5000, and 2000 G, respectively, or a stellar companion of ≈1 R_⊙ with a dipolar surface magnetic field of 60, 15, and 6 G, respectively, at a binary separation of 10¹¹ cm. The choice of our boundary conditions in x and y ensures that the net flux of magnetic field crossing the equatorial plane is constant throughout our simulation.

We started with an isothermal layer and assumed hydrostatic equilibrium to fix the initial vertical density profile. We let the MRI develop, then triggered radiative transfer after 32 local orbits (time was normalized by the quantity 1/Ω₀, thus one orbital period is equivalent to 2π in code units) and let the disk equilibrate and reach a quasi-steady state (if there was one).

For simulations with β ≥ 10⁵, the isothermal temperature T_c0 was set to be approximately the mid-plane temperature found using the vertical structure code of Hameury et al. (1998) for a given surface density Σ₀ and a given effective temperature T_eff, as in the case of ZNF in S18. This code solves the vertical structure equations assuming an α-prescription for the angular momentum transport and associated heating rate. It also uses radiative transfer in the diffusion approximation, and convection described by mixing length theory with a mixing coefficient α_ml = 1.5 (based on models of the Sun). However, when β <  10⁵, to avoid inappropriate box sizes that are due to the influence of the magnetic field, we evaluated the expected mid-plane temperature from test simulations.

2.5. Runs and diagnostics

Table 1 lists the runs that we performed. We partly adopted the notation of Hirose et al. (2014) to label the runs. Σ₀ 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 was ±6H for the hot branch and ±3H for the cold and middle branch. L_x,  L_y, and L_z follow the ratio 1:4:8 on the hot branch and 1:4:4 on the cold branch. The resolution was 32 × 128 × 256, except for test simulations. We chose to use smaller vertical boxes on the cold branch to avoid high Alfvén velocities near the boundaries, resulting in an enhanced numerical diffusion and thus numerical heating.

In the following, we use several averaged quantities denoted in the following way:

$$\begin{array}{cc}& {⟨X⟩}_{\rho }\equiv \frac{1}{\mathrm{\Sigma }{L}_x{L}_y}{\displaystyle {\int }_{-L_z/2}^{L_z/2}{\int }_{-L_y/2}^{L_y/2}{\int }_{-L_x/2}^{L_x/2}\rho X\mathrm{d}x\mathrm{d}y\mathrm{d}z}\hfill \\ \hfill & {⟨X⟩}_{x,y,z}\equiv \frac{1}{L_x{L}_y{L}_z}{\displaystyle {\int }_{-L_z/2}^{L_z/2}{\int }_{-L_y/2}^{L_y/2}{\int }_{-L_x/2}^{L_x/2}X\mathrm{d}x\mathrm{d}y\mathrm{d}z}\hfill \\ \hfill & ⟨X⟩\equiv \frac{1}{L_x{L}_y}{\displaystyle {\int }_{-L_y/2}^{L_y/2}{\int }_{-L_x/2}^{L_x/2}X\mathrm{d}x\mathrm{d}y,}\hfill \end{array}$$

where

$$\begin{array}{c}\hfill \mathrm{\Sigma }={\displaystyle {\int }_{-L_z/2}^{L_z/2}{⟨\rho ⟩}_{x,y}\mathrm{d}z.}\end{array}$$

We use square brackets to denote spatial averages and curly brackets to indicate an average over a time t_avg performed when the simulation reached a quasi-steady state. The averaging timescale t_avg is indicated for each run in the last column of Table 1. We typically averaged over 100 orbits. {Σ}, {T_mid}, and {T_eff} are thus the time-averaged values of the surface density, midplane temperature, and effective temperature, respectively. σ_Tmid, σ_Teff, and σ_α are the standard deviations of the fluctuations with time of these quantities.

We computed T_eff as

$$\begin{array}{c}\hfill T_{\mathrm{eff}}^4=\frac{1}{2{\sigma }_B}(F_{\mathrm{rad}\mathrm{z}}^{+}-F_{\mathrm{rad}\mathrm{z}}^{-}),\end{array}$$

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

$$\begin{array}{c}\hfill {⟨Q^{-}⟩}_{x,y}\equiv \frac{\mathrm{d}}{\mathrm{d}z}{⟨F_{\mathrm{rad}\mathrm{z}}⟩}_{x,y}+\frac{\mathrm{d}}{\mathrm{d}z}{⟨ϵv_z⟩}_{x,y}.\end{array}$$(5)

v_z is the vertical velocity, and the quantity ϵv_z represents the advective flux of internal energy F_adv.

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

$$\begin{array}{c}\hfill \alpha \equiv \frac{\{{⟨W_{\mathit{xy}}⟩}_{x,y,z}\}}{\{{⟨P⟩}_{x,y,z}\}},\end{array}$$

and β, the ratio of the time-averaged midplane pressure to the magnetic pressure due to the mean vertical field as follows:

$$\begin{array}{c}\hfill \beta =\frac{8\pi \{⟨P_{\mathrm{thermal}\mathrm{mid}}⟩\}}{{⟨B_{\mathrm{z}}⟩}^2}·\end{array}$$(6)

The instantaneous stress-to-pressure ratio is $\stackrel{\sim }{\alpha }={⟨W_{\mathit{xy}}⟩}_{x,y,z}/{⟨P⟩}_{x,y,z}$. The turbulent stress W_xy is ρ(u_xu_y − v_Axv_Ay), where $\mathit{u}=\mathit{v}+\frac{3}{2}\mathrm{\Omega }x\widehat{\mathit{y}}$, and v_A is the Alfvén velocity computed from the mean field.

3. Ideal MHD runs

First, we gather all of our simulations and discuss the influence of a net vertical magnetic field on the local properties of MRI in Sect. 3.1 and the properties of MRI-driven outflows in Sect. 3.2. In Sect. 3.3, we discuss the impact of a large-scale poloidal magnetic field with three different magnitudes on the S curves. In Sect. 3.4, we investigate the relative importance of viscous and wind transport of angular momentum on the ensuing accretion rate. Finally, in Sect. 3.5, we discuss the impact of the wind-driven angular momentum transport on observables.

3.1. Local turbulent transport

Local MRI turbulent transport is well known to depend on the magnetization parameter β (Hawley et al. 1995). We plot in Fig. 1 α as a function of β for all of our simulations. Squares are used for cold branch simulations and circles for the hot branch (see Sect. 3.3 for details on the cold and hot branch). The color of the circles gives f_conv, which is defined as

Fig. 1. Square markers denote cold branch simulations, and circle markers denote hot branch simulations. The color indicates the value of fconv. Horizontal dotted gray lines show the value of α in ZNF simulations from S18 in a highly convective simulation of the hot branch (upper line) and a typical cold branch case (bottom line). The linear fit is made using data from the cold branch simulations with β ≲ 2 × 104.

$$\begin{array}{cc}\hfill f_{\mathrm{conv}}& \equiv \frac{1}{2}(\max (0,\max \left(\left\{\frac{F_{\mathrm{adv}}(z>0)}{F_{\mathrm{rad}\mathrm{z}}^{+}}\right\}\right))\hfill \\ \hfill & +\max (0,\max \left(\left\{\frac{F_{\mathrm{adv}}(z<0)}{F_{\mathrm{rad}\mathrm{z}}^{-}}\right\}\right))),\hfill \end{array}$$

where F_adv = ⟨ev_z⟩. f_conv is a dimensionless measure of the importance of convective transport compared to radiative transport.

The gray dashed lines in Fig. 1 correspond to the value of α for the convective run 452O, which shows an enhanced α (top), and for the purely radiative run 434O (bottom) in S18. α ranges between 0.275 and 0.026. We find that α scales as β^−0.56 (with β defined in Eq. (6)) when we only include the simulations from the cold branch. This is consistent with previous works (Hawley et al. 1995; Salvesen et al. 2016).

As in ZNF, convection plays a role in determining the value of α (see Hirose et al. 2014 and S18). For β >  10⁵, simulations from both the hot and cold branch depart from the β^−0.56 relation and tend toward the ZNF case. Moreover, Fig. 1 shows for β  ≈  10⁴ that convective simulations from the hot branch leave the relation, whereas cold branch simulations do not. We note an increasing deviation from the trend with f_conv for these convective simulations. However, in simulations with a high level of magnetization, the disk is thicker because of the magnetic pressure in the midplane and thus pushes the zone where convection could set in out of the box. We doubled the vertical size and resolution of the disk in 456F_bigbox_β3, but only observed one convection episode. The effect on α is not clear from this single event. All we can say is that the vertical equilibrium is dominated by magnetic pressure, hence it would seem reasonable that convection plays a lesser role than that for higher β. We add the caveat that, in this particular simulation, pressure is normalized. This may have an impact on the onset of convection.

S18 showed that α is not correlated with the convective fraction. This result was also confirmed by Coleman et al. (2018), who observed that α appears to be more correlated with the vertical Mach number. However, they added that the speed of the convective eddies is only 10% of the turbulent speed. This indicates that convection is weaker than turbulence and does not significantly drive additional motion that could feed the dynamo. All these characteristics are also true for the net flux. In the absence of any deep insight on how convective MRI leads to an enhancement of α, we plot in Fig. 2, a 2D map of α in the T_eff, β space. We fit our data with Eq. (7), which is similar to Coleman et al. (2016) with an additional term for the dependence in β based on Fig. 1,

$$\begin{array}{c}\hfill \alpha =\{\begin{array}{cc}Aexp(-\frac{{(T_{\mathrm{eff}}-T_0)}^2}{2{\sigma }^2})+Btanh(\frac{T_{\mathrm{eff}}-T_0}{\sigma })+C+D{\beta }^{-E},\hfill & \text{if}\beta <5\times {10}^4\hfill \\ Aexp(-\frac{{(T_{\mathrm{eff}}-T_0)}^2}{2{\sigma }^2})+Btanh(\frac{T_{\mathrm{eff}}-T_0}{\sigma })+C,\hfill & \text{if}\beta >5\times {10}^4\hfill \end{array}.\end{array}$$(7)

We chose to constrain A, T₀, σ, B, and C using only simulations with β >  5 × 10⁴. Then, we fit the remaining parameters in the simulations with β <  10⁴. We plot the fit on the background of Fig. 2. The best-fit parameters are A = 5.35 × 10⁻², T₀ = 6866 K, σ = 853 K, B = 1.65 × 10⁻³, C = 3.65 × 10⁻², D = 1.37, and E = 0.58. This prescription will be useful in a DIM model or any global model that does not resolve the turbulent transport.

Fig. 2. Value of log(α) in our simulations is shown by color. Simulations at β = 109 are ZNF simulations taken from S18 with additional simulations from this article. The background color map shows the result of the fit of α as a function of Teff and β from Eq. (7).

3.2. Outflows and angular momentum transport

Net magnetic flux significantly affects the dynamic of the disk as it allows the formation of magneto-centrifugally driven outflows from MRI (Fromang et al. 2013; Lesur et al. 2013). All of our simulations with a β ≲ 10⁴ clearly exhibit these magneto-centrifugally driven outflows. By computing the Bernoulli invariant along streamlines, we observe that thermal effects also play a significant role in launching the outflow. Our simulations reproduce the main characteristic that have been extensively described in the literature (Blandford & Payne 1982; Ferreira & Pelletier 1995; Fromang et al. 2013; Lesur et al. 2013).

Such outflows extract matter, angular momentum, and energy away from the disk. Since the classical DIM only includes transport of angular momentum through anormalous diffusive transport, we need to quantify the transport that is due to the outflows to include them in a more sophisticated DIM. Following Balbus & Papaloizou (1999), we derived the conservation of angular momentum of an axisymmetric disk in cylindrical coordinates including the surface terms. After some algebra, we can write the surface density equation of evolution in cylindrical coordinates in the following form:

$$\begin{array}{cc}\hfill {\partial }_t\mathrm{\Sigma }& -\frac{2}{R_0}{\partial }_{R_0}(R_0^{1/2}{\partial }_{R_0}\left(\frac{\mathrm{\Sigma }}{{\mathrm{\Omega }}_0}{R}_0^{1/2}\alpha {⟨c_s^2⟩}_{\rho }\right)\hfill \\ \hfill & +\frac{2q}{\beta }⟨P_{\mathrm{thermal}\mathrm{mid}}⟩\frac{R_0}{{\mathrm{\Omega }}_0})+\zeta ⟨{\rho }_{\mathrm{mid}}{c}_{\mathrm{s}\mathrm{mid}}⟩=0.\hfill \end{array}$$(8)

Equation (8) shows three ways for the surface density to change: through radial accretion due to a transport of angular momentum either through turbulence or outflows, and loss of matter from the vertical boundaries. We have studied the radial turbulent transport of angular momentum defined by the α parameter in Sect. 3.1. Outflows lead to the definition of two additional dimensionless parameters: q, which quantifies the strength of the wind torque, and ζ, which is the ratio of the mass loss from the vertical boundaries to the sound speed times the density in the midplane. They are defined as

$$\begin{array}{c}\hfill q=\frac{{[⟨\rho u_{\varphi }{u}_z⟩-⟨B_{\varphi }{B}_z⟩/4\pi ]}_{\mathrm{surface}-}^{\mathrm{surface}+}}{2⟨P_{\mathrm{thermal}\mathrm{mid}}⟩}\beta \end{array}$$(9)

and

$$\begin{array}{c}\hfill \zeta =\frac{{\left[⟨\rho u_z⟩\right]}_{-L_z/2}^{L_z/2}}{2⟨{\rho }_{\mathrm{mid}}{c}_{\mathrm{s}\mathrm{mid}}⟩}·\end{array}$$(10)

Shearing box simulations are very useful for studying the local properties of MRI turbulence. However, they suffer limitations concerning the global geometry of the outflows. The shearing-box symmetries do not ensure a physical configuration of the magnetic field. For example, they can allow a surface torque of the same sign on both sides of the disk; this configuration does not extract angular momentum. Additionally, it is a well-known result that the sign of ⟨B_y⟩ flips regularly in a shearing box (Brandenburg et al. 1995). To cope with these limitations, we defined a function S that gives the sign of ⟨B_y⟩ in each hemisphere. We multiplied the wind torque by S in our time average. This allowed us to compute wind torques that do not depend on the geometry picked up by the shearing box at a given time.

A surface for the disk needs to be defined so that q can be estimated. We took the height where the maximum value of the wind torque is reached as the definition of the disk surface because above this surface, magnetic energy is transferred to kinetic energy to accelerate the outflow. We add the caveat that our vertically modified outflow boundary conditions impose a purely vertical magnetic field that constrains the geometry of the magnetic field lines and therefore imposes q = 0 at the boundaries. Numerically, we defined q as

$$\begin{array}{c}\hfill q_{\text{surface}}=\frac{\text{max}(S(⟨\rho v_{\varphi }{v}_z⟩+⟨B_{\varphi }{B}_z⟩)(z^{+})+S(⟨\rho v_{\varphi }{v}_z⟩+⟨B_{\varphi }{B}_z⟩)(z^{-}))}{2⟨P_{\text{thermal mid}}⟩}\beta .\end{array}$$(11)

We plot q_surface in Fig. 3 as a function of β for all simulations. We find that q_surface is roughly constant for β ≲ 10⁴ with an approximate value of q_surface  ≈  40. For low β, the wind torque is dominated by the magnetic part. This means that, remarkably, across almost three orders of magnitude in β, the MRI gives a quasi-constant ratio B_ϕ/B_z. This value agrees with local simulations from Fromang et al. (2013) and Simon et al. (2013) as well as with the global simulation of Zhu & Stone (2018), where q  ≈  10, according to their Figs. 6 and 7. Bai & Stone (2013a) found a lower value of q, but they measured the wind torque at the external boundaries; our results again agree when we adopt the same definition. When β exceeds 10⁵, we tend toward the ZNF case where the magnetic field becomes mainly azimuthal, and we expect q ∝ β^0.5. We find that a dependence in β^0.6 fits our data well. We propose in Fig. 3 a function that allows a smooth transition between the two asymptotic regimes.

Fig. 3. Squares denote cold branch simulations, and circles denote hot branch simulations. The color indicates the value of fconv.

The amplitude of the wind torque normalized by thermal pressure is given by 2q_surface/β. It varies as β⁻¹ for β ≲ 10⁴ and as β^−0.4 for lower magnetization; angular momentum is more efficiently removed by the wind at low β, as we expected (see Sect. 3.4).

Convective simulations from the hot branch have a slightly higher q_surface than cold branch simulations for β >  10⁵. It is unclear whether convection affects the ϕz component of the stress tensor for higher magnetization. We note that the magnetic component ⟨B_ϕB_z⟩ accounts for ≈70% of the wind torque for β ≲ 10⁵. For β ≳ 10⁵, the Reynolds component can account for up to ≈60% of the total wind torque.

Since the vertical extent of our box is only a few scale heights, part of the matter leaving the box, which in our simulations is definitely lost, is susceptible to falling back onto the disk. This picture is consistent with decreasing mass outflow with increasing box size, as observed in Fromang et al. (2013); this is also observed here. The value of ζ can therefore only be an estimate in a shearing-box simulation (Suzuki & Inutsuka 2009; Fromang et al. 2013). We computed the critical points for the simulation 434F_8G and found that the fast Alfvén point is reached at the vertical boundary exactly as in Fromang et al. (2013), suggesting a numerical artifact of the shearing box. Taking this into account, we conclude that the ζ we measured are maximum values for the mass loss rate through the vertical boundaries.

We show in Fig. 4 the dependence of ζ on β. Crosses show simulations on a cold branch with an extent of 12H in the z-direction, and squares show simulations with an extent of 6H in the z-direction, as for Figs. 1 and 3. The dependence of ζ on the vertical size of the box is clear. Smaller boxes exhibit higher ζ. However, ζ scales approximately as β⁻¹ for the two box sizes in the range β <  10⁵. Like α and q_surface, ζ reaches the asymptotic value given by the zero net-flux case as we increase β. The results are more dispersed than for α since by definition, they depend on the boundary conditions and the relative box size compared to the disk. All these results are consistent with Suzuki & Inutsuka (2009), Bai & Stone (2013a), and Fromang et al. (2013). Moreover, it seems that convection enhances ζ for β >  10⁵, whereas it does not for lower β. This is illustrated by the higher ζ of the convective hot branch simulations.

Fig. 4. Squares denote cold branch simulations with a vertical extent Lz = 6H, and circles denote hot branch simulations with a vertical extent Lz = 12H. Crosses are used for cold branch simulations with the same box size as on the hot branch (results from those simulations are not reported in Table 1). The color indicates the value of fconv. Horizontal dotted gray lines show the value of ζ in ZNF simulations from S18 in a highly convective simulation of the hot branch (upper line) and a typical cold branch case (bottom line). The blue and green solid lines are linear fits to the points from the cold branch simulations with β ≲ 2 × 104. This figure only reports results from simulations with Bz0 = 0.8, 2, and 8 G.

3.3. Thermal equilibrium curves with a constant value of B_z0

We plot in the left (right) panels of Fig. 5 the thermal equilibrium curves in the T_mid(T_eff) versus Σ plane, also known as S-curves, for several net magnetic fields. We emphasize that for a given S-curve, B_z0 is fixed, but not β. The S-curve in ZNF is taken from S18 with additional simulations with outflow z-boundary conditions on the cold branch to compare this with the net flux case. These additional simulations show better agreement than S18 with Hirose et al. (2014) concerning the cold branch.

Fig. 5. Thermal equilibrium curves in the Σ − Tmid(−Teff) plane for different magnetic field configurations. From top to bottom, first panel: result of S18 (circles) with outflow conditions with new additional simulations using a zero net magnetic field configuration, second panel: Bz0 = 0.8 G, third panel: Bz0 = 2 G, last panel: Bz0 = 8 G. The colors of the dot correspond to the value of α. Triangles show the instability points where the patch of disk undergoes critical cooling (downward-pointing triangle) or heating (upward). Stars indicate the simulations from the middle branch that did not converge. Gray vertical lines are plotted at the end of the hot and cold branch to facilitate comparisons between the panels. Finally, the colored S-curves in solid lines shows the DIM model with convection using different values for α.

Returning to the net flux case, we see that similarly to ZNF simulations, there are always two branches of stability: the hot branch, and the cold branch. However, as observed in Hirose (2015) and S18, there is also a middle branch, which is a continuation at higher temperature of the cold branch; these points are indicated as stars in Fig. 5. There is no obvious middle branch in the case of B_z0 = 8 G since the temperature on the cold branch gradually increases from 3005 K to 13751 K. For B_z0 = 0.8, 2 G and ZNF we note that the middle branch does not converge with the box size. The equilibrium temperature strongly depends on the box size. We do not know the origin of this dependence. Given the uncertainties regarding the stability of these points, we did not include these simulations in any other figure (see also Hirose 2015).

The extrema of the S-curves, defined by ${\Sigma }_{\text{crit}}^{+}/{\Sigma }_{\text{crit}}^{-}$ and $T_{\text{effcrit}}^{+}/T_{\text{effcrit}}^{-}$, are of special importance as they determine the shape of the light curves and in particular the presence of reflares (Lasota 2001; Coleman et al. 2016). For B_z0 <  8 G, the S-curves are not very different from the ZNF case. We find ${\Sigma }_{\text{crit}}^{-}\approx 110\text{g}{\text{cm}}^{\text{-2}}$, which is compatible with S18 given the uncertainty on the stability limit. However, for B_z0 = 8 G, we have ${\Sigma }_{\text{crit}}^{-}=80\text{g}{\text{cm}}^{\text{-2}}$, which extends the hot branch to lower densities. The same result applies to the beginning of the middle branch. However, the conclusions concerning Σ_crit⁺ are more difficult to draw from our simulations as we see a box-dependent behavior of the middle branch, as reported above. In the DIM, S-curves with higher α are shifted toward lower densities. Since increasing the magnetization gave higher values of α, our results are consistent with the DIM.

The S-curve with B_z0 = 8 G has another notable feature compared to the lower net magnetic field S-curves and S18: here alone, the maximum value of α is obtained on the cold branch. We postpone the discussion of this point to Sect. 5. For Σ <  116 g cm⁻², magnetization starts to be important on the cold branch with B_z0 = 2 G and α is greater than the ZNF case. Otherwise, in the remaining S-curve and for weaker magnetic field strengths, the values of α are similar to the ZNF case.

Finally, we note that on the hot branch our simulations are well fit by the S-curves from the DIM. However, simulations with the strongest magnetization, which are located on the cold branches of cases B_z0 = 2 and 8 G, clearly do not fit the equilibrium curves from the DIM. As the disk instability model is a purely hydrodynamic code, this is not surprising.

3.4. Vertical and radial transport of angular momentum

In a disk, the removal of angular momentum from turbulent MRI and surface torque leads to an inward flow of matter along the radial direction. We define the mass accretion rate in the disk $\dot{M}$ as

$$\begin{array}{c}\hfill \dot{M}\equiv -2\pi R_0\mathrm{\Sigma }{⟨u_{\mathrm{R}}⟩}_{\rho }.\end{array}$$

From the conservation of angular momentum, we find that

$$\begin{array}{c}\hfill \dot{M}=\frac{4\pi }{{\mathrm{\Omega }}_0}(\frac{1}{R_0}{\partial }_{R_0}(\mathrm{\Sigma }{R}_0^2\alpha {⟨c_s^2⟩}_{\rho })+\frac{2q}{\beta }⟨P_{\mathrm{thermal}\mathrm{mid}}⟩R_0)\end{array}$$

Following Fromang et al. (2013), we decompose $\dot{M}$ into a radial component, ${\dot{M}}_{R\varphi }$, and a vertical component, ${\dot{M}}_{z\varphi }$, of angular momentum transport,

$$\begin{array}{c}\hfill {\dot{M}}_{R\varphi }=\frac{4\pi }{{\mathrm{\Omega }}_0}\frac{1}{R_0}{\partial }_{R_0}(\mathrm{\Sigma }{R}_0^2\alpha {⟨c_s^2⟩}_{\rho }),\\ \hfill {\dot{M}}_{z\varphi }=\frac{8\pi }{\beta }q⟨P_{\mathrm{thermal}\mathrm{mid}}⟩\frac{R_0}{{\mathrm{\Omega }}_0}·\end{array}$$(12)

Because the curvature of the shearing box is absent, the term ${\dot{M}}_{R\varphi }$ measured in our simulation is zero. However, by making a few assumptions on the radial structure of the disk, we can have an estimate of ${\dot{M}}_{R\varphi }$. We assumed the following scalings: Σ  ≈  R^p and ⟨c²_s⟩_ρ  ≈  R^s. Equation (12) then becomes

$$\begin{array}{c}\hfill {\dot{M}}_{R\varphi }=\frac{p+s+2}{{\mathrm{\Omega }}_0}4\pi \mathrm{\Sigma }\alpha {⟨c_s^2⟩}_{\rho }.\end{array}$$

Classical viscous disk theory predicts Σ ∝ R^−3/4 and T_eff ∝ R^−3/4 (Frank et al. 2002). This means that (p + s + 2) is not far from unity. The ratio of the radial to vertical accretion rate is then

$$\begin{array}{c}\hfill \frac{{\dot{M}}_{R\varphi }}{{\dot{M}}_{z\varphi }}=(p+s+2)\frac{\beta }{2q}\frac{\mathrm{\Sigma }{⟨c_s^2⟩}_{\rho }}{⟨P_{\mathrm{thermal}\mathrm{mid}}⟩}\frac{\alpha }{R_0}·\end{array}$$

From our simulations, we find that Σ⟨c_²s⟩_ρ/⟨P_{thermal mid}⟩  ≈  C(β, H)×H, where the factor C(β, H) extends from 7.4 to 1.6 in our simulations. This means that

$$\begin{array}{c}\hfill \frac{{\dot{M}}_{R\varphi }}{{\dot{M}}_{z\varphi }}\approx C\alpha \times \frac{\beta }{2q}\left(\frac{H}{R_0}\right)\end{array}$$

We can use the empirical relations from Figs. 1 and 3, and we find

$$\begin{array}{c}\hfill \frac{{\dot{M}}_{R\varphi }}{{\dot{M}}_{z\varphi }}\approx \{\begin{array}{cc}0.19C\times {\beta }^{0.44}\left(\frac{H}{R_0}\right),\hfill & \text{if}\beta <{10}^4\hfill \\ 0.30C\times {\beta }^{0.4}\left(\frac{H}{R_0}\right),\hfill & \text{if}\beta >{10}^4\hfill \end{array}.\end{array}$$(13)

The values of ${\dot{M}}_{R\varphi }$ and ${\dot{M}}_{z\varphi }$ can be found in Table 1. We plot in Fig. 6 the ratio of the radial to vertical accretion rate as a function of β and H/R. Dashed lines show the limit ${\dot{M}}_{R\varphi }/{\dot{M}}_{z\varphi }=1$ from the empirical formulas given by Eq. (13). Simulations from the hot and cold branch are represented by circles and squares, respectively. Figure 5 shows that the cold branch is entirely dominated by accretion that is due to vertical transport of angular momentum. Only simulations from the denser and hotter part of the hot branch are dominated by viscous radial transport, and even in this case, vertical transport remains non-negligible. This shows that the vision of a purely viscous α-disk is inappropriate to describe an DNe disk if a large-scale poloidal field is present. In such a case, the role of the wind torque in removing angular momentum must be taken into account when computing the disk evolution. This can be done using the fits we provide for q_surface (see Fig. 3). Additionally, the limit ${\dot{M}}_{R\varphi }/{\dot{M}}_{z\varphi }=1$ is well reproduced by empirical formulas from Eq. (13), providing a good estimate of the constant C. As a comparison, in the global simulations from Zhu & Stone (2018), where H/R = 0.1 and β  ≈  10³, the transport is dominated by the viscous torque with a contribution of only 5% from the wind torque. Using our empirical formula with $C=\sqrt{2\pi }$, we find that the wind torque accounts for ≈20% of the total torque. Again, our simple shearing-box estimates of the wind torque are on the same order of magnitude as global simulations.

Fig. 6. Ratio of the accretion rate due to radial transport of angular momentum ${\dot{M}}_{R\varphi }$ to the accretion rate due to vertical transport of angular momentum ${\dot{M}}_{z\varphi }$ as a function of β and H/R. The color shows the value of $log({\dot{M}}_{R\varphi }/{\dot{M}}_{z\varphi })$. Squares indicates cold branch simulations and circles hot branch simulations. Dashed lines show the limit ${\dot{M}}_{R\varphi }/{\dot{M}}_{z\varphi }=1$ from the empirical formulas given by Eq. (13).

3.5. Stability criterion in a disk-wind model

The essence of the DIM is the thermal instability around ionization of hydrogen. In the DIM, the effective temperature is directly related to the mass accretion rate in the disk by

$$\begin{array}{c}\hfill \sigma T_{\mathrm{eff}}^4=\frac{3GM\dot{M}}{4\pi R_0^3}(1-{\left(\frac{R_{\mathrm{in}}}{R_0}\right)}^{1/2}).\end{array}$$(14)

In a stationary system, the mass accretion rate is imposed by the mass transfer rate from the companion. Hence, in a range of $\dot{M}$, part of the disk is ionizing hydrogen and thus unstable; this corresponds to DNe. If the accretion rate is high enough for the annulus at the external radius to be fully ionized, then the disk is stable; this a novae-like system. This simple argument provides an observational test of the DIM. Dubus et al. (2018) systematically examined the relation between $\dot{M}$ and the stability of the system for ≈100 CVs and found that predictions based on the DIM are robust.

We find that T_{eff crit⁻} does not change much between ZNF and net flux simulations. It is determined by the ionization temperature of hydrogen and is ≈7000 K. However, when we include a contribution from the wind in the conservation of energy and angular momentum, we find that the classical relation between $\dot{M}$ and T_eff in the disk changes to

$$\begin{array}{cc}\hfill \sigma T_{\mathrm{eff}}^4& =\frac{3GM\dot{M}}{4\pi R_0^3}(1-{\left(\frac{R_{\mathrm{in}}}{R_0}\right)}^{1/2}-\frac{1}{2\dot{M}{R}_0^{1/2}}{\int }_{R_{\mathrm{in}}}^{R_0}{\dot{M}}_{z\varphi }{r}^{-1/2}\mathrm{d}r)\hfill \end{array}$$(15)

$$\begin{array}{cc}& =\frac{3GM}{8\pi R_0^{7/2}}{\displaystyle {\int }_{R_{\mathrm{in}}}^{R_0}{\dot{M}}_{R\varphi }{r}^{-1/2}\mathrm{d}r.}\hfill \end{array}$$(16)

We can rewrite this equation as

$$\begin{array}{c}\hfill \sigma T_{\mathrm{eff}}^4=\frac{3GM{\dot{\mathcal{M}}}_{\mathrm{vis}}}{4\pi R_0^3}(1-{\left(\frac{R_{\mathrm{in}}}{R_0}\right)}^{1/2}),\end{array}$$(17)

where ${\dot{\mathcal{M}}}_{\mathrm{vis}}$ corresponds to the viscous mass transfer rate, which a purely viscous disk would have to match the dissipation rate of a turbulent, wind-driven disk.

In the presence or absence of wind-driven accretion, heating is only due to viscous accretion. This implies that W_Rϕ must be positive since it is of the same sign as the heating, which, by definition, is positive (see Balbus & Hawley 1998). With the assumption that W_Rϕ(R_in)=0 in a purely viscous disk, ${\dot{M}}_{R\varphi }$ needs to be constant and positive to ensure stationarity. This is no longer true when angular momentum is partly extracted by a wind. We see from Eq. (16) that in a stationary disk, ${\dot{M}}_{R\varphi }$ does not have to be positive at all radii or constant in order to ensure positive heating. This breaks down the classical relation between the viscous mass accretion rate and effective temperature (Pringle 1981) and opens up the possibility for more complex disk structures in which the disk can be viscously excreting and still be in steady state.

Equation (15) implies that a system with a high $\dot{M}$ could be unstable if there is a significant contribution to the angular momentum transport from the wind torque. The term in parentheses on the right-hand side of Eq. (15) is positive and smaller than unity. Hence, the disk heating is weaker than would be deduced from a viscous model with the same total $\dot{M}$. The disk temperature can become very low if most of the transport is due to the wind (with no associated heating), so that even a disk with a high total $\dot{M}$ could reside on the cold branch. Dubus et al. (2018) computed the mass accretion rates from optical luminosities, hence they measured ${\dot{\mathcal{M}}}_{\mathrm{vis}}$, that is, a measure of the viscous dissipation rate, but not of the total $\dot{M}$ or of ${\dot{M}}_{R\varphi }$. The separation they found between stable and unstable systems should still hold since the disk instability criterion is fundamentally based on temperature, hence on ${\dot{\mathcal{M}}}_{\mathrm{vis}}$. However, there could be a large discrepancy between ${\dot{\mathcal{M}}}_{\mathrm{vis}}$ and the total $\dot{M}$, which could be tested if a direct estimate of the mass transfer from the companion $\dot{M}$ were to be available.

4. Resistive MHD runs

Ideal MHD is a very poor approximation on the cold branch. Previous results from S18 showed that in a zero net flux configuration, a threshold in temperature exists below which MRI can no longer sustain turbulence. This transition occurs when the magnetic Reynolds number R_m is lower than ≈5 × 10³ with

$$\begin{array}{c}\hfill R_{\mathrm{m}}=\frac{c_{s}h}{\eta }·\end{array}$$

In the ZNF case, the flow has to sustain the field for the MRI to operate through a dynamo feedback loop. This positive feedback is not required in the net flux case since the field is provided by the environment, such that turbulent motions can be excited at lower R_m. Therefore the presence of a net ⟨B_z⟩ reduces the threshold below which linear MRI cannot grow and offers the possibility for MRI to survive down to lower ionization fractions. This behavior was observed in Fleming et al. (2000), who described a cyclic revival of MRI down to R_m = 50 compared to the case of ZNF, where MRI turbulence disappears for R_m <  10⁴. However, Fleming et al. (2000) considered a uniform fixed vertical profile of resistivity with an isothermal equation of state, which can lead to an artificial quasi-steady state. We wish to determine whether by relaxing this constraint, the patch of disk can keep turbulence active to lower ionization fractions than in ZNF.

Here, resistivity was computed as in Blaes & Balbus (1994):

$$\begin{array}{c}\hfill \eta =230\left(\frac{n_n}{n_e}\right)T^{1/2}{\mathrm{cm}}^2{\mathrm{s}}^{-1}\end{array}$$

with the assumption that this is dominated by electron-neutral collision. The values of η were pre-computed with the equation of state (Saha equilibrium) and saved in a table.

We restarted from previous ideal MHD simulations on the cold branch and added resistivity. We imposed a floor of 50 on the magnetic Reynolds number in order to maintain a reasonably large time step.

4.1. Impact of B_z on the resistive cold branch

R_m,  turb is the critical Reynolds number below which no turbulence is seen in our simulation. It is computed as the midplane value of R_m at the moment where the Maxwell stress ultimately decays to zero.

Following Jin (1996), we find that the linear criterion for stability of resistive MRI in the case of a vertical field and no stratification is

$$\begin{array}{c}\hfill R_{\mathrm{m},\mathrm{linear}}=\frac{\beta }{\sqrt{\frac{3}{2{\pi }^2}\beta -4}}·\end{array}$$

We plot in Fig. 7R_m,  turb as a function of β together with its analytical estimate from linear theory. As expected, the threshold for turbulence decreases with stronger magnetization. Our results agree with the linear criterion plotted as an oblique solid black line on Fig. 7.

Fig. 7. Dots show the magnetic Reynolds number Rm,  turb for which turbulence ceases in our simulation as a function of βturb. The oblique solid line is drawn from Eq. (36), which gives Rm,  linear, the linear criterion for stability of resistive MRI. The oblique dashed line gives the Rm,  cooling below which turbulence is affected by resistivity and will ultimately be suppressed.

Nonetheless, we observe that turbulence is weaker than the pure ideal case even for R_m >  R_m,  turb. In 407F_B0.8G, for example, the Reynolds magnetic number is larger than R_m,  turb for ≈100 orbits after we added resistivity. The turbulence is affected by resistivity, however; α slowly decreases, causing the disk to cool until it reaches the point where R_m <  R_m,  turb and turbulence ceases. This shows that a region exists where MRI is still linearly unstable but turbulence irremediably ceases. The reason is that as T decreases, α decreases faster than Q⁻, leading to a thermally unstable situation. The timescale on which this occurs is inversely proportional to the initial R_m from which we started resistivity. For R_m ≈ 1.5 × 10⁴ and 5 × 10³, we have a time decay for the Maxwell stress of ≈100 and 20 orbits, respectively.

Based on these observations, we define another parameter, R_{m,  cooling}, below which turbulence is still active but not self-sustaining. R_{m,  cooling} does not change much between ZNF where R_{m,  cooling} = 18000 and net flux simulations with B_z0 = 8 G where R_{m,  cooling} = 6000. We summarize these results in Fig. 7, where the dashed black line shows R_{m,  cooling} from our simulations. Three distinct zones are evident: for R_m >  R_{m,  cooling}, turbulence is similar to the ideal case; for R_{m,  cooling} >  R_m >  R_m,  turb, turbulence ceases and the disk cools until it reaches the zone R_m <  R_m,  turb, where MRI is stable.

These results show that ideal turbulence can be maintained to slightly lower ionization levels in the presence of a large-scale magnetic field. This effect is too subtle for B_z0 = 0.8 and 2 G to maintain the cold branch, however. Only for B_z0 = 8 G is the cold branch maintained down to Σ = 70 g cm⁻². Surely, the fact that R_{m,  cooling} = 6000 helps in sustaining the cold branch to such low surface densities, but the main actor is that α takes high values, up to 0.275, because of the low β. This enhanced heating ensures a level of thermal ionization that is higher than in simulations with weaker magnetic fields. In this case, two effects act together in keeping the MRI unstable: first and foremost, the enhanced heating at high magnetization, and secondarily, the effect of net magnetic flux on resistive MRI.

4.2. Transport of angular momentum from the upper ionized layers

It is probable that even with a net magnetic field, there are regions of the disk that are too cold to sustain MRI turbulence. In this last section, we therefore study the effect of X-ray irradiation on the regions where turbulence has ceased. S18 showed that X-ray ionization from the central white dwarf (WD) cannot increase the ionization fraction in the midplane above the threshold required for MRI to be unstable. This is consistent with the picture described in Gammie (1996), where an X-ray illuminated accretion disk of DNe in quiescence was organized into layers with a so-called “dead zone” in the central region. However, S18 neglected the active zones in the upper layers where the ionization fraction suffices to ensure coupling between the field and the plasma. As in protoplanetary disks, this could ensure transport of angular momentum from the wind torque (Bai & Stone 2013b; Béthune et al. 2017). The ensuing supersonic accretion of the upper layers of the disk could then lead to the required accretion rates. Following the procedure developed in Sect. 3.4, we can estimate that for α  ≈  0.035 and Σ  ≈  100 g cm⁻², the typical turbulent accretion rate is ${\dot{M}}_{R\varphi }\approx 5\times {10}^{14}\mathrm{g}{\mathrm{s}}^{-1}$. We wish to estimate the ionization fraction in the upper layers and determine whether we can obtain such an accretion rate.

Following Glassgold et al. (1997a, b), we computed the ionization rate ζ assuming a photon flux

$$\begin{array}{c}\hfill f_0=\mathcal{C}\frac{L_{\mathrm{XR}}}{4\pi R_0^{2}k{T}_{\mathrm{XR}}},\end{array}$$

where $\mathcal{C}$ takes into account the irradiation geometry and the albedo A of the disk. S18 estimated that C = 10⁻². From Byckling et al. (2010), we adopted a maximum X-ray flux L_XR of photons emitted from the white dwarf of 10³¹ erg cm⁻² with a bremsstrahlung spectrum of temperature T_XR = 10 keV. We adopted a characteristic surface density of 100 g cm⁻².

We find that ≈1% of the mass is ionized in the upper active layers. From now on, we assume that only the upper ionized layers of the disk accrete at a constant velocity v_R. The typical accretion rate on the cold branch being $\dot{\mathrm{M}}\approx 2\pi R\mathrm{\Sigma }{⟨v_{\mathrm{R}}⟩}_{\rho }\approx 5\times {10}^{14}\mathrm{g}{\mathrm{s}}^{-1}$, this leads from a surface density of 100 g cm⁻² to v_R  ≈  5 × 10³ cm s⁻¹. From the equation of angular momentum we can obtain the following relation between q and the accretion speed:

$$\begin{array}{c}\hfill v_{\mathrm{R}}=\frac{2q{B}_z^2}{3\mathrm{\Omega }\mathrm{\Sigma }{f}_{\mathrm{ioniz}}},\end{array}$$

which shows that the minimum wind torque that could explain the accretion rate on the cold branch is q_min  ≈  10. From Fig. 1, we find for β <  10⁴ that q_surface takes a constant value of ≈40, which is above q_min. Adding the caveat that q_surface is only an estimate from a shearing box in ideal MHD, we propose sonic accretion in the upper layers as a reasonable option to explain accretion in resistive regions of DNe.

5. Discussion

S18 have confirmed that MRI leads to values of α  ≈  0.1 near the tip of the hot branch when coupled non-linearly with convection. When we add a net magnetic flux, we find that convection also enhances α for magnetizations up to β  ≈  10⁴, showing that convective MRI is a reliable mechanism for enhancing classical MRI transport. The enhancement of α is localized at the tip of the hot branch and cannot explain values of 0.1 in the high-density part of the hot branch, however. This could be obtained if β were high throughout the whole hot branch. However, a strong vertical magnetic field also leads to high values of α on the cold branch. Obviously, the large magnetic field is surely not constant throughout a hysteresis cycle of DNe. Then, it is possible that the advection of the magnetic field leads to a higher magnetization during the eruptive phase than in quiescence. Unfortunately, advection of a large-scale magnetic field is a global phenomenon that cannot be modeled in the shearing-box framework (see Guilet & Ogilvie 2012, 2014 for recents developments on that topic).

We investigated the problem of angular momentum transport on the cold branch and found that resistive MRI does not naturally lead to values of α  ≈  0.01, as required by the DIM models. S18 have indeed shown that when resistivity is included in a ZNF configuration, the cold branch becomes MRI stable and the transport is quenched. The same conclusion applies for this paper in the case of B_z = 0.8 and 2 G. For B_z = 8 G, the disk remains turbulent on most of the cold branch, but exhibits values of α ≳ 0.1.

Some net magnetic flux is unavoidable in the disks of DNe since the white dwarf or companion are likely to have a non-zero dipolar field. If this leads to low values of β, then the mass accretion is mostly due to the torque from the magnetic wind that appears. This will have an effect on the evolution of the disk throughout an outburst cycle. Wind transport will also change the link between light-curve timescales and α. In a purely turbulent case, the light-curve timescales are related to the thermal (t_therm ∝ 1/αΩ) and viscous (t_vis ∼ R²/ν) timescales of the disk. For instance, the propagation time of a cooling front is ∼(t_thermt_vis)^1/2 (Menou et al. 1999; Kotko & Lasota 2012). When a wind is present, however, the observed timescales will become a diffusion or advection time related to α, but also to q. This will have to be borne in mind in future determinations of α from observations.

To build such a disk-wind model, we provided prescriptions for α and the wind torque q estimated from our shearing-box simulations. We compared our results with local simulations from Bai & Stone (2013a), Fromang et al. (2013), and Simon et al. (2013) as well as global simulations from Zhu & Stone (2018), and we found that our estimates of the wind torque and the mass accretion rates agree with those in the literature. However, an uncertainty remains concerning the dependence of q on H/R. This problem needs to be addressed with global simulations.

On a more observational side, we know that there are winds in the eruptive state of DNe that are emitted by the disk (Cordova & Mason 1982), which may be hydromagnetically launched (Cannizzo & Pudritz 1988). However, only one potential indirect detection of wind in quiescence has been reported by Perna et al. (2003). The detection is expected to be more difficult than for a disk in eruption (Drew 1990; Perna et al. 2003), and it should not be concluded from the lack of observational evidence that winds are absent in quiescence.

Mass accretion rates computed through X-ray emission allow access to the total mass accretion rates. Several sources indicate that the X-ray flux in quiescence is higher than what can be explained by the classical DIM (Collins & Wheatley 2010; Mukai 2017; Wheatley et al. 2003). One explanation is that the disk may be truncated by the magnetic field of the WD (Wheatley & West 2003; Balman & Revnivtsev 2012). However, this hypothesis is still debated, and there is also strong evidence of an accretion disk that even reaches the WD (Mukai 2017; Ishida et al. 2009). In view of our results, we suggest that accretion in quiescence is mainly driven by the wind torque. This would reconcile the high X-ray flux coming from the flow of matter $\dot{M}$ onto the white dwarf with the low values of ${\dot{M}}_{R\varphi }$ required by the DIM: in such a case the dissipation in the disk is related to ${\dot{M}}_{R\varphi }$ and not to $\dot{M}$, leading to a darker disk than a turbulence-driven disk with the same $\dot{M}$.

We also propose that a plausible solution for the colder regions of the disk is accretion that is not driven by turbulence, but by removal of angular momentum by the wind torque in the ionized upper layers, much as in protoplanetary disks. We used observational constraints of the X-ray flux in quiescence emitted from the central white dwarf to compute the ionization level in the upper layers. We estimate that it might be sufficient to explain the α that are observed. A layered accretion is likely to have an effect on how outbursts are triggered in the disk, notably on the conditions for which the heating front propagates outward (“inside-out outbursts”) or inward (“outside-in”). There is observational evidence for both in dwarf novae. In this context, the absence of heating from the wind is a major problem. This may prevent the disk from reaching high enough temperatures to ionize H and trigger the thermal-viscous instability. The trigger may then be due to a decrease in the wind torque as matter accumulates in quiescence, a profound modification to the DIM. Moreover, the lack of heating may not enable maintenance of the vertical structure, and some base level of heating is required in the poorly ionized parts of the disk.

6. Conclusion

We have studied the impact of a large net vertical magnetic field on the properties of transport and on the thermal equilibrium curves of DNe. We considered a large-scale magnetic field as could arise from the dipolar component of the primary or the companion, which is constant for a given S-curve. This configuration naturally leads to lower values of β on the cold branch than on the hot branch. For B_z <  2 G, the S-curves, in particular the critical points, are quasi-identical to the ZNF case. For B_z ≳ 2 G, the S-curves are shifted toward low density, and we found higher values of α on the cold branch than on the hot branch.

We investigated the problem of the transport of angular momentum on the cold resistive branch. We find that for B_z = 0.8 and 2 G, MRI turbulence is suppressed on most of the cold branch. However, for B_z = 8 G, the turbulence is sustained to lower surface densities, but the cold branch exhibits values of α >  0.1.

We showed that for R = 1.315 × 10¹⁰ cm, the entire cold branch is wind driven, even for local magnetic fields that are relatively weak. It is therefore essential to include wind-driven transport in the DIM and to review the effect on our understanding of outburst cycles.

Outflows provide an alternative route for extracting angular momentum that has a huge impact on the budget of angular momentum, but does not heat the disk. This means that DNe may be accreting much more than what we infer from the luminosity of the disk, and the high X-ray flux observed in quiescence may support this theory. All in all, we arrive at a point where observations and theory need to work hand in hand again if we wish to proceed in our understanding of accretion in DNe.

Acknowledgments

NS, GL, and GD thank Omer Blaes for fruitful discussions. We thank the anonymous referee for insightful comments. NS acknowledges financial support from the pole PAGE of the Université Grenoble Alpes. GL, GD, and MF 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. M. Flock has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 757957).

References

All Tables

All Figures

Fig. 1. Square markers denote cold branch simulations, and circle markers denote hot branch simulations. The color indicates the value of fconv. Horizontal dotted gray lines show the value of α in ZNF simulations from S18 in a highly convective simulation of the hot branch (upper line) and a typical cold branch case (bottom line). The linear fit is made using data from the cold branch simulations with β ≲ 2 × 104.

In the text

	Fig. 2. Value of log(α) in our simulations is shown by color. Simulations at β = 109 are ZNF simulations taken from S18 with additional simulations from this article. The background color map shows the result of the fit of α as a function of Teff and β from Eq. (7).
In the text

	Fig. 3. Squares denote cold branch simulations, and circles denote hot branch simulations. The color indicates the value of fconv.
In the text

Fig. 4. Squares denote cold branch simulations with a vertical extent Lz = 6H, and circles denote hot branch simulations with a vertical extent Lz = 12H. Crosses are used for cold branch simulations with the same box size as on the hot branch (results from those simulations are not reported in Table 1). The color indicates the value of fconv. Horizontal dotted gray lines show the value of ζ in ZNF simulations from S18 in a highly convective simulation of the hot branch (upper line) and a typical cold branch case (bottom line). The blue and green solid lines are linear fits to the points from the cold branch simulations with β ≲ 2 × 104. This figure only reports results from simulations with Bz0 = 0.8, 2, and 8 G.

In the text

Fig. 5. Thermal equilibrium curves in the Σ − Tmid(−Teff) plane for different magnetic field configurations. From top to bottom, first panel: result of S18 (circles) with outflow conditions with new additional simulations using a zero net magnetic field configuration, second panel: Bz0 = 0.8 G, third panel: Bz0 = 2 G, last panel: Bz0 = 8 G. The colors of the dot correspond to the value of α. Triangles show the instability points where the patch of disk undergoes critical cooling (downward-pointing triangle) or heating (upward). Stars indicate the simulations from the middle branch that did not converge. Gray vertical lines are plotted at the end of the hot and cold branch to facilitate comparisons between the panels. Finally, the colored S-curves in solid lines shows the DIM model with convection using different values for α.

In the text

Fig. 6. Ratio of the accretion rate due to radial transport of angular momentum ${\dot{M}}_{R\varphi }$ to the accretion rate due to vertical transport of angular momentum ${\dot{M}}_{z\varphi }$ as a function of β and H/R. The color shows the value of $log({\dot{M}}_{R\varphi }/{\dot{M}}_{z\varphi })$. Squares indicates cold branch simulations and circles hot branch simulations. Dashed lines show the limit ${\dot{M}}_{R\varphi }/{\dot{M}}_{z\varphi }=1$ from the empirical formulas given by Eq. (13).

In the text

	Fig. 7. Dots show the magnetic Reynolds number Rm,  turb for which turbulence ceases in our simulation as a function of βturb. The oblique solid line is drawn from Eq. (36), which gives Rm,  linear, the linear criterion for stability of resistive MRI. The oblique dashed line gives the Rm,  cooling below which turbulence is affected by resistivity and will ultimately be suppressed.
In the text