© ESO 2020

1. Introduction

Circumstellar disks are a key element in the theory of star and planet formation. They form during the gravitational collapse of pre-stellar cloud cores and disperse as the host stars evolve toward the main sequence. Mass transport in circumstellar disks determines the terminal stellar mass, and dust growth in the disk sets the stage for planet formation.

Among many processes that operate in circumstellar disks, gravitational instability (GI) is of particular interest in the context of protostellar accretion and planet formation. GI is a dominant mechanism of mass transport in the early protostellar stage of disk evolution (Lin & Pringle 1990; Vorobyov & Basu 2009), thus largely setting the rate of mass accretion on the star. Moreover, the character of mass accretion versus time can also be affected by GI, either directly through perturbations caused by global spiral arms (Elbakyan et al. 2016) or indirectly by assisting the magnetorotational instability in the disk innermost regions (e.g., Zhu et al. 2009; Bae et al. 2014). The role of GI in disk evolution is versatile and it can even affect the chemical composition of circumstellar disks (Evans et al. 2015).

The extreme manifestation of GI, disk gravitational fragmentation, is another important phenomenon that can affect the character of mass accretion on the growing protostar. Inward-migrating clumps formed via disk fragmentation bring a large amount of mass to the disk’s inner regions, triggering mass accretion bursts when destroyed by tidal torques (Vorobyov & Basu 2005, 2015; Meyer et al. 2017). Disk fragmentation is also considered a viable and perhaps dominant scenario for the formation of giant planets on wide orbits where disk conditions can be favorable for gravitational fragmentation (Boss 2003; Johnson & Gammie 2003; Mayer et al. 2007; Boley 2009; Vorobyov 2013). When clump migration and dust settling are taken into consideration, disk fragmentation can also account for the formation of icy and rocky planets at closer orbits (Boss 1998; Nayakshin 2010, 2017; Boley et al. 2010; Vorobyov & Elbakyan 2019). It is not yet clear whether disk fragmentation can explain the majority of the known close-in giant planets, for which core accretion remains a preferable scenario (Pollack et al. 1996; Alibert et al. 2005).

Gravitational instability and disk fragmentation have been extensively studied for solar metallicity disks with respect to mass transport and planet formation. When lower metallicity environments are considered, the importance of GI in setting the mass accretion rate and the role of disk fragmentation as a planet formation mechanism becomes less clear. Observations of accretion in low-metallicity galactic and extragalactic star-forming regions show that accretion rates on the low-metallicity stars are on average higher than in their solar metallicity counterparts (Spezzi et al. 2012; De Marchi et al. 2013). On the other hand, Kalari & Vink (2015) found no difference between the mass accretion rates inferred for the low-metallicity (1/5 Z_⊙) open cluster Dolidze 25 and solar metallicity Galactic pre-main sequence stars in the same mass range. The reason for this difference remains to be understood.

Numerical hydrodynamics simulations of primordial (Z = 0) disks indicate that they are vigorously GI unstable and produce many clumps that quickly migrate to the star, triggering accretion and luminosity outbursts (Vorobyov et al. 2013; Sakurai et al. 2016). As the disk metallicity increases to Z = 10⁻⁵ − 10⁻³ Z_⊙, steady-state models also predict that circumstellar disks should be strongly gravitationally unstable (Tanaka & Omukai 2014). However, at subsolar metallicities of Z ∼ 0.1 − 1.0 Z_⊙ studies of circumstellar disks produce conflicting results. For instance, Cai et al. (2006) found that GI is stronger in lower metallicity disks, but disks nevertheless do not fragment. Conversely, Boss (2002) reported that the propensity of circumstellar disks to gravitational fragmentation is largely insensitive to metallicity variations by a factor of 10 on both sides of the solar value because the disk’s thermal balance is mainly determined by the radiative input from the central star. Furthermore, Meru & Bate (2010) also found that reducing metallicity below the solar value promotes disk fragmentation since lower opacities (associated with lower metallicities) allow the disk to cool more quickly.

In this work, we revisit the issue of protostellar accretion and gravitational fragmentation in disks with subsolar metallicities using numerical hydrodynamics simulations in the thin-disk limit. Unlike previous studies, in this paper we consider a wider range of metallicities from Z = 10⁻² Z_⊙ to Z = 1.0 Z_⊙. In addition, we employ a sophisticated thermal scheme (Vorobyov et al. 2020), which allows decoupling of gas and dust temperatures at low metallicities. Finally, the disk hydrodynamics simulations are linked with a stellar evolution code, which allows us to compute the time- and metallicity-dependent stellar radiative input as the disk evolves with time.

The paper is organized as follows. In Sect. 2 we review the main aspects of disk modeling at metallicities lower than solar. In Sect. 3 we present the initial stages of gravitational collapse for different metallicities. In Sect. 4 we discuss the disk evolution in the context of solar metallicity models. In Sect. 5 we present the results of disk modeling and mass accretion rates at low metallicities. In Sect. 6 we summarize our main results.

2. Model description

We use numerical hydrodynamics simulations in the thin-disk limit to model the evolution of protostellar disks with different metallicities. The model was described in detail in Vorobyov et al. (2020). Here we only review the key aspects of disk modeling at low metallicities.

The numerical model takes disk self-gravity and turbulent viscosity into account. The turbulent viscosity is described via the usual Shakura & Sunyaev parameterization with the α-value in the 10⁻⁴ − 10⁻² limits. Unlike many previous studies of disk evolution that make no distinction between gas and dust temperatures, we employ a sophisticated thermal evolution scheme that separates gas and dust temperatures and considers the following energetic processes: radiative continuum cooling (or heating) of gas and dust (including collisional transfer of energy between gas and dust), H₂ and HD line cooling, chemical cooling and heating through considered chemical reactions, metal line cooling, stellar and background irradiation, and viscous heating. Subsolar metallicities are set by scaling down the gas and dust opacities, dust-to-gas mass ratio, and metal content of the solar metallicity disk by the corresponding factor.

The pertinent hydrodynamic equations for mass, momentum, and internal energy density in the thin-disk limit can then be written as

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Sigma }}{\partial t}=-{\mathbf{\nabla }}_p·\left(\mathrm{\Sigma }{\mathit{v}}_p\right),\end{array}$$(1)

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\mathrm{\Sigma }{\mathit{v}}_p\right)+{[\mathbf{\nabla }·(\mathrm{\Sigma }{\mathit{v}}_p\otimes {\mathit{v}}_p)]}_p=-{\mathbf{\nabla }}_{p}P+\mathrm{\Sigma }{\mathit{g}}_p+{(\mathbf{\nabla }·\mathbf{\Pi })}_p,\end{array}$$(2)

$$\begin{array}{c}\hfill \frac{\partial e}{\partial t}+{\mathbf{\nabla }}_p·\left(e{v}_p\right)=-P({\mathbf{\nabla }}_p·v_p)-Q_{\mathrm{tot}}+{\left(\mathbf{\nabla }v\right)}_{p{p}^{\prime }}:{\mathbf{\Pi }}_{p{p}^{\prime }},\end{array}$$(3)

where the subscripts p and p′ refer to the planar components (r, ϕ) in polar coordinates, Σ is the gas mass surface density, P is the vertically integrated gas pressure calculated via the ideal equation of state as P = (γ − 1)e, e is the internal energy per surface area, γ is the ratio of specific heats, ${\mathit{v}}_p=v_{\mathrm{r}}\widehat{\mathit{r}}+v_{\varphi }\widehat{\mathit{\varphi }}$ is the velocity in the disk plane, ${\mathit{g}}_p=g_{\mathrm{r}}\widehat{\mathit{r}}+g_{\varphi }\widehat{\mathit{\varphi }}$ is the gravitational acceleration in the disk plane (including that of the disk and the star), and ${\mathbf{\nabla }}_p=\widehat{\mathit{r}}\partial /\partial r+\widehat{\mathit{\varphi }}{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk. Turbulent viscosity enters the basic equations via the viscous stress tensor Π. The total cooling–heating rate is expressed as

$$\begin{array}{c}\hfill Q_{\mathrm{tot}}=(Q_{\mathrm{cont}}+Q_{{\mathrm{H}}_2}+Q_{\mathrm{HD}}+Q_{\mathrm{chem}}+Q_{\mathrm{metal}})2H,\end{array}$$(4)

where H is the vertical scale height calculated assuming a local hydrostatic equilibrium in the gravitational field of the star and disk (see Vorobyov & Basu 2009), Q_cont is the rate of radiative continuum cooling (or heating) of gas and dust (including collisional transfer of energy between gas and dust), Q_H2 is the H₂ line cooling rate, Q_HD is the HD line cooling rate, Q_chem is the chemical cooling–heating rate (through chemical reactions), and Q_metal is the metal line cooling rate. All constituents of Q_tot are volumetric cooling or heating rates. The dust temperature is determined in the steady-state limit by the energy balance on dust grains due to the thermal emission, absorption, and collision with gas. We note that cooling–heating described by Eq. (4) does not include the adiabatic cooling–heating and viscous heating, which are taken into account by the first and third terms on the right-hand side of Eq. (3).

In addition to hydrodynamic Eqs. (1)–(3), we solve non-equilibrium kinetic equations for H, H₂, H⁺, D, HD, D⁺, and e, while the H⁻ fraction is calculated from the equilibrium assumption. We assume that helium is always neutral and its fractional abundance is y_He = 8.333 × 10⁻². We further assume that our species are collisionally coupled with gas and solve the continuity equation for each considered species using the same third-order-accurate algorithm as for the gas. A more detailed description of the thermal evolution scheme and the solution procedure can be found in Vorobyov et al. (2020).

The central protostar is not only a source of gravity, but also provides radiative heating through irradiation of the disk surface. Stellar characteristics, such as the radius and photospheric luminosity, are calculated in line with the disk evolution using the stellar evolution code STELLAR (Yorke et al. 2008), modified by Hosokawa et al. (2013) to take different metallicities into account. The stellar characteristics are further used to calculate the total luminosity (stellar photospheric plus accretion) and the radiation flux impinging on the surface of the disk. The co-evolution of the disk and the central protostar allows us to accurately follow the evolution of a highly dynamical system with variable accretion and luminosity.

The starting point of numerical simulations is the gravitational collapse of a pre-stellar core. We consider three model realizations with distinct metallicities: Z = 0.01 Z_⊙, Z = 0.1 Z_⊙, and Z = 1.0 Z_⊙. In addition, models with different α-parameters describing the strength of turbulent viscosity in the disk are considered. The initial radial profiles of the gas surface density Σ and angular velocity Ω of the pre-stellar core have the form typical of rotating cloud cores formed via slow expulsion of magnetic field due to ambipolar diffusion, with the angular momentum remaining constant during axially symmetric core compression (Basu 1997):

$$\begin{array}{c}\hfill \mathrm{\Sigma }=\frac{r_0{\mathrm{\Sigma }}_0}{\sqrt{r^2+r_0^2}},\end{array}$$(5)

$$\begin{array}{c}\hfill \mathrm{\Omega }=2{\mathrm{\Omega }}_0{\left(\frac{r_0}{r}\right)}^2[\sqrt{1+{\left(\frac{r}{r_0}\right)}^2}-1].\end{array}$$(6)

Here Σ₀ and Ω₀ are the gas surface density and the angular velocity at the center of the core. The radius of the central plateau is proportional to the thermal Jeans length and is defined as $r_0=A{c}_{\mathrm{s}}/\sqrt{\pi G{\rho }_0}$, where ρ₀ is the central gas volume density, $c_{\mathrm{s}}=\sqrt{\mathcal{R}{T}_{\mathrm{g},0}/\mu }$ is the initial sound speed in the core, μ is the mean molecular weight, T_g, 0 is the initial gas temperature, ℛ is the universal gas constant, and A is the amplitude of the initial positive density perturbation that drives the system out of equilibrium and trigger gravitational collapse.

To define the initial surface densities, temperatures, and chemical fractions of pre-stellar cores at distinct metallicities, we use a thermal model that is similar to that presented by Omukai et al. (2005) for collapsing cores in the one-zone approximation. Their Fig. 1 presents the gas temperature (T_g) evolution of pre-stellar clouds with different metallicities as a function of the gas number density n in the center of the core. For Z = 1.0 Z_⊙ and Z = 0.1 Z_⊙ we choose n = 10⁶ cm⁻³, while for Z = 0.01 Z_⊙ we use n = 2 × 10⁷ cm⁻³. A higher value of n for the latter case is explained by the fact that at lower densities the T_g ∝ n^γ relation is steeper (i.e., γ is higher) than what is allowed for gravitational collapse to occur (γ ≤ 4/3). At lower densities the Z = 0.01 Z_⊙ core would not collapse without strong external compression, which is not considered in our model.

Fig. 1. Radial profiles of the gas surface density (top left), radial infall velocity (top right), gas temperature (bottom left), and dust temperature (bottom right) in model Z1.0. The dashed curves correspond to the initial setup at the onset of gravitational collapse, while the solid curves show the quantities at the time instance immediately preceding disk formation. Curves of different color correspond to models with different metallicities, as shown in the legend.

Once the central number density is set, the corresponding gas and dust temperatures are obtained, and the resulting r₀ and Σ₀ = μm_Hnr₀ are used to define the radial surface density profile of the core, where m_H is the mass of hydrogen. We note that Omukai et al. (2005) neglected the external irradiation heating other than that provided by the cosmic microwave background. Here we also take external stellar irradiation into account with a constant temperature of 10 K. This explains why our solar metallicity core has a higher temperature than the value obtained in Omukai et al. (2005). The outer radius of the core is varied to create objects of different mass. The initial angular velocity profile is set by varying Ω₀ until the desired ratio of rotational to gravitational energy β is obtained.

The initial parameters of the models are detailed in Table 1. We consider two models with the solar metallicity and eight additional models with lower metallicities. The models are dubbed according to their ordinal number and metallicity. For instance, model 2_Z0.1 corresponds to model 2 with metallicity Z = 0.1 Z_⊙. The initial core masses of models with low metallicity vary by a factor of three to produce disks with different masses. The ratio of rotational to gravitational energy β is similar in all models to exclude its effect on the disk evolution. At the same time, β is sufficient to produce gravitationally unstable disks in the solar metallicity case (Vorobyov 2013).

The numerical resolution is 512  ×  512 grid zones for all models except for models 4_Z0.1 and 4_Z0.01 with smallest pre-stellar cores, for which a twice coarser resolution was used. The polar grid (r, ϕ) is logarithmically spaced in the radial direction and has equal spacing in the azimuthal direction. The inner 5 au were cut out and replaced with the sink cell. In this work the mass accretion rates are calculated as the gas mass that crosses the sink–disk interface per unit time. We note that the actual mass accretion rate on the star may be modified by the physical conditions and mechanisms that may operate in the innermost disk regions (such as the magnetorotational instability). The radial size of the innermost grid cell (just outside 5 au) varies in the 0.075–0.14 au limits, depending on the model, and it scales linearly with the radial distance. The Truelove criterion states that the local Jeans length must be resolved by at least four numerical cells to correctly capture disk fragmentation (Truelove et al. 1998). In the thin-disk limit the Jeans length can be expressed as (Vorobyov 2013)

$$\begin{array}{c}\hfill R_{\mathrm{J}}=\frac{c_{\mathrm{s}}^2}{\pi G\mathrm{\Sigma }}.\end{array}$$(7)

Fragments usually condense out of the densest sections of spiral arms at a typical distance of 100 au and then either migrate inward or scatter outward. The typical surface densities and temperatures in spiral arms do not exceed 100 g cm⁻² and 100 K. Adopting these values, the corresponding Jeans length is R_J≈ 20 au. The numerical resolution at 100 au varies in the 1.5–2.8 au limits, thus fulfilling the Truelove criterion. Finally, we note that the inner boundary condition allows matter to flow in both directions through the sink-disk interface, thus greatly reducing a spurious drop in the gas density near the sink that may develop in models with a one-way boundary condition (for details, see Vorobyov et al. 2018).

3. Pre-disk stage

In this section we focus on the initial stages of core collapse for different metallicities. Figure 1 depicts the radial distributions of gas surface density, gas and dust temperatures, and radial infall velocity at the onset of collapse and at the time instance immediately preceding disk formation. Three models with similar core masses but distinct metallicities are shown: model 1_Z1.0, model 2_Z0.1, and model 2_Z0.01. The initial state and the pre-disk state are also shown.

There are several key differences in the pre-stellar evolution of the Z = 0.1 − 1.0 Z_⊙ models, on the one hand, and Z = 0.01 Z_⊙ model, on the other hand. First, the gas and dust temperatures in the latter model are systematically higher than in the former models. At the considered metallicities and densities, the dust continuum emission is a dominant cooling mechanism (Omukai et al. 2005) and the reduction in the dust content leads to an increase in the gas temperature during the collapse phase.

Furthermore, the dust temperature in model 2_Z0.01 decouples notably from that of gas. In particular, in the initial state the decoupling is seen throughout the entire extent of the core, while in the pre-disk phase the decoupling is limited to the outer low-density regions. The effect is also noticeable at higher metallicities, but is strongest for Z = 0.01 Z_⊙. Assuming that the dust radiative cooling balances heating due to collisions with gas, the dust temperature can be expressed as (Vorobyov et al. 2020)

$$\begin{array}{c}\hfill T_{\mathrm{d}}\simeq 120\mathrm{K}{\left(\frac{T_{\mathrm{g}}}{100\mathrm{K}}\right)}^{0.3}{\left(\frac{n}{{10}^{10}{\mathrm{cm}}^{-3}}\right)}^{0.2}.\end{array}$$(8)

By setting T_g = T_d we define the threshold temperature above which gas and dust thermally decouple from each other. This threshold temperature can be written as

$$\begin{array}{c}\hfill T_{\mathrm{crit}}\simeq 130\mathrm{K}{\left(\frac{n}{{10}^{10}{\mathrm{cm}}^{-3}}\right)}^{0.3}.\end{array}$$(9)

Consider now model 2_Z0.01 with n = 2 × 10⁷ cm⁻³ in the initial state. The corresponding T_g = 30.3 K and Eq. (9) yields T_crit = 20 K, meaning that the gas and dust temperatures decouple from each other. This is indeed the case, as is seen in Fig. 1. For model 1_Z1.0 with n = 10⁶ cm⁻³ in the initial state the critical temperature is T_crit = 8.2 K and the gas temperature (set by external stellar irradiation) is 10 K, meaning that decoupling is weak. In the pre-disk state the trend is similar: lower metallicity models show stronger decoupling between gas and dust temperatures.

Finally, we note that the infall velocity v_r in model 2_Z0.01 is appreciably higher (by the absolute value) than in the other two models. This is in agreement with analytical arguments stating that the mass infall rate is related to the gas temperature as (Shu 1977)

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{infall}}=-2\pi r\mathrm{\Sigma }{v}_{\mathrm{r}}\simeq \frac{T_{\mathrm{g}}^{3/2}}{G}.\end{array}$$(10)

As we show later in the text, higher infall rates and temperature decoupling have important consequences for the disk evolution in the lowest metallicity models.

4. Disk evolution and accretion bursts at solar metallicity

In this section we review the disk evolution and protostellar accretion in the solar metallicity case using model 1_Z1.0 as a prototype example. This will help us later to understand the corresponding processes in low-metallicity models. Figure 2 presents the gas surface density distribution in the inner 1200 × 1200 au² box at six time instances after disk formation. Already after 20 kyr the disk becomes gravitationally unstable, as indicated by the presence of two spiral arms. At this stage the disk is nevertheless too small to fragment. Eighty thousand years later the disk is vigorously unstable and hosts several dense gaseous clumps formed via gravitational fragmentation. Disk fragmentation continues to the end of our numerical simulations for at least 400 kyr. We note that the number of fragments changes with time; the increase means disk fragmentation and decrease implies clump destruction. The properties of the formed clumps will be considered in a follow-up study.

Fig. 2. Disk evolution in model 1_Z1.0. Shown is the gas surface density at six evolutionary times since the disk formation instance. The time is counted starting from the instance of disk formation, which occurs at 89 kyr in absolute time elapsed since the onset of gravitational collapse. The scale bar is in g cm−2.

In this analysis we focus on the inward migration of the clumps and on the accretion bursts that follow as the clumps disperse in the innermost disk regions. Figure 3 presents the disk evolution over a short period of time focused on one episode of clump migration at t ≈ 240 kyr. Shown is the inner 400  ×  400 au² box at six times spanning a range of 100 yr. The inward migrating clump is highlighted with the arrows. Initially, the clump is found at a radial distance of 54 au and after 80 yr it moved to 25 au. The quick inward migration is caused by the gravitational interaction between the clumps, as was described in detail in Vorobyov & Elbakyan (2018). As the clump continues its migration toward the star its Hill radius becomes smaller than its radius and the clump starts losing its material. The Hill radius is defined as

$$\begin{array}{c}\hfill R_{\mathrm{H}}=r_{\mathrm{cl}}{\left(\frac{1}{3}\frac{M_{\mathrm{cl}}}{M_{\ast }+M_{\mathrm{cl}}}\right)}^{1/3},\end{array}$$(11)

Fig. 3. Disk evolution in model Z1.0. Shown is the gas surface density at six evolutionary times encompassing a short time period of clump migration at t ≈ 240 kyr. The scale bar is in g cm−2.

where M_cl is the clump mass, M_* is the stellar mass, and r_cl is the radial position of the clump. At t = 80 yr the Hill radius is 4.5 au, while the radius of the clump is approximately 8 au, meaning that the clump has already begun to disintegrate. The hollow structure at about the position expected for the clump at 90 and 100 yrs is related to this burst-like disintegration, perhaps caused by tidal heating of the fast-rotating clump. Part of the released material is accreted through the central sink cell causing an accretion burst. The process of clump migration repeats as long as the disk is massive enough to sustain continual disk fragmentation. This usually occurs in the embedded stage of disk evolution (lasting for several hundred kyr) when mass-loading from the infalling envelope replenishes the disk mass loss via accretion.

The black line in Fig. 4 presents the mass accretion rate as a function of time in model 1_Z1.0, which is calculated as the mass passing through the sink cell per unit time. We note that the radius of the sink (5 au) is much greater than the stellar radius, meaning that the calculated accretion rates may be affected by the disk physics interior to 5 au. For instance, the magnetorotational instability may introduce additional accretion bursts (Zhu et al. 2009) or the infalling clumps themselves may act as a trigger for the magnetorotational instability (Ohtani et al. 2014). The red line shows the mass infall rate on the disk from the envelope. This quantity is calculated as the mass flux at a radial distance of 1000 au from the star. The mass accretion rate is both time-declining and highly variable. High accretion variability reflects a highly dynamical nature of disk evolution. Long-term variations are caused by secular changes in the spiral pattern (Elbakyan et al. 2016), while the high-magnitude short-term spikes are triggered by infalling clumps. The strongest burst at ≈240 kyr corresponds to the clump infall depicted in Fig. 3. The strength of accretion variability diminishes as the mass infall rate on the disk declines with time. Reduced disk mass-loading cannot anymore sustain gravitational instability and fragmentation in the disk, which is the main driving force for accretion variability in our model. The reader is referred to Vorobyov & Basu (2015) for the in-depth discussion of the accretion bursts caused by inward-migrating clumps in the gravitationally unstable solar metallicity disks. It is also worth mentioning that some clumps may not be affected by inward migration and tidal destruction, thus forming giant planets or brown dwarfs on orbits from tens to hundreds of astronomical units (e.g., Cha & Nayakshin 2011; Vorobyov 2013; Stamatellos 2015; Vorobyov & Elbakyan 2018). This is consistent with the detection of giant planets and brown dwarfs on wide orbits via direct imaging (e.g., Marois et al. 2008; Lafrenière et al. 2010; Chauvin et al. 2017).

Fig. 4. Mass accretion rate vs. time in model 1_Z1.0 shown with the black line. The red line presents the mass infall rate on the disk. The time is counted from the onset of gravitational collapse and is shown on the log scale to better resolve the variability at early times. The disk forms at 89 kyr.

5. Disk evolution and accretion burst at lower metallicities

We now discuss the disk evolution in models with metallicities lower than solar. We start with the Z = 0.1 Z_⊙ case and show in Fig. 5 the gas surface densities in the inner 1200 × 1200 au² box for the initial 395 kyr of disk evolution. All four models are presented, starting from the most massive model 1_Z0.1 (left column, M_core = 1.66 M_⊙) and ending with the least massive model 4_Z0.1 (right column, M_core = 0.5 M_⊙). Models 1_Z0.1 and 2_Z0.1 feature vigorous gravitational instability and fragmentation throughout the entire considered period of evolution. The disks in these models have fragmented into highly dynamical clusters that comprise one massive and several less massive clumps. What remained of the original disk is now concentrated around the central star and is characterized by a notably smaller radius. The masses of the most massive clumps are 27.6 and 18.7 M_Jup, indicating the possible formation of low-mass brown dwarfs on wide orbits.

Fig. 5. Disk evolution in models with metallicity ten times lower than solar. Each column shows a particular model as indicated in the top row. Each row corresponds to a specific time since the disk formation instance as indicated in the left column.

On the other hand, models 3_Z0.1 and 4_Z0.1 show disk fragmentation only during the initial 100–200 kyr of evolution. Afterward the disk becomes either weakly unstable or even nearly axisymmetric. This difference in the disk evolution between models 1_Z0.1 and 2_Z0.1, on the one hand, and models 3_Z0.1 and 4_Z0.1, on the other hand, can be explained by two factors. First, more massive cores produce more massive disks and this favors stronger gravitational instability and fragmentation in models 1_Z0.1 and 2_Z0.1. Second, turbulent viscosity represents an additional mechanism for mass transport, thus effectively reducing the disk mass via accretion and spreading the disk outward. Turbulent viscosity also represents an additional source of heating (the last term on the right-hand side of Eq. (3)), thus warming up the disk and reducing the strength of gravitational instability. These factors favor stronger instability in model 2_Z0.1 with α = 10⁻⁴ in comparison to the otherwise identical model 3_Z0.1 but with α = 10⁻². A similar dependence of GI on the strength of turbulent viscosity was also reported by Rice & Nayakshin (2018).

When compared to the solar metallicity case, no quantitative difference is seen in the disk evolution of Z = 0.1 Z_⊙ models. Disks with both metallicities develop gravitational instability and fragmentation, and the strength and longevity of these phenomena are similar when models with similar initial core masses are considered (e.g., model 1_Z1.0 and model 3_Z0.1). This is in agreement with what was found by Boss (2002) and Meru & Bate (2010), but is in stark disagreement with the findings of Cai et al. (2006) who reported no disk fragmentation for subsolar metallicity disks. We note, however, that the last authors did not consider mass-loading from the infalling envelope, which can qualitatively change the disk evolution and drive it to the fragmenting state (Vorobyov & Basu 2005; Kratter et al. 2008). It appears that reducing the metallicity by a factor of 10 does not notably change the overall disk evolution, as was also the case for the pre-disk collapse phase in Fig. 2.

Now we turn to the least metal-rich models of our sample. Figure 6 presents the disk evolution in four models with Z = 0.01 Z_⊙ at the same time instances as for the models with Z = 0.1 Z_⊙. We note that the initial masses and ratios of rotational to gravitational energy in the two sets of models are similar. However, the disk appearance in the Z = 0.01 Z_⊙ models is strikingly different. The disks in all models become virtually axisymmetric after 250 kyr of evolution. GI and fragmentation exist only in the very early stages of disk evolution. Model 4_Z0.01 with the lowest initial core mass (M_core = 0.52 M_⊙) does not show disk fragmentation at all, while in the other models disk fragmentation ends by 100 kyr. A notable exception is the most massive model 1_Z0.01, for which disk fragmentation extends slightly beyond 100 kyr. The low-viscosity model with α = 10⁻⁴ forms a notably more compact disk (no viscous spreading) and its disk appears to be slightly more unstable than its α = 10⁻² counterpart.

Fig. 6. Similar to Fig. 5, but for metallicity 100 times lower than solar.

The reason for the apparently different behavior between the Z = 1.0 − 0.1 Z_⊙ disks and Z = 0.01 Z_⊙ disks lies in the thermal evolution of the considered models. We showed in Fig. 1 that the models with Z = 0.01 Z_⊙ have systematically higher gas temperatures, dust temperatures, and infall velocities than the models with higher metallicity. The immediate consequence is that the mass infall rate on the disk is higher by up to a factor of 5 in the lowest metallicity models, affecting the entire disk evolution. We note that the infall rate is proportional to the gas temperature to the power 3/2 (see Eq. (10)). The comparison of Figs. 3 and 5, on the one hand, and Fig. 6, on the other hand, demonstrates that the disk appears more developed at 20 kyr after its formation in models with the lowest metallicity. With a higher infall rate, the embedded phase in models with Z = 0.01 Z_⊙ lasts much less time (for the same mass reservoir in the pre-stellar core), and this also acts to shorten the gravitationally unstable phase. For instance, the embedded phase in model 3_Z0.1 ends about 320 kyr after disk formation, while in model 3_Z0.01 it ends about 40 kyr after disk formation. When making these calculations, we defined the end of the embedded phase as the time instance when less than 10% of the initial core mass was left in the infalling envelope. We note that a shorter embedded phase in the Z = 0.01 Z_⊙ models does not prevent the disk from fragmenting because disk fragmentation takes place on orbital timescales, which are still much shorter than the time extent of the embedded phase. After the embedded phase no accretion onto the disk takes place and the gas surface density cannot easily increase globally to sustain instability.

Moreover, the temperature decoupling between gas and dust can also contribute to the increased rate of gas infall on the disk in the lowest-metallicity models. As an example, we present in Fig. 7 the spacetime diagrams of the azimuthally averaged gas and dust temperatures in models 2_Z1.0, 4_Z0.1, and 4_Z0.01. The initial 300 kyr of disk evolution is shown and the time is counted from the disk formation instance. These are the least massive models with disks having a regular (not highly fragmented) structure and are therefore best suited to demonstrate the temperature trends. We note that the horizontal high-temperature spikes are caused by accretion bursts discussed in detail in the next section.

Fig. 7. Spacetime diagrams of the gas temperature (left column) and dust temperature (right column) in models with different metallicities, as indicated in each panel. The color bar is in units of Kelvin on the log scale. The time is counted from the disk formation instance.

The comparison between models with different metallicities reveals two key differences in the gas and dust radial distribution. First, the higher metallicity disks are on average hotter in the inner several tens of astronomical units than their low-metallicity counterparts. This trend is caused by a lower optical depth in low-metallicity disks. As a result, the viscous heat released in the disk midplane can more easily be radiated away by dust thermal emission. There are, however, subtle deviations from the described trend. Most notably, the Z = 0.01 Z_⊙ model features a warmer disk in the initial stages of disk evolution (t ≤ 0.05 Myr) than the Z = 0.1 Z_⊙ one. This is caused by a higher accretion rate and associated higher accretion luminosity (stellar irradiation) in the lowest metallicity model (see Fig. 8 below). This high-accretion stage is short-lived, however, and is limited by the embedded phase.

Fig. 8. Mass accretion rate on the star (black line) and mass infall rates on the disk (red lines) vs. time in models with metallicities lower than solar. The left column presents models with metallicity Z = 0.1 Z⊙, while the right column is for the Z = 0.01 Z⊙ metallicity. The time is counted from the onset of gravitational collapse and is shown on the log scale to better resolve the variability at early times.

Second, the gas and dust temperatures decouple from each other in the vicinity of the disk outer edge and in the inner envelope. Moreover, the magnitude of temperature decoupling increases with decreasing metallicity and is remarkable in the Z = 0.01 Z_⊙ case, as can be seen from a local peak in the gas temperature in the 10²–10³ au region (such a peak is absent in the dust temperature distribution). The temperature decoupling occurs in regions with low density (see Eq. (9)) occupied by the outer disk and the inner envelope. The gas temperature is notably higher there than the dust temperature. As a result, the mass infall rate M_infall also increases, and this in turn makes the disk evolve faster for the same mass reservoir in the collapsing core.

In Sect. 4, we demonstrated that gravitationally unstable disks of solar metallicity can be characterized by highly variable mass accretion rates with episodic bursts triggered by inward-migrating clumps. Here we review the character of protostellar mass accretion in models with lower metallicities in the Z = 0.1 − 0.01 Z_⊙ range. The black lines in Fig. 8 present the protostellar mass accretion rates (Ṁ) as a function of time in models with Z = 0.1 Z_⊙ (left column) and Z = 0.01 Z_⊙ (right column). The mass infall rates on the disk (Ṁ_infall) are also plotted with the red lines. This quantity is calculated as a mass flux at a radial distance of 1000 au from the star; all our model disks are smaller than 1000 au in radius. Several interesting features can be seen in the figure.

First, the mass accretion rate is highly variable in all models, but the variability period is limited to the initial stages of disk evolution. The initial smooth behavior of Ṁ typical of the pre-disk phase is followed by a highly variable phase once the disk forms at several tens of thousands of years after the onset of gravitational collapse. The variability in the mass accretion rate diminishes together with a diminishing mass infall rate. As was noted in Sect. 4, mass-loading from the envelope helps to replenish the disk material lost via protostellar accretion and keep the disk in the gravitationally unstable mode. Once the infall is over, the disk quickly attains a stable state that is characterized by a smoothly declining mass accretion rate. Secondly, the variability is stronger and the bursts are more numerous in more massive models. This is again related to the strength of gravitational instability and fragmentation in more massive disks that form in models with initially more massive pre-stellar cores. Thirdly, the burst activity seems to be more pronounced in models with Z = 0.1 Z_⊙ than in models with Z = 0.01 Z_⊙. This is expressed in a longer time period of variability and a higher number of bursts. Finally, both Ṁ and Ṁ_infall are systematically higher in the Z = 0.01 Z_⊙ models than in the Z = 0.1 Z_⊙ ones during the embedded phase of disk evolution (40 and 320 kyr for Z = 0.01 and Z = 0.1 Z_⊙, respectively). A more rigorous statistical analysis of the burst, including their light curves and possible similarity to the classical FU Orionis-type bursts in the solar metallicity environments, will be performed in a follow-up study.

Finally, we want to demonstrate the similarity between accretion burst mechanisms in different metallicity environments caused by disk fragmentation. Figure 9 provides an example of clump infall causing an accretion burst in model 1_Z0.01. Shown is the gas surface density in the inner 200  ×  200 au² box for a limited time period of 110 yr. The inward-migrating clump is highlighted by the yellow arrows. The overall picture is similar to what is shown in Fig. 3 in the context of solar metallicity bursts. The clump spirals down toward the star until it is destroyed by tidal torques in the immediate vicinity of the sink cell. The mass of the clump at this stage is about 10 M_Jup; its radius is about 2 au, while the Hill radius is 1.5 au, meaning that the clump has started to lose its material, which is quickly accreted on the star producing a burst.

Fig. 9. Inward migration and infall of a clump in the gravitationally unstable disk with Z = 0.01 Z⊙. Shown is the gas surface density in log g cm−2. The inward-migrating clump is highlighted by the yellow arrows.

A potential caveat needs to be mentioned with respect to this burst scenario. The numerical resolution of the clump is insufficient to resolve the second collapse (see Sect. 2). This means that part of the clump may avoid tidal destruction if collapsed to protoplanetary densities and sizes (see, e.g., Vorobyov & Elbakyan 2018). In this case the burst might be of a notably smaller magnitude than found here. These topics deserve detailed follow-up studies.

6. Conclusions

We studied numerically the early stages of disk evolution and protostellar accretion in young stellar systems with metallicities from 10 to 100 times lower than solar. For this purpose we employed the numerical hydrodynamics simulations in the thin-disk limit that feature separate gas and dust temperatures and consider disk mass-loading from collapsing parent cloud cores. Low metallicities were set by scaling down the gas and dust opacities, dust-to-gas mass ratio, and metal content of the solar metallicity disk by the corresponding factor. To eliminate the possible uncertainty that may be caused by the initial conditions, we considered the gravitational collapse of pre-stellar cores with similar mass and ratio of rotational-to-gravitational energy, but distinct metallicity. A comparison with the solar metallicity models is also performed. Our findings can be summarized as follows:

– The initial stages of cloud core collapse are distinct in models within the considered metallicity range. In particular, the Z = 0.01 Z_⊙ models are characterized by higher temperatures and infall velocities than their higher metallicity counterparts (Z = 1.0 − 0.1 Z_⊙). The gas and dust temperatures decouple notably in the model with the lowest metallicity.

– The initial stages of disk evolution are characterized by vigorous gravitational instability and fragmentation. The time period that is covered by this unstable stage, however, is much shorter in the Z = 0.01 Z_⊙ models than in their higher metallicity counterparts. This can be explained by elevated mass infall rates on the disk in models with the lowest metallicity, causing the corresponding disks to evolve faster for the same initial mass reservoir in the pre-stellar cores. Once the disk mass-loading from the infalling envelope diminishes (which occurs sooner in lower metallicity models), the disk quickly acquires a gravitationally stable state.

– Protostellar accretion rates are highly variable in the low-metallicity models reflecting the highly dynamical state of the corresponding protostellar disks. Similar to the solar metallicity models, the Z = 0.1 − 0.01 Z_⊙ systems feature short but energetic episodes of mass accretion caused by infall of inward-migrating gaseous clumps that form in the outer disk via gravitational fragmentation. The bursts seem to be more numerous and longer-lasting in the Z = 0.1 Z_⊙ models than in the Z = 0.01 Z_⊙ case.

Acknowledgments

We are thankful to the anonymous referee for useful comments that helped to improve the manuscript. E. I. V. and M. G. acknowledge support from the Austrian Science Fund (FWF) under research grant P31635-N27. V. G. E. acknowledges the Swedish Institute for a visitor Grant allowing to conduct research at Lund University. The simulations were performed on the Vienna Scientific Cluster.

References

All Tables

All Figures

Fig. 1. Radial profiles of the gas surface density (top left), radial infall velocity (top right), gas temperature (bottom left), and dust temperature (bottom right) in model Z1.0. The dashed curves correspond to the initial setup at the onset of gravitational collapse, while the solid curves show the quantities at the time instance immediately preceding disk formation. Curves of different color correspond to models with different metallicities, as shown in the legend.

In the text

	Fig. 2. Disk evolution in model 1_Z1.0. Shown is the gas surface density at six evolutionary times since the disk formation instance. The time is counted starting from the instance of disk formation, which occurs at 89 kyr in absolute time elapsed since the onset of gravitational collapse. The scale bar is in g cm−2.
In the text

	Fig. 3. Disk evolution in model Z1.0. Shown is the gas surface density at six evolutionary times encompassing a short time period of clump migration at t ≈ 240 kyr. The scale bar is in g cm−2.
In the text

	Fig. 4. Mass accretion rate vs. time in model 1_Z1.0 shown with the black line. The red line presents the mass infall rate on the disk. The time is counted from the onset of gravitational collapse and is shown on the log scale to better resolve the variability at early times. The disk forms at 89 kyr.
In the text

	Fig. 5. Disk evolution in models with metallicity ten times lower than solar. Each column shows a particular model as indicated in the top row. Each row corresponds to a specific time since the disk formation instance as indicated in the left column.
In the text

	Fig. 6. Similar to Fig. 5, but for metallicity 100 times lower than solar.
In the text

	Fig. 7. Spacetime diagrams of the gas temperature (left column) and dust temperature (right column) in models with different metallicities, as indicated in each panel. The color bar is in units of Kelvin on the log scale. The time is counted from the disk formation instance.
In the text

Fig. 8. Mass accretion rate on the star (black line) and mass infall rates on the disk (red lines) vs. time in models with metallicities lower than solar. The left column presents models with metallicity Z = 0.1 Z⊙, while the right column is for the Z = 0.01 Z⊙ metallicity. The time is counted from the onset of gravitational collapse and is shown on the log scale to better resolve the variability at early times.

In the text

	Fig. 9. Inward migration and infall of a clump in the gravitationally unstable disk with Z = 0.01 Z⊙. Shown is the gas surface density in log g cm−2. The inward-migrating clump is highlighted by the yellow arrows.
In the text