Free Access
Issue
A&A
Volume 609, January 2018
Article Number A77
Number of page(s) 13
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201731900
Published online 11 January 2018

© 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 solar-type star to a white dwarf. Their light curves show periodic outbursts during which the luminosity typically rises by several magnitudes (Warner 2003). According to the disk instability model (DIM; see Lasota 2001, for a review), these outbursts are caused by a thermal-viscous instability due to the steep temperature dependence of the opacity when hydrogen ionizes around 7000 K. During quiescence, the disk fills up as mass is transferred from the companion, gradually heating up until this critical temperature is reached at some radius. This triggers the propagation of heat fronts through the disk, bringing it into a hot state with a high accretion rate. The disk then empties until the temperature falls below the critical temperature and a cooling front returns the disk to quiescence. The timescales involved are set by the ability of the disk to transport angular momentum, providing an observational handle on α. Outburst-decay 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 magneto-rotational instability (MRI; Balbus & Hawley 1991). The properties of MRI-induced transport have been extensively studied over the past 25 years using local, shearing box simulations. Isothermal, stratified local simulations with zero net magnetic flux show a universal value of α ~ 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 non-stratified shearing box, and showed numerically that the disk is indeed thermally unstable in the conditions of DNe. Their zero net flux simulations give α ~ 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 X-rays from the white dwarf boundary layer may play in the ionization balance of quiescent disks.

2. Methods

We adopt the local shearing-box approximation (Hawley et al. 1995) to simulate a vertically-stratified patch of accretion disk situated at a distance R0 = 1.315 × 1010cm from a 0.6 M white dwarf, giving an angular velocity Ω(R0) = 5.931 × 10-3s-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 co-rotating 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 time-stepping following Flock et al. (2013). In this step, we solve the coupled matter-radiation equations in the flux-limited diffusion approximation where ρ is the density, the fluid velocity vector, P the thermal pressure, B the magnetic field vector, Φ = (−3x2 + z2) × Ω2(R0)/2 is the gravitational potential in the co-rotating frame, η the Ohmic resistivity, J = (c/ 4π)∇∇ × B the current density vector, E = ρϵ + 0.5ρvvv2 + B2/ 8π the total energy, ϵ the internal energy, Pt = P + B2/ 8π the thermal pressure plus the magnetic pressure, c the speed of light, aR = (4σ/c) the radiation constant with σ the Stefan-Boltzmann constant, T the temperature, ER the radiation density energy, κP the Planck opacity and κR the Rosseland opacity. The radiative energy flux, in the flux diffusion approximation, is Frad = ((R) /κR(T)ρ)∇ER. The flux limiter is defined as λ(R) ≡ (2 + R)/(6 + 3R + R2) 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 pre-calculated tables and interpolate linearly between values of the tables. Tables are computed from the Saha equations assuming ionization equilibrium for the solar composition of Grevesse & Sauval (1998) with a hydrogen abundance X = 0.7 and a metallicity Z = 0.02. The adiabatic index Γ, the entropy and the thermal capacity Cv are similarly computed.

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

2.2. Boundary conditions

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

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

2.3. Numerical method

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

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

2.4. Initial conditions

The initial flow is Keplerian with a weak zero net flux magnetic field of the form Bz ∝ 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 quasi-steady state. The isothermal temperature Tc0 is set to the mid-plane temperature found by using the vertical structure code of Hameury et al. (1998) for given surface density Σ0 and effective temperature Teff. 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 quasi-steady state (if there is one).

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

Table 1

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 Hcs(Tc0)/Ω is the pressure scale-height (with cs the sound speed). The horizontal extent of the box is of ±6H for the hot branch and ±3H for the cold and middle branch. Lx,Ly and Lz 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 tavg performed once the simulation has reached a quasi-steady state. The averaging timescale tavg 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). [Σ], [Tmid], [Teff], and [τtot] are thus the time-averaged 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, Σ, Teff as and (12)where is the radiative flux in the vertical direction on the upper/lower boundary of our simulation box. We also define Qx,y the total cooling rate as (13)vz the vertical velocity and the quantity ϵvz represents the advective flux of internal energy Fadv. 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 Wxy is ρ(vxvyvAxvAy) where vA 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 Lx, Ly and Lz on the cold branch and on the hot branch. We note n the number of points per final time-averaged scale height hcs( [Tmid] )/Ω. 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 Tmid 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 Lx/H = 1.5, Ly/H = 6.0 and Lz/H = 8.0 and found a difference on [Tmid] 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 two-fold 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 Lx and Ly 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.

thumbnail Fig. 1

Thermal equilibria in the [Tmid] vs. [Σ] (top) or [Teff] 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 color-coded to the value of α. The color-coded curves are vertical thermal equilibria using an α-prescription. The dashed blue line indicates where the magnetic Reynolds number Rm = 104 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 time-averaged midplane temperature [Tmid] ≡ [ ⟨ T(z = 0) ⟩ x,y] versus the surface density [Σ] (top panel). The bottom panel shows [Teff] 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 degree-of-freedom 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 S-curve) with a hot stable branch and a cold stable branch, providing independent confirmation of the results of Hirose et al. (2014). The S-curve from the simulations is comparable to that predicted by the vertical structure calculations using an α-prescription. These also predict a third stable branch at intermediate temperatures (Tmid ≈ 104 K). Indeed, Fig. 1 shows also a set of three stable runs in an intermediate regime of temperature. A middle-branch in radiative MRI simulations was also reported by Hirose (2015) in the context of protoplanetary disks. The middle branch is not as extended as the other branches and may be seen as a prolongation of the cold branch to higher temperatures. In fact, the middle branch does not overlap with the cold branch, so there is at most two equilibrium states for a given Σ, not three. However, the middle branch runs have a higher opacity than the cold branch and are convectively unstable (Sect. 3.2).

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

We have good agreement between our simulations and the α-prescription calculation concerning the equilibrium temperatures at given α and the critical points for runaway heating/cooling. The equilibrium values follow the calculated curves corresponding to the value of α found in the simulation. Somewhat unexpectedly, the agreement extends to the middle branch even though its location in the ,T) plane depends on the chosen value of the convection mixing length parameter (Cannizzo & Wheeler 1984). Here, we took α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 S-curve implies very high accretion rates that are unlikely to be relevant to dwarf novae. The value of α also increases as we approach the tip of the hot branch, going to a maximum of 0.098 for the lowest-density run on the hot branch (453P). This enhancement is clearly seen when plotting α against temperature (Fig. 2, see also Fig. 1 in Coleman et al. 2016). The enhancement is accompanied by stronger fluctuations as we near the unstable tip of the hot branch (see Sect. 3.2.3).

thumbnail Fig. 2

α as a function of [Tmid] (top panel) or [Teff] (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 Brunt-Väisälä frequency N, defined by (15)where Γ ≡ (lnP/lnρ)s is the adiabatic index and s is the specific entropy. N2 has the sign of the entropy gradient along z, thus a negative Brunt-Väisälä frequency denotes a convectively unstable vertical profile. The Brunt-Väisälä frequency gives a linear stability criterion for the convective stability only in regions where the total pressure is dominated by thermal pressure and does not give direct information about the convection flux.

Figure 3 shows the radiative, advective, and total fluxes (Eq. (13)) as well as the Brunt-Väisälä frequency as a function of height for several illustrative cases. The top panel shows a highly convective episode (averaging over 16 local orbits) of the run 442P (hot branch). The advective flux is the main source of heat transport in the regions where the disk is convectively unstable (N2< 0) as opposed to the upper regions where the radiative flux dominates. Averaging the same run on a longer time (80 local orbits) including convective and non-convective episodes shows the contribution of the advective flux decreases (middle panel). We also see evidence for downward transport of heat between the midplane and H/R0 ~ 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 MRI-driven turbulence. Because motions are approximately adiabatic in this region, this mixing tends to flatten the entropy profile, resulting in a downward transport of heat. Thus, the advective 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.

thumbnail Fig. 3

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

thumbnail Fig. 4

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

3.2.3. Convective/radiative cycles

In the following, we make use of the convective fraction fconv 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 S-curve. 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).

thumbnail Fig. 5

Autocorrelation function fauto 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 fauto.

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 fauto (τ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 fauto (τ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 fconv (Fig. 6). On the middle branch, the cycles are much longer for the same fconv. For instance, run 401O has fconv = 0.14 and τ1 ≈ 111 local orbits, requiring a longer averaging timescale tavg than for the other runs.

thumbnail 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 anti-correlated 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 Tmid are the same for both. This difference can be explained by noting that the temperature gradient, assuming only radiative transport, is κ/T3, 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.

thumbnail Fig. 7

vs. Tmid 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 time-averaged α in the middle branch. Figure 9 shows no trend in the convective runs between the time-averaged values of α and fconv. We conclude that there is no simple relationship between the enhancement of α and the relative strength of the convection measured by fconv. 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.

thumbnail Fig. 8

vs. for the same runs as Fig. 7.

thumbnail Fig. 9

α vs. convective fraction fconv 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 box-averaged Byx,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.

thumbnail Fig. 10

Average azimuthal magnetic field Byx,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 Byx,y, the vertical profile of the azimuthal field. Non-convective MRI simulations are known to form “butterfly wings” as the magnetic field flips sign at each quasi-periodic pulsation and propagates outward (e.g. see Fig. 12). Coleman et al. (2017) reported that the quasi-periodic pulsation does not always lead to a sign flip of Byx,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 Byx,y,z = 0. We also see this behavior in the middle branch during convective episodes.

thumbnail Fig. 11

Vertical profile of azimuthal magnetic field Byx,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 electron-neutral collisions increase. This is quantified by the magnetic Reynolds number (18)with η being the resistivity, cs the sound speed and h = cs/ Ω the local pressure scale-height. For Rm < 104, 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 Rm is of the order of 103 in a quiescent DNe disk.

thumbnail Fig. 12

Time evolution of Rm, By 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.

thumbnail Fig. 13

Same as Fig. 12 for run 462PR (Σ = 174 g cm-2).

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 electron-neutral collision. The values of η are pre-computed 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 Δx2/η where Δx is the cell minimum length. To avoid dramatically small time steps, we impose a minimum value of Rm = 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 , Rmx,y, and Byx,y for those two runs before and after resistivity is switched on.

Turning on resistivity has no influence on run 435PR (Fig. 12). Rm stays well above the critical value of 104 in the midplane. Some regions above the photosphere have a lower Rm, 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 ([Tmid] ≈ 3100 K compared to 3750 K for 435P). Here, Rm starts to drop as soon as resistivity is switched on (Fig. 13). The turbulence is suppressed once Rm decreases below 5000 and the transport of angular momentum ceases. Byx,y diffuses as expected: first near the photosphere and then, as Rm 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 non-ideal effect. The typical values deduced from this run for the density, temperature, ionization fraction and magnetisation vA/cs are 5 × 10-7g cm-3, 3000 K, 10-5.5 and 0.1, respectively. Therefore, in this run, ambipolar diffusion is expected to have a negligible impact while the Hall effect is only an order of magnitude weaker than Ohmic diffusion. We note that these conclusions also hold as the disc cools down when MRI is suppressed.

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

Ambipolar diffusion is characterized by the ambipolar Elsasser number (Balbus & Terquem 2001): (22)where ρi = nimi is the ion density and γi = ⟨ σvni/ (mn + mi) = 2.7 × 1013(41.33mu/mn + mi)cm3 g-1 s-1 is the ion-neutral drag coefficient. Following Balbus & Terquem (2001), we assume a momentum rate coefficient for ion-neutral collision σvni = 1.9 × 10-9cm3 s-1 (Draine 2010), mn = 2.33mu and mi = 39mu (the neutral and ion masses respectively). We obtain ΛA ~ 103 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 two-fold. For strong Hall effect ( ≲ 5), turbulent transport is essentially suppressed and the flow relaxes into a self-organized 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 103 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 ~ 103, 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/Rm does not depend on the temperature, meaning that at the temperature where H ≲ 50, one also expects Rm ≲ 50, well below the critical Rm found in Hall-MRI simulation (Sano & Stone 2002). We conclude that the Hall effect does not affect the critical temperature (or Σ) below which DNe are MRI-stable, nor does it revive the MRI at lower temperatures.

4.3. Influence of X-ray irradiation

The ionization fraction xe 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) xe ≈ 4 × 10-6 in the MRI-stable run 462PR (at midplane) and xe ≈ 6 × 10-6 in the MRI-unstable 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 X-ray emission originating from the accretion boundary layer onto the white dwarf that may provide sufficient ionization to maintain the MRI active in the quiescent disk.

The quiescent X-ray spectra are bremsstrahlung with temperatures kTXR ranging from 1 to 10 keV and luminosities LXR ranging from 1028 to 1032 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 ~103 closer to the X-ray source, so the impinging ionizing flux is much higher.

We calculate the expected ionization fraction xe 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 X-ray emission region is about the size of the white dwarf (Mukai 2017), which is much greater than the disk height. In this case, we may expect (King 1997) (25)with a white dwarf radius RWD = 109 cm.

We obtain xe by solving Eq. (11) from Fromang et al. (2002), which relates xe to the fraction of metals xM = nM/nn, the dissociative recombination rate for molecular ions, the radiative recombination rate for metal atoms and the rate of charge exchange between them. We further assume that the disk is passive, that is, has no internal heating, and the only source of heat is X-ray irradiation. Hence, in steady-state, the disk is isothermal at the temperature set by . The vertical structure is in hydrostatic equilibrium with Σ = 100 g cm-2 to compare to our results with the first MRI-stable run, 462PR.

We find that X-ray ionization is negligible, in agreement with the order-of-magnitude estimate of Gammie & Menou (1998). For example, with and kTXR = 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 × 1017 cm-3, giving xe = 4.1 × 10-10 when no metals are present (xM = 0) and xe = 1.1 × 10-7 when xM = 6.86 × 10-6 (solar abundance). The fraction xe decreases for lower values of or kTXR. Hence, X-ray ionization is unable to provide the required ionization fraction except in the upper layer. However, only a small amount of material is active in this layer: taking the values assumed above, xe > 5 × 10-6 for z > 1.5 × 108 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 well-known S-curve 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 fconv, 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 MRI-driven turbulent transport in the cold, quiescent state of DNe (Gammie & Menou 1998). We verified that Ohmic dissipation is the dominant non-ideal effect and carried out resistive MHD simulations to test its influence on α. We find that the region of the cold branch with Σ < 191 g cm-1 does not maintain the MRI-driven turbulence. Fleming et al. (2000) found in isothermal simulations that the presence of a net vertical magnetic flux, providing external support to the MRI, sustains turbulence down to Rm = 100. Net flux may help push the MRI-inactive 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 β ≤ 105 (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 ~104−5G for the white dwarf, smaller than the typical field of intermediate polars (in which the disk is truncated by the white dwarf magnetosphere), or 102−3G for the companion, somewhat on the high side of the measured magnetic fields in low-mass stars (see Fig. 3 of Donati & Landstreet 2009). Hence, in principle, the binary components might be able to provide enough net magnetic flux to impact transport even though, in practice, how this is achieved is unknown.

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

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

The problem may turn out to be a red herring if the disk actually never samples the regime where MRI is inactive. Our results show that accretion is possible on the high-Σ 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 S-curve. 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 MRI-inactive part of the cold branch. These radial fluxes become important only because the cooling-front size is of order of the vertical height of the disk, a situation that strains the assumption in these calculations that the radial and vertical directions are decoupled. Consequently, there are some uncertainties on how they are taken into account by the DIM and their influence on the light curves (Cannizzo et al. 2010); uncertainties that only a global model of DNe outbursts under MRI-driven transport could resolve.

More prominently, such a global model would also address whether realistic light curves of DNe outbursts can be obtained at all with MRI transport. In this regard, the enhancement of α 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 re-awakened.

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 PHY-1125915). This work was granted access to the HPC resources of IDRIS under the allocation A0020402231 made by GENCI (Grand Equipement National de Calcul Intensif). Some of the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

References

  1. Arzamasskiy, L., & Rafikov, R. R. 2017, ApJ, submitted, ArXiv e-prints [arXiv:1710.01304] [Google Scholar]
  2. Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bai, X.-N., & Stone, J. M. 2013, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
  4. Balay, S., Abhyankar, S., Adams, M. F., et al. 2016, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]
  5. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
  6. Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235 [Google Scholar]
  7. Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 [NASA ADS] [CrossRef] [Google Scholar]
  8. Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] [Google Scholar]
  9. Byckling, K., Mukai, K., Thorstensen, J., & Osborne, J. 2010, MNRAS, 408, 2298 [NASA ADS] [CrossRef] [Google Scholar]
  10. Cannizzo, J. K., & Wheeler, J. C. 1984, ApJS, 55, 367 [NASA ADS] [CrossRef] [Google Scholar]
  11. Cannizzo, J. K., Shafter, A. W., & Wheeler, J. C. 1988, ApJ, 333, 227 [NASA ADS] [CrossRef] [Google Scholar]
  12. Cannizzo, J. K., Still, M. D., Howell, S. B., Wood, M. A., & Smale, A. P. 2010, ApJ, 725, 1393 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cannizzo, J. K., Smale, A. P., Wood, M. A., Still, M. D., & Howell, S. B. 2012, ApJ, 747, 117 [NASA ADS] [CrossRef] [Google Scholar]
  14. Coleman, M. S. B., Kotko, I., Blaes, O., Lasota, J.-P., & Hirose, S. 2016, MNRAS, 462, 3710 [NASA ADS] [CrossRef] [Google Scholar]
  15. Coleman, M. S., Yerger, E., Blaes, O., Salvesen, G., & Hirose, S. 2017, MNRAS, 467, 2625 [Google Scholar]
  16. Donati, J.-F., & Landstreet, J. 2009, ARA&A, 47, 333 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  17. Draine, B. T. 2010, Physics of the interstellar and intergalactic medium (Princeton University Press) [Google Scholar]
  18. Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
  19. Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464 [NASA ADS] [CrossRef] [Google Scholar]
  20. Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gammie, C. F., & Menou, K. 1998, ApJ, 492, L75 [NASA ADS] [CrossRef] [Google Scholar]
  23. Glassgold, A. E., Najita, J., & Igea, J. 1997a, ApJ, 480, 344 [NASA ADS] [CrossRef] [Google Scholar]
  24. Glassgold, A. E., Najita, J., & Igea, J. 1997b, ApJ, 485, 920 [NASA ADS] [CrossRef] [Google Scholar]
  25. Grevesse, N., & Sauval, A. 1998, Space Sci. Rev., 85, 161 [NASA ADS] [CrossRef] [Google Scholar]
  26. Hameury, J.-M., Menou, K., Dubus, G., Lasota, J.-P., & Hure, J.-M. 1998, MNRAS, 298, 1048 [NASA ADS] [CrossRef] [Google Scholar]
  27. Hawley, J. F., & Stone, J. M. 1998, ApJ, 501, 758 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hirose, S. 2015, MNRAS, 448, 3105 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., & Sano, T. 2014, ApJ, 787, 1 [NASA ADS] [CrossRef] [Google Scholar]
  32. Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
  33. Ju, W., Stone, J. M., & Zhu, Z. 2017, ApJ, 841, 29 [NASA ADS] [CrossRef] [Google Scholar]
  34. King, A. R. 1997, MNRAS, 288, L16 [NASA ADS] [CrossRef] [Google Scholar]
  35. King, A., Pringle, J., & Livio, M. 2007, MNRAS, 376, 1740 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kotko, I., & Lasota, J.-P. 2012, A&A, 545, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295 [NASA ADS] [CrossRef] [Google Scholar]
  38. Lasota, J.-P. 2001, New Astron. Rev., 45, 449 [NASA ADS] [CrossRef] [Google Scholar]
  39. Latter, H. N., & Papaloizou, J. C. 2012, MNRAS, 426, 1107 [NASA ADS] [CrossRef] [Google Scholar]
  40. Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Menou, K. 2000, Science, 288, 2022 [NASA ADS] [CrossRef] [Google Scholar]
  42. Menou, K., Hameury, J.-M., & Stehle, R. 1999, MNRAS, 305, 79 [NASA ADS] [CrossRef] [Google Scholar]
  43. Mignone, A. 2009, Mem. Soc. Astron. Ital. Suppl., 13, 67 [NASA ADS] [Google Scholar]
  44. Mukai, K. 2017, PASP, 129, 062001 [NASA ADS] [CrossRef] [Google Scholar]
  45. Ryan, B. R., Gammie, C. F., Fromang, S., & Kestener, P. 2017, ApJ, 840, 6 [NASA ADS] [CrossRef] [Google Scholar]
  46. Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016, MNRAS, 457, 857 [NASA ADS] [CrossRef] [Google Scholar]
  47. Sano, T., & Stone, J. M. 2002, ApJ, 577, 534 [NASA ADS] [CrossRef] [Google Scholar]
  48. Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13 [NASA ADS] [CrossRef] [Google Scholar]
  49. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  50. Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685 [NASA ADS] [CrossRef] [Google Scholar]
  51. Smak, J. 1984, Acta Astron., 34, 161 [NASA ADS] [Google Scholar]
  52. Turner, N., & Stone, J. 2001, ApJS, 135, 95 [Google Scholar]
  53. Warner, B. 2003, Cataclysmic variable stars (Cambridge University Press), 28 [Google Scholar]
  54. 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

Table 1

Initial parameters and results for our ideal MHD simulations.

All Figures

thumbnail Fig. 1

Thermal equilibria in the [Tmid] vs. [Σ] (top) or [Teff] 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 color-coded to the value of α. The color-coded curves are vertical thermal equilibria using an α-prescription. The dashed blue line indicates where the magnetic Reynolds number Rm = 104 based on an isothermal model (see Sect. 4.1).

In the text
thumbnail Fig. 2

α as a function of [Tmid] (top panel) or [Teff] (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
thumbnail Fig. 3

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

In the text
thumbnail Fig. 4

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

In the text
thumbnail Fig. 5

Autocorrelation function fauto 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 fauto.

In the text
thumbnail 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
thumbnail Fig. 7

vs. Tmid for two representative convectively unstable simulations: hot branch run 442O (top) and middle branch run 401O (bottom).

In the text
thumbnail Fig. 8

vs. for the same runs as Fig. 7.

In the text
thumbnail Fig. 9

α vs. convective fraction fconv for all convectively unstable runs (same color and symbol coding as Fig. 2).

In the text
thumbnail Fig. 10

Average azimuthal magnetic field Byx,y,z as a function of time for run 446 (outflow O or periodic P boundary conditions).

In the text
thumbnail Fig. 11

Vertical profile of azimuthal magnetic field Byx,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
thumbnail Fig. 12

Time evolution of Rm, By 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
thumbnail 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.