© N. Scepi et al. 2020

1. Introduction

Magnetic fields in accretion disks are largely unconstrained with regard to observations. However, the magnetic configuration has tremendous importance in achieving an understanding of the observed properties and evolution of astrophysical disks. It is known that active galactic nuclei (AGN; Lister et al. 2009; Blandford et al. 2019), protoplanetary disks (Burrows et al. 1996; Ray et al. 1996; Hirth et al. 1997; Dougados et al. 2000), and low-mass X-ray binaries (LMXBs; Mirabel et al. 1992; Fender et al. 1997) demonstrate the presence of strongly collimated, fast outflows called jets. It was recently proposed that jets may also be present in dwarf novæ (DNe; Körding et al. 2008; Russell et al. 2016; Coppejans et al. 2016; Coppejans & Knigge 2020). The source of these jets may be a strong (near equipartition) large-scale poloidal magnetic field threading the disk that magneto-centrifugally accelerates and collimates the flow (Blandford & Payne 1982; Pudritz & Norman 1983; Königl 1989; Ferreira & Pelletier 1995). More extended, slower ejections, called winds, have also been observed in all the accretion disks cited above (Cordova & Mason 1982; Mauche & Raymond 1987; Crenshaw et al. 2003; Miller et al. 2004; Ponti et al. 2012; Louvet et al. 2018).

Recently, 3D MRI numerical simulations (Fromang et al. 2013; Bai & Stone 2013; Lesur et al. 2013) and self-similar solutions (Jacquemin-Ide et al. 2019) showed that outflows, as per Blandford & Payne (1982), could be extended to lower magnetization to produce slower, more massive solutions compared to the previously computed jet solutions (Ferreira 1997). However, it is not clear that these lower magnetization outflows could fully explain the observed features of winds. Indeed, at least in DNe and LMXBs, it seems that to reproduce the estimated mass loss rates and terminal velocities from observations, there should be some hybrid mechanism involving magnetic and thermal driving or radiative driving at work (Proga 2003; Miller et al. 2006; Luketic et al. 2010; Higginbottom & Proga 2015; Chakravorty et al. 2016; Díaz Trigo & Boirin 2016; Jacquemin-Ide et al. 2019; Trueba et al. 2019). Nonetheless, both jets and winds in accretion disks strongly depend on the presence of a large-scale magnetic field threading the disk.

A large-scale magnetic field does not only impact the ejection mechanisms but also the accretion processes. In fact, the magnetic field is believed to be one of the main drivers of accretion, either through turbulence due to the magneto-rotational instability (MRI; Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991)¹, which transports angular momentum radially, or through magnetic winds or jets (Blandford & Payne 1982; Pudritz & Norman 1983; Königl 1989; Ferreira & Pelletier 1995), which extract angular momentum vertically. Moreover, both the strength of MRI (Hawley et al. 1995; Salvesen et al. 2016; Scepi et al. 2018a) and the properties of the magnetic outflows (Jacquemin-Ide et al. 2019) depend on the large-scale poloidal magnetic field strength.

Large-scale poloidal magnetic fields can have two origins. They could either be generated in situ by dynamo effects (Brandenburg et al. 1995; Tout & Pringle 1996) or could be advected inwards from the surrounding medium (Bisnovatyi-Kogan & Ruzmaikin 1974). It is still not clear whether MRI dynamo can produce an organized poloidal magnetic field on large scales (Brandenburg et al. 1995). Local shearing box simulations tend to show that they do not but state-of-the-art 3D global simulations from Liska et al. (2020) found that loops of poloidal magnetic field, when formed in the outer region and dragged in the inner region of the disk, become “large scale” due to the difference in length scale between the outer and inner region.

In a seminal work, Lubow et al. (1994a) studied the transport of magnetic field in a thin disk where angular momentum transport is only due to turbulence, with turbulence prescribed as an α parameter, where α is the ratio of the radial stress to the thermal pressure (Shakura & Sunyaev 1973). They showed that in an α-disk, it appears impossible to advect a large scale magnetic field inwards. These results were later confirmed by Heyvaerts et al. (1996). Indeed, with a magnetic Prandtl number (ratio of the turbulent viscosity to the magnetic turbulent diffusivity) of unity, as seems to be the case for MRI-driven turbulence (Lesur & Longaretti 2009), the large-scale poloidal magnetic field efficiently diffuses outwards. Lubow et al. (1994a) considered the advection of magnetic field inwards by turbulent accretion and computed the diffusion term by solving the global magnetic structure in the potential field approximation. However, the authors ignored the vertical structure of the disk and assumed that the velocity at which the magnetic flux is advected is the same as the velocity at which the mass is accreted.

Bisnovatyi-Kogan & Lovelace (2007) and Rothstein & Lovelace (2008) proposed that taking into account the vertical structure of the disk could help advect magnetic field inwards. Their idea is based on the fact that the upper layers of the disk, where the MRI is suppressed, have a higher effective conductivity than the bulk of the disk. Although accretion in the upper layers does not play a large role in mass transport because the density is low in those regions, it could help advect the magnetic field inwards faster than it would in the more diffusive midplane of the disk. Evidence for this mechanism is also found in 3D global simulations (Beckwith et al. 2009; Tchekhovskoy et al. 2011; Suzuki & Inutsuka 2014; Zhu & Stone 2018) where funnels of low-density matter above the disk carry the magnetic flux inwards. Guilet & Ogilvie (2012) computed the vertical structure of a thin disk through an asymptotic expansion of the resistive MHD equations. With this semi-local model, they were able to compute the transport rates due to an α turbulent accretion as well as the MHD outflow-driven accretion as a function of magnetization. They showed, in the context of protoplanetary disks, that by taking account the vertical structure of the disk, the magnetic field could be efficiently advected inward for a realistic magnetic Prandtl number due to accretion in the upper low-density layers (Guilet & Ogilvie 2014).

The problem of inefficient dragging of the magnetic field is fundamentally due to the fact that the turbulent magnetic diffusivity is on the same order as that of the turbulent viscosity (Tout & Pringle 1996). This problem could be alleviated if the angular momentum is removed by the magnetic outflow-driven torque, which can be shown to be more efficient than the turbulent torque (Lubow et al. 1994b; Heyvaerts et al. 1996; Ogilvie & Livio 2001; Guilet & Ogilvie 2012). Li & Cao (2019) constructed a global model of thin disks where the removal of angular momentum is done by a magnetic outflow and the magnetic structure is evolved self-consistently. To evolve the magnetic field, they computed the global magnetic structure in the same way as Lubow et al. (1994a). They showed that magnetic field is advected inwards very efficiently and forms a very magnetized zone around the central object even when they start from a relatively moderate external magnetic field. Their work focused on the stationary, global structure of the disk and the outflows.

In this paper, we wish to explore the dynamical evolution of a disk threaded throughout by a large-scale poloidal magnetic field in the context of transient compact objects. We focus, in particular, on the case of dwarf novæ, which are binary systems composed of a white dwarf and a solar type star, where the accretion disk surrounding the white dwarf undergoes thermal-viscous eruptions. Eruptions of dwarf novæ have been modeled for years based on the disk instability model (DIM, Lasota 2001). This model traditionally assumes that the transport of angular momentum is due to turbulence. Dwarf novæ have been used as a testbed to deduce the value of α from observations (Kotko & Lasota 2012; Cannizzo et al. 2012). However, despite many efforts, no physically satisfying solution has been found to produce the values of α that are required with regard to observations (Hirose et al. 2014; Coleman et al. 2016; Scepi et al. 2018b). Recently, Scepi et al. (2019) proposed that the torque produced by a magnetic outflow could help resolve this issue. They showed that the inclusion of the outflow-driven magnetic torque shapes the light curves. In particular, they were able to produce realistic eruptions by assuming the magnetic field B is a function of disk radius B ∝ R⁻³. However, Scepi et al. (2019) fixed the magnetic configuration, whereas it is expected to evolve during the outburst as the magnetic field is advected and becomes diffused in the disk. Here, we wish to take this into account by self-consistently evolving the large-scale magnetic configuration together with the modified disk instability model of Scepi et al. (2019). We emphasize that we do not aim to directly compare our results with observations of DNe but, rather, to provide a first theoretical study of the interplay between magnetic field transport and the eruptions of DNe. We also stress that we are not modeling the magnetospheric interaction of a possible white dwarf dipolar field with the disk. This is usually captured as a truncation of the inner disk at the magnetospheric radius, the radius below which the dipolar field is strong enough to funnel accretion (Frank et al. 2002). Such a truncation has been proposed, in combination with the DIM, to explain the high X-ray flux in quiescence and the lag between the outburst rise in UV and optical in DNe (Livio & Pringle 1992; Schreiber et al. 2003).

In Sect. 2 we expose in detail the method we use to evolve the fields in our model. In Sects. 3.1 and 3.2 we explore the behavior of the disk using two different prescriptions for the magnetic field evolution. We show the impact of this choice on the light curves of DNe. Then in Sect. 4 we discuss the caveats of our model, its consequences for observations, and its application to LMXBs. We present our conclusions in Sect. 5.

2. Methods

In this paper, we consider an axisymetric accretion disk around a white dwarf of 0.6 M_⊙. To study the evolution of the ionized plasma that constitutes the accretion disk, we use the viscous and resistive non-relativistic MHD equations in a frame centered on the white dwarf, along with cylindrical coordinates.

2.1. Induction equation

In an axisymetric disk, we can write the magnetic field B as

$$\begin{array}{c}\hfill \mathit{B}=-\frac{1}{R}{\mathit{e}}_{\varphi }\times \mathrm{\nabla }\psi +B_{\varphi }{\mathit{e}}_{\varphi },\end{array}$$(1)

where ψ is the magnetic flux function and R the cylindrical radius in the disk. The induction equation is written as

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\mathit{v}\times \mathit{B}-\eta \mathrm{\nabla }\times \mathit{B}),\end{array}$$(2)

where v is the velocity and η is the magnetic diffusivity, assumed here to be of turbulent origin.

We can rewrite Eq. (2) using Eq. (1) to find

$$\begin{array}{c}\hfill \frac{\partial \psi }{\partial t}=R[v_z{B}_{\mathrm{R}}-v_{\mathrm{R}}{B}_z-\eta (\frac{\partial B_{\mathrm{R}}}{\partial z}-\frac{\partial B_z}{\partial R})]=0.\end{array}$$(3)

Following the asymptotic expansion in a thin disk of the viscous and resistive MHD equations developed by Guilet & Ogilvie (2012), we write

$$\begin{array}{c}\hfill \psi (R,z,t)={\psi }_0(R,t)+{ϵ}^2{\psi }_2(R,z,t)+o({ϵ}^2),\end{array}$$(4)

where ϵ is a small dimensionless parameter of the order of H/R and H ≡ c_s/Ω is the scale height of the disk, c_s the sound speed, and Ω the angular velocity. According to Guilet & Ogilvie (2012), we assume

$$\begin{array}{c}\hfill \frac{v_z}{v_{\mathrm{R}}}=O(ϵ);\frac{B_{\mathrm{R}}}{B_z}=O(ϵ).\end{array}$$(5)

We can check that indeed

$$\begin{array}{c}\hfill \frac{B_{\mathrm{R}}}{B_z}=\frac{-\frac{1}{R}{\partial }_z({ϵ}^2{\psi }_2)}{\frac{1}{R}{\partial }_{\mathrm{R}}{\psi }_0}=O(ϵ)\end{array}$$

since z/R = O(ϵ) in the disk. We discuss the limitations of these assumptions in Sect. 4.1.

Using Eqs. (1)–(5) we can rewrite the poloidal part of the induction equation in leading order as

$$\begin{array}{c}\hfill \frac{\partial {\psi }_0}{\partial t}+v_{\mathrm{R}}\frac{\partial {\psi }_0}{\partial R}+\eta R(\frac{\partial B_{\mathrm{R}}}{\partial z}-\frac{\partial B_z}{\partial R})=0.\end{array}$$(6)

We then vertically integrate Eq. (6) weighted by the conductivity 1/η (Ogilvie & Livio 2001) to find

$$\begin{array}{c}\hfill \frac{\partial {\psi }_0}{\partial t}+R{B}_z[{\overline{v}}_{\mathrm{R}}+\overline{\eta }(\frac{1}{L_z}\frac{B_{\mathrm{Rs}}}{B_z}-\frac{\partial ln{B}_z}{\partial R})]=0\end{array}$$(7)

with

$$\begin{array}{c}\hfill \frac{1}{\overline{\eta }}=\frac{1}{2L_z}\int \frac{1}{\eta }\mathrm{d}z,\end{array}$$

$$\begin{array}{c}\hfill {\overline{v}}_{\mathrm{R}}=\frac{\overline{\eta }}{2L_z}\int \frac{v_{\mathrm{R}}}{\eta }\mathrm{d}z,\end{array}$$

with L_z as the height up to which we integrate in the vertical direction and B_Rs as the radial magnetic field at the disk surface. We assumed an odd symmetry for the radial magnetic field and integrated all quantities between −L_z and L_z. Equation (7) can be rewritten as

$$\begin{array}{c}\hfill \frac{\partial \psi }{\partial t}+v_{\psi }\frac{\partial \psi }{\partial R}=0,\end{array}$$(8)

where

$$\begin{array}{c}\hfill v_{\psi }\equiv {\overline{v}}_{\mathrm{R}}+\overline{\eta }(\frac{1}{L_z}\frac{B_{\mathrm{Rs}}}{B_z}-\frac{\partial ln{B}_z}{\partial R})\end{array}$$(9)

is the transport velocity of the magnetic flux. We dropped the subscript 0 in the magnetic flux function. From Eq. (8) it appears that the induction equation has become a simple advection problem. However, we should keep in mind that the velocity v_ψ depends on ${\partial }_{\mathrm{R}}{B}_z={\partial }_{\mathrm{R}}\left(\frac{1}{R}\frac{\partial \psi }{\partial R}\right)$ and on $B_{\mathrm{Rs}}=-\frac{1}{R}{\partial }_z{\psi |}_{\mathrm{s}}$, so that Eq. (8) describes a advection/diffusion problem.

2.2. Prescriptions for v_ψ

In the absence of a prescription for v_ψ derived from first principles, we assume that magnetic flux transport is mainly local and can be prescribed as a function of β, with

$$\begin{array}{c}\hfill \beta \equiv \frac{8\pi P_{\mathrm{mid}}}{B_z^2},\end{array}$$

where P_mid is the mid-plane thermal pressure. We write β(z) as $8\pi P(z)/B_z^2$, otherwise β is always computed at the mid-plane. We discuss the limitations of our local approach in Sect. 4.1.

With this model, we aim to study the impact on the dynamics of the disk of different prescriptions for v_ψ. We also tried to use the more sophisticated linear model of Guilet & Ogilvie (2012) to obtain a prescription for v_ψ but we encountered some issues with this approach (see Appendix A).

We decompose v_ψ as

$$\begin{array}{c}\hfill v_{\psi }=-v_{\mathrm{in}}+v_{\mathrm{out}}-v_{D_{B_z}}{D}_{B_z},\end{array}$$(10)

where v_in represents the inwards transport velocity of magnetic flux as the flux is dragged inwards with the accretion flow. In our case, we expect accretion to be due to the turbulent torque or a surface torque due to a magnetized outflow. Inversely, v_out represents the outward transport velocity of magnetic flux. It is expected to be dominated by the vertical diffusion of the radial magnetic field. Finally, v_DBzD_Bz represents the radial diffusion of the vertical magnetic field, with $D_{B_z}\equiv \frac{\partial ln{B}_z}{\partial lnR}$.

For simplicity, we impose v_in and v_out as power laws of β so as to explore a wide range of regime for magnetic flux transport. To reduce the degrees of liberty concerning the choice of v_in and v_out, here we give some expected order of magnitudes for the three velocities.

2.2.1. Constraints on v_DBz

Following our simple approach, we set

$$\begin{array}{c}\hfill v_{D_{B_z}}{D}_{B_z}\equiv \overline{\eta }\frac{\partial ln{B}_z}{\partial R}·\end{array}$$

Assuming that η(z) is dominated by the effective turbulent resistivity and is constant with height, this can be rewritten as

$$\begin{array}{c}\hfill v_{D_{B_z}}=\frac{\alpha }{\mathcal{P}}\frac{H}{R}{c}_{\mathrm{s}},\end{array}$$

where

$$\begin{array}{c}\hfill \alpha =\min (1,15\times {\beta }^{-0.56}+0.03)\end{array}$$(11)

is taken from the 3D shearing box simulations with radiative transfer of Scepi et al. (2018a), adding a cap in the low β limit to mimic the fact that MRI-driven turbulence is quenched when β <  1. We also introduce the magnetic Prandtl number, 𝒫, the ratio of the effective turbulent viscosity to the effective turbulent resistivity, which is set to unity in agreement with Lesur & Longaretti (2009).

Our expression for v_DBz is a crude estimation since it does not take into account the impact of the radial gradient of vertical magnetic field on v_R and assumes a constant resistivity.

2.2.2. Constraints on v_out

From Eq. (9) we can estimate

$$\begin{array}{c}\hfill v_{\mathrm{out}}\sim \frac{\alpha }{\mathcal{P}}\frac{H}{L_z}\frac{B_{\mathrm{Rs}}}{B_z}{c}_{\mathrm{s}}.\end{array}$$(12)

Here, L_z represents the surface of the disk and we can crudely assume that it is of the order of ∼H. Constraining the radial surface magnetic field using the self-similar solutions for disk-magnetized outflow from Jacquemin-Ide et al. (2019) in the regime 1 <  β <  10⁴, we have

$$\begin{array}{c}\hfill \frac{B_{\mathrm{Rs}}}{B_z}=0.62\times {\beta }^{0.4}.\end{array}$$(13)

Using Eqs. (11)–(13) we see that

$$\begin{array}{c}\hfill v_{\mathrm{out}}\sim c_{\mathrm{s}}\end{array}$$

when β ≈ 1.

We note that Eq. (13) is not in agreement with our assumption that B_R/B_z = O(ϵ). We discuss the implication of this in Sect. 4.1.

2.2.3. Constraints on v_in

It is more difficult to obtain a simple order of estimation for v_in as we do not know the vertical structure of v_R a priori. In the model from Guilet & Ogilvie (2012), the transport velocity of magnetic flux inwards, v_in, due to either the turbulent torque or the surface torque from a magnetized outflow, tends towards the mass transport velocity, ${\stackrel{\sim }{v}}_{\mathrm{R}}$, in the highly magnetized case (see their Fig. 8). Indeed, in their model, the height where the surface accretion occurs tends towards the midplane as β decreases and the magnetic flux transport velocity tends towards the mass transport velocity. In the highly magnetized regime (β ≈ 1), the mass transport velocity is dominated by the torque from the magnetized outflow (Ferreira & Pelletier 1995) with

$$\begin{array}{c}\hfill {\stackrel{\sim }{v}}_{\mathrm{R}}=\frac{4q}{\beta }{c}_{\mathrm{s}},\end{array}$$

where q ≡ B_ϕs/B_z is given by the self-similar solutions of Jacquemin-Ide et al. (2019)

$$\begin{array}{c}\hfill q=0.36\times {\beta }^{0.6}.\end{array}$$(14)

We see that when β ≈ 1, the mass transport velocity, ${\stackrel{\sim }{v}}_{\mathrm{R}}$, approaches unity and so

$$\begin{array}{c}\hfill v_{\mathrm{in}}\sim c_{\mathrm{s}}.\end{array}$$

2.2.4. Behavior at large β

We expect v_out to represent the vertical diffusion of the radial magnetic field, where the surface radial magnetic field is set by a magnetized outflow. We argue here that an MHD outflow cannot exist at low magnetization and that subsequently all magnetic flux transport associated with an MHD outflow should decrease with β.

The base of the outflow (or surface of the disk) can be roughly estimated as the height, ζ_B, where β(z) = 1 (Suzuki & Inutsuka 2009; Fromang et al. 2013; Zhu & Stone 2018). So, as β increases, ζ_B increases too. However, in a real disk, the base of the outflow, where angular momentum is removed from the disk, must originate from a resistive region where the matter is able to decouple from the magnetic field to be accreted (Ferreira & Pelletier 1995). In a realistic model, the effective turbulent diffusivity should decrease with height as it is set by the turbulent profile in the disk. Hence, in a real disk, a critical value of β above which an outflow cannot be launched should be established. This critical β corresponds to the magnetization where ζ_B becomes larger than the vertical extent of the resistive region.

In our estimations, we assumed that the turbulent resistivity, η, does not depend on z and so, we do not capture this effect. However, the expected limited vertical extent of the resistive region implies that v_out must decrease with β and be null for high β.

The case of v_in is more complicated since we expect inward advection to be due to both outflow-driven accretion and MRI turbulent driven accretion. In contrast to outflow-driven accretion, turbulent MRI-driven accretion is not expected to fade out as β goes to infinity because of the MRI dynamo (Scepi et al. 2018a). So, if we consider that turbulent driven accretion efficiently advects magnetic flux inward as in Guilet & Ogilvie (2014), v_in does not have to be null at high β. However, if we assume that turbulent driven accretion does not efficiently advect magnetic flux inward as in Lubow et al. (1994a), v_in does have to be null at high β. We explore both options in Sects. 3.1 and 3.2 respectively.

In summary, we have three constraints to build our prescriptions of v_ψ. First, v_in and v_out should be on the order of c_s when β approaches unity. Second, v_DBz is smaller than v_in and v_out by a factor H/R when β approaches unity. Third, v_out (and v_in if we only consider accretion due to a magnetized outflow) decreases when increasing β and goes to zero as β tends to infinity.

2.2.5. A variety of prescriptions

When building prescriptions for v_in and v_out as a function of β, we can set up two different scenarios:

– v_in and v_out are never equal and one dominates the other for every β (cases 1 and 2 in Fig. 1),

Fig. 1. Sketch of the four set of prescriptions obeying the constraints on vin and vout from Sect. 2.2 with a comment on the resulting light curves and disk configuration. Green lines represent the outward transport and blue lines the inward transport.

– or v_in and v_out are equal at some β that we call β_eq (cases 3 and 4 in Fig. 1).

We find that the first two sub-cases are actually of low interest. Indeed, in these cases, the magnetic field is accumulated either in the inner parts (case 1 in Fig. 1) or the outer parts (case 2 in Fig. 1) over a very small range in radii². The range of radii over which the flux accumulates depends on the competition between v_in (or v_out) and v_DBzD_Bz. Since v_DBz is smaller than v_in (or v_out) by a factor of ≈H/R in the magnetized regime, it’s necessary to build a large radial gradient for B_z before reaching a magnetic equilibrium. This explains how the field ends up concentrated within the innermost (outermost) radii only. The disk then becomes unmagnetized, except for a highly magnetized inner (outer) zone. The highly magnetized zone does not evolve with time and the rest of the disk has small eruptions as expected in disks with constant α.

The case where an equilibrium magnetization, β_eq, exists, where the inner transport of magnetic flux is balanced by the outer transport of magnetic flux is of greater interest. We can imagine more complicated scenarios where v_in and v_out cross several times. For simplicity, we omit such cases from our discussion.

Once we have established that there is a β_eq, again we face two choices.

• v_out >  v_in for β >  β_eq (case 3 in Fig. 1),

• or v_in >  v_out for β >  β_eq (case 4 in Fig. 1).

As discussed in Appendix C, we believe that case 4, where v_in >  v_out for β >  β_eq, is the astrophysically interesting case and we mainly discuss this case. A discussion on case 3 is added in Appendix C.

2.3. Modified DIM code with magnetic field transport

To evolve the radial structure of the disk with time we solve the vertically averaged equations of a thin disk:

$$\begin{array}{cc}& \frac{\partial \mathrm{\Sigma }}{\partial t}=\frac{1}{2\pi R}\frac{\partial \dot{M}}{\partial R}\hfill \end{array}$$(15)

$$\begin{array}{cc}& \frac{\partial T_{\mathrm{c}}}{\partial t}=\frac{Q^{+}-Q^{-}}{\mathrm{\Sigma }}\frac{2m_{\mathrm{p}}}{5k_{\mathrm{B}}\mu }\hfill \end{array}$$(16)

and

$$\begin{array}{c}\hfill \frac{\partial \psi }{\partial t}=(v_{\mathrm{in}}-v_{\mathrm{out}}+v_{D_{B_z}}{D}_{B_z})\frac{\partial \psi }{\partial R},\end{array}$$(17)

where Ṁ is the mass accretion rate, T_c the temperature in the midplane of the disk, Q⁺ the heating rate, Q⁻ the cooling rate, m_p the mass of a proton, k_B the Boltzman constant, and μ the mean molecular weight. We note that for the sake of simplicity, we neglect mass losses induced by the wind in Eq. (15). We compute the mass accretion rate as

$$\begin{array}{c}\hfill \dot{M}=\frac{4\pi }{R\mathrm{\Omega }}\frac{\partial }{\partial R}\left(\frac{3}{2}{R}^2\alpha (\beta )\mathrm{\Sigma }{\mathrm{\Omega }}^2{H}^2\right)+q(\beta )B_z^2\frac{R}{\mathrm{\Omega }},\end{array}$$(18)

where we take into account the accretion due to MRI turbulence and the outflow magnetic torque on the surface of the disk. The value of α(β) is taken from (11). The value of q is given by Eq. (14). The heating rates and cooling rates, respectively, are given by

$$\begin{array}{c}\hfill Q^{+}=\frac{9}{4}\alpha {\mathrm{\Omega }}^3{H}^2\mathrm{\Sigma }\end{array}$$(19)

and

$$\begin{array}{c}\hfill Q^{-}=2\sigma T_{\mathrm{eff}}^4.\end{array}$$(20)

The effective temperature T_eff is a function of Σ and κ and is computed using the prescriptions from Latter & Papaloizou (2012). We note that the treatment of the thermal structure in our model is simpler than that in Hameury et al. (1998) and Scepi et al. (2019). Indeed, it does not include cooling due to convection, the change in mean molecular weight and thermal capacity as hydrogen is ionized and recombined, or the radial transport of energy and the heating and cooling by the adiabatic expansion of the fluid. In Scepi et al. (2019) we found that the last two terms, in particular, were only important near sharp transitions in surface density or temperature and heating and cooling fronts (see also Menou et al. 1999a). Comparing our model with the one from Scepi et al. (2019), we find that including these terms do not drastically alter the overall dynamics of the disk, although the eruptions tend to be longer, with a slower rise to the eruptive state, and tend to exhibit reflares. We caution that the simplified cooling function of Faulkner et al. (1983), which inspired the prescription of Latter & Papaloizou (2012) that we adopt here, is known to lead to outbursts with a larger amplitude than more realistic treatments (Lin et al. 1985; Hameury et al. 1998). This is especially true when the inner edge of the disk is large (see Figs. 5 and 6 of Lin et al. 1985 and Sect. 3.1.6). Nonetheless, we believe that our model captures most of the necessary physics to understand the coupling between magnetic field transport and the eruptions. We also note that Eq. (19) overestimates the heating rate. Part of the gravitational energy that is usually extracted from the disk by the turbulent stress tensor is actually used to launch the outflow and goes in the vertical Poynting flux instead of being liberated into thermal energy locally (Ferreira & Pelletier 1995).

The novelty compared to Scepi et al. (2019) is that we evolve the magnetic field at the same time as the surface density and the temperature by solving Eq. (17). Equation (17) is coupled to Eqs. (15) and (16) through the prescriptions for v_in, v_out, and v_DB, which depend on β and H/R. We also use a different prescription for the magnetic outflow’s torque, q, taken from the self-similar solutions of Jacquemin-Ide et al. (2019). While prescriptions from shearing box simulations concerning the turbulent torque are believed to be quite accurate, the prescriptions concerning the magnetic outflow’s torque might be better captured by self-similar solutions that take into account the global curvature of the outflow, which is ignored in local shearing box models. In any case, all the results presented here do not depend on the choice of q(β).

We start from a uniform magnetic field configuration. The disk is isolated, meaning that it cannot lose or gain magnetic flux as we enforce v_ψ to cancel at the inner and outer boundaries. We set up Σ to be small at the inner boundary (equal to 1) and the mass accretion rate to be equal to the external mass accretion rate Ṁ_ext = 10¹⁶ g s⁻¹ at the outer boundary. The internal radius of the disk is R_in = 8 × 10⁸ cm, the external radius of the disk is fixed at R_out = 2 × 10¹⁰ cm, and the mass of the white dwarf is 0.6 M_⊙. We use a static logarithmic grid in radius with a resolution 200 points. We went up to a resolution of 800 points and did not see any meaningful difference concerning the overall dynamics of the disk.

3. Magnetic advection dominates for β >  β_eq

In this section, we study two cases:

– turbulent accretion is efficient in advecting magnetic flux inwards at large β (as in Guilet & Ogilvie 2014) and so, v_in >  v_DBz at large β (blue solid line in Fig. 2 and Sect. 3.1).

Fig. 2. Prescriptions for vin, vout, and vDBz normalized by the sound speed as a function of β. The solid blue line is supposed to represent a disk where the magnetic flux’s advection by turbulent accretion is more efficient than diffusion at large β. The dashed blue line gives vin(β) in the case where the magnetic flux’s advection by turbulent accretion is less efficient than diffusion at large β as studied in Sect. 3.2. The dotted black line gives the value of βeq.

– turbulent accretion is not efficient in advecting magnetic flux inwards at large β (as in Lubow et al. 1994a) and so, v_in <  v_DBz at large β (blue dashed line in Fig. 2 and Sect. 3.2).

3.1. Efficient advection at large β

We start with the case where turbulent accretion is efficient in advecting magnetic flux inwards at large β. For this subsection, we use the following prescriptions for v_in and v_out:

$$\begin{array}{cc}& v_{\mathrm{in}}=3\times {\beta }^{-0.2}{c}_{\mathrm{s}}\hfill \\ \hfill & v_{\mathrm{out}}=4f\times {\beta }^{-0.3}{c}_{\mathrm{s}},\hfill \end{array}$$(21)

where f = 1/(1 + 10⁻¹⁰β²) is used to quench the outflow at large β (see Sect. 2.2.4).

Figure 2 shows the prescriptions that we used. In this case, β_eq ≈ 18. Our results are not dependent on the exact value of β_eq, as shown in Appendix B.

3.1.1. Transitory evolution

In this section, we describe the early dynamics of the magnetic field when we simultaneously evolve the density, temperature, and magnetic flux using the prescriptions from Fig. 2. We start from an initially uniform vertical large-scale magnetic field of 30 G (top panel of Fig. 3) and an external mass transfer rate of 10¹⁶ g cm⁻². The lower panel of Fig. 3 shows the initial distribution of β. Initially, β decreases from 10⁹ at the inner radius to 10⁴ at the outer radius. Since β >  β_eq, the magnetic flux is advected inwards (see Fig. 2).

Fig. 3. Temporal evolution of the magnetic field configuration (upper panel) and plasma β parameter (lower panel) up to the quasi-stationary final configuration. The initial magnetic field is uniform with a value of 30 G. The external mass transfer rate is 1016 g cm−2. The dashed black line in the upper panel gives the theoretical expected slope of the magnetic field, in R−5/4 deduced from Eq. (23). The dotted black line in the lower panel gives the value of βeq.

The magnetic flux is advected inwards until the magnetic configuration reaches an equilibrium at β = β_eq, where diffusion outwards (v_out) compensates advection inwards (v_in). We see from Fig. 3 that after 1.4 days the inner regions have already reached the equilibrium magnetization while the outer disk still advects the magnetic field inwards. As the magnetic flux is accumulated inwards, the inner magnetically dominated region grows with time. At t ≈ 30 days the disk has almost reached a quasi-stationary state. Between, t ≈ 30 days and t ≈ 1200 days the disk is in a stationary state except for small fluctuations in all the quantities. The mass accretion rate is constant and the disk does not undergo thermal eruptions. This would correspond to the case of a nova-like structure. The inner region of the disk concentrates most of the magnetic flux and is stabilized at a β slightly above β_eq (see lower panel of Fig. 3). This small difference is due to the gradient of B_z transporting magnetic flux outwards. The outer region of the disk is emptied of its magnetic flux and reaches the floor in β. We prescribed a ceiling at β = 10¹⁰ for numerical convenience but this does not affect our results.

3.1.2. Timescales

We see from Fig. 3 that advection of the magnetic flux inwards is fast. The timescale for the evolution of the magnetic field, t_mag, is given by Eq. (8) and is

$$\begin{array}{c}\hfill t_{\mathrm{mag}}=\frac{R}{v_{\psi }}·\end{array}$$

In a highly magnetized disk where β ≈ 1 and the accretion is due to a surface torque from a magnetized outflow, v_in and v_out are expected to be on the order of the sound speed c_s (see Sects. 2.2.2 and 2.2.3) so that

$$\begin{array}{c}\hfill t_{\mathrm{mag}}\approx \left(\frac{R}{H}\right)\frac{1}{\mathrm{\Omega }}\end{array}$$

when β ≈ 1. For example, at R = 10⁹ cm, t_mag ≈ 1 h at the beginning of our simulation, which is consistent with the timescales on Fig. 3.

Moreover, the accretion velocity in this case is also ≈c_s, so that the accretion timescale, t_acc is also $\left(\frac{R}{H}\right)\frac{1}{\mathrm{\Omega }}$. This is much smaller than the typical viscous timescale of an α-disk, $t_{\mathrm{vis}}=\frac{1}{\alpha \mathrm{\Omega }}{\left(\frac{R}{H}\right)}^2$, which is ≈10 years at the outer radius at the beginning of our simulation.

When β ≈ 1, we have α ≈ 1 and so the magnetic timescale and the accretion timescale are only an order R/H larger than the thermal timescale, t_th ≡ 1/(αΩ). In these conditions, the magnetic configuration and the density can adapt very rapidly to the thermal structure. The classical thermal-viscous instability may exhibit a different behavior in these conditions.

3.1.3. Magnetic configuration

The magnetic field, after a transition regime, maintains the same configuration from t ≈ 1 day up until t = 1200 days. The inner magnetic field follows a power law with an index of R^−5/4, and the outer magnetic field is almost null. We define R_tr as the transition radius between the inner magnetized region and the outer non-magnetized region; it also corresponds to the radius where the magnetic field deviates from the power law in −5/4. When R_tr increases, we see that the value of the magnetic field at the inner boundary, B_zin, decreases.

We can explain those results using a very simple model. The inner outflow-driven accretion zone has β ≈ β_eq, which is constant. In the inner zone, the accretion is driven by the magnetized outflow’s torque and the mass accretion rate,

$$\begin{array}{c}\hfill \dot{M}=q{B}_z^2\frac{R}{\mathrm{\Omega }},\end{array}$$(22)

is constant since the disk is steady. Since β in the inner zone remains fixed around β_eq, q(β) remains at q(β_eq). Isolating B_z in Eq. (22), we find

$$\begin{array}{c}\hfill B_z\propto R^{-5/4}\end{array}$$(23)

in the inner zone, as in Blandford & Payne (1982). This simple analytical estimate is perfectly reproduced by our model as we can see from the black dashed line on the top panel of Fig. 3.

In an isolated disk (v_ψ(R_out) = v_ψ(R_in) = 0), the magnetic flux is conserved. If we assume that all the magnetic flux is concentrated in the inner part, with B_z ∝ R^−5/4 then the conservation of the magnetic flux,

$$\begin{array}{c}\hfill \mathrm{\Psi }={\int }_{R_{\mathrm{in}}}^{R_{\mathrm{tr}}}{B}_{z}R\mathrm{d}R,\end{array}$$

where Ψ is the total magnetic flux in the disk, can be rewritten as

$$\begin{array}{c}\hfill R_{\mathrm{tr}}=R_{\mathrm{in}}\times {(\frac{3}{4}\frac{\mathrm{\Psi }}{B_{z_{\mathrm{in}}}{R}_{\mathrm{in}}^2}+1)}^{4/3}.\end{array}$$(24)

We can easily constrain the magnetic field at the inner boundary with Eq. (22) and we find

$$\begin{array}{c}\hfill B_{z_{\mathrm{in}}}={\left(\frac{{\dot{M}}_{\mathrm{in}}{\mathrm{\Omega }}_{\mathrm{in}}}{q{R}_{\mathrm{in}}}\right)}^{1/2}.\end{array}$$(25)

The values of R_tr and B_zin are in perfect agreement with our simulations, as we can see from the upper panel of Fig. 4 where we report the temporal and radial evolution of the magnetic field with a dashed black line representing the analytical estimation (24) in which the value of B_zin is set by Eq. (25).

Fig. 4. Magnetic field (top panel), surface density (middle panel), and temperature (bottom panel) as a function of time and radius. The black dashed line in the top panel gives the analytical estimate of the transition radius using Eqs. (24) and (25). The initial magnetic field is uniform with Bz = 10 G. The external mass transfer rate is 1016 g cm−2.

Finally, we point out that in a disk where the magnetic flux is conserved and approximating that R_tr ≫ R_in, the equations above imply that the transition radius evolves as

$$\begin{array}{c}\hfill R_{\mathrm{tr}}\propto {\dot{M}}_{\mathrm{in}}^{-2/3}.\end{array}$$(26)

This behavior is somewhat different from that of the magnetospheric radius, R_mag, which is set by equating the ram pressure of the accretion disk with the magnetic pressure from a possible dipolar magnetic field anchored in the white dwarf. Accretion then proceeds along magnetic field lines, truncating the disk within R_mag. The magnetospheric radius behaves as $R_{\mathrm{mag}}\propto {\dot{M}}_{\mathrm{in}}^{-2/7}$ (Frank et al. 2002). We stress that in our case the disk is not truncated below R_tr but its properties are changed as it becomes highly magnetized.

3.1.4. Temporal evolution

Figure 4 shows the evolution as a function of time and radius of the magnetic field, the surface density, and the temperature. These results were found using an initial constant magnetic field of 3 G and a mass transfer rate of 10¹⁶ g cm⁻¹.

Figure 4 shows that the disk is unstable to the thermal-viscous instability. The only difference with the stable disk shown in Fig. 3 is a lower magnetic field (B_z = 10 G) than in the stable case (30 G). In the high magnetic flux case of Fig. 3 the transition radius R_tr is larger than the radius where the thermal-viscous instability is triggered. This stabilizes the disk and allows it to find a stationary state with a constant mass accretion rate. We recover the instability by decreasing the magnetic flux in the disk.

We see from Fig. 4 that the eruptions are triggered at the transition between the inner magnetized and the outer non magnetized disk. Starting from the quiescent state, the disk accumulates mass until Σ reaches a critical value where the disk becomes thermally unstable and transits to a hotter state, corresponding to the eruptive state. We see two fronts that propagates through the disk: one inwards and the other outwards. The disk stays in eruption for approximately 40 days. Then, a cooling front coming from the outer boundary of the disk signals the end of the outburst as it propagates inwards. We note that the disk has reflares (Menou et al. 1999a) when going back to quiescence. These show up as oscillations in the size of the hot zone as the disk returns to quiescence (bottom panel of Fig. 4). This means that when the cooling front propagates, the density reaches the critical density again and a new heating front appears. Then, another cooling front develops at the outer edge and, if the density in the disk is low enough, the cooling front can propagate all the way down to the transition radius ending the eruption. The disk then returns to quiescence and stays in quiescence for approximately 80 days. The recurrence of these eruptions depends on the initial magnetic field and the parameters of the disk, as we will see in Sect. 3.1.6.

From Fig. 5, we also see that the magnetic configuration changes between the eruptive state and the quiescent state. In fact, the evolution of the transition radius is driven by the thermal-viscous instability operating in the outer unmagnetized disk. We see that R_tr increases when going from eruption to quiescence in response to variations in Ṁ_in as in Eq. (26). Going back to quiescence, β at the transition decreases since the disk temperature drops. When β decreases below β_eq the disk enters a zone where the diffusion is dominant, and the inner magnetic flux diffuses outwards until a β = β_eq is reached again. On the contrary, when going from quiescence to eruption, the value of β increases. At the transition, the magnetic field is then advected inwards since β becomes greater than β_eq. The transition radius moves inwards as the heating front propagates. The evolution of the inner disk is entirely dictated by the eruptions in the outer disk.

Fig. 5. Surface density (first panel), temperature (second panel), total mass accretion rate (third panel), and vertical magnetic field (fourth panel) as a function of radius in eruption (green curves) and quiescence (blue curves). The dashed curves in the second panel represent the effective temperature of the disk. The dashed curves in the third panel represent the outflow-driven mass accretion rate and the dotted curves represent the turbulent-driven mass accretion rate. The initial magnetic field is uniform with Bz = 10 G.

3.1.5. Radial structure

Figure 5 shows the disk’s radial structure in quiescence and in eruption. The disk is radically different in the outer and inner parts during both eruption and quiescence.

The outer disk is emptied of magnetic flux. The value of β in this zone is set by the floor of β = 10¹⁰ that we impose artificially. Thus, the outer parts behave as a standard α-disk with a constant α ≈ 0.03. The accretion is due to turbulence and the accretion energy is deposited locally in the disk. The accretion driven by the torque of the magnetized outflow is very weak in this zone. It is in the outer disk that the thermal-viscous instability takes place. As we see from Fig. 5, the outer disk in outburst is hot, dense, and optically thick. In quiescence, it is cold and has a lower density.

The inner disk, where most of the magnetic flux has accumulated, behaves very differently from an α-disk. It reaches the value β_eq ≈ 18 where accretion is driven by the torque from the magnetized outflow. Because this magnetized accretion is so efficient, that is, v_R ≈ q/βc_s ≥ c_s, the density in the inner disk is low. A significant fraction of the accretion energy is not deposited locally but is lost in the outflow. Thus, the inner disk also has a lower temperature than a viscously-driven disk for the same Ṁ because of the lower heating rate. However, we find that the densities and temperatures in the inner disk are not as low as the values reached in Scepi et al. (2019), where the fixed magnetic field configuration led to very low values of β in the inner disk (Figs. 2 and 5 in Scepi et al. 2019). Surprisingly, both Σ and T_c remain high enough that the thermal-viscous cycle continues to operate in the inner disk. Figure 5 also show the presence of a small hot zone as T_c in the magnetized inner disk, whereas the viscously-driven outer disk is in quiescence. This produces small eruptions on a very short timescale because of the fast accretion timescale in the magnetized inner disk. These eruptions shows up as wiggles during quiescence in panels c and d of Fig. 6.

Fig. 6. Mass-accretion rate at the inner boundary (left panels) and V magnitude (right panels) as a function of days for different magnetic flux. From top to bottom: the initial constant magnetic field Bz0 is equal to 0.1, 1, 10, and 21 G.

3.1.6. Light curves

We showed in Sect. 3.1.4 that the magnetic field configuration is self-consistently changing from the eruptive state to the quiescent state. As we can see from Fig. 6, by changing the value of the initial magnetic flux, we can change the shape of the eruptions.

In a weak initial magnetic field (B_z0 = 0.1 G), accretion driven by the magnetized outflow is negligible as there is not enough magnetic flux in the disk. The disk is mainly turbulent, with α ≈ 0.03, and the eruptions are of small amplitude with short recurrent timescales (panel a of Fig. 6). We attribute the differences in amplitude and the recurrence timescale with the low magnetic field case shown in the top panel of Figs. 3 and 4 in Scepi et al. (2019), where the same disk is modeled to our simplifying assumptions (Sect. 2.3). As we increase the magnetic flux in the disk, the accretion driven by the magnetized outflow starts to be important in the inner disk. The transition radius acts as an effective truncation since the accretion timescale becomes very short below R_tr. Indeed, Fig. 6 shows the light curve behaves as expected from a disk with a progressively higher truncation radius as the magnetic flux is increased. The heating front that triggers the outburst starts close to R_tr for these parameters. Since the critical surface density for the ignition of the viscous-thermal instability is proportional to R^1.1 (Hameury et al. 1998), a truncated disk needs to build more mass between outbursts, increasing the recurrence timescale from panels a–c. The mass accretion rate in quiescence also increases as R_tr increases (left panels of Fig. 6), an effect that has long been identified as a possible solution to the high X-ray flux measured from DNe in quiescence (Mukai 2017, and references therein). The disk becomes “leaky” when its Ṁ_in in quiescence becomes close to the mass transfer rate (Lasota 2001), preventing mass from building up in the disk and further increasing the recurrence timescale as illustrated when going from panels c and d. These results do not depend strongly on the exact value of β_eq. We retrieved similar light curves with different β_eq by changing only the amount of magnetic flux in the disk (see Appendix B).

3.2. Inefficient advection at large β

We now briefly explore the second case where the advection of magnetic flux at large β is inefficient. This case aims at mimicking the results of Lubow et al. (1994a). We use the following expressions for v_in:

$$\begin{array}{c}\hfill v_{\mathrm{in}}=3f\times {\beta }^{-0.2},\end{array}$$(27)

where f = 1/(1 + 10⁻¹⁰β). v_out is still defined according to Eq. (21). The only difference with our previous prescription is that v_in is also quenched at large β as shown by the dashed blue line in Fig. 2. Hence, there is a magnetization above which the magnetic flux cannot be advected anymore because v_in <  v_DBz (see Fig. 2). In Lubow et al. (1994a), the dominant term in the diffusion of the magnetic flux is due to the vertical diffusion of the radial magnetic field, where the radial magnetic field is computed from the global magnetic structure. Our model is local and is unable to capture such an effect. We use a case where v_in <  v_DBz as a proxy for the global method of Lubow et al. (1994a).

We do not discuss this case in as much detail as the one that sees efficient advection at large β. Indeed, the results are very similar except that the light curves obtained with this prescriptions show eruptions of smaller amplitude (≈two magnitudes) as can be seen in Fig. 7. This is due to the fact that the magnetic field at the inner boundary cannot change as much between the eruptive state and the quiescent state as in the case of efficient advection (see a comparison between both cases on Fig. 7). In the case of inefficient advection, the disk can only advect inwards (or diffuse outwards) the magnetic flux from the parts of the disk where β ≲ 10⁵. This leaves the outer parts of the disk out of reach and limits the range over which the magnetic field at the inner boundary can vary. Ultimately, this limits the range over which the transition radius can change (see Fig. 7) and produces smaller eruptions than in Sect. 3.1.6.

Fig. 7. Upper figure: temporal and radial evolution of the magnetic field (upper panels) and of the outflow-driven mass accretion rate (lower panels) for the case of inefficient advection of the magnetic flux at high β (left hand side panels) and for the case of efficient advection of the magnetic flux at high β (right hand side panels). Lower figure: light curves for inefficient advection at large β. The initial constant magnetic field in this case is Bz0 = 20 G.

4. Discussion

4.1. Limitations of our model

Our model assumes that the vertical velocity and the radial magnetic field are small enough so that we can neglect the term v_zB_R compared to v_RB_z in the induction equation. However, the velocity, v_z, can be important, especially in the inner magnetized region where a wind or a jet is formed. In the case where β_eq is 1−100, v_z becomes ≳2 times the local orbital speed (Jacquemin-Ide et al. 2019) and so v_ZB_R/v_RB_Z can no longer be neglected since B_R/B_Z ≳ 1 and v_R ≲ c_s. This could induce another source of magnetic flux transport.

Moreover, we assume that there exists an outflow transporting angular momentum but we do not treat the transport of mass and energy in the outflow. This means that the disk may be even emptier and cooler than it is in the magnetized outflow-dominated inner parts on Fig. 4.

Also, our approach is different from that of Lubow et al. (1994a), Guilet & Ogilvie (2014) or Li & Cao (2019) in that they use a global model, whereas our approach is purely local. We assume that the behavior of the disk only depends on the local magnetization and not on the full global structure. We note that the prescriptions from Jacquemin-Ide et al. (2019) are self-similar and so, by nature, they take into account the global structure of the disk. However, due to the self-similar approximation, the solutions from Jacquemin-Ide et al. (2019) have a uniform magnetization (not dissimilar to global simulations from Zhu & Stone 2018; Mishra et al. 2020) and are not consistent with a realistic disk where the magnetization depends on the radius and is not self-similar. This is especially true close to the fronts.

We would like to emphasize that the global approach followed by Lubow et al. (1994a), Guilet & Ogilvie (2014) and Li & Cao (2019) relies on the approximation that the magnetic field is potential. When an outflow is present, this is clearly not accurate because the outflow matter’s inertia exerts a force on the magnetic field configuration. Our local approximation may be more suited than the potential approximation when an outflow is present, although it still needs a proper support from global simulations. Nevertheless, the diffusion terms in the outer disk, which has a high β and no outflow, may be underestimated as we do not take into account the global structure.

At high magnetizations (β ≲ 10²), other instabilities that we ignored here, such as the magnetic Rayleigh-Taylor instability, might play a role in giving the final magnetic configuration as seen in simulations of McKinney et al. (2012). Taking into account these high magnetization instabilities would require a detailed study providing a simple prescription of the mass, angular momentum (see Marshall et al. 2018), energy, and magnetic field transport that we can use in our model. However, it is not clear that a simple axisymetric model would be able to reproduce the long-term behavior of high-magnetization simulations, which are reported to be highly non-axisymetric with spiral patterns that are possibly related to the properties of transport (McKinney et al. 2012; Mishra et al. 2020).

Beckwith et al. (2009) also expressed concerns regarding the use of a mean field approach to treat the turbulence. In their 3D GRMHD global simulations, turbulent motions reconnect field lines in the surface accretion layer forming closed field loops in the disk and allowing coronal advection of open magnetic field lines on the horizon of the central black hole. A detailed study of whether this effect can be captured using a mean field approach remains to be done.

Finally, we restricted our model to a case where the magnetic flux is conserved. It is not clear whether such a case is relevant for compact binaries where the disk is fed with matter and possibly magnetic flux at the outer radii. We tried to address this problem by continuously feeding the disk with a net vertical magnetic field. However, in the absence of any mechanism to extract magnetic flux from the disk, the magnetic flux keeps accumulating until we stop the simulation. To prevent a dramatic accumulation of magnetic flux, we also tried to feed the disk with a vertical field of switching polarity. In the case of efficient advection at large β, we mostly find modulations of the light curves with an amplitude and a frequency corresponding to those of the injected magnetic flux. It seems unlikely that the eruptions of DNe should actually be driven by the varying incoming magnetic field and so, we did not pursue in this direction. Internal local generation of net vertical magnetic field by the MRI dynamo might also play a role in the long term evolution of DNe (see Begelman & Armitage 2014 for low mass X-ray binaries). However, exploring such effects is clearly beyond the scope of this paper.

4.2. Comparison with Scepi et al. (2019)

One major result of this paper is that the inner disk, where accretion is driven by the torque due to the magnetized outflow, acts exactly as a truncation of the inner disk. The inner disk is passive and responds to the changes in the outer disk; it does not drive the dynamics of the eruptions. Our hybrid disk can then simply be modeled as a truncated α-disk with a varying inner radius, R_in, which follows the relations of Eqs. (24) and (25). To check this assertion, we performed a simulation of such a truncated α-disk and found the same behavior than reported in Fig. 4.

This is different from the results of Scepi et al. (2019), where the inner zone acted only partly as a truncation of the disk. Because the magnetic field configuration had a fixed dipolar configuration, the transition in outflow-driven mass accretion rate was more shallow than in the present work. Hence, the torque of the magnetized outflow also enhanced the accretion in the outer disk (especially during outburst) and did not act as a simple truncation. Indeed, we checked in Scepi et al. (2019) that a simple truncation of the outer α-disk led to weak outburst with short recurrence timescales. Here, the magnetic flux is much more concentrated to the disk center and the transition between magnetized and unmagnetized regions is sharper.

Another difference is that the relation between the transition radius, R_tr, and the mass accretion rate reported here is different from the one of Scepi et al. (2019). In the case of a fixed dipolar magnetic field, R_tr evolves as Ṁ^−2/7 (like the magnetospheric radius) while here it evolves as Ṁ^−2/3 as in Eq. (26). Since a disk can be stabilized if R_tr >  R_inst, where R_inst is the radius at which the thermal instability would be triggered, this modifies the dependency with the orbital period of the stability criterion shown on the Fig. 6 of Scepi et al. (2019) for a given amount of magnetic flux.

4.3. Comparison with global simulations

Because of the difference of timescales, it is difficult to directly compare our work with 3D global (GR-)MHD simulations. Here, we are interested in the secular evolution, that is, taking place over several viscous timescales, whereas state-of-the-art 3D global simulations could typically only reach hundreds of orbital timescales at the outer radius of DNe. Still, Liska et al. (2018) performed a long-duration GRMHD 3D global simulation and found that at late times, the accretion disk is divided into an inner magnetized, low-density region and an outer region sharing the properties of an α-disk. The authors suggested that a large-scale magnetic torque might be responsible for the formation of this inner region.

More generally, GRMHD simulations studying magnetic flux transport report that the magnetic flux is advected very efficiently inwards (Beckwith et al. 2009; McKinney et al. 2012). These simulations typically start with a torus and allow the accretion disk to form. However, the most recent MHD simulations, which initially start with an accretion disk with a constant β, tend to show that the magnetic flux is neither advecting inwards or outwards (Zhu & Stone 2018; Mishra et al. 2020), which is at odds with the GRMHD results. This result, surprisingly, seems to hold whatever value of β between 10² and 10⁴. This conflict between GRMHD and MHD simulations clearly deserves further attention from the community.

4.4. Comparison with observations

Allowing the magnetic field distribution evolve under the assumption of efficient advection concentrates the flux in an inner magnetized disk. The outer disk is purely viscous and behaves dynamically as if it were truncated at the transition radius between the two zones. This has several consequences for observations.

First, this provides an alternative mechanism through which disk truncation can be achieved. Disk truncation has been invoked to explain the delayed rise to outburst of the UV flux, compared to the optical flux, and the relatively high accretion rate onto the white dwarf inferred from the X-ray observations in quiescence, compared to the predictions from the DIM for a disk extending down to the surface of the white dwarf (see Lasota 2001; Mukai 2017, and references therein). Various mechanisms have been proposed to explain this: truncation by the white dwarf magnetic field as in intermediate polars (Livio & Pringle 1992), evaporation of the inner disk by a “coronal siphon flow” (Meyer & Meyer-Hofmeister 1994), or the white dwarf irradiation keeping the inner disk hot (King 1997). We propose that it is the magnetic field carried by the disk itself that causes an apparent truncation of the disk. Such a mechanism would apply equally well to low-mass X-ray binaries, as we discuss further below.

Second, the truncation of the inner disk by the outflow-driven accretion eventually shapes the light curve. The eruptions have a longer recurrence timescale and a larger amplitude than in the case of a purely turbulent disk with a constant α (Fig. 6). Coleman et al. (2016) argued that truncation of the inner disk was essential to match outburst models, using the actual values of α measured in simulations of MRI turbulent transport to observations. They proposed truncation by the white dwarf dipolar field. Our model could produce such a truncation without requiring every DNe to have strong magnetic fields on their white dwarf. We also predict a different evolution of the truncation radius than for the magnetospheric radius.

Third, the magnetic flux in the disk will be one of the two main parameters, with the external mass transfer rate, shaping the light curve. Observationally, there is a large scatter in the recurrence timescale of DNe (Cannizzo et al. 1988) that is usually explained by external mass transfer rate variations at a given orbital period. In our model, we could naturally account for the dispersion in recurrence timescales by varying the initial magnetic flux in the disk. We defer a more detailed comparison of light curves to observations using Eq. (26) to evolve the truncation radius to a future work, together with a more realistic treatment of the radiative transfer of the thermodynamics of the disk (Hameury et al. 1998).

Fourth, a high enough magnetic flux will also affect stable systems. It can stabilize a disk that would be unstable to the classical thermal-viscous instability by truncating the inner hot region, leaving only a stable outer cold region, much like posited for magnetospheric truncation (Lasota et al. 1995; Menou et al. 1999b). This would produce a population of cold, stable systems. Hot, stable systems (novæ likes) could also be impacted. Actually, in novae-like variables, a discrepancy exists between the observed spectrum and the spectrum expected from a classical α-disk extending to the inner radius. On the one hand, the emission coming from the inner radii of the disk seems to be missing (see Nixon & Pringle 2019 for recent work and review on the subject). On the other hand, dwarf novæ in eruption, whose observational properties are expected to be close to those of novæ likes, are well-fitted by an α-disk extending to the inner radius (Hamilton et al. 2007). The discrepancy could simply arise because the disk for the same mass accretion rate can be more or less truncated depending on the available magnetic flux.

Lastly, one observational signature in our model is the presence of magnetized outflows, jets, or winds. Recent observations of radio emission in DNe and novæ-like have been attributed to magnetic jets (Körding et al. 2008; Russell et al. 2016) although more than one mechanism could produce such emission (Coppejans et al. 2015, 2016). If magnetic jets are indeed responsible for the radio emission, we could expect from our model a correlation between strong radio emitters and DNe with a long-recurrence timescale. Also, since DNe with a large magnetic flux are stabilized, we could also expect to find cold, stable systems with strong radio emission. Matching stable, point-like sources in radio surveys to late-type stars in optical surveys may be a way to identify those optically faint systems.

4.5. Application to LMXBs

Low mass X-ray binaries are analogues of DNe where, instead of a white dwarf, the central object is a stellar mass black hole or a neutron star. The spectral and variability properties of LMXBs provide several lines of evidence for a truncated thin disk (see Yuan & Narayan 2014 for a review and Zdziarski & De Marco 2020 for more recent work). These objects have eruptions that can also be explained by the DIM taking into account the X-ray irradiation of the outer region by the inner region and a truncated disk (Lasota et al. 1996; Dubus et al. 2001).

Ferreira et al. (2006) suggested (and Marcel et al. 2018a,b later showed) that all the spectral states of LMXBs could be reproduced using an hybrid disk composed of: (1) an inner jet-emitting disk (JED), with β ≈ 1, where all the angular momentum and part of the energy are transported vertically by a magnetic jet; (2) a standard α accretion disk (SAD). By varying only the inner mass accretion rate and the transition radius between the JED and the SAD, Marcel et al. (2019) were able to reproduce the outburst of GX 339−4.

Our disk configuration is not very different from what has been assumed by Ferreira et al. (2006) and Marcel et al. (2019). The main difference between the inner JED disk and our inner magnetized outflow-driven accretion dominated zone is the degree of magnetization. While the JED solutions require β ≈ 1−10 we have a β that can range from 10 to 10³ in our model without affecting the light curves (see Appendix B). The SAD, outer disk of Marcel et al. (2019) is our outer, non-magnetized, turbulent-driven accretion disk. A hybrid disk with a sharp transition between these two components being a natural outcome of our model, we propose that our model might provide a dynamical justification for the model of Ferreira et al. (2006) and also might be able to explain part of the spectral evolution of LMXBs. This will be the subject of a future work.

5. Conclusions

We studied the evolution of a non-steady, axisymmetric, DNe accretion disk undergoing eruptions due to the thermal-viscous instability. The novelty of this work is that we evolve the large-scale poloidal magnetic configuration during the outburst cycle using a local toy model for magnetic flux transport. The disk is isolated in the sense that it cannot lose or gain magnetic flux. To evolve the density, temperature, and magnetic field, we developed a version of the classical DIM (Hameury et al. 1998) that takes into account the removal of angular momentum by MRI turbulence and by magneto-centrifugal outflows as in Scepi et al. (2019). We assumed that the magnetic flux evolution only depends on the local magnetization and used several different sets of prescriptions for magnetic flux transport to study their impact on the outburst light curves.

The only case where we find realistic light curves is when magnetic flux is advected inwards until it reaches an equilibrium magnetization, β_eq, for which inward magnetic field advection is balanced by outward magnetic field diffusion. This case naturally leads to a hybrid disk configuration, regardless of the exact value of β_eq. The inner region of the disk is highly magnetized and has a uniform β(R) = β_eq. As a consequence, the magnetic field is distributed as R^−5/4 in this region. The angular momentum and accretion energy are mostly transported vertically in the outflow. The outer region of the disk is emptied of magnetic field thus outflow-driven accretion is negligible there. The outer disk resembles a classical, viscously-driven, α-disk.

The transition between the inner outflow-dominated region and the outer turbulent region is very sharp and is of the order of the typical scale height of the disk. We find that our hybrid disk can actually be reduced to a truncated α-disk with an inner boundary evolving as ${\dot{M}}_{\mathrm{in}}^{-2/3}$, which can be easily implemented in DIM codes. Thus, our model provides a physical justification to previous models that require (or observations that suggest) a truncation of the viscous disk in both DNe and LMXBs (see Hameury 2020, and references therein). Here, the inner regions are replaced by a disk with a very fast accretion timescale driven by a strongly magnetized outflow.

Acknowledgments

We thank Jean-Pierre Lasota, the referee, for a thorough and critical report that led to a better understanding and presentation of our results. NS would like to thank Jerôme Guilet for help at the beginning of this project, Grégoire Marcel for many discussions and sharing unpublished results as well as Mitch Begelman and Jason Dexter for reading an early version of this work. NS acknowledges partial financial support from the pôle PAGE of the Université Grenoble Alpes and from a NASA Astrophysics Theory Program grant NNX16AI40G. 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. GL acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 815559 (MHDiscs)). GD acknowledges support from CNES.

References

Appendix A: v_ψ from Guilet & Ogilvie (2012)

In a first approach, we attempted to compute v_ψ as a function of β by obtaining the vertical structure of v_R from the asymptotic expansion of Guilet & Ogilvie (2012) instead of simple power laws. As mentioned in Guilet & Ogilvie (2012), there is a tight connection between solving the vertical structure of the disk in the thin disk asymptotic expansion and solving the eigenvalue problem for viscous, resistive MRI. Indeed, the vertical structure is a linear combination of MRI modes. The growth rates and the eigenvectors of the resistive MRI modes depends on the resistivity and viscosity and so, ultimately, on α. Guilet & Ogilvie (2012) used the value of α necessary to damp all MRI unstable modes to find v_ψ. However, this leads to unrealistic values, for example α ≈ 0.7 for β ≈ 10⁵. In our case, we use the value of α from 3D local MRI simulations with radiative transfer from Scepi et al. (2018a). In doing so, we do not ensure that all MRI modes are stable for the whole range of β we considered. In fact, it is natural that unstable MRI modes with our choice of α exist since the simulations we are based on are fully turbulent. When we use 𝒫 = 1, we find that near the values of β where resistive MRI modes are damped by resistivity, our problem is ill-conditioned and we cannot find one unique solution. This is fundamentally due to the fact that it is physically difficult to find a stationary linear solution to a system that still has unstable modes.

We overcame that problem by enhancing the resistivity to damp all the unstable MRI modes. We found that for 𝒫 ≲ 0.01, v_ψ is well behaved as a function of β. However, depending on the exact value of 𝒫, we can be in different scenarii of Fig. 1. Then, given the uncertainty of our method, we switched to a more simplified method using simple power laws, which allow us to probe a larger range of behavior.

Appendix B: Light curve dependance on β_eq

The aim of this paper is to show which parameters are important to have realistic dwarf novæ light curves. In this appendix, we change the value of β_eq by slightly changing the slope of v_out. We use three different values of β_eq: 18, 320, 1320. We find that the exact choice of β_eq do not affect the shape of the light curves as can be seen from Fig. B.1. Indeed, by adjusting the initial magnetic flux we can find in each case light curves that look alike one another. Of course as β_eq becomes large the magnetized outflow’s torque becomes less important changing the light curves. However, as long as β_eq is ≲1000−10 000 the exact value is not crucial to the overall shape of the light curves.

Fig. B.1. Dwarf novæ light curves as a function of time for different values of βeq. We used Bz0 = 10, 4 and 2.5 G for the panels a–c, respectively.

Appendix C: An outer magnetized disk

In the fourth case of Fig. 1, v_out >  v_in for β >  β_eq. Consequently, we could imagine the opposite magnetic field configuration compared to Sect. 3.1: an outer magnetized disk and an inner weakly magnetized disk. As in Sect. 3.1.3, the outer magnetized disk would have a constant magnetization β_eq and thus a magnetic field distributed as B_z ∝ R^−5/4. Again, similarly to Sect. 3.1.3, we have the following constraint on the outer magnetic field:

$$\begin{array}{c}\hfill B_{z_{\mathrm{out}}}={\left(\frac{{\dot{M}}_{\mathrm{out}}{\mathrm{\Omega }}_{\mathrm{out}}}{q{R}_{\mathrm{out}}}\right)}^{1/2}.\end{array}$$(C.1)

If B_zout is too large then Ṁ(R_out) becomes larger than the mass transfer rate from the companion, Ṁ_ext. In this case, the outer disk empties leading to a decrease in β. The outer disk then falls in the regime β <  β_eq, causing the magnetic flux to be advected inwards to return to the equilibrium magnetization β_eq. However, the outer region acts as an enormous reservoir of magnetic flux. A small magnetic field advected inwards becomes a large magnetic field in the inner region. This also means that as the magnetic field is advected inwards, the outer magnetic field (and so the mass accretion rate) does not vary much. Thus, the disk keeps emptying and ultimately becomes strongly magnetized everywhere. In this case, there is no eruption.

To avoid this situation, we need Ṁ(R_out) to be on the order of Ṁ_ext. For Ṁ_ext ≈ 10¹⁶ g s⁻¹, we can find using Eq. (C.1) that B_zout ≲ 3 G, where we have used q ≈ 1 to set an upper limit for the outer magnetic field. Setting Σ ≈ 10² g cm⁻¹, H ≈ 10⁸ cm, we must have β_eq ≳ 10⁴ to satisfy Ṁ(R_out) ≈ Ṁ_ext. We tried a simulation with β_eq = 10⁵ and Ṁ_ext = 10¹⁶ g s⁻¹ and found small amplitude eruptions (≈2 mag) with very short recurrence timescales (≈1 day).

We could also conceive of an increase of Ṁ_ext to allow lower values of β_eq. However, we are limited by the fact that the disk needs to be thermally unstable. If Ṁ_ext is too large, the disk stays in a stationary hot state, that is, a novæ-like case. We tried a case where Ṁ_ext = 10¹⁸ g s⁻¹ and β_eq = 10³ and found no eruptions.

There may be some other configurations where we could find other behaviors that are yet to be explored. However, since we are not sure of the astrophysical relevance of such a solution, we did not investigate this matter further. Indeed, in protoplanetary disks, the jet (and so the highly magnetized region of the disk) comes from the inner region of the disk and not the outer region (Lee et al. 2017). In DNe, eclipse mapping is also able to reveal an outer disk, which is quite in agreement with a typical viscous α-disk (Horne & Cook 1985).

All Figures

	Fig. 1. Sketch of the four set of prescriptions obeying the constraints on vin and vout from Sect. 2.2 with a comment on the resulting light curves and disk configuration. Green lines represent the outward transport and blue lines the inward transport.
In the text

Fig. 2. Prescriptions for vin, vout, and vDBz normalized by the sound speed as a function of β. The solid blue line is supposed to represent a disk where the magnetic flux’s advection by turbulent accretion is more efficient than diffusion at large β. The dashed blue line gives vin(β) in the case where the magnetic flux’s advection by turbulent accretion is less efficient than diffusion at large β as studied in Sect. 3.2. The dotted black line gives the value of βeq.

In the text

Fig. 3. Temporal evolution of the magnetic field configuration (upper panel) and plasma β parameter (lower panel) up to the quasi-stationary final configuration. The initial magnetic field is uniform with a value of 30 G. The external mass transfer rate is 1016 g cm−2. The dashed black line in the upper panel gives the theoretical expected slope of the magnetic field, in R−5/4 deduced from Eq. (23). The dotted black line in the lower panel gives the value of βeq.

In the text

	Fig. 4. Magnetic field (top panel), surface density (middle panel), and temperature (bottom panel) as a function of time and radius. The black dashed line in the top panel gives the analytical estimate of the transition radius using Eqs. (24) and (25). The initial magnetic field is uniform with Bz = 10 G. The external mass transfer rate is 1016 g cm−2.
In the text

Fig. 5. Surface density (first panel), temperature (second panel), total mass accretion rate (third panel), and vertical magnetic field (fourth panel) as a function of radius in eruption (green curves) and quiescence (blue curves). The dashed curves in the second panel represent the effective temperature of the disk. The dashed curves in the third panel represent the outflow-driven mass accretion rate and the dotted curves represent the turbulent-driven mass accretion rate. The initial magnetic field is uniform with Bz = 10 G.

In the text

	Fig. 6. Mass-accretion rate at the inner boundary (left panels) and V magnitude (right panels) as a function of days for different magnetic flux. From top to bottom: the initial constant magnetic field Bz0 is equal to 0.1, 1, 10, and 21 G.
In the text

Fig. 7. Upper figure: temporal and radial evolution of the magnetic field (upper panels) and of the outflow-driven mass accretion rate (lower panels) for the case of inefficient advection of the magnetic flux at high β (left hand side panels) and for the case of efficient advection of the magnetic flux at high β (right hand side panels). Lower figure: light curves for inefficient advection at large β. The initial constant magnetic field in this case is Bz0 = 20 G.

In the text

	Fig. B.1. Dwarf novæ light curves as a function of time for different values of βeq. We used Bz0 = 10, 4 and 2.5 G for the panels a–c, respectively.
In the text