# Co-evolution of dust grains and protoplanetary disks

Yusuke Tsukamoto1, Masahiro N. Machida2, and Shu-ichiro Inutsuka3
1Graduate Schools of Science and Engineering, Kagoshima University, Kagoshima,
Japan
2Department of Earth and Planetary Sciences, Kyushu University, Fukuoka, Japan
3Department of Physics, Nagoya University, Aichi, Japan

###### Abstract

We propose a new evolutionary process of protoplanetary disks ”co-evolution of
dust grains and protoplanetary disks”, revealed by dust-gas two-fluid non-
ideal magnetohydrodynamics simulations considering the growth of dust and
associated changes in magnetic resistivity. We found that the dust growth
significantly affects disk evolution by changing the coupling between the gas
and magnetic field. Moreover, once the dust grains sufficiently grow, the
physical quantities (e.g., density and magnetic field) of the disk are well
described by nontrivial power laws, regardless of the details of the dust
model. In this disk structure, the radial profile of density is steeper and
the disk mass is smaller than those of the model ignoring dust growth and they
are more consistent with the disk observations. We analytically derive these
power laws from the basic equations of non-ideal magnetohydrodynamics. The
analytical power laws are determined only by observable physical quantities,
e.g., central stellar mass and mass accretion rate, and do not include
difficult-to-determine parameters e.g., viscous parameter $\alpha$. Therefore,
they are applicable to various stages of disk evolution. We believe that the
disk structure provides a new basis for future studies on star and planet
formation.

###### keywords:

star formation – circum-stellar disk – methods: magnetohydrodynamics –
smoothed particle hydrodynamics – protoplanetary disk

##  1 Introduction

Protoplanetary disks are the birthplaces of planets. Planets come into being
as a result of the growth of dust grains in the disk. During their growth,
like the baby kicks, dust grains can influence the evolution of the
protoplanetary disks.

In protoplanetary disks, dust grains are not only the building blocks of the
planets, but also play a key role in determining the ionization degree of the
disk gas by adsorbing charged particles (ions and electrons) in the gas phase.
Since the ionization degree determines the magnetic resistivity, i.e., the
degree of coupling between the gas and magnetic field, the microscopic nature
(${\rm\mu}{\rm m}$ to ${\rm cm}$ scale) of the dust grains is expected to
influence the macroscopic disk evolution ($100$ AU, i.e., $10^{15}{\rm cm}$
scale) via magnetic resistivity (Zhao et al., 2016; Marchand et al., 2020;
Guillet et al., 2020; Tsukamoto & Okuzumi, 2022).

Previous studies show that non-ideal magnetohydrodynamics (MHD) effects
arising from finite resistivity, specifically ambipolar diffusion,
dramatically weaken the coupling between the magnetic field and gas, therby
enabling the formation of a disk (Tomida et al., 2015; Tsukamoto et al., 2015;
Wurster et al., 2016; Masson et al., 2016). Moreover ambipolar diffusion
determines the magnetic flux evolution in protostars (Li, 1998; Tsukamoto et
al., 2020).

Previous studies on the formation and evolution of protoplanetary disks
assumed that dust grains possess the properties (such as the size
distribution) of the interstellar medium (ISM) dust grains. However, in the
disk, the dust growth timescale is $\sim 10^{4}$ years, which is much shorter
than the lifetime of the disks ($\sim 10^{6}$ years). Thus, it is
unsatisfactory to study disk evolution with resistivity assuming ISM dust
(Tsukamoto et al., 2022).

How would the dust growth affect the ionization degree? As dust grains merge
and grow, their total surface area decreases. Therefore, the adsorption of
charged particles by the dust grains becomes ineffective, and the gas-phase
ionization degree is expected to increase and magnetic resistivity to decrease
accordingly. Recent studies on dust growth and associated changes in magnetic
resistivity have shown a decrease in magnetic resistivity (Zhao et al., 2016;
Marchand et al., 2020; Guillet et al., 2020; Tsukamoto & Okuzumi, 2022;
Kawasaki et al., 2022), and some studies have also shown changes in gas
dynamics as a result (Lebreuilly et al., 2023; Marchand et al., 2023).

However, the effect of the dust growth on the evolution of the protoplanetary
disk is still unclear, because the calculations in the aforementioned studies
were performed assuming spherical symmetry (Lebreuilly et al., 2023) or 3D
simulation until the prestellar or first core formation stage (in which the
gas is supported by the pressure gradient force) (Marchand et al., 2023).

In this study, we report simulation results of the formation and evolution of
protoplanetary disks of $\sim 10^{4}$ years after the formation of protostars
considering dust growth inside the disk, the associated change of magnetic
resistivity, and its feedback on the disk dynamics. Moreover we present an
analytical argument that explains the resulting disk structures. Based on
these results, we propose a new evolutionary process for protostars: ”co-
evolution of dust grains and protoplanetary disks”.

##  2 Methods and initial condition

###  2.1 Numerical methods

####  2.1.1 Two-fluid magneto-hydrodynamics simulations

We solve two-fluid magnetohydrodynamics equations for the dust-gas mixture.
The governing equations are given as

$\displaystyle\frac{D\rho}{Dt}$ $\displaystyle=$
$\displaystyle-\rho\nabla\cdot\mathbf{v},$ (1)
$\displaystyle\frac{D\epsilon}{Dt}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho}\cdot\\{\epsilon(1-\epsilon)\rho\Delta\mathbf{v}\\},$
(2) $\displaystyle\frac{D\mathbf{v}}{Dt}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho}\\{\nabla
P-\frac{\mathbf{J}\times\mathbf{B}}{c}\\},$ (3)
$\displaystyle\frac{D\Delta\mathbf{v}}{Dt}$ $\displaystyle=$
$\displaystyle-\frac{\Delta\mathbf{v}}{t_{\rm
stop}}+\frac{1}{\rho_{g}}[-\nabla P+\frac{\mathbf{J}\times\mathbf{B}}{c}],$
(4) $\displaystyle\frac{D\mathbf{B}}{Dt}$ $\displaystyle=$
$\displaystyle-\mathbf{B}(\nabla\cdot\mathbf{v})+(\mathbf{B}\cdot\nabla)\mathbf{v}$
(5) $\displaystyle+$ $\displaystyle
c\nabla\times\\{\eta_{O}\mathbf{J}+\eta_{A}(\mathbf{J}\times{\hat{\mathbf{B}}})\times{\hat{\mathbf{B}}}\\},$

where $\rho_{[g,d]}$ denotes the mass densities and subscripts $[g,d]$ denote
gas and dust components, respectively. $\rho=\rho_{g}+\rho_{d}$ denotes the
total density. $\epsilon=\rho_{d}/\rho$ denotes the dust-to-total-mass ratio,
$\mathbf{v}=(\rho_{g}\mathbf{v}_{g}+\rho_{d}\mathbf{v}_{d})/(\rho_{g}+\rho_{d})$
denotes the barycentric velocity of dust gas mixture where
$\mathbf{v}_{[g,d]}$ denotes the gas and dust velocity,
$\Delta\mathbf{v}=(\mathbf{v}_{d}-\mathbf{v}_{g})$ denotes the the velocity
difference between gas and dust, $P$ denotes the gas pressure. $\mathbf{J}$
denotes the electric current. $c$ denotes the speed of light. $\eta_{O}$ and
$\eta_{A}$ denote the Ohmic and ambipolar resistivities, respectively. The
details of the approximations adopted in the governing equations and numerical
method are described in Tsukamoto et al. (2021a). Our numerical simulations
consider the Ohmic and ambipolar diffusions, but ignore the Hall effect.

####  2.1.2 Dust growth

We consider dust growth with single-size approximation (Sato et al., 2016;
Okuzumi et al., 2016; Tsukamoto et al., 2017). The governing equation of dust
growth is

$\displaystyle\frac{D_{d}a_{d}}{Dt}=A_{\rm gain/loss}\frac{a_{d}}{3t_{\rm
growth}},$ (6)

where $a_{d}$ denotes the representative dust size, $t_{\rm growth}=1/(\pi
a_{d}^{2}n_{d}\Delta v_{\rm dust})$, $n_{d}$ denotes the dust number density,
$\Delta v_{\rm dust}$ denotes the collision velocity between the dust grains,
$D_{d}/Dt=\partial/\partial t+\mathbf{v}_{d}\cdot\nabla$, and

$\displaystyle A_{\rm gain/loss}={\rm min}(1,-\frac{\ln(\Delta v_{\rm
dust}/\Delta v_{\rm frag})}{\ln 5}),$ (7)

which models the collisional mass gain and loss (Okuzumi et al., 2016). We
assume $v_{\rm frag}=30~{}{\rm m}~{}{\rm s}^{-1}$. For the dust relative
velocity $\Delta v_{\rm dust}$, we consider the sub-grid scale turbulence and
Brownian motion.

For the turbulent-induced dust relative velocity $\Delta v_{\rm turb}$, we
adopt the prescription presented by Ormel & Cuzzi (2007),

$\displaystyle\Delta v_{\rm turb}$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{\delta v_{\rm Kol}}{t_{\rm Kol}}(t_{\rm
stop,1}-t_{\rm stop,2})&(t_{\rm stop,1}<t_{\rm Kol})\\\ 1.5\delta
v_{L}\sqrt{\frac{t_{\rm stop,1}}{t_{L}}}&(t_{\rm Kol}<t_{\rm stop,1}<t_{L})\\\
\delta v_{L}\sqrt{\frac{1}{1+t_{\rm stop,1}/t_{L}}+\frac{1}{1+t_{\rm
stop,2}/t_{L}}}&(t_{L}<t_{\rm stop,1})\end{cases}$

where $\delta v_{\rm Kol}=Re_{L}^{-1/4}$ and $t_{\rm Kol}=Re_{L}^{-1/2}t_{L}$
denote the eddy velocity and eddy turn-over timescale at dissipation scale and
$Re_{L}=Lv_{L}/\nu$ denotes the Reynolds number. We set the stopping time
$t_{\rm stop,1}=t_{\rm stop}(a_{d})$ and $t_{\rm stop,2}=1/2~{}t_{\rm
stop}(a_{d})$ referring to Sato et al. (2016), where $t_{\rm stop}(a_{d})$ is
calculated as in our previous study (Tsukamoto et al., 2021a).

We assume sub-grid turbulence of the “$\alpha$ turbulence model ”(Tsukamoto et
al., 2021b), in which we assume

$\displaystyle\delta v_{L}$ $\displaystyle=$ $\displaystyle\sqrt{\alpha_{\rm
turb}}c_{s},$ (9) $\displaystyle t_{L}$ $\displaystyle=$
$\displaystyle\frac{c_{s}}{a_{g}},$ (10) $\displaystyle L$ $\displaystyle=$
$\displaystyle\delta v_{L}t_{L}.$ (11)

$\alpha_{\rm turb}=2\times 10^{-3}$ denotes the dimensionless parameter that
determines the strength of the sub-grid turbulence and $a_{g}$ denotes the
gravitational acceleration. (see (Tsukamoto et al., 2021b) for the underlying
physical assumptions for this turbulence model).

####  2.1.3 Resistivity calculations

For the resistivity model, we adopt the analytical resistivity formula
described in Tsukamoto & Okuzumi (2022) in which we analytically solve the
equations for chemical equilibrium in the gas phase and detailed balance
equations for dust charging. The dust size distribution considered in
resistivity calculations is set to be

$\displaystyle\frac{{\rm d}n_{\rm d}}{{\rm d}a}=A~{}a^{-q}(a_{\rm
min}<a<a_{\rm max}),$ (12)

where $A$ denotes a constant for normalization. Here we assume that the
maximum dust size is $a_{\rm max}=a_{\rm d}$ of equation (6) (which is valid
when $q<4$). The minimum dust size $a_{\rm min}$ and power exponet $q$ are the
parameters of this study.

####  2.1.4 Sink particle

A sink particle is dynamically introduced when the density exceeds $\rho_{\rm
sink}=10^{-12}~{}{\rm g~{}cm}^{-3}$. The sink particle absorbs SPH particles
with $\rho>\rho_{\rm sink}$ within $r_{\rm sink}<1$ AU.

###  2.2 Initial conditions

We adopt the density-enhanced Bonnor-Ebert sphere, which is surrounded by a
medium with a steep density profile used in Tsukamoto et al. (2021b) as a
initial condition.

The radius of the core is $R_{c}=4.8\times 10^{3}$ AU and the enclosed mass
within $R_{c}$ is $M_{c}=1\thinspace M_{\odot}$. We adopt an angular velocity
profile of $\Omega(d)=\Omega_{0}/[\exp[10(d/(1.5R_{c}))-1]+1]$ with
$d=\sqrt{x^{2}+y^{2}}$ and $\Omega_{0}=2.3\times 10^{-13}~{}{\rm s^{-1}}$. We
assume a constant magnetic field $(B_{x},B_{y},B_{z})=(0,0,83~{}\mu G)$.

The parameter $\alpha_{\rm therm}$ ($\equiv E_{\rm therm}/E_{\rm grav}$) is
$0.4$, where $E_{\rm therm}$ and $E_{\rm grav}$ denote the thermal and
gravitational energies of the central core (without surrounding medium),
respectively. The parameter $\beta_{\rm rot}$ ($\equiv E_{\rm rot}/E_{\rm
grav}$) within the core is $0.03$, where $E_{\rm rot}$ denotes the rotational
energy of the core. The mass-to-flux ratio of the core normalized by the
critical value is $\mu/\mu_{\rm crit}=3$.

We adopt a dust density profile of
$\rho_{d}(r)=f_{dg}\rho_{g}(r)/[\exp[10(r/(1.5R_{c}))-1]+1]$, where
$f_{dg}=10^{-2}$ denotes the dust-to-gas mass ratio. The dust density profile
has the same shape with the gas density profile in $r\lesssim 1.5R_{c}$ but is
truncated at $r\geq 1.5R_{c}$. The initial (maximum) dust size is assumed to
be $a_{d}=0.1{\rm\mu}{\rm m}$.

We resolve 1 $\thinspace M_{\odot}$ with $3\times 10^{6}$ SPH particles. The
model names and parameters of the models are listed in Table 1.

Table 1:  Model name, minimum dust size $a_{\rm min}$, power exponent for dust size distribution $q$, and cosmic ray ionization rate $\zeta_{\rm CR}$. ”Y” means that the dust growth is considered in the simulation and ”N” means that the dust growth is not considered.  Model name | $a_{\rm min}[{\rm nm}]$ | $q$ | $\zeta_{\rm CR}[{\rm s^{-1}}]$ | Dust growth
---|---|---|---|---
ModelA100Q25 | 100 | 2.5 | $10^{-17}$ | Y
ModelA100Q35 | 100 | 3.5 | $10^{-17}$ | Y
ModelA5Q25 | 5 | 2.5 | $10^{-17}$ | Y
ModelA5Q35 | 5 | 3.5 | $10^{-17}$ | Y
ModelZeta18 | 100 | 2.5 | $10^{-18}$ | Y
ModelA100Fixed | 100 | - | $10^{-17}$ | N

##  3 Results

###  3.1 Co-evolution of dust grains and protoplanetary disks

Figure 1 shows the time evolution of (a) density, (b) magnetic resistivity,
(c) dust size, and (d) plasma $\beta$ in our fiducial model, ModelA100Q25.
Panel 1-a shows the formation of a disk with a size of $\sim 20$ AU at
$t_{*}=1.1\times 10^{3}$ yr (where $t_{*}$ denotes the time after protostar
formation). At this stage, the total magnetic resistivity
($\eta_{O}+\eta_{A}$) is $\sim 10^{19}~{}{\rm cm}^{2}~{}{\rm s}^{-1}$ and
relatively large (panel 1-b). This large value is attributed to the fact that
the maximum dust size remains $a_{d}<1{\rm\mu}{\rm m}$ within the disk (panel
1-c) and dust adsorption of charged particles is effective. Owing to the
efficient magnetic diffusion, the plasma $\beta$ inside the disk is
significantly high (panel 1-d).

As the time progresses, a remarkable change in the magnetic resistivity
occurs. Panel 2-b shows a significant decrease in the magnetic resistivity
within the disk. This is caused by the dust growth within the disk ($a_{d}$
reaches up to $10{\rm\mu}{\rm m}$) and reduction in the dust adsorption
efficiency. However, at this epoch, the decrease in the magnetic diffusion
does not lead to a decrease in the disk size; instead, the disk continues to
expand. The plasma $\beta$ also remains high.

The decrease in the magnetic resistivity leads to a better coupling between
the magnetic field and gas in the disk. This coupling promotes gas accretion,
which in turn transports the magnetic flux to the center and reduces the
plasma $\beta$ in the disk (panel 3-d). The strong magnetic field amplifies
the magnetic resistivity ($\eta_{A}$) in the central region (panel 3-b). Once
the dust grows sufficiently, $\eta_{A}$ is proportional to the square of the
magnetic field strength even inside the disk. Even at this point, the gas
density map (panel 3-a) indicates the presence of a relatively large disk of
$\sim 50$ AU.

The strong magnetic field in the disk causes efficient magnetic braking and
efficient mass accretion over the entire disk, leading to a decrease in the
disk size from 3-a to 4-a. However, that the density structure of the central
region is very similar in panels 3-a and 4-a. On the other hand, between
panels 2-a and 4-a, the disk size is similar, but the density structure of the
inner region is different. This indicates that the inner disk structure
transits from 2-a to 3-a as the dust grains grow.

In this way, the growth of the dust grains causes a decrease in the magnetic
resistivity and changes the magnetic activity of the disk, and ultimately
determines the evolution of the disk. Conversely, changes in the disk
structures affect the dust growth in the disk (figure 7). Based on these
results, we propose a new evolutionary process of the protoplanetary disk:
“co-evolution of dust grains and protoplanetary disks”.

Figure 1:  Time evolution of the (a) gas density, (b) magnetic resistivity
($\eta_{O}+\eta_{A}$), (c) dust size ($a_{d}$), and (d) plasma $\beta$ on the
$x$-$y$ plane with $100$ AU box. The time after protostar formation are shown
in the upper left. The red and orange arrows show the velocity field of the
dust and gas, respectively, on the $x$-$y$ plane. The red and orange lines
show their respective streamlines. The black lines are the contour of the
quantity of each panel. The contour levels are (a)
$\rho_{g}=10^{-15},10^{-14.5},\cdots,10^{-10}~{}{\rm g~{}cm}^{-3}$, (b)
$\eta_{O}+\eta_{A}=10^{16},10^{16.5},\cdots,10^{21}~{}{\rm cm}^{2}~{}{\rm
s}^{-1}$, (c) $a_{d}=10^{-5},10^{-4.5},\cdots,10^{0}{\rm cm}$, and (d)
$\beta=10^{-2},10^{-1.5},\cdots,10^{3}$.

Figure 2:  Time evolution of the azimuthally averaged radial profile of
(panels a) gas density, (b) vertical magnetic field, (c) radial velocity, and
(d) ambipolar resistivity ($\eta_{A}$). The epoch of each line is the same as
in figure 1. The black dashed lines indicate the analytical solutions
(equations (A) to (A)) in which following values $\rho_{c}=4\times
10^{4}~{}{\rm g~{}cm}^{-3}$, $C\gamma=3\times 10^{-3}{\rm
cm^{\frac{3}{2}}g^{-\frac{1}{2}}s^{-1}}$, $c_{s,ref}=190~{}{\rm m}~{}{\rm
s}^{-1}$, $\dot{M}=2\times 10^{-5}\thinspace M_{\odot}~{}{\rm yr}^{-1}$, and
$M=0.3\thinspace M_{\odot}$ are assumed. The hatched areas are regions within
a factor of three from the solutions.

###  3.2 Radial disk structure and comparison with analytical solutions

Figure 2 shows the azimuthally averaged radial profiles at the midplane of
ModelA100Q25. The red, orange, magenta, and green lines show profiles at the
epochs of panels 1-4 in figure 1, respectively.

In the early evolutionary phase when the dust size is sufficiently small
($t_{*}\lesssim 5\times 10^{3}$ yr; red and orange lines), the gas density has
a relatively shallow profile with power of $D_{\rho_{g}}\sim-1$ (red and
orange lines; where $D_{f}$ denotes the power exponent of a quantity $f$ as
$f(r)\propto r^{D_{f}}$ ). As the dust grains grow, the gas density profile
becomes steeper and appears to converge to a power law with
$D_{\rho_{g}}\sim-2$ (at $t_{*}=9\times 10^{3}$ yr; green line).

In supplementary materials, we analytically derive the steady state solutions
by assuming that the adsorption of charged particles on the dust grains is
negligible, and the ionization degree is determined by cosmic ray ionization
and gas phase recombination. The derived solution of the density profile
predicts $D_{\rho_{g}}=-\frac{9}{5}$ (equation (A)). The black dotted line and
hatched areas indicate the analytical solutions and a region within a factor
of four of the solution, respectively. Here in we have chosen the following
values for the analytical solution: $\rho_{c}=4\times 10^{4}~{}{\rm
g~{}cm}^{-3}$, $C\gamma=3\times 10^{-3}{\rm
cm^{\frac{3}{2}}g^{-\frac{1}{2}}s^{-1}}$, $c_{s,ref}=190~{}{\rm m}~{}{\rm
s}^{-1}$, to be consistent with the simulation setup and $\dot{M}=2\times
10^{-5}\thinspace M_{\odot}~{}{\rm yr}^{-1}$ and $M=0.3\thinspace M_{\odot}$
to be approximately consistent with the values at $t_{*}=9\times 10^{3}$ yr.
Our analytical solution well agrees with the simulation result not only in
terms of the power, but also in terms of the exact value.

Panel (b) shows that in the early evolutionary phase (red and orange lines),
the magnetic field in the disk is almost constant ($D_{B_{z}}\sim 0$), which
is in agreement with previous studies (Masson et al., 2016; Tsukamoto et al.,
2015) in which the dust growth is neglected. In contrast, as the dust grains
grow, the vertical magnetic field profile becomes steeper and converges to a
power law with $D_{B_{z}}\sim-1$. Our analytical solution suggests the power
exponent of $D_{B_{z}}=-\frac{9}{10}$ (equation (A)) and agree with the
simulation results at $t_{*}=9.0\times 10^{3}$ yr (green line). The black
dotted line indicates our analytical solution with the same parameters used in
the density profile and the hatched area are a region within a factor of four
from the solution. The black dotted line confirms that our analytical solution
simultaneously reproduces the density and magnetic field with a single set of
parameters and quantitative agrees with the simulation results.

Panel (c) shows that, when the dust size is small, the absolute value of
radial velocity is $|v_{r}|\lesssim 10~{}{\rm m}~{}{\rm s}^{-1}$ (red line).
As the dust grows, it increases in the order of $100~{}{\rm m}~{}{\rm
s}^{-1}$. The proposed analytical solution suggests a power of
$D_{v_{r}}=\frac{1}{10}$ and value of $\sim 100~{}{\rm m}~{}{\rm s}^{-1}$.
Although $v_{r}$ possesses relatively strong time fluctuation, the simulation
results at $t_{*}=9.0\times 10^{3}$ yr (green line) agrees with the analytical
solutions.

Panel (d) shows that the $\eta_{A}$ profile in the simulation is almost
radially constant once the dust grows (magenta and green lines). On the other
hand, our analytical model predicts slightly positive power law with
$D_{\eta_{A}}=\frac{9}{10}$. The value estimated from the analytical solution
is also slightly larger than the simulation results (upto about factor of
four). This may be due to the assumptions made in the analytical solution such
as $\eta_{A}$ being exactly given by $\eta_{A}=B^{2}/(4\pi C\gamma\rho^{3/2})$
(which is realized when the dust adsorption becomes negligible and ionization
degree is determined by the balance between (cosmic-ray) ionization and gas-
phase recombination (Shu, 1983) or the vertically isothermal disk. They
introduce slight differences between the simulation and analytical solution.
Note that the more the model evolves, the more consistent it becomes with the
analytical solution; see green line in panel (d) of figure 4.

###  3.3 Diversity and universality of disk evolution

As seen in the previous section, the disk evolution is significantly affected
by the dust growth in it. Furthermore, once the dust grains sufficiently
grows, the disk structure of our fiducial model is well described by the power
laws analytically derived in Appendix A. In this section, we examine the
impact of the minimum dust size and dust power exponent on the evolution of
the disk. Moreover we examine the effect of the cosmic-ray ionization rate.
Through these considerations, we discuss the diversity and universality of
protoplanetary disk evolution.

Before discussing the simulation results, we summarize how the minimum dust
size $a_{\rm min}$ and the power exponent $q$ affect the magnetic resistivity
based on our previous studies (Tsukamoto & Okuzumi, 2022). The adsorption of
dust grains, which plays the primary role in determining the resistivity
depends on the dust total cross-section, which depends on $q$ and $a_{\rm
min}$. In the case where $q$ is less than 3, the maximum dust size determines
the total surface area. Conversely, in a case where $q$ is larger than 3, both
the maximum and minimum dust sizes affect the total surface area.
Consequently, when $q$ is large (or the size distribution is steep), the
influence of dust growth tends to be weaker.

Another important factor that influences resistivity is conductivity generated
by the dust itself. If a significant amount of dust with a size of $\lesssim
10{\rm nm}$ is present, the dust itself contributes to conductivity. In such
case, resistivity at a high density is smaller than that when the minimum size
is, for instance, $100~{}{\rm nm}$. This effect is pronounced in a case when
the dust has not grown and the size distribution is steep (i.e., $q$ is
large).

####  3.3.1 Comparison of density structures

Figure 3 shows the density map of all models. Although the disk size is
different among the simulation, the density structure of inner $\sim 20$ AU
region of panels a (ModelA100Q25), b (ModelA100Q35), c (ModelA5Q25), and f
(ModelZeta18) are very similar. In these simulations, the density structures
are consistent with the analytical solution.

On the other hand, figure 3 d (ModelA5Q35) shows the formation of a very small
disk of $\lesssim 10{\rm AU}$ and bubble-like structures around it. This
bubble-like structure is created by magnetic interchange instability
(Krasnopolsky et al., 2012). The large amounts of small dust grains make
ambipolar diffusion (and Ohmic diffusion) ineffective in the high-density
region and leads to the development of interchange instability as a
redistribution mechanism of magnetic flux.

Figure 3 e (ModelA100Fixed) in which the dust growth is artificially neglected
shows that the disk is relatively compact, dense and massive. This massive
disk is consistent with previous theoretical studies; however such massive
disk seems to be inconsistent with the observations (Tsukamoto et al., 2022).

####  3.3.2 Universality of the disk structure

Figure 4 shows the azimuthally averaged radial profiles of the all models. The
dashed lines show the power laws of our analytical solutions (in this figure,
we plot the power laws just as a reference because parameters such as central
star mass or mass accretion rate differ among the models).

Of the five models that consider dust growth, four models have an inner region
of $\sim 20{\rm AU}$ consistent with the analytical solution (red, green,
black, and orange lines). If we regard the disk radius as the radius at which
the density distribution deviates from the power law of the analytical
solution, ModelZeta18 has the largest disk size, and ModelA100Q35,
ModelA100Q25, and ModelA5Q25 have smaller disk sizes in this order (see also
figure 5). This difference may be due to the difference in resistivity in the
regions where dust grains have not grown (outer region of the disk and
envelope). Despite the difference in disk size, the universality of the disk
structure in the inner region is noteworthy.

Note that the disk structure of ModelZeta18 (green lines) is more consistent
with the power laws of the analytical solutions than the fiducial model. For
example, $\eta_{A}$ has a positive power exponent. This is because this model
is more evolved than the fiducial model and approximation on $\eta_{A}$ in the
analytical solution is better validated.

Figure 4 shows that the disk of ModelA5Q25 (yellow) is very small and has a
different structure from the other four models. A large amount of small dust
in this model makes ambipolar diffusion ineffective from the beginning of disk
formation in the high-density region. The disk rapidly shrinks before the dust
grows and $\eta_{A}$ is well described by the analytical form of Shu (Shu,
1983). Thus, the disk evolution of the model varies.

It would be instructive to see the differences in disk structure between the
model ignoring the dust growth (ModelA100Fixed; magenta) and those that well
described by the analytical solutions. In the model without the dust growth,
the density and $\eta_{A}$ are large, and the magnetic field and radial
velocity $v_{r}$ are small. This is because in the absence of the dust growth,
the ambipolar diffusion in the disk is extremely effective. It suppresses
magnetic braking in the disk, resulting in smaller $v_{r}$ and an increased
density due to gas accumulation in the disk. Furthermore, the magnetic flux is
extracted from the disk by the ambipolar diffusion causing the low value of
the disk magnetic field. The density of the disk in ModelA100Fixed is large
and Toomre $Q$ parameter is $O(1)$.

The disk density profile in ModelA100Fixed is approximately $\rho_{g}\propto
r^{-1}$, meaning that the surface density $\Sigma_{g}$ is approximately
$\Sigma_{g}\propto r^{1/4}$ (here we assume $\rho_{g}=\Sigma_{g}/H$ and
$H=c_{s}/\Omega\propto r^{5/4}$). While this shallow profile is consistent
with previous theoretical studies (Tomida et al., 2017), it is inconsistent
with the disk observation ($\Sigma_{g}\propto r^{-0.9}$)(Williams & Cieza,
2011). Our analytical solution and the simulation results with dust growth, on
the other hand, has deeper power law of $\rho_{g}\propto r^{-9/5}$ and is
therefore more consistent with the observations.

Figure 3:  Gas density map of (a) ModelA100Q25, (b) ModelA100Q35, (c)
ModelA5Q25, (d) ModelA5Q35, (e) ModelA100Fixed, (f) ModelZeta18 on the $x$-$y$
plane. The time after protostar formation are shown in upper left. The red and
orange arrows show the velocity field of the dust and gas, respectively, on
the $x$-$y$ plane. The black lines indicate the density contour. The contour
levels are $\rho_{g}=10^{-15},10^{-14.5},\cdots,10^{-10}~{}{\rm g~{}cm}^{-3}$.

Figure 4:  The azimuthally averaged radial profile of (a) gas density, (b)
vertical magnetic field, (c) radial velocity, and (d) ambipolar resistivity
($\eta_{A}$). The red, orange, black, yellow, magenta, and green lines show
the results of ModelA100Q25, ModelA100Q35, ModelA5Q25, ModelA5Q35,
ModelA100Fixed, ModelZeta18, respectively. The thick lines show the models
that are in good agreement with the power laws of the analytical solution. The
time of each line is the same as in figure 3. The black dashed lines indicate
the power laws of the analytical solutions.

###  3.4 Time evolution

In this section, we examine the time evolution of disk angular momentum (and
radius), disk mass, typical dust size and dust abundance in disks. We also
investigate the mass ejection rate by outflows from the disks.

####  3.4.1 Time evolution of disk size

Figure 5 shows the time evolution of centrifugal radius and the angular
momentum of the disk. The angular momentum of disk $J(\rho_{\rm disk})$ is
calculated as

$\displaystyle J(\rho_{\rm disk})\equiv\left|\int_{\rho_{g}>\rho_{\rm
disk}}\rho_{g}(\mathbf{r}\times\mathbf{v})dV\right|.$ (13)

For the density threshold of the disk, we choose $\rho_{\rm
disk}=10^{-13}~{}{\rm g~{}cm}^{-3}$. The centrifugal radius is calculated as

$\displaystyle r_{\rm disk}\equiv r_{\rm cent}=\frac{\bar{j}(\rho_{\rm
disk})^{2}}{GM_{\rm star}}.$ (14)

Here $\bar{j}(\rho_{\rm disk})=J(\rho_{\rm disk})/M_{\rm disk}$, where $M_{\rm
disk}$ denotes the enclosed gas mass within the region $\rho_{g}>\rho_{\rm
disk}$. Comparing the radius of the region with a density of
$\rho_{g}>\rho_{\rm disk}$ (figure 4) with the centrifugal radius, the former
is about 2 times larger than the latter owing to radial density distribution
and temporal density oscillation by the non-axisymmetric structures. In this
study, we consider the centrifugal radius as an estimate of the disk size.

The left panel of figure 5 shows the results of the models with $a_{\rm
min}=100~{}{\rm nm}$. In ModelA100Q25 (red line), the angular momentum of the
disk continues to increase until $t_{*}\sim 7\times 10^{3}$ yr, and then it
shows a sharp decrease. This is because the dust grains grows to $a_{\rm
d}\gtrsim 10{\rm\mu}{\rm m}$ (figure 7), which causes a sudden decrease in
magnetic resistivity and the extraction of angular momentum by magnetic field.
Interestingly, although the angular momentum has decreased by a factor of
$\sim 1/3$ from $t_{*}\sim 7\times 10^{3}$ yr to $t_{*}\sim 9\times 10^{3}$
yr, the centrifugal radius has decreased by only a factor of $\sim 1/2$. This
indicates that the disk mass also decreases rapidly during this period (figure
6). For ModelA100Q35, no such rapid decrease of angular momentum is observed.
This is because $a_{\rm min}$ is also responsible for the total dust surface
area, and thus the decrease in magnetic resistivity is not so drastic
(Tsukamoto & Okuzumi, 2022). As exhibited by ModelZeta-18 (green line), the
low cosmic-ray ionization rate can contribute to maintaining the angular
momentum. This is in agreement with previous studies (Wurster et al., 2018;
Kuffmeier et al., 2020; Kobayashi et al., 2023). The decrease in the disk size
of ModelA100Fixed is caused by the pseudo-disk warp and associated inward
magnetic flux drag (Tsukamoto et al., 2020). In this model, the angular
momentum decreases less than a factor of two, which is not significant
compared to ModelA100Q25.

The right panel of figure 5 shows the results of the models with $a_{\rm
min}=5~{}{\rm nm}$. Interestingly, the relationship between the disk size and
power exponent $q$ is different from the models with $a_{\rm min}=100{\rm
nm}$. the disk size of the model with $q=3.5$ (ModelA5Q35) is significantly
smaller than that in the model with $q=2.5$ (ModelA5Q25). This is because when
there is a large amount of small dust ($\sim{\rm nm}$), dust grains are
responsible for the conductivity and reduce the magnetic resistivity. This
allows magnetic braking to work more effectively in ModelA5Q35 and reduce the
disk size.

####  3.4.2 Time evolution of the disk mass

Figure 6 shows the time evolution of disk mass $M_{\rm disk}$ (solid),
protostar mass $M_{\rm star}$ (dashed), and total mass $M_{\rm star}+M_{\rm
disk}$ (dotted). In ModelA100Q25 (red), the disk mass continues to increase
and reaches $\sim 0.15\thinspace M_{\odot}$ at $t_{*}\sim 7\times 10^{3}$ yr.
Then, it drops sharply to $M_{\rm disk}\sim 0.03\thinspace M_{\odot}$ at
$t_{*}\sim 9\times 10^{3}$ yr. Meanwhile, mass accretion onto protostars is
enhanced and the protostellar mass rapidly increases from $0.1\thinspace
M_{\odot}$ to 0.2 $\thinspace M_{\odot}$ in this period, giving a mass
accretion rate of $\gtrsim 10^{-5}\thinspace M_{\odot}~{}{\rm yr}^{-1}$ for
this model. In ModelZeta18 (green line) and ModelA5Q25 (black line), it
appears that the disk mass also begins to decrease towards the end of the
simulations. However, the time at which the decrease begins is later than in
ModelA100Q25.

In ModelA100Fixed, the disk radius decreases in $t_{*}\gtrsim 6\times 10^{3}$
yr (figure 5) but the disk mass does not significantly change. This suggests
that the disk evolution without dust growth is different from the models that
include dust growth and that are consistent consistent with the analytic
solution. In ModelA5Q35, the mass (and radius) evolution differs from the
other models. This is due to inefficient ambipolar diffusion since disk
formation. Thus, the situation close to the ideal MHD is realized.

####  3.4.3 time evolution of dust size and dust abundance in the disk

Figure 7 shows the time evolution of the dust-to-gas mass ratio and mean dust
size of the disks. The dust mass and mean dust size of the disk is calculated
as

$\displaystyle\bar{a}_{\rm disk}\equiv\frac{1}{M_{\rm
disk}}\int_{\rho_{g}>\rho_{\rm disk}}\rho_{g}a_{d}dV,$ (15)

and

$\displaystyle M_{\rm dust,disk}\equiv\int_{\rho_{g}>\rho_{\rm
disk}}\rho_{d}dV,$ (16)

respectively.

The figure shows that the increase in the dust-to-gas mass ratio occurs later
in the simulation. This increase begins when the average dust size in the disk
exceeds $\sim 100{\rm\mu}{\rm m}$. This is due to the “ash-fall phenomenon”
proposed in our previous study (Tsukamoto et al., 2021b). The largest increase
is observed in ModelA100Q25, where the dust-to-gas mass ratio increases to
1.04% at the end of the simulation. Some readers may think that this value is
small and irrelevant. However, we only considered $\sim 10^{3}$ yr after the
ratio started to increase. If this event continues for, for instance, $10^{5}$
yr (i.e., during the Class 0/I phase), it can cause a significant increase in
the dust abundance.

The increase of the dust-to-gas mass ratio is slower in ModelA100Q35,
ModelA5Q25, and ModelZeta18, compared to ModelA100Q25 owing to the lower
outflow activity (i.e., weaker coupling between the magnetic field and gas in
the upper layers of the disk; figure 8) in these models.

####  3.4.4 Time evolution of the mass ejection rate by outflow

Figure 8 shows the time evolution of the mass ejection rate due to the
molecular outflow. The mass ejection rate is calculated as

$\displaystyle\dot{M}_{\rm out}\equiv\int_{r=100{\rm AU}}\rho v^{+}_{r}dS,$
(17)

where $v_{r}^{+}\equiv\max(\mathbf{v}\cdot\mathbf{r},0)$ denotes the positive
radial velocity, and we perform the surface integral on a sphere with radius
of $100{\rm AU}$.

The mass ejection rate is highly variable and has a peak value of $\gtrsim
10^{-5}\thinspace M_{\odot}~{}{\rm yr}^{-1}$. This is comparable to the mass
accretion rate in the disk. Comparing the mass ejection rate with the epoch
when the dust-to-gas mass ratio begins to increase (at $t_{*}\sim 8\times
10^{3}$ yr; figure 7), a mass ejection rate of $10^{-6}$ to $10^{-5}\thinspace
M_{\odot}~{}{\rm yr}^{-1}$ and the mean dust size of $\sim 100{\rm\mu}{\rm m}$
are required for the ”ash-fall” phenomenon to happen, in which the disk gas is
selectively ejected into interstellar space by the outflow and the dust grains
are resupplied to the disk.

Figure 5:  Time evolution of centrifugal radius (the solid lines and the left
axis) and total angular momentum of disk (the dashed lines and the right
axis). Red, orange, black, yellow, magenta, and green lines show the results
of ModelA100Q25,ModelA100Q35, ModelA5Q25, ModelA5Q35, ModelA100Fixed,
ModelZeta18, respectively. The horizontal axis shows the time after the
protostar formation.

Figure 6:  Time evolution of the mass of the disks (the solid lines), the
protostars (the dashed lines), and the total mass (the dotted lines). Red,
orange, black, yellow, magenta, and green lines show the results of
ModelA100Q25, ModelA100Q35, ModelA5Q25, ModelA5Q35, ModelA100Fixed,
ModelZeta18, respectively. The horizontal axis shows the time after the
protostar formation.

Figure 7:  Time evolution of the dust-to-gas mass ratio (the solid lines and
left axis) and mean dust size in the disk (the dashed lines and right axis).
Red, orange, black, yellow, magenta, and green lines show the results of
ModelA100Q25, ModelA100Q35, ModelA5Q25, ModelA5Q35, ModelA100Fixed,
ModelZeta18, respectively. The horizontal axis shows the time after the
protostar formation.

Figure 8:  Time evolution of the mass ejection rate. Red, orange, black,
yellow, magenta, and green lines show the results of
ModelA100Q25,ModelA100Q35, ModelA5Q25, ModelA5Q35, ModelA100Fixed,
ModelZeta18, respectively. The horizontal axis shows the time after the
protostar formation.

##  4 Discussion

###  4.1 Universality of the disk structure in the disk with grown dust grains

In this study, we propose the new disk evolutionary picture “co-evolution of
dust grains and protoplanetary disks” based on non-ideal dust-gas two-fluid
MHD simulations considering dust growth. The dust growth changes the gas-phase
ionization degree, magnetic resistivity, and evolution of the disk. In the co-
evolution process, the microscopic dust grains of micrometer size couple with
the macroscopic disks of $100$ AU size and they co-evolve. The size scale
difference between two objects is $10^{19}$, which is astounding compared to
the well-known co-evolution of super massive black holes and galaxies (size
scale difference is $\sim 10^{10}$).

Furthermore, once the dust grains grow sufficiently, the structure of
protoplanetary disks is well described by the non-trivial power laws, which we
analytically derive in equations (A) to (A) (figure 2 and 4). From the
assumptions adopted in the analytical solutions, we conclude that the disk
structures will emerge when

  1. 1.

Dust grains grow sufficiently and adsorption of charged particles by the dust
grains becomes negligible,

  2. 2.

The toroidal magnetic field in the disk is determined by the balance between
vertical shear (of the order of $(H/r)^{2}$) and ambipolar diffusion, and

  3. 3.

Angular momentum transport mechanisms other than magnetic braking (such as
turbulent viscosity) are negligible.

We believe that the discovery of this new disk structure is a theoretical
breakthrough for star and planet formation theory. The disk structure is
determined only by observable parameters such as the central star mass, mass
accretion rate, disk temperature, and cosmic-ray ionization rate, without
including difficult-to-determine parameters such as the viscous parameter
$\alpha$. Using the analytcal soluiton, we can study the planet formation
process in the realistic disk and evolution of the magnetic flux during
protostellar evolution. In the future, we will discuss the broad implications
of this disk model for the formation and evolution of protostars and planets.

###  4.2 Assumptions employed in the dust growth model and their uncertainty

Our simulations make several simplifications to the dust growth and dust size
distribution. The largest simplification is the representative size
approximation for dust growth in which we assume that the representative size
corresponds to the peak dust size of the mass distribution (note that the peak
of the dust mass distribution corresponds to the maximum dust size if $q<4$).
Our approximated equation for dust growth can be derived from the coagulation
equation (for derivation, see (Sato et al., 2016)). In this study, we solve
the evolution of the representative size and regard it to be the maximum dust
size $a_{\rm max}$, and set the minimum dust size and power as parameters.
Moreover we implicitly assume that the size distribution can be described by a
single power law. More realistically, the time evolution of the dust size
distribution should be considered, and the validity of the simplifications
employed in this study should be investigated in future more realistic
studies.

However, our claim, ”the disk structure converges to the analytical solution
once the dust has grown sufficiently”, remains valid regardless of the
specific details of the dust distribution and dust growth model. This is
because the essential physics required for the disk to converge to the
analytical solutions is that the total dust surface area becomes sufficiently
small and adsorption of charged particles by the dust grains becomes
negligible. Under these conditions, the ambipolar resistivity converges to
$\eta_{A}\propto B^{2}\rho^{-3/2}$ and we obtain the analytical solutions. In
this sense, our results are universal, regardless of the details of the dust
growth model.

###  4.3 Maximum dust size lifted up by the outflow

Figure 7 as well as our previous study (Tsukamoto et al., 2021b) show that the
dust-to-gas mass ratio increases when the mean dust size in the disk reaches
$a_{d}\gtrsim 100{\rm\mu}{\rm m}$. This is caused by the selective fall of
dust grains from the dust-gas mixture lifted up by the outflow. The
simulations show that, once the dust grows to $\sim 100{\rm\mu}{\rm m}$, ”ash-
fall” phenomenon occurs, in which the dust grains and gas are decoupled in the
outflow, and only the dust falls back into the disk. Hence, the minimum dust
size for the protostellar ash-fall is $\sim 100{\rm\mu}{\rm m}$. Then, how
large is the maximum dust size that can be lifted by the outflow or the
maximum dust size of the ”falling ash”? This is particularly important in
explaining the recent observations of the presence of grown dust in the
envelope (Kwon et al., 2009; Galametz et al., 2019; Valdivia et al., 2019).

The maximum dust size lifted up by the outflow can be estimated from the
following considerations. For the dust grain to be lifted up by the outflow,
the dust grains must couple to the gas at the outflow driving point (or root).
Thus, the stopping time of the dust grains should be less than the orbital
period at the root (otherwise, outflow driving causes dust grains to remain in
the disk and only the gas is ejected).

To estimate the stopping time, we need the density at the outflow root. As
shown in figure 8 and by the observations (Wu et al., 2004), the mass ejection
rate of the outflow in young protostars (e.g., their age is $t\lesssim 10^{5}$
yr) is in the range $\dot{M}_{\rm out}=10^{-6}-10^{-5}\thinspace
M_{\odot}~{}{\rm yr}^{-1}$. Thus, by assuming that the outflow velocity is
comparable to the orbital velocity at the radius of the root (Kudoh & Shibata,
1997), the density at the root of the outflow $\rho_{\rm dp}$ can be estimated
as,

$\displaystyle\rho_{\rm dp}$ $\displaystyle=\frac{\dot{M}_{\rm out}}{v_{\rm
out}\pi r_{\rm disk}^{2}}$ (18) $\displaystyle=1.9\times 10^{-15}$
$\displaystyle\left(\frac{\dot{M}_{\rm out}}{10^{-5}\thinspace
M_{\odot}~{}{\rm yr}^{-1}}\right)\left(\frac{M_{*}}{\rm 0.2\thinspace
M_{\odot}}\right)^{-1/2}\left(\frac{r_{\rm disk}}{50{\rm
AU}}\right)^{-3/2}~{}{\rm g~{}cm}^{-3},$

where we assume that $v_{\rm out}=\sqrt{GM_{*}/r}$ is the Keplerian velocity
at disk outer edge $r_{\rm disk}$. Hence, the stopping time at the root is
estimated as

$\displaystyle t_{\rm stop}$ $\displaystyle=\frac{\rho_{\rm
mat}a_{d}}{\rho_{\rm dp}\sqrt{8/\pi}c_{s}}$ (19) $\displaystyle=3.3\times
10^{2}\left(\frac{\rho_{\rm mat}}{2{\rm~{}g}}\right)\left(\frac{a_{d}}{5{\rm
mm}}\right)$ $\displaystyle\left(\frac{M_{*}}{\rm 0.2\thinspace
M_{\odot}}\right)^{1/2}\left(\frac{r_{\rm disk}}{50{\rm
AU}}\right)^{12/7}\left(\frac{\dot{M}_{\rm out}}{10^{-5}\thinspace
M_{\odot}~{}{\rm yr}^{-1}}\right)^{-1}{\rm yr},$

where the sound velocity and temperature is assumed to be $c_{s}=190(T/10{\rm
K})^{1/2}~{}{\rm m}~{}{\rm s}^{-1}$ and $T=150(r/{\rm AU})^{-3/7}{\rm K}$,
respectively. Then, the ratio of the stopping time to the orbital period at
$r_{\rm disk}$ is calculated as,

$\displaystyle\frac{t_{\rm stop}}{t_{\rm orb}}$
$\displaystyle=0.41\left(\frac{\rho_{\rm mat}}{2{\rm~{}g}}\right)$ (20)
$\displaystyle\left(\frac{a_{d}}{5{\rm mm}}\right)\left(\frac{M_{*}}{\rm
0.2\thinspace M_{\odot}}\right)$ $\displaystyle\left(\frac{r_{\rm
disk}}{50{\rm AU}}\right)^{3/14}\left(\frac{\dot{M}_{\rm
out}}{10^{-5}\thinspace M_{\odot}~{}{\rm yr}^{-1}}\right)^{-1},$

or $t_{\rm stop}\sim t_{\rm orb}$ is realized when

$\displaystyle a_{d}$ $\displaystyle\sim 1.2\left(\frac{\rho_{\rm
mat}}{2{\rm~{}g}}\right)^{-1}\left(\frac{M_{*}}{\rm 0.2\thinspace
M_{\odot}}\right)^{-1}$ (21) $\displaystyle\left(\frac{r_{\rm disk}}{50{\rm
AU}}\right)^{-3/14}\left(\frac{\dot{M}_{\rm out}}{10^{-5}\thinspace
M_{\odot}~{}{\rm yr}^{-1}}\right){\rm cm}.$

This indicates that the dust size of at maximum $a_{d}\sim 1{\rm cm}$ can be
entrained by the outflow with $\dot{M}_{\rm out}\sim 10^{-5}\thinspace
M_{\odot}~{}{\rm yr}^{-1}$ from the disk with a size of $\sim 50$ AU. This
size is larger than the wavelength of sub-millimeter observations such as with
ALMA and may cause the decrease of the spectral index of dust opacity in the
outflow and the envelope. Thus, “ash-fall” can explain the presence of grown
dust in the envelope suggested by the observations.

## Acknowledgments

We thank Dr. Shinsuke Takasao and Mr. Ryoya Yamamoto for the fruitful
discussion. The computations were performed on the Cray XC50 system at CfCA of
NAOJ. This work is supported by JSPS KAKENHI grant number 18H05437, 18K13581,
18K03703.

##  Appendix A Analytic soluiton of steady state disks with magnetic braking
and ambipolar diffusion

In this section, we derive the power laws of steady-state circumstellar disks
which is determined by the angular momentum removal by magnetic braking and
magnetic field structure determined by the balance between gas advection and
ambipolar diffusion. We start from MHD equation with ambipolar diffusion,

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})$
$\displaystyle=$ $\displaystyle 0,$ (22)
$\displaystyle\rho\left(\frac{\partial\mathbf{v}}{\partial
t}+\mathbf{v}\cdot\nabla\mathbf{v}\right)$ $\displaystyle=$
$\displaystyle-\rho\nabla\Phi-\nabla
p+\frac{1}{4\pi}(\nabla\times\mathbf{B})\times\mathbf{B},$
$\displaystyle\frac{\partial\mathbf{B}}{\partial t}$ $\displaystyle=$
$\displaystyle\nabla\times\left(\mathbf{v}\times\mathbf{B}\right)$
$\displaystyle-$
$\displaystyle\nabla\times\left(\frac{\eta_{A}}{|\mathbf{B}|^{2}}((\nabla\times\mathbf{B})\times\mathbf{B})\times\mathbf{B}\right),$
$\displaystyle\nabla\cdot\mathbf{B}$ $\displaystyle=$ $\displaystyle 0.$ (25)

In this appendix, $\rho$ denotes the gas density, $\mathbf{v}$ denotes the gas
velocity, $p$ denotes the gas pressure, $\mathbf{B}$ denotes the magnetic
field. Hereafter, we assume the steady state i.e., $\frac{\partial}{\partial
t}=0$.

Because we focus the structure of the circumstellar disk, we make the
assumptions below following Guilet & Ogilvie (2012, 2013). We introduce a
small dimensionless parameter $\epsilon=O(H/r)$, where $H$ denotes the gas
scale height of the disk. When the disk self-gravity is negligible, the
gravitational potential is expanded as

$\displaystyle\Phi(r,z)=\Phi_{0}(r)+\frac{1}{2}\Phi_{2}(r)z^{2}+O\left(\left(\frac{z}{r}\right)^{4}\right),$
$\displaystyle\therefore\Phi(r,\zeta)=\Phi_{0}(r)+\epsilon^{2}\frac{1}{2}\Phi_{2}(r)\zeta^{2}+O(\epsilon^{4}),$
(26)

where we introduces the rescaled vertical coordinate $\zeta=\epsilon^{-1}z$,
and $\Phi_{0}(r)=-(GM/r)$ and $\Phi_{2}(r)=(GM/r^{3})$. This dependence of the
gravitational potential on $\epsilon$ is the guiding principle that determines
the order of other physical quantities.

We assume that the vertical gravity is balanced by thermal pressure at the
leading order. Thus, the scaling of pressure can be assumed to be

$\displaystyle
p(r,\zeta)=\epsilon^{2}p_{2}(r,\zeta)+\epsilon^{4}p_{4}(r,\zeta)+O(\epsilon^{6}).$
(27)

We assume the scaling of the density to be

$\displaystyle\rho(r,\zeta)=\rho_{0}(r,\zeta)+\epsilon^{2}\rho_{2}(r,\zeta)+O(\epsilon^{4}).$
(28)

Then, the sound velocity $c_{s}$ scales as

$\displaystyle c_{s}(r)=\epsilon
c_{s,1}(r)+\epsilon^{3}c_{s,3}(r)+O(\epsilon^{5}).$ (29)

Here, we assume the disk to be vertically isothermal.

The leading order of the radial velocity is assumed to be $\epsilon c_{s}$,
and thus,

$\displaystyle
v_{r}(r,\zeta)=\epsilon^{2}v_{r,2}(r,\zeta)+\epsilon^{4}v_{r,4}(r,\zeta)+O(\epsilon^{6}).$
(31)

The leading order of the azimuthal velocity is assume to be balanced with the
leading order of $\Phi(r,\zeta)$ and hence,

$\displaystyle
v_{\phi}(r,\zeta)=r\Omega_{0}(r)+\epsilon^{2}v_{\phi,2}(r,\zeta)+\epsilon^{4}v_{\phi,4}(r,\zeta)+O(\epsilon^{6}).$
(32)

The leading order of the vertical velocity is assumed to be smaller than
$v_{r}$,

$\displaystyle
v_{z}(r,\zeta)=\epsilon^{3}v_{z,3}(r,\zeta)+\epsilon^{5}v_{r,5}(r,\zeta)+O(\epsilon^{7}).$
(33)

For the magnetic field, we assume that the vertical magnetic field is the
dominant component and the scaling of magnetic field is assumed to be

$\displaystyle B_{r}(r,\zeta)$ $\displaystyle=$
$\displaystyle\epsilon^{2}B_{r,2}(r,\zeta)+\epsilon^{4}B_{r,4}(r,\zeta)+O(\epsilon^{6}),$
(34) $\displaystyle B_{\phi}(r,\zeta)$ $\displaystyle=$
$\displaystyle\epsilon^{2}B_{\phi,2}(r,\zeta)+\epsilon^{4}B_{\phi,4}(r,\zeta)+O(\epsilon^{6}),$
(35) $\displaystyle B_{z}(r,\zeta)$ $\displaystyle=$ $\displaystyle\epsilon
B_{z,1}(r)+\epsilon^{3}B_{z,3}(r,\zeta)+O(\epsilon^{5}),$ (36)

where the leading term of $B_{z}(r,\zeta)$ does not depend on $\zeta$ because
the leading order of the divergence free condition gives,

$\displaystyle\partial_{\zeta}B_{z,1}=0.$ (37)

The underlying assumption that leads to this ordering is that the leading
order of Alfven velocity ($\propto B_{z}/\sqrt{\rho_{g}}=O(\epsilon)$) is the
order of the sound velocity.

Finally, we assume that $\eta_{A}$ is given as

$\displaystyle\eta_{A}=\frac{\mathbf{B}^{2}}{C\gamma\rho^{3/2}}=\epsilon^{2}\eta_{A,2}+O(\epsilon^{4}),$
(38)

whose leading order is $\epsilon^{2}$ because of the leading order of
$B_{z}(r,\zeta)$ and $\rho(r,\zeta)$ are $O(\epsilon)$ and $O(1)$,
respectively. This implies that the dust grains have sufficiently grow and
their adsorption of charged particles is negligible.

The radial components of the equation of motion at the leading order,

$\displaystyle-\rho_{0}r\Omega_{0}^{2}=-\rho_{0}\partial_{r}\Phi_{0},$ (39)

gives

$\displaystyle\Omega_{0}=\sqrt{\frac{GM}{r^{3}}}.$ (40)

The vertical components of equation of motion at the leading order

$\displaystyle\rho_{0}\Phi_{2}\zeta=-\partial_{\zeta}p_{2},$ (41)

leads

$\displaystyle\rho_{0}=\frac{\Sigma_{0}}{\sqrt{2\pi}H_{1}}\exp\left(-\frac{\zeta^{2}}{2H_{1}^{2}}\right)\equiv\tilde{\rho}\exp\left(-\frac{\zeta^{2}}{2H_{1}^{2}}\right),$
(42)

where $H_{1}(r)=c_{s,1}/\Phi_{2}^{1/2}=c_{s,1}/\Omega_{0}$, where $c_{s,1}$ is
the vertically isothermal sound velocity defined by
$p_{2}=c_{s,1}^{2}\rho_{0}$ and $H_{1}(r)$ is related to the scale height as
$H=\epsilon H_{1}+O(\epsilon^{3})$. $\Sigma_{0}$ is given as
$\Sigma_{0}=\int_{-\infty}^{\infty}\rho_{0}d\zeta$.

The second order of the radial and azimuthal components of the equation of
motion and of the induction equation are written as

$\displaystyle-2\rho_{0}\Omega_{0}v_{\phi,2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\rho_{0}\partial_{r}\Phi_{2}\zeta^{2}$ (43)
$\displaystyle-$
$\displaystyle\partial_{r}\left(p_{2}+\frac{B_{z,1}^{2}}{8\pi}\right)+\frac{B_{z,1}}{4\pi}\partial_{\zeta}B_{r,2},$
$\displaystyle\rho_{0}v_{r,2}\frac{1}{r}\partial_{r}(r^{2}\Omega_{0})$
$\displaystyle=$
$\displaystyle\frac{B_{z,1}}{4\pi}\partial_{\zeta}B_{\phi,2},$ (44)
$\displaystyle 0=B_{z,1}\partial_{\zeta}v_{r,2}$ $\displaystyle+$
$\displaystyle\partial_{\zeta}[\eta_{A,2}(\partial_{\zeta}B_{r,2}-(\partial_{r}B_{z,1})].$
(45) $\displaystyle 0=B_{r,2}r\partial_{r}\Omega_{0}$ $\displaystyle+$
$\displaystyle
B_{z,1}\partial_{\zeta}v_{\phi,2}+\partial_{\zeta}(\eta_{A,2}\partial_{\zeta}B_{\phi,2})$
(46)

These equations correspond to equation (28) to (31) of Guilet & Ogilvie
(2012), if we assume viscous parameter $\alpha$ to be $0$.

Then, we rescale the vertical coordinate with

$\displaystyle\hat{z}\equiv\frac{\zeta}{H_{1}}=\frac{z}{H}.$ (47)

To evaluate the single power law for each physical quantity, we need to
further simplify the equations (43) to (46). Here, we assume that the radial
thermal pressure gradient is much larger than the magnetic pressure gradient,
and the terms with $B_{r,2}$ can be neglected.

Then equations (43) to (46) are rewritten as

$\displaystyle-2\rho_{0}\Omega_{0}v_{\phi,2}$ $\displaystyle=$
$\displaystyle\frac{3}{2r}\rho_{0}\Omega_{0}^{2}H_{1}^{2}\hat{z}^{2}-\partial_{r}(\rho_{0}c_{s,1}^{2}),$
(48) $\displaystyle\frac{1}{2}\rho_{0}v_{r,2}\Omega_{0}$ $\displaystyle=$
$\displaystyle\frac{B_{z,1}}{4\pi H_{1}}\partial_{\hat{z}}B_{\phi,2},$ (49)
$\displaystyle B_{z,1}\partial_{\hat{z}}v_{r,2}$ $\displaystyle=$
$\displaystyle(\partial_{\hat{z}}\eta_{A,2})(\partial_{r}B_{z,1}),$ (50)
$\displaystyle B_{z,1}\partial_{\hat{z}}v_{\phi,2}$ $\displaystyle=$
$\displaystyle\frac{1}{H_{1}}\partial_{\hat{z}}(\eta_{A,2}\partial_{\hat{z}}B_{\phi,2}).$
(51)

Furthermore, from the conservation of the mass, we have

$\displaystyle-2\pi r\int^{H}_{-H}\rho
v_{r}dz=\dot{M}\equiv\epsilon^{3}\dot{M}_{3}+O(\epsilon^{5}).$ (52)

where we approximated the integral range from $-H$ to $H$ instead of from
$-\infty$ to $\infty$ because of the vertical expansion below. $\dot{M}$ is
the mass accretion rate within the disk which is assumed to be constant and
$\dot{M}_{3}$ is its leading term.

Using $\hat{z}$, we consider the vertical expansion. From the analytical form
of $\eta_{A}$, it is expanded as

$\displaystyle\eta_{A,2}$ $\displaystyle=$
$\displaystyle\frac{B_{z,1}(r)^{2}}{C\gamma\rho_{0}(r,\hat{z})^{3/2}}=B_{z,1}^{2}\tilde{\rho}^{-3/2}\exp(\frac{3}{4}\hat{z}^{2})$
(53) $\displaystyle=$
$\displaystyle(C\gamma)^{-1}B_{z,1}^{2}\tilde{\rho}^{-3/2}\left(1+\frac{3}{4}\hat{z}^{2}+O(\hat{z}^{4})\right).$

Because the azimuthal component of the magnetic field should be odd with
respect to the midplane, we assume a vertical dependence of the magnetic field
as

$\displaystyle
B_{\phi,1}=B_{\phi,\hat{z}1}\hat{z}+B_{\phi,\hat{z}3}\hat{z}^{3}+O(\hat{z}^{5}).$
(54)

Note that $B_{r,1}$ has disappeared from the equations and cannot be
determined from our assumptions above.

The radial and azimuthal components of the velocity should be even with
respect to the midplane. Thus, we assume the vertical dependence of the
velocity as

$\displaystyle
v_{r,2}=v_{r,\hat{z}0}+v_{r,\hat{z}2}\hat{z}^{2}+O(\hat{z}^{4}),$ (55)
$\displaystyle
v_{\phi,2}=v_{\phi,\hat{z}0}+v_{\phi,\hat{z}2}\hat{z}^{2}+O(\hat{z}^{4}).$
(56)

Then we assume the power law for the vertical magnetic field $B_{z}$ and
midplane density $\rho_{\rm mid}\equiv\rho(r,z=0)$ with respect to the radius
as

$\displaystyle B_{z}(r)=B_{\rm z,ref}\left(\frac{r}{r_{\rm
ref}}\right)^{D_{B_{z}}},$ (57) $\displaystyle\rho_{\rm mid}(r)=\rho_{\rm
mid,ref}\left(\frac{r}{r_{\rm ref}}\right)^{D_{\rho_{g}}}.$ (58)

To be consistent with our simulations, we assume the polytropic relation for
the sound velocity

$\displaystyle c_{s}(r)=c_{\rm s,ref}\left(\frac{\rho_{\rm
mid}(r)}{\rho_{c}}\right)^{1/3}.$ (59)

By substituting these power laws and taking the leading terms with respect to
the $\epsilon$ of equation (52), we obtain the solutions of the equations (48)
to (52),

$\displaystyle v_{\phi,\hat{z}0}$ $\displaystyle=$ $\displaystyle
2v_{\phi,\hat{z}2},$ (60) $\displaystyle v_{r,\hat{z}0}$ $\displaystyle=$
$\displaystyle 2v_{r,\hat{z}2},$ (61) $\displaystyle B_{\phi,\hat{z}3}$
$\displaystyle=$ $\displaystyle 0,$ (62)

and,

$\displaystyle\rho_{\rm mid}(r)$
$\displaystyle=2\left(\frac{2~{}5^{2}}{3^{4}~{}19^{2}\pi^{2}}\right)^{1/5}\dot{M}^{2/5}\rho_{c}^{2/5}\Omega_{0}^{6/5}$
$\displaystyle(C\gamma)^{-2/5}c_{\rm s,ref}^{6/5}\left(\frac{r}{r_{\rm
ref}}\right)^{-\frac{9}{5}},$ (63) $\displaystyle B_{z}(r)$
$\displaystyle=5\left(\frac{2^{3}~{}5}{3^{17}~{}19^{16}\pi}\right)^{1/30}\dot{M}^{8/15}\rho_{c}^{1/5}\Omega_{0}^{3/5}$
$\displaystyle(C\gamma)^{7/15}c_{\rm s,ref}^{-3/5}\left(\frac{r}{r_{\rm
ref}}\right)^{-\frac{9}{10}},$ (64) $\displaystyle v_{\phi}(r,z=0)-v_{K}$
$\displaystyle=\left(\frac{5^{4}~{}3^{7}}{2^{3}~{}19^{4}\pi^{4}}\right)^{\frac{1}{15}}\dot{M}^{4/15}\rho_{c}^{-2/5}\Omega_{0}^{4/5}$
$\displaystyle(C\gamma)^{-4/15}c_{\rm s,ref}^{1/5}\left(\frac{H_{\rm
ref}}{r_{\rm ref}}\right)\left(\frac{r}{r_{\rm ref}}\right)^{-\frac{7}{10}},$
(65) $\displaystyle v_{r}(r,z=0)$
$\displaystyle=\frac{3}{2}\left(\frac{3~{}5^{7}}{2^{9}~{}19^{7}\pi^{7}}\right)^{\frac{1}{15}}\dot{M}^{7/15}\rho_{c}^{-1/5}\Omega_{0}^{2/5}(C\gamma)^{8/15}$
$\displaystyle c_{\rm s,ref}^{-2/5}\left(\frac{H_{\rm ref}}{r_{\rm
ref}}\right)\left(\frac{r}{r_{\rm ref}}\right)^{-\frac{1}{10}},$ (66)
$\displaystyle B_{\phi}(r,z)$
$\displaystyle=\left(\frac{2^{27}~{}3^{17}\pi}{5~{}19^{14}}\right)^{\frac{1}{30}}\left(\frac{z}{H}\right)\dot{M}^{7/15}\rho_{c}^{-1/5}\Omega_{0}^{7/5}(C\gamma)^{-7/15}$
$\displaystyle c_{\rm s,ref}^{-2/5}\left(\frac{H_{\rm ref}}{r_{\rm
ref}}\right)\left(\frac{r}{r_{\rm ref}}\right)^{-\frac{8}{5}},$ (67)

where $H_{\rm ref}\equiv c_{\rm s,ref}/\Omega(r_{\rm ref})$. The following
estimates are obtained by substituting numerical values for the parameters.

$\displaystyle\rho_{\rm mid}(r)$ $\displaystyle=1.3\times 10^{-12}$
$\displaystyle\left(\frac{\dot{M}}{10^{-5}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{2/5}\left(\frac{\rho_{c}}{10^{-13}~{}{\rm
g~{}cm}^{-3}}\right)^{2/5}$ $\displaystyle\left(\frac{M_{\rm
star}}{0.2\thinspace
M_{\odot}}\right)^{3/5}\left(\frac{C\gamma}{10^{-2}}\right)^{-2/5}$
$\displaystyle\left(\frac{c_{\rm s,ref}}{180~{}{\rm m}~{}{\rm
s}^{-1}}\right)^{6/5}\left(\frac{r}{10{\rm AU}}\right)^{-\frac{9}{5}}~{}{\rm
g~{}cm}^{-3},$ (68) $\displaystyle B_{z}(r)$ $\displaystyle=4.6\times
10^{-1}\left(\frac{\dot{M}}{10^{-5}\thinspace M_{\odot}{\rm
yr}^{-1}}\right)^{8/15}$ $\displaystyle\left(\frac{\rho_{c}}{10^{-13}~{}{\rm
g~{}cm}^{-3}}\right)^{1/5}\left(\frac{M_{\rm star}}{0.2\thinspace
M_{\odot}}\right)^{3/10}\left(\frac{C\gamma}{10^{-2}{\rm
cm^{\frac{3}{2}}g^{-\frac{1}{2}}s^{-1}}}\right)^{7/15}$
$\displaystyle\left(\frac{c_{\rm s,ref}}{180~{}{\rm m}~{}{\rm
s}^{-1}}\right)^{-3/5}\left(\frac{r}{10{\rm AU}}\right)^{-\frac{9}{10}}{\rm
G},$ (69) $\displaystyle v_{\phi}(r,z)$
$\displaystyle=v_{K}-63\left(1+\frac{1}{2}\left(\frac{z}{H}\right)^{2}\right)\left(\frac{\dot{M}}{10^{-5}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{4/15}$
$\displaystyle\left(\frac{\rho_{c}}{10^{-13}~{}{\rm
g~{}cm}^{-3}}\right)^{-2/5}\left(\frac{M_{\rm star}}{0.2\thinspace
M_{\odot}}\right)^{2/5}\left(\frac{C\gamma}{10^{-2}{\rm
cm^{\frac{3}{2}}g^{-\frac{1}{2}}s^{-1}}}\right)^{-4/15}$
$\displaystyle\left(\frac{c_{\rm s,ref}}{180~{}{\rm m}~{}{\rm
s}^{-1}}\right)^{1/5}\left(\frac{r}{10{\rm AU}}\right)^{-\frac{7}{10}}~{}{\rm
m}~{}{\rm s}^{-1},$ (70) $\displaystyle v_{r}(r,z)$
$\displaystyle=-186\left(1+\frac{1}{2}\left(\frac{z}{H}\right)^{2}\right)\left(\frac{\dot{M}}{10^{-5}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{7/15}$
$\displaystyle\left(\frac{\rho_{c}}{10^{-13}~{}{\rm
g~{}cm}^{-3}}\right)^{-1/5}\left(\frac{M_{\rm star}}{0.2\thinspace
M_{\odot}}\right)^{1/5}\left(\frac{C\gamma}{10^{-2}{\rm
cm^{\frac{3}{2}}g^{-\frac{1}{2}}s^{-1}}}\right)^{8/15}$
$\displaystyle\left(\frac{c_{\rm s,ref}}{180~{}{\rm m}~{}{\rm
s}^{-1}}\right)^{-2/5}\left(\frac{r}{10{\rm AU}}\right)^{-\frac{1}{10}}~{}{\rm
m}~{}{\rm s}^{-1},$ (71) $\displaystyle B_{\phi}(r,z)$
$\displaystyle=-1.3\times
10^{-2}\left(\frac{z}{H}\right)\left(\frac{\dot{M}}{10^{-5}\thinspace
M_{\odot}{\rm yr}^{-1}}\right)^{7/15}$
$\displaystyle\left(\frac{\rho_{c}}{10^{-13}~{}{\rm
g~{}cm}^{-3}}\right)^{-1/5}\left(\frac{M_{\rm star}}{0.2\thinspace
M_{\odot}}\right)^{7/10}\left(\frac{C\gamma}{10^{-2}{\rm
cm^{\frac{3}{2}}g^{-\frac{1}{2}}s^{-1}}}\right)^{-7/15}$
$\displaystyle\left(\frac{c_{\rm s,ref}}{180~{}{\rm m}~{}{\rm
s}^{-1}}\right)^{-2/5}\left(\frac{r}{10{\rm AU}}\right)^{-\frac{8}{5}}{\rm
G}.$ (72)

These equations effectively reproduce our simulation results.

So far, we have derived the solution from the basic equations with the
assumptions explicitly stated, which, however, may be intuitively difficult to
understand. The aforemetioned power laws can be derived from a simple and
intuitive extension of the viscous accretion disk model. Besides the numerical
factors, the analytc solutions can also be derived by the following equations:

$\displaystyle-2\pi rv_{r}\Sigma=\dot{M},$ (73) $\displaystyle
v_{\phi}=\sqrt{\frac{GM_{c}}{r}}\equiv r\Omega,$ (74) $\displaystyle
H=\frac{c_{s}}{\Omega},$ (75) $\displaystyle
v_{r}=-\frac{B_{z}B_{\phi,s}}{\pi\Sigma\Omega},$ (76) $\displaystyle
B_{z}v_{r}=-\frac{\eta_{A}}{r}B_{z},$ (77) $\displaystyle
B_{\phi,s}=\left(\frac{H}{r}\right)^{2}\frac{B_{z}H}{\eta_{A}}v_{\phi}.$ (78)

The equations (73) to (76) have a similar form of the standard viscous
accretion disk model (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974),
in which equation (76) is

$\displaystyle v_{r}=-\frac{3}{2}\frac{\alpha c_{s}H}{r},$ (79)

where $\alpha$ is the Shakura-Sunyaev viscous parameter.

In our model, equations (77) and (78) describe the radial and azimuthal
balance between the magnetic field advection by the gas motion and magnetic
field drift by ambipolar diffusion. Note that the right hand side of equation
(78) has the factor of $(H/r)^{2}$, reflecting that the balance between the
vertical shear motion (not Keplerian rotation itself) and magnetic field drift
due to ambipolar diffusion determines the toroidal magnetic field $B_{\phi}$.
This is the key to deriving the disk structure we have identified.

Instead of equation (78), it is often assumed that

$\displaystyle B_{\phi,s}=\frac{B_{z}H}{\eta_{A}}v_{\phi},$ (80)

which corresponds to the assumption that the Keplerian rotation balances the
magnetic field drift by ambipolar diffusion (Krasnopolsky & Königl, 2002;
Braiding & Wardle, 2012; Hennebelle et al., 2016). We found that this relation
leads to considerably large $v_{r}$ (of the order of the Keplerian velocity),
which is inconsistent with the simulation results. The reason why this
relation is inappropriate for the circumstellar disk is that at the leading
order of $\epsilon$, the gravity does not depend on $z$ and is canceled out by
the centrifugal force in the circumstellar disk. Thus, we should consider the
balance between rotation and field drift at the order of $\epsilon^{2}$ to
estimate the toroidal magnetic field in the circumstellar disk (see the
derivations above for details).

The solutions given by equations (A) to (A) well reproduce our three
dimensional simulations. Furthermore, they are specified only by the central
star mass, mass accretion rate, equation of state (or gas temperature), and
ionization and recombination rate, and do not contain the viscous $\alpha$
parameter, which is usually extremely difficult to determine. Therefore, we
believe that the analytical solutions are useful for investigating the long-
term evolution of circumstellar disks.

## References

  * Braiding & Wardle (2012) Braiding, C. R., & Wardle, M. 2012, MNRAS, 422, 261
  * Galametz et al. (2019) Galametz, M., Maury, A. J., Valdivia, V., et al. 2019, A&A, 632, A5
  * Guilet & Ogilvie (2012) Guilet, J., & Ogilvie, G. I. 2012, MNRAS, 424, 2097
  * Guilet & Ogilvie (2013) —. 2013, MNRAS, 430, 822
  * Guillet et al. (2020) Guillet, V., Hennebelle, P., Pineau des Forêts, G., et al. 2020, A&A, 643, A17
  * Hennebelle et al. (2016) Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830, L8
  * Kawasaki et al. (2022) Kawasaki, Y., Koga, S., & Machida, M. N. 2022, MNRAS, 515, 2072
  * Kobayashi et al. (2023) Kobayashi, Y., Takaishi, D., & Tsukamoto, Y. 2023, MNRAS, arXiv:2303.04424
  * Krasnopolsky & Königl (2002) Krasnopolsky, R., & Königl, A. 2002, ApJ, 580, 987
  * Krasnopolsky et al. (2012) Krasnopolsky, R., Li, Z.-Y., Shang, H., & Zhao, B. 2012, ApJ, 757, 77
  * Kudoh & Shibata (1997) Kudoh, T., & Shibata, K. 1997, ApJ, 474, 362
  * Kuffmeier et al. (2020) Kuffmeier, M., Zhao, B., & Caselli, P. 2020, A&A, 639, A86
  * Kwon et al. (2009) Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.-F., & Kemball, A. J. 2009, ApJ, 696, 841
  * Lebreuilly et al. (2023) Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326
  * Li (1998) Li, Z.-Y. 1998, ApJ, 497, 850
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
  * Marchand et al. (2023) Marchand, P., Lebreuilly, U., Mac Low, M.-M., & Guillet, V. 2023, arXiv e-prints, arXiv:2301.01510
  * 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
  * Okuzumi et al. (2016) Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82
  * Ormel & Cuzzi (2007) Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413
  * Sato et al. (2016) Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shu (1983) Shu, F. H. 1983, ApJ, 273, 202
  * Tomida et al. (2017) Tomida, K., Machida, M. N., Hosokawa, T., Sakurai, Y., & Lin, C. H. 2017, ApJ, 835, L11
  * Tomida et al. (2015) Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117
  * Tsukamoto et al. (2015) Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, MNRAS, 452, 278
  * Tsukamoto et al. (2021a) Tsukamoto, Y., Machida, M. N., & Inutsuka, S. 2021a, ApJ, 913, 148
  * Tsukamoto et al. (2021b) Tsukamoto, Y., Machida, M. N., & Inutsuka, S.-i. 2021b, ApJ, 920, L35
  * 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. (2017) Tsukamoto, Y., Okuzumi, S., & Kataoka, A. 2017, ApJ, 838, 151
  * Tsukamoto et al. (2022) Tsukamoto, Y., Maury, A., Commerçon, B., et al. 2022, arXiv e-prints, arXiv:2209.13765
  * Valdivia et al. (2019) Valdivia, V., Maury, A., Brauer, R., et al. 2019, MNRAS, 488, 4897
  * Williams & Cieza (2011) Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67
  * Wu et al. (2004) Wu, Y., Wei, Y., Zhao, M., et al. 2004, A&A, 426, 503
  * Wurster et al. (2018) Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 476, 2063
  * Wurster et al. (2016) Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037
  * Zhao et al. (2016) Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050