Turbulent and winddriven accretion in dwarf novae threaded by a largescale magnetic field
^{1}
Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), 38000
Grenoble, France
email: nicolas.scepi@gmail.com
^{2}
Max Planck Institute for Astronomy, Königstuhl 17, 69117
Heidelberg, Germany
Received:
21
July
2018
Accepted:
30
September
2018
Dwarf novae (DNe) are accreting white dwarfs that show eruptions caused by a thermalviscous instability in the accretion disk. The outburst timescales constrain α, the ratio of the viscous stress to the thermal pressure, which phenomenologically connects to the mechanism of angular momentum transport. The eruptive state has α ≈ 0.1 while the quiescent state has α ≈ 0.03. Turbulent transport that is due to the magnetorotational instability (MRI) is generally considered to be the source of angular momentum transport in DNe. The presence of a largescale poloidal field threading the disk is known to enhance MRIdriven transport. Here, we perform 3D local magnetohydrodynamic (MHD) shearingbox simulations including vertical stratification, radiative transfer, and a net constant vertical magnetic flux to investigate how transport changes between the outburst and quiescent states of DNe. We find that a net vertical constant magnetic field, as could be provided by the white dwarf or by its stellar companion, provides a higher α in quiescence than in outburst, in opposition to what is expected. Including resistivity quenches MRI turbulence in quiescence, suppressing transport, unless the magnetic field is high enough, which again leads to α ≈ 0.1. A major difference between simulations with a net poloidal flux and simulations without a net flux is that angular momentum transport in the former is shared between turbulent radial transport and winddriven vertical transport. We find that winddriven transport dominates in quiescence even for moderately low magnetic fields ∼1 G. This can have a great impact on observational signatures since winddriven transport does not heat the disk. Furthermore, wind transport cannot be reduced to an α prescription. We provide fits to the dependence of α with β, the ratio of thermal to magnetic pressure, and T_{eff}, the effective temperature of the disk, as well as a prescription for the wind torque as a function of β that is in agreement with both local and global simulations. We conclude that the evolution of the thermalviscous instability, and its consequences on the outburst cycles of CVs, needs to be thoroughly revised to take into account that most of the accretion energy may be carried away by a wind instead of being locally dissipated.
Key words: accretion / accretion disks / magnetohydrodynamics (MHD) / turbulence / convection / stars: dwarf novae
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 Xray 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 Xray 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 Rochelobe overflow from a solartype 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 α. Outburstdecay 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 magnetorotational 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 quasiomnipresent conditions in astrophysical disks. However, magnetohydrodynamical 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 MRIdriven 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 selfsustaining version of the MRI turbulence also exists in the absence of an externally imposed magnetic field, known as zero netflux 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 largescale 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 lowdensity part of the hot branch. The actual interpretation is that there exists a nonlinear 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, largescale 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 largescale 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 largescale magnetic field, and to determine whether this field could resolve the discrepancies between MRI transport and observations.
2. Methods
We adopted the local shearingbox approximation (Hawley et al. 1995) to simulate a vertically stratified patch of an accretion disk located at a distance R_{0} = 1.315 × 10^{10} cm from a 0.6 M_{⊙} white dwarf, giving an angular velocity Ω(R_{0})=5.931 × 10^{−3} s^{−1} to facilitate comparisons with Hirose et al. (2014). The simulations include radiative transport in the fluxlimited 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 shearingbox approximation, and the differential Keplerian velocity is modeled as a linear shear flow , where x, y, and z correspond to the radial, azimuthal, and vertical directions, respectively. The set of equations in the corotating frame is
The last three terms of Eq. (1) represent the Coriolis force, the tidal force, and the vertical component of the gravitational force, respectively; and are the unit vectors in the x and z directions.
Radiative transfer is treated separately from the MHD step using an implicit timestepping following the implementation of Flock et al. (2013). In this step, we solve the coupled matterradiation equations in the fluxlimited diffusion approximation
where ρ is the density, v the fluid velocity vector, P the thermal pressure, B the magnetic field vector, Φ = Ω^{2}(R_{0})( − 3x^{2} + z^{2})/2 is the gravitational potential in the corotating frame, η the Ohmic resistivity, J = (c/4π)∇ × B the current density vector, E = ϵ + 0.5ρv^{2} + 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 fluxdiffusion approximation is F_{rad} = (cλ(R)/κ_{R}(T)ρ)∇E_{R}. The flux limiter is defined as λ(R)≡(2 + R)/(6 + 3R + R^{2}) 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:
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 lowtemperature region from 2.7 < log(T)< 4.5, and OPAL (Iglesias & Rogers 1996), which covers the hightemperature 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 zerogradient extrapolation.
2.2. Boundary conditions
We used shearperiodic boundary conditions in the xdirection and periodic boundary conditions in the ydirection. We showed in S18 that the choice of vertical boundary conditions only slightly change the final thermal equilibrium with zero netflux. 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 zerogradient 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 zboundary. 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^{3}, outflows become important. In the simulation 439F with β ≈ 10^{3}, the disk would be empty after seven orbits if no normalization were used. When β ≈ 10^{2}, 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^{2} compared to only 7% when β ≈ 10^{3}. We ran two simulations with β as low as ≈10^{2} 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 Godunovtype code PLUTO (Mignone 2009). We chose a secondorder RungeKutta timeintegration 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^{−3} 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^{−6} and 5 × 10^{−2} of the initial midplane 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^{9} 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^{11} 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/Ω_{0}, thus one orbital period is equivalent to 2π in code units) and let the disk equilibrate and reach a quasisteady state (if there was one).
For simulations with β ≥ 10^{5}, the isothermal temperature T_{c0} was set to be approximately the midplane temperature found using the vertical structure code of Hameury et al. (1998) for a given surface density Σ_{0} 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^{5}, to avoid inappropriate box sizes that are due to the influence of the magnetic field, we evaluated the expected midplane 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. Σ_{0} is the initial surface density and H ≡ c_{s}(T_{c0})/Ω is the pressure scaleheight (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.
Initial parameters and results for our ideal MHD simulations.
In the following, we use several averaged quantities denoted in the following way:
where
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 quasisteady 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 timeaveraged 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
where is the radiative flux in the vertical direction on the upper or lower boundary of our simulation box. We also defined the local cooling rate as
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
and β, the ratio of the timeaveraged midplane pressure to the magnetic pressure due to the mean vertical field as follows:
The instantaneous stresstopressure ratio is . The turbulent stress W_{xy} is ρ(u_{x}u_{y} − v_{Ax}v_{Ay}), where , 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 MRIdriven outflows in Sect. 3.2. In Sect. 3.3, we discuss the impact of a largescale 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 winddriven 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 f_{conv}. 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 × 10^{4}. 

Open with DEXTER 
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^{5}, 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^{4} 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,
We chose to constrain A, T_{0}, σ, B, and C using only simulations with β > 5 × 10^{4}. Then, we fit the remaining parameters in the simulations with β < 10^{4}. We plot the fit on the background of Fig. 2. The bestfit parameters are A = 5.35 × 10^{−2}, T_{0} = 6866 K, σ = 853 K, B = 1.65 × 10^{−3}, C = 3.65 × 10^{−2}, 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 β = 10^{9} 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 T_{eff} and β from Eq. (7). 

Open with DEXTER 
3.2. Outflows and angular momentum transport
Net magnetic flux significantly affects the dynamic of the disk as it allows the formation of magnetocentrifugally driven outflows from MRI (Fromang et al. 2013; Lesur et al. 2013). All of our simulations with a β ≲ 10^{4} clearly exhibit these magnetocentrifugally 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:
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
and
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 shearingbox 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 wellknown 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
We plot q_{surface} in Fig. 3 as a function of β for all simulations. We find that q_{surface} is roughly constant for β ≲ 10^{4} 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 quasiconstant 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^{5}, 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 f_{conv}. 

Open with DEXTER 
The amplitude of the wind torque normalized by thermal pressure is given by 2q_{surface}/β. It varies as β^{−1} for β ≲ 10^{4} 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^{5}. 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^{5}. For β ≳ 10^{5}, 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 shearingbox 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 zdirection, and squares show simulations with an extent of 6H in the zdirection, 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 β^{−1} for the two box sizes in the range β < 10^{5}. Like α and q_{surface}, ζ reaches the asymptotic value given by the zero netflux 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^{5}, 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 L_{z} = 6H, and circles denote hot branch simulations with a vertical extent L_{z} = 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 f_{conv}. 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 × 10^{4}. This figure only reports results from simulations with B_{z0} = 0.8, 2, and 8 G. 

Open with DEXTER 
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 Scurves, for several net magnetic fields. We emphasize that for a given Scurve, B_{z0} is fixed, but not β. The Scurve in ZNF is taken from S18 with additional simulations with outflow zboundary 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 Σ − T_{mid}(−T_{eff}) 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: B_{z0} = 0.8 G, third panel: B_{z0} = 2 G, last panel: B_{z0} = 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 (downwardpointing 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 Scurves in solid lines shows the DIM model with convection using different values for α. 

Open with DEXTER 
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 Scurves, defined by and , 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 Scurves are not very different from the ZNF case. We find , which is compatible with S18 given the uncertainty on the stability limit. However, for B_{z0} = 8 G, we have , 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 boxdependent behavior of the middle branch, as reported above. In the DIM, Scurves with higher α are shifted toward lower densities. Since increasing the magnetization gave higher values of α, our results are consistent with the DIM.
The Scurve with B_{z0} = 8 G has another notable feature compared to the lower net magnetic field Scurves 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^{−2}, 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 Scurve 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 Scurves 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 as
From the conservation of angular momentum, we find that
Following Fromang et al. (2013), we decompose into a radial component, , and a vertical component, , of angular momentum transport,
Because the curvature of the shearing box is absent, the term 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 . We assumed the following scalings: Σ ≈ R^{p} and ⟨c^{2}_{s}⟩_{ρ} ≈ R^{s}. Equation (12) then becomes
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
From our simulations, we find that Σ⟨c_{2s}⟩_{ρ}/⟨P_{thermal mid}⟩ ≈ C(β, H)×H, where the factor C(β, H) extends from 7.4 to 1.6 in our simulations. This means that
We can use the empirical relations from Figs. 1 and 3, and we find
The values of and 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 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 nonnegligible. This shows that the vision of a purely viscous αdisk is inappropriate to describe an DNe disk if a largescale 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 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^{3}, the transport is dominated by the viscous torque with a contribution of only 5% from the wind torque. Using our empirical formula with , we find that the wind torque accounts for ≈20% of the total torque. Again, our simple shearingbox 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 to the accretion rate due to vertical transport of angular momentum as a function of β and H/R. The color shows the value of . Squares indicates cold branch simulations and circles hot branch simulations. Dashed lines show the limit from the empirical formulas given by Eq. (13). 

Open with DEXTER 
3.5. Stability criterion in a diskwind 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
In a stationary system, the mass accretion rate is imposed by the mass transfer rate from the companion. Hence, in a range of , 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 novaelike system. This simple argument provides an observational test of the DIM. Dubus et al. (2018) systematically examined the relation between 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 and T_{eff} in the disk changes to
We can rewrite this equation as
where corresponds to the viscous mass transfer rate, which a purely viscous disk would have to match the dissipation rate of a turbulent, winddriven disk.
In the presence or absence of winddriven 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, 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, 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 could be unstable if there is a significant contribution to the angular momentum transport from the wind torque. The term in parentheses on the righthand 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 . 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 could reside on the cold branch. Dubus et al. (2018) computed the mass accretion rates from optical luminosities, hence they measured , that is, a measure of the viscous dissipation rate, but not of the total or of . The separation they found between stable and unstable systems should still hold since the disk instability criterion is fundamentally based on temperature, hence on . However, there could be a large discrepancy between and the total , which could be tested if a direct estimate of the mass transfer from the companion 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^{3} with
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^{4}. 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 quasisteady 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):
with the assumption that this is dominated by electronneutral collision. The values of η were precomputed 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
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 R_{m, turb} for which turbulence ceases in our simulation as a function of β_{turb}. The oblique solid line is drawn from Eq. (36), which gives R_{m, linear}, the linear criterion for stability of resistive MRI. The oblique dashed line gives the R_{m, cooling} below which turbulence is affected by resistivity and will ultimately be suppressed. 

Open with DEXTER 
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^{4} and 5 × 10^{3}, 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 selfsustaining. 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 largescale 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^{−2}. 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 Xray irradiation on the regions where turbulence has ceased. S18 showed that Xray 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 Xray illuminated accretion disk of DNe in quiescence was organized into layers with a socalled “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^{−2}, the typical turbulent accretion rate is . 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
where takes into account the irradiation geometry and the albedo A of the disk. S18 estimated that C = 10^{−2}. From Byckling et al. (2010), we adopted a maximum Xray flux L_{XR} of photons emitted from the white dwarf of 10^{31} erg cm^{−2} with a bremsstrahlung spectrum of temperature T_{XR} = 10 keV. We adopted a characteristic surface density of 100 g cm^{−2}.
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 , this leads from a surface density of 100 g cm^{−2} to v_{R} ≈ 5 × 10^{3} cm s^{−1}. From the equation of angular momentum we can obtain the following relation between q and the accretion speed:
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^{4} 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 nonlinearly with convection. When we add a net magnetic flux, we find that convection also enhances α for magnetizations up to β ≈ 10^{4}, 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 highdensity 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 largescale magnetic field is a global phenomenon that cannot be modeled in the shearingbox 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 nonzero 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 lightcurve timescales and α. In a purely turbulent case, the lightcurve timescales are related to the thermal (t_{therm} ∝ 1/αΩ) and viscous (t_{vis} ∼ R^{2}/ν) timescales of the disk. For instance, the propagation time of a cooling front is ∼(t_{therm}t_{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 diskwind model, we provided prescriptions for α and the wind torque q estimated from our shearingbox 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 Xray emission allow access to the total mass accretion rates. Several sources indicate that the Xray 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 Xray flux coming from the flow of matter onto the white dwarf with the low values of required by the DIM: in such a case the dissipation in the disk is related to and not to , leading to a darker disk than a turbulencedriven disk with the same .
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 Xray 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 (“insideout outbursts”) or inward (“outsidein”). 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 thermalviscous 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 largescale magnetic field as could arise from the dipolar component of the primary or the companion, which is constant for a given Scurve. This configuration naturally leads to lower values of β on the cold branch than on the hot branch. For B_{z} < 2 G, the Scurves, in particular the critical points, are quasiidentical to the ZNF case. For B_{z} ≳ 2 G, the Scurves 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^{10} cm, the entire cold branch is wind driven, even for local magnetic fields that are relatively weak. It is therefore essential to include winddriven 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 Xray 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 PHY1125915). 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.ujfgrenoble.fr), which is supported by the RhôneAlpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR10EQPX2901) 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
 Bai, X.N., & Stone, J. M. 2013a, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013b, ApJ, 769, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Balay, S., Abhyankar, S., Adams, M. F., et al. 2016, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phy., 70, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Papaloizou, J. C. 1999, ApJ, 521, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Balman, Ş., & Revnivtsev, M. 2012, A&A, 546, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Byckling, K., Mukai, K., Thorstensen, J., & Osborne, J. 2010, MNRAS, 408, 2298 [NASA ADS] [CrossRef] [Google Scholar]
 Cannizzo, J. K., & Pudritz, R. E. 1988, ApJ, 327, 840 [NASA ADS] [CrossRef] [Google Scholar]
 Cannizzo, J. K., Shafter, A. W., & Wheeler, J. C. 1988, ApJ, 333, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Cannizzo, J. K., Smale, A. P., Wood, M. A., Still, M. D., & Howell, S. B. 2012, ApJ, 747, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Coleman, M. S. B., Kotko, I., Blaes, O., Lasota, J.P., & Hirose, S. 2016, MNRAS, 462, 3710 [NASA ADS] [CrossRef] [Google Scholar]
 Coleman, M. S. B., Blaes, O., Hirose, S., & Hauschildt, P. H. 2018, ApJ, 857, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Collins, D. J., & Wheatley, P. J. 2010, MNRAS, 402, 1816 [NASA ADS] [CrossRef] [Google Scholar]
 Cordova, F. A., & Mason, K. O. 1982, ApJ, 260, 716 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, S. 1988, SIAM J. Sci. Stat. Comput., 9, 445 [CrossRef] [Google Scholar]
 Drew, J. 1990, IAU Colloq. (Cambridge: Cambridge University Press), 122, 228 [NASA ADS] [Google Scholar]
 Dubus, G., OtulakowskaHypka, M., & Lasota, J.P. 2018, A&A, 617, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Ferreira, J., & Pelletier, G. 1995, A&A, 295, 807 [NASA ADS] [Google Scholar]
 Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frank, J., King, A., & Raine, D. 2002, Accretion Power in Astrophysics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F., & Menou, K. 1998, ApJ, 492, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Glassgold, A. E., Najita, J., & Igea, J. 1997a, ApJ, 480, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Glassgold, A. E., Najita, J., & Igea, J. 1997b, ApJ, 485, 920 [NASA ADS] [CrossRef] [Google Scholar]
 Grevesse, N., & Sauval, A. 1998, Space Sci. Rev., 85, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Guilet, J., & Ogilvie, G. I. 2012, MNRAS, 424, 2097 [NASA ADS] [CrossRef] [Google Scholar]
 Guilet, J., & Ogilvie, G. I. 2014, MNRAS, 441, 852 [NASA ADS] [CrossRef] [Google Scholar]
 Hameury, J.M., Menou, K., Dubus, G., Lasota, J.P., & Hure, J.M. 1998, MNRAS, 298, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Hameury, J.M., Viallet, M., & Lasota, J.P. 2009, A&A, 496, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690 [NASA ADS] [CrossRef] [Google Scholar]
 Hirose, S. 2015, MNRAS, 448, 3105 [NASA ADS] [CrossRef] [Google Scholar]
 Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., & Sano, T. 2014, ApJ, 787, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Ishida, M., Okada, S., Hayashi, T., et al. 2009, Publ. Astron. Soc. Jpn., 61, S77 [CrossRef] [Google Scholar]
 Jin, L. 1996, ApJ, 457, 798 [NASA ADS] [CrossRef] [Google Scholar]
 King, A., & Ritter, H. 1998, MNRAS, 293, L42 [NASA ADS] [CrossRef] [Google Scholar]
 King, A., Pringle, J., & Livio, M. 2007, MNRAS, 376, 1740 [NASA ADS] [CrossRef] [Google Scholar]
 Kotko, I., & Lasota, J.P. 2012, A&A, 545, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lasota, J.P. 2001, New Astron. Rev., 45, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Latter, H. N., & Papaloizou, J. C. 2012, MNRAS, 426, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Menou, K., Hameury, J.M., & Stehle, R. 1999, MNRAS, 305, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A. 2009, Mem. Soc. Astron. It. Supp., 13, 67 [Google Scholar]
 Minerbo, G. N. 1978, J. Quant. Spectr. Rad. Transf., 20, 541 [NASA ADS] [CrossRef] [Google Scholar]
 Mukai, K. 2017, PASP, 129, 062001 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C., & Terquem, C. 1999, ApJ, 521, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Perna, R., McDowell, J., Menou, K., Raymond, J., & Medvedev, M. V. 2003, ApJ, 598, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Pringle, J. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Ryan, B. R., Gammie, C. F., Fromang, S., & Kestener, P. 2017, ApJ, 840, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016, MNRAS, 457, 857 [NASA ADS] [CrossRef] [Google Scholar]
 Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 [NASA ADS] [CrossRef] [Google Scholar]
 Scepi, N., Lesur, G., Dubus, G., & Flock, M. 2018, A&A, 609, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Bai, X.N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013, ApJ, 775, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Smak, J. 1984, Acta Astron., 34, 161 [NASA ADS] [Google Scholar]
 Starling, R. L., Siemiginowska, A., Uttley, P., & Soria, R. 2004, MNRAS, 347, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Suzuki, T. K., & Inutsuka, S.I. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Tetarenko, B., Lasota, J.P., Heinke, C., Dubus, G., & Sivakoff, G. 2018, Nature, 554, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, N., & Stone, J. 2001, ApJSS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Warner, B. 2003, Cataclysmic Variable Stars (Cambridge: Cambridge University Press), 28 [Google Scholar]
 Wheatley, P. J., & West, R. G. 2003, MNRAS, 345, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Wheatley, P. J., Mauche, C. W., & Mattei, J. A. 2003, MNRAS, 345, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, L. T., & Brent, R. P. 2002, Proc. Fifth International Conference on Algorithms and Architectures for Parallel Processing (IEEE) [Google Scholar]
 Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34 [NASA ADS] [CrossRef] [Google Scholar]
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 f_{conv}. 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 × 10^{4}. 

Open with DEXTER  
In the text 
Fig. 2. Value of log(α) in our simulations is shown by color. Simulations at β = 10^{9} 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 T_{eff} and β from Eq. (7). 

Open with DEXTER  
In the text 
Fig. 3. Squares denote cold branch simulations, and circles denote hot branch simulations. The color indicates the value of f_{conv}. 

Open with DEXTER  
In the text 
Fig. 4. Squares denote cold branch simulations with a vertical extent L_{z} = 6H, and circles denote hot branch simulations with a vertical extent L_{z} = 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 f_{conv}. 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 × 10^{4}. This figure only reports results from simulations with B_{z0} = 0.8, 2, and 8 G. 

Open with DEXTER  
In the text 
Fig. 5. Thermal equilibrium curves in the Σ − T_{mid}(−T_{eff}) 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: B_{z0} = 0.8 G, third panel: B_{z0} = 2 G, last panel: B_{z0} = 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 (downwardpointing 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 Scurves in solid lines shows the DIM model with convection using different values for α. 

Open with DEXTER  
In the text 
Fig. 6. Ratio of the accretion rate due to radial transport of angular momentum to the accretion rate due to vertical transport of angular momentum as a function of β and H/R. The color shows the value of . Squares indicates cold branch simulations and circles hot branch simulations. Dashed lines show the limit from the empirical formulas given by Eq. (13). 

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

Open with DEXTER  
In the text 