Eruptive Behavior of Magnetically Layered Protoplanetary Disks in Low-metallicity Environments

Kundan Kadam; Eduard Vorobyov; Ágnes Kóspál

doi:10.3847/1538-4357/abdab3

1. Introduction

The general picture of low-mass star formation suggests that, as a result of the conservation of angular momentum, a pre-main-sequence star inevitably forms a surrounding flattened accretion disk. Planetary systems are born within the environment of such protostellar (younger and relatively massive) and protoplanetary disks (PPDs). It is widely accepted that the disk accretion in such systems, also called young stellar objects (YSOs), is not steady but time-dependent, with low mass accretion rates punctuated by episodes of powerful eruptions (e.g., see Audard et al. 2014, and references therein). This picture is also supported by the surveys of protostars, which suggest that the luminosity of YSOs is consistently about an order of magnitude lower than that expected from a steady accretion of mass (Eisner et al. 2005; Dunham et al. 2010). Sudden and luminous accretion events known as FUor- (∼100 L_⊙, ∼100 yr, more resembling embedded class I YSOs) and EXor- (∼10 L_⊙, ∼1 yr, spanning both class I and II) type eruptions have also been directly observed in YSOs (Hartmann & Kenyon 1996; Herbig 2008; Audard et al. 2014, and references therein). On average, about 10% of the final stellar mass of a low-mass star is thought to be cumulatively accreted during episodic outbursts, which can reach as high as 35% in extreme cases (Dunham & Vorobyov 2012). Growing evidence suggests that accretion bursts may occur during the formation of massive stars as well (Caratti o Garatti et al. 2017; Meyer et al. 2017; Magakian et al. 2019). The magnitude of these eruptions is large enough to have substantial effects on disk dynamics, chemistry, mineralogy, dust properties, and snow lines, all of which have significant consequences for planet formation.

After its mass, the metallicity of a star (i.e., the content of species heavier than helium) has the most significant effect on its structure and evolution (e.g., Hansen et al. 2004). Hence, it is reasonable to assume that the effects of metal content for a PPD may be similarly profound. However, the observations related to low-metallicity PPDs remain poor and controversial, especially in their early stages. The dust-to-gas ratio in a PPD is difficult to measure observationally, and it is generally considered to be proportional to the metallicity of the host star (Murray et al. 2001; Ercolano & Clarke 2010). The mass accretion rate in young low-metallicity stars is observed to be higher than that in the corresponding solar-metallicity systems (Spezzi et al. 2012; De Marchi et al. 2013). The inner disk (up to a few astronomical units) fraction of low-metallicity YSOs is significantly smaller, suggesting that these disks are dispersed at an earlier stage and have a shorter lifetime due to enhanced photoevaporation—about 1 Myr, as compared to 5 Myr for solar-metallicity counterparts (Ercolano & Clarke 2010; Yasui et al. 2010). The frequency of giant planets is strongly correlated with the host star's metallicity, possibly because of an enhanced rate of planetesimal formation (Gonzalez 1997; Fischer & Valenti 2005; Johansen et al. 2009). Due to the difficulties related to time-domain astronomy at large distances, even less is known about episodic accretion in low-metallicity environments. All of the outbursting YSOs discovered so far have been found in star-forming regions within the immediate solar neighborhood, having approximately solar metallicity. Once the eruption begins, the stellar photosphere is not visible; hence, we cannot measure the metallicity of an ongoing event.

From a theoretical point of view, the episodic accretion in solar-metallicity disks has been studied extensively with the help of several numerical models (Bonnell & Bastien 1992; Bell & Lin 1994; Armitage et al. 2001; Vorobyov & Basu 2005; D'Angelo & Spruit 2010). Most relevant for this paper, Gammie (1996) showed that a typical PPD forms a magnetically "dead zone" at its midplane due to insufficient ionization to sustain magnetorotational instability (MRI) turbulence. As a consequence, the accretion occurs only through the active surface layers, with the dead zone forming an effective bottleneck in the mass and angular momentum transport. Accretion through such a layered disk structure can be intrinsically unsteady, and sudden activation of MRI in the dead zone can give rise to "MRI bursts" (Armitage et al. 2001; Zhu et al. 2010b; Kadam et al. 2020). In the early stages of its evolution, a PPD can also be prone to vigorous gravitational instability (GI), where the disk self-gravity forms large-scale spirals, as well as gravitationally bound clumps or fragments (Kratter & Lodato 2016, Toomre 1964). These clumps usually migrate inward on dynamical timescales, and their accretion onto the central protostar can trigger luminosity outbursts (Vorobyov & Basu 2015; Meyer et al. 2017; Zhao et al. 2018). Both of these mechanisms—MRI bursts and clump accretion—have been shown to be consistent with the observational constraints on the longer-duration, FUor-type outbursts.

When considering low-metallicity environments, numerical hydrodynamic simulations have been carried out primarily to study the effects on disk gravitational fragmentation and possible formation of gas giant planets through this process. Evolution of inviscid models suggested that a factor of 10 increase or decrease in metallicity does not affect the disk evolution significantly (Boss 2002). Cai et al. (2006) supported these findings that the outcome of hydrodynamic simulations is somewhat insensitive to the metallicity. However, in contradiction with Boss (2002), the cooling times in these simulations were too long for the disks to fragment and form clumps. In steady-state α-disk models, the low-metallicity disks tend to be more GI unstable (Tanaka & Omukai 2014). With hydrodynamic simulations, including prestellar cloud collapse and dust physics, Vorobyov et al. (2020) showed that metal-poor disks are also GI unstable. The low-metallicity models in this study showed that the duration of the burst phase of the disk due to clump accretion was much shorter than the solar-metallicity counterpart. In all of the low-metallicity studies so far, the underlying models assumed a fixed Shakura & Sunyaev (1973) α parameter for representing the turbulent viscosity in the disk. The dead zone in the disk was thus neglected, unless the evolution of the entire disk was carried out with a reduced viscosity (α parameter less than the canonical value of 10⁻²). The innermost disk region extending a few astronomical units, which is the most relevant in terms of observed FUor outbursts (Zhu et al. 2007), was also excluded from the computational domain.

In this paper, we present the effects of a low-metallicity environment on the structure and evolution of PPDs with a focus on the dead zone in the inner disk, as well as on the episodic accretion. We use numerical hydrodynamic simulations of PPDs in the thin-disk limit for these investigations. The simulations start with the collapse phase of the molecular cloud core so that the initial conditions and mass loading of the disk are accurately reproduced. The inner boundary of the computational domain is set at 0.42 au so that the sub–astronomical unit scale behavior of the disk can be captured. Note that the extent of our simulations, spanning distances from sub–astronomical unit scale to the radius of the parent molecular cloud, as well as temporal evolution over the first 0.7 Myr, push the limits of 2D hydrodynamic models. The dust grain content of a PPD increases with metallicity, which affects several physical processes. The dust dominates the opacity of a disk and thus controls its cooling properties and thermal equilibrium (Lodato 2008). We model a lower metallicity as a lower dust-to-gas ratio by scaling down the gas and dust opacities in proportion (Boss 2002; Cai et al. 2006). The low metal content offers a lower opacity to the stellar X-ray radiation and can thus allow a larger column density of the disk to be ionized (Hartmann et al. 2006). This can increase the thickness of the MRI-active surface layers. During the core-collapse phase, the decrease in dust continuum emission in the low-metallicity environment can lead to a higher background gas temperature (Vorobyov et al. 2020). We take into account the consequences of reduced metallicity in terms of all three effects: reduced opacity, increased active layer thickness, and increased cloud core temperature (Section 2). We find that the low-metallicity disks form progressively more restrictive dead zones in the inner disk and, as a consequence, are more massive (Section 3.2). The length of the phase during which the disk is prone to outbursts decreases at low metallicities, especially when the effects of inefficient cooling are taken into account (Sections 3.2 and 3.3). The outbursts are relatively infrequent in the disks around lower-mass stars and may be entirely suppressed at the lowest metallicities (Section 3.4). An individual outburst in a low-metallicity disk is typically more luminous and short-rising as compared to its solar-metallicity counterpart (Section 3.5). We conclude that the metal-poor environment has a significant effect on the eruptive behavior of PPDs; the disks tend to be more massive in the innermost regions and exhibit a shorter duration of the burst phase, confining it to the initial embedded phases of the YSO.

2. Model Description and Initial Conditions

In this section, we will describe the hydrodynamic model used for studying the metallicity effects described in Section 3. We will also elaborate on the initial conditions for the simulations and the physical reasons behind the chosen set of model parameters considered in this paper.

The model of PPD formation and evolution is based on the Formation and Evolution Of a Star And its circumstellar Disk (FEoSaD) code, where the full set of numerical hydrodynamic equations is solved in the thin-disk limit with a cylindrical geometry (Vorobyov & Basu 2010, 2015; Vorobyov et al. 2018). FEoSaD has several physics modules that can be turned on depending on the problem at hand in order to find a balance between the gained insights and the computational cost. The following hydrodynamic equations of mass, momentum, and energy transport are solved:

$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}=-{{\boldsymbol{\nabla }}}_{p}\cdot \left({\rm{\Sigma }}{{\boldsymbol{v}}}_{p}\right),\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left({\rm{\Sigma }}{{\boldsymbol{v}}}_{p}\right)+{\left[{\boldsymbol{\nabla }}\cdot \left({\rm{\Sigma }}{{\boldsymbol{v}}}_{p}\otimes {{\boldsymbol{v}}}_{p}\right)\right]}_{p}=-{{\boldsymbol{\nabla }}}_{p}P+{\rm{\Sigma }}{{\boldsymbol{g}}}_{p}+{\left({\boldsymbol{\nabla }}\cdot {\boldsymbol{\Pi }}\right)}_{p},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{{\boldsymbol{\nabla }}}_{{\rm{p}}}\cdot \left({{ev}}_{p}\right)=-P\left({{\boldsymbol{\nabla }}}_{p}\cdot {v}_{p}\right)-{\rm{\Lambda }}+{\rm{\Gamma }}+{\left({\boldsymbol{\nabla }}\cdot v\right)}_{{pp}^{\prime} }\,:{{\boldsymbol{\Pi }}}_{{pp}^{\prime} }.\end{eqnarray} \tag{ 3 }$

Here the subscripts p and $p^{\prime}$ refer to the planar components in polar coordinates (r, ϕ), Σ is the surface mass density, e is the internal energy per unit area, and P is the vertically integrated gas pressure. Among nonscalars, ${{\boldsymbol{v}}}_{p}={v}_{{\rm{r}}}\hat{{\boldsymbol{r}}}+{v}_{\phi }\hat{{\boldsymbol{\phi }}}$ is the velocity in the disk plane, ${{\boldsymbol{g}}}_{p}={g}_{{\rm{r}}}\hat{{\boldsymbol{r}}}+{g}_{\phi }\hat{{\boldsymbol{\phi }}}$ is the gravitational acceleration in the disk plane, and ${{\boldsymbol{\nabla }}}_{{\rm{p}}}=\hat{{\boldsymbol{r}}}\partial /\partial r+\hat{{\boldsymbol{\phi }}}{r}^{-1}\partial /\partial \phi$ is the gradient along the planar coordinates of the disk. The ideal equation of state, P = (γ − 1)e with γ = 7/5, is used for calculating the gas pressure. The viscous stress tensor,

$\begin{eqnarray}&&{\boldsymbol{\Pi }}=2{\rm{\Sigma }}\nu \left({\boldsymbol{\nabla }}v-\displaystyle \frac{1}{3}({\boldsymbol{\nabla }}\cdot v){\boldsymbol{e}}\right),\end{eqnarray} \tag{ 4 }$

accounts for turbulent viscosity in the disk, where ν is the kinematic viscosity, ∇v is a symmetrized velocity gradient tensor, and e is a unit tensor.

The cooling and heating rates Λ and Γ, respectively, in Equation (3) are based on the analytical solution of the radiation transfer equations in the vertical direction (Dong et al. 2016). The cooling function per surface area of the disk is expressed as

$\begin{eqnarray}&&{\rm{\Lambda }}=\displaystyle \frac{8{\tau }_{{\rm{P}}}\sigma {T}_{\mathrm{mp}}^{4}}{1+2{\tau }_{{\rm{P}}}+\tfrac{3}{2}{\tau }_{{\rm{R}}}{\tau }_{{\rm{P}}}},\end{eqnarray} \tag{ 5 }$

where σ is the Stefan–Boltzmann constant, τ_P and τ_R are the Planck and Rosseland optical depths to the disk midplane, and κ_P and κ_R are the Planck and Rosseland mean opacities, respectively. The opacities are calculated for physical conditions typical of PPDs, which include both dust and gas components (Semenov et al. 2003). The heating rate is given by

$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{8{\tau }_{{\rm{P}}}\sigma {T}_{\mathrm{irr}}^{4}}{1+2{\tau }_{{\rm{P}}}+\tfrac{3}{2}{\tau }_{{\rm{R}}}{\tau }_{{\rm{P}}}},\end{eqnarray} \tag{ 6 }$

where T_irr is the irradiation temperature at the disk surface,

$\begin{eqnarray}&&{T}_{\mathrm{irr}}^{4}={T}_{\mathrm{bg}}^{4}+\displaystyle \frac{{F}_{\mathrm{irr}}(r)}{\sigma },\end{eqnarray} \tag{ 7 }$

where T_bg is the temperature of the background blackbody irradiation, and F_irr(r) is the radiation flux absorbed by the disk surface at a radial distance r from the central star. The latter is calculated as

$\begin{eqnarray}&&{F}_{\mathrm{irr}}(r)=\displaystyle \frac{{L}_{* }}{4\pi {r}^{2}}\cos {\gamma }_{\mathrm{irr}},\end{eqnarray} \tag{ 8 }$

where γ_irr is the incidence angle of radiation arriving at the flaring disk surface at a radial distance r. The stellar luminosity L_* is the sum of the accretion and photospheric luminosities. The accretion luminosity, ${L}_{* ,\mathrm{accr}}=(1-\epsilon ){{GM}}_{* }\dot{M}/2{R}_{* }$ , was generated from the accreted gas, where the fraction of accretion energy absorbed by the star () was set to 0.05. The photospheric luminosity, L_*,ph, was due to both the gravitational compression and the deuterium burning in the stellar interior (Vorobyov & Basu 2010). The stellar mass, M_*, and accretion rate onto the star, $\dot{M}$ , are determined using the amount of mass passing through the inner computational boundary. The properties of the forming protostar (L_*,ph and radius R_*) are calculated using the pre-main-sequence stellar evolution tracks of D'Antona & Mazzitelli (1997). The effects of metallicity on stellar evolution were not taken into account in this study.

The MRI is considered to be the primary driver in producing turbulent viscosity in PPDs (Hawley et al. 1995; Turner et al. 2014). The disk material needs to be sufficiently ionized for MRI to operate. The shearing Keplerian motion of the gas coupled with the field causes turbulence and angular momentum transport. In a typical PPD, the midplane temperature outside of about 1 au radius is not high enough to sustain collisional ionization (Armitage 2011). Galactic cosmic rays are considered to be the major source of ionization at a few astronomical units and penetrate an approximately constant column density of the gas (Umebayashi & Nakano 1981). The accretion disk thus forms a layered structure; accretion occurs only through the sufficiently ionized, MRI-active surface layers, and a magnetically dead zone is formed at the midplane (Gammie 1996).

We model the accretion through a layered PPD with an effective and adaptive Shakura & Sunyaev (1973) α-parameter (Bae et al. 2014; Kadam et al. 2019). The kinematic viscosity is parameterized as ν = α_eff c_s H, where c_s is the sound speed and H is the vertical scale height of the disk, which is computed assuming a local hydrostatic equilibrium taking disk self-gravity into account (Vorobyov & Basu 2009). The α_eff is given by

$\begin{eqnarray}&&{\alpha }_{\mathrm{eff}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{a}}}{\alpha }_{{\rm{a}}}+{{\rm{\Sigma }}}_{{\rm{d}}}{\alpha }_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{a}}}+{{\rm{\Sigma }}}_{{\rm{d}}}},\end{eqnarray} \tag{ 9 }$

where Σ_a is the gas surface density of the MRI-active surface layer of the disk and Σ_d is the gas surface density of the magnetically dead layer at the disk midplane. Note that the total gas surface density is then Σ = 2 × (Σ_a + Σ_d). The parameters α_a and α_d are proportional to the strength of turbulent viscosity in the MRI-active and MRI-dead layers of the disk, respectively. We set a canonical value of α_a = 0.01 for the magnetically active region. However, there is some uncertainty in α_a, and we will explore the possibility of a higher value in the future, as suggested by recent magnetohydrodynamic simulations (Zhu et al. 2020). The parameter α_d is defined as

$\begin{eqnarray}&&{\alpha }_{{\rm{d}}}={\alpha }_{\mathrm{MRI},{\rm{d}}}+{\alpha }_{\mathrm{rd}},\end{eqnarray} \tag{ 10 }$

where

$\begin{eqnarray}&&{\alpha }_{\mathrm{MRI},{\rm{d}}}=\left\{\begin{array}{c}{\alpha }_{{\rm{a}}},\,\mathrm{if}\,{T}_{\mathrm{mp}}\gt {T}_{\mathrm{crit}}\\ 0,\,\mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 11 }$

where T_crit is the MRI activation temperature, and T_mp is the disk midplane temperature. The viscosity thus sharply rises to the fully MRI-active value above T_crit. This models the almost exponential increase in the disk ionization fraction due to the thermal effects, such as dust sublimation and ionization of alkali metals. The dead zone in a PPD can have a small but nonzero residual viscosity due to hydrodynamic turbulence driven by the Maxwell stress in the disk active layer (Okuzumi & Hirose 2011; Bae et al. 2014). We set this residual viscosity using

$\begin{eqnarray}&&{\alpha }_{\mathrm{rd}}=\min \left({10}^{-5},{\alpha }_{{\rm{a}}}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{a}}}}{{{\rm{\Sigma }}}_{{\rm{d}}}}\right).\end{eqnarray} \tag{ 12 }$

As the nonzero α_rd is due to turbulence propagating from the active layer down to the disk midplane, the above expression ensures that the accretion in the dead zone cannot exceed that of the active layers.

The numerical simulations start from the gravitational collapse of a starless molecular cloud core. The protostar is formed within the inner computational boundary of the disk, while the envelope continues to accrete inside the centrifugal radius during the embedded phase. The initial surface density and angular velocity profiles of the cloud core are derived from an axisymmetric core compression where the angular momentum remains constant and magnetic fields are expelled due to ambipolar diffusion (Basu 1997),

$\begin{eqnarray}&&{\rm{\Sigma }}=\displaystyle \frac{{r}_{0}{{\rm{\Sigma }}}_{0}}{\sqrt{{r}^{2}+{r}_{0}^{2}}},\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&{\rm{\Omega }}=2{{\rm{\Omega }}}_{0}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}\left[\sqrt{1+{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{2}}-1\right],\end{eqnarray} \tag{ 14 }$

where Σ₀ and Ω₀ are the maximum values at the center of the core, and r₀ is the radius of the central plateau proportional to the thermal Jeans length. The initial cloud cores were constructed such that the ratio of rotational to gravitational energy, β, was approximately 0.14%. This value is consistent with the observations of prestellar cores (Caselli et al. 2002) and was chosen to be relatively low in order to inhibit the formation of self-gravitating clumps during the early, massive phases of the disk evolution (Vorobyov 2013). This was done to avoid possible interference with another accretion burst model caused by clump infall, which was shown to operate in solar and low-metallicity environments (Vorobyov et al. 2020).

The treatment of the inner boundary conditions of the numerical simulations of PPDs needs special attention. A complication may arise if the inner boundary allows for matter to flow in only one direction, i.e., an outflow-only boundary condition from the disk to the sink cell. The wavelike motions near the inner boundary result in a disproportionate flow through the sink–disk interface and cause an artificial drop in the gas density near the inner boundary. We consider a carefully implemented inflow–outflow inner boundary condition for our simulations (Vorobyov et al. 2018). Here the material that passes to the sink cell through the sink–disk interface is redistributed between the central protostar and the sink cell. Depending on the mass surface density and velocity gradients at the boundary, the material is allowed to flow from the sink cell back onto the active disk. This method ensures that the innermost parts of the disk are unaffected by the proximity of the boundary and related numerical artifacts. For all simulations presented in this study, the inner boundary is placed at 0.42 au. This allows us to evolve the innermost parts of the disk with an accuracy sufficient for capturing the sub–astronomical unit scale behavior of the PPD.

For this study of PPDs in low-metallicity environments, we conducted a total of nine hydrodynamic simulations as listed in Table 1. We considered two masses of the parental cloud core—1.152 and 0.576 M_⊙—which will eventually form a solar-type dwarf star and a fully convective low-mass star, respectively. The prescripts "MS" and "ML" describe the final stellar mass. The fiducial models with a postscript "fid" are the standard models, conducted with a composition of solar metallicity against which the low-metallicity results are compared. We consider two values of low metallicities. The postscript "Z0.1" corresponds to models with 10% of the solar metallicity, such as the star-forming regions observed in the outer Galaxy. The extremely low metallicity of 2% of the solar value ("Z0.02" models) may correspond to the young galaxies undergoing their first major episodes of star formation. In accordance with previous numerical studies (e.g., Boss 2002; Cai et al. 2006), we set the subsolar metallicity by scaling down the gas and dust opacities of the solar-metallicity disk by the corresponding factor. Although the thickness of the MRI-active surface layers as determined by the Galactic cosmic rays is uniform at about 100 g cm⁻² (Umebayashi & Nakano 1981), the lower dust content in the low-metallicity PPDs may offer reduced rates of recombination on grain surfaces, resulting in a larger penetration depth (Yasui et al. 2009). We thus demonstrate the effects of the increased MRI efficiency with two simulations (denoted by "Σa" in the model name) with a higher value of 200 g cm⁻². In the case of the lowest-metallicity environment of about Z = 1% Z_⊙, the dominant cooling mechanism of the dust continuum emission is much less effective during the initial stages of the core-collapse phase (Omukai et al. 2005). This results in a higher background gas temperature as compared to the higher-metallicity counterparts and higher disk infall velocities during the pre- and protostellar stages (Vorobyov et al. 2020). To investigate the effects of this phenomenon, we conducted two additional simulations with an increased initial gas temperature of the prestellar cloud core of 25 K, which was also set to the temperature of the background radiation (indicated by "Tc" in the model name). The default value of these temperatures was set to 15 K for the rest of the simulations. The last column of Table 1 summarizes the included physical effects in each model with a self-explanatory notation.

Table 1. List of Simulations

Model Name	${M}_{\mathrm{core}}({M}_{\odot })$	Σ_a (g cm⁻²)	${T}_{\mathrm{core}}\,({\rm{K}})$	Z/Z_⊙	Included Effects
`MS`_`fid`	1.152	100	15	1.0	Fiducial solar-mass model

`MS`_`Z0.1`	1.152	100	15	0.1	κ only
`MS`_`Z0.02`	1.152	100	15	0.02

`MS`_Σ`a`_`Z0.1`	1.152	200	15	0.1
`MS`_Σ`a`_`Z0.02`	1.152	200	15	0.02	κ + Σ_a

`MS`_`Tc`_`Z0.02`	1.152	100	25	0.02	$\kappa +{T}_{\mathrm{core}}$

`ML`_`fid`	0.576	100	15	1.0	Fiducial lower-mass model

`ML`_`Z0.1`	0.576	100	15	0.1	κ only

`ML`_`Tc`_`Z0.02`	0.576	100	25	0.02	$\kappa +{T}_{\mathrm{core}}$

Download table as: ASCII Typeset image

The simulations were conducted for 0.7 Myr, and at the end, the system may be considered to be in the T Tauri phase where the envelope has already accreted onto the star–disk system and most of the burst activity is over. Here the time is measured from the start of the gravitational collapse of the core, i.e., the beginning of a simulation. The radial and azimuthal resolution of the computational grid for the simulations was set to 256 × 256 for all simulations, with a logarithmic spacing in the radial direction and uniform spacing in the azimuthal direction. The maximum grid resolution in the radial direction was 1.6 × 10⁻² au at the innermost cell of the disk (0.417 au), while the poorest resolution near the outer boundary of the cloud (10,210 au for MS simulations) was about 400 au. Note that any two simulations conducted with the same initial conditions and parameters do not produce identical disk evolution because of gravitationally unstable disks; however, the simulation results are always congruent, and we have tested their convergence with the help of higher-resolution models.

3. Results

In this section, we will present our results of disks formed in low-metallicity environments. In Section 3.1, we show the differences in disk structure and evolution on a global scale. In Section 3.2, only the effects of opacity are considered, while in Section 3.2, we elaborate on the added effects of increased active layer thickness and cloud core temperature. In Section 3.4, we consider low-mass cloud cores; finally, in Section 3.5, an individual MRI burst is characterized. Table 2 summarizes some of the key results from all nine simulations for comparison, which will be explained in the upcoming sections.

Table 2. Summary of Simulation Results

Model Name	M_*,F (M_⊙) ^a	M_disk,F (M_⊙)	M_{disk<10 au,F} (M_⊙)	Embedded Phase (Myr)	Burst Phase (Myr)
`MS`_`fid`	0.816	0.227	0.015	0.242	0.552
`MS`_`Z0.1`	0.772	0.282	0.078	0.242	0.414
`MS`_`Z0.02`	0.782	0.285	0.084	0.242	0.274
`MS`_Σ`a`_`Z0.1`	0.821	0.220	0.011	0.242	0.415
`MS`_Σ`a`_`Z0.02`	0.804	0.236	0.048	0.243	0.220
`MS`_`Tc`_`Z0.02`	0.840	0.187	0.082	0.113	0.129
`ML`_`fid`	0.459	0.042	2.34 × 10⁻³	0.095	0.25
`ML`_`Z0.1`	0.458	0.044	3.38 × 10⁻³	0.096	0.18
`ML`_`Tc`_`Z0.02`	0.478	0.022	6.05 × 10⁻³	0.044	⋯

Note.

^aThe subscript "F" denotes final values measured at 0.7 Myr.

Download table as: ASCII Typeset image

3.1. Large-scale Evolution

As we shall see later, a lower metallicity has a significant effect on the innermost parts of a PPD. However, in certain cases, the large-scale structures and global evolution of the disk are also noticeably affected. Figure 1 shows the evolution of the disk in four simulations—MS_fid, MS_Z0.02, MS_Tc_Z0.02, and ML_fid—over a region of 500 × 500 au² during the first 0.5 Myr. Note that the time duration between the snapshots is not uniform. The yellow contours mark the Σ = 1 g cm⁻² level, indicating the approximate extent of the disk. Consider the solar-mass, solar-metallicity model, MS_fid, in the top row. The disk was formed around 0.08 Myr (therefore not observed in the first frame) and then showed viscous spread as it evolved. The initial phase was dominated by vigorous GI instability and associated large-scale spirals. The initial ratio of kinetic to gravitational energy was chosen to be low enough that the disks will typically form no clumps. As seen in the last panel at 0.5 Myr, the disk eventually became weakly unstable and progressively more axisymmetric. This approximate symmetry was maintained until the end of the simulation. Evolution on this large scale for models MS_Z0.1, MS_Z0.02, MS_Σa_Z0.1, and MS_Σa_Z0.02 was very similar to the MS_fid model with two key differences. Even at this scale, the innermost regions of ⪅10 au showed a consistently larger gas surface density (e.g., similar to the MS_Tc_Z0.02 simulation in the last snapshot). The lower-metallicity models were also prone to more GI activity and associated spirals, although no fragmentation was observed. The second row of Figure 1 shows the evolution for model MS_Z0.02, where both of these differences can be noticed. The metal-poor environment offers a lower opacity and hence more efficient cooling of the disk, which can, in turn, favor GI. See Vorobyov et al. (2020) for a detailed study of accretion bursts caused by GI fragments in low-metallicity PPDs.

**Figure 1.** Evolution of the disk gas surface density distribution (in units of log₁₀ g cm⁻²) for the four models—`MS`_`fid`, `MS`_`Z0.02`, `MS`_`Tc`_`Z0.02`, and `ML`_`fid`—over a region of 500 × 500 au, showing the large-scale disk structure. The yellow contours mark Σ = 1 g cm⁻² levels, indicating the approximate outer boundary of the disks.
Download figure:
Standard image High-resolution image

Among the solar-mass simulations, the most significant difference in the large-scale disk evolution was observed for the MS_Tc_Z0.02 simulation. The increased cloud core temperature (25 K as opposed 15 K for the other MS models) affects the thermal evolution of the disk, resulting in larger infall velocities (Vorobyov et al. 2020). As a consequence, the mass infall rate on the disk was higher in the early stages, and the disk showed faster growth, as well as a shorter period of GI activity. As seen in the second row of Figure 1, the disk was already formed at 0.05 Myr and appeared more developed at 0.35 Myr than MS_fid at 0.5 Myr. We will elaborate on this accelerated evolution of MS_Tc_Z0.02 in Section 3.3. The last row of Figure 1 shows the evolution of surface density for the lower-mass solar-metallicity model ML_fid. Disk evolution for a lower-mass cloud core occurs on a faster timescale because of the limited mass reservoir in the parental cloud core. The disk showed some GI activity and associated spirals in the initial stages, but no clumps were formed in any of the low-mass models. The extent of the disk was smaller at later times and showed less accumulation of gas in the innermost regions. The trends for the rest of the lower-metallicity simulations were similar to the higher-mass models. The model ML_Z0.1 showed marginal differences at a large scale in the direction of more vigorous GI, while ML_Tc_Z0.02 showed an accelerating evolution. In metal-poor environments, most notable effects occur at smaller scales at less than 10 au of the inner disk, which we will focus on next.

3.2. Effects of Opacity

In this section, we elaborate on the general features of the fiducial solar-metallicity simulation MS_fid. We compare this to the low-metallicity simulations MS_Z0.1 and MS_0.02 with Z = 0.1 and 0.02 Z_⊙, respectively. Note that in the latter models, a metal-poor environment is emulated by scaling down the gas and dust opacities by the same fraction. In Figure 2, the evolution of the inner disk is shown with the help of spacetime diagrams. The ordinate extends from the inner radius of the computational domain of 0.42–500 au on a logarithmic scale, while the abscissa is time in megayears. The left column represents simulation MS_fid with solar metallicity, and the middle and right columns show models MS_Z0.1 and MS_Z0.02, respectively. The rows compare the azimuthally averaged profiles of quantities: gas mass surface density, effective α, and midplane temperature.

Consider the left column for the simulation MS_fid. In the top row, the green contour plotted at Σ = 1 g cm⁻² shows the approximate extent of the disk. The yellow contour marks the 2 × Σ_a = 200 g cm⁻² level, incorporating both sides of the disk; thus, the region outside this boundary is fully MRI active. The disk formed at about 0.08 Myr, and axisymmetric ring structures soon began appearing in the gas surface density. In Kadam et al. (2019), we have described in detail the evolution of the rings, which form at a distance of a few astronomical units due to viscous torques occurring at the inner boundary of the dead zone. The simulation model1_T1300_S100 in Kadam et al. (2019) was identical to MS_fid but conducted with twice the spatial resolution. The salient features of MS_fid remain consistent with the earlier results, indicating convergence of the numerical simulations. We will demonstrate the congruency across spatial resolutions in detail by comparative analysis of these two simulations in the Appendix. The middle row of Figure 2 compares α_eff, which has a constant value of 0.01 for the fully MRI-active disk in the outer regions. The black contour shows the approximate extent of the dead zone, defined as the region where α_eff is below 80% of this maximum value. The yellow contours mark the 10⁻⁴ level, thus indicating the regions of low viscosity that coincide with the gaseous rings. The bottom row shows the midplane temperature with yellow contours marking the MRI activation temperature, T_crit = 1300 K. The midplane temperature in the innermost regions frequently exceeded this value for a short time, and some of these temperature increases were associated with MRI outbursts.

When comparing the effects of lower metallicity on the disk structure in Figure 2, regions outside of about 10 au were very similar for all three simulations. The green contour in Σ indicating the disk extent and the yellow contour in the same plots marking the fully MRI-active outer disk showed almost identical evolution. However, the inner disk structure was markedly different. At a lower metallicity, clear formation of successive rings was not observed; instead, a progressively broader region of gas accumulation was formed. The surface density of the accumulated gas increased as the metallicity decreased. The α_eff essentially varies inversely to the total Σ (Equation (9)); hence, the lower metallicity resulted in a deeper dead zone in terms of disk viscosity. This is clearly observed with the α_eff = 10⁻⁴ contours in the second row. The dead zone also became broader and moved inward, in general, at low metallicities. With a decrement in metallicity, the midplane temperature of MS_Z0.1 and MS_Z0.02 tended to be proportionally cooler as compared to the fiducial model. The opacity was assumed to be proportional to the dust content and hence the metallicity of the disk. With a lower opacity, the midplane of the disk can cool more efficiently, lowering the midplane temperature. The disk viscosity is not only proportional to α_eff but also to the local sound speed. With a decreasing temperature, the sound speed also decreased, resulting in a further decrement in viscosity. This lower viscosity formed a narrower bottleneck to the mass transport, thus contributing to the observed trends with a lower metallicity—accumulation of more gas, as well as lower α_eff in the inner disk region.

In Figure 3, we compare the outbursting activity of the same three simulations—MS_fid, MS_Z0.1, and MS_Z0.02—with respect to some of the time-dependent quantities. The first row compares the accretion rate onto the central star and infall rate onto the disk from the cloud core ( $\dot{{M}_{* }}$ and ${\dot{M}}_{\mathrm{infall}}$ , respectively). The three boxes (A, B, and C) highlight the three modes of accretion variability that will be discussed later in this section. The $\dot{{M}_{* }}$ is calculated as mass passing through the sink cell of the inflow–outflow inner boundary, with the ratio of mass accreting onto the star to that going into the sink cell set to 95%:5% (Kadam et al. 2019). Our model also assumes that 10% of the mass that passes through the sink cell is lost via jets or outflows before this partitioning. The $\dot{{M}_{\mathrm{infall}}}$ is considered to be the mass flux at a radial distance of 1000 au from the center. The second row shows the total and stellar photospheric luminosities, L_total and L_*, respectively. The total luminosity is the sum of accretion and the stellar photospheric luminosity, where the latter quantity is calculated through the precomputed tracks of D'Antona & Mazzitelli (1997). The third row compares the midplane temperature at the inner boundary of the disk (T_inner). The fourth row shows the evolution of the stellar and disk mass (M_* and M_disk, respectively), as well as the mass of the inner disk less than 10 au (M_{disk,10 au}). In this study, we define the end of the "burst phase" of a PPD as the time up to the last major outburst observed in the simulation. We consider the end of the "embedded phase" of the disk as the time when the infalling envelope retains less than 10% of the initial core mass. Table 2 lists the final masses and durations of the above two phases at the end of all simulations. Although the dead zone may continue to exist, it is unlikely that the systems will show additional large-amplitude, FUor-like outbursts after the termination of the simulations at 0.7 Myr, i.e., after the initial burst phase is over (Zhu et al. 2010a; Bae et al. 2014).

Consider the first row of Figure 3. The overall qualitative behavior of the mass accretion history is consistent with previous 1D and 2D simulations (Zhu et al. 2010a; Bae et al. 2014; Vorobyov & Basu 2015). The outbursts were superimposed on a steadily decreasing background accretion rate that stabilized at a few times 10⁻⁷ M_⊙ yr⁻¹ at later times. The early burst phase was characterized by high accretion variability reflecting the dynamical nature of the disk, primarily the presence of GI-induced spiral waves. As the mass infall rate decremented with time, the mass loading of the disk declined. Thus, vigorous GI activity could not be sustained and the accretion variability also diminished. This variability is best observed in the MS_Z0.02 simulation and highlighted by box C in the figure. Typical FU Orionis–like eruptions are characterized by a sudden increase in $\dot{{M}_{* }}$ by more than an order of magnitude. These large-amplitude eruptions are MRI outbursts related to the sudden triggering of MRI in the inner disk (Kadam et al. 2020). This mode of accretion variability caused by MRI outbursts is highlighted by box B. The arrow marks the individual outburst that will be investigated in detail in the context of the metal-poor environment in Section 3.5. The MRI bursts were clearly distinguishable from the GI-related variability. The outbursts were initially clustered together, superimposed on the accretion variability. The duration between two outbursts increased with time, eventually leading to termination of the burst phase. The simultaneous presence of both MRI and GI in the disk raises the possibility of the recently proposed "spiral wave dynamo," especially in the early times (Riols & Latter 2018, 2019; Deng et al. 2020). The GI spiral waves can initiate a magnetohydrodynamic dynamo in the disk that can amplify the initial magnetic field. The associated magnetic torques can drive mass accretion and may revive the dead zone even before the MRI outburst is triggered. The action of such nonideal magnetohydrodynamic effects should be considered in future investigations.

When considering how the lower metallicity affects accretion rates, two trends are immediately apparent. First, the duration of the burst phase is shortened with a decrease in the metal content. The burst phase ended at approximately 0.55, 0.42, and 0.27 Myr for MS_fid, MS_Z0.1, and MS_Z0.02, respectively. With the duration of the embedded phase remaining constant at about 0.24 Myr, most of the outbursting activity in the lower-metallicity models was hidden from direct view. The second observation is that the amplitude of individual bursts increased with decreasing metallicity. This can also be observed in the luminosity plots in the second row. All of the bursts in the MS_fid simulation were less than 100 M_⊙ in luminosity output and became progressively more powerful, typically above 100 M_⊙ for the lower-metallicity counterparts. We conjecture that the reason behind the shortening of the burst phase with decreasing metal content was the overall decrease in temperature of the innermost parts of the disk. With the inclusion of dead-zone viscosity, the MRI burst is triggered by the viscous heating near the inner edge of the dead zone and propagates outward in an inside-out manner (see Section 3.5). As seen in Figure 2, the metal-poor disks cool very efficiently, which can inhibit such triggering of the outbursts. In the same figure, the spacetime plot of the gas mass surface density shows accumulation of a larger amount of mass at lower metallicities, thus making more material available for accretion during an outburst and explaining the trend in the intensity of the bursts.

The third row of Figure 3 shows the temperature at the inner boundary of the disk. For the lower-metallicity models, the luminosity bursts coincided with the occasional crossings of T_inner above the critical temperature of 1300 K. This is expected for an MRI outburst as the underlying mechanism depends on temperature-dependent viscosity. For the MS_fid simulation, however, the temperature at the inner boundary constantly fluctuated across this critical threshold. Here we mention one caveat of our model. The computational domain of the simulations encompasses a region outside of 0.42 au, and despite a carefully constructed inflow–outflow boundary, the simulation results reflect the conditions at this location in the disk. The real mass accretion rate onto the star may be modified by the physical conditions and mechanisms that may operate inside this innermost disk region, such as magnetospheric accretion in low-mass YSOs. In the case of layered disks, such as modeled here, as we move closer to the central star, the disk changes its state from having a dead zone to a fully MRI-active region due to thermal ionization of the gas. Part of the fluctuations that are observed in T_inner and, consequently, the accretion rate onto the star are due to the radial excursion of the inner boundary of the dead zone across the inner boundary of the computational domain. On the one hand, when the boundary is placed further out, important phenomena occurring in the inner few astronomical units, such as MRI outbursts, are not captured at all. On the other hand, additional complications may arise when the boundary is placed sufficiently close to the star. The latter case seems to hold true for MS_fid, where it is difficult to differentiate between real accretion activity and movement of the inner fully MRI-active zone. This third mode of irregularity in the mass accretion rate is highlighted by box A in Figure 3. The magnitude of this accretion variability was typically a few times smaller than the MRI outbursts. The inner edge of the dead zone moved inward and away from the computational boundary for the metal-poor disks due to the lower temperatures. Thus, this anomalous variability was diminished for the lower-metallicity models. This remains a general drawback of disk simulations that do not model the disk reaching all the way to the stellar surface or magnetosphere. The overall results of this study are not affected by these fluctuations, as our results remain consistent with the previous investigations.

The fourth row of Figure 3 compares the evolution of the stellar and disk masses. For all three simulations, the disk remained substantially massive during the burst phase and was a large fraction of the stellar mass at the end of the simulations at 0.7 Myr. With a lower metallicity, the final stellar mass marginally decreased, while the total disk mass was marginally larger (see Table 2 for numerical values). Although the total final disk masses were comparable, the innermost parts were substantially massive at low metallicities. At the end of the simulations, only 6% of the total disk mass in MS_fid was contained within the inner 10 au region, while about 30% of the disk mass was in the inner 10 au for the lowest-metallicity MS_Z0.02 model. This indicates that the dead zone in a low-metallicity environment formed a narrower bottleneck to the mass transport, affecting the long-term evolution of the system. The observations suggest that the frequency of the disk harboring stars (i.e., disk fraction) decreases significantly at low metallicity, indicating a shorter disk lifetime (Yasui et al. 2010). Our finding of more massive inner disks in metal-poor environments suggests that the disk dispersal via photoevaporation occurring at late stages of these disks may be more efficient than previously estimated (e.g., Ercolano & Clarke 2010). The frequency of close-in super-Earths is observed to be almost independent of the host star's metallicity (Petigura et al. 2018). The massive inner disk may facilitate dust accumulation, its growth and formation of planetesimals, explaining this trend. The postburst phase accretion rates were similar across the three models near the end of the simulations. The transport of accumulated mass in the inner disk at later times may be responsible for the observed trends of higher accretion rates at low metallicities (Spezzi et al. 2012; De Marchi et al. 2013).

3.3. Active Layer Thickness and Cloud Core Temperature

For a magnetically layered PPD, the MRI turbulence depends on the ionization degree of the disk gas, which is essentially determined by the details of the ionization and recombination rates. Galactic cosmic rays are considered to be the major source of ionization at the distance of a few astronomical units from the central star that ionize a nearly uniform gas surface density of about 96 g cm⁻² (Umebayashi & Nakano 1981). However, there is considerable uncertainty, as the strong winds of a protostar or accretion outflows can shield the inner disk (Cleeves et al. 2013). At median T Tauri luminosities, the combined effect of stellar far-UV photons and X-rays is at least an order of magnitude smaller (Bergin et al. 2007; Cleeves et al. 2013). Thus, a canonical value of 100 g cm⁻² is used for the thickness of the MRI-active surface layer in our models presented in Section 3.2. In a low-metallicity environment, we can expect a proportionally lower dust content to be present in the disk. As the recombination on grain surfaces is reduced, this would result in a thicker active layer as compared to solar-metallicity disks (Fromang et al. 2002). However, a low metallicity also results in fewer molecular ions in the disk, which contribute to a reduced ionization fraction (Sano et al. 2000). It is speculated that in this tug of war, the former process should dominate and a thicker active layer would form (Yasui et al. 2009). The exact effects are difficult to calculate due to their dependence on long-term dust, as well as the chemical evolution of the disk, and have not been published in the literature. In this section, we assumed that the active layer thickness increases and approximately doubles in metal-poor disks. Although the magnitude is arbitrary, the comparison will give us insight into how the increase in active layer thickness can affect the disk evolution.

Another effect considered in this section is the increase in the initial temperature of the molecular cloud core in a low-metallicity environment. With global numerical simulations of PPD formation that incorporate a sophisticated heating and cooling treatment with separate gas and dust temperatures, Vorobyov et al. (2020) showed that in the case of extremely metal-poor disks (Z = 0.01 Z_⊙), the gas and dust temperatures show significant decoupling during the core-collapse phase. This is because at low densities, the dust continuum emission is a dominant cooling mechanism (Omukai et al. 2005). The disks with extremely low metallicity cannot cool efficiently, leading to an increase in the gas temperature, as well as greater infall velocities than their higher-metallicity counterparts. Thus, for the lowest-metallicity case, we also conducted simulations with an increased cloud core temperature of 25 K, as opposed to a canonical value of 15 K (see Table 1).

In Figure 4, we present the low-metallicity models with an increased active layer thickness—MS_Σa_Z0.1 and MS_Σa_0.02—and the MS_Tc_Z0.02 simulation, which includes evolution with an increased cloud core temperature. In order to demonstrate the effects on an increased Σ_a only, we compare the first two models with their canonical counterparts, MS_Z0.1 and MS_0.02, respectively, presented in Section 3.2. Consider the top row of Figure 4, showing the spacetime diagram of gas surface density. The yellow contours in the first two models show 2 × Σ_a = 200 g cm⁻² levels. The plots show a lesser accumulation of gas in the vicinity of the dead zones for both the MS_Σa_Z0.1 and MS_Σa_0.02 simulations. When considering the spacetime diagrams of α_eff, MS_Σa_Z0.1 did not show the yellow 10⁻⁴ contours, indicating a shallower depth of the dead zone. The radial extent of the dead zone was also smaller as compared to MS_Z0.1. Similar trends were observed for MS_Σa_0.02, where the dead zone that formed was smaller in terms of both radial and temporal extent. With the doubling of the active layer, the dead zone formed a less restrictive bottleneck for angular momentum and mass transport. From Equation (9), the above trends in both surface density and α_eff are as expected.

The bottom row of Figure 4 shows the midplane temperature, which was marginally increased in the innermost parts, i.e., in the vicinity of 1 au, for both MS_Σa_Z0.1 and MS_Σa_0.02, as compared to the corresponding models with Σ_a = 100 g cm⁻². This can be explained by a general increase in viscous heating caused by a thicker active layer. This indicates a negligible effect of an increased Σ_a on the evolution of the outer disk.

We compare the third simulation, MS_Tc_Z0.02, with MS_Z0.02 for demonstrating the effects of increased core temperature. In Figure 4, the evolution of the gas surface density, the extent of the dead zone as measured by α_eff, and the inner disk temperatures look very similar for the two cases, with one key difference. The disk was formed at an earlier time at approximately 0.02 Myr in MS_Tc_Z0.02, as compared to the canonical counterpart, where the disk was formed at 0.08 Myr. The mass infall rate of the cloud core is proportional to its gas temperature to the power 3/2 (Shu 1977). The higher gas temperature (10 K higher than the standard value of 15 K used for MS_Z0.02) in the inner parts of the collapsing core caused a proportional increase in the mass infall rate. The resulting accelerated evolution of the system manifested as an early formation of the disk.

In Figure 5, we present the quantities related to the episodic accretion of the same three models—MS_Σa_Z0.1, MS_Σa_0.02, and MS_Tc_Z0.02. Consider the first two models, where the overall trend in accretion history was similar to canonical models with Σ_a = 100 g cm⁻² (MS_Z0.1 and MS_0.02). With a decrease in metallicity, the duration of the outbursting phase decreased, while the individual luminosity bursts became more powerful. The correlation of the accretion bursts with an increase in inner disk temperature above the critical value indicates that these are MRI outbursts. The disk mass for the lower-metallicity model was marginally larger at the end of the simulations, with the innermost regions (<10 au) showing a larger difference.

The comparison of MS_Σa_Z0.1 and MS_Σa_0.02 with their counterparts with canonical active layer thickness shows similar evolution. The burst phase ended at about 0.42 Myr for MS_Σa_Z0.1, which was almost equal to that of MS_Z0.1. For MS_Σa_0.02, the burst phase lasted 0.22 Myr, which is about 0.05 Myr shorter as compared with MS_Z0.02. Thus, the doubling in the active layer thickness from Σ_a = 100 g cm⁻² had only a marginal impact on the burst phase. With an increased active layer thickness, the overall bottleneck to angular momentum transport in the inner disk was reduced, and the disk could accrete gas more efficiently. This can be inferred by comparing that the final stellar mass in MS_Σa_Z0.1 and MS_Σa_0.02 was larger as compared to their counterparts with canonical active layer thickness. Most notably, the mass contained within the innermost 10 au was substantially smaller in the case of models with an increased Σ_a, indicating a large impact on the structure of the inner disk. Thus, we conclude that although the increased Σ_a at low metallicity has large effects on the structure of the innermost regions of a PPD, the episodic accretion is only marginally affected.

Consider the evolution of the MS_Tc_Z0.02 simulation in Figure 5. The accelerated evolution due to the increased cloud core temperature can be inferred from the cloud infall rate in the first panel. The embedded phase for this simulation lasted only 0.113 Myr, as compared to 0.242 Myr for the rest of the solar-mass models. This model showed the most extreme shortening of the burst phase, which lasted only 0.13 Myr. Note that the disk formation also occurred at an earlier time. In addition to the shortening, the mass infall rate and the base mass accretion rate were larger by a factor of a few as compared to rest of the low cloud core temperature models during the early times. The enhanced accretion rate was responsible for a larger net accretion of material onto the central protostar at the end of the simulation. As a result, the final stellar mass was larger in this case as compared to the MS_0.02 simulation, despite the similar bottleneck of the dead zone. Near the end of the simulation, the total disk mass was smaller in MS_Tc_Z0.02 as compared to MS_0.02; however, the mass contained within the inner 10 au was similar. Thus, a substantial fraction of the disk mass was contained within the innermost regions.

3.4. Lower Stellar Mass

In this section, we present the results for the lower-mass (ML) models with an initial core gas mass of 0.576 M_⊙, which is half as massive as that of the solar-mass models presented earlier. Such a system would ultimately produce a dwarf star at the boundary of the K and M spectral types, with a fully convective interior. Note that the ratio of kinetic to gravitational energy is kept constant across all of the models. The mass of the cloud core limits the total mass available for disk formation, and a proportionally lower-mass disk is seen at all stages of disk evolution. Lower-mass disks are usually less GI unstable, showing suppression in the associated accretion rate variability. The accretion via the disk also lasts for a shorter duration in such systems, primarily due to the limited size of the mass reservoir (Vorobyov 2010).

In Figure 6, we compare the inner disk structure for the three low-mass models: ML_fid, ML_Z0.1, and ML_Tc_Z0.02. Consider the fiducial low-mass model ML_fid. As compared to the corresponding solar-mass model, MS_fid, several key differences can be noticed. The most significant difference in terms of the disk structure is that the gaseous rings, as well as the dead zone, were less robust and sustained for a much shorter time. After the disk formation at 0.04 Myr, the gas surface density could be maintained above Σ_a because of the mass loading from the collapsing core. This explains the formation of the dead zone and gaseous rings at early times. With the limited mass contained in the cloud core, the accreted mass was not replenished from the outer disk, and the layered structure could not be maintained for an extended period of time. The temperature of the innermost part remained relatively high as long as the mass transfer rate was maintained, and then the disk cooled with time.

Consider the effects of reduced metallicity on the inner disk structure in Figure 6. Note that in simulation ML_Z0.1, only the effects of opacity are considered. The trends across the simulations are similar to those described earlier with the solar-mass models. The outer parts of the disks showed similar behavior, while most differences are confined to the innermost regions. Comparing the spacetime plots of the gas surface density and α_eff, the gas accumulation increased with decreasing metallicity. In addition, the extent of the dead zone increased both in radial direction and in time. As described in Section 3.2, the lower gas temperature was caused by more efficient disk cooling in a low-metallicity environment. This, in turn, decreased the kinematic viscosity, contributing to the accumulation of material in the inner disk region. The inner edge of the dead zone also moved inward due to the lower disk temperatures.

In the last model, ML_Tc_Z0.02, both the effects of opacity and an increased cloud core temperature were taken into account. In this case, we can see an accelerated evolution, with the disk forming at 0.02 Myr, as opposed to 0.04 Myr for the rest of the lower-mass models. This is expected, considering the increased infall velocities resulting from the higher cloud core temperature (see Section 3.3). The dead zone formed in ML_Tc_Z0.02 was more robust compared to the rest of the models, as inferred from the α_eff = 10⁻⁴ contour lines. Again, the progressively lower gas temperature as seen in the bottom row is an indicator of the efficient disk cooling due to reduced opacity.

We now focus on the accretion activity in the lower core mass models. In Figure 7, we compare the time-dependent system properties for the three lower-mass simulations. Consider the fiducial model ML_fid. Here the general features of the mass accretion rate are similar to the higher-mass counterpart but limited to a shorter timescale. The embedded phase in this case lasted about 0.095 Myr, as derived from the mass infall rate. The initial phases are dominated by high accretion variability, partly due to the action of GI spirals and partly because of the excursion of the inner edge of the dead zone across the computational domain. In Section 3.2, we explained that the inner boundary of the computational domain needs to be small enough to capture the burst phenomenon occurring at the scale of the innermost 1–2 au of the disk. However, the disk also becomes fully MRI active in the innermost regions due to thermal ionization, where the disk temperature increases above T_crit. This can be seen in the third row, where the midplane temperature at the inner boundary frequently crossed the critical threshold of 1300 K. For model ML_fid, it was challenging to disentangle the spurious variability due to the interaction of the inner edge of the dead zone with the inner boundary of the disk. However, a few bona fide MRI outbursts, as elaborated later in Section 3.5, did occur during the early time. Thus, the MRI bursts were relatively rare in the case of lower-mass disks, while the burst phase also lasted for a shorter time, up to about 0.25 Myr. The evolution of the stellar and disk masses can be observed in the last panel. Near the end of the simulation, the mass of the protostar approached 0.45 M_⊙, while the disk mass continued to decline below 10% of the stellar mass.

Consider the accretion activity of the low-metallicity, lower-mass simulations—ML_Z0.02 and ML_Tc_Z0.02—presented in Figure 7. The trends in accretion are similar to those discussed earlier for the solar-mass counterparts. In the lower-mass case, the period of accretion variability progressively got shorter with lower metallicity. The accelerated evolution for ML_Tc_Z0.02 can be inferred from its higher magnitude and shorter-duration mass infall rate, as compared to the lower-temperature counterparts. The inner boundary temperature in ML_Z0.1 exceeded T_crit several times and triggered MRI outbursts, which can be noticed in concurrent changes in the mass accretion rate and luminosity. The duration of the burst phase in this case was about 0.18 Myr. However, the inner boundary temperature in model ML_Tc_Z0.02 did not increase above the critical threshold needed for triggering MRI outbursts, and as a consequence, no luminosity bursts were observed. Thus, we conclude that the episodic accretion with the mechanism of MRI bursts was rare in low-mass stars as compared to their solar-mass counterparts. In the case of metal-depleted lower-mass stars, such high-luminosity outbursts may be even more uncommon, or completely absent, if the increased cloud core temperature is taken into account.

3.5. Individual MRI Outburst

In this section, we present the results of an individual representative example of an MRI outburst in a low-metallicity environment. For direct comparison, the analysis is similar to that performed in Kadam et al. (2020) for a solar-metallicity model. We restarted the simulation MS_Z0.1 from a suitable checkpoint at about 0.27 Myr and obtained data outputs of 2D fields with a higher time resolution of 20 yr. The particular representative luminosity burst is highlighted in Figure 3 with an arrow. Figure 8 shows some of the relevant time-dependent quantities over a period of 3000 yr in the vicinity of this outburst. Such an MRI outburst empties the inner disk region and causes discontinuity in the spacetime diagram of the surface density, as can be observed in the top row of Figure 2. The vertical dashed lines in Figure 8 correspond to the 2D snapshots analyzed in Figure 9. The first panel in Figure 8 shows the highly asymmetrical light-curve profile of the outburst. The maximum luminosity reached about 118 L_⊙, and the total duration was about 200 yr. The initial rise time of the outburst was relatively fast, about 9 yr, while the decline after the maximum luminosity was much slower, lasting about 175 yr. The duration of the outburst is somewhat smaller than the associated viscous timescale of about 1000 yr, which indicates that the gravitational torques play an important role due to the transient nonaxisymmetry of the disk. The mass accretion rate, as shown in the second panel, increased from about 10⁻⁶ to about 4 × 10⁻⁵ M_⊙ yr⁻¹. The accretion rate also showed an asymmetrical, sharp increase and a slower decline. The stellar mass as shown in the next panel increased by about 0.008 M_⊙ over the duration of the outburst, which was about 1.5% of the mass of the protostar at the time. The last three panels show the conditions at the inner boundary across the outburst in terms of azimuthally averaged gas surface density, α_eff, and midplane temperature. The sudden increase in Σ_inner by over an order of magnitude indicates the flow of material across the boundary that was accreted onto the central protostar producing the powerful eruption. The fundamental reason behind the eruption was the increase of the midplane temperature above the critical value of 1300 K, when the MRI in the disk was triggered (hence the name "MRI outburst"). The MRI activation can be seen in the plot of α_inner, where it achieves a maximum value of 0.01. The postburst values of Σ_inner, as well as T_inner, were considerably lower as compared to those before the eruption, meaning that the inner disk depleted its mass and became colder. Since the inner disk was emptied of a large amount of material, α_inner increased after the outburst until the gas could be replenished from the outer regions (see Equation (9)).

**Figure 8.** Time-dependent system properties (total luminosity, ${\dot{M}}_{\mathrm{mean}}$ , stellar mass), as well as surface density, α_eff, and midplane temperature at the inner boundary during a typical MRI outburst occurring in low-metallicity simulation `MS`_`Z0.1`. The vertical dashed lines mark the time corresponding to the 2D plots presented in Figure 9. The red dashed lines in the last two panels show α_bg = 10⁻² and T_crit = 1300 K, respectively.
Download figure:
Standard image High-resolution image

**Figure 8.** Time-dependent system properties (total luminosity, ${\dot{M}}_{\mathrm{mean}}$ , stellar mass), as well as surface density, α_eff, and midplane temperature at the inner boundary during a typical MRI outburst occurring in low-metallicity simulation `MS`_`Z0.1`. The vertical dashed lines mark the time corresponding to the 2D plots presented in Figure 9. The red dashed lines in the last two panels show α_bg = 10⁻² and T_crit = 1300 K, respectively.
Download figure:
Standard image High-resolution image

Figure 9 depicts the progression of the outburst in terms of 2D distributions of gas surface density, midplane temperature, and effective α. These snapshots (or frames) are marked with vertical dashed lines in Figure 8 and show the innermost 10 × 10 au² region of the disk. These particular frames were chosen such that the disk structure could be captured during the rise, at the maximum luminosity, and during the decline, as well as before and after the outburst. The first column shows a quiescent disk before an MRI outburst where the gaseous ring was accumulated at about 1 au from the central protostar. The second column captures the disk during the sharp increase in mass transfer rate at the beginning of the outburst. The midplane temperature in the innermost part increased above the critical value, while α_eff reached a maximum value of a fully active disk, indicating MRI activation. The third column shows the disk structure near the maximum of the luminosity burst and also corresponds to the maximum extent of the MRI-activated region, up to about 1.5 au. The MRI-activated region engulfed the former gas ring, which can be seen almost fully disrupted in this frame. The accumulated material accreted through the inner boundary due to the increased viscosity, resulting in the large mass transfer rate. It is clear that this was an inside-out MRI outburst, originating near the inner boundary due to viscous heating and propagating in the outward direction in a snowplow fashion. The fourth column shows the disk structure during the decline. Only a small region where the temperature remained above the critical value remained fully MRI active at this point. Note that the disk structure showed significant transient asymmetries during the outburst. The fifth column shows a quiescent disk shortly after the outburst. As the viscosity channels the bulk of the material inward, the angular momentum is transported outward via a relatively small amount of mass. This can explain the final low surface density ring surrounding the maximum extent of the MRI-active region.

The general properties, as well as the mechanism of an individual MRI outburst in a low-metallicity environment, were same as those in a solar-metallicity disk (Kadam et al. 2020). The MRI was triggered at the inner boundary due to an increase in temperature above the critical value due to viscous heating in the dead zone, resulting in an inside-out burst. However, there were a few key differences. When compared to their solar-metallicity counterparts, the outbursts in low-metallicity disks were more powerful (in this particular example, 118 versus 75 L_⊙) and showed a sharply rising profile (10 versus 100 yr) with a shorter total duration (200 versus 300 yr). As elaborated in Section 3.2, the structure of the innermost region was altered with the lowering of the metallicity. A larger amount of material was accumulated in the gaseous rings, while this material also moved closer to the protostar (see Figure 2). As the outburst progressed, the MRI-activated front encountered this mass reservoir within a shorter time as compared to a solar-metallicity disk. This explains the steeper rise time of the outburst, as well as the larger peak mass transfer rate, making the outburst more luminous. The total duration of the outburst is proportional to the viscous timescale, which was reduced at this shorter radius. Note that the total mass accreted in the low-metallicity outburst was also larger, approximately 0.008 M_⊙, as compared to 0.005 M_⊙ in the solar-metallicity case (Kadam et al. 2020). In Figure 3, we see a continued trend of increased luminosity of individual bursts with a further decrease in metallicity, which can be explained along the same line of reasoning.

4. Conclusions

In this study, we conducted numerical experiments using global hydrodynamic simulations of PPD formation in the thin-disk limit, with the aim of understanding the effects of metal-poor environments—10% and 2% of the solar metallicity—on the young disks. The primary focus was on the episodic accretion occurring during these early stages of evolution, as well as the inner disk structure. The simulations started with the collapse phase of the initial cloud core, and the inner boundary was placed at 0.42 au so that the MRI bursts occurring at the sub–astronomical unit scale could be captured. The dead zone was modeled by an adaptive and effective α, which depended on the disk gas surface density, as well as midplane temperature. The low metallicity was modeled by scaling down the gas and dust opacities. The effects of the increased thickness of the magnetically active surface layer and the temperature of the prestellar cloud core were also studied. At low metallicities, the general structure and behavior of the accreting disks were similar to those at the solar value. The dead zone was formed in the inner approximately 10 au region, with gaseous ringlike structures appearing near the inner edge of the dead zone. The early evolution was dominated by formation of spirals and GI activity, resulting in variability in the mass accretion rate. The disks around solar-mass stars showed episodic accretion via MRI bursts for a limited time during the evolution of the disk. We emphasize that the span of the considered spatial scales (0.4–10,000 au) and the evolution period (0.7 Myr) stretch the boundaries of the thin-disk models to their limit. A simulation in this study typically takes more than a month of computing time on 32 CPU cores, and such calculations are prohibitively expensive in full 3D models.

The salient differences between the solar and lower-metallicity PPD formation simulations are summarized as follows.

1.
With lower metallicity, the ringlike structures that formed in the dead zone were more robust in terms of accumulated gas, as well as low effective α parameter, while the accumulation occurred closer to the central protostar. These effects were due to the effective cooling of the disk at low metallicity, which also resulted in low midplane temperatures in general.
2.
The mass accumulated in the inner regions of metal-poor disks was significantly larger as compared to their solar-metallicity counterparts. This can explain some of the observed trends, e.g., the higher mass accretion rates for low-metallicity YSOs at later stages (Spezzi et al. 2012) and insensitivity of the frequency of short-period super-Earths to the metallicity (Petigura et al. 2018).
3.
The burst phase, when the disk is subject to powerful MRI outbursts, became shorter with decreasing metallicity. This was especially true for the lowest-metallicity models if the increased molecular cloud core temperature due to inefficient line cooling was considered. Lower inner disk temperatures in metal-poor disks also resulted in shortening of the burst phase.
4.
The MRI outbursts were rare for a lower-mass star of a main-sequence mass of ⪅0.5 M_⊙, and they were even more unlikely for their low-metallicity counterparts, when the increased cloud core temperature was included.
5.
An individual MRI outburst in the low-metallicity disk was shorter in duration, more luminous, showed a steep rising luminosity curve, and accreted more mass, as compared to its solar-metallicity counterpart.

Although low-metallicity PPDs are currently difficult to observe, our predictions may be verified in the future by high-resolution infrared spectroscopy and imaging with the next generation of large telescopes, e.g., E-ELT and TMT.

Here we mention some of the limitations of this investigation. The critical temperature for MRI ignition and the thickness of the active layer were considered constant in this study. Magnetohydrodynamics equations that incorporate a detailed ionization balance need to be solved in order to obtain a better estimate of the disk behavior. The evolution of the dust component was not considered in this study, which can alter the disk thermodynamics, especially because of its accumulation in the regions of the pressure maxima. The innermost fully MRI-active region (formed due to collisional ionization) could interfere with the inner boundary of the computational domain, producing spurious variations in the mass accretion rate. This remains a general issue with simulations that terminate at a certain distance from the central protostar. The region between the inner boundary and stellar surface is complex and includes several important phenomena, such as the inner rim of the dust and magnetospheric accretion. Despite these uncertainties, the overall picture and trends presented in this study with respect to decreasing opacity should remain valid.

We thank the anonymous referee for constructive comments that improved the quality of the manuscript. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 716155 (SACCRED). E.V. acknowledges support from the Austrian Science Fund (FWF) under research grant P31635-N27. The simulations were performed on the Vienna Scientific Cluster (VSC-3 and VSC-4).

Appendix: Comparison across Spatial Resolution

When conducting numerical simulations, one usually needs to find a balance between computational cost and a grid resolution sufficient to capture the phenomena of interest. The implementation of the hydrodynamic equations in FEoSaD has been carefully evaluated with several test cases, confirming the ability of our numerical schemes to reproduce the known analytic solutions in the thin-disk limit (Vorobyov & Basu 2006, 2010). As the actual problem of PPD formation and its dynamical evolution is fundamentally nonlinear in nature, the adequate grid resolution can only be found by trial and error. In this appendix, we present the results of a comparative analysis of two identical simulations conducted at different resolutions. The fiducial solar-mass model (MS_fid) is compared with its high-resolution counterpart, model1_T1300_S100 in Kadam et al. (2019). Henceforth, the latter model will be termed MS_fid_HR for convenience. The simulations presented in this study, including MS_fid, were all conducted at a resolution of 256 × 256, while the MS_fid_HR model corresponds to twice this resolution, i.e., 512 × 512. We show that the results of these simulations are congruent, and the lower resolution used for the simulations presented in this study was sufficient for capturing the essential aspects of the evolution of the PPDs.

Even with the advantages of the thin-disk approximation in our calculations, increasing the grid resolution in FEoSaD is computationally demanding for two reasons. Doubling the resolution in both the R- and ϕ-direction increases the number of cells by a factor of 4. In addition, because of the cylindrical geometry of the grid, the Courant et al. (1928) condition near the axis of the grid implies more stringent constraints for the numerical stability of the code. This forces a smaller time step for the hydrodynamical calculations at higher resolution, making the final computation up to 16 times more expensive. As a result, a high-resolution simulation such as MS_fid_HR can run on a computer node with 48 cores (FEoSaD is OpenMP-parallelized) for several months of real time. This longer run time of the high-resolution models was the primary reason for the chosen grid resolution for this study, wherein the entire burst phase needed to be captured.

Figure 10 compares the inner disk structure of model MS_fid with its high-resolution counterpart. The overall evolution of the disk structure was similar with respect to all three quantities: Σ, α_eff, and T_mp. The disks formed at the same time in both models at about 0.08 Myr. The outer boundary of the disk (Σ = 1 g cm⁻² contours), as well as the dead zone (black contours in α_eff), showed almost identical radial extents at all times. The rings formed in the inner disk were similar in strength in terms of both surface density and viscosity, and the latter can be inferred from the α_eff = 10⁻⁴ contours. As elaborated earlier, the discontinuities in the rings typically correspond to MRI bursts, wherein the disk is rapidly accreted through the inner boundary. The evolution of the disk midplane temperature was also consistent across the two simulations. Two key differences between the disk structure were observed in the innermost regions—the rings in the ML_fid_HR simulation showed much more intricate features, while the temperature crossed the critical threshold less frequently. The radial extent of the computational domain in both models was identical, extending from 0.4167 to 10,210 au. The highest grid resolution in the radial direction was 1.6 × 10⁻² au at the innermost cell of the disk for MS_fid, as compared to 8.2 × 10⁻³ au for the MS_fid_HR model. Thus, the high-resolution model was able to spatially resolve the inner disk much better, and the less frequent temperature excursions in the MS_fid_HR model may also be a result of the improved resolution. However, the general picture with respect to the disk structure and its evolution remained unchanged for the two models.

In Figure 11, we compare the low- and high-resolution models with respect to time-dependent properties relevant to episodic accretion. The overall behavior for all four panels was very similar. Note that the simulations are compared up to 0.42 Myr of evolution and not during the entire burst phase. As explained in Section 3.2, the mass accretion rate was a superposition of the three modes of variability: excursion of the inner edge of the dead zone across the computational boundary, MRI bursts, and vigorous GI activity. The temperature at the inner boundary crossed the MRI activation threshold, with the larger-amplitude changes corresponding to the MRI outbursts. The evolution of the stellar and disk masses was almost identical across the resolutions for the two models. One notable difference between the two simulations can be seen in the behavior of the mass accretion rate and hence the total luminosity. In between the larger MRI bursts, the background variability was diminished with the increase in resolution. This was also reflected in the evolution of T_inner, where it crossed the critical threshold of 1300 K less frequently in between the larger-amplitude MRI outbursts for the MS_fid_HR model. As part of the small-scale variability was caused by the movement of the inner edge of the dead zone across the computational boundary, it may be improved with the better-resolved innermost regions. The overall evolution remained congruent for the two simulations, and we conclude that the resolution chosen for this study was sufficient to capture both the inner disk structure and the outbursting behavior of the systems.

Eruptive Behavior of Magnetically Layered Protoplanetary Disks in Low-metallicity Environments

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Model Description and Initial Conditions