#  Evolution of the Water Snow Line in Magnetically Accreting Protoplanetary
Disks

Shoji Mori Department of Astronomy, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo 113-0033, Japan  Astronomical Institute, Tohoku University,
6-3 Aramaki, Aoba-ku, Sendai 980-8578, Japan<EMAIL_ADDRESS>Satoshi
Okuzumi Department of Earth and Planetary Sciences, Tokyo Institute of
Technology, Meguro-ku, Tokyo 152-8551, Japan  Masanobu Kunitomo Department of
Physics, Kurume University, 67 Asahi-machi, Kurume, Fukuoka 830-0011, Japan
Xue-Ning Bai Institute for Advanced Study and Department of Astronomy,
Tsinghua University, Beijing 100084, China

###### Abstract

The low water content of the terrestrial planets in the solar system suggests
that the protoplanets formed within the water snow line. Accurate prediction
of the snow line location moving with time provides a clue to constrain the
formation process of the planets. In this paper, we investigate the migration
of the snow line in protoplanetary disks whose accretion is controlled by
laminar magnetic fields, which have been proposed by various nonideal
magnetohydrodynamic (MHD) simulations. We propose an empirical model of the
disk temperature based on our nonideal MHD simulations, which show that the
accretion heating is significantly less efficient than in turbulent disks, and
calculate the snow line location over time. We find that the snow line in the
magnetically accreting laminar disks moves inside the current Earth’s orbit
within 1 Myr after star formation, whereas the time for the conventional
turbulent disk is much longer than 1 Myr. This result suggests that either the
rocky protoplanets formed in such an early phase of the disk evolution, or the
protoplanets moved outward to the current orbits after they formed close to
the protosun.

##  1 Introduction

The terrestrial planets in the solar system are significantly drier than
solids in the outer solar system. The mass of the present Earth’s ocean is
only 0.02 % of the Earth’s mass (e.g., Charette & Smith, 2010). Even
considering the amount of water that may have been taken into the Earth’s
interior, the Earth’s water content is at most 2 wt% (Nomura et al., 2014; Fei
et al., 2017). In addition, current Venus has a dry atmosphere with a water
content of $\sim$30 ppm (e.g., Basilevsky & Head, 2003), and its interior
would only contain at most 10 % of the Earth’s water content (Elkins-Tanton et
al., 2007). Mars is likely to contain only up to 10-4 Mars’ mass of water
(Kurokawa et al., 2014). For Mercury, the only evidence of water is ice at the
north pole, which is only 10-8 Mercury’s mass (Lawrence et al., 2013).
Therefore, their present water contents would be $\lesssim 1$ wt%. In
contrast, bodies originating in the outer solar system, such as comets and
Neptune, have a higher water content over 10 wt% (e.g., Guillot, 2005; A’Hearn
et al., 2011; Rotundi et al., 2015).

The water content at the time of the formation of the terrestrial planets may
have been higher than the present value, but would still have been $\lesssim
1$ wt%. For the Earth, $\sim 0.1M_{\oplus}$ of water is hardly removed by
stellar irradiation (Machida & Abe, 2010; Hamano et al., 2013) or giant impact
(Genda & Abe, 2005; Schlichting et al., 2015; Biersteker & Schlichting, 2021).
Therefore, the initial Earth would not have contained much more water than it
does today. The high hydrogen isotope ratio of present Venus suggests that
massive dehydration has occurred in the past, but even so, the water content
of early Venus is estimated to be about the same as the Earth (Donahue et al.,
1982). Early Mars may have retained up to 0.05 % of the water in the Mars’
mass (Kurokawa et al., 2014).

The water content of a planet is largely determined by the location of the
water snow line during planet formation. The water snow line is defined as a
radius in protoplanetary disks (PPDs) where water sublimates and condenses. In
typical PPDs, the water snow line lies where the gas temperature is $\sim$
160–170 K (e.g., Hayashi, 1981). Dust inside the snow line is essentially dry,
whereas dust outside the snow line can contain up to $\sim$ 50 wt% of water
ice (Lodders, 2003). Planetesimals forming from icy dust can lose water to
some extent through radiogenic heating by 26Al, but the probability that the
dehydrating planetesimals form rocky planets with a significantly low water
content ($<$1 wt%) is still low with the initial solar abundance of 26Al
(Lichtenberg et al., 2019). Even if rocky planets form outside the snow line
from dehydrated planetesimals, they may acquire a significant amount of water
afterward by capturing icy particles in the gas disk (Sato et al., 2016).
Therefore, it is more natural to consider that the terrestrial planets in the
solar system formed inside the snow line and stayed there until the icy dust
outside the snow line was cleared out.

The position of the snow line is determined by the disk’s thermal structure.
In the simplest case where the disk is optically thin and its internal
temperature is determined by direct stellar irradiation, the snow line lies at
$\approx 3$ au for the present-day solar luminosity (Hayashi, 1981). However,
such a simple model does not apply to young, optically thick PPDs that contain
abundant dust grains. Optically thick disks can receive stellar radiation only
on their surface, and for this reason their interior temperature tends to be
lower than that of optically thin disks (Kusaka et al., 1970; Chiang &
Goldreich, 1997). In an optically thick disk that is passively irradiated by a
Sun-like star, the snow line can lie within 1 au (Sasselov & Lecar, 2000).

However, stellar irradiation is not the only heating mechanism for optically
thick PPDs. They can also be heated internally when the gas accretes toward
the central star and liberates its gravitational energy (Shakura & Sunyaev,
1973; Lynden-Bell & Pringle, 1974). Accretion heating is most efficient when
the energy is liberated deep inside the disk from which it takes a long time
for the heat to escape. This effect is sometimes called the blanketing effect
because an optically thick material acts as a thermal blanket (e.g., Milne,
1921; Chandrasekhar, 1935).

Previous studies for the snow line in optically thick PPDs (e.g., Sasselov &
Lecar, 2000; Garaud & Lin, 2007; Oka et al., 2011; Zhang & Jin, 2015; Xiao et
al., 2017) commonly adopt the classical viscous accretion disk model that
assumes vertically uniform viscosity. In this model, the blanketing effect is
particularly effective, pushing the snow line out to a few au from the central
star for typical values of the disk accretion rate (Sasselov & Lecar, 2000;
Garaud & Lin, 2007; Oka et al., 2011). This indicates that accounting for
accretion heating is essential to infer where the snow line was when the
embryos of the terrestrial planets formed in the solar nebula.

An important consequence of the vertically uniform kinematic viscosity assumed
in the classical viscous model is that the accretion energy is locally
dissipated, and hence is mainly deposited near the midplane, thus maximizing
the thermal blanket effect. The question is then whether such a vertical
heating profile is expected in a realistic model of disk accretion. In fact,
magnetohydrodynamic (MHD) accretion models for PPDs predict that energy
dissipation mainly takes place near the disk surface rather than near the
midplane. Turbulence generated by magnetorotational instability (MRI; Balbus &
Hawley, 1991) has been thought to be the dominant energy dissipation
mechanism. The dissipation profile is affected by nonideal MHD effects (Ohmic
diffusion, Hall effect, and ambipolar diffusion) brought about by the weakly
ionized nature of PPDs. Strong Ohmic diffusion in the inner disk region
suppresses the MRI turbulence around the midplane (Gammie, 1996; Fleming &
Stone, 2003), leading to a dead zone. In that case, the energy dissipates in
the surface MRI-active layers above the dead zone (Hirose & Turner, 2011),
rather than around the midplane. Furthermore, ambipolar diffusion operates in
the lower-density region and hence quenches even the turbulence in the upper
layer (Bai & Stone, 2013a; Bai, 2013; Gressel et al., 2015; Simon et al.,
2015), where the disk accretion is driven by the angular momentum removal by
magnetic disk winds (i.e., wind-driven accretion; Bai & Stone, 2013a).

With the MRI being suppressed, heating is due to Joule dissipation rather than
turbulent dissipation, which is not co-spatial with the location where the
accretion energy is liberated. Moreover, a substantial fraction of the
liberated energy is consumed for driving disk winds rather than leading to
dissipation. Mori et al. (2019) (hereafter MBO19, ) performed MHD simulations
with all the three nonideal MHD effects producing laminar magnetized accreting
disks, and showed that the Joule heating occurs at high altitude and hence the
accretion heating is much less efficient than in the viscous disk
model111Béthune & Latter (2020) also investigated the vertical distribution of
Joule heating with the nonideal MHD effects, but led to a different conclusion
from MBO19. In Section 4.1, based on the results of this paper, we discuss the
differences between the results of MBO19 and Béthune & Latter (2020).. It was
also shown that the Hall effect amplifies and suppresses Joule heating when
the magnetic field threading the disk is aligned and anti-aligned with the
disk rotation axis, respectively (MBO19).

In this paper, we show the temporal evolution of the location of the water
snow line in magnetically accreting PPDs. To do so, we construct a global
model of disk temperature evolution in the terrestrial-planet-forming region
in the disks, based on the results of MBO19. MBO19 suggested that even in the
early phase of PPD evolution, the disk temperature of the inner regions will
be determined by the stellar irradiation. We explore this scenario in much
greater detail. In addition, we also discuss the impact of the migration of
the snow line on the terrestrial planet formation. The water snow line
migrates with the temporal evolution of the disk’s thermal structure (e.g.,
Oka et al., 2011). Thus, considering that protoplanets of the terrestrial
planets should have formed within the snow line, the track of the snow line
constrains their formation time and location. We suggest that the formation
time of rocky protoplanets is strongly constrained, or that the protoplanets
have moved outward from the vicinity of the protosun to the present orbits
after the formation.

This paper is organized as follows. In Section 2, we propose a new disk
temperature model based on MHD simulation results. In Section 3, we show the
time evolution of the water snow line around a 1 $M_{\odot}$ star. In Section
4, we discuss the implications of our results for rock planet formation and
other possible heating mechanisms neglected in this paper.

##  2 MHD Wind-Driven Accretion Disk Model

Figure 1:  Schematic illustration of the expected thermal structure of PPDs.

We divide the inner region of PPDs into two zones in terms of the thermal and
turbulent states as summarized in Figure 1: MRI zones and dead zones. The MRI
zone is the very inner hot region ($T\gtrsim 1000$ K; Desch & Turner, 2015)
where MRI-turbulence is induced by the thermal ionization and heats the gas
with the turbulent viscosity. Outside of the MRI zone is the dead zone, where
the MRI is suppressed without thermal ionization, and cooler because of the
inefficient accretion heating. In the following, we focus on the dead zone
because we are interested in the evolution of the snow line where $T\sim
160~{}\rm K$.

###  2.1 Density Structure and Accretion Rate

We assume that the disk is vertically nearly isothermal and give the vertical
gas density profile as

$\rho(z)=\frac{\Sigma}{\sqrt{2\pi}H}\exp\left(-\frac{z^{2}}{2H^{2}}\right),$
(1)

where $z$ is the distance from the midplane, $H$ is the scale height, and
$\Sigma$ is the gas surface density. The scale height is related to the
isothermal sound speed $c_{\rm s}$ and local Keplerian frequency $\Omega$ as
$H=c_{\rm s}/\Omega$. The gas can indeed be regarded as vertically isothermal
at low altitude below heat sources (see MBO19, ).

Accretion in the dead zone is assumed to be entirely driven by magnetic winds.
In this limit, the mass accretion rate $\dot{M}$ can be written as (Wardle,
2007; Fromang et al., 2013)

$\dot{M}=\frac{4\pi r}{\Omega}w_{\rm wind},$ (2)

where $r$ is the radial distance, and $w_{\rm wind}$ is the $z\phi$ component
of the Maxwell stress exerted on both sides of the disk, which we call the
wind stress.

In principle, $w_{\rm wind}$ depends on the net flux of large-scale magnetic
fields threading the disk (Bai & Stone, 2013a, b; Bai, 2014), which is,
however, highly uncertain. For this reason, we opt for parameterizing $w_{\rm
wind}$ as (Suzuki et al., 2016)

$w_{\rm wind}=P_{\rm mid}\overline{\alpha_{z\phi}},$ (3)

where $P_{\rm mid}$ is the gas pressure at the midplane and
$\overline{\alpha_{z\phi}}$ is a dimensionless parameter that characterizes
the level of the wind stress. According to MHD simulations (Bai, 2017; Mori et
al., 2019), typical values of $\overline{\alpha_{z\phi}}$ range between
$10^{-4}$–$10^{-2}$. We take $\overline{\alpha_{z\phi}}=10^{-3}$ as a default
value and vary $\overline{\alpha_{z\phi}}$ in Section 3.2.3.

The parameter $\overline{\alpha_{z\phi}}$ should not be confused with the
$\alpha$ parameter in the standard viscous disk model (Shakura & Sunyaev,
1973). The former measures the stress vertically transporting angular momentum
to the magnetic wind, whereas the latter measures the radial angular momentum
transport due to the magnetic fields and turbulence within the disk. Moreover,
in general, magnetized wind is much more efficient in driving disk accretion.
This can be found by noting that $\Sigma$ in steady state is related to
$\dot{M}$ as $\Sigma\sim\dot{M}/(\overline{\alpha_{z\phi}}c_{\rm s}r)$ for
wind-driven accretion (see Equation (4) below) and $\Sigma\sim\dot{M}/(\alpha
c_{s}H)$ for viscous accretion. Comparing the two expressions, the two disk
models give the same value of $\Sigma$ for given $\dot{M}$ when
$\alpha\sim(r/H)\overline{\alpha_{z\phi}}$ (Wardle, 2007; Bai & Goodman, 2009;
Fromang, 2013; Bai, 2017).

Assuming wind-driven accretion with $\overline{\alpha_{z\phi}}$ being
constant, one can relate $\Sigma$ to $\dot{M}$. When the disk is isothermal
except at high altitude, the midplane pressure can be written as $P_{\rm
mid}=\Sigma c_{\rm s}^{2}/(\sqrt{2\pi}H)=\Sigma\Omega c_{\rm s}/\sqrt{2\pi}$.
Using this and Equations (2) and (3), we have

$\Sigma=\frac{\dot{M}}{2\sqrt{2\pi}\overline{\alpha_{z\phi}}c_{\rm s}r}.$ (4)

We assume that the disk is in a quasi-steady state, which follows from the
assumption that $\dot{M}$ varies on timescales much longer than the local gas
accretion timescale. We also assume that mass accretion is dominant over mass
loss. Thereby, the mass accretion rate $\dot{M}$ is radially constant (see
Suzuki et al., 2016; Bai, 2016) and is equal to the gas accretion rate onto
the central star. We determine $\dot{M}$ as a function of stellar age $t$
using the empirical relation for the stellar mass $M=1\,\rm M_{\sun}$ by
Hartmann et al. (2016) 222We have rewritten the empirical relation of
$\dot{M}$–$t$ for 0.7 $M_{\sun}$ stars (Equation (12) in Hartmann et al.
(2016)) into that for 1.0 $M_{\sun}$ stars by using the empirical relation of
$\dot{M}$–$M$ (Equation (11) in their paper).,

$\dot{M}=4\times 10^{-8\pm 0.5}\left({\displaystyle\frac{t}{1~{}{\rm
Myr}}}\right)^{-1.07}~{}{\rm M_{\sun}~{}{\rm yr}^{-1}},$ (5)

which is shown in Figure 2. Here, the stellar age $t$ is defined as the time
after star formation is completed, i.e., after the end of the protostellar
accretion phase and the arrival at the stellar birthline.

Figure 2:  Empirical relation between the mass accretion rate $\dot{M}$ and
stellar age $t$ (Equation (5); Hartmann et al. 2016). The shade shows the
observational scatter of 0.5 dex.

###  2.2 Disk Heating

We consider both stellar irradiation and accretion heating and give the
temperature $T$ in the disk interior as

$T=\left(T_{\rm irr}^{4}+T_{\rm acc,\,MHD}^{4}\right)^{1/4},$ (6)

where $T_{\rm irr}$ and $T_{\rm acc,\,MHD}$ represent the contributions from
irradiation and accretion, respectively.

The temperature in the irradiation-dominated disk is approximately given by
(Kusaka et al., 1970; Chiang & Goldreich, 1997)

$\displaystyle T_{\rm irr}$ $\displaystyle=$ $\displaystyle
110\left({\displaystyle\frac{r}{1{\rm\,au}}}\right)^{-3/7}\left({\displaystyle\frac{L}{\rm\,L_{\sun}}}\right)^{2/7}\left({\displaystyle\frac{M}{\rm\,M_{\sun}}}\right)^{-1/7}~{}\rm
K,$ (7)

where $L$ is the stellar luminosity. Equation (7) assumes that stellar light
is absorbed at about 4 $H$ above/below the midplane and half of the star is
hidden by the innermost region.

We take $L$ to be dependent on $t$ using the evolutionary model for a
$1\rm\,M_{\sun}$ star by Feiden (2016), as shown in Figure 3. We note that, in
their model, a $1\rm\,M_{\sun}$ star has the deuterium burning from 0.05 to
0.2 Myr and therefore $L$ is kept high. Here, we neglect the contribution of
accretion luminosity to $L$ and assume that the growth of stellar mass can be
neglected for simplicity (see Kimura et al. (2016) and Section 4.4 for the
uncertainties in the pre-main-sequence evolution).

Figure 3:  Temporal evolutions of stellar luminosity $L$ used in our
calculation, which are based on Feiden (2016) (see text for details).

For $T_{\rm acc,\,MHD}$ at the midplane, we use (see Appendixes A and B for
the derivation)

$T_{\rm acc,\,MHD}=\left[\left(\frac{3\dot{M}\Omega^{2}f_{\rm
heat}}{32\pi\sigma}\right)\left(\tau_{\rm
heat}+\frac{1}{\sqrt{3}}\right)\right]^{1/4},$ (8)

where $\sigma$ is the Stefan-Boltzmann constant, and $\tau_{\rm heat}$ is an
effective optical depth at the midplane and approximately corresponds to the
optical depth from infinity to the bottom of the heating layer (see Equation
(50) in MBO19 and Equation (A3) for the exact expression). The dimensionless
parameter $f_{\rm heat}$ ($\leq 1$) is the fraction of the accretion heat
deposited inside the disk, with $f_{\rm heat}<1$ when the wind material
carries away part of the energy generated by accretion. We set the default
value of $f_{\rm heat}$ to 1 to give a maximum estimate of $T_{\rm
acc,\,MHD}$, where $f_{\rm heat}=1$ is the limiting case that all the
liberated energy is consumed for heating. We discuss the effects of varying
$f_{\rm heat}$ in Section 3.2.3.

We assume that the Rosseland mean opacity $\kappa_{\rm R}$ is constant
throughout the disk and write

$\tau_{\rm heat}=\kappa_{\rm R}\Sigma_{\rm heat},$ (9)

where $\Sigma_{\rm heat}$ represents the mass column depth from infinity to
the bottom of the heating region (see Appendix A and Equation (A7) for the
exact expression of $\Sigma_{\rm heat}$). The remaining task is to evaluate
$\Sigma_{\rm heat}$.

###  2.3 Estimating the Heating Layer Depth from Disk Ionization Structure

Figure 4:  Vertical profile of Am at $r=1$ au of the MHD accretion disk model
at $t=0.6$ Myr. The dashed horizontal line indicates the critical ambipolar
Elsasser number (Am = 0.3) of the current layer.

We estimate $\Sigma_{\rm heat}$ from the vertical resistivity structure of the
disk. We particularly focus on ambipolar diffusion because it is the dominant
non-ideal MHD effect near the disk surface. The strength of ambipolar
diffusion can be measured by the the ambipolar Elsasser number

${\rm Am}=\frac{v_{\rm A}^{2}}{\eta_{\rm A}\Omega},$ (10)

where $v_{\rm A}$ is the Alfvén speed and $\eta_{\rm A}$ is the ambipolar
diffusivity. The ambipolar Elsasser number is in general independent of the
magnetic field strength because both $v_{\rm A}^{2}$ and $\eta_{\rm A}$ scale
quadratically with it. In inner regions of protoplanetary disks, ${\rm Am}$
usually falls below unity at the midplane and increases toward the disk
surface (Wardle, 2007; Salmeron & Wardle, 2008; Bai & Stone, 2013a). The
vertical distribution of Am is shown in Figure 4. MHD simulations show that
electric currents decay where ${\rm Am}\ll 1$ (Bai & Stone, 2013a; Gressel et
al., 2015, MBO19). Therefore, one can expect that a current layer should
always lie above the critical layer where ${\rm Am}$ is around unity. In this
study, we conservatively take the critical Elsasser number to be $0.3$ and
write

$\Sigma_{\rm heat}=f_{\rm depth}\Sigma_{\rm Am=0.3},$ (11)

where the prefactor $f_{\rm depth}$ accounts for the fact that an actual
current sheet can occur above the ${\rm Am}=0.3$ layer (see Figure 4), and
$\Sigma_{\rm Am=0.3}$ is the gas surface density above the layer. In Appendix
C, we demonstrate that all MHD simulations presented in MBO19 satisfy $f_{\rm
depth}\lesssim 1$. The typical ranges of $f_{\rm depth}$, which is affected by
Hall effect, are $\sim$ 0.3–1 and $\sim$ 0.01–0.1 when the vertical magnetic
field is parallel and anti-parallel to the disk rotation axis, respectively.
As in the case of $f_{\rm heat}$, we take the default value of $f_{\rm depth}$
to unity to provide a maximum estimate for $T_{\rm acc,\,MHD}$, but also
quantify the effect of varying the value of $f_{\rm heat}$ in Section 3.2.3.

To determine the height of the ${\rm Am}=0.3$ layer, we compute the vertical
profile of $\eta_{\rm A}$ using the ionization model of MBO19 that includes
charged dust grains (see Section 2.1 of MBO19 for details). The grains
regulate the ionization fraction of the disk gas by capturing charged
particles at a rate roughly proportional to the total surface area of the
grains per unit volume (Bai & Goodman, 2009; Okuzumi, 2009). In this study, we
fix the grain size to $0.1~{}\micron$ and the internal density to
$1.4\rm\,g\,cm^{-2}$, and vary the total grain surface area by changing the
dust-to-gas mass ratio $f_{\rm dg}$. The ionizing sources include galactic
cosmic rays (Umebayashi & Nakano, 2009), stellar X-rays (Igea & Glassgold,
1999; Bai & Goodman, 2009), and radionuclides (Umebayashi & Nakano, 2009). We
represent all ion species with a single species by following Okuzumi (2009).
Using the isothermal density profile Equation (1), the height $z_{\rm Am=0.3}$
of the ${\rm Am=0.3}$ layer is related to the column depth of the layer as

$\Sigma_{\rm Am=0.3}=\frac{\Sigma}{2}{\rm erfc}\left(\frac{z_{\rm
Am=0.3}}{\sqrt{2}H}\right).$ (12)

###  2.4 Numerical Procedures and Parameter Choices

To summarize, Equations (4) and (6) determine the surface density $\Sigma(r)$
and midplane temperature $T(r)$ of our quasi-steady wind-driven accretion disk
model. Because the right-hand sides of these equations involve gas density and
temperature, one solves them iteratively to obtain a self-consistent solution.
We set the temperature tolerance to be $10^{-4}$. We have confirmed that the
converged temperature profile does not depend on the input profile of the
initial guess. At each iteration, we also compute the vertical distribution of
the ionization fraction to determine the column depth $\Sigma_{\rm heat}$ of
the heating layers given by Equation (11). When Am at the midplane is higher
than 0.3, we take $\Sigma_{\rm Am=0.3}$ to be $\Sigma/2$ by assuming that the
heating occurs at the midplane.

Our model involves five parameters $\\{\kappa_{\rm R},f_{\rm
dg},\overline{\alpha_{z\phi}},f_{\rm heat},f_{\rm depth}\\}$. The default
values are taken to be $\kappa_{\rm R}=5\rm\,cm^{2}\,g^{-1}$, $f_{\rm
dg}=0.01$, $\overline{\alpha_{z\phi}}=10^{-3}$, $f_{\rm heat}=1$, and $f_{\rm
depth}=1$. The default value of $\kappa_{\rm R}$ is consistent with the
opacity for 0.1 $\micron$ sized grains with the ISM abundance computed by
Pollack et al. (1985).

In this study, we treat $\kappa_{\rm R}$ and $f_{\rm dg}$ as independent
parameters, differently in MBO19, to investigate their effects on the disk
temperature. Effectively, we give the total extinction and geometric cross
sections of dust grains per unit gas mass independently. Because the ratio
between the two generally depends on the composition and size distribution of
grains, one can interpret our treatment as taking into account uncertainties
in grain composition and size distribution.

##  3 Migration of the Snow Line

In this section, we present the evolution of the water snow line location from
our MHD disk model and compare it with the prediction from the standard
viscous disk model. For each disk model, we calculate the radial temperature
profile from 0.08 au to 20 au and search for the snow line location assuming
the ice sublimation temperature of 170 K. For the viscous model, we express
the viscosity as $\nu=\alpha c_{\rm s}H$ and take the parameter $\alpha$ to be
constant throughout the disk. Assuming steady accretion, the gas surface
density of the viscous disk is $\Sigma=\dot{M}/(3\pi\alpha c_{\rm s}H)$
(Lynden-Bell & Pringle, 1974). The midplane temperature is given by

$T=\left(T_{\rm irr}^{4}+T_{\rm acc,\,visc}^{4}\right)^{1/4},$ (13)

where

$T_{\rm
acc,\,visc}=\left[\left(\frac{9\dot{M}\Omega^{2}}{32\pi\sigma}\right)\left(\frac{\tau_{\rm
mid}}{2}+\frac{1}{\sqrt{3}}\right)\right]^{1/4}\approx\left(\frac{3\kappa_{\rm
R}\dot{M}^{2}\Omega^{3}}{128\pi^{2}\sigma\alpha c_{\rm s}^{2}}\right)^{1/4}$
(14)

represents the contribution from viscous accretion heating with $\tau_{\rm
mid}=\kappa_{\rm R}\Sigma/2$ being the vertical optical depth to the midplane
(Hubeny, 1990; Nakamoto & Nakagawa, 1994; Oka et al., 2011). The effective
optical depth $\tau_{\rm heat}$ for the vertically uniform viscosity is
$\tau_{\rm heat}=\tau_{\rm mid}/2$ (Equation (A8)). The final expression of
Equation (14) assumes $\tau_{\rm heat}\gg 1$, which holds true in our
calculations (see Figure 5). Note that Equation (14) should be taken as an
implicit equation for $T_{\rm acc,\,visc}$ because the right-hand side depends
on $T_{\rm acc,\,visc}$ through $c_{\rm s}^{2}\propto T_{\rm acc,\,visc}$.

For a fair comparison, we take $\alpha$ so that the MHD and viscous disk
models have the similar gas surface density. Since
$\alpha\sim(r/H)\overline{\alpha_{z\phi}}$ (see the discussion in Section
2.1), for a typical turbulent protoplanetary disk with $H/r\sim 0.1$, one has
$\alpha\sim 10\overline{\alpha_{z\phi}}$. Therefore, we set $\alpha=10^{-2}$.
Although this is just an order-of-magnitude estimate, we have confirmed that
the ratio of $\Sigma$ in the MHD model to that in the viscous one is in the
range of $\approx$ 0.5–1.4 at 1 au for the fiducial parameter set.

###  3.1 Fiducial Case

Figure 5:  Upper panels: radial profiles of the midplane temperature from the
MHD accretion disk model (Equation (6); thick solid line) and from the viscous
accretion disk model (Equation (14); thin solid line) at $t$ = 0.6 Myr (left
panel) and $t$ = 10 Myr (right panel). For reference, the dotted, dot-dashed,
and dashed lines show $T_{\rm irr}$, $T_{\rm acc,\,MHD}$, and $T_{\rm
acc,\,visc}$ (Equations (7), (8), and (14)), respectively. The horizontal
solid line marks the sublimation temperature of water ice, $T=170\,{\rm\,K}$.
Lower panels: radial profiles of the effective optical depths $\tau_{\rm
heat}$ at the midplane for the viscous model ($\tau_{\rm heat}=\tau_{\rm
mid}/2$; Equation (A8); dashed line) and for the MHD model (Equation (11);
solid line).

We begin by the default case where $\kappa_{\rm R}=5\rm\,cm^{2}\,g^{-1}$,
$f_{\rm dg}=0.01$, $\overline{\alpha_{z\phi}}=10^{-3}$, $f_{\rm heat}=1$, and
$f_{\rm depth}=1$. The upper panels of Figure 5 show the radial temperature
profiles from the MHD and viscous disk models at two different stellar ages
(see Figure 4 in MBO19 for the typical vertical temperature structure). In the
viscous accretion model, accretion heating is the dominant source of disk
heating at $r\lesssim 10~{}\rm au$ (for $t=$ 0.6 Myr) and $r\lesssim 3~{}\rm
au$ (for $t=$ 10 Myr). In contrast, in the MHD model, accretion heating only
gives a minor contribution, and the midplane temperature is approximately
given by $T_{\rm irr}$ at all radii. For both models, the disk cools with
time, but for different reasons. In the viscous model, the temperature
decrease is mainly due to the decrease of $\dot{M}$ in $T_{\rm acc,\,visc}$
(Equation (14)). In the MHD model, the temperature evolution is rather driven
by the stellar evolution shown in Figure 3.

The inefficient accretion heating in the MHD model can be understood by
looking at the optical depth of the heating layer, $\tau_{\rm heat}$, shown in
the lower panels of Figure 5. In the viscous model, the effective optical
depth $\tau_{\rm heat}$ at the midplane is much higher than unity, with
$\tau_{\rm heat}\approx 500$ (for $t=$ 0.6 Myr) to $\approx 80$ (for $t=$ 10
Myr) at 1 au. Because of the strong blanketing effect, the accretion heating
is the dominant heating mechanism in the inner part of the viscous disk. In
contrast, in the MHD model, the optical depth $\tau_{\rm heat}$ of the heating
layer is only 1–3, and therefore the blanket effect is inefficient. This value
of $\tau_{\rm heat}$ comes from our ionization calculations showing that
$\Sigma_{\rm heat}\approx 0.2$–$0.6\rm\,g\,cm^{-2}$. Interestingly, the depth
of the heating layer is insensitive to the disk age and radial location. This
is because $\Sigma_{\rm heat}$ partly reflects the attenuation depth of
ionizing X-rays. It should be noted, however, that $\Sigma_{\rm heat}$ also
depends on the recombination rate of charged particles and hence on the
abundance of small grains in the disk. We show this in Section 3.2.2.

Figure 6:  Radial location of the snow line at the midplane as a function of
the stellar age in the MHD model (blue solid) and the standard viscous model
(red dotted). The shades show uncertainty coming from observational dispersion
of the mass accretion rate (see Hartmann et al., 2016). The arrows indicate
the stellar ages when the snow lines pass 1 au, which are 0.6 and 10 Myr for
the MHD and viscous models, respectively. The gray dashed line is for the
passive disk with $T=T_{\rm irr}$.

In both disk models, the snow line migrates inward because the disk cools with
time. Figure 6 shows the radial location of the snow line at the midplane in
the two models as a function of time. In the viscous disk model, the snow line
passes 1 au at $t\approx$ 10 Myr. In contrast, in the MHD model, the snow line
arrives 1 au at $t\approx$ 0.6 Myr because of inefficient accretion heating.
For comparison, we also show in Figure 6 the snow line position in a passively
irradiated disk with no internal heating, i.e., $T=T_{\rm irr}$. As expected,
the snow line location in the MHD disk almost agrees with that in the
passively irradiated disk at least until the snow line crosses 1 au.

###  3.2 Parameter Dependence

We here study the dependence of the snow line tracks in the MHD model on the
opacity, ionization fraction, accretion stress, $f_{\rm heat}$, and $f_{\rm
depth}$ to understand when the snow line track is largely affected by those
parameters. The used parameters and obtained results are summarized in Table
1.

Table 1: Summary of used parameters and the time $t_{\rm SL,1au}$ that the snow line reaches 1 au. Model | $\kappa_{\rm R}~{}{\rm[cm^{2}\,g^{-1}]}$ | $f_{\rm dg}$ | $\overline{\alpha_{z\phi}}$ | $f_{\rm heat}$ | $f_{\rm depth}$ | $t_{\rm SL,1au}$ [Myr] | Note
---|---|---|---|---|---|---|---
MHD | 5 | 0.01 | 10^-3 | 1 | 1 | 0.60 | Figure 6; fiducial case
viscous | 5 | 0.01 | 10^-3 | 1 | 1 | 12 | Figure 6
MHD | 50 | 0.01 | 10^-3 | 1 | 1 | 1.8 | Figure 7
MHD | 0.5 | 0.01 | 10^-3 | 1 | 1 | 0.47 | Figure 7
MHD | 5 | 0.1 | 10^-3 | 1 | 1 | 0.50 | Figure 8
MHD | 5 | 0.001 | 10^-3 | 1 | 1 | 2.2 | Figure 8
MHD | 5 | 0.001 | 10^-2 | 1 | 1 | 2.6 | Figure 10
MHD | 5 | 0.001 | 10^-4 | 1 | 1 | 2.1 | Figure 10
MHD | 5 | 0.001 | 10^-3 | 0.1 | 1 | 0.61 | Figure 11
MHD | 5 | 0.001 | 10^-3 | 1 | 0.1 | 0.64 |

####  3.2.1 Variation with the Opacity

Figure 7:  Evolution of the snow line location in the MHD model for different
values of $\kappa_{\rm R}$. The shades show uncertainty coming from
observational dispersion of the mass accretion rate. The arrows indicate the
stellar age when the snow lines pass 1 au. The gray dashed line is for the
passive disk with $T=T_{\rm irr}$.

Since $\tau_{\rm heat}=\kappa_{\rm R}\Sigma_{\rm heat}$, opacity is an
important factor in our model. Figure 7 shows the snow line tracks with
varying the opacity from the fiducial value of $\kappa_{\rm
R}=5\rm\,cm^{2}\,g^{-1}$. When we increase $\kappa_{\rm R}$ to
$50\rm\,cm^{2}\,g^{-1}$, $\tau_{\rm heat}$ is increased to $\approx 20$, and
consequently the arrival of the snow line at 1 au is delayed to $t\sim$ 2 Myr.
This demonstrates that accretion heating can dominate over irradiation heating
even in MHD accretion disks if the opacity is sufficiently high. Decreasing
the opacity from the fiducial value has little effect on the snow line
location because accretion heating is subdominant for such low opacities. The
depth $\Sigma_{\rm heat}$ of the heating layer is insensitive to $\tau_{\rm
heat}$ because non-thermal ionization and recombination reactions depend on
temperature only weakly. In short, we here find that the arrival time of the
snow line is $\gg 1$ Myr when $\tau_{\rm heat}\gg 10$.

####  3.2.2 Variation with the Ionization Fraction

Figure 8:  Same as Figure 7, but for different values of $f_{\rm dg}$. Note
that $f_{\rm dg}$ is the dust-to-gas mass ratio used in the ionization model
and does not affect the disk opacity, which is fixed to $5~{}\rm
cm^{2}~{}g^{-1}$ here.  Figure 9:  Optical depths $\tau_{\rm heat}$ at the
heating layer for different values of $f_{\rm dg}$ at $t=1$ Myr.

In our ionization model, the dust-to-gas mass ratio $f_{\rm dg}$ controls the
disk ionization fraction because grain charging significantly contributes to
the removal of ionized particles from the disk gas. A lower value of $f_{\rm
dg}$ gives a higher ionization fraction, a larger $\Sigma_{\rm heat}$, and
consequently a higher midplane temperature. To see this effect quantitatively,
we show in Figure 8 the snow line tracks with varying $f_{\rm dg}$ from the
fiducial value of $0.01$, with the opacity value fixed. We also show in Figure
9 the radial profile of $\tau_{\rm heat}$ for the cases shown in Figure 8. One
can see that decreasing $f_{\rm dg}$ by a factor of 10 from the fiducial value
leads to an increase of $\tau_{\rm heat}$ by approximately the same factor
(corresponding to $\Sigma_{\rm heat}\approx 6\rm\,g\,cm^{-2}$), resulting in
the arrival of the snow line at $1~{}\rm au$ delayed to 2 Myr. Thus, with a
sufficiently high ionization fraction, accretion heating can become a dominant
heating mechanism even in MHD accretion disks. For $f_{\rm dg}\gtrsim 0.01$,
accretion heating is subdominant, and hence variation with $f_{\rm dg}$ little
affects the snow line evolution.

The optical depth $\tau_{\rm heat}$ depends not only on the opacity but also
on the ionization fraction via $\Sigma_{\rm heat}$, and the both in turn
depend on the dust model. For instance, when the total dust abundance is
increased, the increase in opacity and the decrease in $\Sigma_{\rm heat}$ can
be comparable, and hence $\tau_{\rm heat}$ might not vary significantly. To
properly calculate the disk temperature, it is important to consistently give
the opacity and ionization fraction using an identical dust model.

####  3.2.3 Variation with the MHD Parameters

Figure 10:  Same as Figure 7, but for $f_{\rm dg}=10^{-3}$ with different
values of $\overline{\alpha_{z\phi}}$. The gray line is for the fiducial value
of $\overline{\alpha_{z\phi}}=10^{-3}$.

We have shown in Section 3.2.2 that a high $\kappa_{\rm R}$ and/or a low
$f_{\rm dg}$ makes accretion heating dominant even in the MHD model. The
question is now how strongly the level of the accretion heating depends on the
rather unconstrained parameters involved in the MHD disk model, namely
$\overline{\alpha_{z\phi}}$, $f_{\rm heat}$, and $f_{\rm depth}$.

The magnetic stress $\overline{\alpha_{z\phi}}$ mainly depends on the net flux
of the vertical magnetic field (e.g., Hawley et al., 1995; Bai, 2013) which
evolves with the disk evolution (e.g., Lubow et al., 1994; Guilet & Ogilvie,
2014; Okuzumi et al., 2014; Takeuchi & Okuzumi, 2014; Leung & Ogilvie, 2019).
Thus, $\overline{\alpha_{z\phi}}$ may vary widely.

It can be readily expected that the midplane temperature in the MHD model must
be insensitive to $\overline{\alpha_{z\phi}}$. This parameter determines the
magnitude of the wind accretion stress (see Equation (3)) and controls the
disk surface density $\Sigma$ (see Equation (4)). However, it is the column
depth $\Sigma_{\rm heat}$ of the heating layer, not the total column density
$\Sigma$, that determines $T_{\rm acc,\,MHD}$. The depth of the heating layer
is determined by the ionization structure well above the midplane, and is
therefore insensitive to $\Sigma$. We demonstrate this in Figure 10, where we
plot the snow line tracks for $f_{\rm dg}=10^{-3}$ (see also Figure 8) but
with different values of $\overline{\alpha_{z\phi}}$. It is clear that
changing $\overline{\alpha_{z\phi}}$ has essentially no effect on the snow
line location.

Figure 11:  Same as Figure 7, but for $f_{\rm dg}=10^{-3}$ with different
values of $f_{\rm heat}$. The gray line is for the fiducial value of $f_{\rm
heat}=1$.

The calculations shown above adopt $f_{\rm heat}=f_{\rm depth}=1$. Such
calculations give the maximum estimate of $T_{\rm acc,\,MHD}$. We also vary
the two parameters to illustrate how the choice affects the snow line location
from our model.

As explained in Sections 2.2 and Appendix B, $f_{\rm heat}$ is the fraction of
the energy dissipated inside the disk to the energy liberated by wind-driven
accretion. This number is much smaller than one when the wind carries away a
large fraction of the accretion energy. The value of $f_{\rm heat}$ depends on
the polarity of the net vertical magnetic field as well as $f_{\rm depth}$,
and the range is approximately 10-3 to 1. In Figure 11, we plot the snow line
track for $f_{\rm dg}=10^{-3}$, but now assuming $f_{\rm heat}=0.1$, i.e.,
only 10% of the accretion energy is consumed by the disk’s internal heating.
For such as low value of $f_{\rm heat}$, accretion heating is no longer
significant even with $f_{\rm dg}=10^{-3}$, and the snow line arrives at 1 au
within 1 Myr as in the fiducial model (Figure 6).

The factor $f_{\rm depth}$ is less than unity when Joule heating occurs at an
altitude higher than where ${\rm Am}=0.3$. For instance, as shown in Figure 12
of Appendix C, $f_{\rm depth}$ can be as low as 0.01–0.1 when the magnetic
field threading the disk is anti-aligned with the disk rotation axis. Since
$T_{\rm acc,\,MHD}\propto(f_{\rm heat}f_{\rm depth})^{1/4}$ as long as
$\tau_{\rm heat}\gtrsim 1$, varying $f_{\rm depth}$ from 1 to 0.1 gives the
same effect as varying $f_{\rm heat}$ from 1 to 0.1 (see Figure 11).

##  4 Discussion

###  4.1 Comparison with Béthune & Latter (2020)

Here we compare our results with Béthune & Latter (2020), who also
investigated accretion heating in magnetic laminar PPDs with non-ideal MHD
calculations. Béthune & Latter (2020) showed that Joule heating can
efficiently warm the disk interior, apparently in contradiction with the
results by MBO19. Béthune & Latter (2020) speculated that the discrepancy
arises from the different opacities adopted in the two studies: Béthune &
Latter (2020) used the opacity of the order of $5\rm\,cm^{2}\,g^{-1}$, while
MBO19 adopted a ten times smaller value.

The results presented in this study indicate that the combination of the disk
opacity and ionization state is a key to understanding the discrepancy.
According to Equation (8), the temperature of the disk interior is determined
by the optical depth $\tau_{\rm heat}$ of the heating layer, which is the
product of the opacity and the column depth $\Sigma_{\rm heat}$ of the layer.
The disk ionization state plays an important role here because it determines
$\Sigma_{\rm heat}$. As demonstrated in Section 3.1, the value of $\kappa_{\rm
R}=5~{}\rm cm^{2}~{}s^{-1}$ alone is not sufficient to warm the disk interior.
Efficient heating also requires a low dust abundance producing a high
ionization fraction and hence a high $\Sigma_{\rm heat}$ (Section 3.2.2).
Indeed, the fiducial model of Béthune & Latter (2020) adopted a relatively
high ionization fraction corresponding to the dust-free limit in our model. In
contrast, in the model of MBO19, the opacity was scaled with the dust-to-gas
ratio used in the ionization calculations and therefore decreased whenever
$\Sigma_{\rm heat}$ is increased. These results indicate that treating the
opacity and ionization fraction consistently using an identical dust model is
crucial for properly evaluating the efficiency of Joule heating.

###  4.2 Implications for the Formation of the Inner Solar System

####  4.2.1 Formation of the Dry Earth

The results presented in Section 3 show that the snow line in a magnetically
accreting disk around a Sun-like star passes the current Earth’s orbit as
early as 1 Myr after star formation except for particular parameter choices.
In such an early stage of disk evolution, the region outside the snow line is
always abundant in icy particles, and therefore protoplanets at 1 au would
acquire a significant amount of water ($\gtrsim 1~{}\%$ of the protoplanet
mass) by accreting the icy particles (Sato et al., 2016; Ida et al., 2019).
Therefore, our results strongly indicate that either our Earth did not form at
its current orbit or some mechanism prevented icy particles from migrating
into the 1 au region in the solar nebula.

There are some possible scenarios that may explain the low water content of
the Earth. One is that Jupiter’s core formed earlier than the Earth and carved
a gas gap that blocked inward migrating icy particles (Morbidelli et al.,
2016). If this is the case, our results presented in Section 3 suggest that
the proto-Jupiter formed within 1 Myr after the solar nebula formation. This
is consistent with what is inferred from a recent meteoritic analysis (Kruijer
et al., 2017).

Alternatively, Johansen et al. (2021) proposed a possibility that when the
temperature of planetary atmosphere is high enough, the water vapor sublimated
from icy pebbles accreting to the protoplanet is recycled back to the PPD.
This reduces the water content as the protoplanet grows with the pebble
accretion, and the final water content might match with the current water
content of the Earth.

Another scenario is that the Earth’s embryo formed at a close-in orbit ($\sim
0.1$ au) and then migrated to the current orbit after the icy particles in the
nebula had been depleted. For instance, there are models suggesting that rocky
planetesimals form at the inner boundary of the dead zone where the
temperature is $\approx 1000$ K (e.g., Kretke et al., 2009; Dzyurkevich et
al., 2010; Drążkowska et al., 2013). Ogihara et al. (2015) show that
planetesimals forming in the close-in region move outward when
magnetorotationally driven disk winds create a positive surface density slope.
However, whether the migration of the close-in planetesimals occurs after the
icy particles outside the snow line have been sufficiently depleted would
depend on the efficiency of the wind-driven accretion and wind mass loss in
the inner region.

####  4.2.2 Dichotomy of Planetesimals

The position of the snow line may be a key to understanding a dichotomy
between carbonaceous chondrites (CCs) and noncarbonaceous chondrites (NCs) in
the solar system. Meteoritic analyses of CCs and NCs based on the isotopic
compositions suggest that the parent bodies were born in the spatially
separated reservoirs of CC/NC-like materials and the reservoirs were formed at
different timings (see Kruijer et al., 2017; Kleine et al., 2020). Kruijer et
al. (2017) suggested that the two reservoirs were separated due to the early
forming Jupiter creating a gap in the solar nebula. Lichtenberg et al. (2021)
proposed another scenario that explains the dichotomy in a viscously evolving
disk model. In their scenario, the snow line moves outward as the viscously
heated region expands with disk formation during the Class I phase. Thus, two
types of planetesimals with different formation regions and timings are formed
around the snow line (Drążkowska & Dullemond, 2018) and are responsible for
the dichotomy.

The inefficient accretion heating in the MHD disk model may affect the
scenario in Lichtenberg et al. (2021). In the MHD model, at least in the Class
II phase, the snow line lies further closer to the central star than assumed
in Lichtenberg et al. (2021). For the Class I phase, it is uncertain if this
is the case, e.g., a self-gravitational instability might heat the disk as
much as viscous heating (see Section 4.5.2). If accretion heating is
inefficient only in the Class II phase, the formation region of planetesimals
in the Class II phase should overlap that in the Class I phase. In this case,
numerous planetesimals with compositions intermediate between CCs and NCs may
be formed, inconsistent with the meteorite analyses. If accretion heating is
inefficient in both the Class I and Class II phases, the two types of
planetesimals could be spatially well separated, although they would form in a
more inner region than expected in Lichtenberg et al. (2021).

###  4.3 Outward Migration of the Snow Line

In this paper, the direction of the snow line migration is only inward as in
Figure 6, whereas in Oka et al. (2011), the snow line moves outward after the
inward migration. The outward migration of the snow line occurs when the disk
becomes optically thin and thereby the stellar irradiation gets efficient.
According to Oka et al. (2011), the outward migration occurs when the
accretion rate decreases to less than $\dot{M}\lesssim
10^{-9}\rm\,\rm\,M_{\sun}\,\rm\,yr^{-1}$. At such low accretion rate, the disk
is near dispersal (in our model more than 10 Myrs). We are primarily concerned
with snow line evolution in the bulk disk lifetime, and so the outward
migration of the snow line is not focused on.

###  4.4 Uncertainty of Stellar Luminosity Evolution

We showed in Section 3.1 that the snow line evolution in magnetically
accreting PPDs is primarily determined by the evolution of the stellar
luminosity. In this study, we have adopted a conventional pre-main-sequence
evolutionary model.

We adopt the evolutionary model by Feiden (2016), which is also used by
Hartmann et al. (2016) to determine stellar mass and age. However, we note
that the initial condition of Feiden (2016) is much more luminous than the
birthline by Stahler & Palla (2005) and therefore our stellar luminosity and
thus $T_{\rm irr}$ are overestimated in particular in the early phase
($\lesssim 1\,$Myr; see also Section 2.2 of Hartmann et al., 1998). If time is
defined as the time after the luminosity reaches that of the birthline of
Stahler & Palla (2005), the time when the snow line passes 1 au becomes
earlier to $t\approx 0.2$ Myr for the fiducial case of the MHD model.

The initial luminosity of pre-main-sequence stars may be even lower than that
of Stahler & Palla (2005). They derived the birthline by assuming a large
fraction of accretion energy is injected into the protostar (called high-
entropy accretion). However, recent studies have shown that the evolution of
young stars depends on how much entropy of accreting materials is injected
into the star (e.g. Hartmann et al., 1997; Hosokawa et al., 2011; Baraffe et
al., 2012; Kunitomo et al., 2017). If the gas looses most of the entropy
before it reaches the star (called low-entropy accretion; i.e., most of
accretion energy is radiated away), the luminosity of pre-main-sequence stars
becomes much lower than fiducial models. In this case, the snow line passes 1
au much earlier.

We also note that in the present paper, we neglect the growth of stellar mass
and the accretion luminosity, $L_{\rm acc}$333The accretion luminosity is the
radiation from the accretion shock surface, while the intrinsic luminosity
comes from the stellar photosphere. (see Section 2.2). The accretion still
remains vigorous until $\approx 1\,$Myr (see Figure 3). However, in the case
of high-entropy accretion, $L_{\rm acc}\lesssim 1\,{\rm L_{\odot}}$ after the
star leaves the birthline and always satisfies $L_{\rm acc}\ll L_{\rm int}$,
where $L_{\rm int}$ is the stellar intrinsic luminosity. Hence we can safely
neglect $L_{\rm acc}$. In the case of low-entropy accretion, $L_{\rm acc}$ can
be higher than $L_{\rm int}$, but the total luminosity is comparable to or
less than that of the high-entropy accretion case. Therefore the stellar
luminosity in this paper is unlikely to be underestimated and the conclusion
of this paper (i.e., early arrival of the snow line to the present
terrestrial-planet orbits) should not be affected by the uncertainties of
stellar evolution models. Future work on this issue should include the
protostellar accretion to the stellar evolution model and investigate the
influence of the uncertainty in the stellar evolution model on the snow line
location.

###  4.5 Other Heating Mechanisms

We here discuss the possibility that other heating mechanisms influence the
location of the snow line.

####  4.5.1 Hydrodynamic Turbulence

Hydrodynamic turbulence can be a heat source other than Joule heating.
Hydrodynamic instabilities may generate the turbulence when MRI is fully
suppressed (Klahr et al., 2018; Lyra & Umurhan, 2019; Cui & Bai, 2020). If
part of the accretion energy is converted into heat in optically-thick
regions, the disk temperature can be increased by the blanketing effect. It is
important to note that even if the accretion is driven by the disk wind rather
than the turbulence, the heating by the turbulence can still warm the disk up.

There are at least two conditions for the heating to affect the disk
temperature. The first is the growth of hydrodynamic instabilities in
irradiated disks suggested by the MHD model. Pfeil & Klahr (2019) investigated
the unstable regions for some linear hydrodynamic instabilities (i.e.,
vertical convective instability, vertical shear instability, and convective
overstability). The disk midplane is stable to convective overstability and
vertical convection at 1 au in the absence of viscous heating (see the lowest
$\alpha$ case in their Figure 10). In addition, the vertical shear instability
does not grow in the optically thick region. Therefore, the linear
instabilities may not be a heating source in the optically thick inner region
of MHD accretion disks. On the other hand, some nonlinear instabilities (e.g.,
zombie vortex instability and subcritical baroclinic instability) may grow
even in irradiated disks (Marcus et al., 2015; Lesur & Latter, 2016; Lyra et
al., 2019; Pfeil & Klahr, 2019). The nonlinear instabilities may drive
turbulence if some mechanisms provide sufficient amplitudes.

The second condition is sufficiently strong turbulence around the midplane to
affect the disk temperature compared to irradiation. Hydrodynamic simulations
have shown that the Shakura-Sunyaev $\alpha$ parameter showing the turbulence
strength averaged over the disk height is $\sim 10^{-4}$–$10^{-3}$ (e.g.,
Lyra, 2014; Stoll & Kley, 2014; Marcus et al., 2015; Flock et al., 2020).
Assuming the presence of the uniform turbulent viscosity in the MHD accretion
disk, the $\alpha$ parameter of $\approx 3\times 10^{-4}$ is sufficient to
affect the disk temperature profile around 1 au at 1 Myr. However, it is still
unclear whether the turbulent viscosity is uniform and enough heat is released
inside the disk.

We should note that once any turbulence changes the temperature structure, the
turbulence may remain by sustaining the unstable thermal structure.
Instability growing in the upper region could bring the turbulence near the
midplane (Nelson et al., 2013; Klahr & Hubbard, 2014). Such turbulence may
also transport the energy to around the midplane. Further numerical studies
are necessary for understanding practical criteria for the disk heating.

To summarize this section, hydrodynamic turbulence is a potential disk heating
source but it is still unclear whether such turbulence develops and whether it
releases enough heat around the midplane.

####  4.5.2 Shock Heating of Waves Induced by Gravitational Instability

Disks at early stages may be subject to the gravitational instability. Shock
heating of waves induced by self-gravitational instability can also warm up
the inner disk region (see Boss & Durisen, 2005; Martin & Livio, 2012;
Rafikov, 2016). This instability typically occurs in the outer disk region,
but shock waves excited there may propagate to the inner region (Rafikov,
2016). The shock waves release heat around the midplane (Hirose & Shi, 2017),
and thus can warm up the disk by the blanketing effect.

However, the unstable region at 0.6 Myr is uncertain but is generally expected
to be well beyond $\sim$10 au (see Kratter & Lodato, 2016), which would be too
far away for the wave to propagate to $\sim$ 1 au. Thus, the gravitational
instability would not drive the disk heating when the snow line passes the
rocky-planet-forming region.

####  4.5.3 Shock Heating of Planetary Wakes

Waves induced by a planet can also warm up the inner disk region as well as
the gravitational instability (Lyra et al., 2016; Rafikov, 2016; Ziampras et
al., 2020). Ziampras et al. (2020) showed that in an optically thick disk, the
shock heating of wakes generated by a (sub-)Jupiter mass planet warms up the
inner region and thereby shifts the snow line outward. For instance, when
$\dot{M}=10^{-8}\rm\,\rm\,M_{\sun}\,\rm\,yr^{-1}$ and $\alpha=10^{-3}$, a
Jupiter-mass planet at $r=4$ au provides enough heat to raise the temperature
at 1 au to 300 K (see Figure 6 in their paper).

The important point is that the perturber must have sufficient mass. According
to Ziampras et al. (2020), the perturber with $\gtrsim$ 100 M⊕ provides heat
that increases the temperature to $\gtrsim$ 200 K at $\sim$ 1 au. Therefore,
the proto-Jupiter needs to have formed early in order for this heating to
affect our conclusion that the snow line would have passed 1 au within 1 Myr.
Interestingly, if Jupiter formed early, the oversupply of water to the
terrestrial planets may also be resolved (see Section 4.2.1).

###  4.6 Effects of the Dust Model on Accretion Heating

In this paper, we have assumed that the dust model is as 0.1 $\mu$m-sized
grains and that the opacity is constant. As seen in Sections 3.2.1 and 3.2.2,
the efficiency of accretion heating in the MHD model depends on the ionization
fraction and opacity. We discuss how the dust model can be changed, and how
the change will affect the efficiency of accretion heating.

####  4.6.1 Dust Evolution

Coagulation of dust grains alters the dust spatial and size distributions.
Especially, the abundance of smaller grains in the disk upper layer mainly
determines the opacity and ionization fraction, and thus $\tau_{\rm heat}$.
The dust growth reduces the abundance of smaller grains (e.g., Birnstiel et
al., 2010, 2011). In addition, grains settle toward the midplane, reducing the
dust abundance in the upper layer (Weidenschilling, 1980).

The opacity in the upper layer is decreased by these effects of the dust
evolution. On the other hand, the effects increase the ionization fraction
because the total grain surface area is decreased. The competition between
increasing opacity and decreasing ionization fraction will determine whether
the dust evolution promotes or further suppresses the disk heating. This will
be investigated in our future work.

####  4.6.2 Sublimation of Ice

We have assumed that the opacity is constant value, but the opacity is
decreased at the snow line by sublimation of ice, which affects the radial
temperature profile. Oka et al. (2011) investigated the effect for the
conventional viscously heated disk model by adopting a dust opacity model of
Miyake & Nakagawa (1993). They showed that the opacity is decreased by a
factor of 3 by sublimation of ice, and consequently the temperature inside the
snow line is decreased.

On the other hand, the total surface area of grains is decreased by the
sublimation, resulting in increasing the temperature in the MHD model (Section
3.2.2). Using the mass fraction of ice and silicate in Miyake & Nakagawa
(1993), the total surface area is decreased by a factor of 9. Considering both
effects, if $\tau_{\rm heat}\propto f_{\rm dg}$, as it is when $f_{\rm dg}$ is
changed from 0.01 to 0.001 (Section 3.2.2), $\tau_{\rm heat}$ is increased by
a factor of $\approx 3$ by the ice sublimation. Thus, in the MHD model, the
temperature may increase inside the snow line, which is different from the
conventional models. In this case, the snow line can be shifted outward by 1.3
times when accretion heating determines the snow line location. Nevertheless,
this effect is negligible when irradiation heating is dominant as in the
fiducial case.

##  5 Summary and Conclusions

We have investigated the migration of the water snow line in PPDs whose
accretion is controlled by laminar magnetic fields, which have been proposed
by various nonideal MHD simulations. MBO19 showed that the accretion heating
is much less efficient than in the conventional viscous disk model. This is
because the heating in the MHD model occurs at a high altitude and also
because a substantial fraction of accretion energy released is removed by the
disk wind. We proposed the empirical model of the disk temperature based on
MBO19 to calculate the snow line location over the disk evolution (see Section
2).

The snow line in the MHD disk model reaches the current Earth’s orbit earlier
than in the viscous disk model. For instance, in our fiducial model, the time
where the snow line in the MHD model passes 1 au is 0.6 Myr, whereas that for
the viscous model is 10 Myr (see Section 3.1). In the MHD model, the migration
of the snow line is mainly driven by the temporal evolution of the stellar
luminosity in the pre-main sequence phase. Our parameter study shows that the
important parameters for efficient accretion heating are the disk ionization
level and opacity (see Section 3.2). High opacity and ionization fraction lead
to efficient accretion heating. This indicates that treating the opacity and
ionization fraction consistently using an identical dust model is crucial for
properly evaluating the efficiency of Joule heating. When the opacity is
scaled with the dust-to-gas ratio, effects of these parameters on the
temperature will be cancelled out, and thus the accretion heating is
inefficient, which is consistent with MBO19. Besides, considering that the
present results are based on the maximum estimate of the disk temperature, we
would expect that the accretion heating is not dominant. Thus, the snow line
in the magnetically accreting laminar disks would pass inside the current
Earth’s orbit within 1 Myr after star formation, while the time in the viscous
disk model is much longer than 1 Myr.

The early arrival of the snow line to the current rocky planet orbits
constrains the formation process of the planets (see Section 4.2.1). The
terrestrial planets in the solar system should have formed inside the snow
line because the water contents are significantly lower than that of icy
bodies in the outer solar system. If we assume that the terrestrial planets
formed at the current orbits, the protoplanets should have formed in the early
phase of the disk evolution ($t<1$ Myr). On the other hand, if we consider the
possibility that the protoplanets moved after formation, they could have
formed near the star and then moved outward to the current orbits.

We should note that the present results are based on the assumption that the
magnetic field drives both the disk accretion and heating. If hydrodynamic
instabilities generate the turbulence around the midplane, its energy
dissipation can further warm the disk up. Further research is needed to assess
the validity of this possibility.

The authors thank the anonymous referee for their helpful comments. The
authors also thank Shigenobu Hirose, Ogihara Masahiro, Takahiro Ueda, Yuhito
Shibaike, Shigeru Ida and Ryota Fukai for fruitful discussion and meaningful
comments. This work was supported by JSPS KAKENHI Grant Numbers JP16K17661,
JP17J10129, JP18H05222, JP18H05438, JP19K03926, JP19K03941, JP20H00182,
JP20H00205, JP20H01948, and 21J00086, the National Key R&D Program of China
(No. 2019YFA0405100), and the joint Tohoku-Tsinghua Collaborative Research
Grant.

##  Appendix A Vertical Temperature Structure in Internally Heated Disks

Here we derive the vertical temperature profile of an internally heated gas
disk. For simplicity, we assume disk structure that is symmetric with respect
to $z=0$, and thus half of the heating energy is radiated from one side of the
disk surfaces. Under this assumption, the vertical temperature profile $T_{\rm
acc,\,MHD}(z)$ for given heating rate per unit volume $q(z)$ can be written as
(Hubeny, 1990, MBO19),

$T_{\rm acc}(z)=\left(\frac{3\Gamma}{8\sigma}\right)^{1/4}\left[\tau_{\rm
eff}(z)+\frac{1}{\sqrt{3}}+\frac{2q(z)}{3\rho(z)\kappa_{\rm
P}(z)\Gamma}\right]^{1/4},$ (A1)

where

$\Gamma=\int_{-\infty}^{+\infty}q(z)\,{\rm d}z,$ (A2)

is the heating rate per unit area,

$\tau_{\rm
eff}(z)=\frac{2}{\Gamma}\int_{z}^{+\infty}\rho(z^{\prime})\kappa_{\rm
R}(z^{\prime})\mathcal{F}(z^{\prime})\,{\rm d}z^{\prime},$ (A3)

is the optical depth from $z^{\prime}=z$ to $z^{\prime}=\infty$ weighted by
the radiative energy flux ${\cal F}$, and $\kappa_{\rm R}$ ad $\kappa_{\rm P}$
are the Rosseland and Planck mean opacities, respectively. The radiative
energy flux is related to $q$ as $\partial{\cal F}/\partial z=q$; assuming
that $q(z)$ is symmetric about the disk midplane, we have

$\mathcal{F}(z)=\int_{0}^{z}q(z^{\prime})\,{\rm d}z^{\prime}.$ (A4)

The third term in the brackets of Equation (A1) is negligible except at
locations where $\tau_{\rm eff}(z)$ is small and $q(z)$ is anomalously high
(e.g., at current sheets well above and below the midplane). If the heating
rate scales with density as in the standard viscous model,
$q(z)=(\Gamma/\Sigma)\rho(z)$, and hence the third term becomes
$2/(3\kappa_{\rm P}\Sigma)$. When the Planck-mean optical depth at the
midplane $\kappa_{\rm P}\Sigma$ is $\gg 1$, this term is much smaller than the
second term for all $z$.

Neglecting the third term in the right-hand side of Equation (A1), the
temperature at the midplane can be approximately written as

$T_{\rm acc}(z=0)=\left(\frac{3\Gamma}{8\sigma}\right)^{1/4}\left(\tau_{\rm
heat}+\frac{1}{\sqrt{3}}\right)^{1/4},$ (A5)

where we have defined

$\tau_{\rm heat}\equiv\tau_{\rm eff}(z=0).$ (A6)

If we further assume that $\kappa_{\rm R}$ is constant, $\tau_{\rm heat}$ can
be rewritten as $\tau_{\rm heat}=\kappa_{\rm R}\Sigma_{\rm heat}$, where

$\Sigma_{\rm
heat}\equiv\frac{2}{\Gamma}\int_{0}^{\infty}\rho(z^{\prime})\mathcal{F}(z^{\prime})\,{\rm
d}z^{\prime}.$ (A7)

When accretion heating occurs only above an altitude $z_{\rm heat}$,
$\Sigma_{\rm heat}$ approximately represents the mass column depth to
$z=z_{\rm heat}$ because ${\cal F}(z)=0$ at $|z|<z_{\rm heat}$. In the extreme
case where the region of accretion heating is infinitesimally thin,
$\Sigma_{\rm heat}$ is exactly equal to the depth to the layer (see Appendix A
of MBO19). For these cases, one can show that $T_{\rm acc}(z)$ is constant at
$|z|<z_{\rm heat}$.

When accretion heating is vertically uniform in the sense that
$q(z)=(\Gamma/\Sigma)\rho(z)$, we have ${\cal F}=(\Gamma/\Sigma)\chi(z)$ with
$\chi(z)\equiv\int_{0}^{z}\rho(z^{\prime})dz^{\prime}$, which gives
$\Sigma_{\rm
heat}=(2/\Sigma)\int_{0}^{\infty}\chi(z)\rho(z)dz=(2/\Sigma)\int_{0}^{\Sigma/2}\chi
d\chi=\Sigma/4$ and

$\tau_{\rm heat}=\tau_{\rm mid}/2.$ (A8)

##  Appendix B Energy Balance in Wind-Driven Accretion Disks

The energy balance in wind-driven accretion disks is expressed as (Suzuki et
al., 2016, see the equation above their Equation (B.10))

$\displaystyle L_{z}+\Gamma$ $\displaystyle=$ $\displaystyle r\Omega w_{\rm
wind},$ (B1)

where $L_{z}$ is the energy flux carried away by the wind material and $w_{\rm
wind}$ is the $z\phi$ component of the total Maxwell stress exerted by the
wind on both sides of the disk. Here we assume that the disk accretion is
driven only by the angular momentum transport to the disk wind, and neglect
the energy liberated by the radial stress. The right-hand side of the Equation
shows the total liberated energy in the wind-driven accretion disks, while the
left-hand side shows how the liberated energy is used. The wind stress $w_{\rm
wind}$ is related to the wind-driven accretion rate $\dot{M}$ as

$\dot{M}=\frac{4\pi r}{\Omega}w_{\rm wind}.$ (B2)

Eliminating $w_{\rm wind}$, we obtain

$L_{z}+\Gamma=\frac{\Omega^{2}\dot{M}}{4\pi}.$ (B3)

Now we define

$f_{\rm heat}=\frac{\Gamma}{L_{z}+\Gamma}$ (B4)

as the fraction of heat deposited inside the disk in the total liberated
energy. In wind-driven accretion disks, using Equation (B3), we have

$\Gamma=f_{\rm heat}\frac{\Omega^{2}\dot{M}}{4\pi}.$ (B5)

The fraction $f_{\rm heat}$ depends on the polarity of the net vertical
magnetic field. The ranges of $f_{\rm heat}$ in the MHD simulations of MBO19
are 0.02–0.5 and 5$\times 10^{-4}$–0.1 when the vertical magnetic fields are
aligned and anti-aligned with the disk rotation axis, respectively. Note that
the total liberated energy includes the contribution from the radial stress,
though it is not dominant. How much the energy is dissipated correlates with
the toroidal field strength. Hall effect amplifies the toroidal field for the
aligned field case, while Hall effect damps the toroidal field for the anti-
aligned case.

##  Appendix C Depth of the Heating Layer

Figure 12:  The flux-weighted column density $\Sigma_{\rm heat}$ (see Equation
(A7)), against the column depth $\Sigma_{\rm Am=0.3}$ at the altitude where
${\rm Am}=0.3$, both of which are normalized by the column depth at the
midplane, $\Sigma_{\rm mid}=\Sigma/2$. The data are given from the MHD
simulation (MBO19) and additional simulations for $f_{\rm dg}=0.01$ and
various X-ray models. Circles and triangles correspond to the cases where the
background magnetic field is aligned and anti-aligned with the disk rotation
vector, respectively. The solid, dashed, dash-dotted, and dotted lines show
$f_{\rm depth}=$ 1, 0.1, 0.01, and 0.001, respectively.

We here show that $\Sigma_{\rm Am=0.3}$ gives the maximum $\Sigma_{\rm heat}$.
We use the results of nonideal MHD simulations by MBO19 that take into account
nonideal MHD effects (Ohmic diffusion, Hall effect, and ambipolar diffusion).
The resistivities corresponding to the three nonideal effects are computed
from the ionization balance in the disk (see Section 2.1 in MBO19 and Section
2.3 for the detail of the ionization model). We also use data from new
simulations with $f_{\rm dg}=0.01$ and with different X-ray models, where we
vary the X-ray temperature between $\\{3,5\\}$ keV and the X-ray luminosity
between $\\{0.3,1,3\\}\times 10^{30}$ erg s-1. All simulations are carried out
for two cases where the vertical magnetic field $B_{z}$ threading the disk
(i.e., background magnetic field) is set to parallel/anti-parallel to the disk
rotation axis ($B_{z}>0$ and $B_{z}<0$, respectively), because Hall effect
affects the magnetic field profile depending on the polarity of the background
field.

Figure 12 shows $\Sigma_{\rm heat}$ versus $\Sigma_{\rm Am=0.3}$ for all
simulation data. Both axes are normalized by the midplane column depth
$\Sigma_{\rm mid}=\Sigma/2$. Different data points correspond to simulations
with different parameter sets. We also define $f_{\rm depth}=\Sigma_{\rm
heat}/\Sigma_{\rm Am=0.3}$ (Equation (11)) and plot lines of $f_{\rm
depth}=1,0.1,0.01$ and $0.001$. We find that all except one data point are
below the line $f_{\rm depth}=1$. The outlier corresponds to the simulation
for $r=5$ au and $B_{z}>0$. In this particular simulation, the midplane is
ionized enough to sustain a strong current sheet (Bai, 2015) that gives
substantial heating. In this paper, we mainly focus on the region where
$r\lesssim 5~{}\rm au$, where ${\rm Am}\lesssim 0.1$ at the midplane. For this
case, we can safely use $\Sigma_{\rm Am=0.3}$ as the upper limit of
$\Sigma_{\rm heat}$.

Owing to the property of the Hall effect, the actual value of $f_{\rm depth}$
strongly depends on the direction of the background magnetic field relative to
the disk rotation axis. When $B_{z}>0$, the range of $f_{\rm depth}$ is 1–0.1,
which means that the current layers lie approximately just above the ambipolar
dead zone. This is consistent with the expectation that the ambipolar
diffusion determines the altitude of the current layer (see Section 2.3). When
$B_{z}<0$, $f_{\rm depth}=$ 0.1–0.01. This is because the suppression of the
magnetic field by the Hall effect acts to push the heating layer to higher
altitude. Thus, $f_{\rm depth}=1$ provides the upper limit of the midplane
temperature determined by MHD accretion heating.

## References

  * A’Hearn et al. (2011) A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396
  * Bai (2013) Bai, X.-N. 2013, ApJ, 772, 96
  * Bai (2014) —. 2014, ApJ, 791, 137
  * Bai (2015) —. 2015, ApJ, 798, 84
  * Bai (2016) —. 2016, ApJ, 821, 80
  * Bai (2017) —. 2017, ApJ, 845, 75
  * Bai & Goodman (2009) Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737
  * Bai & Stone (2013a) Bai, X.-N., & Stone, J. M. 2013a, ApJ, 769, 76
  * Bai & Stone (2013b) —. 2013b, ApJ, 767, 30
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  * Baraffe et al. (2012) Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118
  * Basilevsky & Head (2003) Basilevsky, A. T., & Head, J. W. 2003, Reports on Progress in Physics, 66, 1699
  * Béthune & Latter (2020) Béthune, W., & Latter, H. 2020, MNRAS, 494, 6103
  * Biersteker & Schlichting (2021) Biersteker, J. B., & Schlichting, H. E. 2021, MNRAS, 501, 587
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79
  * Birnstiel et al. (2011) Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11
  * Boss & Durisen (2005) Boss, A. P., & Durisen, R. H. 2005, ApJ, 621, L137
  * Chandrasekhar (1935) Chandrasekhar, S. 1935, MNRAS, 96, 21
  * Charette & Smith (2010) Charette, M. A., & Smith, W. H. 2010, Oceanography, 23, 112
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
  * Cui & Bai (2020) Cui, C., & Bai, X.-N. 2020, ApJ, 891, 30
  * Desch & Turner (2015) Desch, S. J., & Turner, N. J. 2015, ApJ, 811, 156
  * Donahue et al. (1982) Donahue, T. M., Hoffman, J. H., Hodges, R. R., & Watson, A. J. 1982, Science, 216, 630
  * Drążkowska & Dullemond (2018) Drążkowska, J., & Dullemond, C. P. 2018, A&A, 614, A62
  * Drążkowska et al. (2013) Drążkowska, J., Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37
  * Dzyurkevich et al. (2010) Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70
  * Elkins-Tanton et al. (2007) Elkins-Tanton, L. T., Smrekar, S. E., Hess, P. C., & Parmentier, E. M. 2007, Journal of Geophysical Research (Planets), 112, E04S06
  * Fei et al. (2017) Fei, H., Yamazaki, D., Sakurai, M., et al. 2017, Science Advances, 3, e1603024
  * Feiden (2016) Feiden, G. A. 2016, A&A, 593, A99
  * Fleming & Stone (2003) Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155
  * Fromang (2013) Fromang, S. 2013, in EAS Publications Series, Vol. 62, EAS Publications Series, ed. P. Hennebelle & C. Charbonnel, 95–142
  * Fromang et al. (2013) Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  * Garaud & Lin (2007) Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606
  * Genda & Abe (2005) Genda, H., & Abe, Y. 2005, Nature, 433, 842
  * Gressel et al. (2015) Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84
  * Guilet & Ogilvie (2014) Guilet, J., & Ogilvie, G. I. 2014, MNRAS, 441, 852
  * Guillot (2005) Guillot, T. 2005, Annual Review of Earth and Planetary Sciences, 33, 493
  * Hamano et al. (2013) Hamano, K., Abe, Y., & Genda, H. 2013, Nature, 497, 607
  * Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385
  * Hartmann et al. (1997) Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, ApJ, 475, 770
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135
  * Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  * Hirose & Shi (2017) Hirose, S., & Shi, J.-M. 2017, MNRAS, 469, 561
  * Hirose & Turner (2011) Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30
  * Hosokawa et al. (2011) Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140
  * Hubeny (1990) Hubeny, I. 1990, ApJ, 351, 632
  * Ida et al. (2019) Ida, S., Yamamura, T., & Okuzumi, S. 2019, A&A, 624, A28
  * Igea & Glassgold (1999) Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848
  * Johansen et al. (2021) Johansen, A., Ronnet, T., Bizzarro, M., et al. 2021, arXiv e-prints, arXiv:2102.08611
  * Kimura et al. (2016) Kimura, S. S., Kunitomo, M., & Takahashi, S. Z. 2016, MNRAS, 461, 2257
  * Klahr & Hubbard (2014) Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21
  * Klahr et al. (2018) Klahr, H., Pfeil, T., & Schreiber, A. 2018, Instabilities and Flow Structures in Protoplanetary Disks: Setting the Stage for Planetesimal Formation, ed. H. J. Deeg & J. A. Belmonte, 138
  * Kleine et al. (2020) Kleine, T., Budde, G., Burkhardt, C., et al. 2020, Space Sci. Rev., 216, 55
  * Kratter & Lodato (2016) Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271
  * Kretke et al. (2009) Kretke, K. A., Lin, D. N. C., Garaud, P., & Turner, N. J. 2009, ApJ, 690, 407
  * Kruijer et al. (2017) Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proceedings of the National Academy of Science, 114, 6712
  * Kunitomo et al. (2017) Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A, 599, A49
  * Kurokawa et al. (2014) Kurokawa, H., Sato, M., Ushioda, M., et al. 2014, Earth and Planetary Science Letters, 394, 179
  * Kusaka et al. (1970) Kusaka, T., Nakano, T., & Hayashi, C. 1970, Progress of Theoretical Physics, 44, 1580
  * Lawrence et al. (2013) Lawrence, D. J., Feldman, W. C., Goldsten, J. O., et al. 2013, Science, 339, 292
  * Lesur & Latter (2016) Lesur, G. R. J., & Latter, H. 2016, MNRAS, 462, 4549
  * Leung & Ogilvie (2019) Leung, P. K. C., & Ogilvie, G. I. 2019, MNRAS, 487, 5155
  * Lichtenberg et al. (2021) Lichtenberg, T., Drążkowska, J., Schönbächler, M., Golabek, G. J., & Hands, T. O. 2021, Science, 371, 365
  * Lichtenberg et al. (2019) Lichtenberg, T., Golabek, G. J., Burn, R., et al. 2019, Nature Astronomy, 3, 307
  * Lodders (2003) Lodders, K. 2003, ApJ, 591, 1220
  * Lubow et al. (1994) Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 267, 235
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
  * Lyra (2014) Lyra, W. 2014, ApJ, 789, 77
  * Lyra et al. (2016) Lyra, W., Richert, A. J. W., Boley, A., et al. 2016, ApJ, 817, 102
  * Lyra & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001
  * Lyra et al. (2019) Lyra, W., Haworth, T., Bitsch, B., et al. 2019, BAAS, 51, 129
  * Machida & Abe (2010) Machida, R., & Abe, Y. 2010, ApJ, 716, 1252
  * Marcus et al. (2015) Marcus, P. S., Pei, S., Jiang, C.-H., et al. 2015, ApJ, 808, 87
  * Martin & Livio (2012) Martin, R. G., & Livio, M. 2012, MNRAS, 425, L6
  * Milne (1921) Milne, E. A. 1921, MNRAS, 81, 361
  * Miyake & Nakagawa (1993) Miyake, K., & Nakagawa, Y. 1993, Icarus, 106, 20
  * Morbidelli et al. (2016) Morbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368
  * Mori et al. (2019) Mori, S., Bai, X.-N., & Okuzumi, S. 2019, ApJ, 872, 98
  * Nakamoto & Nakagawa (1994) Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Nomura et al. (2014) Nomura, R., Hirose, K., Uesugi, K., et al. 2014, Science, 343, 522
  * Ogihara et al. (2015) Ogihara, M., Kobayashi, H., Inutsuka, S.-i., & Suzuki, T. K. 2015, A&A, 579, A65
  * Oka et al. (2011) Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141
  * Okuzumi (2009) Okuzumi, S. 2009, ApJ, 698, 1122
  * Okuzumi et al. (2014) Okuzumi, S., Takeuchi, T., & Muto, T. 2014, ApJ, 785, 127
  * Pfeil & Klahr (2019) Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150
  * Pollack et al. (1985) Pollack, J. B., McKay, C. P., & Christofferson, B. M. 1985, Icarus, 64, 471
  * Rafikov (2016) Rafikov, R. R. 2016, ApJ, 831, 122
  * Rotundi et al. (2015) Rotundi, A., Sierks, H., Della Corte, V., et al. 2015, Science, 347, aaa3905
  * Salmeron & Wardle (2008) Salmeron, R., & Wardle, M. 2008, MNRAS, 388, 1223
  * Sasselov & Lecar (2000) Sasselov, D. D., & Lecar, M. 2000, ApJ, 528, 995
  * Sato et al. (2016) Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15
  * Schlichting et al. (2015) Schlichting, H. E., Sari, R., & Yalinewich, A. 2015, Icarus, 247, 81
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Simon et al. (2015) Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117
  * Stahler & Palla (2005) Stahler, S. W., & Palla, F. 2005, The Formation of Stars (Wiley-VCH), 865
  * Stoll & Kley (2014) Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77
  * Suzuki et al. (2016) Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74
  * Takeuchi & Okuzumi (2014) Takeuchi, T., & Okuzumi, S. 2014, ApJ, 797, 132
  * Umebayashi & Nakano (2009) Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69
  * Wardle (2007) Wardle, M. 2007, Ap&SS, 311, 35
  * Weidenschilling (1980) Weidenschilling, S. J. 1980, Icarus, 44, 172
  * Xiao et al. (2017) Xiao, L., Niu, R., & Zhang, H. 2017, MNRAS, 467, 2869
  * Zhang & Jin (2015) Zhang, Y., & Jin, L. 2015, ApJ, 802, 58
  * Ziampras et al. (2020) Ziampras, A., Ataiee, S., Kley, W., Dullemond, C. P., & Baruteau, C. 2020, A&A, 633, A29


# Non-Relativistic Scenario of a Free Particle in an Infinite Cylindrical
Potential Well

Pratik Adarsh  Department of Physics, Indian Institute of Technology Bhilai,
GEC Campus, Sejbahar, Raipur, Chhattisgarh 492015, India  Sabyasachi Ghosh
Department of Physics, Indian Institute of Technology Bhilai, GEC Campus,
Sejbahar, Raipur, Chhattisgarh 492015, India

###### Abstract

There are various types of infinite potential well problems occurring in
elementary quantum mechanics formalism. The infinite square well (one
dimensional), cubical box and, spherical well are quite common in textbooks.
In this paper, we consider a rather uncommon potential well, infinite
cylindrical well, and try to find its energy eigenvalues and eigenfunctions
using Schrödinger equation. These results are also derived in the paper[1].
But our approach here is to elaborate every step just the way Griffiths did in
his book[2] for hydrogen atom, using spherical potential well. We also plot
some radial wavefunctions and density plots.

##  1 Introduction

In infinite spherical potential well, we readily solve the time-independent
Schrödinger equation to get energy eigenvalues and eigenfunctions. In this
problem of infinite cylindrical well, we proceed with the same approach by
solving the time-independent Schrödinger equation in cylindrical coordinates.
The motivation is to see how cylindrical symmetry differs from spherical one
and to what extend its analytical solutions are possible in the most general
way.

##  2 The Problem

Figure 1: Infinite cylindrical well.

We have a particle inside a cylindrical-shaped well for which the potential is
zero inside the boundaries but infinite at and outside the boundaries. The
dimensional attributes of the well are:
Height of the cylinder = $H$
Radius of the cylinder = $a$
With these we have to solve the time-independent Schrödinger equation:

${\frac{-{\hbar}^{2}}{2m}}{{\nabla}^{2}}{\psi}+{V\psi}=E\psi$

##  3 The Approach

###  3.1 Separation of Variables

The Laplacian in cylindrical polar coordinates is:

${\nabla^{2}}={\frac{1}{r}}{\frac{\partial}{\partial
r}}\left({r}\frac{\partial}{\partial
r}\right)+{\frac{1}{r^{2}}}{\frac{\partial^{2}}{\partial\phi^{2}}}+\frac{\partial^{2}}{\partial
z^{2}}$

The wavefunction is zero outside the well, so the Schrödinger equation
inside($V=0$) the well reads:

$\frac{-\hbar^{2}}{2m}\left[\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial\psi}{\partial
r}\right)+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\phi^{2}}+\frac{\partial^{2}\psi}{\partial
z^{2}}\right]=E\psi$ (3.1.1)

We let the solutions separable into products:

$\psi(r,\phi,z)=\psi_{r\phi}(r,\phi)\psi_{z}(z)$

Substituting this into equation 3.1.1,we have

$\frac{-\hbar^{2}}{2m}\left[\frac{\psi_{z}}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial\psi_{r\phi}}{\partial
r}\right)+\frac{\psi_{z}}{r^{2}}\frac{\partial^{2}\psi_{r\phi}}{\partial\phi^{2}}+\psi_{r\phi}\frac{\partial^{2}\psi_{z}}{\partial
z^{2}}\right]=E\psi_{r\phi}\psi_{z}$ (3.1.2)

Rearranging terms, we get

$\underbrace{\frac{1}{\psi_{r\phi}}\left[\frac{1}{r}\frac{\partial\psi_{r\phi}}{\partial
r}+\frac{\partial^{2}\Psi_{r\phi}}{\partial
r^{2}}+\frac{1}{r^{2}}\frac{\partial^{2}\psi_{r\phi}}{\partial\phi^{2}}\right]}_{\text{depends
on }r\text{ and
}\phi}+\underbrace{\frac{1}{\psi_{z}}\left[\frac{d^{2}\psi_{z}}{dz^{2}}\right]}_{\text{depends
on }z}=-\frac{2mE}{\hbar^{2}}$ (3.1.3)

In the L.H.S of equation 3.1.3, the first term depends only on $r$ and $\phi$.
Whereas the second term depends only on $z$. So we let the first term equals
some factor $A_{1}$(say)and the second term equals to some another factor
$A_{2}$. With these, the equation 3.1.3 reads as:

$A_{1}+A_{2}=\frac{-2mE}{\hbar^{2}}$ (3.1.4)

This equation 3.1.4 is the main equation that would give us the energy
eigenvalues for the system. Our aim would be now to find these factors,
$A_{1}$ and $A_{2}$.
Coming back to equation 3.1.3, we write the first term on the L.H.S as:

$\frac{1}{r}\frac{\partial\psi_{r\phi}}{\partial
r}+\frac{\partial^{2}\Psi_{r\phi}}{\partial
r^{2}}+\frac{1}{r^{2}}\frac{\partial^{2}\psi_{r\phi}}{\partial\phi^{2}}=A_{1}\psi_{r\phi}$
(3.1.5)

This is just another partial differential equation we have to solve for
proceeding in the right direction. We let the solutions separable like
previously:

$\psi_{r\phi}(r,\phi)=R(r)X(\phi)$

Substituting this into equation 3.1.5, we have

$\frac{X}{r}\frac{dR}{dr}+X\frac{d^{2}R}{dr^{2}}+\frac{R}{r^{2}}\frac{d^{2}X}{d\phi^{2}}=A_{1}RX$
(3.1.6)

Rearranging terms, we get

$\underbrace{\left[\frac{r}{R}\frac{dR}{dr}+\frac{r^{2}}{R}\frac{d^{2}R}{dr^{2}}-A_{1}r^{2}\right]}_{\text{depends
on
}r}+\underbrace{\left[\frac{1}{X}\frac{d^{2}X}{d\phi^{2}}\right]}_{\text{depends
on }\phi}=0$ (3.1.7)

The first term on L.H.S of the equation 3.1.7 depends only on $r$ and the
second term depends only on $\phi$. So we let the first term equals to some
factor $B$(say), then obviously the second term equals to $-B$.
With these attributes, we have:

$\frac{r}{R}\frac{dR}{dr}+\frac{r^{2}}{R}\frac{d^{2}R}{dr^{2}}-A_{1}r^{2}=B$
(3.1.8) $\frac{1}{X}\frac{d^{2}X}{d\phi^{2}}=-B$ (3.1.9)

And from equations 3.1.3 and 3.1.4

$\frac{1}{\psi_{z}}\left[\frac{d^{2}\psi_{z}}{dz^{2}}\right]=A_{2}$ (3.1.10)

###  3.2 The Radial Equation

We can rewrite equation 3.1.8 as:

$\frac{d}{dr}\left(rR^{\prime}(r)\right)+\left[-A_{1}-\frac{B}{r^{2}}\right]rR(r)=0\textit{
(Strum-Liouville Equation)}$ (3.2.1)

Its solution is:

$R(r)=C_{1}J_{\scaleto{\sqrt{B}}{5pt}}\left(-i\sqrt{A_{1}}r\right)+C_{2}Y_{\scaleto{\sqrt{B}}{5pt}}\left(-i\sqrt{A_{1}}r\right)$
(3.2.2)

Where $J$ and $Y$ are Bessel functions of first and second kind respectively.
Equation 3.2.2 gives us the radial wavefunction. We have to solve it with
appropriate boundary conditions. We would do this later, now the main task is
to simplify it with this simple view:
For $r\,\to\,0$, equation 3.2.2 reads as:

$R(0)=C_{1}J_{\scaleto{\sqrt{B}}{5pt}}\left(0\right)+C_{2}Y_{\scaleto{\sqrt{B}}{5pt}}\left(0\right)$

But we know that

$Y_{\scaleto{\sqrt{B}}{5pt}}\left(0\right)\,\to\,\infty$
$J_{\scaleto{\sqrt{B}}{5pt}}\left(0\right)=0$

This means that the radial wavefunction, $R(0)$ diverges on the axis of the
well. In other words this implies that it’s a non-normalizable wavefunction.
Therefore we must have $C_{2}=0$.
So now our radial wavefunction, equation 3.2.2 becomes:

$R(r)=C_{1}J_{\scaleto{\sqrt{B}}{5pt}}\left(-i\sqrt{A_{1}}r\right)$ (3.2.3)

We will return to this equation soon for calculating $A_{1}$, to plug in the
equation 3.1.4 and getting the energy eigenvalues. But first, we have to find
what $B$ corresponds and of course, we also have to find $A_{2}$.

###  3.3 The Angular Equation

The equation 3.1.9 is:

$\frac{1}{X}\frac{d^{2}X}{d\phi^{2}}=-B$

Its solution can be written as:

$X(\phi)=C_{3}\exp\left(i\sqrt{B}\phi\right)+C_{4}\exp\left(-i\sqrt{B}\phi\right)$
(3.3.1)

For rotational symmetry, we must have:

$X(\phi+2\pi)=X(\phi)$ $\implies
C_{3}\exp\left(i\sqrt{B}(\phi+2\pi)\right)+C_{4}\exp\left(-i\sqrt{B}(\phi+2\pi)\right)=C_{3}\exp\left(i\sqrt{B}\phi\right)+C_{4}\exp\left(-i\sqrt{B}\phi\right)$

Simplifying and equating real and imaginary parts, we get:

$\left(C_{3}+C_{4}\right)\left[\cos(\sqrt{B}\phi+2\pi\sqrt{B})\right]=\left(C_{3}+C_{4}\right)\cos(\sqrt{B}\phi)$

and

$\left(C_{3}-C_{4}\right)\left[\sin(\sqrt{B}\phi+2\pi\sqrt{B})\right]=\left(C_{3}-C_{4}\right)\sin(\sqrt{B}\phi)$

Therefore $\sqrt{B}$ is an integer, for $C_{3}+C_{4}\neq 0$ and $C_{3}\neq
C_{4}$. Hence we have

$\sqrt{B}=n_{\phi}=0,1,2,...$ (3.3.2)

where $n_{\phi}$ denotes a quantum number associated with azimuthal angle.

###  3.4 The Equation along Z Direction

The equation 3.1.10 is:

$\frac{1}{\psi_{z}}\left[\frac{d^{2}\psi_{z}}{dz^{2}}\right]=A_{2}$

Its solution can be written as:

$\psi_{z}(z)=C_{5}\exp\left(iz\sqrt{\frac{2mE}{\hbar^{2}}+A_{1}}\right)+C_{6}\exp\left(-iz\sqrt{\frac{2mE}{\hbar^{2}}+A_{1}}\right)$
(3.4.1)

Now we have the following boundary conditions:

  1. (i.)

$\psi_{z}(0)=0$

  2. (ii.)

$\psi_{z}(H)=0$

Imposing the first condition on the equation 3.4.1, we have:

$C_{6}=-C_{5}$

With this and by second boundary condition, the equation 3.4.1 becomes:

$C_{5}\exp\left(iz\sqrt{\frac{2mE}{\hbar^{2}}+A_{1}}\right)-C_{5}\exp\left(-iz\sqrt{\frac{2mE}{\hbar^{2}}+A_{1}}\right)=0$
$\implies
2iC_{5}\left[\sin\left(H\sqrt{\frac{2mE}{\hbar^{2}}+A_{1}}\right)\right]=0$
$\implies H\sqrt{\frac{2mE}{\hbar^{2}}+A_{1}}=n_{z}\pi\hskip
5.0pt;n_{z}=1,2,3,...$ (3.4.2)

Where $n_{z}$ denotes a quantum number associated with z-direction.
Using equations 3.4.2 and 3.1.4, we have:

$A_{2}=\frac{-n_{z}^{2}\pi^{2}}{H^{2}}$ (3.4.3)

Therefore we have energy expression redefined as:

$E=\frac{\hbar^{2}(n_{z}^{2}\pi^{2}-A_{1}H^{2})}{2mH^{2}}$ (3.4.4)

###  3.5 Energy Eigenvalues

With some remaining calculations, we are now very near to get the energy
eigenvalues.
It’s now time to return to equation 3.2.3 with the consideration of equation
3.3.2, we have

$R(r)=C_{1}J_{\scaleto{n_{\phi}}{5pt}}\left(-i\sqrt{A_{1}}r\right)$ (3.5.1)

Imposing the boundary condition

$R(a)=0$

We get

$\implies R(a)=C_{1}J_{\scaleto{n_{\phi}}{5pt}}\left(-i\sqrt{A_{1}}a\right)=0$
$\implies A_{1}=\frac{-(j_{\scaleto{n_{\phi},n_{r}}{5pt}})^{2}}{a^{2}}\hskip
5.0pt;n_{r}=1,2,3,...$ (3.5.2)

Where $n_{r}$ denotes a quantum number associated with radial solution and
$j_{\scaleto{n_{\phi},n_{r}}{5pt}}$ denotes $n_{r}^{th}$ root of Bessel
function $J_{\scaleto{n_{\phi}}{5pt}}$
Putting this equation 3.5.2 into equation 3.4.4, we have the final expression
for energy eigenvalues as:

$\boxed{E_{\scaleto{n_{r},n_{\phi},n_{z}}{5pt}}=\frac{\hbar^{2}\left[n_{z}^{2}a^{2}\pi^{2}+(j_{\scaleto{n_{\phi},n_{r}}{5pt}})^{2}H^{2}\right]}{2ma^{2}H^{2}}}$
(3.5.3)

Where

$\displaystyle n_{r}$ $\displaystyle=1,2,3,...$ $\displaystyle n_{\phi}$
$\displaystyle=0,1,2,...$ $\displaystyle n_{z}$ $\displaystyle=1,2,3,...$

###  3.6 The Eigenfunctions

We started with separable solutions of the Schrödinger’s equation 3.1.1 as:

$\psi(r,\phi,z)=\psi_{r\phi}(r,\phi)\psi_{z}(z)=R(r)X(\phi)\psi_{z}(z)$
(3.6.1)

Now we will find $R(r)$, $X(\phi)$ and $\psi_{z}(z)$ separately and put them
back into the equation 3.6.1 to get the complete eigenfunctions.

The Normalization Technique

Before finding the wavefunctions, we emphasize that how normalization can be
done separately for each of them.
In cylindrical coordinates, we have elementary volume as:

$dV=rdrd\phi dz$

So normalization of wavefunction in the view of equation 3.6.1 becomes

$\iiint_{V}\psi(r,\phi,z)\psi^{*}(r,\phi,z){rdrd\phi
dz}=\int_{0}^{a}R(r)R^{*}(r)rdr\int_{0}^{2\pi}X(\phi)X^{*}(\phi)d\phi\int_{0}^{H}\psi_{z}(z)\psi_{z}^{*}(z)dz=1$

We can normalize each wavefunctions separately as follows:

  1. (i.)

$\int_{0}^{a}R(r)R^{*}(r)rdr=1$ (3.6.2)

  1. (ii.)

$\int_{0}^{2\pi}X(\phi)X^{*}(\phi)d\phi=1$ (3.6.3)

  1. (iii.)

$\int_{0}^{H}\psi_{z}(z)\psi_{z}^{*}(z)dz=1$ (3.6.4)

The Radial wavefunctions

From the equations 3.5.1 and 3.5.2, we have:

$R(r)=C_{1}J_{\scaleto{n_{\phi}}{5pt}}\left[\frac{(j_{\scaleto{n_{\phi},n_{r}}{5pt}})r}{a}\right]$
(3.6.5)

To make further calculations easier to follow, we take the following
denotations:

$n_{\phi}=n\text{ and }j_{\scaleto{n_{\phi},n_{r}}{5pt}}=\beta$

With these, the radial equation 3.6.5 becomes

$R(r)=C_{1}J_{n}\left[\frac{\beta{r}}{a}\right]$ (3.6.6)

Normalizing using equation 3.6.2, we have

$C_{1}^{2}\int_{0}^{a}J_{n}^{2}\left[\frac{\beta{r}}{a}\right]rdr=1$ (3.6.7)

The indefinite integral has solution as:

$\int{J_{n}^{2}\left[\frac{\beta{r}}{a}\right]rdr}=\frac{r^{2}}{2}\left[J_{n}^{2}\left[\frac{\beta{r}}{a}\right]-J_{n-1}\left[\frac{\beta{r}}{a}\right]J_{n+1}\left[\frac{\beta{r}}{a}\right]\right]+C$
(3.6.8)

where $C$ is the constant of integration.
With the proper limits of integration, $[0,a]$ and using Bessel recurrence
relation:

$J_{n-1}(x)+J_{n+1}(x)=2\frac{n}{x}J_{n}(x)$ (3.6.9)

we have

$C_{1}=\frac{\sqrt{2}}{aJ_{\scaleto{n+1}{5pt}}\left(\beta\right)}$

Putting back into equation 3.6.5 and restoring the original denotations, we
get the normalized radial wavefunctions as:

$\boxed{R(r)=\frac{\sqrt{2}}{aJ_{\scaleto{n_{\phi}+1}{5pt}}\left[j_{\scaleto{n_{\phi},n_{r}}{5pt}}\right]}J_{\scaleto{n_{\phi}}{5pt}}\left[\frac{(j_{\scaleto{n_{\phi},n_{r}}{5pt}})r}{a}\right]}$
(3.6.10)

The Angular wavefunctions

The angular wavefunctions are given by equation 3.3.1, they are basically two
linearly independent angular solutions. But we can compress it to the form:

$X(\phi)=C_{p}\exp(ip\phi)\text{; where }p=\pm\sqrt{B}=\pm n_{\phi}=0,\pm
1,\pm 2,...$ (3.6.11)

The $p$ is like another quantum number that would be responsible for state
degeneracy. We can normalize 3.6.11 in the view of the equation 3.6.3, we get
the normalized angular wavefunctions as:

$\boxed{X(\phi)=\frac{1}{\sqrt{2\pi}}\exp(ip\phi)\text{ ; }p=\pm
n_{\phi}=0,\pm 1,\pm 2,...}$ (3.6.12)

The wavefunctions along Z-Direction

From subsection 3.4, we have wavefunctions along z-direction as:

$\psi_{z}(z)=2iC_{5}\left[\sin\left(z\sqrt{\frac{2mE}{\hbar^{2}}+A_{1}}\right)\right]$
(3.6.13)

Using the equation 3.4.4, we can rewrite equation 3.6.13 as:

$\psi_{z}(z)=2C_{5}^{\prime}\sin\left(\frac{n_{z}\pi{z}}{H}\right)$ (3.6.14)

Normalizing in the view of equation 3.6.4, we get the normalized wavefunctions
along the z-direction as:

$\boxed{\psi_{z}(z)=\left[\sqrt{\frac{2}{H}}\right]\sin\left(\frac{n_{z}\pi{z}}{H}\right)\text{
; }n_{z}=1,2,3,...}$ (3.6.15)

Now using equation 3.6.1, we stack up all the wavefunctions and get the final
expression as:

$\boxed{\psi_{\scaleto{n_{r},n_{\phi},n_{z},p}{5pt}}=\left[\sqrt{\frac{2}{\pi{Ha^{2}}}}\right]\frac{1}{J_{\scaleto{n_{\phi}+1}{5pt}}\left[j_{\scaleto{n_{\phi},n_{r}}{5pt}}\right]}J_{\scaleto{n_{\phi}}{5pt}}\left[\frac{(j_{\scaleto{n_{\phi},n_{r}}{5pt}})r}{a}\right]\sin\left[\frac{n_{z}\pi{z}}{H}\right]\exp(ip\phi)}$
(3.6.16)

Where

$\displaystyle n_{r}$ $\displaystyle=1,2,3,...$ $\displaystyle n_{\phi}$
$\displaystyle=0,1,2,...$ $\displaystyle n_{z}$ $\displaystyle=1,2,3,...$
$\displaystyle p$ $\displaystyle=\pm n_{\phi}$

##  4 Plots

###  4.1 Radial Wavefunctions

Figure 2: Some Radial Wavefunctions $R_{\scaleto{n_{r},n_{\phi}}{5pt}}$
(3.6.10).

###  4.2 Probability density variation,
${\lvert\psi_{\scaleto{n_{r},n_{\phi},n_{z}}{5pt}}\rvert}^{2}$

Figure 3: ${\lvert\psi_{\scaleto{1,0,1}{5pt}}\rvert}^{2}$

Figure 4: ${\lvert\psi_{\scaleto{2,0,1}{5pt}}\rvert}^{2}$

Figure 5: ${\lvert\psi_{\scaleto{1,1,1}{5pt}}\rvert}^{2}$

Figure 6: ${\lvert\psi_{\scaleto{2,2,2}{5pt}}\rvert}^{2}$

Comments

In 4.2 we don’t have that extra quantum number $p$, as we had in equation
(3.6.16). This is obvious as probability density $=\psi\psi^{*}$.

In fig. 5 and fig. 6, we have zero probability density on the z-axis. This is
evident from the term in (3.6.16):

$J_{\scaleto{n_{\phi}}{5pt}}\left[\frac{(j_{\scaleto{n_{\phi},n_{r}}{5pt}})r}{a}\right]$

On the z-axis, $r=0$. So the term becomes

$J_{\scaleto{n_{\phi}}{5pt}}\left(0\right)$

But

$J_{n}\left(0\right)=0\text{ if }n\neq 0$ $J_{n}\left(0\right)=1\text{ if
}n=0$

That is how plots in fig. 5 and fig. 6 differ from plots in fig. 3 and fig. 4.
These observations can also, be verified with the radial wavefunction plots
(4.1).
In fig. 6, however, we have an additional place (apart from the z-axis) where
probability density is vanishing. This can be verified with radial
wavefunction plot $R_{\scaleto{n_{r},n_{\phi}}{5pt}}=R_{\scaleto{2,2}{5pt}}$
in (4.1).

##  5 Conclusions

We can see from energy eigenvalues expression, equation 3.5.3, that the energy
of a free particle in an infinite cylindrical well depends on three quantum
numbers associated with three independent directions. But the eigenfunctions,
given by equation 3.6.16, depend on an extra quantum number $p$. For each
$n_{\phi}$, there are two values of $p$. Thus all excited states are 2-fold
degenerate.

## References

  * [1] Arkadiy S. Baltenkov and Alfred Z. Msezane.  Electronic quantum confinement in cylindrical potential well.  The European Physical Journal D, 70(4), Apr 2016.  URL: http://dx.doi.org/10.1140/epjd/e2016-60728-2, doi:10.1140/epjd/e2016-60728-2.
  * [2] David J Griffiths.  Introduction to quantum mechanics.  Prentice Hall, second edi, 2010.