{'fraction_non_alphanumeric': 0.03288263926446728, 'fraction_numerical': 0.022841175410131604, 'mean_word_length': 4.360034785969659, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evolution of the Water Snow Line in Magnetically Accreting Protoplanetary Disks

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.


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 ∼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 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 1 wt%. For the Earth, ∼ 0.1M ⊕ 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 ∼ 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 ∼ 50 wt% of water ice (Lodders 2003). Planetesimals forming from icy dust can lose water to some extent through radiogenic heating by 26 Al, 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 26 Al (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 di-rect stellar irradiation, the snow line lies at ≈ 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 cospatial 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 model 1 . 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 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 1 Bé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).  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 1000 K; Desch & Turner 2015) where MRIturbulence 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 ∼ 160 K.

Density Structure and Accretion Rate
We assume that the disk is vertically nearly isothermal and give the vertical gas density profile as where z is the distance from the midplane, H is the scale height, and Σ is the gas surface density. The scale height is related to the isothermal sound speed c s and local Keplerian frequency Ω as H = c s /Ω. 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Ṁ can be written as (Wardle 2007;Fromang et al. 2013) where r is the radial distance, and w wind is the zφ component of the Maxwell stress exerted on both sides of the disk, which we call the wind stress.
In principle, w 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 wind as (Suzuki et al. 2016) where P mid is the gas pressure at the midplane and α zφ 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 α zφ range between 10 −4 -10 −2 . We take α zφ = 10 −3 as a default value and vary α zφ in Section 3.2.3. The parameter α zφ should not be confused with the α 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 Σ in steady state is related toṀ as Σ ∼Ṁ/(α zφ c s r) for wind-driven accretion (see Equation (4) below) and Σ ∼ M/(αc s H) for viscous accretion. Comparing the two expressions, the two disk models give the same value of Σ for giveṅ M when α ∼ (r/H)α zφ (Wardle 2007;Bai & Goodman 2009;Fromang 2013;Bai 2017).
Assuming wind-driven accretion with α zφ being constant, one can relate Σ toṀ. When the disk is isothermal except at high altitude, the midplane pressure can be written as P mid = Σc 2 s /( and (3), we have We assume that the disk is in a quasi-steady state, which follows from the assumption thatṀ 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Ṁ is radially constant (see Suzuki et al. 2016;Bai 2016) and is equal to the gas accretion rate onto the central star. We determineṀ as a function of stellar age t using the empirical relation for the stellar mass M = 1 M by Hartmann et al. (2016) 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.

Disk Heating
We consider both stellar irradiation and accretion heating and give the temperature T in the disk interior as where T irr and T 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) T irr = 110 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 M star by Feiden (2016), as shown in Figure  3. We note that, in their model, a 1 M 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).
For T acc, MHD at the midplane, we use (see Appendixes A and B for the derivation) where σ is the Stefan-Boltzmann constant, and τ 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)  inside the disk, with f heat < 1 when the wind material carries away part of the energy generated by accretion. We set the default value of f heat to 1 to give a maximum estimate of T acc, MHD , where f heat = 1 is the limiting case that all the liberated energy is consumed for heating. We discuss the effects of varying f heat in Section 3.2.3. We assume that the Rosseland mean opacity κ R is constant throughout the disk and write where Σ 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 Σ heat ). The remaining task is to evaluate Σ heat .

Estimating the Heating Layer Depth from Disk Ionization Structure
We estimate Σ 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 where v A is the Alfvén speed and η A is the ambipolar diffusivity. The ambipolar Elsasser number is in general independent of the magnetic field strength because both v 2 A and η A scale quadratically with it. In inner regions of protoplanetary disks, 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 Am 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 Am is around unity. In this study, we conservatively take the critical Elsasser number to be 0.3 and write where the prefactor f depth accounts for the fact that an actual current sheet can occur above the Am = 0.3 layer (see Figure  4), and Σ 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 depth 1. The typical ranges of f depth , which is affected by Hall effect, are ∼ 0.3-1 and ∼ 0.01-0.1 when the vertical magnetic field is parallel and antiparallel to the disk rotation axis, respectively. As in the case of f heat , we take the default value of f depth to unity to provide a maximum estimate for T acc, MHD , but also quantify the effect of varying the value of f heat in Section 3.2.3.
To determine the height of the Am = 0.3 layer, we compute the vertical profile of η 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 µm and the internal density to 1.4 g cm −2 , and vary the total grain surface area by changing the dust-to-gas mass ratio f 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 Am=0.3 of the Am = 0.3 layer is related to the column depth of the layer as

Numerical Procedures and Parameter Choices
To summarize, Equations (4) and (6) determine the surface density Σ(r) and midplane temperature T (r) of our quasisteady wind-driven accretion disk model. Because the righthand 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 Σ heat of the heating layers given by Equation (11). When Am at the midplane is higher than 0.3, we take Σ Am=0.3 to be Σ/2 by assuming that the heating occurs at the midplane.
Our model involves five parameters {κ R , f dg , α zφ , f heat , f depth }. The default values are taken to be κ R = 5 cm 2 g −1 , f dg = 0.01, α zφ = 10 −3 , f heat = 1, and f depth = 1. The default value of κ R is consistent with the opacity for 0.1 µm sized grains with the ISM abundance computed by Pollack et al. (1985).
In this study, we treat κ R and f 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.

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 ν = αc s H and take the parameter α to be constant throughout the disk. Assuming steady accretion, the gas surface density of the viscous disk is Σ =Ṁ/(3παc s H) (Lynden-Bell & Pringle 1974). The midplane temperature is given by where T acc, visc = 9ṀΩ 2 32πσ 128π 2 σαc 2 s 1/4 (14) represents the contribution from viscous accretion heating with τ mid = κ R Σ/2 being the vertical optical depth to the midplane (Hubeny 1990;Nakamoto & Nakagawa 1994;Oka et al. 2011). The effective optical depth τ heat for the vertically uniform viscosity is τ heat = τ mid /2 (Equation (A8)). The final expression of Equation (14) assumes τ heat 1, which holds true in our calculations (see Figure 5). Note that Equation (14) should be taken as an implicit equation for T acc, visc because the right-hand side depends on T acc, visc through c 2 s ∝ T acc, visc . For a fair comparison, we take α so that the MHD and viscous disk models have the similar gas surface density. Since α ∼ (r/H)α zφ (see the discussion in Section 2.1), for a typical turbulent protoplanetary disk with H/r ∼ 0.1, one has α ∼ 10α zφ . Therefore, we set α = 10 −2 . Although this is just an order-of-magnitude estimate, we have confirmed that the ratio of Σ in the MHD model to that in the viscous one is in the range of ≈ 0.5-1.4 at 1 au for the fiducial parameter set.

Fiducial Case
We begin by the default case where κ R = 5 cm 2 g −1 , f dg = 0.01, α zφ = 10 −3 , f heat = 1, and f 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 10 au (for t = 0.6 Myr) and r 3 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 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Ṁ in T 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, τ heat , shown in the lower panels of Figure 5. In the viscous model, the effective optical depth τ heat at the midplane is much higher than unity, with τ heat ≈ 500 (for t = 0.6 Myr) to ≈ 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 τ heat of the heating layer is only 1-3, and therefore the blanket effect is inefficient. This value of τ heat comes from our ionization calculations showing that Σ heat ≈ 0.2-0.6 g cm −2 . Interestingly, the depth of the heating layer is insensitive to the disk age and radial location. This is because Σ heat partly reflects the attenuation depth of ionizing X-rays. It should be noted, however, that Σ 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.
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 ≈ 10 Myr. In contrast, in the MHD model, the snow line arrives 1 au at t ≈ 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 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.

Parameter Dependence
We here study the dependence of the snow line tracks in the MHD model on the opacity, ionization fraction, accretion  (14)), respectively. The horizontal solid line marks the sublimation temperature of water ice, T = 170 K. Lower panels: radial profiles of the effective optical depths τ heat at the midplane for the viscous model (τ heat = τ mid /2; Equation (A8); dashed line) and for the MHD model (Equation (11); solid line).
stress, f heat , and f 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.

Variation with the Opacity
Since τ heat = κ R Σ 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 κ R = 5 cm 2 g −1 . When we increase κ R to 50 cm 2 g −1 , τ heat is increased to ≈ 20, and consequently the arrival of the snow line at 1 au is delayed to t ∼ 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 opac-ities. The depth Σ heat of the heating layer is insensitive to τ 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 1 Myr when τ heat 10.

Variation with the Ionization Fraction
In our ionization model, the dust-to-gas mass ratio f 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 dg gives a higher ionization fraction, a larger Σ heat , and consequently a higher midplane temperature. To see this effect quantitatively, we show in Figure 8 the snow line tracks with varying f dg from the fiducial value of 0.01, with the opacity value fixed. We also show Table 1. Summary of used parameters and the time t SL,1au that the snow line reaches 1 au. in Figure 9 the radial profile of τ heat for the cases shown in Figure 8. One can see that decreasing f dg by a factor of 10 from the fiducial value leads to an increase of τ heat by approximately the same factor (corresponding to Σ heat ≈ 6 g cm −2 ), resulting in the arrival of the snow line at 1 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 dg 0.01, accretion heat- ing is subdominant, and hence variation with f dg little affects the snow line evolution. The optical depth τ heat depends not only on the opacity but also on the ionization fraction via Σ 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 Σ heat can be comparable, and hence τ 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. We have shown in Section 3.2.2 that a high κ R and/or a low f 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 α zφ , f heat , and f depth .

Variation with the MHD Parameters
The magnetic stress α zφ mainly depends on the net flux of the vertical magnetic field (e.g., Hawley et al. 1995;Bai  It can be readily expected that the midplane temperature in the MHD model must be insensitive to α zφ . This parameter determines the magnitude of the wind accretion stress (see Equation (3)) and controls the disk surface density Σ (see Equation (4)). However, it is the column depth Σ heat of the heating layer, not the total column density Σ, that determines T acc, MHD . The depth of the heating layer is determined by the ionization structure well above the midplane, and is therefore insensitive to Σ. We demonstrate this in Figure 10, where we plot the snow line tracks for f dg = 10 −3 (see also Figure 8) but with different values of α zφ . It is clear that changing α zφ has essentially no effect on the snow line location.
The calculations shown above adopt f heat = f depth = 1. Such calculations give the maximum estimate of T 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 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 heat depends on the polarity of the net vertical magnetic field as well as f depth , and the range is approximately 10 −3 to 1. In Figure 11, we plot the snow line track for f dg = 10 −3 , but now assuming f 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 heat , accretion heating is no longer significant even with irr. only f dg = 10 3 f dg = 10 3 , f heat = 0.1 Figure 11. Same as Figure 7, but for f dg = 10 −3 with different values of f heat . The gray line is for the fiducial value of f heat = 1.
f dg = 10 −3 , and the snow line arrives at 1 au within 1 Myr as in the fiducial model ( Figure 6). The factor f depth is less than unity when Joule heating occurs at an altitude higher than where Am = 0.3. For instance, as shown in Figure 12 of Appendix C, f 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 acc, MHD ∝ ( f heat f depth ) 1/4 as long as τ heat 1, varying f depth from 1 to 0.1 gives the same effect as varying f heat from 1 to 0.1 (see Figure 11).

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 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 τ heat of the heating layer, which is the product of the opacity and the column depth Σ heat of the layer. The disk ionization state plays an important role here because it determines Σ heat . As demonstrated in Section 3.1, the value of κ R = 5 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 Σ 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 dustto-gas ratio used in the ionization calculations and therefore decreased whenever Σ 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.

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 ( 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 . 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 (∼ 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 ≈ 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.

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).

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 thaṅ M 10 −9 M 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.

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 irr are overestimated in particular in the early phase ( 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 ≈ 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 premain-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 acc 3 (see Section 2.2). The accretion still remains vigorous until ≈ 1 Myr (see Figure 3). However, in the case of high-entropy accretion, L acc 1 L after the star leaves the birthline and always satisfies L acc L int , where L int is the stellar intrinsic luminosity. Hence we can safely neglect L acc . In the case of low-entropy accretion, L acc can be higher than L 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.

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

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 α 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;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 α parameter showing the turbulence strength averaged over the disk height is ∼ 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 α parameter of ≈ 3 × 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.

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 selfgravitational 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 ∼10 au (see Kratter & Lodato 2016), which would be too far away for the wave to propagate to ∼ 1 au. Thus, the gravitational instability would not drive the disk heating when the snow line passes the rocky-planet-forming region.

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Ṁ = 10 −8 M yr −1 and α = 10 −3 , a Jupitermass 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 100 M ⊕ provides heat that increases the temperature to 200 K at ∼ 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).

Effects of the Dust Model on Accretion Heating
In this paper, we have assumed that the dust model is as 0.1 µ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.

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 τ heat . The dust growth reduces the abundance of smaller grains (e.g., Birnstiel et al. 2010Birnstiel et al. , 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.

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 τ heat ∝ f dg , as it is when f dg is changed from 0.01 to 0.001 (Section 3.2.2), τ heat is increased by a factor of ≈ 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.

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.
where L z is the energy flux carried away by the wind material and w wind is the zφ 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 wind is related to the wind-driven accretion rateṀ aṡ Eliminating w wind , we obtain Now we define f heat = Γ L z + Γ (B12) as the fraction of heat deposited inside the disk in the total liberated energy. In wind-driven accretion disks, using Equation (B11), we have The fraction f heat depends on the polarity of the net vertical magnetic field. The ranges of f heat in the MHD simulations of MBO19 are 0.02-0.5 and 5×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.

C. DEPTH OF THE HEATING LAYER
We here show that Σ Am=0.3 gives the maximum Σ 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 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} × 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/antiparallel 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. . The flux-weighted column density Σ heat (see Equation (A7)), against the column depth Σ Am=0.3 at the altitude where Am = 0.3, both of which are normalized by the column depth at the midplane, Σ mid = Σ/2. The data are given from the MHD simulation (MBO19) and additional simulations for f 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 depth = 1, 0.1, 0.01, and 0.001, respectively. Figure 12 shows Σ heat versus Σ Am=0.3 for all simulation data. Both axes are normalized by the midplane column depth Σ mid = Σ/2. Different data points correspond to simulations with different parameter sets. We also define f depth = Σ heat /Σ Am=0.3 (Equation (11)) and plot lines of f depth = 1, 0.1, 0.01 and 0.001. We find that all except one data point are below the line f 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 5 au, where Am 0.1 at the midplane. For this case, we can safely use Σ Am=0.3 as the upper limit of Σ heat .
Owing to the property of the Hall effect, the actual value of f 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 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 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 depth = 1 provides the upper limit of the midplane temperature determined by MHD accretion heating.