Open Access
Issue
A&A
Volume 669, January 2023
Article Number A84
Number of page(s) 17
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202244408
Published online 16 January 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

In classical pre-main-sequence (pre-MS) evolution, a star forms with a large radius and contracts isotropically along the Hayashi track (Hayashi 1961). During the contraction, the star heats up due to the release of gravitational energy, nuclear reactions are triggered, a radiative core develops, and the star evolves toward the main sequence (MS). Within this standard model, however, the way the star gets its mass is not taken into account.

Within a collapsing cloud, two hydrostatic cores form (e.g., Larson 1969; Bhandare et al. 2018). Starting from the second hydrostatic core (in the context of this work, the stellar seed), material from the parent cloud and the forming accretion disk accretes onto the young star (e.g., Mercer-Smith et al. 1984; Palla & Stahler 1992; Hartmann et al. 1997; Baraffe et al. 2009; Dunham & Vorobyov 2012; Kunitomo et al. 2017; Vorobyov et al. 2017a; Steindl et al. 2021). The nature of the accretion process depends on different parameters (e.g., parent cloud mass, temperature, metallicity, and rotational rate), resulting in a variety of possible accretion rates over time (e.g., Vorobyov & Basu 2010; Vorobyov et al. 2017a, 2020).

The pre-MS evolution starting from the stellar seed is strongly influenced by the amount of energy added to the stellar interior (Ladd) during the accretion process (e.g., Kunitomo et al. 2017; Steindl et al. 2021). If the accreting material loses its energy during the accretion process (Ladd → 0; cold accretion), the stellar radius is significantly smaller (by up to one order of magnitude) compared to the classical pre-MS evolution (e.g., Kunitomo et al. 2017). If, on the other hand, the accreting material can add its energy to the stellar interior (Ladd is large; hot accretion), the stellar radius inflates and can become comparable with classical evolutionary models. After accretion has ended and energy is no longer added to the stellar interior, the star can evolve along the Hayashi track (e.g., Kunitomo et al. 2017). We note that the exact pre-MS evolutionary track, starting from a stellar seed, depends on how much energy can be added to the star and where this energy is added in the stellar interior (e.g., Kunitomo et al. 2017; Steindl et al. 2021). In addition to the cold and warm accretion models, Ladd can also depend on the stellar accretion rate (hybrid accretion; e.g., Vorobyov et al. 2017a). In periods of high accretion rates Ladd is large (similar to hot accretion), and during low accretion rates Ladd stays small (similar to cold accretion).

During the early phases of stellar evolution (Class 0 and Class I), the rate of accretion onto the star can vary significantly over time. During episodic outburst events (e.g., FU Orionis outbursts), the stellar accretion rate can rise by orders of magnitude on a short timescale (on the order of years; e.g., Vorobyov & Basu 2010; Vorobyov et al. 2020). With Ladd coupled to the accretion rate, episodic outburst events during the early disk phase (e.g., FU Orionis outbursts) can lead to radial oscillations (e.g., Bastien et al. 2011; Vorobyov et al. 2017a) that affect the stellar pre-MS evolution. In previous studies, the accretion history on the star is either assumed to be constant over time or precomputed (e.g., Kunitomo et al. 2017; Vorobyov et al. 2017a; Jensen & Haugbølle 2018; Elbakyan et al. 2019; Steindl et al. 2021). In the latter case, the effect of episodic outbursts can be included. For now, however, such calculations only deal with the stellar structure and are not able to include any back-reaction of the evolving star on the accretion disk.

Star–disk interactions furthermore influence the stellar rotation period. Understanding the stellar spin evolution allows the determination of stellar ages (gyrochronology; Barnes 2007), magnetic activity, or high-energy stellar radiation (e.g., Pallavicini et al. 1981; Micela et al. 1985; Pizzolato et al. 2003; Wright et al. 2011; France et al. 2018). During the disk phase, the stellar rotation period can also influence the evolution of the accretion rate. Fast-rotating stars tend to experience more outbursts compared to slow-rotating stars (e.g., Gehrig et al. 2022). According to recent spin evolution models (e.g., Matt et al. 2010, 2012; Gallet et al. 2019; Ireland et al. 2021), the stellar accretion rate is an important factor for determining the stellar rotation rate. A precise calculation of the accretion rate is crucial for understanding the rotational evolution of young stars. Furthermore, stellar metallicity has an effect on the stellar spin evolution (e.g., Amard et al. 2019; Amard & Matt 2020). Amard et al. (2019) present a grid of low-mass (≤1.5 M) isochrones and a spin evolution model. The disk phase, however, is treated in a simplified way. During the disk phase, the stellar period is assumed to be constant (disk-locking) and the disk lifetime is a free parameter that is independent of stellar mass and metallicity. This assumption is in conflict with observations of disk fractions in low-metallicity clusters (e.g., Yasui et al. 2010, 2016, 2021; Guarcello et al. 2021) and theoretical, metallicity-dependent photo-evaporation models (e.g., Nakatani et al. 2018), which indicate short disk lifetimes in low-metallicity environments.

In this study we aim to combine The Adaptive Implicit RHD (TAPIR) code (Ragossnig et al. 2020; Steiner et al. 2021) with the software instrument Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019) to provide a self-consistent numerical description of a star–disk system at a young age. This allows us to study which parameters affect, for example, the disk lifetime, the stellar spin, and the pre-MS evolution, including the back-reaction on the other component. In this first paper, which introduces the new method, we study the effects of metallicity on the disk and stellar spin evolution of T Tauri stars (with ages ≳1 Myr). We seek to understand how stellar metallicity affects the accretion disk and stellar spin evolution of T Tauri stars.

This paper is structured as follows: Sect. 2 describes the disk and stellar evolution models used in this study. Moreover, the combination of the disk, spin, and stellar model is introduced. Our results are presented in Sect. 3 and discussed in Sect. 4. Finally, we draw our conclusions in Sect. 5.

2. Model description

2.1. Hydrodynamic disk evolution with the TAPIR code

Studying the long-term evolution of an accretion disk including the inner disk regions close to the star, we use the implicit TAPIR code, which is described in detail in Ragossnig et al. (2020) and Steiner et al. (2021). The key advantage of our numerical method is the treatment of hydrodynamic equations, rather than a diffusive disk evolution approach (e.g., Pringle 1981; Armitage et al. 2001; Zhu et al. 2007; Schib et al. 2021). We can include the influence of pressure gradients and a stellar magnetic field, which is of particular importance in the inner disk regions (e.g., Romanova et al. 2002; Bessolaz et al. 2008). The following summarizes some key features of the TAPIR code.

Basic equations. In our model we use a time-dependent, vertically integrated, viscous accretion disk formulation (e.g., Shakura & Sunyaev 1973; Armitage et al. 2001; Steiner et al. 2021):

(1)

(2)

(3)

where Σ, u, e, Pgas, and HP stand for the gas column density, gas velocity in the planar components u = (ur, uφ), specific internal energy density per surface area, vertically integrated gas pressure, and the vertical scale height of the gas disk, respectively. The gradient in planar cylindrical coordinates reads ∇ = (∂/∂r, r−1∂/∂φ) with ∂/∂φ = 0 for our axisymmetric model. The Pgas is utilized by the ideal equation of state Pgas = Σe(γ − 1), with the adiabatic coefficient γ = 5/3, and κR denotes the Rosseland-mean opacity. It is composed of a gaseous component, κR, gas, which is based on Ferguson et al. (2005) with the abundances in Caffau et al. (2011; see Fig. 1 in Steiner et al. 2021), and a dust-dominated component κR, dust (based on Pollack et al. 1985); with κR = κR, gas + fdustκR, dust, where fdust is the dust-to-gas mass ratio set to the respective stellar metallicity. In our model, fdust equals the stellar metallicity, ranging from Z = 0.0028 to Z = 0.014. The gravitational potential of the star-disk system is denoted by ψ. The vertical component of the stellar magnetic field Bz is assumed to be a dipole and constant within the disk’s vertical extent. The planar magnetic field components B = (Br, Bφ) are taken at the surface of the disk. Similar to Rappaport et al. (2004) and Kluźniak & Rappaport (2007), we ignore the radial magnetic field component Br = 0. Differential rotation between the stellar magnetic field lines, which are assumed to be in rigid rotation with the star, and the disk’s material generates an angular magnetic field component Bφ (e.g., Rappaport et al. 2004; Kluźniak & Rappaport 2007; Steiner et al. 2021; Gehrig et al. 2022). The radial radiative transport (depicted by (J − S)) is modeled in a radiative diffusion approximation with an Eddington factor of 1/3. We note that an Eddington factor of 1/3 is usually valid within optically thick regions or in optically thin regions in the presence of a uniform radiation field. We assume that the stellar radiation and the ambient temperature roughly meet this requirement. The method, however, has its weaknesses when transitioning between optically thick and thin regions and a variable method would produce better results. Here, J and S stand for the zeroth moment of the radiation field and the source function, respectively (e.g., Ragossnig et al. 2020; Steiner et al. 2021). Finally, Ėrad depict the net heating and cooling rate per unit surface area and Q depict the viscous stress tensor (e.g., Steiner et al. 2021; Gehrig et al. 2022). The kinematic viscosity is formulated according to Shakura & Sunyaev (1973)ν = αcSHP; with the viscous α parameters, the isothermal sound speed cS and the pressure scale height HP. In addition to Eqs. (1)–(3), we utilize an adaptive grid (Dorfi & Drury 1987) to allow a moving inner and outer disk boundary as well as adequate radial grid-point resolution throughout the computational domain.

Computational domain. The innermost disk regions are assumed to be truncated and eventually disrupted by a strong stellar magnetic field (∼1 kG) inside the so-called truncation radius, rtrunc (e.g., Romanova et al. 2002; Bessolaz et al. 2008; Romanova & Kurosawa 2014). With our one-dimensional code, the physical processes inside rtrunc cannot be described. Thus, we chose rtrunc to be the inner boundary of our computational domain. The radial position of rtrunc can be calculated by balancing the stellar magnetic pressure with the maximum from the ram pressure of the accreting disk material and the gas pressure (e.g., Koldoba et al. 2002; Romanova et al. 2002; Bessolaz et al. 2008; Steiner et al. 2021; Gehrig et al. 2022):

(4)

If Pram ≥ Pgas, the truncation radius can be calculated following Hartmann et al. (2016):

(5)

with the correction factor ξ = 0.7, the stellar dipole magnetic field strength B, the stellar radius R, the stellar mass M and the accretion rate on the star . In the case of Pram < Pgas, the position of the truncation radius is calculated by equating the magnetic pressure of the vertical field component and the gas pressure:

(6)

where rtrunc is allowed to move during our simulations according to stellar and disk evolution. Furthermore, the outer boundary of our domain is set to the radial position where the surface density has reached a certain value Σout = 1.0 g cm2 (similar to, e.g., Vorobyov et al. 2020).

Boundary conditions. At the inner boundary, we chose “free” (Neumann) boundary conditions (e.g., ∂Σ/∂r = 0, Romanova et al. 2002, 2004; Steiner et al. 2021). At the outer boundary, we fix the radial (angular velocity) ur (uφ) to zero (the Keplerian value). This indicates no inflow of additional material from an outer envelope or parent cloud, as expected in late Class II star-disk systems. The surface density Σ and the specific internal energy e are again treated with a Neumann boundary condition.

2.2. Stellar evolution with MESA

The stellar evolution calculations in this work are performed with version v12778 of MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019). MESA is an open-source software instrument that solves the fully coupled structure and composition equations simultaneously for a spherically symmetric and hence one-dimensional stellar model (Paxton et al. 2011). Details about the adopted micro-physics are given in Appendix A. The treatment of rotation in MESA is summarized in Appendix B.

The input physics for the stellar evolution calculations performed in this work is identical to Steindl et al. (2021). The initial models used in this study are identical to the ones used in Steindl et al. (2021) with a composition of X1H = 1 − X2H − X3He − X4He − Z, where X2H = 20 ppm, X3He = 85 ppm, X4He = 0.276, and Z = 0.014 for the models with solar metallicity or Z = 0.0028 for models with a fifth of solar metallicity. In the following, we briefly discuss implemented changes from the usual MESA calculations for the treatment of accretion.

We followed the description of Baraffe et al. (2009) for which the accretion onto the protostar proceeds non-spherically, allowing the central object to freely radiate energy across most of the photosphere (Hartmann et al. 1997). The energy budget of the accretion process is given by

(7)

which is a combination of the gravitational (−GM/R) and internal energy (+ϵGM/R) of the accreted material. The geometry of the accretion process is described by the free parameter ϵ. We follow the approach of previous works (Baraffe et al. 2009; Jensen & Haugbølle 2018; Steindl et al. 2021) and set ϵ = 0.5, corresponding to accretion from a thin disk around the equator (Hartmann et al. 1997; Baraffe et al. 2009). Thus, the total energy rate accreted to the star from a thin disk is given as ϵGM/R.

The energy of the accreted material is either added to the star as extra heat (Ladd) or radiated away as accretion luminosity (Lacc). The accretion luminosity is added to the intrinsic stellar luminosity and the additional heating of the disk is considered in Eq. (3). A free parameter controls the corresponding amounts such that

(8)

and

(9)

The value β and its dependence on star and disk properties are free parameters. From earlier studies it is evident that β varies with accretion rate (see, e.g., Baraffe et al. 2012; Jensen & Haugbølle 2018; Elbakyan et al. 2019). In this work, we use a functional dependence first presented by Jensen & Haugbølle (2018) given as

(10)

This describes a step function with a smooth transition between β1 and β2. The functional dependence is controlled by four parameters: β1 = 0.005, β2 = 0.2, a midpoint accretion rate of m = 6.2 × 10−6 M yr−1 and a crossover width of Δ = 5.95 × 10−6M yr−1. For comparability, we use the same values as in Steindl et al. (2022), which were originally adapted from Jensen & Haugbølle (2018) for numerical stability. In our case, this translates to a constant value, β, at the beginning of the evolution that becomes gradually smaller as the accretion rate decreases (see Sect. 3.1). During the calculation of the initial models, the accretion rates and thus β() vary significantly. Given the proportionately small accretion rates, β is basically constant during the joint evolution of MESA and TAPIR.

The dependence of β on the accretion rate can also be physically interpreted. During phases of low accretion rates, the stellar magnetic field can truncate the disk at several stellar radii, and the disk material is funneled on the star (e.g., Bessolaz et al. 2008; Hartmann et al. 2016). The in-falling disk material generates a shock close to the stellar surface and most of the energy is radiated away (e.g., Koenigl 1991). As a consequence, only a small part of the disk material’s energy is assumed to be added to the stellar interior. On the other hand, during phases of high accretion rates, the disk is pushed toward the star, and accretion along magnetic field lines ceases. Without the shocks at the end of the magnetic funnels, a larger amount of energy can be added to the star. The transition between the regimes of low and high β values is expected to occur at ∼10−6 − 10−5M yr−1 (e.g., Vorobyov et al. 2017a; Jensen & Haugbølle 2018), which motivates the choice of m. We note, however, that the description of β in Eq. (10) does not take metallicity into account. We assume that, in case of low accretion rates, the inner part of the disk is sufficiently ionized due to stellar irradiation. Stellar magnetic field lines can couple to the disk material and the accretion geometry is unaffected by metallicity. During phases of high accretion rates, when the disk is pushed toward the stellar surface, differences in disk ionization due to metallicity are assumed to be negligible as well. The transition between the two accretion regimes, however, could be affected by the disk’s ionization fraction, which depends on metallicity. Thus, the value Δ should include a metallicity-dependent factor; however, that is beyond the scope of this study.

The heat added to the stellar model is distributed according to

(11)

which follows the approach of Kunitomo et al. (2017). In this picture, the heat is distributed only in an outer region of fractional mass, Mouter, and increases linearly with the mass coordinate, mr. In our study, we chose Mouter = 0.01, hence injecting the heat only in the outer 1% of the stellar structure.

2.3. Stellar spin model

The stellar spin model used in our simulations is based on Matt et al. (2010) and Gallet et al. (2019) and summarized in detail in Gehrig et al. (2022). In this study we assume solid body rotation for the star. Angular momentum transport within the star is discussed in Sect. 4.3. The stellar angular momentum J = IΩ is influenced by external torques Γext and its temporal derivative reads

(12)

with the stellar angular momentum Ω and moment of inertia I. Rearranging Eq. (12) results in a time-dependent equation for the stellar spin:

(13)

The stellar moment of inertia, I, is calculated by the MESA code. Within Γext the effects of accretion, Γacc, stellar winds, ΓW, in the form of an accretion-powered stellar wind (APSW; Matt & Pudritz 2005), and the magnetic star-disk connection, ΓME, in the form of magnetospheric ejections (MEs; Zanni & Ferreira 2013) are combined. During the process of accretion, the star gains mass and angular momentum from the disk, which increases I and spins up the star. An APSW, on the other hand, is assumed to eject a specific amount of the accreted material (W = WṀ) removing stellar mass and angular momentum. W is expected to be ≲2% (e.g., Cranmer 2008; Pantolmos et al. 2020) and thus, we chose W = 2% in this study. The effect of MEs on the stellar spin evolution depends on the position of rtrunc with respect to the co-rotation radius:

(14)

with the gravitational constant G and the stellar mass M. The MEs increase the stellar angular momentum and spin up the star if . Otherwise, the star spins down. Following Gallet et al. (2019), we used Krot = 0.7 in this study1, and the external torque contributions scale as

(15)

(16)

(17)

with Kacc = 0.4, KME = 0.21, the disk rotation rate at the truncation radius Ωdisk, in, and the stellar Alfvén radius

(18)

where , m = 0.2177, K2 = 0.0506, and K1 = 1.7.

2.4. Numerical combination of the disk, spin, and stellar model

The interaction between the disk (TAPIR, Sect. 2.1), spin (stellar spin model, Sect. 2.3), and stellar model (MESA, Sect. 2.2) is schematically shown in Fig. 1. Starting from the initial or old values at timestep t0, the new disk values, spin values, and stellar values at timestep t1 = t0 + δt are calculated in three steps, which are represented by numbered arrows in Fig. 1. We note that we calculate the stellar spin evolution and the MESA evolution separately. During the MESA calculation, the star is assumed to rotate at a constant rate. In the first step, the TAPIR code calculates the disk evolution from timestep t0 to t1 using the old disk, spin, and stellar values (Arrow 1). Then, the spin values are updated based on the new disk values (Arrow 2). Finally, the MESA code is updated and evolves the stellar parameters from timestep t0 to t1 (Arrow 3).

thumbnail Fig. 1.

Schematic representation of the interaction between TAPIR and MESA. Starting from an initial (old) timestep, t0, the disk, spin, and star values are evolved toward the new timestep, t1 = t0 + δt. The new values at t1 are obtained in three steps, represented by the numbered arrows (see text).

The accuracy and reliability of our method depend on the timestep size δt. Choosing small timesteps results in accurate but long computational times. Large timesteps will reduce computational time but accuracy suffers. We chose δt = 100 yr in this study. During 1 Myr the stellar values are updated 104 times resulting in an acceptable representation of the stellar evolution. We also tested a range of δt from 1 yr to 1000 yr and found no significant change in our results. We note that MESA and TAPIR use multiple (internal) timesteps to evolve the star and the disk from time t0 to t1. For MESA, the internal timesteps range between ∼106 s and ∼108 s. For TAPIR, the timesteps range between ∼104 s and ∼108 s (see Fig. A.1 in Steiner et al. 2021). We note that we do not model for example episodic outburst events in this study. During these outbursts, the accretion rate onto the star can change significantly within short timescales and influence the stellar evolution (e.g., Vorobyov et al. 2017a). In further studies, such outbursts are included in our model and δt has to be adapted to lower values to assure a sufficiently accurate stellar evolution calculation.

3. Results

To study the influence of metallicity on the combined evolution of a T Tauri star and its disk, several stellar and disk parameters have to be discussed. These parameters are defined based on observations, and previous theoretical or numerical studies and, within a reasonable range, initial models for a T Tauri star and its disk are constructed. We chose t0 = 2 Myr as an initial time for our combined simulations. The effects of the early evolution and different values for t0 are discussed in Sect. 4.2. For each model, the stellar metallicity is varied between Z = 0.2 Z (motivated by observations of young, low-metallicity clusters; Yasui et al. 2010, 2016, 2021) and Z = 1 Z.

Disk parameters. One important parameter that defines an accretion disk is the mass accretion rate. For T Tauri star-disk systems (≲10 Myr) accretion rates range within 10−10 and 10−7M yr−1 (e.g., Gullbring et al. 1998; Manara et al. 2016, 2022; Vorobyov et al. 2017a). For t0 = 2 Myr, Manara et al. (2012) and Testi et al. (2022) find disk accretion rates for a ∼1 M star between 10−8M yr−1 and 10−9M yr−1. Furthermore, the viscous α parameter can have a wide range between 0.1 and ∼0.0001 (e.g., Zhu et al. 2007, 2009; Vorobyov & Basu 2009; Mulders & Dominik 2012; Pinte et al. 2016; Yang et al. 2018; Flaherty et al. 2020). In our simulations, we chose the initial (starting) accretion rate start = [1 × 10−9, 5 × 10−9, 1 × 10−8]M yr−1 and set the (constant) viscous α-value to 0.01 unless otherwise stated. This combination results in the expected disk lifetimes of ≲10 Myr (consistent with, e.g., Richert et al. 2018). We stopped our simulations when the disk accretion rate decreased to 10−11M yr−1. This lower limit is motivated by the lowest detectable accretion rate for a 1 M star (e.g., Sicilia-Aguilar et al. 2016). At such low accretion rates, the disk is expected to be dissolved rather quickly due to photo-evaporation, which is currently not included in our model. We note that an adaptive viscosity description will be added in further studies.

Stellar parameters. The inner disk region close to the star is influenced by a stellar magnetic field. The field strength of young T Tauri stars can vary between several hundred to several thousand Gauss (e.g., Johnstone et al. 2014; Lavail et al. 2017, 2019). Unless otherwise stated, the magnetic field strength is chosen to be 2.0 kG. In addition, stellar parameters such as radius, luminosity, and effective temperature at t0 result from the initial stellar evolution (t < t0). The initial stellar parameters are chosen as described in Kunitomo et al. (2017), which follow Stahler (1988) and Hosokawa et al. (2011). We start from a fully convective stellar seed with a mass of 0.01 M (≈10 Jupiter masses) and a radius of 1.5 R (e.g., Kunitomo et al. 2017; Steindl et al. 2021). As pointed out in Hosokawa et al. (2011), the initial entropy of the stellar seed influence the subsequent stellar evolution in the case of cold accretion. Cold accretion, however, is assumed to be unrealistic for most stars (e.g., Kunitomo et al. 2017; Vorobyov et al. 2017a). In the case of hot or warm accretion, the differences due to a variation in the initial stellar seed parameters are expected to be small (e.g., Stahler 1988; Kunitomo et al. 2017). Thus, we do not include a variation of the initial stellar seed parameter in this study and refer to Appendix D in Kunitomo et al. (2017) for further information. Each model accretes material until it has reached a stellar mass of 0.91 M and the desired accretion rate start at time t0 = 2 Myr. The effect of the accretion history is described in Sect. 3.1. During the T Tauri phase, the rotational period can span over an order of magnitude from ≲1 toward ∼10 days (e.g., Serna et al. 2021). To compare the effect of different metallicity values, we start all our models at the same initial stellar rotation period Pstart = 10 days (comparable to the slow models in Amard et al. 2019). During the initial stellar evolution (t < t0), the rotation period is fixed to 10 days. We note that the rotation period does change during the first 2 Myr. The respective effects, however, are beyond the scope of this work and will be treated in subsequent studies.

In order to successfully start a simulation that combines stellar and disk evolution, initial models that represent a T Tauri star and its disk at an age of t0 have to be calculated. The initial models are based on the aforementioned stellar and disk parameters.

3.1. Initial stellar evolution: t < t0

The early evolution of a protostar depends on the temporal evolution of the accretion rate (e.g., Kunitomo et al. 2017; Steindl et al. 2021). Unfortunately, the accretion history of a young star is difficult to measure, hardly constrained, and, thus, treated as a free parameter in our model. Based on the numerical results by Vorobyov et al. (2017a), we start with a constant accretion rate init followed by a decrease toward our desired value of start according to t−2 (see Fig. 2). We note that for the models with start = 1 × 10−9 yr−1, we use t−2.3 to avoid values for init > 2 × 10−5 M yr−1. The reason for this variation is convergence problems in the MESA initial model with higher initial accretion rates (e.g., Sect. 5.1.5 in Steindl et al. 2021). The values of init and the time, at which the accretion rate begins to decrease, tdec, depend on start as well as the desired stellar mass at t0. For a starting accretion rate of start = [1 × 10−9, 5 × 10−9, 1 × 10−8] M yr−1, we find init = [1.87 × 10−5, 1.03 × 10−5, 5.29 × 10−6] M yr−1, and tdec = [28, 45, 87] kyr, respectively. We note that the accretion history of the star is identical for both metallicities. We present our results for the initial stellar evolution in three parts: from the stellar seed to the transition toward the Hayashi track corresponding to an evolution time τH, from τH along the Hayashi track toward t0 and finally reaching t0, which marks the starting point for our combined stellar and disk model.

thumbnail Fig. 2.

Initial stellar evolution for t < t0. Panel a: accretion history of the star. Initially, the star accretes with a constant accretion rate, init. At time tdec, the accretion rate drops according to t−2 (see text). The colored dashed lines show the respective value of start. Panel b: evolution of the stellar mass. Starting from the initial mass of 0.01 M, all models reach their final mass of 0.91 M after 2 Myr. The yellow, red, and blue lines show models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. To show the entire initial stellar evolution, starting from t = 0 yr, we use the symlog representation (Hunter 2007). The vertical dashed line indicates the transition from the linear scale (for t ≤ 104 yr) to the logarithmic scale (for t > 104 yr).

From the stellar seed to τH. Considering the evolution of the star on a Hertzsprung-Russell diagram (HRD) for t < t0, it is noticeable that stellar metallicity and the accretion history have a distinct influence on the star (see Fig. 3). Stars with lower metallicity are usually more compact and hotter compared to higher metallicities and their position on the HRD is moved to the upper left. The accretion history influences the amount of accretion energy added to the stellar interior (see Sect. 2.2 and Eq. (10)). In our model, the fraction of added energy depends on the magnitude of the accretion rate. A higher accretion rate init results in a higher amount of energy added to the star2. We note that a higher init results in a lower start (see Fig. 2). The additional energy causes the star to expand and the stellar luminosity increases with an increasing amount of energy added to the star. On the HRD, the stars are shifted in a vertical direction dependent on their accretion history (see Fig. 3). In other words, for a given effective temperature and luminosity, a star is younger if less energy is added to the star during its formation (e.g., Kunitomo et al. 2017). At a few thousand years, deuterium burning sets in (τD, marked with diamonds in Fig. 3) and additionally increases the stellar luminosity. The onset times range between 1.31 kyr and 2.72 kyr (see Table 1). The point of maximum luminosity (τH, marked with crosses in Fig. 3) indicates a transition. Before τH, the pre-MS star derives most of its energy from accretion and deuterium burning and afterward from contraction on the Kelvin-Helmholtz timescale (e.g., Stahler & Palla 2004). The young star turns toward the Hayashi track.

thumbnail Fig. 3.

Evolution on the HRD for t < t0 for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The values at time t = 0 Myr (t = 2 Myr) are marked with circles (triangles) and the yellow, red, and blue lines show models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. The onset of deuterium burning (t = τD) is marked with a diamond. The transition from expansion toward the Hayashi track (t = τH) is marked with a cross.

Table 1.

Onset of deuterium burning (τD), the stellar mass (MH), and evolution time (τH) when the star turns onto the Hayashi track.

Evolution on the Hayashi track. As the accretion rate decreases over time (see Fig. 2), the amount of energy added to the stellar interior decreases as well and the expansion of the radius stops. The star contracts subsequently on the Hayashi track, similar to the classical stellar evolution model. In Fig. 3, this transition toward the Hayashi track is marked with a cross (“x”). The respective stellar mass (MH) and evolution time (τH) is summarized in Table 1. The stellar masses (MH) range between 0.66 and 0.83 M (with slightly lower values for Z = 0.2 Z). The time τH depends, as expected, on the accretion history. Shortly after the accretion rate decreases from init toward start the amount of energy added to the star is insufficient to expand the radius and counteract contraction. Our models are comparable with the results of Kunitomo et al. (2017) for warm (hot) accretion models.

Reaching t0 = 2 Myr. When reaching an age of t0 = 2 Myr, stellar radii and temperatures (see Fig. 4) are used to calculate the initial disk model (see Sect. 3.2). While the stellar radius differs due to the accretion history, the effective temperature stays approximately constant. When connecting the colored triangles in Fig. 3 for Z = 1 Z, they are located on a vertical line indicating evolution on the Hayashi track (stellar luminosity decreasing with age). For Z = 0.2 Z however, the stars at t0 already turn toward the Henyey track (stellar luminosity and effective temperature increase with age) indicating the development of a radiative core.

thumbnail Fig. 4.

Stellar radii and temperatures at t0 = 2 Myr for start = 5 × 10−9 (red triangles) and 1 × 10−8 (blue triangles). Solar (subsolar) metallicity values are connected by a solid (dashed) line.

3.2. Initial disk model

Starting a simulation with an implicit method requires an initial model that already solves Eqs. (1)–(3); (e.g., Dorfi et al. 2006). Following Steiner et al. (2021), we constructed a steady-state accretion disk that serves as the initial model in our simulation. The mass flow rate throughout the disk corresponds to the respective value of start. In Fig. 5 the radial surface density (panel a), disk midplane temperature (panel b), and the optical depth, τ (panel c) profile of the initial disk models are shown.

thumbnail Fig. 5.

Radial profile of the initial disk model for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. Panel a shows the disk surface density, Σ, panel b the disk midplane temperature, Tdisk, and panel c the optical depth, τ. The thin vertical lines in panel a indicate the respective position of Rc. The horizontal line in panel b indicates the ambient temperature Tamb = 20 K. The horizontal dashed line in panel c symbolizes τ = 1, which divides optically thin regions (τ ≪ 1) from optically thick regions (τ ≫ 1).

For larger initial accretion rates, the disks are more massive, have a higher temperature, and are pushed closer toward the star (e.g., Romanova et al. 2002; Bessolaz et al. 2008; Romanova & Kurosawa 2014; Steiner et al. 2021). Due to the smaller stellar radii at t0 (see Fig. 4), the absolute values of the inner disk edge rtrunc is smaller for Z = 0.2 Z. The radial temperature profile can be divided into three regions. Close to the star the disk is directly heated by stellar irradiation r ≲ 10−1. For low metallicities, stellar luminosity is larger and, thus, the disk temperature is higher. From r ≲ 10−1 AU to ∼1 AU, the midplane of the disk is shielded from stellar irradiation and the disk is optically thick (panel c in Fig. 5). This region is referred to as the optically thick zone, Rthick, from now on. In this region, the disk midplane temperature ranges from ∼100 to 1000 K and is dominated by dust opacities (e.g., Semenov et al. 2003). In our model, we assume that the disk gas-to-dust ratio corresponds to the stellar metallicity, fdust = 0.014 and 0.0028 for Z = 1 Z and 0.2 Z, respectively. Consequently, larger stellar metallicity results in higher optical depths and temperatures in the disk midplane compared to lower stellar metallicities. We note that in the models with start = 1 × 10−9 M yr−1, these effects are almost negligible due to the low surface density values. Outside Rthick, the disk temperature is again dominated by stellar irradiation and smaller stellar metallicities result in higher disk temperatures. The disk temperature decreases with increasing disk radius until the ambient temperature of Tamb = 20 K is reached (see panel b in Fig. 5).

To compare the initial disk models to observations and other theoretical models, we summarized the initial disk sizes, given by Rc, and masses in Table 2. The characteristic radius of the disk Rc includes 2/3 of the disk’s mass. The values of Rc, ranging between 20 AU and 128 AU, and mdisk, ranging from 0.001 M to 0.055 M are within the limits of observed disk radii and masses in Ophiuchus II (Andrews et al. 2010) and Lupus star-forming regions (Ansdell et al. 2018).

Table 2.

Characteristic radius, Rc, and disk mass, mdisk, of the initial disk model.

3.3. Disk evolution: t > t0

Starting from the stellar and disk initial models, we evolve the star-disk system until the accretion rate onto the star decreases below 10−11M yr−1. In Fig. 6, the evolution of the disk accretion rate is shown for both metallicities Z = 1 Z and Z = 0.2 Z. Independent of the accretion history of the star, a lower stellar metallicity results in a shorter disk lifetime3.

thumbnail Fig. 6.

Accretion rate onto the star over time for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. The dotted black line shows the termination criterion used in our simulations of = 10−11 M yr−1.

Observations of disk fractions in low-metallicity clusters (Yasui et al. 2010) and statistical studies (Elsender & Bate 2021) suggest that the disk lifetime is indeed shorter for lower stellar metallicities. Nakatani et al. (2018) explain this behavior with a metallicity-dependent photo-evaporation rate. As we did not include photo-evaporation in our model, we are dealing with a different mechanism affecting the disk lifetime.

Owing to the higher stellar luminosity at lower metallicities, the innermost disk region (inside Rthick) and the outer disk regions (outside Rthick) are hotter compared to higher metallicities (see Fig. 5). With higher disk temperatures, the disk’s viscosity, ν, increases. Thus, the viscous timescale, τν = r2/ν, decreases. Since the viscous timescale is a measure of the time, in which the disk material is accreted onto the star, the disk material accretes faster onto the star for lower metallicities. As shown in Fig. C.1, in which the evolution of the stellar parameters for t ≥ t0 are summarized, the differences in intrinsic luminosity between solar and subsolar metallicity increase for older ages. As a consequence, the older a star-disk system is, the greater the difference in stellar irradiation and disk heating between solar and subsolar metallicities. Usually, the (global) viscous timescale is measured in the outer (∼10 AU) disk regions (e.g., Armitage 2011). In the few past years, however, studies have highlighted the importance of the inner disk, the sub-AU region (e.g., Vorobyov et al. 2019; Hennebelle et al. 2020). The accretion rate onto the star is eventually regulated by the interaction between the star and the disk.

To test which region (or perhaps even both) influences the disk lifetime, we artificially shielded a certain part of the disk from 25% of the stellar irradiation. This resulted in two test models, one with a lower disk temperature in the inner disk region (< 1 AU; model Linlow) and the other with a lower temperature in the outer disk region (> 1 AU; Loutlow; see panel b of Fig. 7). We compare the disk evolution of the models Linlow and Loutlow to a reference case (ref) without shielding certain regions from stellar irradiation. The evolution of the accretion rate (panel a in Fig. 7) clearly shows that the different disk temperature in the outer region influences the disk lifetime. We note that in our current model, short-term effects such as episodic outbursts are not included. The number and strength of these events depend strongly on the inner disk region (e.g., Vorobyov et al. 2020; Steiner et al. 2021) and can also affect the long-term disk evolution. Including a more detailed physical model, could increase the importance of the inner disk region. In the simulation shown in Fig. 7, the stellar evolution has been switched off to preserve computational resources. The qualitative statement of the comparison should not be impaired by this choice.

thumbnail Fig. 7.

Comparison between models Linlow (yellow line) and Loutlow (red line). The reference model (ref, blue line) corresponds to a Z = 1 Z star with start = 1 × 10−8 M yr−1. The stellar luminosity is artificially reduced by 25% in the inner 1 AU (Linlow) and outside 1 AU (Loutlow). Panel a shows the evolution of the accretion rate. The dotted black line shows the termination criterion used in our simulations of = 10−11 M yr−1. Panel b shows the temperature structure of the initial disk model. For model Linlow the inner disk region is cooler compared to the reference model, and for model Loutlow the outer disk region is cooler. The dotted black line shows the ambient temperature of 20 K.

3.4. Influence of metallicity on stellar rotation

The influence of metallicity on stellar rotation has been studied for older (≳1 Gyr) stellar populations based on the dependence of wind-breaking formulations on stellar metallicity (Amard & Matt 2020; Amard et al. 2020). Low-metallicity stars spin faster compared to their solar-metallicity counterparts. During the disk phase, however, mechanisms besides wind-braking influence the stellar rotation rate (see Sect. 2.3), and these studies assume a constant stellar period (disk-locking).

In Fig. 8 the evolution of the stellar rotation period is shown over time for both metallicities Z = 1 Z and Z = 0.2 Z. All simulations are started at a stellar period at t0 of P0 = 10 days. We expect a variation in P0 to not change the qualitative result shown in Fig. 8, but rather shift the respective lines in the vertical direction. Independent of the stellar accretion history, stars with lower stellar metallicities rotate faster compared to solar metallicity values. We can identify three reasons for this development during the disk phase.

thumbnail Fig. 8.

Stellar rotational periods over time for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. For low metallicities, stellar spin-up due to contraction (post-disk contraction, pdc) until the disk lifetime of the solar metallicity is reached is indicated by colored dotted lines.

Stellar radius. First of all, the metallicity-dependent stellar initial values have an effect on spin evolution. Low metallicities result in smaller stellar radii at t0. According to the stellar spin model described in Sect. 2.3, the two mechanisms that can remove angular momentum from the star (APSW and MEs) both scale with the stellar radius. Consequently, smaller stellar radius results in less effective stellar spin-down mechanisms and low-metallicity stars spin faster.

Stellar contraction. Contraction of the pre-MS star causes the star to spin up. Dependent on the stellar metallicity and the accretion history of the star, the evolutionary stage at a certain time and, thus, the stellar contraction can differ (see Fig. 3). Panel a in Fig. C.1 shows the evolution of the stellar radii over time for both metallicities Z = 1 Z and Z = 0.2 Z. For low metallicities, the star contracts more slowly. A slower stellar contraction rate counteracts the effect of smaller stellar radii with respect to stellar rotation (see above).

Different disk lifetime. For low metallicities, the lifetime of an accretion disk is, within the scope of our models, ∼1 Myr shorter (see Sect. 3.3). After the disk is dissolved, there are no mechanisms that can remove angular momentum and the star spins up due to contraction4. Stellar spin-up due to contraction, without the interference of a disk, starts at an earlier age. As a result, the difference in rotation period between low and solar-metallicity stars increases further.

3.5. Robustness of the results

The results presented in this section are based on certain, fixed stellar and disk parameters, for example, the stellar magnetic field strength B and the disk’s viscous α parameter. We want to verify the robustness of our results by changing these parameters. For our models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, the stellar magnetic field is reduced to 0.5 kG (model MC1), the initial rotation period is set to 2 days (model MC2), and the viscous α parameter is reduced to 0.005 (model MC3), respectively.

Panel a of Fig. 9 shows the evolution of the accretion disk (see also Fig. 6) and panel b the stellar spin evolution (see also Fig. 8) for models MC1, MC2, and MC3. Similar to the results shown in Sects. 3.3 and 3.4, the trends regarding disk lifetime and stellar rotation period with respect to the metallicity are clearly visible. This indicates that our results are valid over a wider range of parameters.

thumbnail Fig. 9.

Influence of the stellar magnetic field strength, B, initial rotation period, and viscous α parameter on the results for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The cyan, magenta, and orange lines show models MC1, MC2, and MC3, respectively. Panel a: same as Fig. 6. Panel b: same as Fig. 8 but without the post-disk contraction.

4. Discussion

4.1. Comparing our results to observed rotation periods

At ages, ≳10 Myr, the rotation periods of young stars can directly affect the high-energy radiation that interacts with planetary atmospheres (e.g., France et al. 2018). Fast rotating stars usually have stronger high-energy emissions compared to slow rotating stars (e.g., Johnstone et al. 2021). Our results indicate that low-metallicity stars rotate faster than their solar-metallicity counterparts after the disk phase (with differences of up to 2 days; see Sect. 3.4).

We compare our results to observed periods in young clusters to check the plausibility of the faster rotation periods for low-metallicity stars. Unfortunately, the age estimates of young clusters are subjected to large uncertainties, depending on the stellar evolutionary tracks used to fit the data. Additionally, within one cluster, the different stellar populations can coexist expanding the age span within one cluster up to several Myr. With these limitations in mind, we chose six young clusters with ages between ∼2 Myr and ≲5 Myr. The clusters used in this comparison are NGC 6530, NGC 2362, NGC 2264, the Orion star-forming cluster (OSFC), σ Ori, and the Taurus star-forming region, with observed rotation periods taken from Henderson & Stassun (2012), Irwin et al. (2008), Affer et al. (2013), Serna et al. (2021), Cody & Hillenbrand (2010), and Rebull et al. (2020), respectively. The observed ages, mass ranges, and metallicities of these clusters are summarized in Table E.1. We note that our models are not included in the mass range of observed periods in σ Ori and low-mass stars tend to rotate faster compared to higher stellar masses (e.g., Irwin et al. 2008). Owing to the small number of available data, we chose to include σ Ori in this comparison.

The observed rotation periods are combined in two bins: the cluster metallicity is assumed to be solar or above (PZhigh) and the metallicity is assumed to be below the solar value (PZlow). In Fig. 10, the ranges of observed rotation periods from the 25th to the 90th percentile are shown for PZhigh (light gray) and PZlow (light blue). The mean values () are indicated by horizontal lines in the respective color. Consistent with our results, the range of observed periods is shifted to faster values for lower metallicities. The difference in rotation period due to lower metallicity values, however, is small, compared to other stellar parameters, for example, the stellar magnetic field strength. We note that the lower stellar mass range of observed periods in σ Ori could further reduce the difference even further. Thus, the faster rotation periods due to lower metallicity values alone, will probably not affect the stellar high-energy output and, in consequence, the planetary atmosphere significantly. If, on the other hand, the magnetic field strength (that is assumed to be constant during our simulations) can be related to metallicity, the importance of metallicity could increase. For such relations, however, additional theoretical and observational work is needed.

thumbnail Fig. 10.

Stellar rotational periods over time for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The colored lines represent the respective model. The ranges of the observed periods of young clusters with ages between ∼3 Myr and ≲5 Myr (see text) are indicated by the light gray (light blue) area for PZhigh (PZlow). The mean values () are indicated by horizontal lines in the respective color.

4.2. Early evolution of the star–disk system: Episodic accretion at t < t0

In this study, we start the combined star-disk simulation at t0 = 2 Myr. For ages t < t0 fixed parameters are used to define the accretion history on the star (see Fig. 2). The accretion rate for t < t0, however, can vary significantly. Episodic accretion events (outbursts) triggered by gravitational instabilities and disk fragmentation, thermal and magneto-rotational instabilities, or a combination of those, can increase the accretion rate by orders of magnitude over short periods of time (e.g., Nayakshin & Lodato 2012; Bae et al. 2014; Vorobyov & Basu 2015; Vorobyov et al. 2017a).

The amount of energy added to the stellar interior depends on the accretion rate. Stronger and more numerous outbursts result in more energy added to the star. The number and intensity of outbursts also depend on the stellar mass (e.g., Vorobyov et al. 2017a; Kadam et al. 2020). High stellar masses show more numerous and luminous outbursts compared to low-mass protostars.

Furthermore, the stellar rotation period is also affected by early evolution. The principle contributions summarized in Sect. 2.3 can influence the star starting from its formation. It is interesting to note that during accretion outbursts the additional energy can cause the star to inflate (e.g., Vorobyov et al. 2017a). While the additional angular momentum carried by the accreted disk material spins up the star, the inflation of the stellar radius causes a spin-down effect. We plan, for future studies, a starting point t0 reasonably close to the formation of the stellar seed, allowing the consideration of accretion outbursts and early stellar spin evolution.

4.3. Angular momentum transport within the star

In this study, we assume solid body rotation for the star. This simplification is justified as long as the star is fully convective without the presence of a radiative core (e.g., Amard et al. 2019). When a radiative core develops together with contraction on the pre-MS, a strong meridional circulation (flowing from north to south) transports angular momentum from the core toward the surface (e.g., Amard et al. 2016, 2019).

In general, during its presence, a radiative core and the convective envelope above, which is still assumed to rotate as a solid body, have to be treated individually. The convective envelope is still affected by external torques plus the angular momentum flux induced by the meridional circulation. Additionally, the core and envelope can exchange angular momentum due to shear-induced turbulence (e.g., Amard et al. 2019). The extent of this angular momentum transfer between the core and the surface of the star depends on the stellar rotation period and stellar mass. Based on the calculations of Amard et al. (2019), either stars rotate close to solid-body rotation (valid for fast rotating stars during the first ≲108 yr) or the core rotates faster compared to the convective envelope.

Observational confirmation of such a model poses problems. Asteroseismic data reveal only a few estimates of the core angular velocity of low-mass, MS stars (Benomar et al. 2015). In contrast to the model results presented in Amard et al. (2019), these observations suggest slowly rotating stellar cores, motivating efforts to advance theoretical models and obtain more observational data.

On the other hand, helioseismic data suggest a fast-rotating solar core from g-mode observations (Fossat et al. 2017). In addition, precise Transiting Exoplanet Survey Satellite (TESS) observations of the massive star HD 192575 reveal differential rotation between the near-core region and envelope (Burssens et al. 2021). Based on their results, the near-core region rotates at a higher frequency. To what extent these results can be applied to our model, however, is unclear due to the different structures and ages of the observed stars. Furthermore, recent theoretical works of Buldgen & Eggenberger (2022) and Eggenberger et al. (2022) present progress in combining the theory of internal stellar angular momentum transport and observational constraint, for example, the solar surface lithium abundance.

4.4. Evolution of B

The stellar magnetic field strength is a highly variable parameter that can change its magnitude within a short period of time (several years; e.g., Johnstone et al. 2014). The values used in this study for B must be understood as averaged over time. Furthermore, the field topology is assumed to depend on the evolutionary stage (e.g., Zaire et al. 2017; Emeriau-Viard & Brun 2017). With the presence of a radiative core, the magnetic field topology can become more complex, which can result in a weaker large-scale dipole component. In our model, the magnetic field strength B represents the large-scale, dipole component of the magnetic field, which usually dominates the star–disk interaction region and the stellar spin evolution (Finley & Matt 2018).

For low metallicities, the star evolves further and even begins its turn toward the Henyey track, indicating the development of a radiative core. During this time, the large-scale, dipole field strength can decrease, which affects the stellar spin evolution. Referring to Fig. 9, such a decrease in B for low-metallicity stars could be visualized by a gradual transition from the blue-dashed line (where B = 2.0 kG is assumed) toward the cyan-dashed line (B = 0.5 kG in model MC1). As a result, the spread in rotational periods for different metallicity values increases further.

4.5. Disk lifetimes in low-metallicity clusters

Our models predict a dependence of disk lifetimes on stellar metallicity. The disks around low-metallicity stars dissolve faster compared to their solar metallicity counterparts (see Fig. 6). This correlation can be explained (within the scope of our model) by a hotter outer disk region enhancing viscosity, reducing the viscous timescale, and, thus, resulting in a shorter disk lifetime.

Over the last decade several observational studies can confirm this trend between metallicity and disk lifetime (e.g., Yasui et al. 2010, 2016, 2021; Guarcello et al. 2021). In general, the observations show a disk fraction in low-metallicity disks that is lower compared to solar metallicity. For example, at an age of 1 Myr, the disk fraction in a Z = 0.2 Z cluster is ∼10% (e.g., Yasui et al. 2010) and considerably lower compared to solar metallicity values with disk fractions at 1 Myr of 60 − 80%.

Different mechanisms affect the dispersal of the accretion disk in low-metallicity environments. The ionization fraction in low-metallicity disks can be higher compared to solar metallicities. Smaller amounts of dust result in a slower recombination rate of free charges on dust grain surfaces. A higher ionization fraction in the disk increases the efficiency of magneto-rotational instabilities (Balbus & Hawley 1991), which drive accretion (e.g., Hartmann 2009). Alternatively, the efficiency of photo-evaporation also depends on metallicity. Lower metallicities and the resulting lower disk opacities allow a deeper penetration of high-energy radiation into the disk (e.g., Ercolano & Clarke 2010). As a result, a larger part of the disk is heated and the removal of gas and dust by photo-evaporation is facilitated. More recently, numerical simulations (e.g., Nakatani et al. 2018) confirm the metallicity dependence of photo-evaporation. For subsolar metallicities, the photo-evaporative mass-loss rates are higher and the accretion disk disperses at a faster rate. It is worth noting that the observed disk lifetimes in low-metallicity clusters cannot be conclusively explained by one mechanism alone (e.g., Nakatani et al. 2018) and results, more likely, from the combination of the aforementioned aspects. In a more recent study, Elsender & Bate (2021) propose that increased stellar multiplicity in low-metallicity environments and resulting increased dynamical interactions might cause accretion disks to be smaller and have shorter lifetimes compared to those with higher metallicities.

4.6. Differences between our results and observations

The observations shown in Yasui et al. (2010) suggest disk lifetimes of ∼3 Myr for low-metallicity clusters, matching our results for start = 1 × 10−9 M yr−1. In this context, cluster ages can be used to estimate disk lifetimes. The disk lifetimes for models with start = 5 × 10−9 yr−1 and 1 × 10−8M yr−1, on the other hand, are longer compared to observations by a factor of ∼3. We want to address three possible reasons for this discrepancy. First, our models do currently not take photoevaporation into account. Including photoevaporation in future models will reduce our model’s lifetime toward the observed values. Second, the cluster ages presented by Yasui et al. (2010) are based on models that do not take magnetic fields into account. Cluster ages based on magnetic pre-MS stellar evolution models (e.g., Feiden et al. 2015) are up to 2 Myr older compared to previous models (e.g., Richert et al. 2018). Disk lifetimes for low-metallicity stars could increase up to ∼5 Myr. Third, the initial conditions in our models are independent of metallicity. It is, however, likely that the accretion rate, stellar rotation period, and stellar magnetic field strength depend on metallicity. Moving t0 to earlier evolutionary stages in combination with additional observational constraints for low-metallicity clusters will improve the quality of our results in future work.

4.7. Previous spin evolution studies that include metallicity

Many new formulations have emerged over the last decade to explain the stellar spin evolution (e.g., Matt et al. 2015; Finley & Matt 2018; Gallet et al. 2019; Ireland et al. 2021). These studies, however, are often based on solar properties and composition. Thus, an explicit dependence on stellar metallicity is not given.

Amard et al. (2019) presents, for the first time, a comprehensive set of stellar evolutionary tracks including rotation and different metallicities, using the STAREVOL code. They find that metallicity does influence the stellar spin evolution during the pre-MS and MS. Low-metallicity stars are more compact and evolve faster compared to their solar-metallicity counterparts. Thus, during the pre-MS, they spin up faster and retain faster rotation periods during the MS. Internally, low-metallicity stars experience less differential rotation compared to solar metallicities. A more recent study Amard & Matt (2020) provides an in-depth analysis of the effects of metallicity on different wind-braking formulations. They also conclude that low-metallicity stars rotate faster than their solar-metallicity counterparts. During the disk phase, a constant stellar period (disk-locking) is assumed for a fixed timescale independent of stellar mass and metallicity (Amard et al. 2019; Amard & Matt 2020). This assumption is in conflict with observations of metallicity-dependent disk lifetimes (e.g., Yasui et al. 2010, 2016, 2021; Guarcello et al. 2021), photo-evaporation models (e.g., Nakatani et al. 2018) and the results presented in Sect. 3.3.

Starting from the T Tauri phase, our star-disk evolution model does include the effects of, for example, metallicity, stellar magnetic field strength, stellar spin, and mass. Over a wide range of stellar and disk parameters, we can provide realistic, metallicity-dependent disk lifetimes and stellar periods after the disk phase. Those values can then be used as input parameters for the aforementioned spin evolution studies.

4.8. Impact on disk chemistry

Various studies confirm the influence of metallicity on the chemical evolution of disks. For example, by observing the chemical abundances of HII regions in 12 nearby dwarf galaxies, Berg et al. (2016) suggest that there are different chemistry cycles in non-solar-metallicity environments. This evidence hints at the fact that disk evolution and planet formation in protoplanetary disks may proceed differently in diverse metallicity environments. Certain aspects of the physical conditions considered by our model could change the way metallicity influences chemistry.

One of the important features is the usage of passive or active disks. For passive disks, the only heating source is stellar radiation and cosmic rays. Therefore, passive disks exhibit an increase in temperature with decreasing metallicity. In the case of active disks, the hydrodynamical processes provide additional heating via viscous and compressional heating, which operate predominantly in the disk midplane and depend on the assumed metallicity. A decrease in metallicity leads to a lower temperature in the midplane (see Fig. 5).

In a recent study, Guadarrama et al. (2022) showed that snow lines of chemical species are pushed back as metallicity decreases. Assuming a passive disk, their model does not take viscous heating into account resulting in an increasing disk temperature throughout the disk with decreasing metallicity. This assumption might be appropriate for outer disk regions ≳10 AU, where the disk temperature is dominated by stellar irradiation. Within Rthick, however, the disk midplane temperature decreases with decreasing metallicity (see Sect. 3.2) and snow lines of species with sublimation temperatures ≳50 K (e.g., H20 and CO2) are shifted inward. As a consequence, dust grows and the radial chemical composition of the disk is affected. We note that our model starts at an already advanced evolutionary disk stage. To study the full impact of the aforementioned temperature structures for different metallicities on dust growth and planetesimal formation, the starting point of our model t0 has to be shifted toward earlier stages.

4.9. Model limitations

Although several important aspects of stellar and disk evolution are combined in this study, the results are still preliminary due to certain limitations of our model:

One-dimensional model. The main limitation of our disk model is the one-dimensionality. Due to timestep limitations, long-term evolution models that include the innermost disk regions are currently restricted to one dimension (e.g., Steiner et al. 2021). Unfortunately, several important stellar and disk properties require (at least) an extension in the vertical direction. Accretion via funnels, the interaction between stellar magnetic fields and the disk, vertical disk structure, and heating of the disk are intrinsic multi-dimensional processes. In the current version of our model, these aspects are included in a simplistic way. Nevertheless, the future goal can only be to expand our model by (at least) another dimension.

Viscous accretion disk. The transport of mass and angular momentum in our model is described by a simple (constant α) viscosity prescription. α most certainly varies with respect to time and location (e.g., Gammie 1996; King et al. 2007; Vorobyov et al. 2020), and other mechanisms such as disk magneto-centrifugal winds, photo-evaporation and the combination of both also affect the disk evolution (e.g., Königl et al. 2011; Roquette et al. 2021; Kunitomo et al. 2020). These are currently not included but will be added in some future versions.

Stellar multiplicity and close encounters. In the current model, only single stars are considered. Binary or even multiple companions can also influence the accretion disk (e.g., Messina et al. 2017). Additionally, a fly-by or close encounter with a stellar object can impact the disk evolution (e.g., Vorobyov et al. 2017b). These effects, however, are naturally non-axisymmetric events and are currently not feasible with our model.

Variable dust component. The dust-to-gas mass ratio, fdust, in this study is kept constant at the respective stellar metallicity value trough out the simulation and over the disk’s radial extent. While this assumption might be reasonable in the first step, the dust fraction can differ within the disk over time due to a great number of processes, such as coagulation, drift, fragmentation, and vertical settling (e.g., Testi et al. 2014; Vorobyov et al. 2022). Our results indicate a direct connection between the amount of dust in the disk and its evolution and a more sophisticated dust evolution model will improve our results.

5. Conclusion

In this study we combine three important aspects of a young star: the stellar evolution with MESA, the disk evolution with TAPIR, and the stellar spin evolution based on Gallet et al. (2019). With reference to the influence of stellar metallicity on the evolution of the accretion disk and stellar spin during the star-disk phase, the main conclusions are summarized as follows:

  • The accretion history and metallicity can influence the evolutionary stage of a star. A low-metallicity star is more compact and luminous. The higher the initial accretion rate, init, the more energy is added to the stellar interior during the early stellar evolution, t < t0, and the larger the stellar radius. On the HRD, the star appears less evolved compared to when there are small values of init (see Figs. 2 and 3).

  • Larger luminosities of low-metallicity stars affect the radial temperature profile of the accretion disk. The disk midplane temperature of the innermost disk region, ≲0.1 AU, and the disk beyond ≳1 AU are dominated by stellar irradiation and thus hotter compared to solar-metallicity disks. Within Rthick, however, the disk midplane is shielded from stellar irradiation and the temperature is dominated by dust opacities. High stellar metallicities result in higher opacity values, and the temperatures in the disk midplane are higher compared to lower stellar metallicities (see Fig. 5).

  • The lifetime of an accretion disk is influenced by metallicity. The increased temperatures in the outer disk regions (see Fig. 5) result in a more effective (faster) accretion onto the star, and the lifetimes are 0.6−1.4 Myr shorter for low metallicities, within the scope of our model (see Fig. 6). We present an additional explanation for the observed short disk lifetimes in young, low-metallicity clusters (Yasui et al. 2010). A metallicity-dependent photo-evaporation description (e.g., Nakatani et al. 2018) arrives at similar results. Their lifetimes, however, are still too long to match the observations. A combination of our star-disk evolution model, coupled with a metallicity-dependent photo-evaporation, could be a key step toward understanding disk lifetimes in different metallicity environments.

  • Similar to their older counterparts (e.g., Amard & Matt 2020; Amard et al. 2020), young, low-metallicity stars rotate faster compared to those with solar metallicities (see Fig. 8). We can identify three mechanisms that affect stellar spin evolution. For low metallicities, stellar radii are smaller, and thus the mechanisms that can remove angular momentum from the star are less effective (see Sect. 2.3). Additionally, the stellar contraction rate depends on the stellar evolutionary stage. The stellar contraction rate is smaller if a star is more evolved. As a consequence, low-metallicity stars contract more slowly (see Fig. C.1) and the spin-up due to contraction is reduced. Finally, due to the shorter lifetime of low-metallicity disks, stellar spin-up due to contraction, without the interference of a disk, starts at an earlier age. Thus, the difference in rotation period between solar- and subsolar-metallicity stars increases further.

  • Our results, which take into account the effects of metallicity, stellar magnetic field strength, and stellar mass in the calculation of the disk lifetime and stellar rotation period during the disk phase, can be used as input parameters for other recent spin evolution studies (e.g., Amard et al. 2019; Amard & Matt 2020; Gossage et al. 2021).

Our model clearly indicates the potential of jointly studying the hydrodynamic evolution of the star and the disk. We recommend a step-by-step inclusion of additional physical aspects (e.g., photo-evaporation and magneto-centrifugal disk winds) and studies that cover wider parameter spaces to improve our understanding of these young star-disk systems.


1

We note that in more recent studies of MEs (e.g., Pantolmos et al. 2020; Ireland et al. 2021) smaller values for Krot are found. Krot = 0.7 can be understood as an upper limit (Ireland et al. 2021).

2

The different behavior at the beginning of the evolution (the first ∼100 yr) can be explained by the outer boundary condition of the MESA code that reacts differently on different values of Ladd, which depends on the accretion rate.

3

In this context, disk lifetime is the time at which the accretion rate has dropped to 10−11M yr−1.

4

We note that there is also a stellar wind mechanism acting after the disk phase. Compared to stellar contraction, however, its effect is small and can be neglected during a relatively short time period of ∼1 Myr. After contraction has ended and the star evolves toward and on the MS, stellar winds slow down the star on longer timescales (e.g., Gallet & Bouvier 2013).

Acknowledgments

The authors would like to thank the anonymous referee, who provided very detailed comments that helped to improve the quality of the manuscript. Furthermore, the work of Bill Paxton and his collaborators on the stellar evolution code MESA is gratefully acknowledged. E.I.V. and R.G. acknowledge support by the Austrian Science Fund (FWF) under research grant P31635-N27.

References

  1. Affer, L., Micela, G., Favata, F., Flaccomio, E., & Bouvier, J. 2013, MNRAS, 430, 1433 [Google Scholar]
  2. Amard, L., & Matt, S. P. 2020, ApJ, 889, 108 [Google Scholar]
  3. Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, A&A, 587, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Amard, L., Roquette, J., & Matt, S. P. 2020, MNRAS, 499, 3481 [NASA ADS] [CrossRef] [Google Scholar]
  6. Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241 [Google Scholar]
  7. Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [Google Scholar]
  8. Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [NASA ADS] [CrossRef] [Google Scholar]
  9. Armitage, P. J. 2011, ARA&A, 49, 195 [Google Scholar]
  10. Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61 [NASA ADS] [CrossRef] [Google Scholar]
  12. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
  13. Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJ, 702, L27 [NASA ADS] [CrossRef] [Google Scholar]
  14. Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118 [NASA ADS] [CrossRef] [Google Scholar]
  15. Barnes, S. A. 2007, ApJ, 669, 1167 [Google Scholar]
  16. Bastien, F. A., Stassun, K. G., & Weintraub, D. A. 2011, AJ, 142, 141 [NASA ADS] [CrossRef] [Google Scholar]
  17. Benomar, O., Takata, M., Shibahashi, H., Ceillier, T., & García, R. A. 2015, MNRAS, 452, 2654 [Google Scholar]
  18. Berg, D. A., Skillman, E. D., Henry, R. B. C., Erb, D. K., & Carigi, L. 2016, ApJ, 827, 126 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bessolaz, N., Zanni, C., Ferreira, J., Keppens, R., & Bouvier, J. 2008, A&A, 478, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Bhandare, A., Kuiper, R., Henning, T., et al. 2018, A&A, 618, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Buchler, J. R., & Yueh, W. R. 1976, ApJ, 210, 440 [NASA ADS] [CrossRef] [Google Scholar]
  22. Buldgen, G., & Eggenberger, P. 2022, in Proceedings of the Sixteenth Marcel Grossman meeting - July 5–10 2021, - Session on “Rotation in Stellar Evolution” [Google Scholar]
  23. Burssens, S., Bowman, D. M., Michielsen, M., Simón-Díaz, S., & Aerts, C. 2021, https://doi.org/10.5281/zenodo.5126788 [Google Scholar]
  24. Caffau, E., Ludwig, H. G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255 [Google Scholar]
  25. Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094 [NASA ADS] [CrossRef] [Google Scholar]
  26. Cauley, P. W., Johns-Krull, C. M., Hamilton, C. M., & Lockhart, K. 2012, ApJ, 756, 68 [NASA ADS] [CrossRef] [Google Scholar]
  27. Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cody, A. M., & Hillenbrand, L. A. 2010, ApJS, 191, 389 [NASA ADS] [CrossRef] [Google Scholar]
  29. Cranmer, S. R. 2008, ApJ, 689, 316 [NASA ADS] [CrossRef] [Google Scholar]
  30. Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240 [NASA ADS] [CrossRef] [Google Scholar]
  31. Damiani, F., Prisinzano, L., Micela, G., & Sciortino, S. 2019, A&A, 623, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. D’Orazi, V., Randich, S., Flaccomio, E., et al. 2009, A&A, 501, 973 [CrossRef] [EDP Sciences] [Google Scholar]
  33. D’Orazi, V., Biazzo, K., & Randich, S. 2011, A&A, 526, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Dorfi, E., & Drury, L. 1987, J. Comput. Phys., 69, 175 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dorfi, E., Pikall, H., Stökl, A., & Gautschy, A. 2006, Comput. Phys. Commun., 174, 771 [NASA ADS] [CrossRef] [Google Scholar]
  36. Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52 [NASA ADS] [CrossRef] [Google Scholar]
  37. Eggenberger, P., Buldgen, G., Salmon, S. J. A. J., et al. 2022, Nat. Astron., 6, 788 [NASA ADS] [CrossRef] [Google Scholar]
  38. Elbakyan, V. G., Vorobyov, E. I., Rab, C., et al. 2019, MNRAS, 484, 146 [NASA ADS] [CrossRef] [Google Scholar]
  39. Elsender, D., & Bate, M. R. 2021, MNRAS, 508, 5279 [NASA ADS] [CrossRef] [Google Scholar]
  40. Emeriau-Viard, C., & Brun, A. S. 2017, ApJ, 846, 8 [Google Scholar]
  41. Ercolano, B., & Clarke, C. J. 2010, MNRAS, 402, 2735 [Google Scholar]
  42. Feiden, G. A., Jones, J., & Chaboyer, B. 2015, 18th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 18, 171 [NASA ADS] [Google Scholar]
  43. Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [Google Scholar]
  44. Finley, A. J., & Matt, S. P. 2018, ApJ, 854, 78 [Google Scholar]
  45. Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109 [NASA ADS] [CrossRef] [Google Scholar]
  46. Fossat, E., Boumier, P., Corbard, T., et al. 2017, A&A, 604, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. France, K., Arulanantham, N., Fossati, L., et al. 2018, ApJS, 239, 16 [Google Scholar]
  48. Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1 [NASA ADS] [CrossRef] [Google Scholar]
  49. Gallet, F., & Bouvier, J. 2013, A&A, 556, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Gallet, F., Zanni, C., & Amard, L. 2019, A&A, 632, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Gammie, C. F. 1996, ApJ, 457, 355 [Google Scholar]
  52. Gehrig, L., Steiner, D., Vorobyov, E. I., & Güdel, M. 2022, A&A, 667, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Gossage, S., Dotter, A., Garraffo, C., et al. 2021, ApJ, 912, 65 [NASA ADS] [CrossRef] [Google Scholar]
  54. Guadarrama, R., Vorobyov, E. I., Rab, C., & Güdel, M. 2022, A&A, 667, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Guarcello, M. G., Biazzo, K., Drake, J. J., et al. 2021, A&A, 650, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323 [Google Scholar]
  57. Hartmann, L. 2009, Accretion Processes in Star Formation: Second Edition (Cambridge: Cambridge University Press) [Google Scholar]
  58. Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, ApJ, 475, 770 [NASA ADS] [CrossRef] [Google Scholar]
  59. Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135 [Google Scholar]
  60. Hayashi, C. 1961, PASJ, 13, 450 [NASA ADS] [Google Scholar]
  61. Henderson, C. B., & Stassun, K. G. 2012, ApJ, 747, 51 [Google Scholar]
  62. Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140 [NASA ADS] [CrossRef] [Google Scholar]
  64. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  65. Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752 [Google Scholar]
  66. Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
  67. Ireland, L. G., Zanni, C., Matt, S. P., & Pantolmos, G. 2021, ApJ, 906, 4 [Google Scholar]
  68. Irwin, J., Hodgkin, S., Aigrain, S., et al. 2008, MNRAS, 384, 675 [CrossRef] [Google Scholar]
  69. Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411 [NASA ADS] [CrossRef] [Google Scholar]
  70. Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 474, 1176 [Google Scholar]
  71. Johnstone, C. P., Jardine, M., Gregory, S. G., Donati, J. F., & Hussain, G. 2014, MNRAS, 437, 3202 [Google Scholar]
  72. Johnstone, C. P., Bartel, M., & Güdel, M. 2021, A&A, 649, A96 [EDP Sciences] [Google Scholar]
  73. Kadam, K., Vorobyov, E., Regály, Z., Kóspál, Á., & Ábrahám, P. 2020, ApJ, 895, 41 [NASA ADS] [CrossRef] [Google Scholar]
  74. King, A. R., Pringle, J. E., & Livio, M. 2007, MNRAS, 376, 1740 [NASA ADS] [CrossRef] [Google Scholar]
  75. Kluźniak, W., & Rappaport, S. 2007, ApJ, 671, 1990 [CrossRef] [Google Scholar]
  76. Koenigl, A. 1991, ApJ, 370, L39 [Google Scholar]
  77. Koldoba, A. V., Lovelace, R. V. E., Ustyugova, G. V., & Romanova, M. M. 2002, AJ, 123, 2019 [NASA ADS] [CrossRef] [Google Scholar]
  78. Königl, A., Romanova, M. M., & Lovelace, R. V. E. 2011, MNRAS, 416, 757 [NASA ADS] [Google Scholar]
  79. Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A, 599, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Kunitomo, M., Suzuki, T. K., & Inutsuka, S.-I. 2020, MNRAS, 492, 3849 [CrossRef] [Google Scholar]
  81. Langanke, K., & Martínez-Pinedo, G. 2000, Nucl. Phys. A, 673, 481 [NASA ADS] [CrossRef] [Google Scholar]
  82. Larson, R. B. 1969, MNRAS, 145, 271 [Google Scholar]
  83. Lavail, A., Kochukhov, O., Hussain, G. A. J., et al. 2017, A&A, 608, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Lavail, A., Kochukhov, O., & Hussain, G. A. J. 2019, A&A, 630, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Manara, C. F., Robberto, M., Rio, N. D., et al. 2012, ApJ, 755, 154 [NASA ADS] [CrossRef] [Google Scholar]
  86. Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2022, in Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [Google Scholar]
  88. Matt, S., & Pudritz, R. E. 2005, ApJ, 632, L135 [Google Scholar]
  89. Matt, S. P., Pinzón, G., de la Reza, R., & Greene, T. P. 2010, ApJ, 714, 989 [NASA ADS] [CrossRef] [Google Scholar]
  90. Matt, S. P., Pinzón, G., Greene, T. P., & Pudritz, R. E. 2012, ApJ, 745, 101 [NASA ADS] [CrossRef] [Google Scholar]
  91. Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23 [Google Scholar]
  92. Mercer-Smith, J. A., Cameron, A. G. W., & Epstein, R. I. 1984, ApJ, 279, 363 [NASA ADS] [CrossRef] [Google Scholar]
  93. Messina, S., Lanzafame, A. C., Malo, L., et al. 2017, A&A, 607, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Micela, G., Sciortino, S., Serio, S., et al. 1985, ApJ, 292, 172 [CrossRef] [Google Scholar]
  95. Mulders, G. D., & Dominik, C. 2012, A&A, 539, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018, ApJ, 857, 57 [NASA ADS] [CrossRef] [Google Scholar]
  97. Nayakshin, S., & Lodato, G. 2012, MNRAS, 426, 70 [Google Scholar]
  98. Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, At. Data Nucl. Data Tab., 56, 231 [NASA ADS] [CrossRef] [Google Scholar]
  99. Palla, F., & Stahler, S. W. 1992, ApJ, 392, 667 [NASA ADS] [CrossRef] [Google Scholar]
  100. Pallavicini, R., Golub, L., Rosner, R., et al. 1981, ApJ, 248, 279 [NASA ADS] [CrossRef] [Google Scholar]
  101. Pantolmos, G., Zanni, C., & Bouvier, J. 2020, A&A, 643, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  103. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  104. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  105. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  106. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  107. Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [Google Scholar]
  108. Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Pollack, J. B., McKay, C. P., & Christofferson, B. M. 1985, Icarus, 64, 471 [NASA ADS] [CrossRef] [Google Scholar]
  110. Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964 [Google Scholar]
  111. Potekhin, A. Y., & Chabrier, G. 2010, Contrib. Plasma Phys., 50, 82 [NASA ADS] [CrossRef] [Google Scholar]
  112. Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
  113. Ragossnig, F., Dorfi, E. A., Ratschiner, B., et al. 2020, Comput. Phys. Commun., 256, 107437 [NASA ADS] [CrossRef] [Google Scholar]
  114. Rappaport, S. A., Fregeau, J. M., & Spruit, H. 2004, ApJ, 606, 436 [NASA ADS] [CrossRef] [Google Scholar]
  115. Rebull, L. M., Stauffer, J. R., Cody, A. M., et al. 2020, AJ, 159, 273 [Google Scholar]
  116. Richert, A. J. W., Getman, K. V., Feigelson, E. D., et al. 2018, MNRAS, 477, 5191 [NASA ADS] [CrossRef] [Google Scholar]
  117. Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
  118. Romanova, M., & Kurosawa, R. 2014, in 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013), eds. N. V. Pogorelov, E. Audit, & G. P. Zank, ASP Conf. Ser., 488, 127 [Google Scholar]
  119. Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2002, ApJ, 578, 420 [Google Scholar]
  120. Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2004, ApJ, 610, 920 [Google Scholar]
  121. Roquette, J., Matt, S. P., Winter, A. J., Amard, L., & Stasevic, S. 2021, MNRAS, 508, 3710 [NASA ADS] [CrossRef] [Google Scholar]
  122. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  123. Schib, O., Mordasini, C., Wenger, N., Marleau, G. D., & Helled, R. 2021, A&A, 645, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Serna, J., Hernandez, J., Kounkel, M., et al. 2021, ApJ, 923, 177 [NASA ADS] [CrossRef] [Google Scholar]
  126. Shakura, N., & Sunyaev, R. 1973, A&A, 24, 337 [Google Scholar]
  127. Sherry, W. H., Walter, F. M., Wolk, S. J., & Adams, N. R. 2008, AJ, 135, 1616 [NASA ADS] [CrossRef] [Google Scholar]
  128. Sicilia-Aguilar, A., Banzatti, A., Carmona, A., et al. 2016, PASA, 33, e059 [NASA ADS] [CrossRef] [Google Scholar]
  129. Stahler, S. W. 1988, ApJ, 332, 804 [NASA ADS] [CrossRef] [Google Scholar]
  130. Stahler, S. W., & Palla, F. 2004, The Formation of Stars (Wiley-VCH) [CrossRef] [Google Scholar]
  131. Steindl, T., Zwintz, K., Barnes, T. G., Müllner, M., & Vorobyov, E. I. 2021, A&A, 654, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  132. Steindl, T., Zwintz, K., & Vorobyov, E. 2022, Nat. Commun., 13, 5355 [NASA ADS] [CrossRef] [Google Scholar]
  133. Steiner, D., Gehrig, L., Ratschiner, B., et al. 2021, A&A, 655, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  134. Tadross, A. 2003, New Astron., 8, 737 [NASA ADS] [CrossRef] [Google Scholar]
  135. Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339 [Google Scholar]
  136. Testi, L., Natta, A., Manara, C. F., et al. 2022, A&A, 663, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
  138. Vorobyov, E. I., & Basu, S. 2009, MNRAS, 393, 822 [CrossRef] [Google Scholar]
  139. Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896 [NASA ADS] [CrossRef] [Google Scholar]
  140. Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115 [NASA ADS] [CrossRef] [Google Scholar]
  141. Vorobyov, E. I., Elbakyan, V., Hosokawa, T., et al. 2017a, A&A, 605, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  142. Vorobyov, E. I., Steinrueck, M. E., Elbakyan, V., & Guedel, M. 2017b, A&A, 608, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  143. Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  144. Vorobyov, E. I., Khaibrakhmanov, S., Basu, S., & Audard, M. 2020, A&A, 644, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Vorobyov, E. I., Skliarevskii, A. M., Molyarova, T., et al. 2022, A&A, 658, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  146. Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48 [Google Scholar]
  147. Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27 [NASA ADS] [CrossRef] [Google Scholar]
  148. Yasui, C., Kobayashi, N., Tokunaga, A. T., Saito, M., & Tokoku, C. 2010, ApJ, 723, L113 [NASA ADS] [CrossRef] [Google Scholar]
  149. Yasui, C., Kobayashi, N., Saito, M., & Izumi, N. 2016, AJ, 151, 115 [NASA ADS] [CrossRef] [Google Scholar]
  150. Yasui, C., Kobayashi, N., Saito, M., Izumi, N., & Skidmore, W. 2021, AJ, 161, 139 [NASA ADS] [CrossRef] [Google Scholar]
  151. Zaire, B., Guerrero, G., Kosovichev, A. G., Smolarkiewicz, P. K., & Landin, N. R. 2017, in Living Around Active Stars, eds. D. Nandy, A. Valio, & P. Petit, 328, 30 [NASA ADS] [Google Scholar]
  152. Zanni, C., & Ferreira, J. 2013, A&A, 550, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  153. Zhu, Z., Hartmann, L., Calvet, N., et al. 2007, ApJ, 669, 483 [NASA ADS] [CrossRef] [Google Scholar]
  154. Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009, ApJ, 701, 620 [CrossRef] [Google Scholar]

Appendix A: MESA microphysics

The MESA equation of state is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), PTEH (Pols et al. 1995), HELM (Timmes & Swesty 2000), and PC (Potekhin & Chabrier 2010) equations of state.

Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005) and the high-temperature, Compton-scattering dominated regime by Buchler & Yueh (1976). Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are a combination of rates from NACRE (Angulo et al. 1999) and JINA REACLIB (Cyburt et al. 2010), plus additional tabulated weak reaction rates (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996).

Appendix B: Treatment of rotation in MESA

In MESA, rotation is treated as a modification of the stellar equations due to centrifugal acceleration and nonspherical symmetry, when the star is rotating (see Sect. 6 in Paxton et al. 2013). Angular momentum transport due to rotation (rotational mixing) is implemented as a diffusion term, including five rotationally induced mixing processes: dynamical shear instability, Solberg-Høiland instability, secular shear instability, Eddington-Sweet circulation, and the Goldreich-Schubert-Fricke instability (Paxton et al. 2013).

Appendix C: Stellar evolution: t > t0

Fig. C.1 shows the evolution of the stellar radius (Panel a), intrinsic luminosity (Panel b), and mass (Panel c) for Z = 1 Z and Z = 0.2 Z.

thumbnail Fig. C.1.

Evolution of the stellar radius (Panel a), intrinsic luminosity (Panel b), and mass (Panel c) for Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively.

Appendix D: Time evolution of radial disk structure

The radial profiles for the model with start = 1.0 × 10−8 M yr−1 at the starting point of the simulation, t0 = 2 Myr, as well as at 5.0 Myr and 8.0 Myr, are shown in Fig. D.1.

thumbnail Fig. D.1.

Radial profiles for the model with start = 1.0 × 10−8 M yr−1 at the starting point of the simulation, t0 = 2 Myr (left column), as well as at 5.0 Myr (center column) and 8.0 Myr (right column). The solid (dashed) lines correspond to a metallicity of Z = 1 Z (Z = 0.2 Z). Panel (a) shows the disk surface density, Σ, panel (b) the disk midplane temperature, Tdisk, panel (c) the optical depth, τ, and panel (d) the mass flux flux. The thin vertical lines in panel (a) indicate the respective position of rc. The horizontal lines in panel (b) indicate the ambient temperature Tamb = 20 K. The horizontal dashed lines in panel (c) symbolize τ = 1, dividing optically thin regions (τ ≪ 1) from optically thick regions (τ ≫ 1).

Appendix E: Age, stellar mass range, and metallicity of young clusters

Table E.1.

Age, stellar mass range, and metallicity of young (∼2 − 5 Myr) clusters used in Sect. 4.1.

All Tables

Table 1.

Onset of deuterium burning (τD), the stellar mass (MH), and evolution time (τH) when the star turns onto the Hayashi track.

Table 2.

Characteristic radius, Rc, and disk mass, mdisk, of the initial disk model.

Table E.1.

Age, stellar mass range, and metallicity of young (∼2 − 5 Myr) clusters used in Sect. 4.1.

All Figures

thumbnail Fig. 1.

Schematic representation of the interaction between TAPIR and MESA. Starting from an initial (old) timestep, t0, the disk, spin, and star values are evolved toward the new timestep, t1 = t0 + δt. The new values at t1 are obtained in three steps, represented by the numbered arrows (see text).

In the text
thumbnail Fig. 2.

Initial stellar evolution for t < t0. Panel a: accretion history of the star. Initially, the star accretes with a constant accretion rate, init. At time tdec, the accretion rate drops according to t−2 (see text). The colored dashed lines show the respective value of start. Panel b: evolution of the stellar mass. Starting from the initial mass of 0.01 M, all models reach their final mass of 0.91 M after 2 Myr. The yellow, red, and blue lines show models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. To show the entire initial stellar evolution, starting from t = 0 yr, we use the symlog representation (Hunter 2007). The vertical dashed line indicates the transition from the linear scale (for t ≤ 104 yr) to the logarithmic scale (for t > 104 yr).

In the text
thumbnail Fig. 3.

Evolution on the HRD for t < t0 for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The values at time t = 0 Myr (t = 2 Myr) are marked with circles (triangles) and the yellow, red, and blue lines show models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. The onset of deuterium burning (t = τD) is marked with a diamond. The transition from expansion toward the Hayashi track (t = τH) is marked with a cross.

In the text
thumbnail Fig. 4.

Stellar radii and temperatures at t0 = 2 Myr for start = 5 × 10−9 (red triangles) and 1 × 10−8 (blue triangles). Solar (subsolar) metallicity values are connected by a solid (dashed) line.

In the text
thumbnail Fig. 5.

Radial profile of the initial disk model for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. Panel a shows the disk surface density, Σ, panel b the disk midplane temperature, Tdisk, and panel c the optical depth, τ. The thin vertical lines in panel a indicate the respective position of Rc. The horizontal line in panel b indicates the ambient temperature Tamb = 20 K. The horizontal dashed line in panel c symbolizes τ = 1, which divides optically thin regions (τ ≪ 1) from optically thick regions (τ ≫ 1).

In the text
thumbnail Fig. 6.

Accretion rate onto the star over time for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. The dotted black line shows the termination criterion used in our simulations of = 10−11 M yr−1.

In the text
thumbnail Fig. 7.

Comparison between models Linlow (yellow line) and Loutlow (red line). The reference model (ref, blue line) corresponds to a Z = 1 Z star with start = 1 × 10−8 M yr−1. The stellar luminosity is artificially reduced by 25% in the inner 1 AU (Linlow) and outside 1 AU (Loutlow). Panel a shows the evolution of the accretion rate. The dotted black line shows the termination criterion used in our simulations of = 10−11 M yr−1. Panel b shows the temperature structure of the initial disk model. For model Linlow the inner disk region is cooler compared to the reference model, and for model Loutlow the outer disk region is cooler. The dotted black line shows the ambient temperature of 20 K.

In the text
thumbnail Fig. 8.

Stellar rotational periods over time for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively. For low metallicities, stellar spin-up due to contraction (post-disk contraction, pdc) until the disk lifetime of the solar metallicity is reached is indicated by colored dotted lines.

In the text
thumbnail Fig. 9.

Influence of the stellar magnetic field strength, B, initial rotation period, and viscous α parameter on the results for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The cyan, magenta, and orange lines show models MC1, MC2, and MC3, respectively. Panel a: same as Fig. 6. Panel b: same as Fig. 8 but without the post-disk contraction.

In the text
thumbnail Fig. 10.

Stellar rotational periods over time for both metallicities Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The colored lines represent the respective model. The ranges of the observed periods of young clusters with ages between ∼3 Myr and ≲5 Myr (see text) are indicated by the light gray (light blue) area for PZhigh (PZlow). The mean values () are indicated by horizontal lines in the respective color.

In the text
thumbnail Fig. C.1.

Evolution of the stellar radius (Panel a), intrinsic luminosity (Panel b), and mass (Panel c) for Z = 1 Z (solid lines) and Z = 0.2 Z (dashed lines). The yellow, red, and blue lines indicate models with start = [1 × 10−9, 5 × 10−9, and 1 × 10−8] M yr−1, respectively.

In the text
thumbnail Fig. D.1.

Radial profiles for the model with start = 1.0 × 10−8 M yr−1 at the starting point of the simulation, t0 = 2 Myr (left column), as well as at 5.0 Myr (center column) and 8.0 Myr (right column). The solid (dashed) lines correspond to a metallicity of Z = 1 Z (Z = 0.2 Z). Panel (a) shows the disk surface density, Σ, panel (b) the disk midplane temperature, Tdisk, panel (c) the optical depth, τ, and panel (d) the mass flux flux. The thin vertical lines in panel (a) indicate the respective position of rc. The horizontal lines in panel (b) indicate the ambient temperature Tamb = 20 K. The horizontal dashed lines in panel (c) symbolize τ = 1, dividing optically thin regions (τ ≪ 1) from optically thick regions (τ ≫ 1).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.