Issue 
A&A
Volume 609, January 2018



Article Number  A77  
Number of page(s)  13  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201731900  
Published online  11 January 2018 
Impact of convection and resistivity on angular momentum transport in dwarf novae
^{1} Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), 38000 Grenoble, France
email: nicolas.scepi@gmail.com
^{2} Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California 91109, USA
Received: 5 September 2017
Accepted: 16 October 2017
The eruptive cycles of dwarf novae are thought to be due to a thermalviscous instability in the accretion disk surrounding the white dwarf. This model has long been known to imply enhanced angular momentum transport in the accretion disk during outburst. This is measured by the stress to pressure ratio α, with α ≈ 0.1 required in outburst compared to α ≈ 0.01 in quiescence. Such an enhancement in α has recently been observed in simulations of turbulent transport driven by the magnetorotational instability (MRI) when convection is present, without requiring a net magnetic flux. We independently recover this result by carrying out PLUTO magnetohydrodynamic (MHD) simulations of vertically stratified, radiative, shearing boxes with the thermodynamics and opacities appropriate to dwarf novae. The results are robust against the choice of vertical boundary conditions. The thermal equilibrium solutions found by the simulations trace the wellknown Scurve in the densitytemperature plane that constitutes the core of the disk thermalviscous instability model. We confirm that the high values of α ≈ 0.1 occur near the tip of the hot branch of the Scurve, where convection is active. However, we also present thermally stable simulations at lower temperatures that have standard values of α ≈ 0.03 despite the presence of vigorous convection. We find no simple relationship between α and the strength of the convection, as measured by the ratio of convective to radiative flux. The cold branch is only very weakly ionized so, in the second part of this work, we studied the impact of nonideal MHD effects on transport. Ohmic dissipation is the dominant effect in the conditions of quiescent dwarf novae. We include resistivity in the simulations and find that the MRIdriven transport is quenched (α ≈ 0) below the critical density at which the magnetic Reynolds number R_{m} ≤ 10^{4}. This is problematic because the Xray emission observed in quiescent systems requires ongoing accretion onto the white dwarf. We verify that these Xrays cannot selfsustain MRIdriven turbulence by photoionizing the disk and discuss possible solutions to the issue of accretion in quiescence.
Key words: accretion, accretion disks / convection / magnetohydrodynamics (MHD) / turbulence / stars: dwarf novae
© ESO, 2018
1. Introduction
A fundamental, yet challenging issue in accretion theory is the transport of angular momentum. Indeed, for matter to accrete it needs to transfer its angular momentum outward. Historically, the transport of angular momentum has been parametrized by the dimensionless parameter α, the ratio of the fluid stress (responsible for the transport) to the local thermal pressure (Shakura & Sunyaev 1973). Thus, this “αprescription” reduces the physics to setting a phenomenological value for α.
Dwarf novae (DNe) provide the best observational constrains on α (King et al. 2007). DNe are binary systems where matter is transferred by Roche lobe overflow from a solartype star to a white dwarf. Their light curves show periodic outbursts during which the luminosity typically rises by several magnitudes (Warner 2003). According to the disk instability model (DIM; see Lasota 2001, for a review), these outbursts are caused by a thermalviscous instability due to the steep temperature dependence of the opacity when hydrogen ionizes around 7000 K. During quiescence, the disk fills up as mass is transferred from the companion, gradually heating up until this critical temperature is reached at some radius. This triggers the propagation of heat fronts through the disk, bringing it into a hot state with a high accretion rate. The disk then empties until the temperature falls below the critical temperature and a cooling front returns the disk to quiescence. The timescales involved are set by the ability of the disk to transport angular momentum, providing an observational handle on α. Outburstdecay timescales imply α ~ 0.1 (Kotko & Lasota 2012) whereas recurrence timescales imply α ~ 0.01 in quiescence (Cannizzo et al. 1988, 2012). Although it has long been known that the DIM requires transport to be more efficient in outburst than in quiescence (Smak 1984), the physical reason for this change in α has remained elusive.
It is now widely accepted that angular momentum transport in disks is due to turbulence driven by the development of the magnetorotational instability (MRI; Balbus & Hawley 1991). The properties of MRIinduced transport have been extensively studied over the past 25 years using local, shearing box simulations. Isothermal, stratified local simulations with zero net magnetic flux show a universal value of α ~ 0.03 (Hawley et al. 1996; Simon et al. 2012), comparable to the value required for quiescent DNe. Higher values of α ~ 0.1 require an external net vertical magnetic flux (Hawley et al. 1995), whose origin in DNe is unclear. Isothermal simulations are convenient idealisations but do not account for turbulent heating or radiative losses, and thus cannot be used to investigate the thermal equilibrium states of DNe. Latter & Papaloizou (2012) carried out more realistic simulations by using an analytic approximate local cooling law in a nonstratified shearing box, and showed numerically that the disk is indeed thermally unstable in the conditions of DNe. Their zero net flux simulations give α ~ 0.01 in both cold and hot states. Hirose et al. (2014) performed the first simulations including radiative transfer, vertical stratification, and the realistic thermodynamics appropriate to DNe. They found that convection increases α to 0.1 in the hot state near the critical temperature, in the absence of net magnetic flux, providing a tantalizing solution to the change in α in DNe (Coleman et al. 2016).
Another explanation for the difference in transport efficiency between hot and cold states was proposed by Gammie & Menou (1998). In the quiescent state of DNe, the plasma is expected to be largely neutral and thus the magnetic field decouples from the disk. The conductivity drops as electrons become scarce and increasingly collide with neutrals. With the electron fraction being a strong function of the temperature, they pointed out that the MRI may not be able to grow or, at least, sustain fully developed turbulence in a quiescent DNe disk.
In light of these results, we have carried out numerical simulations (Sect. 2) to assess the impact of convection (Sect. 3) and resistivity (Sect. 4) on the transport of angular momentum in DNe in the hot, outburst and cold quiescent states (respectively). In particular, we present the first local magnetohydrodynamics (MHD) shearing box simulations with realistic thermodynamics and radiative transfer that include Ohmic diffusion and conclude on the role that Xrays from the white dwarf boundary layer may play in the ionization balance of quiescent disks.
2. Methods
We adopt the local shearingbox approximation (Hawley et al. 1995) to simulate a verticallystratified patch of accretion disk situated 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}. These values are identical to those chosen by Hirose et al. (2014) to facilitate comparisons. The simulations include radiative transport in the flux diffusion approximation and thermodynamic quantities appropriate to the temperature and density regime sampled by DNe.
2.1. Basic equations
Curvature terms are not taken into account in the shearing box approximation and the differential Keplerian velocity is modeled as a linearized shear flow , where the x, y, and z directions correspond to the radial, azimuthal, and vertical directions, respectively. The set of equations in the corotating frame is: The last three terms of Eq. (2) represent, respectively, the Coriolis force, the tidal force and the vertical component of the gravitational force; and are the units vectors in the x and z directions.
Radiative transfer is treated separately from the MHD step using an implicit timestepping following Flock et al. (2013). In this step, we solve the coupled matterradiation equations in the fluxlimited diffusion approximation where ρ is the density, the fluid velocity vector, P the thermal pressure, B the magnetic field vector, Φ = (−3x^{2} + z^{2}) × Ω^{2}(R_{0})/2 is the gravitational potential in the corotating frame, η the Ohmic resistivity, J = (c/ 4π)∇^{∇}∇ × B the current density vector, E = ρϵ + 0.5ρv^{v}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 StefanBoltzmann 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^{2}) with R ≡  ∇E  /(κ_{R}ρE) (Turner & Stone 2001). We do not take into account radiation pressure as it is negligible for the temperatures reached by DNe; Hirose et al. (2014) found that it contributes at most 6% of the gas+radiation pressure.
To close our set of equations we use the following equation of state (EOS) and internal energy function: where μ is the mean molecular weight. To compute μ and ϵ, we use precalculated tables and interpolate linearly between values of the tables. Tables are computed from the Saha equations assuming ionization equilibrium for the solar composition of Grevesse & Sauval (1998) with a hydrogen abundance X = 0.7 and a metallicity Z = 0.02. The adiabatic index Γ, the entropy and the thermal capacity C_{v} are similarly computed.
We use 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 cover the hightemperature region from 3.75 < log (T) < 8.7. We use a linear interpolation to connect the two and extend the resulting table, where necessary, using a zerogradient extrapolation. Our tables of opacities and thermodynamic quantities are fully consistent with those used by Hirose et al. (2014).
2.2. Boundary conditions
We use shearperiodic conditions in the xdirection, periodic conditions in the ydirection, and either periodic or modified outflow conditions for the zdirection to test their influence on the results (Hirose et al. 2014, used only outflow conditions). Outflow conditions usually assume a zerogradient extrapolation to the ghost cell if material is going outside of the box and prevent matter from outside from coming into the simulation box. Following Brandenburg et al. (1995), we also impose the magnetic field to point in the zdirection at the interface with the ghost cells. This circumvents a spurious increase in the magnetic field intensity that we observed, but did not further investigate, when using pure outflow conditions.
The mass in the shearing box may decrease due to the outflow boundaries. To avoid this, we normalize the total mass to the initial mass at each time step by multiplying the density by a corrective factor. For the run 442O, the mass flux is of the order of ~0.1 gcm^{2} per orbital period. The timeaveraged advective flux at the vertical boundary [⟨ ϵv_{z} ⟩ ] ( ± z_{max}) remains negligible in outflow conditions (the averaging procedure is described in Sect. 2.5).
2.3. Numerical method
We solve the MHD equations on a 3D Cartesian grid with the conservative, Godunovtype code PLUTO (Mignone 2009). We choose a second order RungeKutta time integration method. Constrained transport ensures that ∇·B = 0. The Riemann solver is HLLD except where the pressure difference between two adjacent cell exceeds five times the local pressure, in which case we allow for shock flattening by using the more diffusive solver HLL. To solve the radiative transfer equations, we follow the same implicit scheme as Flock et al. (2013) except that we use the 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 use floors of 10^{6} and 5 × 10^{2} of the initial midplane values for the density and the temperature, respectively. We find that floors are activated over ≲2% of the box and, thus, have no significant impact on the thermodynamical equilibrium of the simulations.
2.4. Initial conditions
The initial flow is Keplerian with a weak zero net flux magnetic field of the form B_{z} ∝ sin(2πx). We usually start with an isothermal layer and assume hydrostatic equilibrium to fix the initial vertical density profile. We let the MRI develop and then trigger radiative transfer after 48 local orbits (time is normalized by the quantity 1/Ω_{0} thus one orbital period is equivalent to 2π in code units). This time is sufficient for the MRI to saturate and for the disk to reach a quasisteady state. The isothermal temperature T_{c0} is set to the midplane temperature found by using the vertical structure code of Hameury et al. (1998) for given surface density Σ_{0} and effective temperature T_{eff}. This code solves the vertical structure equations assuming an αprescription for the angular momentum transport and associated heating rate, radiative transfer in the diffusion approximation, and convection described by mixing length theory with a mixing coefficient α_{ml} = 1.5 (based on models of the Sun). We then activate radiative transfer and let the disk equilibrate and reach a quasisteady state (if there is one).
Some runs are initialized from another simulation by changing the surface density manually instead of starting from an isothermal layer. This has proven useful close to the hydrogen ionization regime, where the disk easily undergoes critical heating/cooling. In these cases, starting from an isothermal state can cause the disk to miss the equilibrium state. Starting from a nearby quasi steadystate allows for a smoother transition to capture the thermal equilibrium.
Initial parameters and results for our ideal MHD simulations.
2.5. Runs and diagnostics
Table 1 lists the runs that we performed. We adopt the notation of Hirose et al. (2014) for labeling the runs, with the addition of a final O to indicate simulations with vertical outflow boundary conditions or a final P for a periodic vertical boundary. Σ_{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 is of ±6H for the hot branch and ±3H for the cold and middle branch. L_{x},L_{y} and L_{z} follow the ratio 1:2:4 (see Table 1).
The horizontal average of a quantity f is defined as: (9)where the integration in x and y is done over the whole simulation box. We use ⟨⟩ brackets to denote spatial averages and [] brackets to indicate an average over a time t_{avg} performed once the simulation has reached a quasisteady state. The averaging timescale t_{avg} is indicated for each run in the last column of Table 1. We typically average over a hundred orbits except for some cases where the variability occurs on long timescales (e.g., 401O, see Sect. 3.2.3). [Σ], [T_{mid}], [T_{eff}], and [τ_{tot}] are thus the timeaveraged values of the surface density, midplane temperature, effective temperature and total optical depth, respectively. σ_{Tmid}, σ_{Teff} and σ_{α} are the standard deviations of the fluctuations with time of these quantities.
We compute τ_{tot}, Σ, T_{eff} as and (12)where is the radiative flux in the vertical direction on the upper/lower boundary of our simulation box. We also define ⟨ Q^{−} ⟩ _{x,y} the total cooling rate as (13)v_{z} the vertical velocity and the quantity ϵv_{z} represents the advective flux of internal energy F_{adv}. The integration over z extends over the full box.
We define α, the ratio of stress to pressure, as (14)and the instantaneous stress to pressure ratio as . The turbulent stress W_{xy} is ρ(v_{x}v_{y}−v_{Ax}v_{Ay}) where v_{A} is the Alfvén velocity.
2.6. Numerical convergence
We take 32 points in the x direction, 128 points in the y direction and 256 points in the z direction in all our simulations. This number of points must be sufficient to resolve the most unstable wavelength of the MRI, which is of order H (as the root mean square turbulent magnetic field is of the order of 0.1–0.2 for all simulations), while still allowing for a significant vertical extension. Hence, we use different values of L_{x}, L_{y} and L_{z} on the cold branch and on the hot branch. We note n the number of points per final timeaveraged scale height h ≡ c_{s}( [T_{mid}] )/Ω. On the high Σ part of the hot branch n is of the order of 20 (Table 1), enough to resolve the MRI. But at lower Σ some of our simulations are restarted from previous runs with higher temperatures, decreasing n to 15 in run 452P. We checked those cases with additional runs with a more appropriate box size and resolution and found that the results are quantitatively the same. For instance, we find a difference of ~5% in T_{mid} and ~2% in α between the run 442O started with an isothermal layer (with n = 21) or initialized from run 437O (with n = 13).
High Alfvénic velocities are allowed near the boundaries when the box height is large, about 12H, increasing the numerical diffusion due to the HLLD solver and providing a numerical source of heating. On the hot branch, the contribution of this numerical coronal heating with respect to the total numerical heating is negligible: we performed the simulation 437O with L_{x}/H = 1.5, L_{y}/H = 6.0 and L_{z}/H = 8.0 and found a difference on [T_{mid}] of only 2% for the hot branch. The effect is much stronger on the middle and cold branch, heating being weaker, leading us to take smaller box sizes. Consequently, n is approximately two times larger in the cold branch than on the hot branch.
According to Ryan et al. (2017), isothermal stratified simulations are not fully resolved since α decreases as n^{− 1/3} up to the highest achievable resolutions. Therefore, a caveat when comparing runs on the cold and hot branch is that the twofold increase in resolution may lower the value of α by ~20% in the cold branch simulations. We ran the simulation 434O with half the resolution on each direction and found α = 0.019. This difference is of the same order as the estimated error on α. Hence we cannot conclude about consistency with Ryan et al. (2017). A proper convergence study would be useful to know if this result holds for our simulations including realistic thermodynamics and radiative transfer.
Finally, Simon et al. (2012) showed that the box size L_{x} and L_{y} must be ≥2H to have converged diagnostics of the MRI. Again, these results were obtained using an isothermal equation of state and we do not know for sure if they hold here. According to this criterion, our box size is appropriate on the hot branch but may be too small in x on the cold branch (Table 1). In particular, convective eddies in our simulations could be limited by small box sizes. We doubled the size in x and y of the run 442O to check for this effect and found no notable difference.
Fig. 1
Thermal equilibria in the [T_{mid}] vs. [Σ] (top) or [T_{eff}] vs. [Σ] (bottom) plane. Squares and circles are for periodic and outflow runs, respectively. Triangular dots represent runs with runaway cooling (triangle facing down) or heating (triangle facing up). Error bars represent the standard deviation of the temperature fluctuations. The symbols are colorcoded to the value of α. The colorcoded curves are vertical thermal equilibria using an αprescription. The dashed blue line indicates where the magnetic Reynolds number R_{m} = 10^{4} based on an isothermal model (see Sect. 4.1). 
3. Ideal MHD runs
3.1. Thermal equilibrium curve
Figure 1 shows the thermal equilibria reached by the simulations listed in Table 1. The top panel shows the timeaveraged midplane temperature [T_{mid}] ≡ [ ⟨ T(z = 0) ⟩ _{x,y}] versus the surface density [Σ] (top panel). The bottom panel shows [T_{eff}] versus [Σ] (bottom panel). The intensity of α is indicated for each run by the color of the dot. The colored curves in Fig. 1 are thermal equilibria curves obtained by using the Hameury et al. (1998) code, which assumes an αprescription and allows for thermal convection (Sect. 2.4). Whereas the latter requires choosing a value for α, the simulations do not have this degreeoffreedom since the value of α is set by the turbulent angular momentum transport generated by the MRI.
The simulations trace an equilibrium thermal curve in the shape of an S (an Scurve) with a hot stable branch and a cold stable branch, providing independent confirmation of the results of Hirose et al. (2014). The Scurve from the simulations is comparable to that predicted by the vertical structure calculations using an αprescription. These also predict a third stable branch at intermediate temperatures (T_{mid} ≈ 10^{4} K). Indeed, Fig. 1 shows also a set of three stable runs in an intermediate regime of temperature. A middlebranch in radiative MRI simulations was also reported by Hirose (2015) in the context of protoplanetary disks. The middle branch is not as extended as the other branches and may be seen as a prolongation of the cold branch to higher temperatures. In fact, the middle branch does not overlap with the cold branch, so there is at most two equilibrium states for a given Σ, not three. However, the middle branch runs have a higher opacity than the cold branch and are convectively unstable (Sect. 3.2).
The hot and cold/middle branches are terminated at, respectively, low and high Σ by unstable runs displaying runaway cooling or heating (noted R in Table 1 and indicated by arrows in Fig. 1). We have followed these runaways to the stable equilibrium on the opposite branch using appropriate box sizes. Periodic simulations seem to be slightly more stable around points of critical heating/cooling as 453O is unstable whereas 453P is not. Restarting simulations from previous runs by changing the surface density allows us to follow the hot branch to lower surface densities than in Hirose et al. (2014) where 442O is their last stable run.
We have good agreement between our simulations and the αprescription calculation concerning the equilibrium temperatures at given α and the critical points for runaway heating/cooling. The equilibrium values follow the calculated curves corresponding to the value of α found in the simulation. Somewhat unexpectedly, the agreement extends to the middle branch even though its location in the (Σ,T) plane depends on the chosen value of the convection mixing length parameter (Cannizzo & Wheeler 1984). Here, we took α_{ml} = 1.5 as in the solar convection zone (Sect. 2.4).
3.2. Convection and transport of angular momentum
3.2.1. Enhancement of α
We see an enhancement of α in the low Σ part of the hot branch, albeit with slightly lower maximum values of α ≈ 0.098 than the α ≈ 0.121 found by Hirose et al. (2014). On the cold and middle branch, we find values of α ranging from 0.029 to 0.042 typical of a zero net vertical flux stratified MRI simulation (Simon et al. 2012). The values of α are comparable on the hot branch for 174 g cm^{2} ≤ Σ ≤ 1030 g cm^{2}. For higher Σ the value of α increases up to 0.066; we did not further investigate this trend as this part of the Scurve implies very high accretion rates that are unlikely to be relevant to dwarf novae. The value of α also increases as we approach the tip of the hot branch, going to a maximum of 0.098 for the lowestdensity run on the hot branch (453P). This enhancement is clearly seen when plotting α against temperature (Fig. 2, see also Fig. 1 in Coleman et al. 2016). The enhancement is accompanied by stronger fluctuations as we near the unstable tip of the hot branch (see Sect. 3.2.3).
Fig. 2
α as a function of [T_{mid}] (top panel) or [T_{eff}] (bottom panel). Error bars represent the standard deviations of the fluctuations around the mean values. Blue, green, and red colors are for the cold, middle and hot branches, respectively. Squares and circles indicate periodic and outflow runs. The shaded area corresponds to the thermally unstable region. 
3.2.2. Convective transport of heat
Hirose et al. (2014) attributed the enhanced α to convection. We also observe that the enhanced α runs are convectively unstable and further note that the value of α does not depend on the chosen vertical boundary conditions. We measure the stability of a run against convection from the BruntVäisälä frequency N, defined by (15)where Γ ≡ (∂lnP/∂lnρ)_{s} is the adiabatic index and s is the specific entropy. N^{2} has the sign of the entropy gradient along z, thus a negative BruntVäisälä frequency denotes a convectively unstable vertical profile. The BruntVäisälä frequency gives a linear stability criterion for the convective stability only in regions where the total pressure is dominated by thermal pressure and does not give direct information about the convection flux.
Figure 3 shows the radiative, advective, and total fluxes (Eq. (13)) as well as the BruntVäisälä frequency as a function of height for several illustrative cases. The top panel shows a highly convective episode (averaging over 16 local orbits) of the run 442P (hot branch). The advective flux is the main source of heat transport in the regions where the disk is convectively unstable (N^{2}< 0) as opposed to the upper regions where the radiative flux dominates. Averaging the same run on a longer time (80 local orbits) including convective and nonconvective episodes shows the contribution of the advective flux decreases (middle panel). We also see evidence for downward transport of heat between the midplane and H/R_{0} ~ 2.3 × 10^{2}. This downward transport is also present in the top panel but is less extended in height. Convection is expected to drive heat upward, where the entropy is lower, hence this downward transport is not due to convective motions. In fact, run 434F on the cold branch is convectively stable yet also shows this downward advective transport (bottom panel). We attribute this transport to vertical mixing by the MRIdriven turbulence. Because motions are approximately adiabatic in this region, this mixing tends to flatten the entropy profile, resulting in a downward transport of heat. Thus, the advective flux includes contributions from both convective motions (when present) and turbulent mixing.
Besides the region of enhanced α near the unstable tip of the hot branch, we find convection also plays a major role in the middle branch runs, notably run 441O, with the fluxes behaving as in the top panel of Fig. 3.
Fig. 3
Timeaveraged vertical profiles of the radiative, advective, and total flux for different runs. Top panel: run 442P during a highly convective episode. Middle panel: same run but averaging over convective and non convective episodes. Bottom panel: run 434P where there is no convection. The vertical line indicates the height above which magnetic pressure becomes larger than thermal pressure. The dashed line shows N^{2}/ Ω^{2} where N is the BruntVäisälä frequency. 
Fig. 4
Evolution of as a function of time for runs located on different parts of the Scurve: (top panel) 434P is a cold branch run, 468P a nonconvective hot branch run, and (bottom panel) 442 a convective hot branch run. For the latter we show the simulation run with outflow (442O) and periodic boundary conditions (442P). We also plot, for run 442O, the contribution of the Maxwell stress to as a dashed line. 
3.2.3. Convective/radiative cycles
In the following, we make use of the convective fraction f_{conv} of the flux, defined as (16)We also use an instantaneous convective fraction obtained by smoothing over a time window with a width of 1.6 local orbits. We note that can have values >1 due to the fluctuations remaining in the vertical profile.
Figure 4 shows the evolution of with time for runs 434P, 468P, 442P and 442O. The first two are convectively stable and show only weak fluctuations in . However, run 442 is located near the tip of the hot branch and is convectively unstable. This run shows has large fluctuations, going through cycles with bursts of strong angular momentum transport. We see from Fig. 4 that the main contribution to is the Maxwell part. The maximum and the recurrence of these bursts in increase towards the unstable tip of the Scurve. For comparison, there is one cycle every ~20 orbits in run 437O (Σ = 174 g cm^{2}) and one every ~10 orbits in run 442O (Σ = 113 g cm^{2}).
Fig. 5
Autocorrelation function f_{auto} as a function of time for the run 442O. τ_{0} is the parabolic decay time of the autocorrelation function (with the parabolic fit shown as a dashed line). τ_{1} is the location of the first peak of f_{auto}. 
To be more quantitative about the frequency of the cycles of , we compute its autocorrelation function (17)Figure 5 shows this function for run 442O. We extract two characteristic times. The first is the parabolic decay time τ_{0} of f_{auto} (τ_{0} ≈ 2 local orbits in Fig. 5); it corresponds to the time on which small fluctuations of remain correlated. The second, τ_{1}, corresponds to the first peak of f_{auto} (τ_{1} ≈ 9 local orbits in Fig. 5); this corresponds roughly to the period of the convective cycles.
We do not see any clear evolution of τ_{0} as the convective fraction increases on the hot branch. However, τ_{1} clearly decreases with f_{conv} (Fig. 6). On the middle branch, the cycles are much longer for the same f_{conv}. For instance, run 401O has f_{conv} = 0.14 and τ_{1} ≈ 111 local orbits, requiring a longer averaging timescale t_{avg} than for the other runs.
Fig. 6
Period τ_{1} of the cycles in α as a function of the convective fraction for convectively unstable simulations on the hot branch. 
The stronger variability near the unstable tip of the hot branch is also manifest as increased temperature fluctuations (error bars in Figs. 1 and 2). These fluctuations in temperature are anticorrelated with on the hot branch (Fig. 7, top panel). However, the fluctuations are in phase at the highest temperatures on the middle branch (Fig. 7, bottom panel). The behavior on the hot and middle branch is very different although the fractional amplitudes ([max – min]/[max + min]) of the fluctuations in and T_{mid} are the same for both. This difference can be explained by noting that the temperature gradient, assuming only radiative transport, is ∝κ/T^{3}, and that convection requires this gradient to be greater than the adiabatic gradient of the gas. On the hot branch, the opacity decreases when the temperature rises, flattening the temperature profile and quenching convection. Upward fluctuations in temperature lead the disk from a convective regime to a radiative regime and vice versa. The opposite correlation is expected on the middle branch since the opacity rises steeply with T, triggering a convective regime when fluctuations raise the temperature.
Fig. 7
vs. T_{mid} for two representative convectively unstable simulations: hot branch run 442O (top) and middle branch run 401O (bottom). 
This difference in behavior between the convective parts of the hot branch and the middle branch also extends to α. Run 442O on the hot branch shows a trend for the instantaneous to increase with the instantaneous convective fraction (Fig. 8, green markers), but this trend is absent in run 401O on the middle branch (Fig. 8, blue markers). Furthermore, the maximum value of is lower than on the hot branch run even though reaches higher values. The runs also show a difference in duty cycle. The hot branch run spends roughly equal time in the low and high parts of the diagram. In contrast, the middle branch run spends most of its time in the low part of the diagram, with only temporary excursions to the higher values. This will naturally lead to a smaller timeaveraged α in the middle branch. Figure 9 shows no trend in the convective runs between the timeaveraged values of α and f_{conv}. We conclude that there is no simple relationship between the enhancement of α and the relative strength of the convection measured by f_{conv}. Finding the physical origin for this enhancement is likely to require a detailed study of the impact of vertical convection on the MRI dynamo, which is beyond the scope of this paper.
Fig. 9
α vs. convective fraction f_{conv} for all convectively unstable runs (same color and symbol coding as Fig. 2). 
3.3. Impact of the vertical boundary conditions
We tested outflow and periodic vertical boundary conditions (Sect. 2.4) because of their potential impact on the interplay between the MRI dynamo and convection. In periodic runs, the boxaveraged ⟨ B_{y} ⟩ _{x,y,z} = 0 by construction whereas it fluctuates over time with outflow boundary conditions (Fig. 10). The flips in sign are related to the dynamo of stratified MRI turbulence and are still poorly understood (Brandenburg et al. 1995). We found no quantitative difference between periodic and outflow runs so these flips have no impact on α or on the thermal equilibrium values.
Fig. 10
Average azimuthal magnetic field ⟨ B_{y} ⟩ _{x,y,z} as a function of time for run 446 (outflow O or periodic P boundary conditions). 
The similar behavior between outflow and periodic runs extends to the “butterfly” diagram that emerges from maps of the time evolution of ⟨ B_{y} ⟩ _{x,y}, the vertical profile of the azimuthal field. Nonconvective MRI simulations are known to form “butterfly wings” as the magnetic field flips sign at each quasiperiodic pulsation and propagates outward (e.g. see Fig. 12). Coleman et al. (2017) reported that the quasiperiodic pulsation does not always lead to a sign flip of ⟨ B_{y} ⟩ _{x,y} when convection is active. We confirm that this is also the case in our convective outflow simulations (Fig. 11, top panel) and in our periodic boundary simulations (bottom panel), even though the magnetic field must adopt a symmetric configuration in the periodic simulations to ensure that ⟨ B_{y} ⟩ _{x,y,z} = 0. We also see this behavior in the middle branch during convective episodes.
Fig. 11
Vertical profile of azimuthal magnetic field ⟨ B_{y} ⟩ _{x,y} (in code units) as a function of time with outflow (top) and periodic (bottom) boundary conditions for run 446 (hot branch). 
4. Resistive MHD runs
Ideal MHD may not apply to the cold branch due to the very low ionization fractions. The development of the MRI can be suppressed as the resistivity due to electronneutral collisions increase. This is quantified by the magnetic Reynolds number (18)with η being the resistivity, c_{s} the sound speed and h = c_{s}/ Ω the local pressure scaleheight. For R_{m} < 10^{4}, it is expected that diffusion of the magnetic field becomes too important for the disk to sustain MHD turbulence (Hawley et al. 1996; Fleming et al. 2000). Gammie & Menou (1998) pointed out that R_{m} is of the order of 10^{3} in a quiescent DNe disk.
Fig. 12
Time evolution of ⟨ R_{m} ⟩, ⟨ B_{y} ⟩ and α for run 435PR (Σ = 191 g cm^{2}). The vertical black line marks the time at which the simulation is restarted from 435P with resistivity turned on. 
We investigate whether or not our cold branch runs can maintain turbulence when resistivity is included. The resistivity is computed as in Blaes & Balbus (1994): (19)with the assumption that we are dominated by electronneutral collision. The values of η are precomputed with the equation of state (Saha equilibrium) and saved in a table. The time step Δt associated with the diffusive term is of the order of Δx^{2}/η where Δx is the cell minimum length. To avoid dramatically small time steps, we impose a minimum value of R_{m} = 50 in our runs. This floor is largely under the limit for MRI turbulence to be suppressed and does not affect our results. We find Δt ~ 10^{3} to be of the order of the time step of the ideal MHD run. The resistive runs are restarted from the periodic runs in ideal MHD and labeled PR to differentiate them.
4.1. Turbulence decay due to resistivity
We find that the resistivity has a critical impact by suppressing turbulence on the cold branch below some critical density located between Σ = 174 g cm^{2} (run 462P) and Σ = 191g cm^{2} (run 435). Figures 12, 13 show the temporal evolution of , ⟨ R_{m} ⟩ _{x,y}, and ⟨ B_{y} ⟩ _{x,y} for those two runs before and after resistivity is switched on.
Turning on resistivity has no influence on run 435PR (Fig. 12). R_{m} stays well above the critical value of 10^{4} in the midplane. Some regions above the photosphere have a lower R_{m}, which is without consequences as most of the heating due to MRI happens in the densest part of the disk.
Run 462P has a lower Σ and temperature ([T_{mid}] ≈ 3100 K compared to 3750 K for 435P). Here, R_{m} starts to drop as soon as resistivity is switched on (Fig. 13). The turbulence is suppressed once R_{m} decreases below ≈5000 and the transport of angular momentum ceases. ⟨ B_{y} ⟩ _{x,y} diffuses as expected: first near the photosphere and then, as R_{m} decreases, in the whole box.
4.2. Influence of ambipolar diffusion and Hall effect
Let us first compare the order of magnitude of Ohmic (O), ambipolar (A) and the Hall effect (H). From Balbus & Terquem (2001), we have We take 462PR as a fiducial run (as it is the first one to be altered by Ohmic diffusion) to compare the importance of each nonideal effect. The typical values deduced from this run for the density, temperature, ionization fraction and magnetisation v_{A}/c_{s} are 5 × 10^{7}g cm^{3}, 3000 K, 10^{5.5} and 0.1, respectively. Therefore, in this run, ambipolar diffusion is expected to have a negligible impact while the Hall effect is only an order of magnitude weaker than Ohmic diffusion. We note that these conclusions also hold as the disc cools down when MRI is suppressed.
However, comparing the amplitude of each nonideal effect might not be enough to conclude on the dynamics since they affect the MRI in different ways. Let us therefore analyze the impact of ambipolar diffusion and of the Hall effect individually.
Ambipolar diffusion is characterized by the ambipolar Elsasser number (Balbus & Terquem 2001): (22)where ρ_{i} = n_{i}m_{i} is the ion density and γ_{i} = ⟨ σv ⟩ _{ni}/ (m_{n} + m_{i}) = 2.7 × 10^{13}(41.33m_{u}/m_{n} + m_{i})cm^{3} g^{1} s^{1} is the ionneutral drag coefficient. Following Balbus & Terquem (2001), we assume a momentum rate coefficient for ionneutral collision ⟨ σv ⟩ _{ni} = 1.9 × 10^{9}cm^{3} s^{1} (Draine 2010), m_{n} = 2.33m_{u} and m_{i} = 39m_{u} (the neutral and ion masses respectively). We obtain Λ_{A} ~ 10^{3} for our typical values. Hawley & Stone (1998) and Bai & Stone (2011) found that angular momentum transport is not significantly impacted when Λ_{A}> 100. Their result is based on numerical simulations with a net magnetic flux and should still hold when there is no net flux: a net magnetic flux implies currents outside of our domain of simulation that are not affected by resistivity and help sustain MRI turbulence. We conclude that the influence that ambipolar diffusion would have on run 462PR is negligible compared to Ohmic diffusion.
The Hall effect is characterized by the Hall Lundquist number (23)The impact of the Hall effect on the MRI saturation level is twofold. For strong Hall effect (ℒ_{ℋ} ≲ 5), turbulent transport is essentially suppressed and the flow relaxes into a selforganized state (Kunz & Lesur 2013). For weaker Hall effect (ℒ_{ℋ} ≳ 10) the system produces sustained turbulent transport provided that the magnetic Reynolds number is sufficiently large, similarly to the pure Ohmic case. This critical Reynolds number for sustained turbulence is reduced to 10^{3} for ℒ_{ℋ} ≃ 30 while it is identical to the pure Ohmic case for ℒ_{ℋ} ≳ 60 (Sano & Stone 2002), indicating that the Hall effect becomes negligible above this threshold.
Our typical values imply ℒ_{H} ~ 10^{3}, so the Hall effect does not affect the MRI close to the region where Ohmic diffusion becomes problematic. As the discs cools down in the MRI stable region, ℒ_{H} decreases, making the Hall effect stronger and potentially reviving the MRI. However, the ratio ℒ_{H}/R_{m} does not depend on the temperature, meaning that at the temperature where ℒ_{H} ≲ 50, one also expects R_{m} ≲ 50, well below the critical R_{m} found in HallMRI simulation (Sano & Stone 2002). We conclude that the Hall effect does not affect the critical temperature (or Σ) below which DNe are MRIstable, nor does it revive the MRI at lower temperatures.
4.3. Influence of Xray irradiation
The ionization fraction x_{e} is critical to determine whether or not the coupling between fluid and magnetic field is strong enough for MRI turbulence to be active. In the resistive runs described above (Sect. 4.1) x_{e} ≈ 4 × 10^{6} in the MRIstable run 462PR (at midplane) and x_{e} ≈ 6 × 10^{6} in the MRIunstable run 435PR. Therefore, a very small ionization fraction of order 5 × 10^{6} can be sufficient to maintain turbulence. These ionization fractions are thermal and do not take into account external sources of ionizations. In particular, quiescent DNe show hard Xray emission originating from the accretion boundary layer onto the white dwarf that may provide sufficient ionization to maintain the MRI active in the quiescent disk.
The quiescent Xray spectra are bremsstrahlung with temperatures kT_{XR} ranging from 1 to 10 keV and luminosities L_{XR} ranging from 10^{28} to 10^{32} erg s^{1} (Byckling et al. 2010), properties akin to those of a protostar. The typical surface densities under consideration are similar to those of a protostellar disk at 1 AU but the DNe disk is ~10^{3} closer to the Xray source, so the impinging ionizing flux is much higher.
We calculate the expected ionization fraction x_{e} by following studies of protoplanetary disks. We first compute the ionization rate ζ as in Glassgold et al. (1997a,b), assuming a photon flux (24)where is a factor taking into account the irradiation geometry and the albedo A of the disk. The Xray emission region is about the size of the white dwarf (Mukai 2017), which is much greater than the disk height. In this case, we may expect (King 1997) (25)with a white dwarf radius R_{WD} = 10^{9} cm.
We obtain x_{e} by solving Eq. (11) from Fromang et al. (2002), which relates x_{e} to the fraction of metals x_{M} = n_{M}/n_{n}, the dissociative recombination rate for molecular ions, the radiative recombination rate for metal atoms and the rate of charge exchange between them. We further assume that the disk is passive, that is, has no internal heating, and the only source of heat is Xray irradiation. Hence, in steadystate, the disk is isothermal at the temperature set by . The vertical structure is in hydrostatic equilibrium with Σ = 100 g cm^{2} to compare to our results with the first MRIstable run, 462PR.
We find that Xray ionization is negligible, in agreement with the orderofmagnitude estimate of Gammie & Menou (1998). For example, with and kT_{XR} = 10 keV, we find ζ = 3.4 × 10^{9} s^{1} at the disk midplane. The isothermal atmosphere has T = 1700 K and a midplane density n = 3.7 × 10^{17} cm^{3}, giving x_{e} = 4.1 × 10^{10} when no metals are present (x_{M} = 0) and x_{e} = 1.1 × 10^{7} when x_{M} = 6.86 × 10^{6} (solar abundance). The fraction x_{e} decreases for lower values of or kT_{XR}. Hence, Xray ionization is unable to provide the required ionization fraction except in the upper layer. However, only a small amount of material is active in this layer: taking the values assumed above, x_{e} > 5 × 10^{6} for z > 1.5 × 10^{8} cm (about 2H) and the column density of the active layer represents less than 3% of Σ.
5. Conclusions
We have performed stratified, radiative, ideal and resistive MHD, shearing box simulations in conditions appropriate to DNe. We find that the thermal equilibrium solutions found by the simulations trace the wellknown Scurve derived from αprescription models, including a middle branch that extends the cold branch to higher temperatures and is characterized by vigorous convection. We confirm that α increases to ≈0.1 near the unstable tip of the hot branch as reported by Hirose et al. (2014). This increase is thus robust against the choice of numerical code and, as we investigated, against the choice of outflow or periodic vertical boundary conditions. Although convection plays a major role in transporting heat in the runs with an enhanced α, we find no clear relationship between α and f_{conv}, the average fraction of the flux carried by convection. Notably, α is not enhanced on the middle branch although it is strongly convectively unstable. Detailed studies of the interplay between convection and the MRI dynamo will be required to understand the exact mechanism behind the enhancement of α towards the tip of the hot branch.
Ohmic dissipation can prevent MRIdriven turbulent transport in the cold, quiescent state of DNe (Gammie & Menou 1998). We verified that Ohmic dissipation is the dominant nonideal effect and carried out resistive MHD simulations to test its influence on α. We find that the region of the cold branch with Σ < 191 g cm^{1} does not maintain the MRIdriven turbulence. Fleming et al. (2000) found in isothermal simulations that the presence of a net vertical magnetic flux, providing external support to the MRI, sustains turbulence down to R_{m} = 100. Net flux may help push the MRIinactive region to slightly lower Σ but the dynamical consequences might be difficult to evaluate without a global model to study the evolution of the net flux. A net magnetic field is also expected to impact substantially the saturated value of α when the corresponding β ≤ 10^{5} (Bai & Stone 2013; Salvesen et al. 2016), that is, B ≥ 20G on the hot branch at a Σ of 174g cm^{2} and B ≥ 4G on the cold branch at a Σ of 93g cm^{2}. Using a dipole approximation, this corresponds to a surface magnetic field of ~10^{4−5}G for the white dwarf, smaller than the typical field of intermediate polars (in which the disk is truncated by the white dwarf magnetosphere), or 10^{2−3}G for the companion, somewhat on the high side of the measured magnetic fields in lowmass stars (see Fig. 3 of Donati & Landstreet 2009). Hence, in principle, the binary components might be able to provide enough net magnetic flux to impact transport even though, in practice, how this is achieved is unknown.
A quiescent disk will not accrete if it enters the MRIinactive region. This is problematic because the Xray emission observed in quiescent DNe requires ongoing accretion onto the white dwarf. The Xray emission region is constrained to the immediate vicinity of the white dwarf by eclipses, as expected from the accretion boundary layer (see Mukai 2017, and references therein). We verified that this emission is unable to selfsustain MRIdriven accretion by photoionizing the disk.
If MRI remains inactive, angular momentum transport may be due to an entirely different mechanism in quiescence. Spiral shocks are often invoked as a means to accrete. They have been proposed to explain the observed anticorrelation between outburst recurrence time and binary mass ratio in some subtypes of dwarf novae (Cannizzo et al. 2012; Menou 2000). However, other parameters, such as the masstransfer rate from the secondary, also depend on the mass ratio and therefore impact accretion, rendering the interpretation of this anticorrelation difficult. Besides, global simulations studying the propagation of spiral shocks in hydrodynamics (Savonije et al. 1994; Lesur et al. 2015; Arzamasskiy & Rafikov 2017) and MHD (Ju et al. 2017) showed that when the ratio H/R is of the order of 10^{2}, as is the case in the cold state (especially in the outer region which is the region of interest as spiral shocks propagate inward), spiral shock dissipate rapidly and do not lead to any accretion in the inner parts of the disk.
The problem may turn out to be a red herring if the disk actually never samples the regime where MRI is inactive. Our results show that accretion is possible on the highΣ part of the cold branch, and this may be enough for the DNe cycle to operate. DIM models of DNe outburst cycles with an αprescription show disk rings follow complex tracks in the (Σ,T) plane when a cooling front propagates through. The ring does not simply cool on a thermal timescale at constant Σ from the tip of the hot branch to the cold branch of the Scurve. Radial heat fluxes cause it to follow a sideways trajectory to a higher Σ on the cold branch (see e.g., Fig. 9 in Menou et al. 1999), and this would ensure it never samples the MRIinactive part of the cold branch. These radial fluxes become important only because the coolingfront size is of order of the vertical height of the disk, a situation that strains the assumption in these calculations that the radial and vertical directions are decoupled. Consequently, there are some uncertainties on how they are taken into account by the DIM and their influence on the light curves (Cannizzo et al. 2010); uncertainties that only a global model of DNe outbursts under MRIdriven transport could resolve.
More prominently, such a global model would also address whether realistic light curves of DNe outbursts can be obtained at all with MRI transport. In this regard, the enhancement of α up to the values required by the DIM is both a major advance; a set back. A set back because the high values of α are limited to the tip of the hot branch, whereas DIM models assumed α to be increased over the whole hot branch. The difference has important consequences on the light curves predicted by the DIM when α(T) is taken from the MRI simulations: Coleman et al. (2016) obtained very jagged light curves that, in our opinion, compare less favorably to the observations than those of the initial DIM (see their Fig. 9). Yet, this is also undisputedly a major advance because this at last provides a physical explanation of the phenomenological change in α between hot and cold state, first noticed over 30 yr ago; and because, after a certain lull, interest in dwarf novae, the source of much of the basic understanding of accretion physics, has been reawakened.
Acknowledgments
We thank the referee for valuable comments and suggestions. N.S. acknowledges financial support from the pole PAGE of the Université Grenoble Alpes. G.L., G.D. and M.F. are grateful to the participants of the KITP 2017 program on Confronting MHD Theories of Accretion Disks with Observations for the many useful conversations pertaining to this work (and thus this research was supported in part by the National Science Foundation under Grant No. NSF 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.
References
 Arzamasskiy, L., & Rafikov, R. R. 2017, ApJ, submitted, ArXiv eprints [arXiv:1710.01304] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2011, ApJ, 736, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 767, 30 [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., & Terquem, C. 2001, ApJ, 552, 235 [Google Scholar]
 Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 [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., & Wheeler, J. C. 1984, ApJS, 55, 367 [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., Still, M. D., Howell, S. B., Wood, M. A., & Smale, A. P. 2010, ApJ, 725, 1393 [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., Yerger, E., Blaes, O., Salvesen, G., & Hirose, S. 2017, MNRAS, 467, 2625 [Google Scholar]
 Donati, J.F., & Landstreet, J. 2009, ARA&A, 47, 333 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Draine, B. T. 2010, Physics of the interstellar and intergalactic medium (Princeton University Press) [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [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]
 Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18 [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]
 Hameury, J.M., Menou, K., Dubus, G., Lasota, J.P., & Hure, J.M. 1998, MNRAS, 298, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., & Stone, J. M. 1998, ApJ, 501, 758 [NASA ADS] [CrossRef] [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]
 Ju, W., Stone, J. M., & Zhu, Z. 2017, ApJ, 841, 29 [NASA ADS] [CrossRef] [Google Scholar]
 King, A. R. 1997, MNRAS, 288, L16 [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]
 Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295 [NASA ADS] [CrossRef] [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., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Menou, K. 2000, Science, 288, 2022 [NASA ADS] [CrossRef] [Google Scholar]
 Menou, K., Hameury, J.M., & Stehle, R. 1999, MNRAS, 305, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A. 2009, Mem. Soc. Astron. Ital. Suppl., 13, 67 [NASA ADS] [Google Scholar]
 Mukai, K. 2017, PASP, 129, 062001 [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., & Stone, J. M. 2002, ApJ, 577, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13 [NASA ADS] [CrossRef] [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]
 Smak, J. 1984, Acta Astron., 34, 161 [NASA ADS] [Google Scholar]
 Turner, N., & Stone, J. 2001, ApJS, 135, 95 [Google Scholar]
 Warner, B. 2003, Cataclysmic variable stars (Cambridge University Press), 28 [Google Scholar]
 Yang, L. T., & Brent, R. P. 2002, Proc. of the Fifth International Conference on Algorithms and Architectures for Parallel Processing, IEEE [Google Scholar]
All Tables
All Figures
Fig. 1
Thermal equilibria in the [T_{mid}] vs. [Σ] (top) or [T_{eff}] vs. [Σ] (bottom) plane. Squares and circles are for periodic and outflow runs, respectively. Triangular dots represent runs with runaway cooling (triangle facing down) or heating (triangle facing up). Error bars represent the standard deviation of the temperature fluctuations. The symbols are colorcoded to the value of α. The colorcoded curves are vertical thermal equilibria using an αprescription. The dashed blue line indicates where the magnetic Reynolds number R_{m} = 10^{4} based on an isothermal model (see Sect. 4.1). 

In the text 
Fig. 2
α as a function of [T_{mid}] (top panel) or [T_{eff}] (bottom panel). Error bars represent the standard deviations of the fluctuations around the mean values. Blue, green, and red colors are for the cold, middle and hot branches, respectively. Squares and circles indicate periodic and outflow runs. The shaded area corresponds to the thermally unstable region. 

In the text 
Fig. 3
Timeaveraged vertical profiles of the radiative, advective, and total flux for different runs. Top panel: run 442P during a highly convective episode. Middle panel: same run but averaging over convective and non convective episodes. Bottom panel: run 434P where there is no convection. The vertical line indicates the height above which magnetic pressure becomes larger than thermal pressure. The dashed line shows N^{2}/ Ω^{2} where N is the BruntVäisälä frequency. 

In the text 
Fig. 4
Evolution of as a function of time for runs located on different parts of the Scurve: (top panel) 434P is a cold branch run, 468P a nonconvective hot branch run, and (bottom panel) 442 a convective hot branch run. For the latter we show the simulation run with outflow (442O) and periodic boundary conditions (442P). We also plot, for run 442O, the contribution of the Maxwell stress to as a dashed line. 

In the text 
Fig. 5
Autocorrelation function f_{auto} as a function of time for the run 442O. τ_{0} is the parabolic decay time of the autocorrelation function (with the parabolic fit shown as a dashed line). τ_{1} is the location of the first peak of f_{auto}. 

In the text 
Fig. 6
Period τ_{1} of the cycles in α as a function of the convective fraction for convectively unstable simulations on the hot branch. 

In the text 
Fig. 7
vs. T_{mid} for two representative convectively unstable simulations: hot branch run 442O (top) and middle branch run 401O (bottom). 

In the text 
Fig. 8
vs. for the same runs as Fig. 7. 

In the text 
Fig. 9
α vs. convective fraction f_{conv} for all convectively unstable runs (same color and symbol coding as Fig. 2). 

In the text 
Fig. 10
Average azimuthal magnetic field ⟨ B_{y} ⟩ _{x,y,z} as a function of time for run 446 (outflow O or periodic P boundary conditions). 

In the text 
Fig. 11
Vertical profile of azimuthal magnetic field ⟨ B_{y} ⟩ _{x,y} (in code units) as a function of time with outflow (top) and periodic (bottom) boundary conditions for run 446 (hot branch). 

In the text 
Fig. 12
Time evolution of ⟨ R_{m} ⟩, ⟨ B_{y} ⟩ and α for run 435PR (Σ = 191 g cm^{2}). The vertical black line marks the time at which the simulation is restarted from 435P with resistivity turned on. 

In the text 
Fig. 13
Same as Fig. 12 for run 462PR (Σ = 174 g cm^{2}). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.