# Co-evolution of dust grains and protoplanetary disks II: structure and
evolution of protoplanetary disks; an analytical approach

Yusuke Tsukamoto1
1Graduate Schools of Science and Engineering, Kagoshima University, Kagoshima,
Japan

###### Abstract

In our previous study (Tsukamoto et al., 2023a), we investigated formation and
early evolution of protoplanetary disks with 3D non-ideal magnetohydrodynamics
simulations considering dust growth, and found that the modified equations of
the conventional steady accretion disk model which consider the magnetic
braking, dust growth and ambipolar diffusion reproduce the disk structure
(such as density and vertical magnetic field) obtained from simulations very
well. The disk structure predicted from the disk model is described by
observationally or experimentally constrainable parameters such as the mass of
central star, mass accretion rate, recombination rate, temperature, and
ionization rate. In this paper, as a sequel of the our previous study, we
analytically investigate the structure and evolution of protoplanetary disks
corresponding to Class 0/I young stellar objects using the modified steady
accretion disk model combining an analytical model of envelope accretion. We
estimate that the disk radius is several AU at disk formation epoch and
increases to several 100 AU at the end of the accretion phase. The disk mass
is estimated to be $0.01\thinspace M_{\odot}\lesssim M_{\rm disk}\lesssim
0.1\thinspace M_{\odot}$ for a disk with radius of several 10 AU and mass
accretion rate of $\dot{M}_{\rm disk}\sim 10^{-6}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$. These estimates seems to be consistent with
recent observations. We also found that, with typical disk ionization rates
($\zeta\gtrsim 10^{-19}{\rm~{}s^{-1}}$) and moderate mass accretion rate
($\dot{M}_{\rm disk}\gtrsim 10^{-8}\thinspace M_{\odot}\thinspace{\rm
yr^{-1}}$), magneto-rotational instability is suppressed in the disk because
of low plasma $\beta$ and efficient ambipolar diffusion. We argue that the
radial profile of specific angular momentum (or rotational velocity) at the
disk outer edge should be continuously connected to that of the envelope if
the disk evolves by magnetic braking, and should be discontinuous if the disk
evolves by internal angular momentum transport process such as gravitational
instability or magneto-rotational instability. Future detailed observations of
the specific angular momentum profile around the disk outer edge are important
for understanding the angular momentum transport mechanism of protoplanetary
disks.

###### keywords:

star formation – circum-stellar disk – methods: magnetohydrodynamics –
protoplanetary disk

##  1 Introduction

Predicting the structure and evolution of protoplanetary disks such as surface
density and magnetic field profile with observationally or experimentally
constrainable parameters is a major challenge in theoretical astrophysics.
Toward this end, many studies have been carried out using multidimensional
magnetohydrodynamics (MHD) simulations that take into account the non-ideal
MHD effects (Li et al., 2011; Machida et al., 2011; Tomida et al., 2015;
Tsukamoto et al., 2015b, a; Wurster et al., 2016; Masson et al., 2016;
Tsukamoto et al., 2017a; Zhao et al., 2018; Wurster & Bate, 2019; Tsukamoto et
al., 2020; Zhao et al., 2020; Xu & Kunz, 2021). The consensus from these
studies is that non-ideal MHD effects, especially ambipolar diffusion, play a
crucial role in disk evolution. Because the strength of non-ideal MHD effects,
magnetic resistivity strongly depends on the microscopic physics (such as the
cosmic-ray ionization, the gas-phase recombination, and the adsorption rate of
charged particles on the dust grains), various studies have also been carried
out to investigate the impact of the microscopic physics on the disk evolution
(Zhao et al., 2016; Wurster et al., 2018; Tsukamoto et al., 2020; Kuffmeier et
al., 2020; Marchand et al., 2021, 2023; Kobayashi et al., 2023).

While it is wonderful that such multidimensional, microphysics-aware
simulations are now possible and thereby reveal the formation and evolution of
protoplanetary disks, the high computational costs of such simulations make it
difficult to study the long-term evolution of protoplanetary disks over $\sim
10^{6}$ yr after protostar formation, through the Class 0 and Class I
evolutionary stages, to the Class II stage. Furthermore, their microscopic and
macroscopic complexities make it seem extremely difficult to achieve the
ultimate goal described above, i.e., to concisely describe the disk evolution
and disk structure in terms of a few parameters.

Recently, Tsukamoto et al. (2023a) (here after TMI23) investigated the
formation and evolution of protoplanetary disks using dust-gas two-fluid non-
ideal MHD simulations with dust growth. The simulations also take into account
the changes of magnetic resistivity associated with dust growth. This study
proposed the concept of a ”co-evolution of dust grains and protoplanetary
disks”, a process in which the global evolution of the disk is significantly
affected by dust growth through the changes in the adsorption efficiency of
charged particles onto the dust grains.

TMI23 also found that when the dust grains becomes sufficiently large (broadly
speaking, the maximum dust size becomes $100~{}{\rm\mu}{\rm m}$ to $1~{}{\rm
mm}$; see Tsukamoto & Okuzumi (2022) for more detailed calculations) and the
adsorption of ions and electrons on the dust grains becomes negligible, the
radial profiles of the disk such as density, magnetic field, radial velocity,
and ambipolar resistivity are well described by the power laws. TMI23 also
showed that a simple extension of the steady accretion disk model reproduces
the disk structure surprisingly well not only the power law exponent but also
the absolute value.

The disk structure of TMI23 is described in terms of observationally
constrainable (or experimentally measurable, such as the gas phase
recombination rate) parameters, and is conducive to the ultimate goal
described above. Furthermore, because it is written analytically, it can be an
important tool for understanding long term evolution of the disks until the
end of Class I. Therefore, it is important to clarify the predictions obtained
from the model and give suggestions for future observations.

However, the disk structure derived in TMI23 is not suitable for the purpose
of comparison with observations because the model employed (simulation-based)
assumptions such as barotropic equation of state to allow the detailed
comparisons with the simulation results. Therefore, in this paper, we derive a
disk structure suitable for comparison to observations by making simple
assumptions on the temperature structure and gas-phase recombination rate for
the disk model proposed in TMI23. Then we discuss its properties and, by
combining an analytical model of envelope accretion and the disk model, we
make observationally testable predictions. We think that some of the
nontrivial predictions for disk mass, radius, and specific angular momentum
profile presented in this paper can be verified by future observations.

##  2 Governing equations and the disk structure

###  2.1 Governing equations of disk model

In TMI23, we suggested that the disk structures obtained by non-ideal MHD
simulations are well described by a steady accretion disk model of,

$\displaystyle-2\pi rv_{r}\Sigma=\dot{M}_{\rm disk},$ (1) $\displaystyle
v_{\phi}=\sqrt{\frac{GM_{\rm star}}{r}}\equiv r\Omega,$ (2) $\displaystyle
H=\frac{c_{s}}{\Omega},$ (3)
$\displaystyle\Sigma(r\Omega)v_{r}=-r\frac{B_{z}B_{\phi,s}}{\pi},$ (4)
$\displaystyle B_{z}v_{r}=-\frac{\eta_{A}}{r}B_{z},$ (5)
$\displaystyle\eta_{A}\frac{B_{\phi,s}}{H}=\left(\frac{H}{r}\right)^{2}B_{z}v_{\phi}.$
(6)

where $r$ is the radius from the central protostar, $\Sigma$ is the surface
density of the disk, $v_{r}$ is the radial velocity of the gas, $\dot{M}_{\rm
disk}$ is the mass accretion rate in the disk which is assumed to be constant,
$v_{\phi}$ is the azimuthal velocity of the gas, $\Omega$ is the angular
velocity. $M_{\rm star}$ is the protostar mass, $c_{s}$ is the sound speed,
$H$ is the scale height, $B_{z}$ is the vertical magnetic field, $B_{\phi,s}$
is the azimuthal magnetic field at $z=H$, $\eta_{A}$ is the ambipolar
resistivity. Among these quantities, $\dot{M}_{\rm disk}$ and $M_{\rm star}$
are constants and do not depend on the radius and other quantities are
functions of the radius. $G$ is the gravitational constant.

The equation (1) to equation (3) are the same as the standard viscous
accretion disk model and our extension is expressed in (4) to (6). The
equation (4) shows that the magnetic braking balances the radial angular
momentum advection. The equation (5) shows that the radial magnetic flux
advection balances the radial magnetic field drift by the ambipolar diffusion.
The equation (6) shows that the azimuthal magnetic field generation by the
vertical shear motion balances the azimuthal magnetic field drift by the
ambipolar diffusion.

Strictly speaking, TMI23 showed that the solutions of equations (A29) to (A33)
in TMI23 which are derived from perturbation theory reproduce the simulation
results well. However, we have confirmed that the disk structures predicted by
the above equations and those predicted by the equations used in TMI23 are
essentially the same in absolute values as well as in power exponents (for
example, the differences in values of $\rho$ and $B_{z}$ are about 20%, and
the all power exponents of the physical quantities are completely the same).
So, we use the equations above in this paper because they are simpler and more
intuitive.

To obtain the solution of the equations (1) to (6), we need a model of
$\eta_{A}$ and the temperature profile. We showed in Tsukamoto & Okuzumi
(2022) that, as dust growth proceeds, $\eta_{A}$ tends to decrease and
converges to the analytic power law of

$\displaystyle\eta_{A}$ $\displaystyle=\frac{\mathbf{B}^{2}}{4\pi
C\gamma\rho^{3/2}}\sim\frac{B_{z}^{2}}{4\pi C\gamma\rho_{\rm mid}^{3/2}},$ (7)

even in the disk, which is determined by the balance between gas phase
ionization and gas phase recombination (Shu, 1983). Here, $\rho_{\rm mid}$ is
the midplane density of the disk. To derive the right-hand side expression, we
assume that the vertical magnetic field is dominant and the contribution of
the radial and azimuthal magnetic field on the field strength are negligible.
We also assume that the density variation in the vertical direction is
negligible.

Note that, using this formula for $\eta_{A}$, we model that the dust grains
have sufficiently grown in the disk. Broadly speaking, this formula is
justified if the maximum dust size becomes $a_{\rm max}\gtrsim 100{\rm\mu}{\rm
m}$ in the disk (see, Tsukamoto & Okuzumi, 2022; Lebreuilly et al., 2023, for
more detail). The dust growth to several hundred microns in size can occur
even in young disks. The growth timescale of dust grains in the disks is
estimated to be the order of $10^{3}$ yr and is much shorter than the age of
Class 0/I young stellar objects (YSOs) (see Tsukamoto & Okuzumi, 2022, for
more detailed estimate). Recent simulations actually show that dust growth and
associated resistivity change can proceed in very young disks (or even in the
first core; Tsukamoto et al., 2021; Marchand et al., 2023; Tsukamoto et al.,
2023a). In addition, signs of dust growth to ${\rm mm}$ in size even in very
young disks have been revealed by recent high-resolution observations (Harsono
et al., 2018; Carrasco-González et al., 2019). Therefore, there are
theoretical and observational motivations for using equation (7).

$C$ is given as

$\displaystyle C=\sqrt{\frac{m_{i}^{2}\zeta}{m_{g}\beta_{r}}},$ (8)

where $m_{i}$ and $m_{g}$ are the mass of ion and neutral particles. In this
paper, We assume that the major ion is HCO+ and $m_{i}=29m_{p}$ where $m_{p}$
is the proton mass. We also assume $m_{g}=2.34m_{p}$. $\zeta$ is the
ionization rate.

$\beta_{r}$ is the recombination rate and assumed to be

$\displaystyle\beta_{r}=\beta_{r,0}\left(\frac{T}{300~{}{\rm K}}\right)^{-1}$
(9)

where $\beta_{r,0}=1.1\times 10^{-7}~{}{\rm cm^{3}s^{-1}}$ which is the rate
of HCO+ and taken from UMIST database (McElroy et al., 2013). Note that we
adopted older value from RATE99 by Smith & Adams (1984) to make power law
simpler (latest recombination rate (RATE12) of UMIST gives $\beta_{r}\propto
T^{-0.69}$).

$\gamma$ is given as

$\displaystyle\gamma=\frac{\langle\sigma v\rangle_{in}}{(m_{g}+m_{i})},$ (10)

where $\langle\sigma v\rangle_{in}$ is the rate coefficient for collisional
momentum transfer between ions and neutrals. We assume $\langle\sigma
v\rangle_{in}=1.3\times 10^{-9}~{}{\rm cm^{3}s^{-1}}$ which is calculated from
the Langevin rate (Pinto & Galli, 2008).

For temperature, we assume the simple power law,

$\displaystyle T(r)=T_{0}\left(\frac{r}{\rm AU}\right)^{-1/2},$ (11)

and vertically isothermal just for simplicity. However, it is straightforward
to use a more complex temperature structures.

We assume the sound velocity to be,

$\displaystyle c_{\rm s}=1.9\times 10^{4}\left(\frac{T}{10{\rm
K}}\right)^{1/2}~{}{\rm cm}~{}{\rm s}^{-1}.$ (12)

###  2.2 Assumptions of the model

In our governing equations, we (implicitly) assume that the internal angular
momentum transport mechanisms such as magneto-rotational instability (MRI),
gravitational instability (GI) are negligible, and magnetic braking is the
mechanism which determines angular momentum evolution of the disk. The
justification for ignoring MRI and GI will be discussed after the actual disk
structure is obtained.

Another important assumption is that the toroidal electric current is
determined by radial gradient of $B_{z}$ as

$\displaystyle J_{\phi}$
$\displaystyle=\frac{4\pi}{c}(\nabla\times\mathbf{B})_{\phi}=\frac{4\pi}{c}(\frac{\partial
B_{r}}{\partial z}-\frac{\partial B_{z}}{\partial r})$
$\displaystyle\sim-\frac{4\pi}{c}\frac{\partial B_{z}}{\partial r}.$ (13)

This assumption is adopted in e.g., Basu & Mouschovias (1994); Dapp & Basu
(2010); Dapp et al. (2012) and leads to right hand side term of equation (5).

On the other hand, another possibility is

$\displaystyle J_{\phi}$ $\displaystyle\sim\frac{4\pi}{c}\frac{\partial
B_{r}}{\partial z}.$ (14)

If we further assume $B_{r}\sim B_{z}$ at $z=H$, the equation (5) is replaced
as

$\displaystyle B_{z}v_{r}=-\frac{\eta_{A}}{H}B_{z}.$ (15)

We also investigated the solutions with this equation and found a few times
weaker vertical magnetic field at $10$ AU and slightly shallower profile of
$B_{z}\propto r^{-29/24}$ because of the stronger radial diffusion for given
magnetic field. On the other hand, the density and radial velocity do not
depend on the choice of the azimuthal electric current. We do not use the
equation (15) in this paper because the resulting magnetic field strength does
not reproduce the simulation results in TMI23.

We assume that radial velocity for the magnetic flux advection ($v_{r}$ in the
left hand side term of equation (5)) equals to the (vertically density
weighted) radial gas velocity. We acknowledge that some previous studies show
that the radial velocity at several scale height is much faster than the
radial velocity at the midplane (so called layered accretion; e.g., Gammie,
1996; Fleming & Stone, 2003; Takasao et al., 2018; Zhu & Stone, 2018). And
there are the discussions that $v_{r}$ in the equation (5) (magnetic field
transport velocity) should be different from gas radial velocity (Guilet &
Ogilvie, 2012, 2014). On the other hand, in the previous simulation studies
started from cloud cores (Machida, 2014; Machida & Basu, 2019), we do not
observe such a significant inward accretion at several scale height rather the
gas is significantly outflowing from several scale height. We speculate that
the difference of flow pattern at several scale height is due to the
significant difference of the magnetic field structure/strength in and around
the disk among the studies, and/or whether non-ideal MHD effects are
considered. For example, the layered accretion becomes less pronounced or
disappears when ambipolar diffusion are taken into account (Gressel et al.,
2015; Bai, 2017).

The causes of the difference of flow pattern should be investigated in future
studies. At the same time, we think the validity of the assumptions made in a
model should ultimately be verified observationally. In this study, we provide
the quantitative predictions of the radial profiles, mass, radius of the disk
predicted from our model. This will help to validate our disk model by future
observational studies.

###  2.3 Disk profiles

By solving the equations (1) to (6) with the model for $\eta_{A}$ (equation
(7)) and temperature (equation (11)) assuming the physical quantities have
power law $f(r)\propto r^{q}$, we can derive the physical quantities of disk.

From equation (5) and (6), $v_{r}$ and $B_{\phi,s}$ can be written as

$\displaystyle v_{r}$ $\displaystyle=-\frac{\eta_{A}}{r},$ $\displaystyle
B_{\phi,s}$
$\displaystyle=\left(\frac{H^{3}}{\eta_{A}r^{2}}\right)B_{z}v_{\phi}.$

By substituting these expressions into equations (1) and (4), we obtain

$\displaystyle 2\pi\eta_{A}\Sigma$ $\displaystyle=\dot{M}_{\rm disk},$
$\displaystyle H^{3}B_{z}^{2}$ $\displaystyle=\pi\eta_{A}^{2}\Sigma.$

Using equations (2), (3), (7), (11), and $\Sigma=\sqrt{2\pi}H\rho_{\rm mid}$
these equations are the simultaneous equations for midplane density and
vertical magnetic field. By solving the simultaneous equations and using the
solutions for $\rho_{\rm mid}$ and $B_{z}$, the disk radial profiles are given
as,

$\displaystyle\rho_{\rm mid}(r)$ $\displaystyle=1.6\times 10^{-12}$
$\displaystyle\left(\frac{\dot{M}_{\rm disk}}{10^{-6}~{}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{\frac{2}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{-\frac{4}{3}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace M_{\odot}}\right)\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{-\frac{1}{3}}$ $\displaystyle\left(\frac{r}{10~{}{\rm
AU}}\right)^{-\frac{7}{3}}~{}{\rm g~{}cm}^{-3},$ (17) $\displaystyle B_{z}(r)$
$\displaystyle=6.2\times 10^{-2}$ $\displaystyle\left(\frac{\dot{M}_{\rm
disk}}{10^{-6}~{}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{\frac{2}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{-\frac{1}{3}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{\frac{1}{6}}$ $\displaystyle\left(\frac{r}{10~{}{\rm
AU}}\right)^{-\frac{4}{3}}{\rm G},$ (18) $\displaystyle v_{r}(r)$
$\displaystyle=-1.4\times 10^{3}$ $\displaystyle\left(\frac{\dot{M}_{\rm
disk}}{10^{-6}~{}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{\frac{1}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{\frac{5}{6}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace
M_{\odot}}\right)^{-\frac{1}{2}}\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{\frac{1}{3}}$ $\displaystyle\left(\frac{r}{10~{}{\rm
AU}}\right)^{\frac{1}{12}}~{}{\rm cm}~{}{\rm s}^{-1},$ (19) $\displaystyle
B_{\phi}(r,z)$ $\displaystyle=-1.2\times 10^{-2}\left(\frac{z}{H}\right)$
$\displaystyle\left(\frac{\dot{M}_{\rm disk}}{10^{-6}~{}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{\frac{1}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{\frac{1}{3}}$ $\displaystyle\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{-\frac{1}{6}}$ $\displaystyle\left(\frac{r}{10~{}{\rm
AU}}\right)^{-\frac{7}{6}}{\rm G},$ (20) $\displaystyle\eta_{A}(r)$
$\displaystyle=2.1\times 10^{17}$ $\displaystyle\left(\frac{\dot{M}_{\rm
disk}}{10^{-6}~{}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{\frac{1}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{\frac{5}{6}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace
M_{\odot}}\right)^{-\frac{1}{2}}\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{\frac{1}{3}}$ $\displaystyle\left(\frac{r}{10~{}{\rm
AU}}\right)^{\frac{13}{12}}{\rm cm^{2}s^{-1}},$ (21)

where we assume that $B_{\phi}$ is linear function of $z/H$.

Then, the surface density is given as follows,

$\displaystyle\Sigma(r)$ $\displaystyle=\sqrt{2\pi}H\rho_{\rm mid}=4.8\times
10^{1}$ $\displaystyle\left(\frac{\dot{M}_{\rm disk}}{10^{-6}~{}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{\frac{2}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{-\frac{5}{6}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{-\frac{1}{3}}$ $\displaystyle\left(\frac{r}{10~{}{\rm
AU}}\right)^{-\frac{13}{12}}~{}{\rm g~{}cm}^{-2}.$ (22)

Plasma $\beta$ and ”Elsasser number” for ambipolar diffusion Am at midplane
are given as

$\displaystyle\beta(r)$ $\displaystyle=\frac{P_{\rm gas}}{P_{\rm
mag}}=3.0\times 10^{1}$ $\displaystyle\left(\frac{\dot{M}_{\rm
disk}}{10^{-6}~{}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{-\frac{2}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{\frac{1}{3}}$ $\displaystyle\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{-\frac{2}{3}}\left(\frac{r}{10~{}{\rm
AU}}\right)^{-\frac{1}{6}},$ (23)

and

$\displaystyle{\rm Am}(r)$
$\displaystyle=\frac{v_{A}^{2}}{\eta_{A}\Omega}=2.7\times 10^{-1}$
$\displaystyle\left(\frac{\dot{M}_{\rm disk}}{10^{-6}~{}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{\frac{1}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{-\frac{1}{6}}$ $\displaystyle\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{\frac{1}{3}}\left(\frac{r}{10~{}{\rm
AU}}\right)^{\frac{1}{12}}.$ (24)

Here, $P_{\rm gas}=\gamma\rho_{\rm mid}c_{s}^{2}$ and $P_{\rm
mag}=B_{z}^{2}/(8\pi)$ are the gas and magnetic pressure, respectively (where
$\gamma=5/3$ is specific heat ratio). $v_{A}=B_{z}/\sqrt{4\pi\rho_{\rm mid}}$
is the Alfvén velocity.

Plasma $\beta$ and Am have a simple relation of,

$\displaystyle\beta({\rm Am})$ $\displaystyle=2.1{\rm Am}^{-2}.$ (25)

Figure 1 shows the radial profile of midplane density (equation (2.3)) and
surface density of the model (equation (2.3)). It shows that the surface
density is more or less consistent with the the Minimum Mass Solar Nebula
(MMSN) model (Hayashi, 1981) when $\dot{M}_{\rm disk}\sim 10^{-6}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$. We also plot the surface density of
gravitationally unstable disk $\Sigma_{\rm GI}$ in the figure calculated from
the condition that Toomre’s $Q$ value is $Q=1$ and the temperature profile of
equation (11),

$\displaystyle\Sigma_{\rm GI}(r)$ $\displaystyle=\frac{c_{s}\Omega}{\pi
GQ}=6.8\times 10^{2}$ $\displaystyle\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{\frac{1}{2}}\left(\frac{M_{\rm star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}$ $\displaystyle\left(\frac{r}{10~{}{\rm
AU}}\right)^{-\frac{7}{4}}~{}{\rm g~{}cm}^{-2}.$ (26)

With large mass accretion rate of $\dot{M}_{\rm disk}\sim 10^{-5}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$ which is almost upper limit of mass
accretion rate for low-mass protostar formation, the surface density becomes
$\Sigma_{\rm GI}<\Sigma$ in $r\gtrsim 100$ AU, and gravitational instability
may develop there even with magnetic braking. On the other hand, once mass
accretion rate becomes $\dot{M}_{\rm disk}\lesssim 10^{-6}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$, the surface density satisfies $\Sigma_{\rm
GI}>\Sigma$ in r$<100$ AU. Although we will neglect gravitational instability
for our disk model in the rest sections of this paper, it is important to
recognize that gravitational instability can develop if the disk becomes
larger than $\sim 100$ AU at a stage where the mass accretion rate is
$\dot{M}_{\rm disk}\sim 10^{-5}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$.

Figure 1:  Radial profiles of midplane density and surface density. The
parameters of $M_{\rm star}=0.3\thinspace M_{\odot}$, $T_{0}=150$ K,
$\zeta=10^{-18}{\rm~{}s^{-1}}$ are used. Long dashed (magenta), solid (black),
and short dashed (red) lines show the profile with $\dot{M}_{\rm
disk}=10^{-5},10^{-6},10^{-7}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$,
respectively. dotted-dashed (green) and dotted (orange) lines show the surface
density of gravitationally unstable disk and the Minimum Mass Solar Nebula
(MMSN) model ($\Sigma_{\rm MMSN}=1.7\times 10^{3}(r/{\rm AU})~{}{\rm
g~{}cm}^{-2}$) (Hayashi, 1981), respectively.

Figure 2 shows the radial profiles of vertical magnetic field (equation (2.3))
and plasma $\beta$ (equation (2.3)). Our model predict that, as the protostar
evolves and the mass accretion rate decreases, the vertical magnetic field
decreases as $\propto\dot{M}^{2/3}_{\rm disk}$ and midplane plasma $\beta$
increases as $\propto\dot{M}^{-2/3}_{\rm disk}$. With the $\dot{M}_{\rm
disk}\sim 10^{-6}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$ which is the
typical value of the Class 0/I YSOs, the vertical magnetic field and $\beta$
at $10$ AU are $B_{z}\sim 60~{}{\rm mG}$ and $\beta\sim 30$, respectively.
This indicates that the disk needs this level of the magnetic field strength
to cause the sufficient mass accretion of $\dot{M}_{\rm disk}\sim
10^{-6}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$ only by magnetic braking
without resorting to other mass accretion mechanisms such as gravitational
instability. Note that the magnetic field strength are consistent with the
disk formation simulations starting from collapsing cloud cores (e.g., Masson
et al., 2016).

On the other hand, with the $\dot{M}_{\rm disk}\sim 10^{-8}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$ which may be suitable for late Class I YSOs
to Class II objects, the vertical magnetic field and $\beta$ at $10$ AU are
$B_{z}\sim 3~{}{\rm mG}$ and $\beta\sim 650$, respectively. This magnetic
field profile may be suitable for the initial conditions of simulations of
isolated star-disk systems just after the envelope is depleted.

Note that our disk model predicts a stronger magnetic field than the typical
initial magnetic field strength of the simulations that examine the impact of
magnetic field in isolated star-disk systems (typically $\beta>10^{3}$; Suzuki
& Inutsuka, 2014; Gressel et al., 2015; Bai, 2017; Zhu & Stone, 2018). By
examining the new parameter space predicted from our disk model with future
high-resolution simulations, a new disk evolution picture is expected to
emerge.

Figure 2:  Radial profiles of vertical magnetic field and plasma $\beta$. The
parameters of $M_{\rm star}=0.3\thinspace M_{\odot}$, $T_{0}=150$ K,
$\zeta=10^{-18}{\rm~{}s^{-1}}$ are used. Long dashed (magenta), solid (black),
short dashed (red), and dotted (blue) lines show the profile with
$\dot{M}_{\rm disk}=10^{-5},10^{-6},10^{-7},10^{-8}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$, respectively.  Figure 3:  Plasma $\beta$ of
the disk as a function of Am. The parameters of $M_{\rm star}=0.3\thinspace
M_{\odot}$, $T_{0}=150$ K, $\zeta=10^{-18}{\rm~{}s^{-1}}$ are used. black,
red, blue, and cyan lines show the profile of the disk with $\dot{M}_{\rm
disk}=10^{-6},10^{-7},10^{-8},10^{-9}$ and $10^{-10}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$, respectively. Upper left end of each line
corresponds to $r=1$ AU and Lower right end of each line corresponds to
$r=100$ AU. The red (or black) point on each line corresponds to $r=10$ AU.
The dashed line shows the boundary between whether MRI operates or not, and in
the region below this line, hatched in light red, MRI does not operate (Bai &
Stone, 2011).

Bai & Stone (2011) performed three-dimensional shearing-box simulations over a
wide range of parameter space, and revealed the parameter space on
$\beta$-${\rm Am}$ plane in which magnetorotational instability (MRI)
operates, corresponding to the requirement that the most unstable vertical
wavelength of MRI should be less than the disk scale height. With the aid of
their conditions, we can discuss whether MRI develops in our disk model, and
whether it is reasonable to neglect MRI.

Figure 3 shows plasma $\beta$ (equation (2.3)) of the disk as a function of Am
(equation (2.3)). In this plot, we consider the disk radius of $1~{}{\rm
AU}<r<100~{}{\rm AU}$. We derived radius of given Am from equation (2.3) and
calculate the plasma $\beta$ at the obtained radius. Upper left end of each
line corresponds to $r=1$ AU, the symbol on the each line corresponds to
$r=10$ AU, and lower right end of each line corresponds to $r=100$ AU. The
different lines represent the different mass accretion rate from
$10^{-6}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$ to $10^{-10}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$. The figure shows that with our fiducial
parameters ($M_{\rm star}=0.3\thinspace M_{\odot}$, $T_{0}=150$ K,
$\zeta=10^{-18}{\rm~{}s^{-1}}$), MRI becomes active in $r\lesssim 100$ AU when
the mass accretion rate satisfies $\dot{M}_{\rm disk}\lesssim
10^{-10}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$. In this paper, we are
interested in Class 0/I YSOs and the mass accretion rate in this phase is
expected to be larger than $10^{-9}\thinspace M_{\odot}\thinspace{\rm
yr^{-1}}$ and neglecting MRI in our disk model is justified.

This figure indicates that MRI is difficult to develop during the early
evolutionary stage of the disk such as those with $\dot{M}_{\rm disk}\gtrsim
10^{-9}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$. For MRI to develop, Am
and $\beta$ must be simultaneously large. However, when Am is large, the
outward transport of magnetic flux by ambipolar diffusion becomes inefficient
and the vertical magnetic field of the disk increases due to the radial
advection, and hence $\beta$ decreases. On the other hand, when ambipolar
diffusion works efficiently, the vertical magnetic field is extracted by
outward diffusion and $\beta$ increases. Because of this self-regulative
behavior of radial advection and outward diffusion of vertical magnetic field,
the disk evolves towards the upper left direction on the $\beta$-Am plane as
the mass accretion rate decreases.

The possible scenario in which MRI develops even with larger $\dot{M}_{\rm
disk}$ is the case that cosmic ray are strongly shielded and the ionization
rate in the disk is lowered. Since Am and $\beta$ depend on the mass accretion
rate and the ionization rate as ${\rm Am}\propto(\dot{M}_{\rm
disk}\zeta)^{1/3}$ and $\beta\propto(\dot{M}_{\rm disk}\zeta)^{-2/3}$,
respectively, decreasing the ionization rate by one order of magnitude and
decreasing the mass accretion rate by one order of magnitude have the same
impact for Am and $\beta$. Thus, for example, if $\zeta$ decreases to
$10^{-22}{\rm~{}s^{-1}}$ in a disk due to low cosmic-ray ionization rate and
lack of short-lived radionuclides, MRI can develop in the disk even when
$\dot{M}_{\rm disk}\sim 10^{-6}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$
because the orange line can also be regarded as a disk with $\dot{M}_{\rm
disk}=10^{-6}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$ and
$\zeta=10^{-22}{\rm~{}s^{-1}}$. In this case, outer disk would be
gravitationally unstable because $\Sigma\propto\zeta^{-1/3}$ and it may become
$\Sigma\gtrsim\Sigma_{\rm GI}$ in $r\gtrsim$ 10 AU (see figure 1) The another
possible scenario in which MRI develops is the case the electric current is
given as $J_{\phi}\sim B_{z}/H$ (as discussed in §2.2). In this case, the
magnetic field becomes a few times smaller and plasma $\beta$ increases about
one order of magnitude, and disk can enter the MRI permitted region.

###  2.4 Mass, magnetic flux, and mass-to-flux ratio of disk

By radially integrating $\Sigma(r)$ and $B_{z}(r)$, we can obtain the disk
mass, disk magnetic flux and mass-to-flux ratio of the protostar-disk system.

$\displaystyle M_{\rm disk}$ $\displaystyle=\int_{0}^{r_{\rm disk}}2\pi
r\Sigma(r)dr=3.0\times 10^{-2}$ $\displaystyle\left(\frac{\dot{M}_{\rm
disk}}{10^{-6}~{}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{\frac{2}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{-\frac{5}{6}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{-\frac{1}{3}}$ $\displaystyle\left(\frac{r_{\rm
disk}}{100~{}{\rm AU}}\right)^{\frac{11}{12}}\thinspace M_{\odot},$ (27)

where $r_{\rm disk}$ is the disk radius. Note that the radial power index of
$\Sigma$ is larger than $-2$ and lower limit of integral does not contribute
to the disk mass and we set it to be zero.

the magnetic flux of the disk is calculated as

$\displaystyle\Phi_{\rm disk}$ $\displaystyle=\int_{0}^{r_{\rm disk}}2\pi
rB_{z}(r)dr=2.7\times 10^{2}$ $\displaystyle\left(\frac{\dot{M}_{\rm
disk}}{10^{-6}~{}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{\frac{2}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{-\frac{1}{3}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{\frac{1}{6}}$ $\displaystyle\left(\frac{r_{\rm
disk}}{100~{}{\rm AU}}\right)^{\frac{2}{3}}{\rm G~{}AU^{2}}.$ (28)

By assuming that the protostar has the most of the mass and the disk has the
most of the magnetic flux of the system, the mass-to-flux ratio of the system
can be estimated as,

$\displaystyle\mu$ $\displaystyle\equiv\left(\frac{M_{\rm star}}{\Phi_{\rm
disk}}\right)/\left(\frac{M}{\Phi}\right)_{\rm crit}=2.0\times 10^{1}$
$\displaystyle\left(\frac{\dot{M}_{\rm disk}}{10^{-6}~{}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{-\frac{2}{3}}\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{\frac{1}{3}}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}\left(\frac{\zeta}{10^{-18}~{}{\rm
s^{-1}}}\right)^{-\frac{1}{6}}\left(\frac{r_{\rm disk}}{100~{}{\rm
AU}}\right)^{-\frac{2}{3}},$ (29)

where we assume $\left(\frac{M}{\Phi}\right)_{\rm
crit}=\frac{0.53}{3\pi}\left(\frac{5}{G}\right)^{1/2}$ (Mouschovias & Spitzer,
1976).

Figure 4:  Disk mass and mass-to-flux ratio as a function of mass accretion
rate. The parameters of $M_{\rm star}=0.3\thinspace M_{\odot}$, $T_{0}=150$ K,
$\zeta=10^{-18}{\rm~{}s^{-1}}$ are used. Solid (black), dashed (red), and
dotted (magenta) lines show the profile with $r_{\rm disk}=100,30,10$ AU,
respectively. The gray horizontal lines in the disk mass show the mass of GI
disks (equation (2.4)).

Figure 4 shows the disk mass and mass-to-flux ratio $\mu$ as a function of
mass accretion rate. The disk mass is $0.01\thinspace M_{\odot}\lesssim
M\lesssim 0.1\thinspace M_{\odot}$ for mass accretion rates $\dot{M}_{\rm
disk}\gtrsim 10^{-6}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$ and several
10 AU sized disks. This value is consistent to the disk mass suggested from
observations (Tobin et al., 2020; Tychoniec et al., 2020). In the figure, we
also plot the mass of gravitationally unstable disk,

$\displaystyle M_{\rm GI,disk}$ $\displaystyle=\int_{0}^{r_{\rm disk}}2\pi
r\Sigma_{\rm GI}(r)dr=3.4\times 10^{-1}$
$\displaystyle\left(\frac{T_{0}}{150~{}{\rm
K}}\right)^{\frac{1}{2}}\left(\frac{M_{\rm star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}$ $\displaystyle\left(\frac{r_{\rm
disk}}{100~{}{\rm AU}}\right)^{\frac{1}{4}}\thinspace M_{\odot}.$ (30)

The mass of gravitationally unstable disk is $M_{\rm GI,disk}\gtrsim
0.1\thinspace M_{\odot}$. This value seems not to be consistent to the disk
mass estimated from observations although the disk mass estimate from the
observations also has uncertainties such as the uncertainty of dust opacity or
the optical thickness (see Tsukamoto et al., 2017b, 2023b).

It is worth noting that the mass of gravitationally unstable disk depends
weakly on the radius. Since the disk size has a large dispersion (Maury et
al., 2019), it may be possible to determine whether the disk is
gravitationally unstable or not by investigating whether the disk mass has a
correspondingly large dispersion.

The mass-to-flux ratio is $\mu\gtrsim 10$ unless high mass accretion rate
($\dot{M}_{\rm disk}\sim 10^{-5}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$)
and large disk $r\sim 100$ AU are realized simultaneously. This value is
sufficiently larger than the mass-to-flux ratio of the molecular cloud core
($\mu\sim 1$). In other words, the magnetic flux brought into the disk by mass
accretion from envelope is sufficient to realize $B_{z}$ shown in equation
(2.3). More precisely, the disk with size of $r\lesssim 100$ AU cannot hold
the all magnetic flux of the cloud core in it and the magnetic flux piles up
around the disk which may cause the ion-neutral drift around the disk (Li et
al., 2011; Masson et al., 2016; Tsukamoto et al., 2020; Marchand et al., 2020;
Zhao et al., 2021). Although the vertical magnetic field strength of the disk
given by equation (2.3) may seem very large, it turns out to be very small
compared to the magnetic flux possessed by the molecular cloud core.

To connect the disk size with the angular momentum accretion rate from the
envelope, we estimate the angular momentum removal rate by the magnetic
torque, which is calculated as,

$\displaystyle\dot{J}_{\rm out}$ $\displaystyle=\int^{r_{\rm
disk}}_{0}r\frac{B_{\phi,s}B_{z}}{4\pi H}~{}2H~{}2\pi rdr=\dot{M}_{\rm
disk}\sqrt{GMr_{\rm disk}}$ $\displaystyle=1.5\times
10^{40}\left(\frac{\dot{M}_{\rm disk}}{10^{-6}~{}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)$ $\displaystyle\left(\frac{M_{\rm star}}{0.3\thinspace
M_{\odot}}\right)^{\frac{1}{2}}\left(\frac{r_{\rm disk}}{100~{}{\rm
AU}}\right)^{\frac{1}{2}}{\rm g~{}cm^{2}~{}s^{-2}}.$ (31)

This simple form of $\dot{J}_{\rm out}$ is immediately derived from equations
(1) and (4), i.e., the angular momentum removal should balance with the
angular momentum supply by radial advection at each radius.

If the disk evolves by the magnetic braking and is in a steady state, the
angular momentum removal from the disk should balance with the angular
momentum supply from the envelope. Thus, we can assume $\dot{J}_{\rm
out}=\dot{J}_{\rm env}~{}(\equiv\dot{M}_{\rm env}j_{\rm env})$ where
$\dot{M}_{\rm env}$ and $j_{\rm env}$ are the mass accretion rate from the
envelope to disk and specific angular momentum of the accretion flow,
respectively. By further assuming $\dot{M}_{\rm disk}=\dot{M}_{\rm env}$ i.e.,
the mass accretion rate from the envelope to disk equals to the mass accretion
rate within the disk, we can estimate the disk size as

$\displaystyle r_{\rm disk}=160\left(\frac{j_{\rm
env}}{10^{-3}{\rm~{}km~{}pc~{}s^{-1}}}\right)^{2}\left(\frac{M_{\rm
star}}{0.3\thinspace M_{\odot}}\right)^{-1}{\rm AU},$ (32)

which is actually showing that the disk radius is equal to the centrifugal
radius of the accretion flow,

$\displaystyle r_{\rm disk}=r_{\rm cent}\equiv\frac{j_{\rm env}^{2}}{GM_{\rm
star}}.$ (33)

This result may seem obvious at first glance. Note, however, that this does
not hold for disks that evolve by internal angular momentum transport within
them (e.g., viscous accretion due to MRI or GI).

If the disk evolves by internal viscous accretion, the angular momentum is
reserved in the disk and the radius of the disk is determined by the total
angular momentum having brought into the disk. Thus, the disk radius should
depend on the entire accretion history, not on the instantaneous angular
momentum accretion from the envelope. In other words, $r_{\rm disk}$ should be
calculated from the equation

$\displaystyle\int^{r_{\rm disk}}_{0}\Sigma~{}r^{2}\Omega~{}2\pi
rdr=\int^{t}_{0}\dot{M}_{\rm env}(t)j_{\rm env}(t)dt.$ (34)

In this case, the specific angular momentum at the outer edge of the disk is
in general different from the specific angular momentum of the envelope
accretion at the envelope-to-disk boundary (see, figure 2 of Takahashi et al.,
2016, for example). This provides an observational test to determine whether
the disk is evolving by magnetic braking or by internal angular momentum
transport. We discuss this point at the end of the next section.

##  3 Disk evolution under the mass and angular momentum accretion from the
envelope

In this section, we discuss the time evolution of protoplanetary disks in
Class 0/I phase by combining the disk model constructed in the previous
section with a model of an envelope accretion.

Figure 5 shows the schematic picture of envelope-disk system considered in
this section. We assume that mass and angular momentum accrete to the outer
edge of the disk. This is based on the observations that Class 0/I YSOs
commonly have bipolar outflows in the normal direction of the disk, making the
accretion from the normal direction difficult. Furthermore, we assume that the
total mass of the disk and protostar is equal to the enclosed mass inside a
certain radius of the collapsing envelope. By relating the mass and angular
momentum accretion, and enclosed mass of envelope models to the mass accretion
rate in the disk, the angular momentum extraction rate, and the mass of the
central protostar and disk of our disk model, we can investigate the time
evolution of the disk with the mass accretion from the envelope.

Figure 5:  The schematic figure of the system we consider in this section.

###  3.1 Envelope accretion model

For the envelope accretion model, we consider the pressure free spherical
symmetric collapse of Vorobyov & Basu (2005).

The equation of motion of spherical shell is

$\displaystyle\frac{dv_{r}}{dt}=-G\frac{M(r)}{r^{2}},$ (35)

where $v_{r}$ and $r$ are the radial velocity and radius of the spherical
shell. $M(r)$ is the enclosed mass within the $r$. By solving this equation
with the assumption that $v_{r}=0$ at $t=0$, we have

$\displaystyle\sqrt{2G\frac{M(r_{0})}{r_{0}^{3}}}t$
$\displaystyle=\arccos\left(\sqrt{\frac{r}{r_{0}}}\right)+\frac{1}{2}\sin\left(2\arccos\left(\sqrt{\frac{r}{r_{0}}}\right)\right),$
(36)

where $r_{0}$ is the radius at $t=0$. By solving this equation for given $t$
and $r$, we have $r_{0}(r,t)$. Then, the radial velocity, density, and mass
accretion, rate, and specific angular momentum can be obtained as,

$\displaystyle v_{r}(r,t)$
$\displaystyle=\sqrt{2GM(r_{0})\left(\frac{1}{r}-\frac{1}{r_{0}}\right)},$
(37) $\displaystyle\rho(r,t)$
$\displaystyle=\rho_{0}(r_{0})\frac{r_{0}^{2}}{r^{2}}\frac{dr_{0}}{dr},$ (38)
$\displaystyle M(r,t)$ $\displaystyle=M(r_{0}),$ (39)
$\displaystyle\dot{M}_{\rm env}(r,t)$ $\displaystyle=4\pi
r^{2}\rho(r,t)v_{r}(r,t),$ (40) $\displaystyle j_{\rm env}(r,t)$
$\displaystyle=\epsilon_{\rm mag}f_{J}r_{0}^{2}\Omega_{0}(r_{0}).$ (41)

where $\Omega_{0}(r_{0})$ is the angular velocity at the equator of the
spherical shell. $f_{J}$ is the correction factor that arises when the angular
momentum of the rotating spherical shell is described by the angular velocity
at the equator of the spherical shell. For example $f_{J}=2/3$ if the shell is
rigidly rotating (note that the moment of inertia is $I=2/3Mr^{2}$ for the
rigidly rotating shell), and $f_{J}=0.72$ if $\Omega\propto r^{-1/2}$. For
arbitrary angular velocities of $\Omega\propto r^{\alpha}$,
$f_{J}=\sqrt{\pi}\Gamma(2+\alpha/2)/(2\Gamma(5/2+\alpha/2))$ for $\alpha>-4$.
$\epsilon_{\rm mag}$ is a factor to mimic the magnetic braking in the
envelope. We assume $\epsilon_{\rm mag}=1/3$ just for simplicity. Note that
the early studies have suggest the $\epsilon_{\rm mag}\sim 1/3$ (Basu &
Mouschovias, 1994; Tomisaka, 2002); but it may change in the late accretion
phase or magnetic field geometry (e.g., Joos et al., 2012; Marchand et al.,
2020). We believe that determining the value of $\epsilon_{\rm mag}$ precisely
is one of the major issues for future research.

Here,

$\displaystyle\frac{dr}{dr_{0}}$ $\displaystyle=\frac{r}{r_{0}}$
$\displaystyle-\sqrt{\frac{G}{2M(r_{0})r_{0}}}\left(\frac{dM(r_{0})}{dr_{0}}-\frac{3M(r_{0})}{r_{0}}\right)\sin\left(2\arccos\sqrt{\frac{r}{r_{0}}}\right).$
(42)

By specifying the initial density and rotation profile, we obtain
$M(r,t),\dot{M}_{\rm env}(r,t),j_{\rm env}(r,t)$. In this paper, we consider
the Bonnor-Ebert sphere as the initial density profile,

$\displaystyle\rho_{0}(r)=\varrho_{0}\Psi_{\rm BE}(r/a)~{}{\rm for},$ (43)

and

$\displaystyle a=c_{\rm s,iso}\left(\frac{1}{4\pi G\varrho_{0}}\right)^{1/2}.$
(44)

where $\Psi_{\rm BE}$ is non-dimensional density profile of the critical
Bonnor-Ebert sphere, $\varrho_{0}$ is the central density, and $R_{c}=1.82a$
is the characteristics radius of the Bonnor-Ebert sphere.

$\Psi_{\rm BE}(\xi)\equiv\exp(-\psi(\xi))$ where $\xi\equiv r/a$ is calculated
from the Lane-Emden equation of isothermal gas with the boundary condition of
$\psi(\xi)=d\psi(\xi)/d\xi=0$ at $\xi=0$,

$\displaystyle\frac{d}{d\xi}\left(\xi^{2}\frac{d\psi}{d\xi}\right)=\xi^{2}\exp(-\psi).$
(45)

A Bonnor-Ebert sphere is determined by specifying central density
$\varrho_{0}$, the ratio of the central density to density at $R_{c}$,
$\varrho_{0}/\rho_{0}(R_{c})$. In this study, we adopted the values of
$\varrho_{0}=2.8\times 10^{-18}~{}{\rm g~{}cm}^{-3}$,
$\varrho_{0}/\rho_{0}(R_{c})=14$. The enclosed mass within $R_{c}=5.3\times
10^{3}$ AU and $10^{4}$ AU are $0.53\thinspace M_{\odot}$ and $1\thinspace
M_{\odot}$, respectively. We assume a rigid rotation with the angular velocity
of $\Omega_{0}=10^{-13}{\rm~{}s^{-1}}$, to be consistent with the
observational results (e.g., Tobin et al., 2011).

Figure 6 shows the radial profile of density, enclosed mass, mass accretion
rate, and specific angular momentum of our envelope model. We can see that
enclosed mass, mass accretion rate, and specific angular momentum are almost
radially constant in $r\sim 100$ AU in late accretion phase (thick solid
lines). Thus, we use enclosed mass at $r=100$ AU as the total mass of central
protostar and disk

$\displaystyle M(r={\rm 100AU},t)=M_{\rm tot}$ (46)

the mass accretion rate at $r=100$ AU as the envelope-to-disk mass accretion
rate,

$\displaystyle\dot{M}_{\rm env}(r={\rm 100AU},t)=\dot{M}_{\rm env}$ (47)

and the specific angular momentum at $r=100$ AU as the specific angular
momentum of accretion flow to the disk.

$\displaystyle j(r={\rm 100AU},t)=j_{\rm env}$ (48)

Figure 7 shows the time evolution of enclosed mass, mass accretion rate, and
specific angular momentum of our envelope model at $r=100$ AU. In our collapse
model, it takes $2\times 10^{5}$ yr for the total mass of the central star and
disk $M_{\rm tot}$ to reach $\sim 1\thinspace M_{\odot}$. During this period,
the mass accretion rate reaches $\sim 10^{-5}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$ at $M_{\rm tot}=0.1\thinspace M_{\odot}$ and
decreases to $\sim 4\times 10^{-6}\thinspace M_{\odot}\thinspace{\rm yr^{-1}}$
at $M_{\rm tot}=1\thinspace M_{\odot}$. On the other hand, the specific
angular momentum of the accretion flows continues to increase monotonically
and reaches $j_{\rm env}\sim 10^{21}~{}{\rm cm}^{2}~{}{\rm s}^{-1}$ at $M_{\rm
tot}=1\thinspace M_{\odot}$.

Figure 6:  Radial profiles of density, enclosed mass, mass accretion rate, and
specific angular momentum of the collapse model. Thick solid lines show
profiles of late accretion phase and black, red, blue, and cyan lines show the
profiles at $t=4\times 10^{4},5\times 10^{4},1\times 10^{5},2\times 10^{5}$
yr, respectively. Thin solid black line in density profile shows initial
density. Dotted lines in density profile show the profiles of prestellar
collapse phase and black, red, and blue lines show the profiles at
$t=3.0\times 10^{4},3.1\times 10^{4},3.6\times 10^{4}$ yr, respectively.

Figure 7:  Time evolution of enclosed mass, and specific angular momentum of
the collapse model at $r=100$ AU. The mass accretion rate and specific angular
momentum are shown as a function of enclosed mass.

###  3.2 Connection between the envelope model and the disk model

Here, we explicitly specify the relationship between the quantities of the
envelope model and those of the disk model. To connect the accretion model
described in the previous subsection with our disk model, we assume that the
total mass of central protostar and disk is given as

$\displaystyle M_{\rm tot}=M_{\rm disk}+M_{\rm star},$ (49)

that the mass accretion rate from envelope to be disk mass accretion rate is
given as,

$\displaystyle\dot{M}_{\rm disk}=\dot{M}_{\rm env},$ (50)

and that the angular momentum accretion rate from envelope is equal to the
angular momentum removal rate,

$\displaystyle\dot{J}_{\rm out}=\dot{J}_{\rm env}=\dot{M}_{\rm env}j_{\rm
env}.$ (51)

Hereafter, we refer to this disk-envelope model as ”co-evolution disk model”.

On the other hand, we also consider a disk that evolves by internal angular
momentum transport due to gravitational instability for comparison. Hereafter,
we refer to this model as ”GI disk model”. In this GI disk model, we assume

$\displaystyle\int^{r_{\rm disk}}_{0}\Sigma_{\rm GI}~{}r^{2}\Omega~{}2\pi
rdr=\int^{t}_{0}\dot{M}_{\rm env}(t)j_{\rm env}(t)dt,$ (52)

instead of equation (51) because total angular momentum in the disk is not
removed by internal viscous accretion process.

To obtain disk mass and protostar mass, we numerical solve equations (51) (for
co-evolution disk model) or (52) (for GI disk model) and equation (49) to
obtain $M_{\rm disk}$ and $M_{\rm star}$. Although the angular velocity of the
disk can be different from the angular velocity of the Keplerian rotation when
the disk mass is not negligible compared to the central protostar mass, we
simply assume

$\displaystyle\Omega=\sqrt{\frac{GM_{\rm star}}{r^{3}}},$ (53)

even for the gravitationally unstable disk.

###  3.3 Evolution of disk radius

Figure 8 shows the disk radius as a function of the total mass, which can be
regarded as the time evolution of the disk radius. The co-evolution model
predicts the disk radius of several AU at $M_{\rm tot}\sim 10^{-1}\thinspace
M_{\odot}$ to several 100 AU in $M_{\rm tot}\sim 1\thinspace M_{\odot}$. This
relatively compact disk seems to be consistent with the recent observations of
Class 0/I YSOs (Yen et al., 2017; Maury et al., 2019).

On the other hand, radius of GI disk model is almost 10 times larger than the
co-evolution model and several $10$ AU at $M_{\rm tot}\sim 10^{-1}\thinspace
M_{\odot}$ and reaches several 1000 AU at $M_{\rm tot}\sim 1\thinspace
M_{\odot}$. This radius especially around $M_{\rm tot}\sim 1\thinspace
M_{\odot}$ is very large and seems to be inconsistent with the observations.
Note that, the radius of GI disk is the lower limit among disks which evolves
with the internal angular momentum transport, because the surface density of
GI disk is the upper limit and given total angular momentum, smaller surface
density causes larger disk radius. Thus, we conclude in general that disk
evolution models in which disks evolve with internal angular momentum
transport mechanisms predict too large disk radius and at least a certain
amount of angular momentum must be extracted from the disk during Class 0/I
phase.

Figure 8:  Time evolution of disk radius of the co-evolution model (black) and
GI disk model (red). The temperature for GI disk model is assumed to be
$T_{0}=150$ K.

###  3.4 Evolution of disk mass

Figure 9 shows the evolution of the mass of disk and protostar. The disk mass
of co-evolution model predicts the disk mass of several $10^{-3}\thinspace
M_{\odot}$ at $M_{\rm tot}\sim 10^{-1}\thinspace M_{\odot}$ and mass of $\sim
10^{-1}\thinspace M_{\odot}$ at $M_{\rm tot}\sim 1\thinspace M_{\odot}$ (which
depends on the ionization rate as $\zeta^{-1/3}$). The disk mass is much
smaller than the central protostar mass. On the other hand, the disk mass of
GI disk model is very large and it is larger than the central protostar mass,
$M_{\rm disk}/M_{\rm star}\gtrsim 1$. Although $M_{\rm disk}/M_{\rm
star}\gtrsim 1$ is actually obtained in the hydrodynamics simulations starting
from collapsing cloud core (e.g., Tsukamoto & Machida, 2011; Machida et al.,
2011), such a extremely massive disk seems to be inconsistent with the
observations.

Figure 9:  Time evolution of mass of disk (solid lines) and central protostar
(dashed lines) of the co-evolution model (black, blue) and GI disk model
(red). The temperature is assumed to be $T_{0}=150$ K. $\zeta$ for black and
blue lines is $\zeta=10^{-18}{\rm~{}s^{-1}}$, $\zeta=10^{-17}{\rm~{}s^{-1}}$,
respectively.

###  3.5 Evolution of disk magnetic-flux

Figure 10 shows the time evolution of the mass-to-flux ratio $\mu$ of co-
evolution model. With this figure, we investigate whether the magnetic flux
provided by envelope accretion is sufficient to realize the vertical magnetic
field of co-evolution disk model. If the magnetic flux of the disk is larger
than that of the initial molecular cloud core, the vertical magnetic field of
our disk cannot be realized and our model becomes inappropriate.

The figure shows that $\mu\sim 30$ at $M_{\rm tot}\sim 10^{-1}\thinspace
M_{\odot}$. On the other hand, the mass-to-flux ratio of the cloud core is
$\mu=O(1)$ (Crutcher, 2012). Thus, it indicate that sufficient magnetic flux
is supplied to the disk in the early evolutionary phase. $\mu$ monotonically
decreases mainly due to the increase of the disk radius. When $M_{\rm tot}\sim
1\thinspace M_{\odot}$ (and $r_{\rm disk}\sim 300$ AU with our envelope
model), $\mu\sim 3$ which is the same order of the mass-to-flux ratio of the
cloud core. Thus, if the disk grows to several $100$ AU, the disk can reserve
the same level of magnetic flux of the cloud core. Furthermore, if $\mu$ of
the core is $\mu\gtrsim 3$, the magnetic flux of the cloud core is not
sufficient to realize vertical magnetic field of equation (2.3). Therefore, we
must be careful about applying the disk model in this paper to model the outer
region of several hundred AU sized disks.

Figure 10:  Time evolution of mass-to-flux ratio of the disk-protostar system
normalized by $\left(\frac{M}{\Phi}\right)_{\rm
crit}=\frac{0.53}{3\pi}\left(\frac{5}{G}\right)^{1/2}$ of co-evolution model.
The temperature and $\zeta$ are assumed to be $T_{0}=150$ K and
$\zeta=10^{-18}{\rm~{}s^{-1}}$, respectively.

###  3.6 Specific angular momentum profile from disk to envelope

So far, we have compared the co-evolution disk model with the GI disk model,
and we pointed out that the co-evolution disk model predicts the radius and
mass of the disk more consistent with observations than the GI disk model.

Is there any way to more rigorously distinguish these two disk evolution
models from observations? (A more generalized question would be a way to
distinguish the disk models in which the angular momentum removal balances the
angular momentum supply and disk models which evolve by the internal angular
momentum transport mechanism.) We propose that these two (types of) models can
be distinguished by examining the radial profile of specific angular momentum
(or rotational velocity) from the envelope to the disk.

For disks that evolve by angular momentum transfer within the disk, such as
viscous accretion disks, the disk size is determined by the total angular
momentum brought in during the accretion history, since the total angular
momentum of the disk is conserved (equation (52)). In this case, the specific
angular momentum of the accreting flow at a given instant does not generally
coincide with the specific angular momentum at the outer edge of the disk.
Therefore, there must be a jump in the specific angular momentum profile at
the envelope-disk boundary.

On the other hand, in a steady accretion disk in which the angular momentum
extraction from the disk is balanced by the angular momentum supply to the
disk, the specific angular momentum of the accretion gas and that at the outer
edge of the disk do coincide, as seen in equation (33). In this case, the
angular momentum profiles are continuously connected at the envelope-disk
boundary. Therefore, we can determine whether the disk is evolving according
to viscous evolution by MRI, GI, etc. or by the magnetic braking from the
existence of jumps in the specific angular momentum profiles.

We show the radial profile of the specific angular momentum (and rotation
velocity) of co-evolution disk model (solid lines) and GI disk model (dashed
lines) in figure 11. In the co-evolution disk model, the specific angular
momentum continuously transit from $j\propto r^{0}$ to $j\propto r^{1/2}$ (or
equivalently, the rotation velocity continuously transit from
$v_{\rm\phi}\propto r^{-1}$ to $v_{\rm\phi}\propto r^{-1/2}$) without jump. On
the other hand, in the GI disk, the specific angular momentum increases about
a factor of two at the boundary between the disk and the envelope throughout
the entire accretion phase. The difference of the disk rotation profiles
between two models at the same epoch are due to the difference of the
protostar masses (we assume $\Omega=\sqrt{GM_{\rm star}/r^{3}}$). Because the
specific angular momentum profile in the disk depends weakly on radius as
$j\propto r^{1/2}$, the jump if it exists would be observable with the
observations in which the disk is spatially resolved.

Note that the magnitude of the jump of GI disk in the figure is a conservative
estimate. In the GI disk, the mass of the central star is much smaller than
the mass of the disk, but we only consider the mass of the central star to
calculate $\Omega$. In reality, at the outer edge of the GI disk, the radial
self-gravity of the disk may cause the faster angular velocity (or shallower
radial profile of $\Omega$ than Keplerian rotation profile) (see, e.g.,
Tsukamoto et al., 2015c). In this case, the magnitude of the jump is expected
to be larger. It is noteworthy that there is a factor of two increase even
with such a conservative estimate.

Figure 11:  Radial profiles of specific angular momentum and rotation velocity
of co-evolution model (solid) and GI disk model (dashed) at $t=6\times 10^{4}$
(black;co-evolution model, gray;GI model), $8\times 10^{4}$ (red;co-evolution
model, purple;GI model), $1\times 10^{5}$ yr (green;co-evolution model,
yellow-green;GI model), respectively. The temperature and $\zeta$ are assumed
to be $T_{0}=150$ K and $\zeta=10^{-18}{\rm~{}s^{-1}}$, respectively.

##  4 Summary and Discussion

###  4.1 Summary

In this paper, we analytically investigate the disk structure predicted from
”co-evolution disk model”. In our disk model, we assume that the dust grains
have grown sufficiently and ambipolar resistivity is determined by the balance
between ionization and gas-phase recombination. We also assume that the disk
evolves with magnetic braking rather than internal viscous accretion. We
utilize the equations (1) to (6) which can be regarded as the modified steady
accretion disk model, and the equation (7) which describes the sufficient dust
growth in the disk, to obtain the disk structure. We have confirmed in our
previous paper TMI23 that these equations describe the disk structures
obtained by 3D non-ideal MHD simulations with dust growth very well (within a
factor of three).

Our disk structure is shown by equation (2.3) to (2.3), and is described as a
function of central stellar mass, mass accretion rate, gas phase recombination
rate, gas temperature, and cosmic ray ionization rate, and does not include
free-parameter such as viscous parameter $\alpha$. With this disk model, we
discussed whether GI and MRI can develop or not during Class 0/I phase. Then,
by combining an analytical model of envelope accretion with the disk model, we
made predictions for the disk mass and radius during Class 0/I phase.

Our main findings are summarized as follow.

  * •

Disk radius is estimated to be several AU to $\sim 100$ AU during the
accretion phase which is about 10 times smaller than the disk evolving with
angular momentum conservation.

  * •

Disk mass is estimated to be several $10^{-3}\thinspace M_{\odot}$ to $\sim
10^{-1}\thinspace M_{\odot}$ during the accretion phase which is consistent
with the estimate from the observations.

  * •

With typical disk ionization rates ($\zeta\gtrsim 10^{-18}{\rm~{}s^{-1}}$) and
moderate mass accretion rate ($\dot{M}_{\rm disk}>10^{-8}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$), magneto-rotational instability is
prohibited in the disk.

  * •

Plasma $\beta$ at midplane depends on the mass accretion rate (i.e., disk
evolutionary stage) as $\beta\propto\dot{M}^{-2/3}_{\rm disk}$, and $\beta\sim
30$ at $\dot{M}_{\rm disk}\sim 10^{-6}\thinspace M_{\odot}\thinspace{\rm
yr^{-1}}$ and $\beta\sim 700$ at $\dot{M}_{\rm disk}\sim 10^{-8}\thinspace
M_{\odot}\thinspace{\rm yr^{-1}}$ at $r=10$ AU.

  * •

Specific angular momentum (or rotation velocity) profile at the envelope to
disk boundary is continuous if the disk angular momentum evolves only with
magnetic braking.

###  4.2 Discussion

####  4.2.1 The relation of the disk size and the specific angular momentum of
accretion flow

An important prediction obtained from our analysis of the disk model is that
the disk radius coincides with the centrifugal radius of the accreting flow
when the disk evolves with magnetic braking rather than viscous accretion.
This follows just from the assumption of steadiness of the disk (i.e.,
constant $\dot{M}_{\rm disk}$) and the assumption that angular momentum
extraction from the disk is balanced with angular momentum supply to the disk,
and does not depend on the detail of the assumptions on $\eta_{A}$.

This prediction can be verified by detailed observations of the specific
angular momentum (or rotation velocity) profile around the disk-envelope
boundary. If the specific angular momentum is continuous at the disk-envelope
boundary, it suggests that disk angular momentum evolution is determined by
the angular momentum extraction due to the magnetic braking. On the other
hand, if the specific angular momentum is discontinuous (or has the jump) at
the disk-envelope boundary, it suggests that disk angular momentum evolution
is determined by the internal (viscous) angular momentum transport (see figure
11).

There have been several observations that investigate the radial profile of
specific angular momentum or rotation velocity of Class 0/I YSOs. Recent ALMA
velocity profile observations seem to show that rotational velocity profile
from the envelope to disk of the Class 0/I phase seems to be continuous and no
jump at the disk-envelope boundary (Ohashi et al., 2014; Aso et al., 2015; Yen
et al., 2017). If future high-resolution observations confirm the continuous
connection of the specific angular momentum of the envelope and the outer edge
of the disk, this indicates that the evolution of the angular momentum and
radius of the disk is determined by the extraction of angular momentum from
the disk (i.e., magnetic braking) rather than by internal viscous evolution of
the disk.

####  4.2.2 Caveat of co-evolution disk model and Future Prospects

In the present study, ambipolar diffusion is assumed to be the most efficient
non-ideal MHD effect. When dust growth proceeds, it is shown that Ohmic
dissipation and Hall effects are weaker than ambipolar diffusion for cases
where the magnetic field is relatively strong (Tsukamoto & Okuzumi, 2022).
However, when the magnetic field becomes weak, these effects may possibly play
the role. In future studies, we will examine the disk structure where Ohm
dissipation or the Hall effect play the dominant role. Note, however that,
when Ohm dissipation is dominant, there may not be a unique steady structure
because $\eta_{O}$ is independent of the magnetic field.

In this study, we assumed that the dust grans were well grown over the entire
disk. However, at the edge of the disk, it is expected that there is a supply
of ${\rm\mu}{\rm m}$ sized dust grains from the envelope. In that case, more
efficient magnetic diffusion would be expected at the edge of the disk, and
our approximation for $\eta_{A}$ may not be justified. It is important to
study the physical state of the outer edge of the disk in detail in the future
studies which consider the dust growth and envelope-to-disk accretion (e.g.,
Tsukamoto et al., 2017b). Note, however, that the above discussion of the
specific angular momentum jump at the outer edge of the disk is not affected
by the dust size as long as it evolves with magnetic braking, since it does
not depend on the details of $\eta_{A}$.

We also note that the detailed modeling of dust coagulation and fragmentation
and its impact on $\eta_{A}$ is important because the dust fragmentation can
cause a variety of dust size distributions (Birnstiel et al., 2011). As shown
in Lebreuilly et al. (2023), dust fragmentation can changes $\eta_{A}$ profile
especially in the high density region. Once a model of $\eta_{A}$ that takes
into account the effects of dust coagulation/fragmentation is obtained, it
would be possible to predict the disk structure according to the dust growth
model and resulting $\eta_{A}$ using the equations in this paper.

Although there are issues to be solved in the future as described above, we
believe that our co-evolution disk model provides a new perspective on the
structure and evolution of protoplanetary disks.

## Acknowledgments

I thank Dr. Shinsuke Takasao, Dr. Satoshi Okuzumi, and Mr. Ryoya Yamamoto for
the fruitful discussions. I also thank Prof. Shu-ichiro Inutsuka and Prof.
Masahiro Machida for their continuous encouragement and the comments. We also
thank the referee, Dr. Pierre Marchand for his helpful and insightful
comments. This work is supported by JST FOREST.

## References

  * Aso et al. (2015) Aso, Y., Ohashi, N., Saigo, K., et al. 2015, ApJ, 812, 27
  * Bai (2017) Bai, X.-N. 2017, ApJ, 845, 75
  * Bai & Stone (2011) Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144
  * Basu & Mouschovias (1994) Basu, S., & Mouschovias, T. C. 1994, ApJ, 432, 720
  * Birnstiel et al. (2011) Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11
  * Carrasco-González et al. (2019) Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71
  * Crutcher (2012) Crutcher, R. M. 2012, ARA&A, 50, 29
  * Dapp & Basu (2010) Dapp, W. B., & Basu, S. 2010, A&A, 521, L56
  * Dapp et al. (2012) Dapp, W. B., Basu, S., & Kunz, M. W. 2012, A&A, 541, A35
  * Fleming & Stone (2003) Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  * Gressel et al. (2015) Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84
  * Guilet & Ogilvie (2012) Guilet, J., & Ogilvie, G. I. 2012, MNRAS, 424, 2097
  * Guilet & Ogilvie (2014) —. 2014, MNRAS, 441, 852
  * Harsono et al. (2018) Harsono, D., Bjerkeli, P., van der Wiel, M. H. D., et al. 2018, Nature Astronomy, 2, 646
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  * Joos et al. (2012) Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128
  * Kobayashi et al. (2023) Kobayashi, Y., Takaishi, D., & Tsukamoto, Y. 2023, MNRAS, arXiv:2303.04424
  * Kuffmeier et al. (2020) Kuffmeier, M., Zhao, B., & Caselli, P. 2020, A&A, 639, A86
  * Lebreuilly et al. (2023) Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326
  * Li et al. (2011) Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180
  * Machida (2014) Machida, M. N. 2014, ApJ, 796, L17
  * Machida & Basu (2019) Machida, M. N., & Basu, S. 2019, ApJ, 876, 149
  * Machida et al. (2011) Machida, M. N., Inutsuka, S., & Matsumoto, T. 2011, PASJ, 63, 555
  * Marchand et al. (2021) Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2021, A&A, 649, A50
  * Marchand et al. (2023) Marchand, P., Lebreuilly, U., Mac Low, M. M., & Guillet, V. 2023, A&A, 670, A61
  * Marchand et al. (2020) Marchand, P., Tomida, K., Tanaka, K. E. I., Commerçon, B., & Chabrier, G. 2020, ApJ, 900, 180
  * Masson et al. (2016) Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32
  * Maury et al. (2019) Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76
  * McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36
  * Mouschovias & Spitzer (1976) Mouschovias, T. C., & Spitzer, Jr., L. 1976, ApJ, 210, 326
  * Ohashi et al. (2014) Ohashi, N., Saigo, K., Aso, Y., et al. 2014, ApJ, 796, 131
  * Pinto & Galli (2008) Pinto, C., & Galli, D. 2008, A&A, 484, 17
  * Shu (1983) Shu, F. H. 1983, ApJ, 273, 202
  * Smith & Adams (1984) Smith, D., & Adams, N. G. 1984, ApJ, 284, L13
  * Suzuki & Inutsuka (2014) Suzuki, T. K., & Inutsuka, S.-i. 2014, ApJ, 784, 121
  * Takahashi et al. (2016) Takahashi, S. Z., Tomida, K., Machida, M. N., & Inutsuka, S.-i. 2016, MNRAS, 463, 1390
  * Takasao et al. (2018) Takasao, S., Tomida, K., Iwasaki, K., & Suzuki, T. K. 2018, ApJ, 857, 4
  * Tobin et al. (2011) Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2011, ApJ, 740, 45
  * Tobin et al. (2020) Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130
  * Tomida et al. (2015) Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117
  * Tomisaka (2002) Tomisaka, K. 2002, ApJ, 575, 306
  * Tsukamoto et al. (2015a) Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015a, ApJ, 810, L26
  * Tsukamoto et al. (2015b) —. 2015b, MNRAS, 452, 278
  * Tsukamoto & Machida (2011) Tsukamoto, Y., & Machida, M. N. 2011, MNRAS, 416, 591
  * Tsukamoto et al. (2021) Tsukamoto, Y., Machida, M. N., & Inutsuka, S.-i. 2021, ApJ, 920, L35
  * Tsukamoto et al. (2023a) —. 2023a, PASJ, 75, 835
  * Tsukamoto et al. (2020) Tsukamoto, Y., Machida, M. N., Susa, H., Nomura, H., & Inutsuka, S. 2020, ApJ, 896, 158
  * Tsukamoto & Okuzumi (2022) Tsukamoto, Y., & Okuzumi, S. 2022, ApJ, 934, 88
  * Tsukamoto et al. (2017a) Tsukamoto, Y., Okuzumi, S., Iwasaki, K., Machida, M. N., & Inutsuka, S.-i. 2017a, PASJ, 69, 95
  * Tsukamoto et al. (2017b) Tsukamoto, Y., Okuzumi, S., & Kataoka, A. 2017b, ApJ, 838, 151
  * Tsukamoto et al. (2015c) Tsukamoto, Y., Takahashi, S. Z., Machida, M. N., & Inutsuka, S. 2015c, MNRAS, 446, 1175
  * Tsukamoto et al. (2023b) Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023b, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 317
  * Tychoniec et al. (2020) Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19
  * Vorobyov & Basu (2005) Vorobyov, E. I., & Basu, S. 2005, MNRAS, 360, 675
  * Wurster & Bate (2019) Wurster, J., & Bate, M. R. 2019, MNRAS, 486, 2587
  * Wurster et al. (2018) Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, arXiv:1802.04872
  * Wurster et al. (2016) Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037
  * Xu & Kunz (2021) Xu, W., & Kunz, M. W. 2021, MNRAS, 508, 2142
  * Yen et al. (2017) Yen, H.-W., Koch, P. M., Takakuwa, S., et al. 2017, ApJ, 834, 178
  * Zhao et al. (2018) Zhao, B., Caselli, P., Li, Z.-Y., & Krasnopolsky, R. 2018, MNRAS, 473, 4868
  * Zhao et al. (2020) Zhao, B., Caselli, P., Li, Z.-Y., et al. 2020, MNRAS, 492, 3375
  * Zhao et al. (2021) —. 2021, MNRAS, 505, 5142
  * Zhao et al. (2016) —. 2016, MNRAS, 460, 2050
  * Zhu & Stone (2018) Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34