© ESO 2021

1. Introduction

During the collapse of a molecular cloud, a protostar and its surrounding protostellar disk are formed. Observations of such star–disk systems reveal that the luminosity of the central object varies with time (e.g. Herbig 1989; Contreras Peña et al. 2017). Some of those young variable stars can be classified as FU Orionis (FU Ori) objects, which show a rise in luminosity over several orders of magnitudes for a short period of time (e.g., Hartmann & Kenyon 1996; Audard et al. 2014). Spectroscopic observations by Zhu et al. (2007) and Eisner & Hillenbrand (2011) indicate that such high-luminosity phases (bursts) are associated with an increased accretion rate of disk material onto the protostar. The frequency of FU Ori-type eruptions in the solar neighborhood suggests that a protostellar accretion burst is not a single but a recurring phenomenon (Hartmann & Kenyon 1996). Moreover, numerical simulations indicate that an FU Ori-like star–disk system undergoes at least 10 to 20 such eruptive events (Hartmann & Kenyon 1996).

Thermal instability (TI) of the inner disk provides a possible explanation for the aforementioned eruptive episodic accretion events. Ionization of most of the hydrogen in the very inner parts of a disk (10 R_⊙ ≲ r ≲ 0.1 AU) leads to enhanced viscosity caused by runaway heating of the optically thick inner disk regions and therefore to an increased accretion rate (e.g., Bell & Lin 1994). Although their model is able to predict bursts, the explanation of the burst duration is unsatisfying because of the narrow radial region in which the disk can become sufficiently ionized and hence may become thermally unstable (Zhu et al. 2009a). To expand the potentially thermally unstable regions of a disk, a model with vertical layers of different viscosity has been proposed (e.g., Gammie 1996). A warm, irradiated surface layer is sufficiently ionized for the magneto-rotational instability (MRI) to efficiently operate, but also acts as a shield for the deeper disk layer close to the midplane. Hence, MRI cannot develop in the deep layer and the subsequent reduced viscosity leads to a pile-up of material from the outer disk (dead zone). Eventually, the gas temperature in the dead zone rises above a certain threshold, caused by enhanced viscous heating of the accumulated mass. The deep disk layer then becomes MRI unstable, which may be sufficient to trigger the TI in the dead zone, effectively broadening the radial region where the disk can become thermally unstable. Recent work (e.g., Armitage et al. 2001; Zhu et al. 2009a, 2010) assumed gravitational instabilities (GIs) to act as an effective mechanism to enhance the mass transport rate from the outer disk towards the inner regions. However, observations of low-luminosity accretion events –identified as FU Ori events (see e.g., Kóspál et al. 2016; Hillenbrand et al. 2018)– show that low-mass protostellar star–disk systems with masses lower than 0.01 M_⊙ can also undergo episodic accretion. Those systems are, contrary to FU Ori for example, not in the embedded (Class I) but in the late Class II phase (cf. Kóspál et al. 2016) and additionally have disks with masses M_disk ≈ 0.01 M_⊙, which are too low for GI to contribute significantly (even in the outer disk) to the averaged mass transport rate.

As TI is likely to develop in the very inner disk parts, especially for low-mass disks, the position of the inner disk rim is a crucial parameter. Magneto-hydrodynamic (MHD) simulations in 3D show that protoplanetary disks are eventually disrupted by strong protostellar magnetic fields (e.g., Romanova et al. 2004; Bouvier et al. 2007; Romanova & Kurosawa 2014; Romanova & Owocki 2015). The accretion flow is then funneled along magnetic field lines onto the star. This is especially important in the case of episodic accretion, as then the inner disk rim is moving towards the star (e.g., Hartmann et al. 2016; Zhu et al. 2020) and consequently the inner disk dynamics is altered. To be able to consider the influence of stellar magnetic torques on the disk evolution and simultaneously calculate the radial position of the inner disk radius self-consistently, it is necessary to solve the full set of hydrodynamic equations.

The work of Bell & Lin (1994), Armitage et al. (2001) and Zhu et al. (2007, 2008, 2009a) for example is able to explain FU Ori-like outbursts by either purely thermally unstable inner disks or by including a layered viscosity model. Due to time-step limitations of explicit methods, especially in disk regions close to the star (r ≤ 0.5 AU) (Appendix A), neither the influence of protostellar magnetic fields nor the effects of local steep pressure gradients could be reproduced. However, the findings of Bouvier et al. (2007) and Romanova & Kurosawa (2014), for example, indicate that the very inner regions of disks are crucial for understanding episodic accretion. Recent hydrodynamic models of MRI and thermal instability bursts (Bae et al. 2013; Kadam et al. 2020; Vorobyov et al. 2020) have strong limitations in modeling the innermost disk regions, therefore necessitating the development of numerical codes that can realistically treat the star–disk interface in global disk simulations. Combining the requirements of (i) self-consistent treatment of both the gas dynamics in the inner disk including the effects of a stellar magnetic field and the position of the inner disk radius, (ii) the inclusion of an energy equation to investigate the effects of TI-induced episodic accretion, and (iii) maintaining the ability to resolve the very inner disk parts of a few stellar radii and carry out long-term calculations up to a few million years, we choose to use an implicit 1+1D RHD code (TAPIR, see Stoekl & Dorfi 2014; Ragossnig et al. 2020) that enables us to meet all the requirements without being time-step-limited by the Courant-Friedrichs-Levy (CFL) condition (Appendix A). Until now, no comparable numerical code has combined the above requirements. We also compare our simulations to similar work (e.g., Zhu et al. 2010, 2020) and point out the differences due to a time-dependent and self-consistent treatment of the component u_φ(r, t).

In Sect. 2 we briefly describe the physical equations and the layered viscosity model. Details of our method and the implementation of the boundary conditions used are explained in Sect. 3. The position of the inner disk edge as well as the importance of a consistent treatment of the inner disk region is outlined in Sect. 4. Our results are presented and discussed in Sect. 5 and are then summarized in Sect. 6.

2. Physical setup

Protoplanetary disks and their evolution in time can be described by adopting the equations of hydrodynamics using the turbulent viscosity prescription introduced by for example Shakura & Sunyaev (1973) or Balbus & Hawley (1991) and an energy equation. This set of equations is briefly presented in Sect. 2.1. Additionally, the stellar magnetic field (Sect. 2.2), a description of the energy equation (Sect. 2.3), the adapted viscosity model (Sect. 2.4), and the gas and dust opacities (Sect. 2.5) are presented.

2.1. Theoretical framework

The following equations describe the time evolution of a viscous accretion disk in the thin-disk limit (pressure scale-height H_p ≪ radius r) around a protostar. The energy equation contains radial and vertical radiation transport in an approximated way (see Sect. 2.3) and incorporates irradiation from the central object. Because of the thin-disk approximation, the equations can be vertically integrated and read as,

$$\begin{array}{cc}& \frac{\partial }{\partial t}\mathrm{\Sigma }+\mathrm{\nabla }·(\mathrm{\Sigma }\mathit{u})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial }{\partial t}(\mathrm{\Sigma }\mathit{u})+\mathrm{\nabla }·(\mathrm{\Sigma }\mathit{u}:\mathit{u})-\frac{B_z\mathit{B}}{2\pi }\hfill \\ \hfill & +\mathrm{\nabla }{P}_{\mathrm{gas}}+\mathrm{\nabla }·Q+\mathrm{\Sigma }\mathrm{\nabla }\psi +H_{\mathrm{p}}\mathrm{\nabla }\left(\frac{B_z^2}{4\pi }\right)=0,\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial }{\partial t}(\mathrm{\Sigma }e)+\mathrm{\nabla }·(\mathrm{\Sigma }\mathit{u}e)+P_{\mathrm{gas}}\mathrm{\nabla }·\mathit{u}\hfill \\ \hfill & +Q:\mathrm{\nabla }\mathit{u}-4\pi \mathrm{\Sigma }\kappa (J-S)+{\dot{E}}_{\mathrm{rad}}=0,\hfill \end{array}$$(3)

where Σ, P_gas, and u denote the gas column density, the vertically integrated gas pressure, and the gas velocity vector, respectively. Q represents the viscous pressure tensor and ψ marks the gravitational potential of the system (star and disk). In most of our calculations, we focus on late class II systems, in which the gravitational force from the protostar dominates. The protostellar magnetic field is presented by B = (B_r, B_φ, B_z)^T and is modeled via a prescribed field in polodial and torodial directions (see Sect. 2.2). The term B_z B/2 π corresponds to magnetic stress exerted on the disk gas, whereas $\nabla (B_z^2/4\pi )$ describes the magnetic pressure force density. In Eq. (3), e denotes the specific internal energy density, J stands for the zeroth moment of the radiation field and S represents the source function. The (J − S) in Eq. (3) is treated approximately as described in Sect. 2.3. To describe the internal structure of the ideal gas in radial and vertical directions, we utilize the ideal equation of state (EOS),

$$\begin{array}{c}\hfill P_{\mathrm{gas}}=(\gamma -1)e,\end{array}$$(4)

where γ denotes the adiabatic coefficient.

2.2. Stellar magnetic field

Livio & Pringle (1992) used an analytic approach to describe a stellar magnetic field which threads an accretion disk,

$$\begin{array}{cc}& B_z(r)=B_{\star }{\left(\frac{R_{\star }}{r}\right)}^3,\hfill \end{array}$$(5)

$$\begin{array}{cc}& B_{\phi }^{\mathrm{I}}(r)\simeq B_z(r)(1-\frac{\mathrm{\Omega }(r)}{{\mathrm{\Omega }}_{\star }}),\hfill \end{array}$$(6)

where B_z is the vertical magnetic field component and is assumed to have a positive sign. The stellar magnetic field is modeled in Eqs. (5) and (6) as a dipole field and a toroidal component $B_{\phi }^{\mathrm{I}}$ due to the torque exerted on the threaded protoplanetary disk by the stellar magnetic field, respectively. However, Eq. (6) has the problem that |B_φ|≫|B_z| very close to the star. Therefore, various authors (e.g., Rappaport et al. 2004; Kluźniak & Rappaport 2007) use a slightly modified version of Eq. (6), which reads,

$$\begin{array}{cc}& {\alpha }_{\mathrm{cor}}=\{\begin{array}{cc}+1\hfill & \mathrm{if}r<r_{\mathrm{cor}}\hfill \\ -1\hfill & \mathrm{if}r>r_{\mathrm{cor}}\hfill \end{array},\hfill \end{array}$$(7)

$$\begin{array}{cc}& B_{\phi }^{\mathrm{II}}(r)\simeq -{\alpha }_{\mathrm{cor}}{B}_z(r)[1-{\left(\frac{{\mathrm{\Omega }}_{\star }}{\mathrm{\Omega }(r)}\right)}^{{\alpha }_{\mathrm{cor}}}]·\hfill \end{array}$$(8)

The superscripts I and II are applied to differentiate between the two models for B_φ, whereas α_cor takes account of changing the field prescription to inside and outside the corotation radius. The model described by $B_{\phi }^{\mathrm{II}}$ confers the advantage that |B_φ|≈|B_z| for radii not close to r_cor. This is also physically motivated, as larger |B_φ| are not possible because a too tightly wound up toroidal field would lead to an opening up of the field configuration because of the increased magnetic energy (see e.g., Rappaport et al. 2004). Typical values for the protostellar magnetic field are of the order of several kG (see e.g., Bouvier et al. 2007; Kurosawa et al. 2008).

2.3. Energy equation

The energy equation Eq. (3) is described by Ragossnig et al. (2020) and accounts for viscous heating, radiative cooling, and radial diffusion of inner energy. Furthermore, a vertical radiation transport is included to also account for stellar irradiation on to the surface of the disk. In cylindrical coordinates, Eq. (3) reads,

$$\begin{array}{cc}& Q=\frac{\mu }{2}[\mathrm{\nabla }\mathit{u}+{(\mathrm{\nabla }\mathit{u})}^{\mathrm{T}}-\frac{2}{3}(\mathrm{\nabla }·\mathit{u})\mathbb{1}],\hfill \end{array}$$(9)

$$\begin{array}{cc}& {ϵ}_Q\equiv 2Q:\mathbf{\nabla }\mathit{u},\hfill \end{array}$$(10)

$$\begin{array}{cc}& \frac{\partial }{\partial t}(\mathrm{\Sigma }e)+\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}_{\mathrm{r}}\mathrm{\Sigma }e\right)+P\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}_{\mathrm{r}}\right)+{ϵ}_Q+{\dot{E}}_{\mathrm{rad}}\hfill \\ \hfill & -4\pi \frac{\sigma }{3r\pi }\frac{\partial }{\partial r}\left(\frac{r}{{\kappa }_{\mathrm{R}}\rho }\frac{\partial T_0^4}{\partial r}\right)=0,\hfill \end{array}$$(11)

where the second and third terms in Eq. (11) represent advection due to accretion and work done by pressure, respectively. Here, ϵ_Q denotes the viscous energy dissipation (e.g., Tscharnuter & Winkler 1979). The last term in Eq. (11) describes radial radiative transport as a diffusion approximation in the Eddington limit (e.g., Ragossnig et al. 2020), whereas Ė_rad depicts the net radiation heating or cooling rate per unit surface area of the protoplanetary disk. This latter can be obtained by solving the vertical, stationary radiation transfer equation in the optically thick limit (τ → ∞) assuming local thermal equilibrium and utilizing an Eddington factor f_edd = 1/3,

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{rad}}=\frac{8\sigma }{3\tau }(T_0^4-T_{\mathrm{surf}}^4).\end{array}$$(12)

Ė_rad is a balance of protostellar irradiation Ė_irr, radiative cooling Ė_cool, and irradiation from the ambient medium Ė_amb,

$$\begin{array}{cc}& {\dot{E}}_{\mathrm{rad}}+{\dot{E}}_{\mathrm{irr}}+{\dot{E}}_{\mathrm{amb}}-{\dot{E}}_{\mathrm{cool}}=0,\hfill \end{array}$$(13)

$$\begin{array}{cc}& {\dot{E}}_{\mathrm{irr}}=\frac{L_{\star }}{r^2\pi }{f}_{\mathrm{irr}}max[\mathrm{\Delta }\left(\frac{H_{\mathrm{p}}-H_{\star }}{r}\right),0],\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\dot{E}}_{\mathrm{amb}}=2\sigma T_{\mathrm{amb}}^4,\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\dot{E}}_{\mathrm{cool}}=2\sigma T_{\mathrm{surf}}^4,\hfill \end{array}$$(16)

where L_⋆ and H_⋆ denote stellar luminosity and an effective stellar radius, respectively. Furthermore, black-body irradiation is assumed for the cooling term and the ambient radiation. The final form reads (Ragossnig et al. 2020),

$$\begin{array}{cc}& {\dot{E}}_{\mathrm{rad}}=\sigma \frac{1}{1+{\tau }^{\prime }}(T_0^4-T_{\mathrm{amb}}^4)-\frac{L_{\star }}{r^2\pi }{f}_{\mathrm{irr}}max[\mathrm{\Delta }\left(\frac{H_{\mathrm{p}}-H_{\star }}{r}\right),0],\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\tau }^{\prime }\equiv \frac{3}{4}\tau ,\hfill \end{array}$$(18)

where f_irr determines how much of the irradiation can be processed by the gas and which fraction gets reflected.

2.4. Viscosity model

The kinematic viscosity ν is defined as (Shakura & Sunyaev 1973)

$$\begin{array}{c}\hfill \nu =\alpha c_{\mathrm{S}}{H}_{\mathrm{P}},\end{array}$$(19)

where α is the viscous parameter, c_S the isothermal sound speed, and H_p the pressure scale height. We apply the thin-disk approximation (Armitage 2010) which allows a vertical integration of the physical disk quantities via the EOS,

$$\begin{array}{c}\hfill P_{\mathrm{gas},0}={\rho }_0{c}_{\mathrm{S}}^2,\end{array}$$(20)

where P_gas, 0 is the midplane gas pressure. The vertically integrated dynamical viscosity μ is then written as

$$\begin{array}{c}\hfill \mu =\alpha c_{\mathrm{S}}{H}_{\mathrm{P}}\mathrm{\Sigma }.\end{array}$$(21)

For the viscosity parameter α we adopt the layered-viscosity model (Gammie 1996), which is a sum of contributions of vertically stratified disk layers,

$$\begin{array}{c}\hfill \alpha ={\alpha }_{\mathrm{base}}+{\alpha }_{\mathrm{surf}}+{\alpha }_{\mathrm{deep}}+{\alpha }_{\mathrm{grav}},\end{array}$$(22)

where α_base is a base value for the viscosity (e.g., within a dead-zone) and α_surf is the MRI viscosity in the permanently active surface layer, which is ionized by cosmic rays and stellar irradiation for example. α_deep accounts for the viscosity in parts of the disk, where the temperature exceeds the threshold for thermal ionization (e.g., Zhu et al. 2009a). The parameter α_grav denotes GIs acting as effective viscosity.

The separation of the surface and deep layer is controlled by the local surface density. If the surface density at an arbitrary radius exceeds a value Σ₀, the disk develops a deep disk layer at that radius. To handle this, we introduce a density switch that is s_Σ = 1 if Σ ≥ Σ₀ and s_Σ = 0 if Σ ≤ Σ₀. Parts of the deep layer then become MRI active if the local gas temperature T₀ is higher than a constant temperature threshold T_active for thermal ionization. Hence, the viscosity parameter for the surface layer reads

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{surf}}(r)={\alpha }_{\mathrm{MRI}}[s_{\mathrm{\Sigma }}\frac{{\mathrm{\Sigma }}_0}{\mathrm{\Sigma }(r)}+(1-s_{\mathrm{\Sigma }})],\end{array}$$(23)

where α_MRI (e.g., Zhu et al. 2010; Hartmann & Bae 2018) is a constant value for MRI-active regions. We adopt Σ₀ = 100 g cm⁻² for all calculations in this work. Additionally, the viscosity parameter for the deep layer is temperature dependent, that is,

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{deep}}(r)={\alpha }_{\mathrm{MRI}}{s}_{\mathrm{\Sigma }}(1-\frac{{\mathrm{\Sigma }}_0}{\mathrm{\Sigma }(r)})s(T_0),\end{array}$$(24)

where s(T₀) is a temperature switch,

$$\begin{array}{c}\hfill s(T_0)=\frac{1}{2}[1+tanh\left(\frac{T_0-T_{\mathrm{active}}}{T_{\mathrm{width}}}\right)],\end{array}$$(25)

which controls whether the deep layer is MRI-active or not (see e.g., Flock et al. 2016). For faster convergence of the Newton-Raphson iteration method (see Sect. 3), Eq. (25) is utilized to establish a smooth transition between active and inactive parts within the deep layer, where T_width = 10 − 50 K is a smoothing width for the temperature. This smooth transition creates a temperature zone below T_active in which viscosity is already increased and the disk becomes MRI unstable. This MRI active zone starts at ≈T_active − 2 ⋅ T_width.

We define the viscosity parameter for gravitationally unstable regions as

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{grav}}={\alpha }_{\mathrm{GI}}{s}_{\mathrm{G}}(\frac{Q_{\mathrm{T},\mathrm{crit}.}^2}{Q_{\mathrm{T}}^2}-1),\end{array}$$(26)

where α_GI = 0.01 is a constant value for gravitationally unstable regions within the disk, and Q_T is the Toomre parameter (Toomre 1964). The critical value for the Toomre parameter at which GIs occur is set to Q_T, crit = 1.0. The quantity s_G is a switch for the GI and is s_G = 1 if the disk is unstable (Q_T < Q_T, crit.) and s_G = 0 otherwise. In this paper we have chosen disks with masses low enough (M_disk < 0.1 M_⋆) to not exceed the Toomre criterion, thus avoiding GIs within the disk.

2.5. Gas and dust opacities

For the gas opacity description κ_R, gas we use the Rosseland-mean opacity tables created by Ferguson et al. (2005), who compiled low-temperature opacity tables by employing opacity sampling methods with a considerably higher sampling resolution than the previous opacity tables provided by Alexander & Ferguson (1994). In this work we use opacities for solar abundance (X = 0.7, Z = 0.02) based on Caffau et al. (2011) in the temperature range from 500 K to 30 000 K.

Opacities below 500 K are dust-dominated, and therefore we use the dust opacity model introduced by Pollack et al. (1985) for very low temperatures. Dust grains are believed to consist of silicates, iron, troilite, organics, and ice (see e.g., Pollack et al. 1994; Henning & Stognienko 1996), which can have different metal abundances and shapes. Pollack et al. (1985) divided silicates into iron-poor, iron-rich, and normal species, which have Fe/(Fe + Mg) ratios of 0.0, 0.4, and 0.3, respectively. These latter authors also considered different shapes as spherical and aggregate grains, as well as a distinction between homogeneous, composite, and porous dust grains. For our calculations, the dust-to-gas ratio f_dust is kept constant in the whole disk and in time, which is an oversimplification. However, we do not treat gas and dust separately in our simulations and cannot calculate a time-dependent and radially changing f_dust. The total Rosseland-mean opacity κ_R for various midplane gas densities ρ₀(r) is shown in Fig. 1 and can be calculated by combining gas and dust opacities,

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{R}}={\kappa }_{\mathrm{R},\mathrm{gas}}+f_{\mathrm{dust}}{\kappa }_{\mathrm{R},\mathrm{dust}}.\end{array}$$(27)

Fig. 1. Rosseland-mean opacities κR for various midplane gas densities ρ0.

3. Method description

The reasons for using a 1D implicit integration method in this work are briefly explained in Sect. 3.1 and the method is described in Sect. 3.2. Some vital details of our numerical method are briefly presented in Sect. 3.3, whereas the boundary conditions used in our simulation are presented in Sect. 3.4. An implicit numerical integration scheme needs an initial model to start the simulation with. The construction of such a model is briefly outlined in Sect. 3.5.

3.1. Justification of the 1D approach

In recent years, the advent of 2D and 3D models of protoplanetary disks (e.g., Zhu et al. 2009b, 2020; Vorobyov & Basu 2007; Vorobyov & Pavlyuchenkov 2017; Kadam et al. 2019) has lead to a detailed understanding of the complicated interactions of equally important physical processes abundant in such star–disk systems. Hence, the question must be asked, what new insights can be added by another 1D approach.

The study of magnetically truncated, self-consistently simulated protoplanetary disks cannot (at the time of writing) be combined with an investigation of the influence of magnetic fields on the long-term disk behavior. This is because of severe time-step limitations for simulations of the inner disk. An implicit integration scheme (see Sect. 3.2) is not limited to those restrictions (see Appendix A).

We understand that global protoplanetary disk simulations outside the very inner regions (r > 0.5 AU) should be at least 2D for studies of the gravitationally unstable phases of disk evolution (e.g., Vorobyov et al. 2019), and 3D for a detailed understanding of complex features such as turbulence for example (e.g., Flock et al. 2016). However, in the region within approximately 0.5 AU, those models are confined to coarse simplifications (e.g., to the viscous diffusion ansatz or an inner “smart” sink-cell approach) and cannot include important features such as the effects of torques caused by a protostellar magnetic field. Besides time-step limitations imposed by the CFL condition, multi-dimensional models are either confined to a limited, specific simulation region compared to the whole disk dimensions (e.g., Vorobyov et al. 2019) or tend to evolve towards a quasi-stationary state after several orbital periods of the inner disk (e.g., Romanova et al. 2004; Zhu et al. 2020). Both cases are not ideal for investigating the long-term effects over several million years in the innermost disk region. Additionally, the increased shearing forces towards smaller radii tend to smooth out instabilities in angular directions (see e.g., Lesur 2020) in the very inner disk.

Using an implicit numerical method allows consistent treatment of the inner disk by including the effects of a stellar magnetic field while simultaneously maintaining the ability to carry out long-term evolution studies. Additionally, the self-consistent treatment of large-scale magnetic fields during protoplanetary disk evolution as well as the interaction of the stellar magnetic field and the inner disk makes it necessary to solve the toroidal velocity component u_φ in a time-dependent manner in order to cover magnetic torques (see e.g., Lubow et al. 1994; Guilet & Ogilvie 2012, 2014). Disk winds are not included in our simulations but will be addressed separately in a follow-up paper. Apart from a more realistic handling of the radial gas flow, this further supports the choice of solving all hydrodynamic quantities for long-term evolution studies of protoplanetary accretion disks.

The hydrodynamic equations for protoplanetary disks are formulated as a boundary value problem. It is equally true for explicit and implicit methods that boundary conditions for the column density Σ, the gas velocity u and the internal energy e imposed on the disk at the inner and outer boundaries determine the radial structure and influence its evolution. An implicit time integration on the other hand will not converge towards a solution at a new time, if the boundary conditions are chosen such that no physical radial structure for Σ, u, and e can be found for the equations involved. Therefore, an implicit method, despite being more difficult to set up, provides feedback with respect to the physical correctness of the boundary conditions used (see Dorfi & Drury 1987; Stoekl & Dorfi 2014, for further details).

Consequently, a 1D implicitly integrated simulation of the inner disk can indeed lead to a better understanding of important features and could for example be used as the inner boundary for more sophisticated 2D and 3D models by providing a physically more accurate description of the very inner regions (e.g., Crida et al. 2007).

3.2. Implicit integration scheme

A numerical method utilizing an implicit integration scheme (e.g., Dorfi & Drury 1987; Stoekl & Dorfi 2014) requires a substantial amount of additional work in development compared to explicit methods. This is because, for each iteration in time Δt, a Jacobian has to be constructed, which consists of the derivations of each variable at a certain grid point with respect to each variable at every grid point in the computational domain. This matrix must then be subsequently inverted as often as a multidimensional Newton-Raphson iterator needs to converge towards a new solution of the problem at a new time t + Δt. Another complication is the requirement of an initial model that already solves the discretized set of algebraic equations. The additional work for an implicit integration scheme has the advantage of not being time-step-limited by the Courant-Friedrichs-Levy (CFL) condition. This translates into a condition that information is not transported over more than one grid cell during a single time-step. This constraint effectively limits the time-step of every hydrodynamic disk simulation with an explicit time-integration scheme. Comparing the time-step of our implicit method to the constraint dictated by the CFl condition (see Appendix A), we can confirm that the simulation time would be increased by several orders of magnitude to achieve the same radial resolution (2000 radial points are used throughout the simulations in this paper) in the inner regions of the disk with an explicit method. Consequently, full-blown 3D MHD calculations are very time-consuming and can only be conducted for between a few hundred and a maximum of one thousand orbital periods (see e.g., Romanova & Kurosawa 2014; Zhu et al. 2020), which would result in an inner disk radius of r_in = 0.06 AU for a covered simulation time-period of t_total ≤ 15 years.

3.3. Numerical method details

The equations are formulated in one-dimension in the radial direction, are axisymmetric and conservative, and are integrated using a finite-volume method. Also, a van Leer advection scheme (van Leer 1977) is used at the cell boundaries. This ensures that the method is second order in space, while the same accuracy in time is achieved by using time-centered variables (for details see e.g., Dorfi 1998; Dorfi et al. 2006).

The equations are discretized by using a staggered mesh (see Ragossnig et al. 2020), where scalar entities as column density Σ and internal energy e are discretized on a scalar mesh while vector-like quantities as the gas velocity component u_r are discretized on the vector mesh, which is spatially displaced from the scalar mesh by half a grid cell. The velocity component u_φ is of vectorial nature, because for the degeneration of grid cell faces in the angular direction for a 1D radial treatment, u_φ has to be discretized at the scalar mesh. In addition to the numerical treatment of the physical equations, we use an adaptive numerical grid. We apply the grid point distribution introduced by Dorfi (1998), which avoids the need for coping with artificially introduced perturbations due to grid adaptations, because the mesh configuration is simultaneously solved with the equations. The specific details of our numerical method as well as the detailed discretizations performed for the physical set of equations presented in equations Eqs. (1)–(3) and the implementation of the grid equation can be studied in more detail in Ragossnig et al. (2020).

3.4. Numerical boundary conditions

Protoplanetary disk quantities such as surface density structure Σ(r), disk mass M_disk, and accretion flow Ṁ(r) sensitively depend on the choice of the inner and outer boundary conditions imposed on the disk. At the inner boundary, magnetic torques caused by the stellar magnetic field brake the disk in the angular direction, and consequently the radial flow starts to accelerate towards the star. At a certain point, the magnetic field dominates the disk dynamics and the accretion transitions from a radial drift into an accretion stream along magnetic funnels. The simulation of this magnetically dominated region is intrinsically 3D and cannot be tackled by a 1D approach (see e.g., Bouvier et al. 2007; Romanova & Kurosawa 2014). Hartmann et al. (2016) argued that the disk approximation fails at the magnetic truncation radius r_trunc and is therefore chosen as our inner boundary r_in as defined in Eq. (47).

A self-consistent treatment of the inner boundary requires the angular velocity to adapt to the magnetic torque of the stellar magnetic field. This is best represented with a zero gradient condition,

$$\begin{array}{cc}& {\frac{\partial u_{\phi }}{\partial r}|}_{r=r_{\mathrm{in}}}=0,\hfill \end{array}$$(28)

$$\begin{array}{cc}& {\frac{\partial u_{\mathrm{r}}}{\partial r}|}_{r=r_{\mathrm{in}}}=0.\hfill \end{array}$$(29)

The gas temperature structure at the inner boundary is chosen such that the sum of contributions of viscous heating and stellar irradiation determines the gas temperature at the inner boundary, which is mapped by employing a Van Neumann boundary condition for the midplane gas temperature T₀ at the inner boundary,

$$\begin{array}{c}\hfill {\frac{\partial T_0}{\partial r}|}_{r=r_{\mathrm{in}}}=0.\end{array}$$(30)

At the inner boundary, the surface density Σ_in must be able to adapt to different mass transport rates Ṁ_in, for example during an episodic accretion event. Hence, we choose a Van Neumann boundary condition for Σ_in,

$$\begin{array}{c}\hfill {\frac{\mathit{\partial }\mathbf{\Sigma }}{\mathit{\partial }\mathit{r}}|}_{\mathit{r}\mathbf{=}{\mathit{r}}_{\mathbf{in}}}\mathbf{=}\mathbf{0}.\end{array}$$(31)

At the outer boundary, a mass flux Ṁ_out is assumed to act as an external mass reservoir for the disk.

$$\begin{array}{cc}& {\frac{\partial \mathrm{\Sigma }}{\partial r}|}_{r=r_{\mathrm{out}}}=0,\hfill \end{array}$$(32)

$$\begin{array}{cc}& u_{\mathrm{r}}(r_{\mathrm{out}})=\frac{{\dot{M}}_{\mathrm{out}}}{2\pi r_{\mathrm{out}}\mathrm{\Sigma }(r_{\mathrm{out}})}·\hfill \end{array}$$(33)

The angular velocity u_φ is assumed to be Keplerian and a Van Neumann boundary condition is used for the internal energy e at the outer boundary condition, which results in

$$\begin{array}{cc}& {\frac{\partial e}{\partial r}|}_{r=r_{\mathrm{out}}}=0,\hfill \end{array}$$(34)

$$\begin{array}{cc}& u_{\phi }(r_{\mathrm{out}})={\mathrm{\Omega }}_{\mathrm{K}}(r_{\mathrm{out}})r_{\mathrm{out}},\hfill \end{array}$$(35)

where Ω_K(r_out) is the Keplerian velocity at the outer boundary. In this study, we focus on the star–disk interaction, and therefore an outer radius of 30 AU is chosen for all simulations.

3.5. Stationary initial model

An implicit numerical method requires an initial model to start with, which already solves the discretized equations along with the grid equation described by Dorfi (1998). Although, in principle, every solution to the set of equations can be used as a starting model, we aim to construct stationary solutions to start our simulations. Such solutions can be used for verification of the numerical method, because analytical results exist for the assumption of time-independence (see e.g., Armitage 2010).

We construct these stationary solutions in a two-stage process. At first, we choose a Keplerian disk with no radial velocity u_r (see Fig. 2). This disk is in local thermal equilibrium (LTE) in the vertical and radial directions. Second, we introduce a time-independent mass flow across the outer boundary Ṁ_out. Collectively, a certain outer mass flux, other boundary conditions, stellar parameters (mass M_⋆, luminosity L_⋆, magnetic field B_⋆, rotation period P, and radius R_⋆), and viscosity model parameters (cf. Sect. 2.4) fully determine the internal radial structure of the disk.

Fig. 2. Snapshots in time, which show the convergence towards a stationary initial model for (a) the surface density Σ(r), (b) the radial velocity ur(r), and (c) the mass transport rate Ṁ. The black dashed line shows the initial Keplerian disk, and the black solid line denotes the final stationary model.

The transition from the starting model to the stationary model occurs on the viscous timescale t_visc,

$$\begin{array}{c}\hfill t_{\mathrm{visc}}(r)=\frac{r^2}{\nu (r)}·\end{array}$$(36)

In Fig. 2 this transition can be seen as mass is transported through the disk, fed by Ṁ_out. After t ⪆ t_visc(r_out), the disk becomes stationary and has adjusted to the boundary conditions and the mass transport rate of the stationary initial model Ṁ_init is constant throughout the whole disk, and therefore Ṁ_init is equal to Ṁ_out (cf. panel c of Fig. 2). We also note that r_in adjusts accordingly for a changing value of Ṁ (cf. Fig. 2 and Sect. 4.3).

We want our simulations to trigger thermal instability and to develop an episodic accretion onto the protostar. During construction of the initial model, this behavior is undesired, and therefore outbursts have to be prevented. We achieve this by increasing the activation temperature T_active during the first stage to a value that the disk gas temperature T₀ does not surpass at any radius. As the model becomes stationary, the time-step can rise due to the implicit nature of the numerical method. Contrary to explicit numerical methods, time-steps can become larger than the viscous timescale t_visc, because they are not restricted by the CFL condition. This effectively subdues the thermal instability, which develops on a much shorter thermal timescale. Therefore, after Δt > t_visc(r_out), the initial model construction enters the second stage, where we smoothly decrease T_active to the desired value and let the disk model become stationary. We note that this procedure is not physical but is utilized to generate a quasi-stationary initial model, which solves the equations of hydrodynamics.

4. The importance of the inner disk region

The limitations imposed on explicit numerical methods by the CFL condition (see Sect. 3.3) are mitigated in Bell & Lin (1994) and subsequent works (e.g., Armitage et al. 2001; Zhu et al. 2009a) by using a diffusion ansatz with a low spatial resolution (compared to our approach) for describing the radial disk evolution. In Sect. 4.1 the differences between a full hydrodynamic simulation and the diffusion ansatz are investigated. The importance of the pressure force is evaluated in Sect. 4.2 and the position of the inner disk boundary is discussed in Sect. 4.3.

4.1. Revisiting the diffusion equation

The radial evolution equation reads (e.g., Pringle 1981),

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Sigma }}{\partial t}=-\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{(r^2\mathrm{\Omega }{)}^{\prime }}\frac{\partial }{\partial r}(\nu \mathrm{\Sigma }{r}^3{\mathrm{\Omega }}^{\prime })\right],\end{array}$$(37)

where Ω(r) corresponds to the angular velocity. In the derivation of Eq. (37), pressure gradients are not included in the angular momentum equation (e.g., Pringle 1981). If Ω(r)∝Ω_K(r), with Ω_K(r) depicting the Keplerian velocity, Eq. (37) becomes

$$\begin{array}{cc}& \frac{\partial \mathrm{\Sigma }}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[\sqrt{r}\frac{\partial }{\partial r}(\nu \mathrm{\Sigma }\sqrt{r})\right],\hfill \end{array}$$(38)

$$\begin{array}{cc}& u_{\mathrm{r}}=-\frac{3}{\mathrm{\Sigma }\sqrt{r}}\frac{\partial }{\partial r}(\nu \mathrm{\Sigma }\sqrt{r}),\hfill \end{array}$$(39)

where the continuity equation Eq. (1) is used to obtain a formulation for the radial drift velocity u_r. Torques due to for example protostellar magnetic fields cannot be treated in a self-consistent way with Eq. (38), because this would involve braking (or acceleration) of the angular velocity component in a way in which Ω(r) is no longer proportional to Ω_K(r). In order to see the impact of a perturbation in the angular velocity profile Ω(r), one can model the angular velocity as Keplerian with a small linear perturbation Ω₁,

$$\begin{array}{c}\hfill \mathrm{\Omega }(r)={\mathrm{\Omega }}_{\mathrm{K}}(r)+{\mathrm{\Omega }}_1(r).\end{array}$$(40)

After linearization of Eq. (37) with respect to Ω₁, the following expression is obtained:

$$\begin{array}{cc}& f_{\mathrm{corr}}=\frac{2(r^2{\mathrm{\Omega }}_1{)}^{\prime }}{r{\mathrm{\Omega }}_{\mathrm{K}}}+\frac{2(\nu \mathrm{\Sigma }{r}^3{\mathrm{\Omega }}_1^{\prime }{)}^{\prime }}{3(\nu \mathrm{\Sigma }{r}^2{\mathrm{\Omega }}_{\mathrm{K}}{)}^{\prime }},\hfill \end{array}$$(41)

$$\begin{array}{cc}& \frac{\partial \mathrm{\Sigma }}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[\sqrt{r}\frac{\partial }{\partial r}(\nu \mathrm{\Sigma }\sqrt{r})(1-f_{\mathrm{corr}})\right].\hfill \end{array}$$(42)

$$\begin{array}{cc}& u_{\mathrm{r}}=-\frac{3}{\mathrm{\Sigma }\sqrt{r}}\frac{\partial }{\partial r}(\nu \mathrm{\Sigma }\sqrt{r})(1-f_{\mathrm{corr}}),\hfill \end{array}$$(43)

where a primed quantity A′ denotes its radial derivative ∂A/∂r. The first term of the correction factor f_corr depends on the radial gradient of Ω₁; hence a radially localized, sharp perturbation in Ω(r) can lead to a non-negligible effect. During the onset of an FU Ori-like outburst, an ionization front is propagating radially inwards and outwards through the disk (for a detailed discussion see Bell & Lin 1994). Along these waves, the disk is adapting itself to a changed disk temperature and viscosity. Various authors (e.g., Bell & Lin 1994; Armitage 2010; Zhu et al. 2009a) have shown that those waves occur in the column density Σ(r), the internal energy profile e(r), and in the radial mass transport rate Ṁ(r). However, the ionization front is also visible as a wave in the radial angular velocity profile Ω(r), which then modifies the accretion rate and subsequently the density profile. The second term in Eq. (41) takes into account that a modified relative angular velocity between two neighboring radial disk annuli also changes the shear that those disk rings are experiencing. f_corr can be interpreted as the deviation from the viscous contribution to the diffusion equation. In Eq. (43) a correction factor f_corr > 1 leads to a radially localized, outward-bound mass flow, which is caused by redistribution of angular momentum in sharp angular velocity perturbations Ω₁(r).

In Sect. 5.1 we provide a detailed description of the perturbations during the onset of a FU Ori-like outburst, whose influence on f_corr is shown in Fig. 3. An angular velocity perturbation Ω₁ induced by sharp peaks of MRI efficiency due to density waves causes a radially localized peak in |f_corr|. Four snapshots are taken for four different times t₀ < t₁ < t₂ < t₃ immediately after the onset of TI and shows the inward-bound ionization front (cf. Sect. 5.1). The first term of Eq. (41) is plotted separately in Fig. 3 to show the direct influence of a nonKeplerian velocity u_φ compared to the additional influence of altered shear between two neighboring disk annuli. For a yet unpronounced perturbation at time t₀, the modified shear clearly dominates f_corr, but at later times t₁ to t₃ the direct influence of Ω₁ is also non-negligible. While the diffusion equation approach Eq. (38) appears to be well suited to describing the disk evolution in the absence of such waves, without solving the momentum equation in angular direction, there are, at least locally, substantial deviations to be expected.

Fig. 3. a: Ω1 waves for four different, arbitrarily chosen times t0 to t3 during the onset of a FU Ori-like outburst. b: correction factors as defined in Eq. (41) for the same times as in (a). The dashed-dotted lines correspond to the first term in Eq. (41), whereas the full lines show the total value of fcorr.

The deviations in Fig. 3 represent a snapshot during time evolution at a certain time t, which can affect the short term behaviour of the disk. A more qualitative way to investigate how strongly a certain disk region is affected by f_corr over a longer period of time is shown in Fig. 4. During the duration of an outburst (Δt_burst ∼ 20 years) in our fiducial model, we have added up the times in which the absolute value of f_corr is larger than a certain threshold (0.05, 0.5 and 1.0). In the innermost disk (≲0.1 AU) f_corr can exceed 0.05 over 90% of Δt_burst. Even values of f_corr = 0.50 and f_corr = 1.00 are exceeded during 80% of Δt_burst in some regions of the inner disk.

Fig. 4. Time relative to the burst duration in the fiducial model Δtburst, in which a region in the inner 0.5 AU exceeds a certain fcorr threshold (0.05, 0.5 and 1.0). The radial range is divided in 200 equally spaced logarithmic bins.

4.2. Importance of the pressure gradient

In a steady-state protoplanetary disk the gas rotation velocity u_φ (neglecting radial viscous forces, magnetic fields, and radial advection of momentum) is essentially given by

$$\begin{array}{cc}\hfill \frac{u_{\phi }^2}{r}& \approx \frac{G{M}_{\star }}{r^2}+\frac{1}{\rho }\frac{\mathrm{d}P}{\mathrm{d}r}·\hfill \end{array}$$(44)

Therefore, in order to have a subKeplerian steady-state accretion disk, the pressure gradient has to yield dP/dr < 0. For all radii but the very innermost parts (r < 0.1 AU), this condition is usually fulfilled. Furthermore, the pressure gradients are of minor order compared to the centrifugal forces, and at the inner boundary a zero-torque boundary condition is usually applied for such models in the absence of a stellar magnetic field. For example the IBC applied by Bell & Lin (1994),

$$\begin{array}{c}\hfill \mathrm{\Sigma }(r=r_{\mathrm{in}})\approx 0,\end{array}$$(45)

models a disk extending to the stellar surface and therefore accounts for the transition of a protoplanetary disk to u_r = 0. However, a magnetically braked inner disk has an increased accretion rate and therefore a lower column density Σ and diminished gas pressure P_gas, 0 close to the inner rim (see Fig. 5). Hence, in the very innermost regions (r ⪅ 0.2 AU) the pressure gradient force −∇P (cf. Eq. (2)) helps the gravitational force to push material inwards, which leads to a super-Keplerian u_φ according to Eq. (44). This can also be seen in Fig. 3, where Ω₁ > 0 (except for the declining slopes of the ionization fronts). The pressure gradient force also exceeds the viscous force contribution in the region of the ionization front (see Sect. 5.1). In the very inner regions, the diffusion equation ansatz of Eq. (38) is less suitable for describing the time evolution of the disk. This is especially true for thermally unstable disks (cf. Figs. 3 and 5), where due to localized perturbations of u_φ the deviations from the diffusion equation Eq. (38) are non-negligible (see Fig. 4).

Fig. 5. a: midplane gas pressure Pgas, 0 between 0.04 AU and 2 AU for two different times during an FU Ori like outburst. b: midplane pressure gradient force −∇Pgas, 0 relative to the gravitational force is plotted. The radially localized changes in sign occur at the positions of the ionization fronts running inwards and outwards. The deviation from the Keplerian angular velocity is shown in (c).

4.3. The position of the inner disk boundary

The inner disk is coupled to the star via the stellar magnetic field, which means that, inside the corotation radius r_cor, the angular momentum of the disk is transferred to the protostar, whereas outside of r_cor the disk is accelerated. This star–disk connection results in mutual transfer of angular momentum between the star and the inner disk (e.g., Matt et al. 2010). However, in the present study, the spin-up or spin-down of the star due to gain or loss of angular momentum is not considered, as our focus is on the influence of stellar magnetic torque on the inner disk. The inclusion of angular momentum transfer introduces a further complexity which makes it more difficult to separate the various influences altering the long-term evolution of a protoplanetary disk and has to be investigated in further studies.

The corotation radius is defined according to Kepler’s third law,

$$\begin{array}{c}\hfill r_{\mathrm{cor}}={\left(\frac{G{M}_{\star }{P}^2}{4{\pi }^2}\right)}^{1/3},\end{array}$$(46)

where P is the stellar rotational period. Following Herbst et al. (2001), we choose the stellar rotation period P = 6 days such that r_cor ≈ 0.06 AU for all simulations in this work. Magnetic torques tend to slow down the disk inside r_cor and the disk’s toroidal velocity u_φ becomes subKeplerian. Consequently, the disk material inside r_cor starts to radially accelerate towards the star until it reaches the magnetic truncation radius r_trunc. Inside r_trunc, the magnetic pressure of the stellar dipole field exceeds the ram pressure P_ram of the infalling material, and therefore the disk gets disrupted and accretes along magnetic funnels onto the star. Consequently, r_trunc is a natural choice for the inner disk edge in the case of strong accretion, where P_ram exceeds the gas pressure P_gas. Consequently, r_in is given as follows (for details cf. Hartmann et al. 2016):

$$\begin{array}{cc}\hfill r_{\mathrm{in}}(P_{\mathrm{gas}}<P_{\mathrm{ram}})\approx & 18\xi R_{\odot }{\left(\frac{B_{\star }}{{10}^3\mathrm{G}}\right)}^{4/7}{\left(\frac{R_{\ast }}{2R_{\odot }}\right)}^{12/7}{\left(\frac{M_{\star }}{0.5M_{\odot }}\right)}^{-1/7}\hfill \\ \hfill & \times {\left(\frac{{\dot{M}}_{\star }}{{10}^{-8}{M}_{\odot }/\mathrm{yr}}\right)}^{-2/7},\hfill \end{array}$$(47)

where B_⋆, R_⋆, and M_⋆ are the stellar magnetic field at the stellar surface, the protostellar radius, and its mass, respectively, and ξ is a correction factor that accounts for the rather complicated details of disk–star interactions and is usually set to ξ < 1 according to Hartmann et al. (2016). In this work, we choose ξ = 0.65. Further, Ṁ_⋆ denotes the accretion flow over the inner boundary onto the star and is in general time-dependent. As the inner radius r_in is dependent on the accretion rate Ṁ_⋆, it is likely to move inwards during increased accretion, for example during a FU Ori-like burst. The magnetic field cannot withstand the high accretion rate during an outburst event, and therefore the stellar magnetosphere gets squashed closer to the star. The movement of the inner boundary according to the time-dependent mass accretion rate Ṁ_⋆(t) for our fiducial model (see Table 1 and Sect. 5.1) during an outburst is shown in Fig. 6. During an outburst, r_in is pushed towards the stellar radius R_⋆ (cf. Eq. (47)). We want to emphasize that our method is able (from a numerical point of view) to cover the case of r_in ∼ R_⋆ in combination with long-term simulations.

Fig. 6. a: accretion rate for a full TI-triggered outburst for our fiducial model fid (cf. Table 1). b: variation of the inner boundary during the burst is shown in (a). Here, t0 corresponds to the onset of the first burst. The two dashed lines denote the protostellar radius R⋆ and the corotation radius rcor.

In regions where the ram pressure P_ram exceeds the gas pressure P_gas, Eq. (47) describes the movement of r_in well. If P_gas is larger than P_ram, the protostellar magnetic field is able to disrupt the disk as soon as the magnetic pressure $P_{\text{mag}}=B_z^2/(8\pi )$ becomes comparable to P_gas, which for a dipole field approximation occurs at

$$\begin{array}{c}\hfill r_{\mathrm{in}}(P_{\mathrm{gas}}>P_{\mathrm{ram}})=\frac{{B_{\star }}^{1/3}{R}_{\star }}{{P_{\mathrm{gas}}}^{1/6}{(8\pi )}^{1/6}}·\end{array}$$(48)

In this work, we find two scenarios in thermally unstable disks where the gas pressure P_gas exceeds the accretion ram pressure P_ram. The first is after an outburst, when the accretion rate decreases until the disk evolves towards the next burst (cf. Sect. 5). The second is towards the end of a disk’s lifetime when the accretion rate decreases, which yields a lower ram pressure P_ram (see Sect. 5.2).

5. Results and discussion

The goal of our simulations is to show the influence of a protostellar magnetic field on the bursting behavior of MRI- and TI-unstable disks. Therefore, in Sect. 5.1 the bursting behavior and the long-term evolution of such disks is investigated. In Sect. 5.2 we perform long-term runs with different magnetic field strengths B_⋆ to investigate the effect of magnetic torques on the inner disk dynamics and its consequences for the long-term behavior. Additionally, the effect of different MRI activation temperatures is discussed in Sect. 5.3.

5.1. Episodic accretion and long-term evolution of magnetically truncated low-mass protoplanetary disks

Armitage et al. (2001) and Zhu et al. (2010) argue that a protoplanetary disk cannot sustain a steady mass transport rate in the radial direction from r ≈ 100 AU to its inner edge at a few stellar radii, because GI in the outer disk feeds too much mass to the inner disk where it piles up. This leads to viscous heating in this region of enhanced density and eventually to the activation of the MRI. Enhanced radial mass transport Ṁ in the MRI-active disk yields more viscous heating and hence facilitates triggering of the thermal instability (e.g., Bell & Lin 1994).

In Table 1 the parameters for our fiducial model fid are stated. Ṁ_init corresponds to the constant mass transport rate throughout the disk, which is obtained for the stationary initial model constructed as described in Sect. 3.5. The protostellar parameters (mass M_⋆, radius R_⋆, magnetic dipole field strength B_⋆, and the corotation radius r_cor), together with the disk parameters (α-viscosity parameter due to MRI α_MRI, the ratio α_MRI/α_base, and the thermal instability activation temperature T_active) and the opacity prescription (cf. Fig. 1), determine the initial disk mass M_disk, init for the stationary initial model. We choose our parameters such that M_disk, init ≈ 0.01 M_⊙ to ensure a gravitationally stable disk (see Sects. 2.1 and 2.4). For all simulation runs in Table 1, no external mass reservoir is assumed to resemble the configuration of a late class II star–disk system, and therefore Ṁ(r_out) = 0.

The column density Σ, midplane gas temperature T_gas, scale height H_p, the radial mass transport rate Ṁ, the relative deviation from the Keplerian velocity u_K, and the midplane gas pressure P_gas, 0 of our fiducial model fid are shown in Fig. 7 during onset of a TI-triggered burst. Zhu et al. (2010) observed short phases during outbursts in their simulations, where the thermal instability could not be sustained (they termed this phenomenon drop-outs), whose cause they interpreted to be the lack of radial advection in their simulations. We include radial advection (cf. Sect. 2.1) in our models and do not observe such drop-outs, which confirms their interpretation. Shortly before t_onset, 1, a TI is triggered at r ≈ 15 R_⊙, which can be seen by a steep rise in T_gas, 0 (panel b of Fig. 7). Simultaneously, the TI strongly increases local viscosity, to which the disk reacts by redistributing material to both the inner and outer disk annuli (e.g., Bell & Lin 1994). This increases Σ at both neighboring disk annuli above the threshold at which TI is also triggered. Hence, two ionization fronts are starting to travel inwards and outwards, which can be seen most clearly in panels a, b, and e of Fig. 7 for different points in time from t_onset, 1 to t_onset, 4. Both ionization fronts heat up the disk (see Fig. 7b) which leads to an increased scale height H_p (Fig. 7c). This in turn yields a region where the disk is shadowed from protostellar irradiation. Our models incorporate geometric shadowing caused by local elevations in H_p (cf. Sect. 2.3), which can be seen in Fig. 7b as a drop in surface temperature T_surf just behind (from the perspective of the protostar) the outward-traveling ionization front. As soon as the inward-bound wave reaches the rim of the inner disk, heated (and hence more viscous) material can no longer be redistributed further inwards, and the disk starts to transition into the outburst phase. The gas temperature at the inner boundary T_gas, 0(r_in) rises very quickly and is accompanied by a sharp increase in viscosity and hence in Ṁ_⋆ (cf. panel g of Fig. 7).

Fig. 7. Various disk quantities for certain times tonset, 1 − tonset, 4 describing an outburst onset phase. Panels a–c: radial surface density structure Σ and the pressure scale-height Hp during the onset of an outburst, respectively. Panel b: midplane gas temperature T0 (solid) and the surface temperature Tsurf (dashed). The dotted line marks the activation temperature Tactive. Panel d: radial mass transport rate Ṁ, where the dotted line corresponds to Ṁ = 0. The deviation from a Keplerian accretion disk is shown in panel e, whereas panel f: gas pressure at the midplane Pgas, 0. Finally, panel g maps the various snapshots in time to the stages during the onset of an outburst. The dashed black line in panels a–f denote the stationary initial model as described in Sect. 3.5.

The deviation from the diffusion approach investigated in Sect. 4.1 can be seen clearly in panels e and f of Fig. 7. The relative difference of u_φ compared to the Keplerian velocity u_K peaks at ∼0.02, which is an order of magnitude larger than the predicted u_φ ≈ 0.996 u_K for a viscous protoplanetary accretion disk (see Armitage 2010). Additionally, the rotational velocity u_φ not only becomes subKeplerian but also superKeplerian (see panel e in Fig. 7). This feature is strongest along the propagating ionization fronts, as there the surface density has localized peaks, which in turn leads to an altered (with respect to a quasiKeplerian accretion disk) viscous net shear on the disk annuli at those radii. The relative influence of gas pressure P_gas on a certain annulus compared to the viscous force is shown in Fig. 8. The region at around 10 solar radii in Fig. 8 corresponds to the location of the inward-bound ionization front, whereas the outer ionization front is located at approximately 24 R_⊙. The magnetic braking of the inner disk inside r_cor leads to an increased radial velocity u_r and hence to a drop in Σ towards r_in, which leads to the pressure force changing sign at around 22 R_⊙. Furthermore, it can be seen in Fig. 8 that the gas pressure alters the disk dynamics close to the inner disk radius, which is because that is where the pressure gradient becomes steepest (apart from local disturbances like ionization fronts). The influence of the pressure gradient even exceeds the viscous contributions, which is, once more, an indication that the full set of hydrodynamic equations are required to properly model the inner regions of a disk. The sharp decrease in Ṁ_⋆ between t_onset, 2 and t_onset, 3 in Fig. 7 occurs because of the small-scale structure of the inbound wave when it reaches r_in (for a thorough analysis of the thermal instability see Bell & Lin 1994).

Fig. 8. Net contributions of pressure Δfp and viscous force Δfν on an disk annulus are shown for the inner disk. The forces are in units of the net gravitational force Δfg at its corresponding radius. The dashed lines outline the various force contributions for the stationary initial model. The dotted horizontal line shall help to visualize where the force contributions change sign.

Figure 9 shows the whole disk during a TI-induced outburst. The outward-bound ionization front as described in detail by Bell & Lin (1994) stalls at t ≈ t_burst, 4 at a radius where the swept-up material is no longer capable of increasing Σ sufficiently for the disk annulus to become thermally unstable (cf. b in Fig. 9). The increased accretion rate yields a higher accretion luminosity, which in turn causes a higher temperature at the disk surface and a puffing-up of the irradiated outer disk (panel c of Fig. 9). The detailed effects of the increased scale height, relating to shadowing of the disk behind this bump, will be studied in another paper of this series. The ongoing depletion of the inner disk inside the ionization front due to strongly increased viscosity eventually leads to a drop in gas temperature T_gas, 0. At approximately t_burst, 3, hydrogen is able to recombine again and a cooling wave starts to propagate inwards (radially outer green vertical dashed line and arrow in panel b of Fig. 9). The cooler gas then leads to a drop in optical depth (cf. Fig. 1), which results in termination of the TI. This inward-bound perturbation becomes clearly visible at t_burst, 4 and t_burst, 5 in panels a to c of Fig. 9, where the position and direction of the cooling front is marked with a vertical dashed line and an arrow, respectively. As this wave propagates towards the inner disk radius r_in, the lower temperature leads to a decrease in viscosity and hence the outburst enters the decaying phase and the disk inside the ionization front empties on a viscous timescale t_visc, which for our model fid is t_visc ≲ 100 years in the MRI-active region.

Fig. 9. Panels a–c, as well as (e), are the same disk properties as in Fig. 7, but for five points in time tburst, 1 to tburst, 5 during an outburst. Panel d: same quantity as (e) in Fig. 7. The dashed black line in panels a–d denote the stationary initial model as described in Sect. 3.5. The radial range in this plot is chosen from rin to 10 AU. The color-coded, dashed vertical lines and corresponding arrows mark the position of the cooling waves and their propagation direction, respectively.

In panel d the influence of the protostellar magnetic field (see also Sect. 5.2) is revealed in the form of a considerably sub-Keplerian rotational velocity u_φ towards the inner radius. This is due the protostellar magnetic torque acting on the disk. As a further consequence, the radial velocity of the disk starts to increase inside the corotation radius r_cor. However, before the radial velocity u_r is even close to the free-fall velocity, the protostellar magnetic field dominates the disk dynamics (e.g., Bessolaz et al. 2008) and the accretion flow transitions into accreting funnels. This transition zone at the inner rim leads to a drop in surface density Σ and as a consequence yields a lower gas temperature T_gas, 0. Additionally, because of the decreasing accretion rate Ṁ_⋆, the inner radius is moving outwards (cf. panels a to d in Fig. 9), which in combination with the lower temperature leads to an additional cooling wave starting at the inner edge and propagating outwards (color-coded vertical dashed lines and corresponding radially outward-oriented arrows close to the inner rim in panel b of Fig. 9). At t_burst, 5, the two waves have almost met each other (cf. color-coded vertical lines corresponding to t_burst, 5 in panel b of Fig. 9), resulting in the turnoff of the TI almost everywhere in the inner region, eventually marking the end of the outburst. Hence, compared to models without the influence of a stellar magnetic field, an additional cooling wave traveling outwards has an influence on the time needed for the inner disk to become thermally stable again. This because the outer cooling front does not travel inwards until it reaches the inner disk rim but is met by the outward-bound cooling wave earlier. Material is then fed to the inner regions from the outer regions and accumulates until a disk annulus becomes thermally unstable again, therefore completing a burst cycle.

Such burst cycles occur as long as the inner disk can be fed with material from the outer regions and consequently is capable of heating up the inner regions sufficiently through viscous dissipation. In Fig. 10 this repeated episodic accretion can be seen in panel a. This process is repeated 49 times before the disk mass can no longer provide a sufficiently strong accretion flow of gas to heat the inner disk above the MRI activation temperature T_active once more. For our fiducial model fid, the disk becomes quiescent at ∼2.35 × 10⁵ years and below a disk mass of ∼7 × 10⁻³ M_⊙. Panel b of Fig. 10 shows the variation of the inner radius r_in as it reacts to a changing accretion rate. After the bursts have ceased and the disk has become quiescent, the accretion rate decreases, as the disk continually loses mass over its lifetime (cf. panel c of Fig. 10). The constantly decreasing accretion rate also has an impact on the location of the inner radius. The magnetic field is able to disrupt the disk farther out, as Ṁ_⋆ becomes weaker. The apparent push of r_in outside the initial corotation radius r_cor, init can be explained with a super-Keplerian rotation velocity caused by the positive pressure gradient in the innermost disk in the vicinity of r_cor, init (cf. Fig. 8) and therefore a slight change in r_cor towards larger radii.

Fig. 10. Long-term evolution for our fiducial model fid (see Table 1). Panels a and b: accretion rate Min at the inner boundary and its corresponding movement of the inner radius rin. The evolving total disk mass Mdisk is presented in panel c. The inset in panel c shows the step-like decrease in disk mass Mdisk caused by outbursts.

5.2. Influence of the stellar magnetic field

In Sect. 5.1 we analyze and describe our fiducial model fid. It can be clearly seen that a protostellar magnetic field is influencing the disk close to the inner boundary r_in. In agreement with Vorobyov et al. (2019), we argue that the inner boundary has an effect on the global disk structure and therefore on the long-term evolution of a protoplanetary disk. We emphasize in Sect. 3.1 that a 1D approach outside of r ≈ 0.5 AU has certain limitations in describing the disk structure appropriately. However, our focus in this section is to investigate how different magnetic field strengths (see Table 1) are effecting the bursting behavior of the disk. The MRI and TI ignition point lies well below 0.5 AU, which is an essentially axisymmetric region of the disk (e.g., Lesur 2020), and therefore the assumptions of our model are valid approximations. Additionally, we restrict our studies in this work to low-mass disks where no gravitational instabilities are expected to break their axisymmetry.

In Fig. 11 the initial magnetic field strengths for different stationary models (cf. Table 1) are shown. We can determine two effects of the stellar magnetic field. First, the inner boundary of the disk is pushed outward for a higher field strength due to the stronger magnetic pressure (see Eq. (47)). Additionally, a stronger field also has a stronger propelling effect on the disk gas just outside of the corotation radius r_cor (cf. panel a of Fig. 11). Material in this area starts to get pushed outwards, and the disk starts to pile up material because disk gas is still accreting from the outer disk inwards (see panel c of Fig. 11). This leads to increased viscous heating and therefore to a higher midplane gas temperature T_gas, 0 (see panel b of Fig. 11). If the magnetic field is chosen to be stronger than a certain threshold field strength (the exact value depends on the disk and stellar parameters), then the aforementioned pile-up of material is sufficient to trigger the MRI. For the disk parameters chosen in this work, this limiting value is B_⋆, up ≈ 4.5 kG. Consequently, the higher T_gas, 0 is sufficient for deeper disk layers to be become MRI-active, which in turn leads to even higher viscosity and temperature (cf. panel b of Fig. 11). Eventually, this triggers TI at some ignition radius, which then starts an outburst cycle.

Fig. 11. Comparison of (a) the toroidal protostellar magnetic field component Bφ as stated in Eq. (8), (b) the midplane gas temperature Tgas, 0, and (c) the surface density Σ for the different stationary initial models in Table 1 of varying protostellar magnetic field strengths, ranging from B⋆ = 1.5 kG to 5 kG. The MRI active zone starts at little over 1400 K for our choice of Tactive = 1530 K and Twidth = 50 K. The vertical dashed lines denote the position of the stellar radius R⋆ and the initial corotation radius rcor, init.

Also notable is the effect of a relatively weak magnetic field on the disk. A lower magnetic field strength yields an inner radius r_in closer to the protostar. Consequently, the disk is heated up to higher temperatures and the MRI/TI can be triggered. As a result, for all stellar magnetic fields B_⋆ lower than a lower limit field strength B_⋆, low, a disk also undergoes episodic accretion. For the parameters in Table 1, this lower limit is at B_⋆, low ≈ 2 kG, which we use for our fiducial model fid (cf. panel b of Fig. 11).

We also conduct long-term runs for all models in Table 1. In Fig. 12, the accretion rate on the star Ṁ_⋆ and the movement of the inner radius r_in are plotted in panels a and b, respectively. The models 1p5kG (blue solid) and fid (orange solid) become thermally unstable due to the proximity of the inner radius r_in to the star, which can be seen in panel b of Fig. 12. The burst frequency (cf. panel a of Fig. 12) depends on how fast the gas temperature T_gas, 0 exceeds the activation temperature T_active. For disks with stellar magnetic field strengths B_⋆ < B_⋆, low the burst frequency increases for weaker magnetic fields B_⋆, because the weaker the stellar field B_⋆, the further inwards the disk stretches. This leads to even higher gas temperatures and therefore to a disk prone to TI. Equally importantly, for B_⋆ > B_⋆, up, the stronger B_⋆, the more effectively the propelling effect can fling material outwards, and hence the more material piles up outside of r_cor. The gas temperature then also increases with higher B_⋆ and TI is triggered faster. Models 1p5kG and fid satisfy the constraint B_⋆ < B_⋆, low and 5p0kG satisfies B_⋆ > B_⋆, up. For intermediate magnetic fields of B_⋆, low < B_⋆ < B_⋆, up, no outbursting behavior develops (models 3p0kG and 4p0kG). We note that the actual threshold values for B_⋆ depend on the choice of stellar and disk parameters and have to be determined anew for every different parameter set (cf. Fig. 16). For weaker fields, the inner disk radius r_in is pushed very close to the stellar radius (cf. Eq. (47) and panel b of Fig. 12); in the case of our model 1p5kG, the inner radius r_in ≈ 1.5 R_⋆.

Fig. 12. Comparison of the models in Table 1. Panel a: accretion rate Ṁ⋆ and (b) the corresponding inner radius rin, respectively. The colors correspond to the same models as in Fig. 11. The horizontal dashed lines denote the position of the stellar radius R⋆ and the initial corotation radius rcor, init.

In Fig. 13 the full disk lifetime is plotted. In panels a and b it can be seen that the transition of the disk into its quiescent state occurs later for models 1p5kG and 5p0kG, which is due to their higher maximum temperature T_gas, 0 (compared to fid) which can also be seen in panel b of Fig. 11. The physical explanation for the different rest disk masses in panel c of Fig. 13 is that below a certain disk mass threshold, the protostellar field is effectively preventing further accretion onto the protostar (propeller regime) and hence the disk mass stabilizes at a certain value. The strength of the magnetic field then determines the lower disk mass M_disk, below which further accretion is quenched by B_φ.

Fig. 13. Same models as in Fig. 12, but for the full disk lifetime. The onset of the first burst of every model has been aligned in time to provide a better comparability regarding burst frequency and burst strength. Panel b: the dashed lines denote the initial corotation radius Rcor, init and the protostellar radius R⋆.

Towards the end of the disk lifetime, the disk mass and consequently the accretion rate Ṁ_⋆ drop, resulting in r_in shifting towards larger radii (cf. Fig. 13 and Eq. (47)). With decreasing Ṁ_⋆ the gas pressure P_gas eventually exceeds the ram pressure P_ram and controls the movement of the inner disk radius r_in. For all models, this transition occurs inside the corotation radius r_cor. For a sufficiently strong stellar magnetic field, the propelling effect just outside r_cor piles up material (cf. Fig. 11) and increases P_gas until it exceeds P_ram, thus keeping r_in inside r_cor. We want to note here that current 2D simulations of the propelling regime show a more complex behavior (e.g., Ustyugova et al. 2006; Romanova et al. 2009, 2018). A specific part of the piled up matter is accreted onto the star while another is ejected in winds. Inclusion of disk winds in our model, for example, which are shown to play an important role for the inner boundary (e.g., Königl et al. 2011), would reduce pile up of matter (see panel a in Fig. 11) and consequently increase B_⋆, up. Thus, the results presented here have to be considered as a limiting case.

In addition to the long-term evolution of the disk for different values of the stellar magnetic field, we compare the duration of bursts and the accreted disk mass onto the star ΔM_⋆ for each individual outburst (cf. Fig. 14). For a given stellar magnetic field strength, the burst duration as well as ΔM_⋆ remain approximately constant until the disk has depleted to a degree where the stellar disk can no longer feed the inner region with sufficient mass to become MRI- or TI-unstable. The small variability of the duration and the accreted mass of consecutive bursts is due to the fact that the ignition radius for the MRI/TI remains effectively unaltered for the disk during its bursting phase, which then yields a very similar maximum radius for the outward-bound ionization front and consequently a similar burst duration and ΔM_⋆. Relaxing the assumption of axisymmetry probably adds more variability to the results (e.g., Vorobyov et al. 2020) and shall be reviewed in further studies. However, for an increasing stellar magnetic field strength, the burst duration and ΔM_⋆ increase. Because of the stronger magnetic pressure with an increasing stellar magnetic field, the inner disk radius r_in as well as the MRI ignition point are pushed outwards (cf. panel b of Fig. 11). Additionally, more mass is piled up directly outside the MRI ignition point (cf. panel c of Fig. 11), which also increases the disk temperature in this region. As a result, the outward-bound ionization front reaches further out, which increases the burst duration and the piled up mass outside the MRI ignition point increases ΔM_⋆.

Fig. 14. Time evolution of the outbursting behavior of the disk, shown for all consecutive bursts during disk evolution. Panel a: mass accretion onto the protostar during an outburst ΔM⋆, whereas panel b: burst duration for a certain burst.

Figure 15 specifically depicts the thermally stable zone for B_⋆ between roughly 2.5 kG and 4.5 kG. The final disk mass towards the end of a disk’s lifetime seems to behave like a power-law (cf. panel c of Fig. 15). A possible explanation for this is that an almost dissolved disk with little mass cannot counteract the propelling forces of the stellar magnetic field and hence cannot accrete onto the star. The stellar field is modeled as a dipole and the resulting centrifugal force would also yield a power law; nevertheless a more thorough analysis is planned in further studies to investigate this behavior in more detail. Panel b of Fig. 15 depicts the point in time where the disk ceases to burst (burst offset), whereas panel a shows the number of bursts for a certain model. Figure 15 indicates that stellar magnetic fields, which are weaker than the lower limit or stronger than the upper limit of B_⋆, will lead to more bursts, depending on the size of the deviations from those limits.

Fig. 15. Comparison of (a) the number of bursts, (b) the burst offset time of the disk, and (c) the final disk mass for all models in Table 1 and additionally for models with B⋆ = 2.5 kG, 3.5 kG and 4.5 kG. The color coding is the same as for Fig. 11, whereas the black dashed line linearly connects the data points and has the sole purpose of emphasizing the burst-free region in panels a and b as well as indicating the power-law for the final mass in (c).

5.3. Variation of activation temperature and mass transport rate

To test the robustness of the results of Sect. 5.2, we varied the MRI activation temperature T_active together with the initial mass transport rate Ṁ_init. For all previous simulations, T_active is fixed to 1530 K (see Table 1). A different choice of T_active will affect the radius at which the deep layer of the disk becomes MRI-active and TI is triggered. Choosing a lower activation temperature will result in more models becoming thermally unstable. The value range between B_⋆, low and B_⋆, up becomes smaller. A higher activation temperature T_active, on the other hand, will increase this region as there will be fewer models in which T_gas can exceed T_active. We can reproduce our main results for different values of T_active by adapting the mass transport rate Ṁ throughout the disk. Choosing T_active = 1300 K and Ṁ_init of 0.97 × 10⁻⁸ M_⊙/yr⁻¹, the range of the stellar magnetic field B_⋆ in which outbursts are suppressed extends from 2.5 to 3.5 kG (see Fig. 16). For choosing T_active = 1400 K and Ṁ_init of 1.0 × 10⁻⁸ M_⊙/yr⁻¹, outbursts are suppressed between 2.0 and 4.5 kG. The adapted values of Ṁ_init, which correspond to the quiescent accretion rate between two outbursts, are all in agreement with Armitage et al. (2001). We would like to emphasize that the exact interpretation of these results requires a more detailed parameter study and is left for future work.

Fig. 16. Bursting behavior depending on the stellar magnetic field B⋆ for different Tactive. A distinction is made between models that produce outbursts (B) and models in which outbursts are suppressed (0). Tactive is set to 1300 K and 1400 K in panels a and b, respectively. The accretion rate Ṁinit is given in units of 10−8 M⊙/yr−1.

6. Conclusion

In this study we conduct fully implicit 1+1D hydrodynamic simulations of protoplanetary disks with a focus on the inner regions. The TAPIR code (see Ragossnig et al. 2020) is capable of solving the full set of hydrodynamic equations, which enables us to include additional physics, such as for example the effects of local pressure gradients, protostellar magnetic fields, or features in the angular velocity.

Previous studies (e.g., Bell & Lin 1994; Armitage et al. 2001; Zhu et al. 2009a) wherein 1D long-term simulations were also performed used the viscous diffusion equation for the temporal evolution of the disk. We demonstrate in Sect. 4 that this approach shows significant deviations –at least locally close to the star– from a full solution of the hydrodynamic equations. Furthermore, we also adopt a description of the magnetic truncation radius based on the equality of magnetic pressure and ram pressure (see e.g., Hartmann et al. 2016). Complementary to that approach, we incorporate the gas pressure influence in our models, covering P_gas > P_ram.

In the first paper of this series, we study the effect of protostellar magnetic torques on the long-term behavior of protoplanetary disks. These torques slow down or accelerate the disk depending on the radial location compared to the corotation radius r_cor, and thus necessitate a self-consistent treatment of the momentum equation in the angular direction. Diffusion approaches are not well suited to modeling an inner disk threaded by a stellar magnetic field.

Layered-viscosity models as proposed by Armitage et al. (2001) can, for certain parameters, become thermally unstable and repeatedly experience FU Ori-like episodic accretion events (e.g., Bell & Lin 1994). Interestingly, the burst strength and burst frequency are strongly dependent on the stellar magnetic field strength B_⋆.

On the one hand, a stronger B_⋆ corresponds to a larger inner disk radius r_in, which leads to less radiative heating of the surface of the inner disk and consequently to a cooler disk. On the other hand, a stronger toroidal magnetic field component B_φ also acts as a propeller outside of the corotation radius r_cor, effectively flinging material outwards. This pushed out disk gas combines with the inward-directed accretion flow, resulting in a larger radially localized surface density Σ. This leads to a higher gas temperature T_gas, which counteracts the cooling effect of a farther out r_in. We find that there exists a limit value for B_⋆, above which a disk gets sufficiently heated to become thermally unstable.

A weaker B_⋆ leads to a small inner radius r_in, and therefore the inner disk gets hotter due to the proximity to the protostar. We find that if the protostellar magnetic field B_⋆ is weaker than a certain upper limit, then the inner disk becomes hot enough to also become thermally unstable.

Summarizing these results, a region of protostellar magnetic field of intermediate strength exists where (for a certain set of disk and star parameters) no bursts occur. Below and above this region the models can become thermally unstable. These results are later reproduced with different value sets of the MRI activation temperature T_active and mass transport rate Ṁ. Additionally, the burst strength and burst frequency are also altered by a varying B_⋆, which changes the duration of the thermally unstable phase of such disks. We find that the burst duration and mass accretion rate caused by a single outburst only vary slightly for consecutive bursts, which is due to the almost constant radial position of the MRI/TI ignition point during the bursting phase of a disk. Moreover, for stronger stellar magnetic fields, the burst duration and accretion mass during an outburst increase because the MRI ignition point gets pushed radially outwards. We conducted long-term simulations to study the influence of the protostellar field and found that disks tend to have more frequent but weaker bursts if B_⋆ is either far below or far above the aforementioned intermediate zone. Our simulations show that for weaker stellar magnetic fields, the inner disk radius r_in is pushed very close to the stellar radius R_⋆. For even weaker fields or stronger bursts the magnetic field could be squashed to the stellar surface and could therefore suppress the accretion funnels as well as the resulting accretion shock (e.g., Hartmann et al. 2016) and thus provide an explanation for the lack of UV excess in FU Ori-like star–disk systems (e.g., Bell & Lin 1994).

In this work, we show that a detailed model of the inner disk regions is essential for the simulations of protoplanetary disks, especially during episodic outbursts. We agree in this regard with Vorobyov et al. (2019). for example, and are confident that the TAPIR code (Ragossnig et al. 2020) is well suited to representing the inner disk region in global 2D long-term simulations.

6.1. Limitations of this work

While the TAPIR-code is capable of treating the inner disk self-consistently, there are also limitations to this approach. Our 1D model is mainly suited to simulation of the inner disk. This is because in the inner regions, viscous shear becomes stronger and therefore angular perturbations cannot be sustained, even for short times (for a conclusive discussion see e.g., Lesur 2020). Therefore the assumption of axisymmetry is a valid approximation in the inner regions.

Furthermore, because of the assumption of vertical thermal equilibrium, radiative transport in the z-direction is only included approximately. Photoevaporation of the inner disk definitely has an effect on the inner radius of the disk, as well as on the star–disk-interaction as gas is removed, especially from the inner part of the disk. However, these effects will not effect the qualitative statements of this work and will be investigated in upcoming studies. Furthermore, we use a simple dust model with a constant dust-to-gas ratio. The dust is included in the opacity calculation, but the back-reaction onto the gas or dust dynamics is not considered yet. This approximation has only a minor effect on the main findings from this study. However, a more realistic time-dependent dust model would increase the accuracy of our code.

We also neglect some details of the star–disk interaction. In our model, B_⋆ is prescribed as a dipole field, which is an approximation, especially in the case of very high accretion rates. Nevertheless, for an qualitative study of the effects of magnetic torque on the disk evolution, the details of the field topology are of minor order. Additionally, the transfer of angular momentum from the star to the disk and vice versa has not been modeled in this work, as we focused on the effect of an existing magnetic field on disk evolution. We agree that this effect should be investigated in further studies, because a time-dependent stellar angular velocity also results in the corotation radius r_cor moving radially inwards and outwards, which has a potential effect on the outbursting behavior of disks.

6.2. Further work

In upcoming studies of this series, we will study the effects of magnetic fields and stellar radiation on the pressure scale height, which leads to a shadowed disk region, where less stellar radiation can heat the disk surface. Another study will involve a model for the star–disk interaction and the inclusion of higher-order magnetic field moments like quadrupole and octopole contributions to the stellar magnetic field topology. We will use the altered magnetic field structure to study the effect on the outbursting behavior of disks.

Equally interesting is the inclusion of large-scale disk fields. Those fields have been investigated in 1D models (e.g., Guilet & Ogilvie 2012, 2014). The incorporation of such magnetic field models in the TAPIR code enables us to self-consistently solve the inner disk and the disk field topology simultaneously. This further allows estimates of how much mass and angular momentum can be transported away by magneto-centrifugally driven outflows during for example an FU Ori-like outburst. Also important is the inclusion of photo-evaporative winds, especially during an outburst, because this is when the accretion luminosity is highly elevated. The effects of outflows due to photo-evaporation coupling with a large-scale disk magnetic field are equally interesting, as then not only mass is transported away but also angular momentum, which further modifies the disk dynamics.

Acknowledgments

E. I. V. acknowledge support of Ministry of Science and Higher Education of the Russian Federation under the grant 075-15-2020-780 (N13.1902.21.0039; Sect. 6).

References

Appendix A: Implicit time-step

Implicit numerical models are not bound to the CFL condition and can benefit from larger time-steps and consequently shorter simulation times compared to explicit models. To quantify this difference, we compare the maximum possible time-step– Δt ≤ Δr/u_i –allowed by the CFL condition; with the radial grid size Δr and a velocity u_i with which information propagates. In an accretion disk, there are several velocities that transport information. For this order-of-magnitude approximation, we want to compare the effects of the speed of sound c_S with those of the radial velocity u_r. Additionally, the time-step limitation of diffusion in the disk– Δt ≤ Δr²/D –is taken into account. Here, D is the diffusion coefficient and is connected to the viscosity ν via the Schmidt number Sc = ν/D. Following Armitage (2011), Sc is of the order of unity and Sc = 1 is used. In Fig. A.1, the time-step used in the implicit TAPIR code is compared to the theoretical time-step limitations imposed by the CFL condition for an explicit method. During an outburst, the implicit time-step exceeds the explicit time-step by up to two orders of magnitude and during the quiescent phase, by up to six orders of magnitude. Our model also incorporates an adaptive time-step control, which increases or decreases Δt until the Newton-Raphson iteration (cf. Sect. 3.2) converges towards a solution. This results in shorter time-steps during bursts and longer time-steps in quiet phases, because in the latter case the disk quantities change slowly over time and a solution can still be found for longer time-steps. We can conclude that an explicit model with the same radial resolution would increase the computational time from days to months or even years, which again justifies this implicit approach.

Fig. A.1. Comparison of the implicit time-step Δtimp to time-step limitations imposed by the CFL condition with respect to the speed of sound ΔtcS, the radial velocity Δtur, and the diffusion Δtdiff for our fiducial model; (a) during an outburst and (b) during the quiescent phase between two outbursts. The decrease in Δtimp at the right side denotes the beginning of the next outburst.

All Tables

All Figures

	Fig. 1. Rosseland-mean opacities κR for various midplane gas densities ρ0.
In the text

	Fig. 2. Snapshots in time, which show the convergence towards a stationary initial model for (a) the surface density Σ(r), (b) the radial velocity ur(r), and (c) the mass transport rate Ṁ. The black dashed line shows the initial Keplerian disk, and the black solid line denotes the final stationary model.
In the text

	Fig. 3. a: Ω1 waves for four different, arbitrarily chosen times t0 to t3 during the onset of a FU Ori-like outburst. b: correction factors as defined in Eq. (41) for the same times as in (a). The dashed-dotted lines correspond to the first term in Eq. (41), whereas the full lines show the total value of fcorr.
In the text

	Fig. 4. Time relative to the burst duration in the fiducial model Δtburst, in which a region in the inner 0.5 AU exceeds a certain fcorr threshold (0.05, 0.5 and 1.0). The radial range is divided in 200 equally spaced logarithmic bins.
In the text

	Fig. 5. a: midplane gas pressure Pgas, 0 between 0.04 AU and 2 AU for two different times during an FU Ori like outburst. b: midplane pressure gradient force −∇Pgas, 0 relative to the gravitational force is plotted. The radially localized changes in sign occur at the positions of the ionization fronts running inwards and outwards. The deviation from the Keplerian angular velocity is shown in (c).
In the text

	Fig. 6. a: accretion rate for a full TI-triggered outburst for our fiducial model fid (cf. Table 1). b: variation of the inner boundary during the burst is shown in (a). Here, t0 corresponds to the onset of the first burst. The two dashed lines denote the protostellar radius R⋆ and the corotation radius rcor.
In the text

Fig. 7. Various disk quantities for certain times tonset, 1 − tonset, 4 describing an outburst onset phase. Panels a–c: radial surface density structure Σ and the pressure scale-height Hp during the onset of an outburst, respectively. Panel b: midplane gas temperature T0 (solid) and the surface temperature Tsurf (dashed). The dotted line marks the activation temperature Tactive. Panel d: radial mass transport rate Ṁ, where the dotted line corresponds to Ṁ = 0. The deviation from a Keplerian accretion disk is shown in panel e, whereas panel f: gas pressure at the midplane Pgas, 0. Finally, panel g maps the various snapshots in time to the stages during the onset of an outburst. The dashed black line in panels a–f denote the stationary initial model as described in Sect. 3.5.

In the text

	Fig. 8. Net contributions of pressure Δfp and viscous force Δfν on an disk annulus are shown for the inner disk. The forces are in units of the net gravitational force Δfg at its corresponding radius. The dashed lines outline the various force contributions for the stationary initial model. The dotted horizontal line shall help to visualize where the force contributions change sign.
In the text

Fig. 9. Panels a–c, as well as (e), are the same disk properties as in Fig. 7, but for five points in time tburst, 1 to tburst, 5 during an outburst. Panel d: same quantity as (e) in Fig. 7. The dashed black line in panels a–d denote the stationary initial model as described in Sect. 3.5. The radial range in this plot is chosen from rin to 10 AU. The color-coded, dashed vertical lines and corresponding arrows mark the position of the cooling waves and their propagation direction, respectively.

In the text

	Fig. 10. Long-term evolution for our fiducial model fid (see Table 1). Panels a and b: accretion rate Min at the inner boundary and its corresponding movement of the inner radius rin. The evolving total disk mass Mdisk is presented in panel c. The inset in panel c shows the step-like decrease in disk mass Mdisk caused by outbursts.
In the text

Fig. 11. Comparison of (a) the toroidal protostellar magnetic field component Bφ as stated in Eq. (8), (b) the midplane gas temperature Tgas, 0, and (c) the surface density Σ for the different stationary initial models in Table 1 of varying protostellar magnetic field strengths, ranging from B⋆ = 1.5 kG to 5 kG. The MRI active zone starts at little over 1400 K for our choice of Tactive = 1530 K and Twidth = 50 K. The vertical dashed lines denote the position of the stellar radius R⋆ and the initial corotation radius rcor, init.

In the text

	Fig. 12. Comparison of the models in Table 1. Panel a: accretion rate Ṁ⋆ and (b) the corresponding inner radius rin, respectively. The colors correspond to the same models as in Fig. 11. The horizontal dashed lines denote the position of the stellar radius R⋆ and the initial corotation radius rcor, init.
In the text

	Fig. 13. Same models as in Fig. 12, but for the full disk lifetime. The onset of the first burst of every model has been aligned in time to provide a better comparability regarding burst frequency and burst strength. Panel b: the dashed lines denote the initial corotation radius Rcor, init and the protostellar radius R⋆.
In the text

	Fig. 14. Time evolution of the outbursting behavior of the disk, shown for all consecutive bursts during disk evolution. Panel a: mass accretion onto the protostar during an outburst ΔM⋆, whereas panel b: burst duration for a certain burst.
In the text

Fig. 15. Comparison of (a) the number of bursts, (b) the burst offset time of the disk, and (c) the final disk mass for all models in Table 1 and additionally for models with B⋆ = 2.5 kG, 3.5 kG and 4.5 kG. The color coding is the same as for Fig. 11, whereas the black dashed line linearly connects the data points and has the sole purpose of emphasizing the burst-free region in panels a and b as well as indicating the power-law for the final mass in (c).

In the text

	Fig. 16. Bursting behavior depending on the stellar magnetic field B⋆ for different Tactive. A distinction is made between models that produce outbursts (B) and models in which outbursts are suppressed (0). Tactive is set to 1300 K and 1400 K in panels a and b, respectively. The accretion rate Ṁinit is given in units of 10−8 M⊙/yr−1.
In the text

	Fig. A.1. Comparison of the implicit time-step Δtimp to time-step limitations imposed by the CFL condition with respect to the speed of sound ΔtcS, the radial velocity Δtur, and the diffusion Δtdiff for our fiducial model; (a) during an outburst and (b) during the quiescent phase between two outbursts. The decrease in Δtimp at the right side denotes the beginning of the next outburst.
In the text