11institutetext:  Department of Astrophysics, University of Vienna,
Türkenschanzstrasse 17, A-1180 Vienna, Austria 22institutetext: Research
Institute of Physics, Southern Federal University, Rostov-on-Don 344090,
Russia

# What governs the spin distribution of very young ¡ 1 Myr low mass stars

L. Gehrig and E. I. Vorobyov 111122

(Received ….; accepted ….)

###### Abstract

Context. The origin of the stellar spin distribution at young ages is still
unclear. Even in very young clusters ($\sim 1$ Myr), a significant spread is
observed in rotational periods ranging from $\lesssim 1$ to $\sim 10$ days.

Aims. We aim to study the parameters that might govern the spin distribution
of low mass ($\lesssim 1.0\leavevmode\nobreak\ \mathrm{M_{\odot}}$) stars
during the first million years of their evolution.

Methods.  We compute the evolution and rotational periods of young stars,
using the MESA code, starting from a stellar seed, and take protostellar
accretion, stellar winds, and the magnetic star-disk interaction into account.
Furthermore, we add a certain fraction of the energy of accreted material into
the stellar interior as additional heat and combine the resulting effects on
stellar evolution with the stellar spin model.

Results.  For different combinations of parameters, stellar periods at an age
of 1 Myr range between 0.6 days and 12.9 days. Thus, during the relatively
short time period of 1 Myr, a significant amount of stellar angular momentum
can already be removed by the interaction between the star and its accretion
disk. The amount of additional heat added into the stellar interior, the
accretion history, and the presence of disk and stellar winds have the
strongest impact on the stellar spin evolution during the first million years.
The slowest stellar rotations result from a combination of strong magnetic
fields, a large amount of additional heat, and effective winds. The fastest
rotators combine weak magnetic fields and ineffective winds or result from a
small amount of additional heat added to the star. Scenarios that could lead
to such configurations are discussed. Different initial rotation periods of
the stellar seed, on the other hand, quickly converges and do not affect the
stellar period at all.

Conclusions.  Our model matches up to 90% of the observed rotation periods in
six young ($\lesssim 3$ Myr) clusters. Based on these intriguing results, we
motivate to combine our model with a hydrodynamic disk evolution code to self-
consistently include several important aspects such as episodic accretion
events, magnetic disk winds, internal, and external photo-evaporation. Such a
combined model could replace the widely-used disk-locking model during the
lifetime of the accretion disk and provide valuable insights into the origin
of the rotational period distribution of young clusters.

###### Key Words.:

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

##  1 Introduction

The origin of the stellar spin and angular momentum (AM) distribution of young
stars still poses problems for stellar evolution models. Observations of
stellar rotation in young clusters over the last decade have provided valuable
insights for the evolution of the stellar spin and AM (e.g., Herbst et al.,
2002; Rebull et al., 2006; Irwin et al., 2008; Irwin & Bouvier, 2009;
Hartmann, 2009; Affer et al., 2013; Gallet & Bouvier, 2013; Rebull et al.,
2020; Serna et al., 2021). Especially the early distribution ($\lesssim 10$
Myr) of stellar spin and AM is an important benchmark for gyrochronology
studies (e.g., Barnes, 2007) and can affect stellar high energy radiation of
low mass stars (e.g., Pallavicini et al., 1981; Micela et al., 1985; Wright et
al., 2011; France et al., 2018).

We want to focus on the first Myr of stellar evolution, in which the star is
embedded in an envelope and surrounded by an accretion disk. During the
accretion disk phase, it is assumed that the accretion disk plays an important
role in regulating the stellar spin distribution. Mass and AM can be exchanged
between the star and the disk via magnetic field lines, influencing the
stellar rotation period (e.g., Armitage & Clarke, 1996; Matt & Pudritz, 2005b;
Matt et al., 2010a). Early models propose that the star rotates at the orbital
period of the inner disk boundary (disk-locking, e.g., Ghosh & Lamb, 1979;
Koenigl, 1991). However, because of a limited region, in which the magnetic
field lines can connect to the accretion disk, the disk-locking mechanisms
cannot remove enough AM from the star to counteract accretion and contraction
and does not explain the existence of slowly rotating stars (e.g., Matt &
Pudritz, 2005a, b; Zanni & Ferreira, 2013). Since then additional mechanisms
have been proposed that can remove AM from the star during the accretion disk
phase: the process of accretion itself can drive accretion-powered stellar
winds (APSW, Matt & Pudritz, 2005a) and disk material can be ejection along
magnetic field lines (magnetospheric ejections, MEs, e.g., Zanni & Ferreira,
2013) removing AM from the star.

Even at a very young age of $\sim 1$ Myr, young clusters already show a
significant spread in stellar rotation periods. Observations in the Orion
Nebular Cluster (ONC) (e.g., Herbst et al., 2002; Serna et al., 2021), NGC
2264 (e.g., Venuti et al., 2017), Taurus (Rebull et al., 2020), NGC 6530
(Henderson & Stassun, 2012), and other young clusters (YSOVAR, Rebull et al.,
2014) show rotation periods ranging from $<1$ to $\sim 10$ days. We note that
there is a small number of stars rotating even faster or slower than this
range at a young age (e.g., Rebull et al., 2014). To date, the parameters
responsible for this large spread are not constrained. Bouvier et al. (2014)
argues that processes during the cloud collapse or the embedded phase might
govern the rotational distribution at $\sim 1$ Myr. Later studies have not
solved for the very early stellar spin evolution and adopted the observed
spread from $\sim 1$ day to $\sim 10$ days as initial conditions (e.g., Gallet
et al., 2019; Roquette et al., 2021).

To explore the parameter and processes responsible for the rotational
distribution of low-mass stars ($\lesssim 1\leavevmode\nobreak\
\mathrm{M_{\odot}}$), we combine Modules for Experiments in Stellar
Astrophysics (MESA Paxton et al., 2011, 2013, 2015, 2018, 2019) with a recent
stellar spin evolution model presented by Ireland et al. (2021). Similar to
current stellar evolution models, we start with a stellar seed
($M_{\mathrm{\star}}=0.05\leavevmode\nobreak\ \mathrm{M_{\odot}}$) that
accretes mass from the surrounding disk (e.g., Baraffe et al., 2012; Vorobyov
et al., 2017b; Kunitomo et al., 2017; Steindl et al., 2021). A fraction of the
gravitational energy of the accreted material is deposited into the stellar
interior as additional heat, which can lead to an inflation of the stellar
radius and affect the stellar evolution significantly (e.g., Vorobyov et al.,
2017b; Kunitomo et al., 2017). The amount of additional heat depends on the
fraction of energy deposited in the stellar interior, the region, in which the
energy is deposited, and the accretion rate onto the star. The effect of those
parameters on the stellar spin evolution, however, is not clear and will be
tested in this study. Our model considers the influence of stellar inflation
on the spin evolution of stars. It also takes into account different rotation
rates of the stellar seed for the spin evolution of these stars. Furthermore,
we can test a wide parameter space including the stellar magnetic field
strength, the initial angular velocity distribution of the stellar seed, winds
originating from the star and the disk, and the initial radii of the stellar
seed.

Our results suggest that certain mechanisms and parameters have a major impact
on the early stellar spin evolution, while others are virtually unimportant.
Distinct slow and fast rotating models differ in the rotation period at an age
of 1 Myr by up to two orders of magnitude. A comparison with observed
rotational periods of six young clusters allows us to test these model
results. This study is structured as follows: Sec. 2 summarizes the model
highlights and defines the parameter space used in our simulations. The
results are presented in Sec. 3, discussed, and compared to observations in
Sec. 4. Finally, we summarize our conclusions in Sec. 5.

##  2 Model and parameter space

In this study, we combine version v12778 of MESA with the stellar spin
evolution model of Ireland et al. (2021). In the following, we highlight
important aspects of the respective models and provide the respective
references for additional information. The adopted micro-physics and the
treatment of rotation in MESA are summarized in Appendix A and Appendix B,
respectively. The MESA input files for re-creating our results are publicly
available111will be uploaded to ZENODO after acceptance.

###  2.1 Stellar evolution with MESA

Similar to recent stellar evolution studies (e.g., Kunitomo et al., 2017;
Steindl et al., 2021), we start our simulations with a stellar seed with an
initial stellar mass of $M_{\mathrm{init}}=0.05\leavevmode\nobreak\
\mathrm{M_{\odot}}$222We note that Kunitomo et al. (2017) and Steindl et al.
(2021) start with $M_{\mathrm{init}}=0.01\leavevmode\nobreak\
\mathrm{M_{\odot}}$. The combination of high accretion rates and fast stellar
rotation leads to convergence issues with MESA for stellar masses
$<0.05\leavevmode\nobreak\ \mathrm{M_{\odot}}$, solar metallicity,
$Z_{\star}=0.014$, and a Deuterium abundance of 20 ppm. During the initial
stages of stellar evolution, disk material is accreted onto the star. The
associated energy of the accreted material is partly radiated away in the form
of an accretion luminosity

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

and partly added to the stellar interior as additional heat

$\dot{E}_{\mathrm{add}}=\beta\epsilon\frac{GM_{\star}\dot{M}_{\star}}{R_{\star}}\,,$
(2)

with the stellar mass $M_{\star}$, the stellar radius $R_{\star}$, the
gravitational constant $G$ and the accretion rate onto the star
$\dot{M}_{\star}$. The geometry of the accretion process is described by
$\epsilon$. Accretion from a thin accretion disk is described with
$\epsilon\leq 0.5$ (Hartmann et al., 1997) and we choose $\epsilon=0.5$ in
this study (e.g., Baraffe et al., 2009; Steindl et al., 2021). The amount of
additional energy added to the star is controlled by $\beta$ and is assumed to
depend on the accretion rate (e.g., Baraffe et al., 2012; Vorobyov et al.,
2017b; Kunitomo et al., 2017; Elbakyan et al., 2019; Steindl et al., 2021).
During phases of low accretion rates, the disk is truncated at several stellar
radii by magnetic field lines anchored in the star (e.g., Bessolaz et al.,
2008; Hartmann et al., 2016). Disk material is channeled through so-called
funnels and is accreted onto the star at nearly free-fall speeds. Close to the
stellar surface, a shock front develops, and most of the material’s energy is
radiated away (e.g., Koenigl, 1991). As a result, only a small fraction of the
accreted energy is added to the stellar interior and $\beta$ is small. During
phases of high accretion rates, the disk material is pushed toward the star
and the stellar magnetic field is no longer able to channel the disk material
into the funnel flows. Without the radiative loss of energy at the shocks
close to the stellar surface a larger amount of energy can be added to the
star. Based on previous studies, accretion rates of $\dot{M}_{\star}\gtrsim
10^{-5}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$
result in high values of $\beta=\beta_{\mathrm{up}}$ and accretion rates of
$\dot{M}_{\star}\lesssim 10^{-7}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ result in low values of
$\beta=\beta_{\mathrm{low}}$ (e.g., Baraffe et al., 2012; Vorobyov et al.,
2017b; Jensen & Haugbølle, 2018). In between the value of $\beta$ transitions
smoothly between $\beta_{\mathrm{up}}$ and $\beta_{\mathrm{low}}$ according to

$\beta(\dot{M}_{\star})=\frac{\beta_{\mathrm{up}}+\beta_{\mathrm{low}}}{2}+\frac{\beta_{\mathrm{up}}-\beta_{\mathrm{low}}}{2}\tanh{\frac{\log{\dot{M}_{\star}}+6}{0.5}}\,,$
(3)

with the accretion rates in units of $\mathrm{M_{\odot}\leavevmode\nobreak\
yr^{-1}}$. In this study, we choose $\beta_{\mathrm{low}}=0.005$ and vary
$\beta_{\mathrm{up}}$ between 0.005 (corresponding to the case of cold
accretion, $\beta=\mathrm{const.}=0.005$) and 0.5 (corresponding to the case
of warm or hybrid accretion). $\dot{E}_{\mathrm{add}}$ is distributed linearly
starting from the stellar surface, down to a relative mass coordinate
$M_{\mathrm{add}}$ (e.g., Kunitomo et al., 2017; Steindl et al., 2021). We
choose $M_{\mathrm{add}}$ to be in the range of $0.01-0.1$, meaning that
$\dot{E}_{\mathrm{add}}$ is distributed within the outer $1-10$% of the
stellar mass. $\dot{E}_{\mathrm{add}}$ can cause the star to inflate,
depending on the value of $\beta$ and the position in the stellar interior, up
to which this energy is added $M_{\mathrm{add}}$ (e.g., Baraffe et al., 2009,
2012; Vorobyov et al., 2017b; Kunitomo et al., 2017; Steindl et al., 2021).

The outer, atmospheric boundary condition is chosen according to Kunitomo et
al. (2017). For stellar radii larger than $0.7\leavevmode\nobreak\
\mathrm{R_{\odot}}$, we choose the “photospheric” table with an optical depth
of $\tau_{\mathrm{s}}=2/3$, a surface temperature
$T_{\mathrm{s}}=T_{\mathrm{eff}}$ and the surface pressure $P_{\mathrm{s}}$
obtained from PHOENIX (Hauschildt et al., 1999a, b) and ATLAS9 (Castelli &
Kurucz, 2003) model atmospheres. For smaller stellar radii, we adopt the
“$\tau_{\mathrm{s}}=100$” model, in which $P_{\mathrm{s}}$ and
$T_{\mathrm{s}}$ are calculated at $\tau_{\mathrm{s}}=100$ from ATLAS9 and
COND (Allard et al., 2001) model atmospheres.

###  2.2 Stellar spin evolution

In the scope of our model, the star rotates as a solid body with the stellar
AM, $J_{\star}=I_{\star}\Omega_{\star}$; $I_{\star}$ and $\Omega_{\star}$ are
the stellar moment of inertia and angular velocity $\Omega_{\star}$,
respectively. External torques $\Gamma_{\mathrm{ext}}$ change $J_{\star}$
according to
$\dot{J}_{\star}=I_{\star}\dot{\Omega}_{\star}+\dot{I}_{\star}\Omega_{\star}=\Gamma_{\mathrm{ext}}$.
This results in a stellar spin evolution according to

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

The moment of inertia $I_{\star}$ is calculated by the MESA code and the
external torque contributions $\Gamma_{\mathrm{ext}}$ are calculated according
to the description provided by Ireland et al. (2021). The contributions to
$\Gamma_{\mathrm{ext}}$ is composed of the influence of accretion
$\Gamma_{\mathrm{acc}}$, magnetospheric ejections $\Gamma_{\mathrm{me}}$ and
stellar winds $\Gamma_{\mathrm{w}}$. The internal contributions due to
accretion, mass loss, and changes in the stellar radius are combined in
$\Gamma_{\mathrm{int}}=-\dot{I}_{\star}\Omega_{\star}$. An increase in stellar
mass or the stellar radius results in a positive $\dot{I}_{\star}$ and the
star spins down, for example, during phases of stellar radius inflation. On
the other hand, the star spins up if the star contracts. According to Ireland
et al. (2021), the external torque contributions depend on different stellar
parameters and disk parameters. We want to highlight the general dependencies
and specific trends of the individual torque values. The applicability of the
spin model presented in Ireland et al. (2021) for this study is discussed in
Sec. 4.1. For additional information, we refer to the original article.

Accretion torque $\Gamma_{\mathrm{acc}}$: Material that accretes onto the star
usually adds the AM of the disk material at the inner disk boundary to the
star (e.g., Matt et al., 2010b; Gallet et al., 2019; Ireland et al., 2021),
increasing the stellar AM, and causing a stellar spin-up according to

$\Gamma_{\mathrm{acc}}=|\dot{M}_{\star}|(\Omega(R_{\mathrm{t}})-\Omega_{\star})R_{\mathrm{t}}^{2}\,,$
(5)

with the disk rotation rate at the truncation radius $R_{\mathrm{r}}$,
$\Omega(R_{\mathrm{t}})$333We note that the formulation of
$\Gamma_{\mathrm{acc}}$ in Ireland et al. (2021) do not include the
$-\Omega_{\star}$ term. They take into account stellar rotation periods much
smaller compared to the rotation rates at the inner disk boundaries. At young
ages, however, fast rotation rates close to break-up should be considered and
the additional term should be included.. The truncation radius in units of the
stellar radius can be written as

$\frac{R_{\mathrm{t}}}{R_{\star}}=K_{\mathrm{t}}\left(\frac{B_{\star}^{2}R_{\star}^{2}}{4\pi|\dot{M}_{\star}|v_{\mathrm{esc}}}\right)^{m_{\mathrm{t}}}\,,$
(6)

with $K_{\mathrm{t}}=0.756$, $m_{\mathrm{t}}=0.34$, the stellar dipole field
strength $B_{\star}$, and the escape velocity from the stellar surface
$v_{\mathrm{esc}}=(2GM_{\star}/R_{\star})^{1/2}$. The disk at the truncation
radius is assumed to rotate slower compared to the Keplerian velocity (e.g.,
Bessolaz et al., 2008; Zanni & Ferreira, 2013) due to the removal of AM from
the disk due to the magnetic star-disk interaction (e.g., Zanni & Ferreira,
2013) or disk winds (e.g., Zhang et al., 2017). Thus, $\Omega(R_{\mathrm{t}})$
can be written as

$\Omega(R_{\mathrm{t}})=K_{\mathrm{acc}}\left(\frac{R_{\mathrm{t}}}{R_{\star}}\right)^{m_{\mathrm{acc}}}\Omega_{\mathrm{K}}(R_{\mathrm{t}})\,,$
(7)

with the correction factor for the non-Keplerian rotation
$K_{\mathrm{acc}}=0.775$, $m_{\mathrm{acc}}=-0.147$, and the Keplerian
rotation rate at the truncation radius $\Omega_{\mathrm{K}}(R_{\mathrm{t}})$.

Magnetic star-disk interaction: The influence of the magnetic star-disk
interaction in Ireland et al. (2021) is modeled as magnetospheric ejections
(MEs, Matt & Pudritz, 2005b). MEs can transfer angular momentum between the
star and the disk due to the twist of magnetic field lines that are rooted at
the stellar surface and coupled to the disk material (e.g., Armitage & Clarke,
1996; Matt & Pudritz, 2005b; Zanni & Ferreira, 2013; Ireland et al., 2021).
The amount of angular momentum transferred to or away from the star depends on
the relation between the truncation radius $R_{\mathrm{t}}$ and the co-
rotation radius

$R_{\mathrm{co}}=\left[K_{\mathrm{acc}}\left(\frac{R_{\mathrm{t}}}{R_{\star}}\right)^{m_{\mathrm{acc}}}\right]^{2/3}\left(\frac{GM_{\star}}{\Omega^{2}}\right)^{1/3}\,,$
(8)

which includes the effect of non-Keplerian rotation. Disk material, which is
coupled to the stellar magnetic field lines inside (outside) $R_{\mathrm{co}}$
rotates faster (slower) compared to the star and increases (decreases) the
stellar angular momentum. This process is highly dynamic and complex. Thus,
the following formulation should be treated with caution and interpreted as a
time-averaged approximation. Following Ireland et al. (2021), the total torque
due to MEs can be written as

$\Gamma_{\mathrm{ME}}=K_{\mathrm{ME}}B_{\star}^{2}R_{\star}^{3}\left(\frac{R_{\mathrm{t}}}{R_{\star}}\right)^{m_{\mathrm{ME}}}\left[1-\left(\frac{R_{\mathrm{t}}}{R_{\mathrm{co}}}\right)^{3/2}\right]\,,$
(9)

with numerical fit parameters $K_{\mathrm{ME}}=0.00772$ and
$m_{\mathrm{ME}}=-2.54$. If $R_{\mathrm{t}}<R_{\mathrm{co}}$
($R_{\mathrm{t}}>R_{\mathrm{co}}$), the total torque due to MEs is positive
(negative) and the star spins up (spins down).

Accretion-powered stellar wind (APSW): In the picture of an accretion-powered
stellar wind, a certain fraction $W$ of the accretion rate is ejected by the
star; $\dot{M}_{\mathrm{w}}=W\dot{M}_{\star}$. According to Weber & Davis
(1967), the torque due to such a stellar outflow can be written as

$\Gamma_{\mathrm{w}}=-|\dot{M}_{\mathrm{w}}|\Omega_{\star}r_{\mathrm{A}}^{2}\,,$
(10)

with the Alfvén radius $r_{\mathrm{A}}$. The Alfvén radius can be related to
the open stellar magnetic flux (Réville et al., 2015) and based on numerical
best-fit parameters (Ireland et al., 2021), the torque due to an APSW can be
written as

$\Gamma_{\mathrm{w}}=-K_{\mathrm{w1}}|\dot{M}_{\star}|(GM_{\star}R_{\star})^{1/2}W^{m_{\mathrm{w1}}}f^{m_{\mathrm{w2}}}\left[1+\left(\frac{f}{K_{\mathrm{w2}}}\right)^{2}\right]^{-m_{\mathrm{w3}}}\left(\frac{R_{\mathrm{t}}}{R_{\star}}\right)^{m_{\mathrm{w4}}}\,,$
(11)

with best-fit parameters $K_{\mathrm{w1}}=57.92$, $m_{\mathrm{w1}}=0.254$,
$m_{\mathrm{w2}}=1.092$, $K_{\mathrm{w2}}=0.0356$, $m_{\mathrm{w3}}=0.373$,
$m_{\mathrm{w4}}=0.195$, and the fraction of the break-up rotation rate
$f=\Omega_{\star}/\Omega_{\mathrm{crit}}$. The torque due to stellar winds
$\Gamma_{\mathrm{w}}$ always removes AM from the star.

###  2.3 Effects of different stellar rotation rates

Based on this spin evolution model presented above, we can predict the
behavior of a star that rotates at a certain rate $\Omega_{\mathrm{old}}$ and
spins up to $\Omega_{\mathrm{new}}>\Omega_{\mathrm{old}}$, with all other
parameters unchanged. First, the accretion torque $\Gamma_{\mathrm{acc}}$ is
reduced (see Eq. 5). Second, the torque due to the magnetic star-disk
interaction $\Gamma_{\mathrm{ME}}$ is smaller as the co-rotation radius moves
closer to the truncation radius for fast rotation rates (see Eq. 9). Third,
the rate of angular momentum that is removed by the APSW,
$\Gamma_{\mathrm{w}}$, increases (see Eq. 10). An increasing stellar rotation
rate ($\Omega_{\mathrm{new}}>\Omega_{\mathrm{old}}$) would thus result in a
spin-down torque and the star would return to $\Omega_{\mathrm{old}}$. Based
on the same principles, a decreasing stellar rotation rate
($\Omega_{\mathrm{new}}<\Omega_{\mathrm{old}}$) would result in a spin-up
torque. Hence, $\Omega_{\mathrm{new}}$ evolves back to
$\Omega_{\mathrm{old}}$, and the timescale, on which this change occurs is
studied in Sec. 3.2.

###  2.4 Accretion rate onto the star

In this study, we start from an initial mass of $M_{\mathrm{init}}=0.05$ and
consider three different stellar masses at an age of 1 Myr;
$M_{\mathrm{final}}=0.3,\leavevmode\nobreak\ 0.7,\leavevmode\nobreak\
\mathrm{and}\leavevmode\nobreak\ 1.0\leavevmode\nobreak\ \mathrm{M_{\odot}}$.
The accretion rate is defined in two regimes. Initially, we assume a high,
constant accretion rate of $\dot{M}_{\mathrm{init}}=5\times
10^{-6}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$,
which is equivalent to the mass infall rate $\dot{M}_{\mathrm{infall}}\propto
c_{\mathrm{s}}^{3}/G$ with the isothermal sound speed $c_{\mathrm{s}}$ and an
external temperature of $\approx 21$ K (e.g., Vorobyov et al., 2017b). Thus,
an initial stellar mass of $M_{\mathrm{init}}=0.05$ corresponds to an age of
$\sim 10^{4}$ years, assuming a roughly constant accretion rate during this
time. After a specific time $t_{\mathrm{acc}}$, the accretion rate decreases
according to a power law with an exponent of -1.5, $\dot{M}_{\star}(t)\propto
t^{-1.5}$. A stellar mass of 0.3, 0.7, and 1.0 $\mathrm{M_{\odot}}$ can be
obtained at 1 Myr for $t_{\mathrm{acc}}=21.4$, 54.5 and 81.5 kyr,
respectively. Fig. 1 shows the evolution of the stellar accretion rate (Panel
a) and the stellar mass (Panel b) for a final stellar mass of 0.3, 0.7 and 1.0
$\mathrm{M_{\odot}}$ (solid colored lines). To study the influence of a
different accretion history, we increase the initial accretion rate to
$\dot{M}_{\mathrm{init}}=2\times 10^{-5}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ and decrease the power-law
exponent to $-2$ ($\dot{M}_{\mathrm{high}}$ dashed line). The final stellar
mass at 1 Myr is also 0.7 $\mathrm{M_{\odot}}$ in case of using
$\dot{M}_{\mathrm{high}}$.

Figure 1:  Panel (a) shows the accretion rate and Panel (b) shows the stellar
mass over time. The yellow, black, and blue lines represent a final stellar
mass of 0.3, 0.7, and 1.0 $\mathrm{M_{\odot}}$, respectively. The black dashed
line ($\dot{M}_{\mathrm{high}}$) represents a different accretion history,
which also results in a final stellar mass of 0.7 $\mathrm{M_{\odot}}$ at 1
Myr.

###  2.5 Parameter space

Finding reasonable ranges for stellar and accretion parameters for very young,
embedded objects is challenging due to theoretical and observational
uncertainties and a limited number of observations. The ranges for the free
parameters presented below are based on previous studies, theoretical models,
and observational results.

Magnetic field strengths: One important stellar parameter that influences the
spin evolution model presented above is the magnetic field strength. It is
assumed that the large-scale dipole field components spin evolution of young
stars (e.g., Finley & Matt, 2018). Small-scale field components, on the other
hand, play only a minor role. We note that the results presented in Finley &
Matt (2018) do not cover very young, embedded objects with high accretion
rates (see Sec. 4.1). Due to observational limitations, there is only a
limited number of magnetic field estimates for Class I (0) objects. Based on
the observation of magnetically-sensitive line broadening (Zeeman effect), a
total magnetic field strength of $2.9\pm 0.4$ kG for WL 17 (Johns-Krull et
al., 2009) and $1.36\pm 0.06$ kG for V347 Aur (Flores et al., 2019) can be
inferred. Unfortunately, the total field strength combines the small and
large-scale features of the magnetic field and a distinction is only possible
with additional (more precise) measurement (e.g., Zeeman-Doppler imaging, ZDI,
Johnstone et al., 2014). Such observations are currently not available for
embedded Class I (0) objects. In addition, the small number of available
measurements does not have to be representative. There are, however, estimates
on the magnetic field strength of Class I (0) objects ranging from $\sim 0.1$
kG to $\sim 1$ kG based on flux conservation during the core collapse
(Tsukamoto et al., 2022), near-infrared K-band spectra (Laos et al., 2021),
and laboratory plasma experiments (Burdonov et al., 2021). These estimates
match the inferred field strengths for WL 17 and V347 Aur. Furthermore, Laos
et al. (2021) argue that Class I (0) objects should have similar accretion
mechanisms and magnetic field strengths to Class II objects. For these
slightly further evolved objects, the large-scale field component ranges from
$\sim 0.1$ kG to $\gtrsim 1$ kG (e.g., Johnstone et al., 2014). Thus, we
choose a range for the large-scale magnetic field strength, $B_{\star}$, from
0.1 kG to 2.0 kG. These values are kept constant throughout the simulations in
Sec. 3. In addition, we explore the possibility of mass-dependent magnetic
field strengths as presented in Browning et al. (2016) (see Sec. 4.5).

Initial rotation rates: While there are rotational studies of young ($\sim 1$
Myr) T Tauri stars in NGC 2264 (Venuti et al., 2017), ONC (Rebull et al.,
2006), $\sigma$ Ori (Cody & Hillenbrand, 2010), and Taurus (Rebull et al.,
2020) , there are no large-scale surveys of stellar rotation periods of Class
I (0) objects. Similar to magnetic fields, there are only singular estimates
for individual stars. The Class I object Elias 29 has a rotation period of
$\sim 2$ days (Pillitteri et al., 2019). Furthermore, White et al. (2007)
suggest that Class I objects rotate slightly faster than their Class II
counterparts. This would result in rotation periods of $\sim 1$ day during the
Class I phase. The initial rotation periods at the beginning of the Class I
(0) phase, however, are not well constrained. Thus, we explore a wide range of
initial rotation rates from very slow-rotating stars with
$\Omega_{\mathrm{init}}=0.01\Omega_{\mathrm{crit}}$ to fast-rotating stars
with $\Omega_{\mathrm{init}}=0.9\Omega_{\mathrm{crit}}$, with
$\Omega_{\mathrm{crit}}$ being the critical or brake-up rotation rate, which
is usually estimated as (e.g., Maeder, 2009)

$\Omega_{\mathrm{crit}}=\sqrt{\frac{GM_{\star}}{(1.5R_{\star})^{3}}}\,.$ (12)

Initial stellar radius: Studies of early stellar evolutionary stages that
start shortly after the second collapse and the formation of the second
hydrostatic core (with an initial stellar mass of $0.01\leavevmode\nobreak\
\mathrm{M_{\odot}}$) assume stellar radii of $0.25\leavevmode\nobreak\
\mathrm{R_{\odot}}$ to $3.0\leavevmode\nobreak\ \mathrm{R_{\odot}}$ (e.g.,
Hosokawa et al., 2011; Kunitomo et al., 2017; Steindl et al., 2021). According
to Kunitomo et al. (2017), initial stellar radii of $\gtrsim
1.5\leavevmode\nobreak\ \mathrm{R_{\odot}}$ are required to explain the most
luminous stars and small initial radii can explain the low luminosities, and
the observed luminosity spread in young clusters (Hillenbrand, 2009; Soderblom
et al., 2014). Although we start with slightly further evolved initial stellar
models (with an initial mass of $0.05\leavevmode\nobreak\
\mathrm{M_{\odot}}$), we assume similar initial stellar radii ranging from
$0.6\leavevmode\nobreak\ \mathrm{R_{\odot}}$ to $2.5\leavevmode\nobreak\
\mathrm{R_{\odot}}$.

Accretion powered stellar wind efficiency: In this study, winds originating
from the star are represented by accretion-powered stellar winds (APSW). Such
winds extract a certain amount of angular momentum from the star and cause a
spin-down (see Eq. 10). The exact driving mechanisms of these winds are still
under debate (see, e.g., the discussion in Gallet et al., 2019). We adapt the
driving mechanisms proposed by Decampli (1981) and Matt & Pudritz (2005a).
Angular momentum can be removed from the star by Alfvén waves, excited by the
accreting material at the stellar surface and traveling outwards. The mass-
loss rates of such winds are assumed to be in the order of one percent of the
accretion rate onto the star ($W\sim 1$%, e.g., Cranmer, 2008; Pantolmos et
al., 2020). Thus, the APSW efficiency value $W$ ranges from 0% to 5% in this
study. We note that there are also other wind-driving mechanisms affecting the
star-disk system (e.g., jets launched at the inner disk). We discuss the
possible impact of such winds or jets on the rotational evolution of the young
star in Sec. 3.4.

###  2.6 Reference model

We define a reference model, to qualitatively show how each parameter affects
the rotational evolution of young stars. In the following sections, the
reference values are marked in boldface. The reference parameters are
summarized in Table 1.

Table 1: Reference model parameters used in this study. We note that the stellar mass corresponds to an age of 1 Myr. $M_{\star}$ | $0.7\leavevmode\nobreak\ \mathrm{M_{\odot}}$ | $\Omega_{\mathrm{init}}$ | $0.1\leavevmode\nobreak\ \mathrm{\Omega_{crit}}$
---|---|---|---
$B_{\star}$ | 1.0 kG | $R_{\mathrm{init}}$ | $1.1\leavevmode\nobreak\ \mathrm{R_{\odot}}$
$W$ | 2% | $\beta_{\mathrm{up}}$ | 0.2
$M_{\mathrm{add}}$ | 0.1 |  |

##  3 Results

###  3.1 Evolution of the reference value

Before going into detail on the individual parameters, we want to show the
evolution of the reference model (see Fig. 2). The choice of initial
conditions leads to an initial rotation period of 6 days. During the first few
thousand years, the star spins up to $\approx 4$ days. This initial spin-up
can be explained by adjusting the rotation period to the accretion rate when
starting the simulation. Starting from the initial stellar radius of
$R_{\mathrm{init}}=1.1\leavevmode\nobreak\ \mathrm{R_{\odot}}$, the star
inflates due to the energy of the accreting material added into the stellar
interior and the onset of Deuterium burning. Up to an age of $\approx 18$ kyr,
the inflation of the star (combined with the APSW) can counteract the spin-up
influence of accretion and the star-disk interaction and the star spins down
towards 6 days ($\Gamma_{\mathrm{tot}}<0$, in Panel c of Fig. 2). With the
increasing stellar radius, the spin-up influence of the MEs increases as well,
and the stellar period decreases toward 4 days up to an age of 80 kyr. At
$\sim 80$ kyr, the stellar inflation increases for a short period of time,
resulting in a brief spin-down. Afterward, the star contracts for the
remaining simulation, reaching a final stellar radius of
$1.84\leavevmode\nobreak\ \mathrm{R_{\odot}}$ after 1 Myr. During contraction,
$\Gamma_{\mathrm{int}}$ increases the stellar spin up and the star reaches a
period of 1.6 days after 1 Myr.

Figure 2:  Evolution of the reference model. The stellar rotation period and
radius are shown in Panels (a) and (b), respectively. The different torque
contributions are shown in Panel (c). The accretion torque,
$\Gamma_{\mathrm{acc}}$, the influence of the magnetic star-disk interaction,
$\Gamma_{\mathrm{ME}}$, the APSW torque, $\Gamma_{\mathrm{W}}$, and the
internal contributions due to variation of the stellar moment of inertia,
$\Gamma_{\mathrm{int}}$, are shown as a solid blue line, dashed blue line,
dotted blue line, and dash-dotted blue line, respectively. The total torque,
$\Gamma_{\mathrm{tot}}$, defined as the sum of the other torque contributions
is shown as a solid black line. A grey dashed line marks the change between
the spin-up and spin-down regime.

###  3.2 Effects of individual parameters

Now we examine how the stellar rotational evolution is affected by the initial
rotation rate $\Omega_{\mathrm{init}}$, the stellar magnetic field strength
$B_{\star}$, the APSW efficiency $W$, the initial stellar radius
$R_{\mathrm{init}}$, the final stellar mass at 1 Myr $M_{\star}$, and the
amount of energy of the infalling material added into the star,
$\beta_{\mathrm{up}}$ (see Fig. 3). For a clear description of different
effects, we show, in addition to the stellar rotation period, the evolution of
the stellar radius in Fig. 4, the individual torque contributions in Appendix
C, and the HRD for each model in Fig. 15.

The initial rotation rate: First, we want to show the effects of different
initial rotation rates, $\Omega_{\mathrm{init}}$ (Panel a in Fig. 3 and Fig.
4). We vary $\Omega_{\mathrm{init}}$ from 1% to 90% of the critical (or break-
up) velocity. The black, blue, and yellow lines symbolize the reference, fast,
and slow rotating models. Within the first $\sim 5$ kyr, the initial
difference in the stellar periods of $\approx 59$ days converges towards a
common value (see Sec. 2.3). As a consequence of the fast initial rotation and
strong centrifugal forces, the initial stellar radius for a star with
$\Omega_{\mathrm{init}}=0.9\leavevmode\nobreak\ \mathrm{\Omega_{crit}}$ is
slightly larger compared to the other models. However, this difference
disappears after a few thousand years and the stellar evolution toward 1 Myr
is completely unaffected by the initial rotation rate
$\Omega_{\mathrm{init}}$.

Stellar magnetic field strength: The stellar magnetic field strength does
affect all external torque contributions, such as $\Gamma_{\mathrm{acc}}$,
$\Gamma_{\mathrm{ME}}$, and $\Gamma_{\mathrm{W}}$. The torque due to MEs,
however, is affected most strongly (see Eq. 9 in Sec. 2.2). As shown in Panel
(b) of Fig. 3, a strong magnetic field causes the star to spin up during the
first $\sim 100$ kyr. After $\sim 100$ kyr, a strong magnetic field causes a
spin down. The strong and weak magnetic field strengths are symbolized by the
blue and yellow lines, respectively. This behavior can be explained with Eq.
(9). During the initial $\sim 100$ kyr, the star is still heavily embedded and
the accretion rates can be assumed to be high (see Fig. 1). The truncation
radius, $R_{\mathrm{t}}$, is pushed inwards towards the star (regardless of
the magnetic field strengths used in this study), and the ratio between
$R_{\mathrm{t}}$ and the co-rotation radius, $R_{\mathrm{co}}$, is small.
Thus, the torque due to the magnetic star-disk interaction,
$\Gamma_{\mathrm{ME}}$, is positive, and a stronger magnetic field strength
results in a stronger spin-up torque. After $\sim 100$ kyr, the accretion rate
decreases, and a stronger magnetic field pushes the truncation radius outwards
closer to the co-rotation radius. The disk’s rotation at the truncation
radius, $\Omega(R_{\mathrm{t}})$, and the difference between
$\Omega(R_{\mathrm{t}})-\Omega_{\star}$ decreases. Consequently, stronger
magnetic fields result in a weaker accretion spin-up torque (see Eq. 5), a
weaker torque due to the magnetic star-disk interaction,
$\Gamma_{\mathrm{ME}}$ (see Eq. 9), and a slower rotation period compared to
weak magnetic field strengths. The differences in the rotation periods are
relatively small ($\sim 1$ day) and the stellar radius evolution is almost
unaffected (see Penal b of Fig. 4).

APSW efficiency $W$: In the current spin evolution model, a stellar wind
removes angular momentum from the star and causes a spin-down. The amount of
AM that is removed scales with the mass loss, which is proportional to the
APSW efficiency $W$ (see Eq. 10). Thus, a higher efficiency value results in a
stronger spin-down torque, and the star slows down compared to a weaker or no
APSW (see Panel d in Fig. 3). The difference, especially during the early
evolution ($\lesssim 10^{5}$ years) can be significant. With no angular
momentum removed by an APSW ($W=0$%), the star spins up to $81.7$% of the
critical rotation rate (see the dotted lines in Fig. 3). The stellar rotation
period differs by more than 7 days, equivalent to a factor of $\approx 4$. The
stellar radius is also affected by the different rotation periods. A fast-
rotating star shows a larger radius (Panel c in Fig. 4) and a lower effective
temperature (see Panel c in Fig. 15). Thus, the first nuclear reactions ignite
at a later age (see the shifted radius evolution in Panel c of Fig. 4). While,
toward an age of 1 Myr, the rotation periods still show a difference of a
factor of $\sim 2$, the stellar radii converge to a similar value.

Maximum amount of additional heat $\dot{E}_{\mathrm{add}}$: The amount of
additional heat can be controlled by the free parameter
$\beta_{\mathrm{up}}$444We expect that the stellar evolution is largely
independent of $\beta_{\mathrm{low}}$, as long the value is chosen to be small
enough, $\beta_{\mathrm{low}}\ll 0.1$.. The value of $\beta_{\mathrm{up}}$
does have a significant impact on stellar evolution. While the evolution of
most stars can be explained by a substantial amount of additional heat added
to the star ( corresponding to $\beta_{\mathrm{up}}\gtrsim 0.1$,e.g., Vorobyov
& Pavlyuchenkov, 2017; Kunitomo et al., 2017; Steindl et al., 2021), objects
with very low luminosity require $\beta\sim\beta_{\mathrm{low}}$ (Kunitomo et
al., 2017; Vorobyov et al., 2017a). The amount of additional heat affects the
evolution of the stellar radius. The more energy added to the stellar
interior, the stronger the star inflates (see Panel d in Fig. 4). In the cold-
accretion scenario ($\beta=\beta_{\mathrm{low}}$), the stellar radius does not
inflate and remains approximately constant at a value of $\sim
1\leavevmode\nobreak\ \mathrm{R_{\odot}}$. These different evolution tracks
influence stellar rotation (see Panel d in Fig. 3). The strong initial
inflation to $R_{\star}>3\leavevmode\nobreak\ \mathrm{R_{\odot}}$ for
$\beta_{\mathrm{up}}=0.5$ results in a spin-down to 9.2 days. With larger
radii, MEs become more effective and start to spin up the star slowly. After 1
Myr, the stellar period has reached $\sim 2$ days. Without the initial
inflation, the star continuously spins up during its evolution (excluding the
short period of inflation due to the first nuclear reactions). After 1 Myr,
the stellar period has reached $\sim 0.57$ days. Within the parameter ranges
studied so far, a variation in $\beta_{\mathrm{up}}$ marks the largest impact
on the stellar rotation period.

Stellar mass at 1 Myr: To reach different stellar masses at 1 Myr, we vary the
age, at which the initially constant accretion rate decreases (see Fig. 1).
During the first $\sim 10^{5}$ years, smaller final stellar masses result in
faster-rotating stars (see Panel e in Fig. 3). To explain this behavior, the
amount of additional heat has to be considered. For a decreasing accretion
rate, $\beta_{\mathrm{up}}$ and consequently $\dot{E}_{\mathrm{add}}$ decrease
as well. Thus, the stellar inflation stops and the star starts contracting,
resulting in a stellar spin-up. After $10^{5}$ years, the influence of the
additional heat diminishes as accretion rates decrease further. Now, a lower
accretion rate results in less accretion torque, $\Gamma_{\mathrm{acc}}$, and
a smaller stellar radius results in less effective MEs. As a consequence, a
smaller accretion rate causes a slower rotation compared to high accretion
rates.

Initial stellar radii: Finally, we show the impact of different initial
stellar radii on the spin evolution. For a given stellar mass, a larger radius
results in a cooler temperature throughout the star. Starting at a large
initial stellar radius, gravity prevails over the inflation due to the
additional heat, and the star contracts slightly (see Panel f in Fig. 4).
After $\sim 20$ kyr, the inflation of the stellar radius has balanced the
initial differences. The temperature deep in the stellar interior, however, is
still cooler in case of an initially larger radius and the first nuclear
reactions start slightly later (maximum radii are reached later). Reaching 1
Myr, the stellar radii converge toward a common value and the initial
differences are erased. A similar picture emerges for the stellar rotation
periods. Initial differences converge fast and throughout the evolution, only
small differences arise due to different stellar radii. After 1 Myr, the
rotation period evolved independently of the initial stellar radius. In our
reference case, we choose $\beta_{\mathrm{up}}=0.2$, which corresponds to a
warm accretion scenario. Kunitomo et al. (2017), however, pointed out that the
importance of the initial stellar radius increases, when decreasing the amount
of additional heat. Hence, additional simulations have been carried out with
different initial radii and $\beta_{\mathrm{up}}=0.005$ (see Appendix D).
Without the additional heat and the resulting stellar inflation, the stellar
radii maintain a significant difference throughout their evolution and after 1
Myr, the stellar radii differ by a factor of $\sim 2$ (see Panel b in Fig.
14). As a result, the stellar rotation period does not converge and shows a
dependence on the initial stellar radius. Smaller stellar radii result in
faster-rotating stars and after 1 Myr, the rotation periods differ by a factor
of 2.6 (see Panel a in Fig. 14).

Figure 3:  Evolution of the stellar rotation period, $P_{\star}$, for each
parameter variation: the initial rotation rate $\Omega_{\mathrm{init}}$ (Panel
a), the stellar magnetic field strength $B_{\star}$ (Panel b), the APSW
efficiency $W$ (Panel c), the amount of energy of the infalling material added
into the star, $\beta_{\mathrm{up}}$ (Panel d), the final stellar mass at 1
Myr, $M_{\star}$ (Panel e), and the initial stellar radius,
$R_{\mathrm{init}}$ (Panel f). The reference values are marked in boldface.
For comparison, the critical rotation period is shown for each evolution as
dotted lines.

Figure 4:  Evolution of the stellar radius, $R_{\star}$, for each parameter
variation (same as in Fig. 3). The reference values are marked in boldface.

###  3.3 Other factors affecting $\dot{E}_{\mathrm{add}}$

The additional heat, $\dot{E}_{\mathrm{add}}$, is not only controlled by
$\beta$. The region, in which $\dot{E}_{\mathrm{add}}$ is added, and the
accretion history can influence the stellar evolution. We reduce the region,
in which the additional heat, $\dot{E}_{\mathrm{add}}$, is distributed to
$M_{\mathrm{add}}=0.01$. With more additional heat injected close to the
stellar surface, the initial stellar inflation increases (see the sharp
increase of the stellar radius in Column a of Fig. 5). During the whole
evolution, the stellar radius remains slightly larger compared to the
reference model. Reaching 1 Myr, the stellar radius converges toward the
reference value as the influence of $\dot{E}_{\mathrm{add}}$ becomes less
important as the accretion rate decreases. As discussed before, larger
inflated stars tend to rotate more slowly. During the strong initial
inflation, the difference to the reference model is the largest ($\sim 3$
days). Similar to the stellar radius, reaching an age of 1 Myr, the stellar
period converges toward the reference value.

Next, we change the accretion history. Starting from an initial accretion rate
of $2\times 10^{5}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\
yr^{-1}}$, which is factor 4 larger compared to the reference value, the
accretion rate drops earlier and faster to reach the same stellar mass
($0.7\leavevmode\nobreak\ \mathrm{M_{\odot}}$) compared to the reference model
(see Fig. 1). The strong accretion rate injects a large amount of additional
heat into the stellar interior and the star inflates rapidly to
$>4\leavevmode\nobreak\ \mathrm{R_{\odot}}$ (see Column a of Fig. 5). Due to
this immense expansion, the star slows down to over 10 days. After the
accretion rate drops, the stellar expansion stops and the star contracts
toward the reference value, reaching the end of the simulation. In contrast to
the previous model with $M_{\mathrm{add}}=0.01$, the smaller accretion rate at
ages $>10^{5}$ years reduces the accretion torque and the torque due to MEs
and the star spins down. At the age of 1 Myr, the stellar rotation period
differs by a factor of 1.9.

While the region, in which the additional heat is injected, plays only a minor
role in the stellar spin evolution, the accretion history significantly
affects the stellar period. In our model, the accretion rate is assumed to be
smooth, without short-term variations. This idealistic assumption, however,
should merely be considered as an average over time. In hydro-dynamic disk
simulations, especially the embedded phase ($<10^{5}$ years), show episodic
outbursts with peak accretion rates of $>10^{5}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ (e.g., Vorobyov & Basu, 2015;
Vorobyov et al., 2017b; Hsieh et al., 2019; Kadam et al., 2021). During such
an outburst, a huge amount of additional heat is injected into the stellar
interior, the star inflates, and the rotation period increases. The nature of
these outbursts depends, amongst other things, on the stellar mass. More
massive stars undergo more and stronger outbursts over a larger period of
their early life. If each outburst spins down the star slightly, more massive
stars should tend to rotate more slowly compared to low-mass stars. This trend
is, in fact, reported in Henderson & Stassun (2012) for young ($\lesssim 5$
Myr) clusters.

###  3.4 Disk winds

In the presence of a large-scale magnetic field in the accretion disk, disk
winds, launched by magneto-rotational processes in the inner disk regions, can
remove a significant amount of angular momentum from the disk (e.g., Zhang et
al., 2017). The exact regions, in which the wind is launched, and the total
amount of mass and angular momentum transported away, depend on a variety of
factors (e.g., Guilet & Ogilvie, 2014, Steiner et al. in prep). A detailed
description of this process would go beyond the scope of this work. However,
we can approximate the effect of a strong disk wind on the stellar rotation.
The loss of angular momentum in the inner disk region causes the disk material
to slow down and reduces the accretion torque. For our simulation, we extend
Eq. 7 by a factor of $K_{\mathrm{DW}}=0.5$. We note that this factor is most
likely not constant and depend on different parameters, but for showing the
qualitative effect of a disk wind, this assumption should be viable. The
rotational evolution including a disk wind is shown in Column (b) of Fig. 5
(yellow line). As expected, the reduced accretion torque leads to a rapid
spin-down during the radius inflation. Reaching 1 Myr, the rotational period
is still slower by a factor of 2.0 compared to the reference model, which is
similar to the effect of the model with $\dot{M}_{\mathrm{high}}$.

###  3.5 Stellar metallicity

Finally, stellar metallicity can affect the rotational evolution of the star
during the pre-main sequence and the main sequence. Low-metallicity stars are
expected to rotate faster compared to their solar-metallicity counterparts
(e.g., Amard et al., 2019; Amard & Matt, 2020; Amard et al., 2020; Gehrig et
al., 2023). As low-metallicity stars are more compact compared to solar
metallicities, mechanisms that remove angular momentum from the star and scale
with the stellar radius are less effective in a low-metallicity environment
during the Class II phase and beyond the main sequence (e.g., Amard & Matt,
2020; Gehrig et al., 2023). In addition, the disk lifetime of low-metallicity
stars is shorter compared to solar metallicities, and the pre-main sequence
spin-up due to contraction starts at an earlier age (Yasui et al., 2016, 2021;
Guarcello et al., 2021; Gehrig et al., 2023). We find similar results for the
evolution during the first million years for a metallicity of
$Z_{\star}=0.1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (see Column b in Fig.
5, blue line). The initial inflation is weaker and the star is more compact
during the whole simulation time compared to the reference model. Thus, the
star spins up during the embedded phase. Reaching the age of 1 Myr, the APSW
is less effective in removing angular momentum from the star due to the
smaller stellar radius and the star continues to spin up compared to the
reference model.

Figure 5:  Evolution of the stellar rotation period, $P_{\star}$, and the
stellar radius, $R_{\star}$, for different variations from the reference
model. Column (a): region, in which $\dot{E}_{\mathrm{add}}$ is distributed,
$M_{\mathrm{add}}=0.01$ (blue line) and using a different accretion history,
$\dot{M}_{\mathrm{high}}$ (see Fig. 1, yellow line). Column (b): low
metallicity, $Z_{\star}=0.1\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (blue
line), and including the influence of disk winds (model DW, yellow line). The
reference model is marked in boldface.

###  3.6 Summarizing the effects on $P_{\star}$

In this work, we have studied the influence of various stellar and disk
parameters on the stellar spin evolution during the first million years of
stellar evolution. We want to summarize our results in an illustrative way in
Table 2. We divide the stellar evolution into two phases: the embedded phase
$\lesssim 10^{5}$ year and the later ($\gtrsim 10^{5}$ years) phase, which
corresponds to the Class II or T Tauri phase. For each parameter variation, we
indicate if the star rotates faster (red), slower (blue), or approximately
unchanged (grey) compared to the reference model.

Table 2: Summary of the influence of different stellar and disk parameters on the stellar spin evolution. For each parameter, we distinguish between 3 different effects: the star rotates faster (red), slower (blue), or approximately unchanged (grey) compared to the reference model. | Evolution phase
---|---
| Embedded | Class II
fast $\Omega_{\mathrm{init}}$ |  gray!25unchanged |  gray!25unchanged
strong $B_{\star}$ |  red!25fast |  blue!25slow
high $W$ |  blue!25slow |  blue!25slow
high $\beta_{\mathrm{up}}$ |  blue!25slow |  blue!25slow
low $M_{\star}$ at 1 Myr |  red!25fast |  blue!25slow
high $R_{\mathrm{init}}$ ($\beta_{\mathrm{up}}=0.2$) |  gray!25unchanged |  gray!25unchanged
high $R_{\mathrm{init}}$ ($\beta_{\mathrm{up}}=0.005$) |  blue!25slow |  blue!25slow
high $\dot{M}_{\star}$ |  blue!25slow |  red!25fast
$M_{\mathrm{add}}=0.01$ |  blue!25slow |  gray!25unchanged
disk winds (DW) |  blue!25slow |  blue!25slow
low $Z_{\star}$ |  red!25fast |  red!25fast

##  4 Discussion

###  4.1 Applicability of the spin model for embedded stars

In this study, we use the stellar spin evolution model presented in Ireland et
al. (2021) that depends strongly on best-fit parameters and proportionality
constants (see Eqs. 5 \- 11). These parameters are based on models with
smaller accretion rates $\dot{M}_{\mathrm{acc}}\lesssim
10^{-8}\leavevmode\nobreak\ \mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$.
Such accretion rates are typically found in older, further evolved systems
with ages of $\sim 1$ Myr, which marks the end of our simulations. In their
work, Ireland et al. (2021) recommend the use of their formulation for
$R_{\star}\leq R_{\mathrm{t}}\leq R_{\mathrm{co}}$, which is the case for all
our models. However, the large accretion rates can result in effects not
considered by Ireland et al. (2021).

Disk truncation radius: The disk truncation radius, $R_{\mathrm{t}}$, depends
on the balance between stellar magnetic pressure and ram pressure of the
infalling material (e.g., Hartmann et al., 2016). For accretion rates of
$\dot{M}_{\mathrm{acc}}\lesssim 10^{-8}\leavevmode\nobreak\
\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$, the truncation radius is
usually located at several stellar radii and the dipolar magnetic field
component dominates the magnetic field strength at this position (e.g., Finley
& Matt, 2018; Ireland et al., 2021). During phases of higher accretion rates,
however, the truncation radius is pushed toward the star and can reach
$R_{\mathrm{t}}=R_{\star}$. Close to the stellar surface, the importance of
high-order magnetic field components increases. Quadro- or octopole components
can have field strengths of several kilogauss (e.g., Kochukhov et al., 2020),
which could also affect the position of the truncation radius close to the
star. On the other hand, trends in the evolution of the stellar magnetic field
suggest that young, fully-convective stars have axisymmetric fields dominated
by a dipole component (e.g., Folsom et al., 2016). Thus, we assume that the
calculation of the truncation radius as presented in Ireland et al. (2021)
provides a good approximation for young, fully-convective stars.

Accretion torque: The accretion torque, $\Gamma_{\mathrm{acc}}$, depends on
the ’Keplerianity’ of the disk material at the truncation radius. Angular
momentum removed by MEs or stellar magnetic torques results in sub-Keplerian
rotation of the disk material (e.g., Bessolaz et al., 2008; Matt et al., 2012;
Ireland et al., 2021). In the current formulation, these effects depend on the
proportionality constant $K_{\mathrm{acc}}=0.775$ and the position of the
truncation radius (see Eq. 7). While the formulation should provide reasonable
estimates for the accretion torque for $R_{\mathrm{t}}>R_{\star}$,
$\Gamma_{\mathrm{acc}}$ could be underestimated in case the truncation radius
would be located inside the stellar radius, due to strong accretion rates or,
for example, weak magnetic field strengths. In such cases, we set
$R_{\mathrm{t}}=R_{\star}$. As a result, the ram pressure of the accreting
material exceeds the stellar magnetic pressure, and the accreting material
could be close to Keplerian rotation.

Magnetospheric ejections: The torque due to magnetospheric ejections,
$\Gamma_{\mathrm{ME}}$, depends, similarly to the truncation radius, on the
stellar magnetic field strength and the location of the truncation radius (see
Eq. 9). For truncation radii close to the stellar surface, small-scale
magnetic field components could affect the magnitude of
$\Gamma_{\mathrm{ME}}$. Assuming primarily axisymmetric fields dominated by a
dipole component (e.g., Folsom et al., 2016), the current formulation should
provide reasonable estimates for $\Gamma_{\mathrm{ME}}$.

APSW torque: In the current model, accretion-powered stellar winds (APSW)
remove angular momentum from the star (see Eq. 11 and, e.g., Matt & Pudritz,
2005a; Gallet et al., 2019; Ireland et al., 2021). These studies, however, do
not consider the very early, embedded stages in their formulations. It is
unclear how winds emerging from the stellar surface interact with an envelope
falling onto the disk and how much angular momentum can be transported away
from the star under these circumstances. In case the envelope around the young
star has little or no influence on the APSW, the current formulation with
efficiency rates of $W\sim 2\%$ provides a realistic estimate for the torque
due to APSW. If, on the other hand, the infalling envelope interferes with the
APSW, the efficiency rate could be reduced significantly, which is
demonstrated in this study by $W=0\%$.

In summary, the model presented in Ireland et al. (2021) should provide
reasonable estimates for the torques acting on the young protostar during the
embedded phase with high accretion rates. The individual relations between
stellar and disk parameters are assumed to be valid. We note, however, that
the proportionality constants used in Eqs. 5 \- 11 could change when repeating
the calculations presented in Ireland et al. (2021) with higher accretion
rates.

###  4.2 Magnitude and distribution of $\beta_{\mathrm{up}}$

The evolution of young stars strongly depends on the amount of additional
heat, $\dot{E}_{\mathrm{add}}$, added into the stellar interior (e.g., Baraffe
et al., 2012; Vorobyov et al., 2017a, b; Kunitomo et al., 2017). Usually,
$\dot{E}_{\mathrm{add}}$ is assumed to be a certain fraction, $\beta$, of the
energy of the accreted material. In the hybrid accretion scenario (e.g.,
Vorobyov et al., 2017b), $\beta$ depends on the accretion rate on the star,
ranging from small values, $\beta_{\mathrm{low}}$, during phases of low
accretion rates to high values, $\beta_{\mathrm{up}}$, during phases of high
accretion rates (e.g., Jensen & Haugbølle, 2018; Elbakyan et al., 2019;
Steindl et al., 2021). According to Kunitomo et al. (2017), relatively large
values of $\beta\gtrsim 0.1$ are necessary to explain the evolution of a
majority of stars in young clusters555We note that for a higher Deuterium
abundance lower values of $\beta$ suffice to explain the evolution of a
majority of stars (Kunitomo et al., 2017). Since a variation in the Deuterium
abundance has similar effects on the stellar evolution compared to a variation
in $\beta_{\mathrm{up}}$, we forego the simulation of specific models.. Low
values of $\beta$ (e.g., $\beta=\beta_{\mathrm{low}}$), on the other hand,
result in stellar luminosities of $\lesssim 0.1\leavevmode\nobreak\
\mathrm{L_{\odot}}$ and can explain very low-luminosity objects (VeLLOs, e.g.,
Vorobyov et al., 2017a). To explain the luminosity spread in a young cluster,
a distribution of $\beta_{\mathrm{up}}$ has to be assumed. The exact processes
that constrain the magnitude of $\beta_{\mathrm{up}}$ are still unknown and
more detailed modeling of the accretion process and the distribution of the
additional heat is required to formulate reasonable limits. Our results show
that the amount of additional heat affects the stellar spin evolution
significantly. Low values of $\beta_{\mathrm{up}}$ result in very fast
rotating stars (see Panel d in Fig. 3). It would be interesting to test if
low-luminosity stars tend to rotate faster compared to objects with higher
luminosities. Unfortunately, the available observational data do not allow any
conclusion to be drawn so far.

###  4.3 External influences

The previous results are valid for a single star without any external or
environmental influences. Recently, Roquette et al. (2021) studied the
influence of external photo-evaporation caused by massive stars in young
regions. The presence of a high-mass star can remove mass from the disk due to
high-energy radiation and reduce the disk’s lifetime. The timescale, on which
the star and the disk can exchange angular momentum is reduced and the stars
start to spin up due to contraction at an early age. Thus, such stars show
faster rotation periods after several million years. During the initial
stellar evolution, however, Roquette et al. (2021) assume disk-locking,
meaning that the stellar rotation period is kept constant as long as the
accretion disk is present. The effect of external photo-evaporation on the
stellar spin evolution during the first million years can be summarized as
follows: During the embedded phase the disk material is effectively shielded
by the envelope and external photo-evaporation is assumed to have little
effect on the disk evolution (e.g., Winter et al., 2020). Without the
envelope, on the other hand, the high-energy radiation removes material from
the disk, and the disk accretion rate decreases for $\gtrsim 10^{5}$ years. As
shown in the previous models with reduced accretion rates at $\gtrsim 10^{5}$
years (e.g., model $\dot{M}_{\mathrm{high}}$, see Sec. 4.2), the star spins
down. While assuming disk-locking, this effect is not seen in Roquette et al.
(2021). To model the extent of the spin-down due to lower accretion rates and
the spin-up due to shorter disk lifetimes in a high-energy radiation
environment, a combination of our early spin evolution model and a
hydrodynamic disk model including the effects of external photo-evaporation
will provide valuable insights in the early spin evolution.

Another environmental factor is the tidal interaction between different stars
in close encounters, binary, or multiple systems. Tidal interactions in close
encounters can remove disk material and reduce the disk accretion rate and
lifetime (e.g., Winter et al., 2018). In contrast to external photo-
evaporation, these close encounters can act on the star-disk system during all
phases of its evolution. The amount of mass removed, depends on the mass of
the ’intruding object’ (e.g., Vorobyov et al., 2020) and the number of
encounters (e.g., Winter et al., 2018). After the embedded phase, a reduced
accretion rate results in a stellar spin-down (similar to the effect of
external photo-evaporation). During the embedded phase, a reduced accretion
rate will reduce $\dot{E}_{\mathrm{add}}$ and the stellar inflation. As a
consequence, the star rotates faster (see our model with
$M_{\star}=0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$, Panel e in Fig. 3). On
the other hand, close encounters can ignite outburst events (e.g., Vorobyov et
al., 2021) that result in a stellar spin-down. Which of the two mechanisms
predominates is unclear and requires a combined model.

Binaries or multiple stellar companions can have different effects on the
disk’s lifetime and thus, on the stellar spin evolution. Accretion disks,
truncated by a wide binary companion, are assumed to have shorter lifetimes as
mass is removed from the accretion disk (e.g., Kraus et al., 2011). As a
consequence, the star spins up during the embedded phase and spins down
afterward compared to a single star. On the other hand, circumbinary disks
have longer lifetimes compared to single stars (e.g., Alexander, 2012). In
multiple stellar systems, certain configurations can also result in long disk
lifetimes (e.g., Ronco et al., 2021).

###  4.4 Rotational distributions in young clusters

Our results show the impact of different parameters on stellar spin evolution.
Now, we want to overview the observed rotational distributions of young
clusters and estimate how fast (slow) a $\sim 1$ Myr old star probably spins.
In the next step, we infer the parameters necessary to reproduce the observed
distribution and discuss the physical consequences of these choices (see Sec.
4.5). We have selected the rotational period distribution of six young
($\lesssim 3$ Myr) clusters: NGC 6530 (Henderson & Stassun, 2012), Orion SFR
(Serna et al., 2021), Taurus (Rebull et al., 2020), NGC 2264 (Venuti et al.,
2017), $\sigma$ Ori (Cody & Hillenbrand, 2010), and ONC (Rebull et al., 2006).
In Fig. 6, we show the 5. (left-facing triangle), 50. (diamond), and 95.
(right-facing triangle) percentile for each cluster. For all six clusters, the
95. percentiles indicating the slow-rotating stars are located at $\gtrsim 10$
days with a mean value of 13.0 days (vertical dashed line). The fast-rotating
boundary is at $\lesssim 1$ day with the exception of Taurus, which contains
faster-rotating stars with the 5. percentile at 0.15 days. The mean value of
the 5. percentile is 0.76 days (vertical dotted line). One explanation for the
rapidly rotating stars in Taurus could be its age and the evolutionary stage
of the star-disk systems. With an age of $\sim 3$ Myr, the sample of Taurus
consists of disk-bearing stars and stars without an accretion disk (Rebull et
al., 2020). In addition, the ages of these young clusters are subjected to
large scatters in the order of up to $\sim 1$ Myr (e.g., Testi et al., 2022).
Some stars might be significantly older than the final age of our models (1
Myr). During this additional time, the accretion disk could be dissolved due
to accretion, winds, or photo-evaporative effects and the star starts to spin
up due to contraction. This theory is supported by a comparison of stars with
and without an accretion disk in Rebull et al. (2020). Stars, which are still
surrounded by an accretion disk, tend to rotate slower than stars without a
disk. The fastest, disk-bearing stars have rotation periods slightly below 1
day, which matches the 5. percentiles of the other five clusters in our
selection.

Figure 6:  Comparison between the rotational distribution in young clusters
and rotation period range at 1 Myr of our fast and slow-rotating models (blue
area). The diamonds, left-pointing triangles, and right-pointing triangles are
the 50., 5., and 95. percentile of the rotational distribution for the
respective cluster. The mean values of the 95. (5.) percentiles are indicated
by the vertical dashed (dotted) line. See text for references.

###  4.5 Slow and fast rotating configurations

Based on the rotational distribution of young star-forming regions (see Sec.
4.4), 90% of the stars in young clusters show rotation periods between 0.76
days and 13.0 days. While the rotation periods of our models match the fast-
rotating stars in observed clusters (especially in the case of $\beta=0.005$
and $W=0\%$, see Fig. 3), observed slow-rotating stars have rotation periods
that are higher compared to our model results by a factor of $\gtrsim 4$. We
want to explore possible scenarios that can lead to these rotation periods.

A rising number of observed rotational distributions of young clusters allows
for deriving empirical relations between the stellar rotation period and other
parameters. For accreting T Tauri stars, a positive correlation between the
rotation period and the stellar magnetic field strength has been reported
(e.g., Vidotto et al., 2014). For Class II stars, a strong magnetic field
result in a stellar spin down. This relation agrees with our results (see
Panel b in Fig. 3) and other studies (e.g., Matt et al., 2010b; Ireland et
al., 2021; Gehrig et al., 2022). Since the magnetic field strength alone does
not spin down the star sufficiently to match the observed slow-rotating
stellar periods (see Fig. 3), the question arises if strong magnetic fields
(in the star or the disk) might affect other parameters. Bu & Mosallanezhad
(2018) report that disk winds in a strongly magnetized environment are more
powerful compared to a weakly magnetized outflow. Similarly, Bai et al. (2016)
explored the relation between angular momentum loss in protoplanetary disks
due to the disk winds and the magnetic field strength. The higher the magnetic
field strength, the more effective is the wind in removing angular momentum
from the disk and the higher the accretion rate induced by the wind. A
stronger accretion rate increases the amount of additional heat,
$\dot{E}_{\mathrm{add}}$, which is added to the star. In addition, a strong
magnetic field affects the accretion-powered stellar wind. The Alfvén radius,
$r_{\mathrm{A}}$, which acts as the lever arm of the APSW grows with an
increasing stellar magnetic field strength (e.g., Gallet et al., 2019; Ireland
et al., 2021). The effect of the magnetic field strength on the mass-loss rate
of the APSW, however, is poorly constrained by both, theory and observations.
To summarize, a strong magnetic field strength in the star and the disk can
result in an effective disk wind, APSW, and can increase the accretion rate
onto the star. On the other hand, a weak magnetic field strength would be
accompanied by weak winds (APSW and disk winds).

To test the effects of these configurations on the stellar spin evolution, two
models slow and fast are constructed. Model slow is characterized by a
magnetic field strength of $B_{\star}=2.0$ kG, the inclusion of disk winds
removing angular momentum from the inner disk (see Sec. 3.4), and
$\dot{M}_{\mathrm{high}}$ (see Sec. 3.3). All other values match the reference
model. Model fast combines a weak magnetic field strength, $B_{\star}=0.1$ kG,
with an ineffective wind (no disk winds and $W=0\%$). In Fig. 7, the evolution
of the rotation periods for the model slow (blue line) and fast (yellow line)
are shown in comparison to the reference model (black line). The slow-rotating
model spins down during the initial inflation to 22 days. After spinning up to
$\sim 10$ days during the embedded phase, the strong magnetic field strength
spins down the star, reaching 1 Myr with a period of 12.9 days. The fast-
rotating model spins up initially to a rotation period of 1.7 days due to the
weak APSW wind. This period is kept roughly constant during the embedded
phase. During the Class II phase, the weak magnetic field strength results in
a further spin-up to a rotation period of 0.68 days at an age of 1 Myr. Thus,
the effect of the magnetic fields in the young star-disk system on other
parameters like the effectiveness of winds or the accretion rate can result in
rotation periods that agree with the observed range in young clusters (blue
shaded are in Fig. 7). We note that the influence of magnetic fields on the
amount of additional energy controlled by $\beta$ is unclear. Additionally,
the magnetic field strength itself is a free parameter in our model. The
generation of the magnetic field and its evolution is still unclear. One
possible evolution scenario for the magnetic field strength is given by
Browning et al. (2016). The magnetic energy is scaled with the convective
energy in the stellar interior and a positive correlation between the stellar
mass and the magnetic field strength is proposed. Further investigation is
needed to constrain the evolution of the free parameters used in our model.

Another empirical relation connects the stellar rotation period and the
stellar mass. In Herbst et al. (2007), a higher stellar mass is associated
with longer rotation periods in $\sim 1$ Myr old clusters. We note that in the
context of this work, a high-mass star corresponds to a solar-mass star. On
the other hand, a weak and even negative correlation between the stellar
rotation period and the stellar mass is reported in Henderson & Stassun (2012)
for cluster ages of $\sim 1$ Myr. Our results presented in Fig. 3 agree with
the negative correlation between stellar period and mass. However, the disks
surrounding high-mass stars, which are more massive compared to low-mass stars
(e.g., Vorobyov, 2011; Andrews et al., 2013), could launch stronger disk winds
as these winds scale with the density in the disk (e.g., Bai et al., 2016). In
addition, the stellar mass scales with the accretion rate (e.g., Testi et al.,
2022; Somigliana et al., 2022) and a high-mass star could be injected with a
larger amount of additional heat, $\dot{E}_{\mathrm{add}}$. Thus, a
combination of high stellar masses, strong accretion rates, and an effective
disk wind can also result in a slow-rotating star. To test the relation
between the stellar mass and the rotation period in more detail and constrain
the mentioned parameters, a combined stellar and disk evolution model,
including disk winds, and further observations are required.

Figure 7:  Evolution of the stellar rotation period, $P_{\star}$ for the slow
(blue line), reference (black line), and fast rotating models (yellow line).
For comparison, the range of the mean values of the 5. and 95. percentile is
highlighted as a blue area, $P_{\mathrm{\star\leavevmode\nobreak\ obs}}$ (see
Sec. 4.4).

##  5 Conclusion

In this work, we present the evolution of stellar rotation periods for young
($\leq 1$ Myr), low-mass ($\leq 1\leavevmode\nobreak\ \mathrm{M_{\odot}}$)
stars. We start our simulation from a stellar seed ($0.05\leavevmode\nobreak\
\mathrm{M_{\odot}}$) and take protostellar accretion, stellar winds, the
magnetic star-disk interaction as well as the internal stellar evolution into
account. Our model allows us to test the influence of a variety of parameters
on the stellar rotation period. Based on our results, we can draw the
following conclusions:

  * •

The stellar rotation period is significantly influenced by the amount of
additional heat, $\dot{E}_{\mathrm{add}}$, (controlled by
$\beta_{\mathrm{up}}$) that is injected into the stellar interior, the
accretion history, and the presence of strong winds (APSW and disk winds). The
magnetic field strength, which is dominant during later evolutionary stages,
plays a noticeable but smaller role during the first million years.

  * •

Differences in the initial rotational periods of the stellar seed,
$\Omega_{\mathrm{init}}$, are quickly forgotten ($\lesssim 10^{4}$ years) and
do not impact the rotational distribution at later stellar ages. Similarly,
the initial stellar radius, $R_{\mathrm{init}}$, (assuming
$\beta_{\mathrm{up}}\sim 0.2$) and the extent of the region, in which the
additional heat is distributed, $M_{\mathrm{add}}$, do not affect the rotation
period noticeably.

  * •

We explore the possibility of fast and slow-rotating models. The magnetic
field strength in the star-disk system can affect the amount of additional
heat, $\dot{E}_{\mathrm{add}}$, and the effectiveness of both, stellar and
disk winds. A strong magnetic field can lead to an enhanced accretion rate and
effective winds (APSW and disk winds), resulting in slow-rotating stars. Weak
magnetic fields or a small amount of additional heat ($\beta=0.005$) lead to
the fastest rotating stars. The stellar rotation periods in our models cover a
range from $0.6-12.9$ days.

  * •

Our models match the rotational distribution of six young ($\lesssim 3$ Myr)
clusters. The rotation periods of up to 90% of all observed stars are located
within our results.

  * •

Motivated by these intriguing results, our model, extended over the whole disk
lifetime, can be used as an alternative to the widely used disk-locking model
during the presence of an accretion disk. In combination with a hydrodynamic
disk model, including effects such as external photo-evaporation and episodic
accretion outbursts, valuable insights into the origin of the rotational
distribution of young clusters can be obtained.

###### Acknowledgements.

We thank the anonymous referee for the constructive feedback that helped to
improve the manuscript. The authors gratefully acknowledge the work of Bill
Paxton and his collaborators on the stellar evolution code MESA. Furthermore,
the authors thank Thomas Steindl and Konstanze Zwintz for providing valuable
insights into ”handling” MESA and its parameter files. E.I.V. acknowledges
support by the Ministry of Science and Higher Education of the Russian
Federation (State assignment in the field of scientific activity 2023). The
simulations were performed on the VSC Vienna Scientific Cluster.

## References

  * Affer et al. (2013) Affer, L., Micela, G., Favata, F., Flaccomio, E., & Bouvier, J. 2013, MNRAS, 430, 1433
  * Alexander (2012) Alexander, R. 2012, The Astrophysical Journal Letters, 757, L29
  * Allard et al. (2001) Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357
  * Amard & Matt (2020) Amard, L. & Matt, S. P. 2020, ApJ, 889, 108
  * 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. (2013) Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129
  * Angulo et al. (1999) Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A, 656, 3
  * Armitage & Clarke (1996) Armitage, P. J. & Clarke, C. J. 1996, MNRAS, 280, 458
  * Bai et al. (2016) Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, The Astrophysical Journal, 818, 152
  * 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
  * Bessolaz et al. (2008) Bessolaz, N., Zanni, C., Ferreira, J., Keppens, R., & Bouvier, J. 2008, A&A, 478, 155
  * Bouvier et al. (2014) Bouvier, J., Matt, S. P., Mohanty, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 433
  * Browning et al. (2016) Browning, M. K., Weber, M. A., Chabrier, G., & Massey, A. P. 2016, ApJ, 818, 189
  * Bu & Mosallanezhad (2018) Bu, D.-F. & Mosallanezhad, A. 2018, A&A, 615, A35
  * Buchler & Yueh (1976) Buchler, J. R. & Yueh, W. R. 1976, ApJ, 210, 440
  * Burdonov et al. (2021) Burdonov, K., Bonito, R., Giannini, T., et al. 2021, A&A, 648, A81
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094
  * Castelli & Kurucz (2003) Castelli, F. & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20
  * 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
  * Decampli (1981) Decampli, W. M. 1981, ApJ, 244, 124
  * Elbakyan et al. (2019) Elbakyan, V. G., Vorobyov, E. I., Rab, C., et al. 2019, MNRAS, 484, 146
  * 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
  * Flores et al. (2019) Flores, C., Connelley, M. S., Reipurth, B., & Boogert, A. 2019, ApJ, 882, 75
  * Folsom et al. (2016) Folsom, C. P., Petit, P., Bouvier, J., et al. 2016, Monthly Notices of the Royal Astronomical Society, 457, 580
  * 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
  * Gehrig et al. (2023) Gehrig, L., Steindl, T., Vorobyov, E. I., Guadarrama, R., & Zwintz, K. 2023, A&A, 669, A84
  * Gehrig et al. (2022) Gehrig, L., Steiner, D., Vorobyov, E. I., & Güdel, M. 2022, A&A, 667, A46
  * Ghosh & Lamb (1979) Ghosh, P. & Lamb, F. K. 1979, ApJ, 234, 296
  * Guarcello et al. (2021) Guarcello, M. G., Biazzo, K., Drake, J. J., et al. 2021, A&A, 650, A157
  * Guilet & Ogilvie (2014) Guilet, J. & Ogilvie, G. I. 2014, Mon. Not. R. Astron. Soc. Vol. 441, Issue 1, p.852-868, 441, 852
  * 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
  * Hauschildt et al. (1999a) Hauschildt, P. H., Allard, F., & Baron, E. 1999a, ApJ, 512, 377
  * Hauschildt et al. (1999b) Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., & Alexander, D. R. 1999b, ApJ, 525, 871
  * Henderson & Stassun (2012) Henderson, C. B. & Stassun, K. G. 2012, The Astrophysical Journal, 747, 51
  * Herbst et al. (2002) Herbst, W., Bailer-Jones, C. A. L., Mundt, R., Meisenheimer, K., & Wackermann, R. 2002, A&A, 396, 513
  * Herbst et al. (2007) Herbst, W., Eislöffel, J., Mundt, R., & Scholz, A. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 297
  * Hillenbrand (2009) Hillenbrand, L. A. 2009, in The Ages of Stars, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, Vol. 258, 81–94
  * Hosokawa et al. (2011) Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, The Astrophysical Journal, 738, 140
  * Hsieh et al. (2019) Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, The Astrophysical Journal, 884, 149
  * 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 & Bouvier (2009) Irwin, J. & Bouvier, J. 2009, in The Ages of Stars, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, Vol. 258, 363–374
  * 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
  * Johns-Krull et al. (2009) Johns-Krull, C. M., Greene, T. P., Doppmann, G. W., & Covey, K. R. 2009, ApJ, 700, 1440
  * Johnstone et al. (2014) Johnstone, C. P., Jardine, M., Gregory, S. G., Donati, J. F., & Hussain, G. 2014, MNRAS, 437, 3202
  * Kadam et al. (2021) Kadam, K., Vorobyov, E., & Ágnes Kóspál. 2021, The Astrophysical Journal, 909, 31
  * Kochukhov et al. (2020) Kochukhov, O., Hackman, T., Lehtinen, J. J., & Wehrhahn, A. 2020, A&A, 635, A142
  * Koenigl (1991) Koenigl, A. 1991, ApJ, 370, L39
  * Kraus et al. (2011) Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2011, The Astrophysical Journal, 745, 19
  * Kunitomo et al. (2017) Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A, 599, A49
  * Langanke & Martínez-Pinedo (2000) Langanke, K. & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481
  * Laos et al. (2021) Laos, S., Greene, T. P., Najita, J. R., & Stassun, K. G. 2021, ApJ, 921, 110
  * Maeder (2009) Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars
  * Matt & Pudritz (2005a) Matt, S. & Pudritz, R. E. 2005a, ApJ, 632, L135
  * Matt & Pudritz (2005b) Matt, S. & Pudritz, R. E. 2005b, MNRAS, 356, 167
  * Matt et al. (2010a) Matt, S. P., Pinzón, G., de la Reza, R., & Greene, T. P. 2010a, ApJ, 714, 989
  * Matt et al. (2010b) Matt, S. P., Pinzón, G., de la Reza, R., & Greene, T. P. 2010b, 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
  * Micela et al. (1985) Micela, G., Sciortino, S., Serio, S., et al. 1985, ApJ, 292, 172
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231
  * 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
  * Pillitteri et al. (2019) Pillitteri, I., Sciortino, S., Reale, F., et al. 2019, A&A, 623, A67
  * 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
  * Rebull et al. (2014) Rebull, L. M., Cody, A. M., Covey, K. R., et al. 2014, AJ, 148, 92
  * Rebull et al. (2020) Rebull, L. M., Stauffer, J. R., Cody, A. M., et al. 2020, AJ, 159, 273
  * Rebull et al. (2006) Rebull, L. M., Stauffer, J. R., Megeath, S. T., Hora, J. L., & Hartmann, L. 2006, The Astrophysical Journal, 646, 297
  * Réville et al. (2015) Réville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015, ApJ, 798, 116
  * Rogers & Nayfonov (2002) Rogers, F. J. & Nayfonov, A. 2002, ApJ, 576, 1064
  * Ronco et al. (2021) Ronco, M. P., Guilera, O. M., Cuadra, J., et al. 2021, The Astrophysical Journal, 916, 113
  * 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
  * Serna et al. (2021) Serna, J., Hernandez, J., Kounkel, M., et al. 2021, The Astrophysical Journal, 923, 177
  * Soderblom et al. (2014) Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 219
  * Somigliana et al. (2022) Somigliana, A., Toci, C., Rosotti, G., et al. 2022, MNRAS, 514, 5927
  * Steindl et al. (2021) Steindl, T., Zwintz, K., Barnes, T. G., Müllner, M., & Vorobyov, E. I. 2021, A&A, 654, A36
  * 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
  * Tsukamoto et al. (2022) Tsukamoto, Y., Maury, A., Commerçon, B., et al. 2022, arXiv e-prints, arXiv:2209.13765
  * Venuti et al. (2017) Venuti, L., Bouvier, J., Cody, A. M., et al. 2017, A&A, 599, A23
  * Vidotto et al. (2014) Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361
  * Vorobyov (2011) Vorobyov, E. I. 2011, ApJ, 729, 146
  * Vorobyov & Basu (2015) Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115
  * Vorobyov et al. (2017a) Vorobyov, E. I., Elbakyan, V., Dunham, M. M., & Guedel, M. 2017a, A&A, 600, A36
  * Vorobyov et al. (2017b) Vorobyov, E. I., Elbakyan, V., Hosokawa, T., et al. 2017b, A&A, 605, A77
  * Vorobyov et al. (2021) Vorobyov, E. I., Elbakyan, V. G., Liu, H. B., & Takami, M. 2021, A&A, 647, A44
  * Vorobyov & Pavlyuchenkov (2017) Vorobyov, E. I. & Pavlyuchenkov, Y. N. 2017, AAP, 606, A5
  * Vorobyov et al. (2020) Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2020, A&A, 635, A196
  * Weber & Davis (1967) Weber, E. J. & Davis, Leverett, J. 1967, ApJ, 148, 217
  * White et al. (2007) White, R. J., Greene, T. P., Doppmann, G. W., Covey, K. R., & Hillenbrand, L. A. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 117
  * Winter et al. (2018) Winter, A. J., Clarke, C. J., Rosotti, G., et al. 2018, Monthly Notices of the Royal Astronomical Society, 478, 2700
  * Winter et al. (2020) Winter, A. J., Kruijssen, J. M. D., Chevance, M., Keller, B. W., & Longmore, S. N. 2020, MNRAS, 491, 903
  * Wright et al. (2011) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48
  * 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
  * Zanni & Ferreira (2013) Zanni, C. & Ferreira, J. 2013, A&A, 550, A99
  * Zhang et al. (2017) Zhang, Q., Claus, B., Watson, L., & Moran, J. 2017, The Astrophysical Journal, 837, 53

##  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 rotational-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 Detailed evolution of individual torque components

We show each component affecting the stellar spin evolution: the external
torques due to the accretion process $\Gamma_{\mathrm{acc}}$, APSW
$\Gamma_{\mathrm{W}}$ and the magnetic star-disk interaction
$\Gamma_{\mathrm{ME}}$, the internal contribution $\Gamma_{\mathrm{int}}$ and
the resulting total torque $\Gamma_{\mathrm{tot}}$.

Figure 8:  Evolution of the different external torque components
($\Gamma_{\mathrm{acc}}$, $\Gamma_{\mathrm{me}}$ and $\Gamma_{\mathrm{w}}$;
solid blue, dashed blue, and dotted blue lines, respectively), the internal
contribution to the stellar spin evolution $\Gamma_{\mathrm{int}}$ (blue dash-
dotted line) and the combination of these mechanisms $\Gamma_{\mathrm{tot}}$
(black solid line) for different initial rotation rates,
$\Omega_{\mathrm{init}}$. The value of the reference model is marked in
boldface. The horizontal dashed black line divides the regions, in which the
star spins up and down.

Figure 9:  Same as Fig. 8 for different magnetic field strengths, $B_{\star}$.

Figure 10:  Same as Fig. 8 for different APSW efficiencies,$W$.

Figure 11:  Same as Fig. 8 for different initial stellar radii,
$R_{\mathrm{init}}$.

Figure 12:  Same as Fig. 8 for different stellar masses, $M_{\star}$, at an
age of 1 Myr.

Figure 13:  Same as Fig. 8 for different maximum energy fractions that can be
added to the stellar interior, $\beta_{\mathrm{up}}$.

##  Appendix D Influence of the initial stellar radius for
$\beta_{\mathrm{up}}=0.005$

Evolution of the stellar rotation period and radius for
$\beta_{\mathrm{up}}=0.005$ and different initial stellar radii.

Figure 14:  Evolution of the stellar rotation periods, $P_{\star}$ (Panel a),
and the stellar radius (Panel b) for different initial stellar radii,
$R_{\mathrm{init}}$. For these simulations we choose
$\beta_{\mathrm{up}}=0.005$.

##  Appendix E HRD for individual simulations

Here, we show the influence of each parameter on the HRD.

Figure 15:  Intrinsic stellar luminosity, $L_{\star}$, over effective
temperature, $T_{\mathrm{eff}}$, for each model presented in Sec. 3 and Sec.
4. The reference model is marked in boldface.