EVOLUTION OF SNOW LINE IN OPTICALLY THICK PROTOPLANETARY DISKS: EFFECTS OF WATER ICE OPACITY AND DUST GRAIN SIZE

A geometrically thin axisymmetric disk revolving around a central star is considered. A plane symmetry with respect to the midplane is assumed. Cylindrical coordinates (R, ϕ, Z) are adopted. The origin and the Z = 0 plane are placed at the central star and the disk midplane.

The hydrodynamical and thermal structures of the disk are modeled in almost the same way as in Dullemond et al. (2002) and Inoue et al. (2009) except for the disk's gas surface density Σ(R). The radiative cooling and two heating sources, irradiation by the central star and viscous dissipation, are considered. The temperature is determined with the radiative equilibrium condition. To obtain the whole disk structure, the 1+1D approach is adopted. The direct irradiation from the central star at the disk inner edge is ignored, and the disk structure is assumed to change gradually in the radial direction. This approach is justified if the snow line location is sufficiently far from the inner edge. The 2D structure of the inner disk region (the puffed-up inner rim and the adjacent shadowed region) seems to be important only in the innermost part of the disk around a T Tauri star, judging from results for Herbig Ae/Be stars of Dullemond (2002). Although the shadowed region extends to slightly large heliocentric distances, it does not matter for the scope of this study because it only lowers the disk temperature and does not shift the snow line location outward. Details are given in the Appendix.

The radial distribution of the gas surface density is modeled assuming that the disk is in the steady accretion state, i.e., the mass accretion rate $\dot{M}$ is constant along the radial direction R. The gas surface density in a region where R ≫ R_* (the stellar radius) is given as (e.g., Pringle 1981)

$$\begin{equation} \Sigma (R) = \frac{\dot{M}}{3\pi \nu _{t}(R)}, \end{equation} \tag{ 1 }$$

where ν_t(R) is the viscosity of the disk. The assumption of steady accretion is justified in the region around the snow line because the snow line is usually located within 10 AU, and the viscous diffusion timescale there is sufficiently smaller than the entire disk evolution timescale.

The viscosity ν_t is modeled with the α-prescription (Shakura & Sunyaev 1973) as ν_t = αc²_s/Ω_K, where α is a non-dimensional parameter, $c_{s}=\ssty\sqrt{kT / (\bar{\mu } m_{u})}$ is the isothermal sound velocity of gas, k is the Boltzmann constant, T is the gas temperature, $\bar{\mu }$ is the mean molecular weight, m_u is the atomic mass unit, $\Omega _{{K}} =\ssty\sqrt{GM_{*} / R^{3}}$ is the Kepler angular velocity, G is the gravitational constant, and M_* is the stellar mass. The value of α ranges from 0.001 to 0.1 according to ideal MHD simulations (Hawley et al. 1995, 1996). The value of α = 0.01 is adopted as a fiducial one throughout the disk in this study, except in Section 3.3. In the evaluation of ν_t, the sound velocity c_s at the midplane (Z = 0) is used. A value of $\bar{\mu } = 2.3$ is used for the mean molecular weight.

The frequency-dependent radiative transfer is solved numerically. Absorption and scattering by dust particles are taken into consideration as the opacity source in the disk. Isotropic scattering is assumed when the dust grain size is sufficiently smaller than the radiation wavelength. When the dust grain size is larger than the wavelength, the scattering coefficient is set to zero because forward scattering dominates in that case.

The steady-state disk structure is solved self-consistently as a function of the mass accretion rate. Although density, viscosity, and temperature are related to one another, they are numerically obtained consistently by iterative calculations.

Only dust particles suspended in the disk are considered to be the opacity source because the gas opacity is negligibly small. Two kinds of dust particles, silicate and pure water ice particles, are assumed to be present separately. Effects of icy mantles on silicate cores will be discussed in Section 5.2. It is assumed that dust particles are uniformly sized spheres and mass ratios of the silicate and ice particles to the gas are constant throughout the disk. Radial transport, vertical mixing, and sedimentation of dust particles are neglected. Also, complete thermal coupling between dust particles and gas is assumed; this is valid in the disk considered in this study (Kamp & Dullemond 2004).

Absorption and scattering coefficients of dust particles and mass ratios to the gas are taken from Miyake & Nakagawa (1993); silicate and (maximum) ice mass ratios to the gas are ζ_sil = 0.0043 and ζ_ice = 0.0094, respectively, which are consistent with the solar elemental abundance given by Anders & Grevesse (1989). With this water ice abundance and the mean molecular weight of 2.3, the H₂O mole fraction in the disk gas becomes 1.2 × 10⁻³. Figure 1 shows the absorption and scattering coefficients of a unit mass disk medium.

Zoom In Zoom Out Reset image size

The ice abundance in the disk gas is determined with the thermal equilibrium condition under which the partial pressure of water vapor is limited by the saturated vapor pressure (Bauer et al. 1997):

$$\begin{equation} P_{\mathrm{sat}}(T) = \exp (-6070\;\mathrm{K} / T + 30.86)\, \mathrm{dyn \, cm^{-2}}. \end{equation} \tag{ 2 }$$

The partial pressure of water vapor $P_{\mathrm{H_{2}O}}$ is given approximately by the product of the water vapor mole fraction $X_{\mathrm{H_{2}O}}$ and the total gas pressure:

$$\begin{equation} P_{\mathrm{H_{2}O}} = X_{\mathrm{H_{2}O}} P. \end{equation} \tag{ 3 }$$

The water vapor mole fraction $X_{\mathrm{H_{2}O}}$ is limited to 1.2 × 10⁻³, so excess water molecules condense as water ice. Thus, $X_{\mathrm{H_{2}O}}$ is given by

$$\begin{equation} X_{\mathrm{H_{2}O}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1.2 \times 10^{-3} & (1.2 \times 10^{-3} P \le P_{\mathrm{sat}}(T)), \\ P_{\mathrm{sat}}(T)/ P & (1.2 \times 10^{-3} P > P_{\mathrm{sat}}(T)). \end{array}\right. \end{equation} \tag{ 4 }$$

Then, the mass ratio of ice to the total water molecules in the disk gas, x_ice, is given by

$$\begin{equation} x_{\mathrm{ice}} = 1 - X_{\mathrm{H_{2}O}} {/} (1.2 \times 10^{-3}). \end{equation} \tag{ 5 }$$

Using x_ice, the absorption and scattering coefficients of the disk medium, κ_abs and κ_sca, are given as

$$\begin{eqnarray} \kappa _{\mathrm{abs}} &=& \kappa _{\mathrm{sil, abs}} + x_{\mathrm{ice}} \kappa _{\mathrm{ice, abs}}, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \kappa _{\mathrm{sca}} &=& \kappa _{\mathrm{sil, sca}} + x_{\mathrm{ice}} \kappa _{\mathrm{ice, sca}}, \end{eqnarray} \tag{ 7 }$$

where κ_{sil, abs}, κ_{sil, sca}, κ_{ice, abs}, and κ_{ice, sca} represent the absorption and scattering coefficients of the disk medium attributed to silicate and ice dust particles, respectively. Equations (6) and (7) are held for each frequency of the radiative transfer. The disk temperature depends on x_ice; hence, the temperature and x_ice are solved consistently by iterative calculations.

Numerical results shown in this subsection are obtained by calculations with the following input parameters: T_eff = 3000 K, R_* = 2 R_☉, and M_* = 0.5 M_☉. The dust grain size is 0.1 μm.

3.1.1. Disk Structure

Figures 2 and 3 show vertical temperature profiles at various radii with $\dot{M} = 10^{-8}\,M_{\odot }\;\mathrm{yr}^{-1}$ and 10⁻¹⁰ M_☉ yr⁻¹, respectively, as typical cases for high and low mass accretion rates. In Figure 2, the temperature in the disk with ice opacity is higher than that in the one without it. The dominant heating source in this case is viscous dissipation, and the blanket effect of dust particles is enhanced due to the ice in the upper layers. On the other hand, in Figure 3, the midplane temperatures are not considerably different from each other, though ice condenses. In this case, the dominant heating source is the incident stellar radiation, so the midplane temperature does not depend so much on dust opacity. The temperature in the disk is lowered due to the ice condensation. The low temperature causes a slightly low altitude of the scale height and H_s as seen in Figure 3 (H_s is the height from the midplane where the incident stellar radiation energy flux is reduced by a factor of e).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The disk temperature is altered only a little by considering dust scattering, because a single scattering albedo dust particles of size 0.1 μm is small. Though small, the midplane temperature is increased slightly when viscous heating is dominant (Figure 2); the scattering by dust particles enhances the optical depth in the disk and increases the blanket effect. On the other hand, the midplane temperature is lowered slightly when stellar irradiation is the dominant heating source (Figure 3), since a large fraction of incident radiation is reflected toward space (Dullemond & Natta 2003; Inoue et al. 2009).

3.1.2. Evolution of the Ice-Condensing Region

Figure 4 shows the ice-condensing region in the (R, Z) sectional plane with various mass accretion rates.¹ A decrease in the mass accretion rate can be regarded as the time evolution of the disk.

Zoom In Zoom Out Reset image size

When $\dot{M}$ is high (⩾10⁻⁹ M_☉ yr⁻¹), the condensation front has a two-branched structure, as found by Cassen (1994) and Davis (2005). The condensation front in the upper layer extends closer to the central star than to the snow line (hereafter, the snow line is defined as the condensation front at the midplane) irrespective of the presence of ice. The lower branch of the condensation front is formed by a high temperature caused by viscous dissipation, while the upper branch is formed by a high temperature caused by direct stellar radiation and a low partial pressure of water vapor. The snow line shifts outward by a factor of 1.3 due to the ice opacity independent of the mass accretion rate.

On the other hand, when $\dot{M}$ is low (⩽10⁻¹⁰ M_☉ yr⁻¹), the condensation front is a smooth curve and the snow line location is not shifted significantly by the ice opacity. This can be attributed to two reasons: when the disk is mainly heated by stellar irradiation (1) the midplane temperature does not depend so much on the opacity of the disk medium and (2) icy dust particles do not condense in the upper layer above the snow line due to the temperature increasing with height.

3.1.3. Evolution of the Snow Line

Figure 5 shows the snow line locations obtained with and without ice opacity as a function of the mass accretion rate. In both cases, the snow line migrates inward when $\dot{M}$ is high (≳ 10⁻¹⁰ M_☉ yr⁻¹) and outward when $\dot{M}$ is low (≲ 10⁻¹⁰ M_☉ yr⁻¹). An inward migration is caused by the decrease of the midplane temperature due to the reduction of viscous heating and the decrease of the blanket effect. The inward migration ceases when stellar irradiation starts to dominate the heating around the midplane. Then, an outward migration is started by the decrease in the disk's optical thickness, which results in a higher temperature due to a deeper penetration of the stellar radiation into the disk interior. The snow line in this study migrates outward more slowly than the results of Garaud & Lin (2007) (see Section 3.3).

Zoom In Zoom Out Reset image size

When $\dot{M}$ is high (≳ 10⁻¹⁰ M_☉ yr⁻¹), the ratio of the snow line distance with ice opacity (R_{SL, sil + ice}) to that without ice opacity (R_{SL, sil}), f_SL = R_{SL, sil + ice}/R_{SL, sil}, is about 1.3 and almost constant. This constant ratio of the snow line shift factor f_SL originates from the assumptions that the disk is the steady accretion disk and the mass fraction of ice dust particles to silicate dust particles is uniform throughout the disk (see Section 4.1). On the other hand, when $\dot{M}$ is low (≲ 10⁻¹⁰ M_☉ yr⁻¹), the snow line location is not substantially affected by ice opacity (i.e., f_SL ≃ 1). This small difference in snow line location is caused by a small difference in the disk structure around the snow line.

The stellar parameters in this subsection are set as L_* = 1 L_☉, M_* = 1 M_☉, and T_eff = 4000 K for comparison with the solar system. Figure 6 shows the heliocentric distance of the snow line as a function of the mass accretion rate. The different curves correspond to the different dust grain sizes assumed in each calculation. Calculations are stopped when the optical depth for stellar radiation in the region around the snow line becomes less than unity. When the grain size is ⩾100 μm, the scattering coefficient is set to zero.

Zoom In Zoom Out Reset image size

During the inward migration of the snow line, the snow line location is almost independent of the dust grain size if the grain size is smaller than 10 μm. The wavelength of radiation from the disk interior at the condensation temperature (T ≃ 170 K) is longer than 10 μm, so the opacity of the dust particle is in the Rayleigh regime, dependent only on the total mass of the dust particles, and dominated by absorption. On the other hand, when the size is larger than 10 μm, i.e., in the geometrical optics regime, the opacity depends on the total cross section of the dust particles and decreases with the dust grain size; hence, the temperature enhancement by viscous heating becomes inefficient, although both absorption and scattering contribute in this case.

When the dust grain size is larger than 10 μm, the evolutional track of the snow line levels off at the minimum value, because the opacity of large grains is almost independent of the wavelength emitted from the disk. This implies that the snow line necessarily passes the heliocentric distance of 1 AU, unless the dust grain size exceeds a certain large value. This point will be considered in Section 4.4.

The shift ratio of the snow line f_SL is qualitatively the same as in Figure 5. The snow line shifts to a larger heliocentric distance in its inward migration phase, whereas it does not shift substantially in its outward migration and leveling-off phases.

The shift ratio of the snow line f_SL depends on the dust grain size: f_SL ≃ 1.3 for dust grain sizes ≲10 μm, and f_SL ≃ 1.6 for grain sizes ≳100 μm. When the dust grain size is small, the total dust opacity increases in proportion to the total mass of the dust, while the total dust opacity increases in proportion to the total cross section of dust particles when the dust grain size is large. This different dependence of the opacity on the total amount of dust causes the different dependence of the shift ratio.

The critical mass accretion rate, that is the mass accretion rate with which the snow line starts to migrate outward, also depends on dust grain size. The snow line starts to migrate outward when the disk region around the snow line becomes transparent. When the dust grain size is sufficiently larger than the wavelength of the stellar radiation (∼1 μm), the total opacity of the dust is inversely proportional to the grain size. Thus, the critical mass accretion rate is proportional to the dust grain size. Indeed, when the dust grain size is 10 μm, 100 μm, and 1 mm, the snow line starts to migrate outward at about 10⁻¹¹ M_☉ yr⁻¹, 10⁻¹⁰ M_☉ yr⁻¹, and 10⁻⁹ M_☉ yr⁻¹, respectively (Figure 6).

The snow line evolution obtained from our numerical simulations is compared with that obtained from the semi-analytical model of Garaud & Lin (2007). The stellar parameters are L_* = 1 L_☉, M_* = 1 M_☉, and R_* = 1 R_☉, and the viscosity parameter is α = 0.001. Other settings are almost the same as in Garaud & Lin (2007).

Garaud & Lin (2007) used the Planck mean opacity; thus, a new frequency-dependent opacity model of the mixture of disk gas and dust particles is developed so that its Planck mean agrees with that of Garaud & Lin (2007). The Planck mean opacity is proportional to temperature, so the absorption opacity is modeled as a linear function of frequency: κ^abs_ν = Cν, where C is a constant determined by

$$\begin{equation} \frac{ \int _{0}^{\infty } \!\kappa _{\nu }^{\mathrm{abs}} B_{\nu }(T_{\mathrm{eff}})\,d\nu }{ \int _{0}^{\infty } \! B_{\nu }(T_{\mathrm{eff}})\,d\nu } = \kappa _{V}, \end{equation} \tag{ 8 }$$

where κ_V is the Planck mean opacity with the stellar effective temperature, and κ_V is set to be 2 cm² g⁻¹ according to Garaud & Lin (2007).

The snow line location calculated as a function of $\dot{M}$ is shown in Figure 7, as is the result of Garaud & Lin (2007). Our result agrees well with that of Garaud & Lin (2007) for higher $\dot{M}$ (>10⁻⁹ M_☉ yr⁻¹), though our result does not show an abrupt outward shift of the snow line at $\dot{M} \sim 10^{-9.5}\,M_{\odot }\;\mathrm{yr}^{-1}$. Garaud & Lin (2007) adopted a semi-analytical method to solve the temperature in the optically thin region. Although semi-analytical methods are much easier to use than numerical models, they may not be accurate enough for the evaluation of the disk midplane temperature in optically thin regions. (They are accurate enough and useful in optically thick regions.) The discontinuous evolution found by Garaud & Lin (2007) may come from the simplified treatment of the radiative transfer. Hereafter, the discontinuity is ignored and the overall features of the snow line evolution are discussed.

Zoom In Zoom Out Reset image size

In order to obtain the midplane temperature, the estimation of the height of the superheated layer H_s is a key factor, because its radial dependence determines the grazing angle of the disk for irradiation. When the midplane temperature in a region is increased, the surface density there is lowered according to the steady accretion assumption (Equation (1)). Then, the grazing angle decreases, the midplane temperature tends to fall, and the surface density tends to increase: a negative feedback takes place. Therefore, in our calculations, the disk does not become optically so thin and the snow line migrates outward more slowly than in the result of Garaud & Lin (2007), in which H_s is assumed to be proportional to the disk pressure scale height.

The minimum heliocentric distance of the snow line obtained from our calculations is slightly smaller than that of Garaud & Lin (2007). This may be attributed to the two-layer model (Chiang & Goldreich 1997) they adopted. Inoue et al. (2009) showed that the three-layer approach is necessary to obtain an accurate midplane temperature.

Figures 8 and 9 show the midplane temperature and the surface density profiles from our calculation and from Garaud & Lin (2007). These figures clearly show that the midplane temperatures of both models agree well in the inner optically thick region, whereas our result is lower than that of Garaud & Lin (2007) in the outer optically thin region. Similarly, the surface density distributions in both models agree in the inner region, whereas there is a difference in the outer region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We conclude that if we are concerned with the detailed evolution of the snow line in the optically thin phases in which planet formation may proceed, we need to calculate the temperature numerically. (If we are concerned with optically thick phases or qualitative features of the snow line evolution in entire phases, semi-analytical calculations are sufficient.)

In the early phase in which the snow line migrates inward, the main heating source at the midplane around the snow line is viscous dissipation (Figure 2). Thus, the temperature is determined by viscous heating and radiative cooling. The diffusive radiative energy transfer gives the midplane temperature, T_c, as (e.g., Nakamoto & Nakagawa 1994)

$$\begin{equation} T_{c} = \left[ \left(\frac{3}{2}+\frac{9\tau _{c}}{16} \right) \left(\frac{GM_{*}\dot{M}}{4\pi \sigma R^{3}} \right) \right]^{\frac{1}{4}}, \end{equation} \tag{ 9 }$$

where τ_c is the optical depth for the midplane and σ is the Stefan–Boltzmann constant. Substituting τ_c = κΣ/2 (κ is the Rosseland mean opacity of the disk medium) and noticing that τ_c is much larger than unity, we obtain

$$\begin{eqnarray} &&T_{c} \simeq \left[ \left(\frac{9\kappa \Sigma }{32} \right) \left(\frac{GM_{*}\dot{M}}{4\pi \sigma R^{3}} \right) \right]^{\frac{1}{4}} \propto \left(\frac{\kappa \Sigma }{R^{3}} \right)^{\frac{1}{4}}. \end{eqnarray} \tag{ 10 }$$

The surface density of the steady accretion disk is given by $\Sigma =\dot{M}/(3\pi \alpha c_{s}^{2}/\Omega)\propto 1/(T_{c}R^{3/2})$. So, we have

$$\begin{equation} T_{c}^{5} \propto \kappa R^{-\frac{9}{2}}. \end{equation} \tag{ 11 }$$

In the annulus located at the snow line, both silicate and ice particles contribute to the opacity κ (see Figure 4). This Rosseland mean opacity is denoted as κ_{sil + ice}, whereas the Rosseland mean opacity due to silicate particles only is expressed as κ_sil. Since the ice condensation temperature at the midplane is almost independent of the heliocentric distance, the right-hand side of Equation (11) can be regarded as constant, and a relation between the snow line locations with and without water ice opacity, R_{SL, sil + ice} and R_{SL, sil}, is obtained as

$$\begin{equation} \kappa _{\mathrm{sil+ice}} R_{\mathrm{SL, sil+ice}}^{-\frac{9}{2}} = \kappa _{\mathrm{sil}} R_{\mathrm{SL, sil}}^{-\frac{9}{2}}. \end{equation} \tag{ 12 }$$

The shift factor of the snow line location due to the additional ice is then written as

$$\begin{equation} f_{\rm SL} = \frac{R_{\mathrm{SL, sil+ice}}}{R_{\mathrm{SL, sil}}} = \left(\frac{\kappa _{\mathrm{sil+ice}}}{\kappa _{\mathrm{sil}}} \right)^{\frac{2}{9}}. \end{equation} \tag{ 13 }$$

Adopting the dust opacity model of Miyake & Nakagawa (1993), we have κ_{sil + ice}/κ_sil ≃ 3 for the thermal radiation at 170 K (the midplane temperature at the snow line) when the dust grain size is 0.1 μm. Then, f_SL is evaluated as $f_{\rm SL} = 3^{\frac{2}{9}} = 1.28$, which is consistent with our numerical results. Equation (13) implies that the snow line location is shifted outward in accordance with the water abundance around the snow line.

The snow line location is considerably affected by α, because α relates to the disk's optical thickness and energy generation rate by viscous dissipation. In the early phase of the disk evolution, the midplane temperature in the inner disk region is obtained in the same way as in Section 4.1. Then, we have T⁵_c∝α⁻¹R^−9/2. We then find the dependence of the snow line location R_snow on α as

$$\begin{equation} R_{\mathrm{snow}} \propto \alpha ^{-\frac{2}{9}}. \end{equation} \tag{ 14 }$$

Taking into account the range of α (0.001 ≲ α ≲ 0.1), the snow line location shifts inward or outward by at most 67% from the results with α = 0.01 in the inward migration phase. Note that the minimum heliocentric distance of the snow line hardly changes, because it is determined not by viscous heating but by stellar radiation heating.

In the later phase of the disk evolution, a change in α changes the correlation between R_snow and $\dot{M}$. In the outward migration phase, R_snow is a function of Σ, and $\dot{M}$ is proportional to the product of ν_t and Σ. The viscosity ν_t is proportional to α, so $\dot{M}$ corresponding to each R_snow varies in proportion to α. Thus, R_snow in the later phase is shifted in proportion to α (Figure 6). This means that the $\dot{M}$ at which the disk becomes optically thin to the stellar radiation also decreases in proportion to α.

The dust-to-gas mass ratio, ζ_d, changes the opacity of the disk medium. Hence, it changes the disk temperature in the same way as ice opacity. In the inward migration phase, assuming that the opacity of the disk medium varies in proportion to ζ_d, a proportional relation between R_snow and ζ_d is obtained in the same way as in Sections 4.1 and 4.2:

$$\begin{equation} R_{\mathrm{snow}} \propto \zeta _{d}^{\frac{2}{9}}. \end{equation} \tag{ 15 }$$

In the outward migration phase, the correlation between R_snow and $\dot{M}$ depends on ζ_d. In this phase, R_snow is a function of solid mass surface density (ζ_dΣ), and $\dot{M}$ is proportional to Σ, which is the product of solid mass surface density and the inverse of ζ_d. Thus, $\dot{M}$ corresponding to each R_snow varies in inverse proportion to ζ_d. Hence, R_snow is shifted in inverse proportion to ζ_d (Figure 6). The mass accretion rate with which the disk becomes optically thin to the stellar radiation increases in inverse proportion to ζ_d.

According to the results shown in Section 3, the snow line comes inside the Earth's current orbit in the course of the disk evolution. However, if planetesimals in the terrestrial planet region were icy, this is not consistent with current solar system bodies. Here, discussing the snow line evolution track on the $(\dot{M}, R_{\mathrm{snow}})$-plane in Figure 6, we examine to what extent that conclusion is certain.

From Figure 6, we can see that if (α/ζ_d) is increased by more than a factor of 100 and if the dust grain size is 1 mm, the snow line would not come inside Earth's orbit. However, there are some caveats. Noticing that the range of α is from 0.001 to 0.1 and the fiducial value used in our calculations is α = 0.01, we realize that ζ_d should be decreased by a factor of at least 10 from the value found by Miyake & Nakagawa (1993) to increase the value of (α/ζ_d) by more than a factor of 100. When α  <  0.1, ζ_d should be much smaller depending on the size of α. On the other hand, as the dust grain size becomes small, (α/ζ_d) should be large (or ζ_sil should be small). Thus, it seems that the snow line necessarily comes inside Earth's orbit if the dust grain size is smaller than 1 mm, ζ_d is higher than a tenth of that of the solar abundance, and α ranges from 0.001 to 0.1.

Another possibility is an increase in ice abundance; this can prevent the snow line from coming inside Earth's orbit. An increase in ice abundance lowers the value of $\dot{M}$ at which the snow line migration levels off at the minimum distance or starts to move outward. For example, if the ratio of water ice mass to disk gas, ζ_ice, is increased by a factor of 10⁴, it is expected that the snow line would not come inside Earth's orbit even in the case in which α = 0.01, the dust grain size is 1 mm, and ζ_sil = 0.0043 (the solar abundance). When the dominant dust grain size or (α/ζ_d) is small, a greater increase in ice abundance is needed.

In our solar system, the water distribution shows a drastic change at the asteroid belt; this may be a clue to the snow line. Terrestrial planets are thought to form from planetesimals that have formed inside the snow line, otherwise a large amount of water is inevitably accumulated into planets. Then, the timing of planetesimal formation may be restricted by the snow line evolution.

First, we consider the possibility that planetesimals form after the snow line migrates outward from the terrestrial planet region. When we look at the solid mass, planetesimal formation with dust grains of large size is favored because the gas surface density at which the snow line changes its direction of motion and moves outward increases as the dust grain size increases (Figure 6). The most favorable case from our numerical results is for dust grains of size 1 mm. Our numerical simulation shows that the surface density of solid material at the heliocentric distance of 1 AU is 1.1 × 10⁻² g cm⁻² when the outwardly moving snow line reaches 1 AU. This is much smaller than that of the minimum mass solar nebula model (about 10 g cm⁻²). To match the minimum mass solar nebula model, or to match the current terrestrial planet's mass, the size of the dust particles would necessarily be larger than 1 m. This would be unrealistic, partly because complete depletion of fine dust particles smaller than 1 m would be difficult, and partly because dust particles of size 1 m are eliminated efficiently by accretion toward the central star due to gas drag. Thus, planetesimal formation during the snow line's outward migration would be unrealistic.

Next, let us consider the possibility of planetesimal formation during the snow line's inward migration. If planetesimal formation completes before the inwardly migrating snow line reaches the terrestrial planet region, planetesimals devoid of water can be formed. As for the solid mass, planetesimal formation in this phase is favorable. However, as the grain size grows, the snow line quickly migrates inward and eventually comes inside Earth's orbit; in Figure 6, at a fixed $\dot{M}$, we can see that R_snow is small if the grain size is large. Hence, the dust grain growth and planetesimal formation should be completed in a short period of time compared to the timescale with which the snow line moves inward. If not, a large amount of ice would necessarily accumulate deep within planetesimals. We are not sure if planetesimal formation is completed in such a short period of time.

In conclusion, the snow line location during disk evolution does not seem to match the current solar system. This inconsistency may originate from the following assumptions adopted in this study: (1) the continuous boundary condition at the inner disk edge, (2) uniform dust grain size and uniform dust-to-gas mass ratio throughout the disk, and (3) neglecting the evolution of solid material and water distribution.

Here, we consider possible scenarios by relaxing the above assumptions.

• 1.  
  Continuous disk boundary. The presence of the inner disk edge, where the disk is heated by direct stellar radiation and becomes hot, is ignored in this assumption. If the dust grain grows at the inner edge and the inner edge moves outward after planetesimal formation, dust grain growth would always proceed in a hot environment. Then the resulting planetesimals would be devoid of water.
• 2.  
  Uniform dust grain size and dust-to-gas mass ratio. In a real disk, the dust grain size and the dust-to-gas mass ratio would be non-uniform. Dust particles grow mainly by collisions. Growth tends to proceed from the inner disk region, so the dust grain size in the inner disk is likely to be larger than that in the outer disk. When the dust grain is large enough, the height of the disk surface is lowered. Consequently, a disk region that is located just outside the dust grain region developed when dust grain growth is about to start can receive more stellar radiation flux. Then, dust grain growth may proceed at a hot temperature and water-deficit planetesimals may form.
• 3.  
  Solid material and water distribution evolution. According to Ciesla & Cuzzi (2006) and Garaud (2007), large amounts of water vapor and fine dust particles are transported from the outer to the inner part of the disk during disk evolution. If fine dust particles or water vapor are supplied to the terrestrial planet region, and the snow line location is kept farther from the terrestrial planet region until planetesimal formation is completed, dust grain growth may proceed in a water-ice-free environment; consequently, water-devoid planetesimals may be formed. Whether or not this mechanism works may depend on the initial condition of the disk.

Yet another possible scenario may be the elimination of water from icy planetesimals after the snow line migrates outward from the terrestrial planet region. Future work on these points is needed.

In this study, it is assumed that ice condenses as pure ice particles and that no ice mantle forms on silicate particles. It would be more realistic, however, if H₂O molecules condense as ice mantles on silicate particles. Here, we discuss how ice mantle formation affects the snow line location. We consider its effects only in the inward migration phase, since ice opacity does not change the snow line location in the outward migration phase.

To evaluate the effects of the ice mantle, we consider to what extent the dust opacity changes due to the ice mantle. Dust opacity properties change as the ratio between the dust grain size and the wavelength of the thermal radiation (λ ∼ 10–20 μm) changes. When the dust grain size is sufficiently smaller than 10–20 μm (the Rayleigh regime), the total dust opacity is hardly affected by the ice mantle. On the other hand, when the dust grain size is sufficiently larger than 10–20 μm (the geometrical optics regime), the ice mantle changes the total dust opacity.

When the ice condenses, the opacity of dust particles generally increases. However, the opacity enhancement factor g_ice can be different between cases in which the ice condenses as independent ice particles and those in which the ice condenses as ice mantles on silicate cores, even though the amount of condensed ice is the same. In the geometrical optics regime, g_ice for the mantle case becomes smaller than that for the independent case, since the cross section of ice mantles is smaller than that of the independent particles. According to Miyake & Nakagawa (1993), the material densities of silicate and ice are 3.3 and 0.92 g cm⁻³, respectively, and the mass fractions of silicate and water ice in the disk medium are 0.0043 and 0.0094, respectively. This means that ice has a volume about six times larger than that of silicate per unit mass at disk medium. As a result, if only pure ice particles are formed, the total cross section of the dust particles is increased by a factor of 7.9, while if ice mantles are formed on silicate cores, the total cross section of dust particles is increased only by a factor of 3.7. Thus, ice mantle formation reduces the opacity enhancement factor by ice condensation. Consequently, the snow line's shift ratio by ice condensation (f_SL) would be reduced from 1.54 to 1.33 due to ice mantle formation (see Equation (13)).

Finally, we consider the non-uniformity of the dust grain size. When independent ice particles are formed, the opacity enhancement depends on the size distribution of ice particles. If all the ice particles are sufficiently small compared to the wavelength of the thermal radiation (the Rayleigh regime), the opacity enhancement depends only on the total mass of ice and is largest. When the dominant size of ice particles is larger than 10–20 μm, and as the dominant size increases, the opacity enhancement decreases. On the other hand, if ice mantles are formed on silicate particles, the opacity enhancement would be limited to a certain range. The real situation should be somewhere between the Rayleigh regime and the geometrical optics regime, so the enhancement of the dust opacity by ice condensation would be around 30%.

In summary, it is likely that the enhancement of the dust opacity would be limited to a range around 30% when ice mantles are formed around silicate dust particles. If pure ice particles are formed, the size distribution of ice particles is needed to evaluate the enhancement accurately.

The hydrostatic equilibrium along the vertical direction is expressed as

$$\begin{equation} \frac{\partial P(R,Z)}{\partial Z} = -\rho (R,Z)\frac{GM_{*}}{R^{3}}Z, \end{equation} \tag{ A1 }$$

where P is the gas pressure and ρ is the mass density of the gas. The equation of the state of the gas is given as $P = {\rho kT \over \bar{\mu } m_{u}}= c_{s}^{2} \rho$, where $\bar{\mu }$ is the mean molecular weight, m_u is the atomic mass unit, and c_s is the sound velocity. Integrating the density over Z, the surface density Σ is given as Σ = ∫^∞_{− ∞}ρ dZ.

Two heating sources are taken into consideration: viscous dissipation and irradiation of the central star. Balancing the heating sources and the radiative cooling, the following equilibrium temperature is obtained:

$$\begin{eqnarray} &&q_{\mathrm{vis}}\, {+}\, q_{\mathrm{irr}}\, {+}\, \int _{0}^{\infty }\!\rho \kappa _{\nu }(1-\varpi _{\nu })\, d\nu \oint I_{\nu }(\bm {\Omega })\,d\Omega\nonumber\\ &&\quad = 4\pi \int _{0}^{\infty }\!\rho \kappa _{\nu }(1-\varpi _{\nu })B_{\nu }(T)\,d\nu, \end{eqnarray} \tag{ A2 }$$

where q_vis is the heating rate by viscous dissipation, q_irr is the heating rate by irradiation of the central star, κ_ν is the extinction coefficient, ϖ_ν is the single scattering albedo, $I_{\nu }(\bm {\Omega })$ is the diffuse radiation field in the disk (intensity without including the direct radiation from the central star), and B_ν(T) is the Planck function. Quantities with subscript ν are a function of frequency ν. The integral ∮dΩ on the left-hand side of Equation (A2) represents the integral over the entire solid angle. The viscous heating rate q_vis is given by (e.g., D'Alessio et al. 1998) $q_{\mathrm{vis}} = \frac{9}{4}\rho \nu _{t}\Omega _{{K}}^{2}.$

The diffuse component of the intensity $I_{\mu,\nu }(\bm {\Omega })$ and q_irr are obtained by solving the radiative transfer. Assuming a plane-parallel structure along the Z-direction, the radiative transfer equation is written as

$$\begin{equation} \mu \frac{dI_{\mu,\nu }(Z)}{dZ} = \rho \kappa _{\nu } (S_{\nu }(Z) - I_{\mu,\nu }(Z)), \end{equation} \tag{ A3 }$$

where μ is the cosine of the propagation direction of radiation and S_ν is the source function. The source function S_ν is given by

$$\begin{eqnarray} S_{\nu }(Z) &=& (1-\varpi _{\nu }(Z)) B_{\nu }(T) + \varpi _{\nu }(Z) \nonumber\\ &&\times \frac{1}{2}\int _{-1}^{1} I_{\mu,\nu }d\mu + \varpi _{\nu }(Z)\frac{F_{\mathrm{irr},\nu }(Z)}{4\pi }, \end{eqnarray} \tag{ A4 }$$

where F_{irr, ν}(Z) is the radiation energy flux from the central star. Therein, isotropic scattering is assumed. The energy flux F_{irr, ν}(Z) is evaluated using the so-called grazing angle recipe (e.g., Chiang & Goldreich 1997; Dullemond et al. 2002):

$$\begin{equation} F_{\mathrm{irr},\nu }(Z) = \frac{L_{\nu }}{4\pi R^{2}} \exp \left(-\frac{\tau _{\nu }(R,Z)}{\beta } \right), \end{equation} \tag{ A5 }$$

where β, L_ν, and τ_ν(R, Z) are the cosine of the penetrating angle of the radiation from the central star into the disk surface, the luminosity of the central star, and the optical depth between a point (R, Z) and infinity along the vertical direction (R, ∞), respectively. The optical depth τ_ν is defined as τ_ν(R, Z) = ∫^∞_Zκ_νρ(R, Z)dZ. The central star is assumed to emit blackbody radiation with effective temperature T_eff, so the luminosity is given by L_ν = 4π²R²_*B_ν(T_eff). β is evaluated by (e.g., Kusaka et al. 1970; Chiang & Goldreich 1997; Dullemond et al. 2002)

$$\begin{equation} \beta (R) = 0.4\frac{R_{*}}{R}+ \xi (R)\frac{H_{s}}{R}, \end{equation} \tag{ A6 }$$

where H_s is the height from the midplane where the radiation from the central star loses its energy by 1 − (1/e) and ξ is the so-called flaring index defined as $\xi = \frac{d \ln (H_{s}/R)}{d \ln R}.$ Then, the heating rate by irradiation of the central star in Equation (A2) is given by

$$\begin{equation} q_{\mathrm{irr}} = \int _{0}^{\infty }\rho (Z)\kappa _{\nu }(1-\varpi _{\nu })F_{\mathrm{irr},\nu }(Z)\,d\nu. \end{equation} \tag{ A7 }$$

The boundary condition for Equation (A3) is

$$\begin{eqnarray} &&I_{\mu, \nu }(+\infty) = 0 \quad (\mu <0), \end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray} &&I_{\mu, \nu }(0) = I_{-\mu, \nu }(0) \quad (\mu \ge 0). \end{eqnarray} \tag{ A9 }$$

This condition represents the mirror boundary at the Z = 0 plane and no incident diffuse radiation at the disk surface.

Vertical profiles of the diffuse radiation $I_{\mu,\nu }(\bm {\Omega })$ and the temperature T are obtained by solving Equations (A2)–(A4) iteratively. A straightforward way to solve these equations converges very slowly, so the variable Eddington factor (VEF) method, described in Dullemond et al. (2002) and Inoue et al. (2009), is employed. In this study, the VEF method is modified to deal with the heating by viscous dissipation.