11institutetext:  Department of Astrophysics, University of Vienna,
Türkenschanzstrasse 17, A-1180 Vienna, Austria 22institutetext: Institut für
Astro- und Teilchenphysik, Universität Innsbruck, Technikerstraße 25, A-6020
Innsbruck, Austria

# The influence of metallicity on a combined stellar and disk evolution

L. Gehrig 11 T. Steindl 22 E. I. Vorobyov 11 R. Guadarrama 11 K. Zwintz 22

(Received ….; accepted ….)

###### Abstract

Context. The effects of an accretion disk are crucial to understanding the
evolution of young stars. During the combined evolution, stellar and disk
parameters influence each other, motivating a combined stellar and disk model.
This makes a combined numerical model, evolving the disk alongside the star,
the next logical step in the progress of studying early stellar evolution.

Aims. We aim to understand the effects of metallicity on the accretion disk
and the stellar spin evolution during the T Tauri phase.

Methods. We combine the numerical treatment of a hydrodynamic disk with
stellar evolution, including a stellar spin model, allowing a self-consistent
calculation of the back-reactions between the individual components.

Results.  We present the self-consistent theoretical evolution of T-Tauri
stars coupled to a stellar disk. We find that disks in low metallicity
environments are heated differently and have shorter lifetimes, compared to
their solar metallicity counterparts. Differences in stellar radii, the
contraction rate of the stellar radius, and the shorter disk lifetimes result
in faster rotation of low metallicity stars.

Conclusions. We present an additional explanation for the observed short disk
lifetimes in low metallicity clusters. A combination of our model with
previous studies (e.g., a metallicity-based photo-evaporation) could help to
understand disk evolution and dispersal at different metallicities.
Furthermore, the stellar spin evolution model includes several important
effects, previously ignored (e.g., the stellar magnetic field strength and a
realistic calculation of the disk lifetime) and we motivate to include our
results as initial or input parameters for further spin evolution models,
covering the stellar evolution towards and during the main sequence.

###### Key Words.:

protoplanetary disks – accretion, accretion disks – stars: protostars – stars:
rotation – stars: formation

##  1 Introduction

In the 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 on 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, or 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 ($L_{\mathrm{add}}$) 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
($L_{\mathrm{add}}\rightarrow 0$, cold accretion), the stellar radius is
significantly (up to one order of magnitude) smaller 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 in the stellar interior
($L_{\mathrm{add}}$ large, hot accretion), the stellar radius inflates and can
become comparable with classical evolutionary models. After accretion has
finished and the energy added into the stellar interior ceases, 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). Besides
the cold and warm accretion models, $L_{\mathrm{add}}$ can also depend on the
stellar accretion rate (hybrid accretion, e.g., Vorobyov et al. 2017a). In
periods of high accretion rates, $L_{\mathrm{add}}$ is large (similar to hot
accretion) and during low accretion rates $L_{\mathrm{add}}$ stays small
(similar to cold accretion).

During the early (Class 0 and Class I) phases of stellar evolution, the
accretion rate on the star can vary significantly over time. During episodic
outburst events (e.g., FU Orionis outbursts), the stellar accretion rate can
rise over orders of magnitude during a short timescale ($\sim$ years, e.g.,
Vorobyov & Basu 2010; Vorobyov et al. 2020). With $L_{\mathrm{add}}$ 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), affecting the stellar pre-MS evolution. In
previous studies, the accretion history on the star is either assumed to be
constant over time or pre-computed (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.
Up to 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, L. 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 to determine 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). In their work, Amard et al. (2019) present a
grid of low-mass ($\leq 1.5\leavevmode\nobreak\ \mathrm{M_{\odot}}$)
isochrones including 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 and 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. 2010a, 2016, 2021; Guarcello et al. 2021) and theoretical,
metallicity-dependent photo-evaporation models (e.g., Nakatani et al. 2018),
indicating 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, e.g., the disk lifetime, the stellar spin, and the pre-MS evolution,
including the back-reaction on the other component. In this first paper,
introducing the new method, we want to study the effects of metallicity on the
disk and stellar spin evolution of T Tauri stars (with ages $\gtrsim 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: Sec. 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 Sec. 3 and
discussed in Sec. 4. Finally, we draw our conclusions in Sec. 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).

$\displaystyle\frac{\partial}{\partial t}\,\Sigma$
$\displaystyle+\nabla\cdot(\Sigma\,\@vec{u})=0\;,$ (1)
$\displaystyle\frac{\partial}{\partial t}(\Sigma\,\@vec{u})$
$\displaystyle+\nabla\cdot(\Sigma\,\@vec{u}:\@vec{u})-\frac{B_{\mathrm{z}}\@vec{B}}{2\pi}$
$\displaystyle+\nabla P_{\mathrm{gas}}+\nabla\cdot
Q+\Sigma\,\nabla\psi+H_{\mathrm{P}}\nabla\left(\frac{B_{\mathrm{z}}^{2}}{4\pi}\right)=0\;,$
(2) $\displaystyle\frac{\partial}{\partial t}(\Sigma\,e)$
$\displaystyle+\nabla\cdot(\Sigma\,\@vec{u}\,e)+P_{\mathrm{gas}}\,\nabla\cdot\@vec{u}$
$\displaystyle+Q:\nabla\@vec{u}-4\pi\,\Sigma\,\kappa_{\mathrm{R}}\left(J-S\right)+\dot{E}_{\mathrm{rad}}=0\;,$
(3)

where $\Sigma$, $\@vec{u}$, $e$, $P_{\mathrm{gas}}$ and $H_{\mathrm{P}}$ stand
for the gas column density, gas velocity in the planar components
$\@vec{u}=(u_{\mathrm{r}},u_{\mathrm{\varphi}})$, 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 $\nabla=(\partial/\partial
r,r^{-1}\partial/\partial\varphi)$ with $\partial/\partial\varphi=0$ for our
axisymmetric model. $P_{\mathrm{gas}}$ is utilized by the ideal equation of
state $P_{\mathrm{gas}}=\Sigma e(\gamma-1)$, with the adiabatic coefficient
$\gamma=5/3$ and $\kappa_{\mathrm{R}}$ denotes the Rosseland-mean opacity,
which is composed of a gaseous component $\kappa_{\mathrm{R,gas}}$, which are
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
$\kappa_{\mathrm{R,dust}}$ (based on Pollack et al. 1985); with
$\kappa_{\mathrm{R}}=\kappa_{\mathrm{R,gas}}+f_{\mathrm{dust}}\kappa_{\mathrm{R,dust}}$,
where $f_{\mathrm{dust}}$ is the dust-to-gas mass ratio set to the respective
stellar metallicity. In our model, $f_{\mathrm{dust}}$ equals the stellar
metallicity, ranging from $Z=0.0028$ to $Z=0.014$. The gravitational potential
of the star-disk system is denoted by $\psi$. The vertical component of the
stellar magnetic field $B_{\mathrm{z}}$ is assumed to be a dipole and constant
within the disk’s vertical extent. The planar magnetic field components
$\@vec{B}=(B_{\mathrm{r}},B_{\mathrm{\varphi}})$ 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 $B_{\rm r}=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_{\rm\varphi}$ (e.g. Rappaport et al. 2004;
Kluźniak & Rappaport 2007; Steiner et al. 2021; Gehrig, L. 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 0th moment of the radiation
field and the source function, respectively (e.g. Ragossnig et al. 2020;
Steiner et al. 2021). Finally, $\dot{E}_{\mathrm{rad}}$ and $Q$ depict the net
heating/cooling rate per unit surface area and the viscous stress tensor,
respectively (e.g. Steiner et al. 2021; Gehrig, L. et al. 2022). The kinematic
viscosity is formulated according to Shakura & Sunyaev (1973) $\nu=\alpha
c_{\mathrm{S}}H_{\mathrm{P}}$; with the viscous $\alpha$ parameters, the
isothermal sound speed $c_{\mathrm{S}}$ and the pressure scale height
$H_{\mathrm{P}}$. 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 gridpoint resolution throughout the computational domain.

Computational domain: The innermost disk regions are assumed to be truncated
and eventually disrupted by a strong ($\sim$ kG) stellar magnetic field inside
the so-called truncation radius $r_{\rm trunc}$ (e.g. Romanova et al. 2002;
Bessolaz et al. 2008; Romanova & Kurosawa 2014). With our one-dimensional
code, the physical processes inside $r_{\rm trunc}$ can not be described.
Thus, we choose $r_{\rm trunc}$ to be the inner boundary of our computational
domain. The radial position of $r_{\rm trunc}$ 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, L. et
al. 2022)

$P_{\mathrm{magn}}(r_{\mathrm{trunc}})=\text{max}\left[P_{\mathrm{ram}}(r_{\mathrm{trunc}}),P_{\mathrm{gas}}(r_{\mathrm{trunc}})\right]\;.$
(4)

In case $P_{\mathrm{ram}}\geq P_{\mathrm{gas}}$, the truncation radius can be
calculated following Hartmann et al. (2016)

$\displaystyle r_{\mathrm{trunc}}(P_{\mathrm{ram}}\geq
P_{\mathrm{gas}})\approx 18\,\xi\,R_{\odot}$
$\displaystyle\,\left(\frac{B_{\star}}{10^{3}\,G}\right)^{4/7}\left(\frac{R_{\star}}{2\,R_{\odot}}\right)^{12/7}\left(\frac{M_{\star}}{0.5\,M_{\odot}}\right)^{-1/7}$
$\displaystyle\left(\frac{\dot{M}_{\star}}{10^{-8}\,M_{\odot}/\mathrm{yr}}\right)^{-2/7}\;,$
(5)

with the correction factor $\xi=0.7$, the stellar dipole magnetic field
strength $B_{\star}$, the stellar radius $R_{\star}$, the stellar mass
$M_{\star}$ and the accretion rate on the star $\dot{M}_{\star}$. In the case
of $P_{\mathrm{ram}}<P_{\mathrm{gas}}$, the position of the truncation radius
is calculated by equating the magnetic pressure of the vertical field
component and the gas pressure

$\displaystyle
r_{\mathrm{trunc}}(P_{\mathrm{ram}}<P_{\mathrm{gas}})=\frac{{B_{\star}}^{1/3}\,R_{\star}}{{P_{\mathrm{gas}}}^{1/6}\,(8\pi)^{1/6}}\;.$
(6)

$r_{\rm trunc}$ 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
$\Sigma_{\rm out}=1.0\leavevmode\nobreak\ {\rm g/cm^{2}}$ (similar to, e.g.,
Vorobyov et al. 2020).

Boundary conditions: At the inner boundary, we choose ”free”/Neumann boundary
conditions (e.g. $\partial\Sigma/\partial r=0$, Romanova et al. 2002, 2004;
Steiner et al. 2021). At the outer boundary, we fix the radial (angular
velocity) $u_{\rm r}$ ($u_{\rm\varphi}$) to zero (the Keplerian value),
respectively. 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 $\Sigma$ 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 Modules for Experiments in Stellar Astrophysics (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
$X_{{}^{1}\mathrm{H}}=1-X_{{}^{2}\mathrm{H}}-X_{{}^{3}\mathrm{He}}-X_{{}^{4}\mathrm{He}}-Z$,
where $X_{{}^{2}\mathrm{H}}=20$ ppm, $X_{{}^{3}\mathrm{He}}=85$ ppm,
$X_{{}^{4}\mathrm{He}}=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 follow 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

$\frac{\mathrm{d}E_{\rm
acc}}{\mathrm{d}t}=(\epsilon-1)\frac{GM_{\star}\dot{M}_{\star}}{R_{\star}}\,,$
(7)

which is a combination of the gravitational ($-GM_{\star}/R_{\star}$) and
internal energy ($+\epsilon GM_{\star}/R_{\star}$) of the accreted material.
The geometry of the accretion process is described by the free parameter
$\epsilon$. We follow the approach of previous works (Baraffe et al. 2009;
Jensen & Haugbølle 2018; Steindl et al. 2021) and set $\epsilon=0.5$,
corresponding to accretion from a thin disc 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 $\epsilon
GM_{\star}\dot{M_{\star}}/R_{\star}$.

The energy of the accreted material is either added to the star as extra heat
($L_{\rm add}$) or radiated away as accretion luminosity ($L_{\rm acc}$). 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

$L_{\rm add}=\beta\epsilon\frac{GM_{\star}\dot{M}_{\star}}{R_{\star}}$ (8)

and

$L_{\rm acc}=(1-\beta)\epsilon\frac{GM_{\star}\dot{M}_{\star}}{R_{\star}}.$
(9)

The value $\beta$ and its dependence on star and disk properties are free
parameters. From earlier studies it is evident that $\beta$ 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 dependency first
presented by Jensen & Haugbølle (2018) given as

$\beta(\dot{M_{\star}})=\frac{\beta_{1}\exp\left(\frac{\dot{M}_{\rm
m}}{\Delta}\right)+\beta_{2}\exp\left(\frac{\dot{M}_{\star}}{\Delta}\right)}{\exp\left(\frac{\dot{M}_{\rm
m}}{\Delta}\right)+\exp\left(\frac{\dot{M}_{\star}}{\Delta}\right)}$ (10)

This describes a step function with a smooth transition between
$\beta_{\mathrm{1}}$ and $\beta_{\mathrm{2}}$. The functional dependency is
controlled by four parameters: $\beta_{1}=0.005$, $\beta_{2}=0.2$, a midpoint
accretion rate of $\dot{M}_{\rm m}=6.2\cdot 10^{-6}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ and a crossover width of
$\Delta=5.95\cdot 10^{-6}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ 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 $\beta$ at the beginning of the evolution which
becomes gradually smaller as the accretion rate decreases (see Figure 2).
During the calculation of the initial models, the accretion rates and thus
$\beta(\dot{M_{\star}})$ vary significantly. Given the proportionately small
accretion rates, $\beta$ is basically constant during the joint evolution of
MESA and TAPIR.

The dependence of $\beta$ 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 towards 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 $\beta$ values is expected to occur at $\sim
10^{-6}\leavevmode\nobreak\ -\leavevmode\nobreak\ 10^{-5}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ (e.g, Vorobyov et al. 2017a;
Jensen & Haugbølle 2018), which motivates the choice of $\dot{M}_{\rm m}$. We
note, however, that the description of $\beta$ 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 towards 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
$\Delta$ should include a metallicity-dependent factor, that would exceed the
frame of this study.

The heat added to the stellar model is distributed according to

$l_{\rm extra}=\frac{L_{\rm add}}{M_{\star}}{\rm max}\left\\{0,\frac{2}{M_{\rm
outer}^{2}}\left(\frac{m_{r}}{M_{\star}}-1+M_{\rm outer}\right)\right\\}.$
(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 $M_{\rm outer}$
and increases linearly with the mass coordinate $m_{r}$. In our study, we
choose $M_{\rm outer}=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, L. et al. (2022).
In this study, we assume solid body rotation for the star. Angular momentum
transport within the star is discussed in Sec. 4.3. The stellar angular
momentum $J_{\star}=I_{\star}\Omega_{\star}$ is influenced by external torques
$\Gamma_{\rm ext}$ and its temporal derivative reads

$\dot{J}_{\star}=I_{\star}\dot{\Omega}_{\star}+\dot{I}_{\star}\Omega_{\star}=\Gamma_{\mathrm{ext}}\,,$
(12)

with the stellar angular momentum $\Omega_{\star}$ and moment of inertia
$I_{\star}$. Rearranging Eq. (12) results in a time-dependent equation for the
stellar spin

$\dot{\Omega}_{\star}=\frac{\Gamma_{\mathrm{ext}}}{I_{\star}}-\frac{\dot{I}_{\star}}{I_{\star}}\Omega_{\star}\,.$
(13)

The stellar moment of inertia $I_{\star}$ is calculated by the MESA code.
Within $\Gamma_{\rm ext}$ the effects of accretion $\Gamma_{\mathrm{acc}}$,
stellar winds $\Gamma_{\mathrm{W}}$ (in form of an accretion powered stellar
wind APSW, Matt & Pudritz 2005) and the magnetic star-disk connection
$\Gamma_{\mathrm{ME}}$ (in 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 that increases $I_{\star}$ and spins
up the star, respectively. An APSW, on the other hand, is assumed to eject a
specific amount of the accreted material ($\dot{M}_{\rm W}=W\dot{M}$) removing
stellar mass and angular momentum. $W$ is expected to be $\lesssim 2$ % (e.g.
Cranmer 2008; Pantolmos et al. 2020) and thus, we choose $W=2$ % in this
study. The effect of MEs on the stellar spin evolution depends on the position
of $r_{\rm trunc}$ with respect to the co-rotation radius

$r_{\mathrm{cor}}=\left(\frac{G\,M_{\star}}{\Omega_{\star}^{2}}\right)^{1/3}\;,$
(14)

with the gravitational constant $G$ and the stellar mass $M$. MEs increase the
stellar angular momentum and spin up the star if
$r_{\mathrm{trunc}}<K_{\mathrm{rot}}^{2/3}r_{\mathrm{cor}}$. Otherwise, the
star spins down. Following Gallet et al. (2019), we use $K_{\mathrm{rot}}=0.7$
in this study111We note that in more recent studies of MEs (e.g. Pantolmos et
al. 2020; Ireland et al. 2021) smaller values for $K_{\rm rot}$ are found.
$K_{\rm rot}=0.7$ can be understood as an upper limit (Ireland et al. 2021).
and the external torque contributions scale as:

$\Gamma_{\mathrm{acc}}=K_{\mathrm{acc}}\dot{M}_{\mathrm{\star}}r_{\mathrm{trunc}}^{2}\Omega_{\mathrm{disk,in}}\,$
(15)
$\Gamma_{\mathrm{W}}=\dot{M}_{\mathrm{W}}r_{\mathrm{A}}^{2}\Omega_{\star}\,\,\mathrm{and}$
(16)
$\Gamma_{\mathrm{ME}}=K_{\mathrm{ME}}\frac{B_{\star}^{2}R_{\star}^{6}}{r_{\mathrm{trunc}}^{3}}\left[K_{\mathrm{rot}}-\left(\frac{r_{\mathrm{trunc}}}{r_{\mathrm{cor}}}\right)^{3/2}\right]\,,$
(17)

with $K_{\mathrm{acc}}=0.4$, $K_{\mathrm{ME}}=0.21$, the disk rotation rate at
the truncation radius $\Omega_{\mathrm{disk,in}}$, and the stellar Alfvén
radius

$r_{\mathrm{A}}=K_{\mathrm{1}}\left[\frac{B_{\star}^{2}R_{\star}^{2}}{\dot{M}_{\mathrm{wind}}\sqrt{K_{\mathrm{2}}^{2}v_{\mathrm{esc}}^{2}+\Omega_{\star}^{2}R_{\star}^{2}}}\right]^{m}R_{\star}\,,$
(18)

where $v_{\mathrm{esc}}=\sqrt{2\mathrm{G}M_{\star}/R_{\star}}$, $m=0.2177$,
$K_{\mathrm{2}}=0.0506$, and $K_{\mathrm{1}}\leavevmode\nobreak\
=\leavevmode\nobreak\ 1.7$.

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

The interaction between the disk (TAPIR, Sec. 2.1), spin (stellar spin model,
Sec. 2.3), and stellar model (MESA, Sec. 2.2) is schematically shown in Fig.
1. Starting from the initial or old values at timestep $t_{\mathrm{0}}$, the
new disk values, spin values, and stellar values at timestep
$t_{\mathrm{1}}=t_{\mathrm{0}}+\delta t$ are calculated in 3 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
$t_{\mathrm{0}}$ to $t_{\mathrm{1}}$ 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 $t_{\mathrm{0}}$ to $t_{\mathrm{1}}$ (Arrow 3).

The accuracy and reliability of our method depend on the timestep size $\delta
t$. Choosing small timesteps results in accurate but long computational times.
Large timesteps will reduce computational time but accuracy suffers. We choose
$\delta t=100$ yr in this study. During 1 Myr the stellar values are updated
$10^{4}$ times resulting in an acceptable representation of the stellar
evolution. We also tested a range of $\delta 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
$t_{\mathrm{0}}$ to $t_{\mathrm{1}}$. For MESA, the internal timesteps range
between $\sim 10^{6}$ s and $\sim 10^{8}$ s. For TAPIR, the timesteps range
between $\sim 10^{4}$ s and $\sim 10^{8}$ 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 time scales and influence the stellar
evolution (e.g. Vorobyov et al. 2017a). In further studies, such outbursts are
included in our model and $\delta t$ has to be adapted to lower values to
assure a sufficiently accurate stellar evolution calculation.

Figure 1: Schematic representation of the interaction between TAPIR and MESA.
Starting from an initial/old timestep $t_{0}$, the disk, spin and star values
are evolved towards the new timestep $t_{1}=t_{0}+\delta t$. The new values at
$t_{1}$ are obtained in three steps, represented by the numbered arrows (see
text).

##  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 choose $t_{\rm 0}=2$ Myr as an
initial time for our combined simulations. The effects of the early evolution
and different values for $t_{\rm 0}$ are discussed in Sec. 4.2. For each
model, the stellar metallicity is varied between $Z=0.2\leavevmode\nobreak\
\mathrm{Z_{\odot}}$ (motivated by observations of young, low metallicity
clusters Yasui et al. 2010a, 2016, 2021) and $Z=1\leavevmode\nobreak\
\mathrm{Z_{\odot}}$.

Disk parameters: One important parameter that defines an accretion disk is the
mass accretion rate. For T Tauri star-disk systems ($\lesssim 10$ Myr)
accretion rates range within $10^{-10}$ and $10^{-7}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ (e.g. Gullbring et al. 1998;
Manara et al. 2016; Vorobyov et al. 2017a; Manara et al. 2022). For $t_{\rm
0}=2$ Myr, Manara et al. (2012) and Testi et al. (2022) find disk accretion
rates for a $\sim 1\leavevmode\nobreak\ \mathrm{M_{\odot}}$ star between
$10^{-8}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$
and $10^{-9}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\
yr^{-1}}$. Furthermore, the viscous $\alpha$-parameter can have a wide range
between $0.1$ and $\sim 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 choose the initial (starting) accretion rate
$\dot{M}_{\rm start}=[1\times 10^{-9},5\times 10^{-9},1\times
10^{-8}]\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$
and set the (constant) viscous $\alpha$-value to 0.01 unless otherwise stated.
This combination results in the expected disk lifetimes of $\lesssim 10$ Myr
(consistent with e.g. Richert et al. 2018). We stop our simulations when the
disk accretion rate has decreased to $10^{-11}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$. This lower limit is motivated
by the lowest detectable accretion rate for a $1\leavevmode\nobreak\
\mathrm{M_{\odot}}$ 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 $t_{\mathrm{0}}$ result from
the initial stellar evolution ($t<t_{\mathrm{0}}$). 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\leavevmode\nobreak\ \mathrm{M_{\odot}}$
($\approx$ 10 Jupiter masses) and a radius of $1.5\leavevmode\nobreak\
\mathrm{R_{\odot}}$ (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$ $\mathrm{M_{\odot}}$ and the desired accretion rate
$\dot{M}_{\mathrm{start}}$ at time $t_{\rm 0}=2$ Myr. The effect of the
accretion history is described in Sec. 3.1. During the T Tauri phase, the
rotational period can span over an order of magnitude from $\lesssim 1$
towards $\sim 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 $P_{\mathrm{start}}=10$ days (comparable to the slow
models in Amard et al. 2019). During the initial stellar evolution
($t<t_{\mathrm{0}}$), 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 combining stellar and disk
evolution, initial models, which represent a T Tauri star and its disk at an
age of $t_{\mathrm{0}}$ have to be calculated. The initial models are based on
the aforementioned stellar and disk parameters.

###  3.1 Initial stellar evolution: $t<t_{\rm 0}$

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 $\dot{M}_{\rm init}$ followed by a decrease towards our desired
value of $\dot{M}_{\rm start}$ according to $\dot{M}_{\star}\propto t^{-2}$
(see Fig. 2). We note that for the models with $\dot{M}_{\rm start}=1\times
10^{-9}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$,
we use $\dot{M}_{\star}\propto t^{-2.3}$ to avoid values for $\dot{M}_{\rm
init}>2\times 10^{-5}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$. The reason for this variation
is convergence problems in the MESA initial model with higher initial
accretion rates (e.g., Sec. 5.1.5 in Steindl et al. 2021). The values of
$\dot{M}_{\rm init}$ and the time, at which the accretion rate begins to
decrease, $t_{\rm dec}$, depend on $\dot{M}_{\rm start}$ as well as the
desired stellar mass at $t_{\rm 0}$. For a starting accretion rate of
$\dot{M}_{\rm start}=[1\times 10^{-9},5\times 10^{-9},1\times
10^{-8}]\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$,
we find $\dot{M}_{\rm init}=[1.87\times 10^{-5},1.03\times 10^{-5},5.29\times
10^{-6}]\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$
and $t_{\rm dec}=[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 towards the Hayashi track corresponding to an evolution
time $\tau_{\mathrm{H}}$, from $\tau_{\mathrm{H}}$ along the Hayashi track
towards $t_{\mathrm{0}}$ and finally reaching $t_{\mathrm{0}}$, which marks
the starting point for our combined stellar and disk model.

Figure 2: Initial stellar evolution for $t<t_{\rm 0}$. Panel (a): Accretion
history of the star. Initially, the star accretes with a constant accretion
rate $\dot{M}_{\rm init}$. At time $t_{\rm dec}$, the accretion rate drops
according to $\dot{M}_{\star}\propto t^{-2}$ (see text). The colored dashed
lines show the respective value of $\dot{M}_{\rm start}$. Panel (b): The
evolution of the stellar mass. Starting from the initial mass of
$0.01\leavevmode\nobreak\ \mathrm{M_{\odot}}$, all models reach their final
mass of $0.91\leavevmode\nobreak\ \mathrm{M_{\odot}}$ after 2 Myr. The yellow,
red and blue lines show models with $\dot{M}_{\rm start}=[1\times
10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ 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\leq 10^{4}$ yr) to the logarithmic
scale (for $t>10^{4}$ yr).

From the stellar seed to $\tau_{\mathrm{H}}$: Considering the evolution of the
star on a Hertzsprung-Russell-Diagram (HRD) for $t<t_{\rm 0}$, 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 Sec. 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 $\dot{M}_{\rm init}$ results in a
higher amount of energy added to the star.222The different behavior at the
beginning of the evolution (the first $\sim 100$ years) can be explained by
the outer boundary condition of the MESA code that reacts differently on
different values of $L_{\mathrm{add}}$, which depends on the accretion rate.
We note that a higher $\dot{M}_{\rm init}$ results in a lower $\dot{M}_{\rm
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 kyr, Deuterium
burning sets in ($\tau_{\mathrm{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
($\tau_{\mathrm{H}}$, marked with crosses in Fig. 3) indicates a transition.
Before $\tau_{\mathrm{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 towards
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 towards the Hayashi track is marked with a cross (’x’). The
respective stellar mass ($M_{\mathrm{H}}$) and evolution time
($\tau_{\mathrm{H}}$) is summarized in Table 1. The stellar masses
($M_{\mathrm{H}}$) range between 0.66 and 0.83 $\mathrm{M_{\odot}}$ (with
slightly lower values for $Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$). The
time $\tau_{\mathrm{H}}$ depends, as expected, on the accretion history.
Shortly after the accretion rate decreases from $\dot{M}_{\mathrm{init}}$
towards $\dot{M}_{\mathrm{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 $t_{\mathrm{0}}=2$ Myr: When reaching an age of $t_{\mathrm{0}}=2$
Myr, stellar radii and temperatures (see Fig. 4) are used to calculate the
initial disk model (see Sec. 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\leavevmode\nobreak\
\mathrm{Z_{\odot}}$, they are located on a vertical line indicating evolution
on the Hayashi track (stellar luminosity decreasing with age). For
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ however, the stars at
$t_{\mathrm{0}}$ already turn towards the Henyey track (stellar luminosity and
effective temperature increase with age) indicating the development of a
radiative core.

Table 1: Onset of Deuterium burning ($\tau_{\mathrm{D}}$), stellar mass ($M_{\mathrm{H}}$) and evolution time ($\tau_{\mathrm{H}}$) the star turns onto the Hayashi track $\dot{M}_{\mathrm{start}}$ [$\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$] | $Z\leavevmode\nobreak\ \mathrm{[Z_{\odot}]}$ | $M_{\mathrm{H}}\leavevmode\nobreak\ \mathrm{[M_{\odot}]}$ | $\tau_{\mathrm{H}}\leavevmode\nobreak\ \mathrm{[Myr]}$ |  $\tau_{\mathrm{D}}$ [kyr]
---|---|---|---|---
$1\times 10^{-9}$ | 1 | 0.82 | 0.09 | 1.45
$5\times 10^{-9}$ | 1 | 0.71 | 0.10 | 2.31
$1\times 10^{-8}$ | 1 | 0.69 | 0.16 | 2.72
$1\times 10^{-9}$ | 0.2 | 0.81 | 0.09 | 1.31
$5\times 10^{-9}$ | 0.2 | 0.66 | 0.08 | 1.84
$1\times 10^{-8}$ | 0.2 | 0.66 | 0.15 | 2.59

Figure 3:  Evolution on the HRD for $t<t_{\rm 0}$ for both metallicities
$Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines) and
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (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 $\dot{M}_{\rm start}=[1\times
10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, respectively. The onset of
Deuterium burning ($t=\tau_{\mathrm{D}}$) is marked with a diamond. The
transition from expansion towards the Hayashi track ($t=\tau_{\mathrm{H}}$) is
marked with a cross.

Figure 4:  Stellar radii and temperatures at $t_{\mathrm{0}}=2$ Myr for
$\dot{M}_{\mathrm{start}}=5\times 10^{-9}$ (red triangles) and $1\times
10^{-8}$ (blue triangles). Solar (sub-solar) metallicity values are connected
by a solid (dashed) line, respectively.

###  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 construct 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 $\dot{M}_{\mathrm{start}}$. In Fig. 5,
the radial surface density (Panel a), disk midplane temperature (Panel b), and
the optical depth $\tau$ (Panel c) profile of the initial disk models are
shown. The stellar metallicities are represented by solid lines
($Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$) and dashed lines
($Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$) respectively. The yellow, red
and blue lines symbolize models with $\dot{M}_{\rm start}=[1\times
10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, respectively.

For larger initial accretion rates, the disks are more massive, have a higher
temperature, and are pushed closer towards 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 $t_{\mathrm{0}}$ (see Fig. 4), the
absolute values of the inner disk edge $r_{\mathrm{trunc}}$ is smaller for
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$. The radial temperature profile
can be divided into three regions. Close to the star the disk is directly
heated by stellar irradiation $r\lesssim 10^{-1}$. For low metallicities,
stellar luminosity is larger and, thus, the disk temperature is higher. From
$r\lesssim 10^{-1}$ AU to $\sim 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 optical thick zone, $\mathrm{R_{thick}}$,
from now on. In this region, the disk midplane temperature ranges from $\sim
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, $f_{\mathrm{dust}}=0.014$ and $0.0028$ for
$Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ and $0.2\leavevmode\nobreak\
\mathrm{Z_{\odot}}$, 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
$\dot{M}_{\mathrm{start}}=1\times 10^{-9}\mathrm{M_{\odot}\leavevmode\nobreak\
yr^{-1}}$, these effects are almost negligible due to the low surface density
values. Outside $R_{\mathrm{thick}}$, 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 $T_{\mathrm{amb}}=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 $R_{\mathrm{c}}$, and
masses in Table 2. The characteristic radius of the disk $R_{\mathrm{c}}$
includes $2/3$ of the disk’s mass. The values of $R_{\mathrm{c}}$, ranging
between 20 AU and 128 AU, and $m_{\mathrm{disk}}$, ranging from 0.001
$\mathrm{M_{\odot}}$ to 0.055 $\mathrm{M_{\odot}}$ 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 $R_{\mathrm{c}}$ and disk mass $m_{\mathrm{disk}}$ of the initial disk model  $\dot{M}_{\mathrm{start}}$ [$\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$] | $Z\leavevmode\nobreak\ \mathrm{[Z_{\odot}]}$ |  $R_{\mathrm{c}}$ [AU] |  $m_{\mathrm{disk}}$ [$\mathrm{M_{\odot}}$]
---|---|---|---
$1\times 10^{-9}$ | 1 | 23 | 0.003
$1\times 10^{-9}$ | 0.2 | 20 | 0.001
$5\times 10^{-9}$ | 1 | 78 | 0.019
$5\times 10^{-9}$ | 0.2 | 72 | 0.017
$1\times 10^{-8}$ | 1 | 128 | 0.055
$1\times 10^{-8}$ | 0.2 | 121 | 0.051

Figure 5:  The radial profile of the initial disk model for both metallicities
$Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines) and
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (dashed lines). The yellow, red
and blue lines indicate models with $\dot{M}_{\rm start}=[1\times
10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, respectively. Panel (a) shows
the disk surface density $\Sigma$, Panel (b) the disk midplane temperature
$T_{\mathrm{disk}}$, and Panel (c) the optical depth $\tau$, respectively. The
thin vertical lines in Panel (a) indicate the respective position of
$R_{\mathrm{c}}$. The horizontal line in Panel (b) indicates the ambient
temperature $T_{\mathrm{amb}}=20$ K. The horizontal dashed line in Panel (c)
symbolizes $\tau=1$, dividing optical thin regions ($\tau<<1$) from optical
thick regions ($\tau>>1$).

###  3.3 Disk evolution: $t>t_{\rm 0}$

Starting from the stellar and disk initial models, we evolve the star-disk
system until the accretion rate onto the star decreases below
$10^{-11}\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$. In Fig. 6, the
evolution of the disk accretion rate is shown for both metallicities
$Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines) and
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (dashed lines). The red and
blue lines indicate models with $\dot{M}_{\rm start}=[1\times 10^{-9},5\times
10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, respectively. Independent of
the accretion history of the star, a lower stellar metallicity results in a
shorter disk lifetime333In this context, disk lifetime is the time, at which
the accretion rate has dropped to
$10^{-11}\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$..

Observations of disk fractions in low metallicity clusters (Yasui et al.
2010a) 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 $R_{\mathrm{thick}}$) and the outer disk regions (outside
$R_{\mathrm{thick}}$) are hotter compared to higher metallicities (see Fig.
5). With higher disk temperatures, the disk’s viscosity $\nu$ increases. Thus,
the viscous timescale $\tau_{\mathrm{\nu}}=r^{2}/\nu$ 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. 11, in which the evolution of the
stellar parameters for $t\geq t_{\mathrm{0}}$ are summarized, the differences
in intrinsic luminosity between solar and sub-solar 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 sub-solar
metallicities. Usually, the (global) viscous timescale is measured in the
outer ($\sim 10$ AU) disk regions (e.g., Armitage 2011). During the past
years, however, recent studies highlighted the importance of the inner disk,
the sub-AU region (e.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 shield a certain part of the disk from 25 % of the stellar
irradiation, resulting in two test models with a lower disk temperature in the
inner disk region $<1$ AU (model $\mathrm{L_{inlow}}$) and the outer disk
region $>1$ AU ($\mathrm{L_{outlow}}$), respectively (see Panel b of Fig. 7).
We compare the disk evolution of the models $\mathrm{L_{inlow}}$ and
$\mathrm{L_{outlow}}$ 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.

Figure 6:  Accretion rate onto the star over time for both metallicities
$Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines) and
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (dashed lines). The yellow, red
and blue lines indicate models with $\dot{M}_{\rm start}=[1\times
10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, respectively. The black
dotted line shows the termination criterion used in our simulations of
$\dot{M}=10^{-11}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\
yr^{-1}}$.

Figure 7:  Comparison between models $\mathrm{L_{inlow}}$ (yellow line) and
$\mathrm{L_{outlow}}$ (red line). The reference model (ref, blue line)
corresponds to a $Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ star with
$\dot{M}_{\mathrm{start}}=1\times 10^{-8}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$. The stellar luminosity is
artificially reduced by 25 % in the inner 1 AU ($\mathrm{L_{inlow}}$) and
outside 1 AU ($\mathrm{L_{outlow}}$). Panel (a) shows the evolution of the
accretion rate. The black dotted line shows the termination criterion used in
our simulations of $\dot{M}=10^{-11}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$. Panel (b) shows the
temperature structure of the initial disk model. For model
$\mathrm{L_{inlow}}$, the inner disk region is cooler and for model
$\mathrm{L_{outlow}}$, the outer disk region is cooler compared to the
reference model. The black dotted 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
($\gtrsim 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 Sec. 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\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines)
and $Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (dashed lines). The yellow,
red and blue lines indicate models with $\dot{M}_{\rm start}=[1\times
10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, respectively. All simulations
are started at a stellar period at $t_{\mathrm{0}}$ of $P_{\mathrm{0}}=10$
days. We expect that a variation in $P_{\mathrm{0}}$ does not change the
qualitative result shown in Fig. 8, but rather shifts 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.

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 $t_{\mathrm{0}}$. According to the stellar spin model described in
Sec. 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. 11 shows the evolution
of the stellar radii over time for both metallicities $Z=1\leavevmode\nobreak\
\mathrm{Z_{\odot}}$ (solid lines) and $Z=0.2\leavevmode\nobreak\
\mathrm{Z_{\odot}}$ (dashed lines). 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, $\sim 1$ Myr shorter (see Sec. 3.3).
After the disk is dissolved, there are no mechanisms that can remove angular
momentum and the star spins up due to contraction. 444We 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 $\sim 1$ Myr. After contraction has ended and
the star evolves towards and on the MS, stellar winds slow down the star on
longer timescales (e.g., Gallet & Bouvier 2013).  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.

Figure 8:  Stellar rotational periods over time for both metallicities
$Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines) and
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (dashed lines). The yellow, red
and blue lines indicate models with $\dot{M}_{\rm start}=[1\times
10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ 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.

###  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_{\star}$
and the disk’s viscous $\alpha$-parameter. We want to verify the robustness of
our results by changing these parameters. For our models with $\dot{M}_{\rm
start}=[1\times 10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ 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 $\alpha$-parameter is reduced to 0.005 (model
MC3), respectively.

Fig. 9 shows the evolution of the accretion disk (Panel a, see Fig. 6) and the
stellar spin evolution (Panel b, see Fig. 8) for models MC1, MC2, and MC3.
Similar to the results shown in Sec. 3.3 and Sec. 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.

Figure 9:  Influence of the stellar magnetic field strength $B_{\star}$,
initial rotation period, and viscous $\alpha$-parameter on the results for
both metallicities $Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines)
and $Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (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 without the post-disk
contraction (pdc).

##  4 Discussion

###  4.1 Comparing our results to observed rotation periods

At ages, $\gtrsim 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 up
to 2 days, see Sec. 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
co-exist expanding the age span within one cluster up to several Myr. With
these limitations in mind, we choose 6 young clusters with ages between $\sim
2$ Myr and $\lesssim 5$ Myr. The clusters used in this comparison are NGC
6530, NGC 2362, NGC 2264, Orion star-forming cluster (OSFC), $\sigma$ Ori and
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 3. We note that our models are not included in the mass range of
observed periods in $\sigma$ 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 choose to include $\sigma$ Ori in this
comparison.

The observed rotation periods are combined in two bins: the cluster
metallicity is assumed to be solar or above ($\mathrm{P_{Zhigh}}$) and the
metallicity is assumed to be below the solar value ($\mathrm{P_{Zlow}}$). In
Fig. 10, the ranges of observed rotation periods from the $25^{\mathrm{th}}$
to the $90^{\mathrm{th}}$ percentile are shown for $\mathrm{P_{Zhigh}}$ (light
grey) and $\mathrm{P_{Zlow}}$ (light blue). The mean values
($\mathrm{\bar{P}}$) 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, e.g., the stellar magnetic field strength. We note that the lower
stellar mass range of observed periods in $\sigma$ 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.

Figure 10:  Stellar rotational periods over time for both metallicities
$Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (solid lines) and
$Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (dashed lines). The colored
lines represent the respective model. The ranges of the observed periods of
young clusters with ages between $\sim 3$ Myr and $\lesssim 5$ Myr (see text)
are indicated by the light grey (light blue) area for $\mathrm{P_{Zhigh}}$
($\mathrm{P_{Zlow}}$), respectively. The mean values ($\mathrm{\bar{P}}$) are
indicated by horizontal lines in the respective color.

###  4.2 Early evolution of the star-disk system: Episodic accretion during
$t<t_{\mathrm{0}}$

In this study, we start the combined star-disk simulation at
$t_{\mathrm{0}}=2$ Myr. For ages $t<t_{\mathrm{0}}$ fixed parameters are used
to define the accretion history on the star (see Fig. 2). The accretion rate
for $t<t_{\mathrm{0}}$, 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 Sec. 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
$t_{\mathrm{0}}$ 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
towards 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), stars rotate either close to solid-body rotation (valid for
fast rotating stars during the first $\lesssim 10^{8}$ years) 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, main
sequence 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 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_{\star}$

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_{\star}$ 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_{\star}$ 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
towards 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 of
$B_{\star}$ for low metallicity stars could be visualized by a gradual
transition from the blue-dashed line (where $B_{\star}=2.0$ kG is assumed)
towards the cyan-dashed line ($B_{\star}=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. 2010a, 2016;
Guarcello et al. 2021; Yasui 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\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ cluster is $\sim 10$ % (e.g.
Yasui et al. 2010a) 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
deeper penetration of high energy radiation in 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 sub-solar 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) proposes 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 higher metallicities.

###  4.6 Differences between our results and observations

The observations shown in Yasui et al. (2010b) suggest disk lifetimes of $\sim
3$ Myr for low-metallicity clusters, matching our results for $\dot{M}_{\rm
start}=1\times 10^{-9}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$. In this context, cluster ages
can be used to estimate disk lifetimes. The disk lifetimes for models with
$\dot{M}_{\rm start}=5\times 10^{-9}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ and $1\times
10^{-8}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$,
on the other hand, are longer compared to observations by a factor of $\sim
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 towards the
observed values. Second, the cluster ages presented by Yasui et al. (2010b)
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 $\sim 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
$t_{\mathrm{0}}$ 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 including 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. 2010a, 2016, 2021;
Guarcello et al. 2021), photo-evaporation models (e.g., Nakatani et al. 2018)
and the results presented in Sec. 3.3.

Starting from the T-Tauri phase, our star-disk evolution model does include
the effects of, e.g., 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) suggests 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, R. et al. (2022) showed that snowlines 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 $\gtrsim 10$ AU,
where the disk temperature is dominated by stellar irradiation. Within
$R_{\mathrm{thick}}$, however, the disk midplane temperature decreases with
decreasing metallicity (see Sec. 3.2) and snowlines of species with
sublimation temperatures $\gtrsim 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 $t_{\mathrm{0}}$ has
to be shifted towards 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 $\alpha$) viscosity prescription.
$\alpha$ 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 $f_{\mathrm{dust}}$ 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.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 $\dot{M}_{\mathrm{init}}$, the more energy is added to
the stellar interior during the early stellar evolution $t<t_{\mathrm{0}}$ and
the larger the stellar radius. On the HRD, the star appears less evolved
compared to small values of $\dot{M}_{\mathrm{init}}$ (see Fig. 2 and Fig. 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 $\lesssim 0.1$ AU as well as the disk beyond $\gtrsim 1$ AU are
dominated by stellar irradiation and thus hotter compared to solar metallicity
disks. Within $R_{\mathrm{thick}}$, however, the disk midplane is shielded
from stellar irradiation and the temperature is dominated by dust opacities.
High stellar metallicities results 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. 2010a). 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 towards the
understanding of 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 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 Sec. 2.3). Additionally, the stellar contraction rate depends
on the stellar evolutionary stage. The stellar contraction rate is smaller if
a star is further evolved. As a consequence, low metallicity stars contract
more slowly (see Fig. 11) 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 sub-solar
metallicity stars increases further.

  * •

Including 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, we motivate that our results 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 the joint hydrodynamic evolution
of the star and the disk. We motivate the step-by-step inclusion of additional
physical aspects (e.g., photo-evaporation and magneto-centrifugal disk winds)
and studies covering wider parameter spaces to improve our understanding of
these young star-disk systems.

###### Acknowledgements.

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

  * Affer et al. (2013) Affer, L., Micela, G., Favata, F., Flaccomio, E., & Bouvier, J. 2013, MNRAS, 430, 1433
  * Amard & Matt (2020) Amard, L. & Matt, S. P. 2020, ApJ, 889, 108
  * Amard et al. (2016) Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, A&A, 587, A105
  * Amard et al. (2019) Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77
  * Amard et al. (2020) Amard, L., Roquette, J., & Matt, S. P. 2020, MNRAS, 499, 3481
  * Andrews et al. (2010) Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241
  * Angulo et al. (1999) Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A, 656, 3
  * Ansdell et al. (2018) Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21
  * Armitage (2011) Armitage, P. J. 2011, Annual Review of Astronomy and Astrophysics, 49, 195
  * Armitage et al. (2001) Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705
  * Bae et al. (2014) Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61
  * Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, Astrophys. J., 376, 214
  * Baraffe et al. (2009) Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJ, 702, L27
  * Baraffe et al. (2012) Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118
  * Barnes (2007) Barnes, S. A. 2007, ApJ, 669, 1167
  * Bastien et al. (2011) Bastien, F. A., Stassun, K. G., & Weintraub, D. A. 2011, AJ, 142, 141
  * Benomar et al. (2015) Benomar, O., Takata, M., Shibahashi, H., Ceillier, T., & García, R. A. 2015, Monthly Notices of the Royal Astronomical Society, 452, 2654
  * Berg et al. (2016) Berg, D. A., Skillman, E. D., Henry, R. B. C., Erb, D. K., & Carigi, L. 2016, The Astrophysical Journal, 827, 126
  * Bessolaz et al. (2008) Bessolaz, N., Zanni, C., Ferreira, J., Keppens, R., & Bouvier, J. 2008, A&A, 478, 155
  * Bhandare et al. (2018) Bhandare, A., Kuiper, R., Henning, T., et al. 2018, A&A, 618, A95
  * Buchler & Yueh (1976) Buchler, J. R. & Yueh, W. R. 1976, ApJ, 210, 440
  * Buldgen & Eggenberger (2022) Buldgen, G. & Eggenberger, P. 2022, arXiv e-prints, arXiv:2202.10086
  * Burssens et al. (2021) Burssens, S., Bowman, D. M., Michielsen, M., Simón-Díaz, S., & Aerts, C. 2021, Internal rotation and mixing in the massive star HD192575
  * Caffau et al. (2011) Caffau, E., Ludwig, H. G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094
  * Cauley et al. (2012) Cauley, P. W., Johns-Krull, C. M., Hamilton, C. M., & Lockhart, K. 2012, The Astrophysical Journal, 756, 68
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028
  * Cody & Hillenbrand (2010) Cody, A. M. & Hillenbrand, L. A. 2010, The Astrophysical Journal Supplement Series, 191, 389
  * Cranmer (2008) Cranmer, S. R. 2008, ApJ, 689, 316
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240
  * Damiani et al. (2019) Damiani, F., Prisinzano, L., Micela, G., & Sciortino, S. 2019, A&A, 623, A25
  * D’Orazi et al. (2011) D’Orazi, V., Biazzo, K., & Randich, S. 2011, A&A, 526, A103
  * D’Orazi et al. (2009) D’Orazi, V., Randich, S., Flaccomio, E., et al. 2009, A&A, 501, 973
  * Dorfi & Drury (1987) Dorfi, E. & Drury, L. 1987, J. Comput. Phys., 69, 175
  * Dorfi et al. (2006) Dorfi, E., Pikall, H., Stökl, A., & Gautschy, A. 2006, Comput. Phys. Commun., 174, 771
  * Dunham & Vorobyov (2012) Dunham, M. M. & Vorobyov, E. I. 2012, ApJ, 747, 52
  * Eggenberger et al. (2022) Eggenberger, P., Buldgen, G., Salmon, S. J. A. J., et al. 2022, Nature Astronomy, 6, 788
  * Elbakyan et al. (2019) Elbakyan, V. G., Vorobyov, E. I., Rab, C., et al. 2019, MNRAS, 484, 146
  * Elsender & Bate (2021) Elsender, D. & Bate, M. R. 2021, MNRAS, 508, 5279
  * Emeriau-Viard & Brun (2017) Emeriau-Viard, C. & Brun, A. S. 2017, ApJ, 846, 8
  * Ercolano & Clarke (2010) Ercolano, B. & Clarke, C. J. 2010, MNRAS, 402, 2735
  * Feiden et al. (2015) Feiden, G. A., Jones, J., & Chaboyer, B. 2015, in Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, Vol. 18, 18th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 171–176
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  * Finley & Matt (2018) Finley, A. J. & Matt, S. P. 2018, ApJ, 854, 78
  * Flaherty et al. (2020) Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109
  * Fossat et al. (2017) Fossat, E., Boumier, P., Corbard, T., et al. 2017, A&A, 604, A40
  * France et al. (2018) France, K., Arulanantham, N., Fossati, L., et al. 2018, ApJS, 239, 16
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1
  * Gallet et al. (2019) Gallet, Zanni, C., & Amard, L. 2019, A&A, 632, A6
  * Gallet & Bouvier (2013) Gallet, F. & Bouvier, J. 2013, A&A, 556, A36
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  * Gehrig, L. et al. (2022) Gehrig, L., Steiner, D., Vorobyov, E. I., & Güdel, M. 2022, A&A, 667, A46
  * Gossage et al. (2021) Gossage, S., Dotter, A., Garraffo, C., et al. 2021, ApJ, 912, 65
  * Guadarrama, R. et al. (2022) Guadarrama, R., Vorobyov, E. I., Rab, Christian, & Güdel, M. 2022, A&A submitted
  * Guarcello et al. (2021) Guarcello, M. G., Biazzo, K., Drake, J. J., et al. 2021, A&A, 650, A157
  * Gullbring et al. (1998) Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323
  * Hartmann (2009) Hartmann, L. 2009, Accretion Processes in Star Formation: Second Edition
  * Hartmann et al. (1997) Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, ApJ, 475, 770
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, Annual Review of Astronomy and Astrophysics, 54, 135
  * Hayashi (1961) Hayashi, C. 1961, PASJ, 13, 450
  * Henderson & Stassun (2012) Henderson, C. B. & Stassun, K. G. 2012, The Astrophysical Journal, 747, 51
  * Hennebelle et al. (2020) Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67
  * Hosokawa et al. (2011) Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Iglesias & Rogers (1993) Iglesias, C. A. & Rogers, F. J. 1993, ApJ, 412, 752
  * Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
  * Ireland et al. (2021) Ireland, L. G., Zanni, C., Matt, S. P., & Pantolmos, G. 2021, ApJ, 906, 4
  * Irwin et al. (2008) Irwin, J., Hodgkin, S., Aigrain, S., et al. 2008, MNRAS, 384, 675
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411
  * Jensen & Haugbølle (2018) Jensen, S. S. & Haugbølle, T. 2018, MNRAS, 474, 1176
  * Johnstone et al. (2021) Johnstone, C. P., Bartel, M., & Güdel, M. 2021, A&A, 649, A96
  * Johnstone et al. (2014) Johnstone, C. P., Jardine, M., Gregory, S. G., Donati, J. F., & Hussain, G. 2014, MNRAS, 437, 3202
  * Kadam et al. (2020) Kadam, K., Vorobyov, E., Regály, Z., Kóspál, Á., & Ábrahám, P. 2020, ApJ, 895, 41
  * King et al. (2007) King, A. R., Pringle, J. E., & Livio, M. 2007, Monthly Notices of the Royal Astronomical Society, 376, 1740
  * Kluźniak & Rappaport (2007) Kluźniak, W. & Rappaport, S. 2007, ApJ, 671, 1990
  * Koenigl (1991) Koenigl, A. 1991, ApJ, 370, L39
  * Koldoba et al. (2002) Koldoba, A. V., Lovelace, R. V. E., Ustyugova, G. V., & Romanova, M. M. 2002, AJ, 123, 2019
  * Königl et al. (2011) Königl, A., Romanova, M. M., & Lovelace, R. V. E. 2011, MNRAS, 416, 757
  * Kunitomo et al. (2017) Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A, 599, A49
  * Kunitomo et al. (2020) Kunitomo, M., Suzuki, T. K., & Inutsuka, S.-i. 2020, MNRAS, 492, 3849
  * Langanke & Martínez-Pinedo (2000) Langanke, K. & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481
  * Larson (1969) Larson, R. B. 1969, MNRAS, 145, 271
  * Lavail et al. (2019) Lavail, A., Kochukhov, O., & Hussain, G. A. J. 2019, A&A, 630, A99
  * Lavail et al. (2017) Lavail, A., Kochukhov, O., Hussain, G. A. J., et al. 2017, A&A, 608, A77
  * Manara et al. (2022) Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2022, arXiv e-prints, arXiv:2203.09930
  * Manara et al. (2012) Manara, C. F., Robberto, M., Rio, N. D., et al. 2012, The Astrophysical Journal, 755, 154
  * Manara et al. (2016) Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3
  * Matt & Pudritz (2005) Matt, S. & Pudritz, R. E. 2005, ApJ, 632, L135
  * Matt et al. (2015) Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23
  * Matt et al. (2010) Matt, S. P., Pinzón, G., de la Reza, R., & Greene, T. P. 2010, ApJ, 714, 989
  * Matt et al. (2012) Matt, S. P., Pinzón, G., Greene, T. P., & Pudritz, R. E. 2012, The Astrophysical Journal, 745, 101
  * Mercer-Smith et al. (1984) Mercer-Smith, J. A., Cameron, A. G. W., & Epstein, R. I. 1984, ApJ, 279, 363
  * Messina et al. (2017) Messina, S., Lanzafame, A. C., Malo, L., et al. 2017, A&A, 607, A3
  * Micela et al. (1985) Micela, G., Sciortino, S., Serio, S., et al. 1985, ApJ, 292, 172
  * Mulders & Dominik (2012) Mulders, G. D. & Dominik, C. 2012, A&A, 539, A9
  * Nakatani et al. (2018) Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018, ApJ, 857, 57
  * Nayakshin & Lodato (2012) Nayakshin, S. & Lodato, G. 2012, MNRAS, 426, 70
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231
  * Palla & Stahler (1992) Palla, F. & Stahler, S. W. 1992, ApJ, 392, 667
  * Pallavicini et al. (1981) Pallavicini, R., Golub, L., Rosner, R., et al. 1981, ApJ, 248, 279
  * Pantolmos et al. (2020) Pantolmos, G., Zanni, C., & Bouvier, J. 2020, A&A, 643, A129
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25
  * Pizzolato et al. (2003) Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147
  * Pollack et al. (1985) Pollack, J. B., McKay, C. P., & Christofferson, B. M. 1985, Icarus, 64, 471
  * Pols et al. (1995) Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964
  * Potekhin & Chabrier (2010) Potekhin, A. Y. & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82
  * Pringle (1981) Pringle, J. E. 1981, ARAA, 19, 137
  * Ragossnig et al. (2020) Ragossnig, F., Dorfi, E. A., Ratschiner, B., et al. 2020, Computer Physics Communications, 256, 107437
  * Rappaport et al. (2004) Rappaport, S. A., Fregeau, J. M., & Spruit, H. 2004, ApJ, 606, 436
  * Rebull et al. (2020) Rebull, L. M., Stauffer, J. R., Cody, A. M., et al. 2020, AJ, 159, 273
  * Richert et al. (2018) Richert, A. J. W., Getman, K. V., Feigelson, E. D., et al. 2018, MNRAS, 477, 5191
  * Rogers & Nayfonov (2002) Rogers, F. J. & Nayfonov, A. 2002, ApJ, 576, 1064
  * Romanova & Kurosawa (2014) Romanova, M. & Kurosawa, R. 2014, in Astronomical Society of the Pacific Conference Series, Vol. 488, 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013), ed. N. V. Pogorelov, E. Audit, & G. P. Zank, 127
  * Romanova et al. (2002) Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2002, ApJ, 578, 420
  * Romanova et al. (2004) Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2004, ApJ, 610, 920
  * Roquette et al. (2021) Roquette, J., Matt, S. P., Winter, A. J., Amard, L., & Stasevic, S. 2021, MNRAS, 508, 3710
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713
  * Schib et al. (2021) Schib, O., Mordasini, C., Wenger, N., Marleau, G. D., & Helled, R. 2021, A&A, 645, A43
  * Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
  * Serna et al. (2021) Serna, J., Hernandez, J., Kounkel, M., et al. 2021, The Astrophysical Journal, 923, 177
  * Shakura & Sunyaev (1973) Shakura, N. & Sunyaev, R. 1973, $\backslash$aap, 24, 337
  * Sherry et al. (2008) Sherry, W. H., Walter, F. M., Wolk, S. J., & Adams, N. R. 2008, AJ, 135, 1616
  * Sicilia-Aguilar et al. (2016) Sicilia-Aguilar, A., Banzatti, A., Carmona, A., et al. 2016, PASA, 33, e059
  * Stahler (1988) Stahler, S. W. 1988, ApJ, 332, 804
  * Stahler & Palla (2004) Stahler, S. W. & Palla, F. 2004, The Formation of Stars
  * Steindl et al. (2021) Steindl, T., Zwintz, K., Barnes, T. G., Müllner, M., & Vorobyov, E. I. 2021, A&A, 654, A36
  * Steindl et al. (2022) Steindl, T., Zwintz, K., & Vorobyov, E. 2022, Nature Communications, 13, 5355
  * Steiner et al. (2021) Steiner, D., Gehrig, L., Ratschiner, B., et al. 2021, A&A, 655, A110
  * Tadross (2003) Tadross, A. 2003, New Astronomy, 8, 737
  * Testi et al. (2014) Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339
  * Testi et al. (2022) Testi, L., Natta, A., Manara, C. F., et al. 2022, A&A, 663, A98
  * Timmes & Swesty (2000) Timmes, F. X. & Swesty, F. D. 2000, ApJS, 126, 501
  * Vorobyov & Basu (2009) Vorobyov, E. I. & Basu, S. 2009, MNRAS, 393, 822
  * Vorobyov & Basu (2010) Vorobyov, E. I. & Basu, S. 2010, The Astrophysical Journal, 719, 1896
  * Vorobyov & Basu (2015) Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115
  * Vorobyov et al. (2017a) Vorobyov, E. I., Elbakyan, V., Hosokawa, T., et al. 2017a, A&A, 605, A77
  * Vorobyov et al. (2020) Vorobyov, E. I., Khaibrakhmanov, S., Basu, S., & Audard, M. 2020, A&A, 644, A74
  * Vorobyov et al. (2019) Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627, A154
  * Vorobyov et al. (2022) Vorobyov, E. I., Skliarevskii, A. M., Molyarova, T., et al. 2022, A&A, 658, A191
  * Vorobyov et al. (2017b) Vorobyov, E. I., Steinrueck, M. E., Elbakyan, V., & Guedel, M. 2017b, A&A, 608, A107
  * Wright et al. (2011) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27
  * Yasui et al. (2016) Yasui, C., Kobayashi, N., Saito, M., & Izumi, N. 2016, The Astronomical Journal, 151, 115
  * Yasui et al. (2021) Yasui, C., Kobayashi, N., Saito, M., Izumi, N., & Skidmore, W. 2021, The Astronomical Journal, 161, 139
  * Yasui et al. (2010a) Yasui, C., Kobayashi, N., Tokunaga, A. T., Saito, M., & Tokoku, C. 2010a, ApJ, 723, L113
  * Yasui et al. (2010b) Yasui, C., Kobayashi, N., Tokunaga, A. T., Saito, M., & Tokoku, C. 2010b, ApJ, 723, L113
  * Zaire et al. (2017) Zaire, B., Guerrero, G., Kosovichev, A. G., Smolarkiewicz, P. K., & Landin, N. R. 2017, in Living Around Active Stars, ed. D. Nandy, A. Valio, & P. Petit, Vol. 328, 30–37
  * Zanni & Ferreira (2013) Zanni, C. & Ferreira, J. 2013, A&A, 550, A99
  * Zhu et al. (2007) Zhu, Z., Hartmann, L., Calvet, N., et al. 2007, ApJ, 669, 483
  * Zhu et al. (2009) Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009, ApJ, 701, 620

##  Appendix A MESA microphysics

The MESA EOS 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) EOSes.

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), 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 non-spherical symmetry, when the star is rotating
(see Sec. 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>t_{\mathrm{0}}$

Fig. 11 shows the evolution of the stellar radius (Panel a), intrinsic
luminosity (Panel b), and mass (Panel c) for $Z=1\leavevmode\nobreak\
\mathrm{Z_{\odot}}$ (solid lines) and $Z=0.2\leavevmode\nobreak\
\mathrm{Z_{\odot}}$ (dashed lines). The yellow, red, and blue lines indicate
models with $\dot{M}_{\rm start}=[1\times 10^{-9},5\times 10^{-9},1\times
10^{-8}]\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$,
respectively.

Figure 11:  Evolution of the stellar radius (Panel a), intrinsic luminosity
(Panel b), and mass (Panel c) for $Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$
(solid lines) and $Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (dashed
lines). The yellow, red, and blue lines indicate models with $\dot{M}_{\rm
start}=[1\times 10^{-9},5\times 10^{-9},1\times 10^{-8}]\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, respectively.

##  Appendix D Time evolution of radial disk structure

The radial profiles for the model with $\dot{M}_{\mathrm{start}}=1.0\times
10^{-8}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ at
the starting point of the simulation $t_{\mathrm{0}}=2$ Myr, as well as at 5.0
Myr and 8.0 Myr are shown in Fig. 12. The solid (dashed) lines correspond to a
metallicity of $Z=1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$
($Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$), respectively.

Figure 12:  Radial profiles for the model with
$\dot{M}_{\mathrm{start}}=1.0\times 10^{-8}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ at the starting point of the
simulation $t_{\mathrm{0}}=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\leavevmode\nobreak\ \mathrm{Z_{\odot}}$
($Z=0.2\leavevmode\nobreak\ \mathrm{Z_{\odot}}$), respectively. Panels (a)
shows the disk surface density $\Sigma$, panels (b) the disk midplane
temperature $T_{\mathrm{disk}}$, panels (c) the optical depth $\tau$ and
panels (d) the mass flux $\dot{M}_{\mathrm{flux}}$, respectively. The thin
vertical lines in panels (a) indicate the respective position of
$r_{\mathrm{c}}$. The horizontal lines in panels (b) indicate the ambient
temperature $T_{\mathrm{amb}}=20$ K. The horizontal dashed lines in panels (c)
symbolizes $\tau=1$, dividing optical thin regions ($\tau<<1$) from optical
thick regions ($\tau>>1$).

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

Table 3:  Age, stellar mass range, and metallicity of young ($\sim 2-5$ Myr)
clusters used in Sec. 4.1.

Name | Age [Myr] | Stellar mass range [$\mathrm{M_{\odot}}$] | Metallicity
---|---|---|---
NGC 6530 |  $\lesssim 5$(1)(1)footnotemark: $(1)$ |  $0.2-1.9$(2)(2)footnotemark: $(2)$ | [Fe/H] = 0.11(3)(3)footnotemark: $(3)$
NGC 2362 |  $\sim 5$(4)(4)footnotemark: $(4)$ |  $0.1-1.2$(4)(4)footnotemark: $(4)$ | [Fe/H] $\approx 0$(4)(4)footnotemark: $(4)$
NGC 2264 |  $\sim 3$(5)(5)footnotemark: $(5)$ |  $0.2-2.5$(5)(5)footnotemark: $(5)$ | [Fe/H] ¡ 0(6)(6)footnotemark: $(6)$
OSFC |  $4$(7)(7)footnotemark: $(7)$ |  $0.2-1.9$(7)(7)footnotemark: $(7)$ | [Fe/H] $\approx 0$(8)(8)footnotemark: $(8)$
$\sigma$ Ori |  $2-3$(9)(9)footnotemark: $(9)$ |  $0.2-0.5$(10)(10)footnotemark: $(10)$ | [Fe/H] = -0.27(9)(9)footnotemark: $(9)$
Taurus |  $\sim 3$(11)(11)footnotemark: $(11)$ |  $\lesssim 2$(11)(11)footnotemark: $(11)$ | [Fe/H] $\approx 0$(12)(12)footnotemark: $(12)$

(1) Damiani et al. (2019); (2) Henderson & Stassun (2012); (3) Tadross (2003);
(4) Irwin et al. (2008); (5) Affer et al. (2013); (6) Cauley et al. (2012), we
note that no exact metallicity value is given; (7) Serna et al. (2021); (8)
D’Orazi et al. (2009); (9) Sherry et al. (2008); (10) Cody & Hillenbrand
(2010); (11) Rebull et al. (2020); (12) D’Orazi et al. (2011)