Gravitoviscous Protoplanetary Disks with a Dust Component. V. The Dynamic Model for Freeze-out and Sublimation of Volatiles

In this section, we describe the dust growth scheme, which is an updated version of the scheme first presented in Vorobyov et al. (2018). In the FEOSAD code, small and grown dust components are represented by their surface densities Σ_d,sm and Σ_d,gr, respectively. Each dust population has the size distribution N(a) described by a simple power-law function N(a) = Ca^−p with a fixed exponent p = 3.5 and a normalization constant C. For small dust, the minimum size is ${a}_{\min }=5\times {10}^{-7}$ cm and the maximum size is a_* = 10⁻⁴ cm. For grown dust, a_* is the minimum size and ${a}_{\max }$ is the maximum size, which can vary due to dust coagulation and fragmentation.

The dynamics of small and grown dust grains is described by the following continuity and momentum equations (note that small dust is assumed to be strictly dynamically coupled to gas):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}}{\partial t}+{{\rm{\nabla }}}_{p}\cdot \left({{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}{{\boldsymbol{v}}}_{p}\right)=-S({a}_{\max }),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}}{\partial t}+{{\rm{\nabla }}}_{p}\cdot \left({{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}{{\boldsymbol{u}}}_{p}\right)=S({a}_{\max }),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left({{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}{{\boldsymbol{u}}}_{p}\right)+{\left[{\rm{\nabla }}\cdot \left({{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}{{\boldsymbol{u}}}_{p}\otimes {{\boldsymbol{u}}}_{p}\right)\right]}_{p}={{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}\,{{\boldsymbol{g}}}_{p}\\ \quad +\,{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}{{\boldsymbol{f}}}_{p}+S({a}_{\max }){{\boldsymbol{v}}}_{p},\end{array}\end{eqnarray} \tag{ 8 }$$

where u _p are the planar components (r, ϕ) of the grown dust velocity.

The quantity $S({a}_{\max })$ is the rate of small-to-grown dust conversion per disk surface area (in g s⁻¹ cm⁻²). The general idea behind our dust conversion scheme is illustrated in Figure 1, showing the dust size distributions at n and n + 1 time steps with red and blue lines, respectively. The blue area schematically represents the amount of small dust (per surface area) converted into grown dust during one hydrodynamic time step. We note that the area highlighted in orange is transferred above ${a}_{\max }^{n}$, but it does not change the mass of grown dust. The change in the surface density of small dust due to conversion into grown dust ${\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}={{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}^{n+1}-{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}^{n}$ can be expressed as

$$\begin{eqnarray}&&{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}={{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}^{n}\displaystyle \frac{{I}_{1}\left({C}_{\mathrm{sm}}^{n+1}{C}_{\mathrm{gr}}^{n}\,{I}_{2}-{C}_{\mathrm{sm}}^{n}{C}_{\mathrm{gr}}^{n+1}{I}_{3}\right)}{\left({C}_{\mathrm{sm}}^{n+1}{I}_{1}+{C}_{\mathrm{gr}}^{n+1}{I}_{3}\right)\left({C}_{\mathrm{sm}}^{n}{I}_{1}+{C}_{\mathrm{gr}}^{n}{I}_{2}\right)},\end{eqnarray} \tag{ 9 }$$

where ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}^{n}={{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}^{n}+{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}^{n}$ is the total dust surface density, C_sm and C_gr are the normalization constants for small and grown dust size distributions at the current (n) and next (n + 1) time steps, and the integrals I₁, I₂, and I₃ are defined as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{1} & = & {\int }_{{a}_{\min }}^{{a}_{* }}{a}^{3-{\rm{p}}}{da},\\ {I}_{2} & = & {\int }_{{a}_{* }}^{{a}_{\max }^{n}}{a}^{3-{\rm{p}}}{da},\\ {I}_{3} & = & {\int }_{{a}_{* }}^{{a}_{\max }^{n+1}}{a}^{3-{\rm{p}}}{da}.\end{array}\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

By introducing different normalization constants for small and grown dust, we implicitly assume that the dust size distribution can be discontinuous at a_*. In the original paper (see Equation (12) in Vorobyov et al. 2018), we set ${C}_{\mathrm{sm}}^{n}={C}_{\mathrm{gr}}^{n}$ and ${C}_{\mathrm{sm}}^{n+1}={C}_{\mathrm{gr}}^{n+1}$, effectively suggesting that the dust grows in such a manner that the distribution is continuous across a_*. However, the dust surface density in a given computational cell can change not only due to growth (the right-hand side term in Equations (6) and (7)) but also due to dust flow through the cell (the second term on the left-hand side). Because of different dynamics of small and grown dust, a discontinuity may develop at a_* after we account for advection.

To account for this effect, we modify our approach in the following manner. We suggest that the dust size distribution can develop a discontinuity at a_* due to differential dust dynamics. However, dust growth due to the S term smooths out the discontinuity each time it could occur after the advection step. Physically, this assumption corresponds to the dominant role of dust evolution over the dust flow in setting the fixed shape of the dust size distribution. This can be achieved by setting ${C}_{\mathrm{sm}}^{n+1}={C}_{\mathrm{gr}}^{n+1}$ in Equation (9) while keeping ${C}_{\mathrm{sm}}^{n}$ and ${C}_{\mathrm{gr}}^{n}$. The resulting expression is then

$$\begin{eqnarray}&&{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}={{\rm{\Sigma }}}_{\mathrm{tot}}^{n}\displaystyle \frac{{I}_{1}\left({C}_{\mathrm{gr}}^{n}{I}_{2}-{C}_{\mathrm{sm}}^{n}{I}_{3}\right)}{{I}_{4}\left({C}_{\mathrm{sm}}^{n}{I}_{1}+{C}_{\mathrm{gr}}^{n}{I}_{2}\right)},\end{eqnarray} \tag{ 11 }$$

where the integral I₄ is defined as

$$\begin{eqnarray}&&{I}_{4}={\int }_{{a}_{\min }}^{{a}_{\max }^{n+1}}{a}^{3-{\rm{p}}}{da},\end{eqnarray} \tag{ 12 }$$

and the normalization constants can then be found as

$$\begin{eqnarray}&&{C}_{\mathrm{sm}}^{n}=\displaystyle \frac{3{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}^{n}{\rm{\Delta }}s}{4\pi {\rho }_{{\rm{s}}}{I}_{1}},\,\,{C}_{\mathrm{gr}}^{n}=\displaystyle \frac{3{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}^{n}{\rm{\Delta }}s}{4\pi {\rho }_{{\rm{s}}}{I}_{2}},\end{eqnarray} \tag{ 13 }$$

where ρ_s = 3 g cm⁻³ is the material density of dust grains and Δs is the total dust surface area in a given computational cell. Substituting the normalization constants in Equation (11), we finally obtain

$$\begin{eqnarray}&&{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}^{n}{I}_{1}-{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}^{n}{I}_{3}}{{I}_{4}}.\end{eqnarray} \tag{ 14 }$$

The rate of small-to-grown dust conversion during one time step Δt is then written as

$$\begin{eqnarray}&&S({a}_{\max })=-{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}/{\rm{\Delta }}t.\end{eqnarray} \tag{ 15 }$$

To finalize the calculation of $S({a}_{\max })$, the maximum radius of grown dust ${a}_{\max }$ must be computed at each time step and in each computational cell. The evolution of ${a}_{\max }$ is described as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {a}_{\max }}{\partial t}+({u}_{{\rm{p}}}\cdot {{\rm{\nabla }}}_{p}){a}_{\max }={ \mathcal D },\end{eqnarray} \tag{ 16 }$$

where the growth rate ${ \mathcal D }$ accounts for the change in ${a}_{\max }$ due to coagulation, and the second term on the left-hand side accounts for the change of ${a}_{\max }$ due to dust flow through the cell. We write the source term ${ \mathcal D }$ as

$$\begin{eqnarray}&&{ \mathcal D }=\displaystyle \frac{{\rho }_{{\rm{d}}}{\text{}}{v}_{\mathrm{rel}}}{{\rho }_{{\rm{s}}}},\end{eqnarray} \tag{ 17 }$$

where ρ_d is the total dust volume density, ρ_s is the material density of dust grains, and v_rel is the dust-to-dust collision velocity. We consider Brownian and turbulence-induced particle velocities but not the drift velocity, which is of less importance for the grain sizes relevant for our study (see, e.g., Figure 3 in Testi et al. 2014). The adopted approach is similar to the monodisperse model of Stepinski & Valageas (1997) and described in more detail in Vorobyov et al. (2018). The value of ${a}_{\max }$ is capped by the fragmentation barrier (Birnstiel et al. 2012), defined as

$$\begin{eqnarray}&&{a}_{\mathrm{frag}}=\displaystyle \frac{2{{\rm{\Sigma }}}_{{\rm{g}}}{v}_{\mathrm{frag}}^{2}}{3\pi {\rho }_{{\rm{s}}}\alpha {c}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 18 }$$

where v_frag is the fragmentation velocity (see Equation (41)). Whenever ${a}_{\max }$ exceeds a_frag, the growth rate ${ \mathcal D }$ is set to zero. We note that the assumption that any discontinuity in the slope of the dust size distribution (which may develop as a result of dust drift) smooths out owing to dust growth also implies that grown-to-small dust conversion can occur even if ${a}_{\max }\lt {a}_{\mathrm{frag}}$.

As the dust grows, the span in the dust sizes and corresponding Stokes numbers covered by the grown component may become substantial. However, we calculate the stopping time using ${a}_{\max }$, meaning that our model tracks the dynamics of the upper end of the dust size distribution where most of the dust mass reservoir is located if p = 3.5. We introduced an effective Stokes number for the grown dust population in the bidisperse approach in Akimkin et al. (2020). A more rigorous approach to studying dust dynamics requires the use of multiple bins with narrower ranges of dust sizes, and this model is currently under development.

The number of major volatiles in the gas phase and on the icy mantles of the two dust populations is defined by their surface densities: ${{\rm{\Sigma }}}_{s}^{\mathrm{gas}}$ for the gas phase, ${{\rm{\Sigma }}}_{s}^{\mathrm{sm}}$ for the ice on the surface of small dust, and ${{\rm{\Sigma }}}_{s}^{\mathrm{gr}}$ for the ice on the surface of grown dust (s stands for the index of a species). We assume that all of the volatiles are initially on the surfaces of small grains (which are the only grains present at the onset of prestellar core collapse).

The surface densities of the considered species in the gas and ice phases can change via three main physical processes that are considered at each time step and for every grid cell. First, the surface densities ${{\rm{\Sigma }}}_{s}^{\mathrm{gas}}$, ${{\rm{\Sigma }}}_{s}^{\mathrm{sm}}$, and ${{\rm{\Sigma }}}_{s}^{\mathrm{gr}}$ are updated, taking into account freeze-out and sublimation of volatile species for both dust populations. Second, dust evolution (coagulation and fragmentation) is considered, which redistributes the solids (including ices) between small and grown dust fractions. Finally, the surface densities are updated to take into account the transport of volatile species with gas, small dust grains, and grown dust grains. The advection of volatiles is calculated using the same third-order-accurate piecewise parabolic method as for the gas and dust components (Vorobyov et al. 2018). The processes concerning the first and second steps are described in detail below.

Expanding the method from Vorobyov et al. (2013), we describe the phase transitions of volatiles using the following equations:

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Sigma }}}_{s}^{\mathrm{gas}}}{{dt}}=-{\lambda }_{s}{{\rm{\Sigma }}}_{s}^{\mathrm{gas}}+{\eta }_{s}^{\mathrm{sm}}+{\eta }_{s}^{\mathrm{gr}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Sigma }}}_{s}^{\mathrm{sm}}}{{dt}}={\lambda }_{s}^{\mathrm{sm}}{{\rm{\Sigma }}}_{s}^{\mathrm{gas}}-{\eta }_{\mathrm{sm}}^{s},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Sigma }}}_{s}^{\mathrm{gr}}}{{dt}}={\lambda }_{s}^{\mathrm{gr}}{{\rm{\Sigma }}}_{s}^{\mathrm{gas}}-{\eta }_{s}^{\mathrm{gr}},\end{eqnarray} \tag{ 21 }$$

where ${\lambda }_{s}^{\mathrm{sm}}$ and ${\lambda }_{s}^{\mathrm{gr}}$ (s⁻¹) are the adsorption rates onto small and grown dust grains and ${\lambda }_{s}={\lambda }_{s}^{\mathrm{sm}}+{\lambda }_{s}^{\mathrm{gr}}$. We also define ${\eta }_{s}={\eta }_{s}^{\mathrm{sm}}+{\eta }_{s}^{\mathrm{gr}}$ (g cm⁻² s⁻¹), where ${\eta }_{s}^{\mathrm{sm}}$ and ${\eta }_{s}^{\mathrm{gr}}$ are the desorption rates from small and grown dust populations, respectively. Each of the desorption rates is a sum of thermal and photodesorption rates ${\eta }_{s}^{\mathrm{sm}}={\eta }_{s}^{\mathrm{th},\mathrm{sm}}+{\eta }_{s}^{\mathrm{ph},\mathrm{sm}}$ and ${\eta }_{s}^{\mathrm{gr}}\,={\eta }_{s}^{\mathrm{th},\mathrm{gr}}+{\eta }_{s}^{\mathrm{ph},\mathrm{gr}}$. We note that the desorption rates do not depend on ${{\rm{\Sigma }}}_{s}^{\mathrm{sm}}$ and ${{\rm{\Sigma }}}_{s}^{\mathrm{gr}}$, thus implying the zeroth-order desorption model. This approach implies that only the upper layer of the icy mantle can sublimate, which is reasonable for multilayer mantles. Laboratory studies show that in simulated interstellar conditions, the desorption of various molecules is rather described by zeroth-order kinetics (Fraser et al. 2001; Öberg et al. 2005; Bisschop et al. 2006).

For the initial composition of volatiles, we use the ice abundances measured with Spitzer in protostellar cores (Öberg et al. 2011). The initial surface densities of volatile species relative to the initial surface density of gas (${{\rm{\Sigma }}}_{{\rm{g}}}^{\mathrm{init}}$) in the models are shown in Table 2. Although there are estimates of the ice mass being comparable to the rock mass (Pontoppidan et al. 2014), we adopt a total ice mass from Öberg et al. (2011), which is an order of magnitude lower than that of rock, so that we can suggest that icy mantles are thin and do not affect the dust size. As the total surface density of volatiles at the starting moment is only about 8.6% of the dust surface density, this approximation is likely valid along the disk evolution. We will check if at some point the mass of ice notably exceeds that of rock and discuss the possible effects in Section 3.5.

Table 2. Binding Energies, Molecular Weights, and Initial Abundances for the Considered Volatiles

Species	Eb/kB	msp	fice	${{\rm{\Sigma }}}_{s}^{\mathrm{sm}}/{{\rm{\Sigma }}}_{{\rm{g}}}^{\mathrm{init}}$
	(K)	(amu)	(%)
H2O	5770	18	100	3.90 × 10−4
CO2	2360	44	29	2.77 × 10−4
CO	850	28	29	1.76 × 10−4
CH4	1100	16	5	1.74 × 10−5

Note. Binding energies for H₂O, CO₂, and CO are based on the experimental data from Cuppen et al. (2017) for desorption from crystalline water ice, and the binding energy for methane is taken from Aikawa et al. (1996). Initial surface densities of ice on small dust relative to total gas surface density ${{\rm{\Sigma }}}_{s}^{\mathrm{sm}}/{{\rm{\Sigma }}}_{{\rm{g}}}^{\mathrm{init}}$ are based on ice abundances relative to water f_ice determined for low-mass protostars (Öberg et al. 2011).

Download table as:  ASCIITypeset image

The system of Equations (19)–(21) has an analytic solution. We use the explicit expressions for the solution, restricting ice abundances to positive values. The analytic solution, its asymptotic behavior, and the results of tests using a single-point model are presented in Appendix A.

Volatiles can be adsorbed from the gas phase to become icy mantles on the surface of dust grains; they can also evaporate from these mantles by thermal or photodesorption. Here we describe the parameterization of these processes adopted in our model. For simplicity, we drop species index s in this subsection.

2.3.1. Thermal Desorption

The rate coefficient k_td (s⁻¹) for thermal desorption of species from the dust surface can be written as (Tielens & Allamandola 1987)

$$\begin{eqnarray}&&{k}_{\mathrm{td}}=\sqrt{\displaystyle \frac{2{N}_{\mathrm{ss}}{E}_{{\rm{b}}}}{{\pi }^{2}{m}_{\mathrm{sp}}}}\exp \left(-\displaystyle \frac{{E}_{{\rm{b}}}}{{k}_{{\rm{B}}}T}\right),\end{eqnarray} \tag{ 22 }$$

where N_ss (cm⁻²) is the surface density of binding sites, E_b (erg) is the binding energy of the species to the surface, m_sp (g) is the species mass, k_B is the Boltzmann constant, and T is the midplane temperature. In the case of the zeroth-order desorption, the mass desorption rate per disk unit area η_td (g cm⁻² s⁻¹) can be approximated as

$$\begin{eqnarray}&&{\eta }_{\mathrm{td}}=\chi {n}_{\mathrm{act}}{\widetilde{\sigma }}_{\mathrm{tot}}{N}_{\mathrm{ss}}{m}_{\mathrm{sp}}{k}_{\mathrm{td}},\end{eqnarray} \tag{ 23 }$$

where χ is the fractional abundance of the species in the icy mantle, n_act is the number of actively evaporating monolayers, and ${\widetilde{\sigma }}_{\mathrm{tot}}$ (cm² cm⁻²) is the total surface area of dust grains (small or grown) per disk unit area. Typically, two to four upper monolayers can evaporate freely. And, as all considered species have comparable abundances, we assume χ n_act = 1 as a typical constant value across the disk. For the surface density of binding sites, we take N_ss = 10¹⁵ cm⁻², which is a characteristic value for amorphous water ice (Cuppen et al. 2017). The adopted values of binding energies (see Table 2) are based on the experimental data for sublimation from an icy surface (Cuppen et al. 2017). As most of the gas and dust mass lies in high-density midplane regions, we assume equal gas and dust temperatures. The thermal balance includes heating by the stellar and background radiation, dust radiative cooling, and gas adiabatic and viscous terms (Vorobyov et al. 2018). As we track volatiles on small and grown dust, the corresponding values of ${\widetilde{\sigma }}_{\mathrm{tot}}$ are taken for these two populations. For the power-law size distribution (p = 3.5) and dust size not dependent on the vertical height a ≠ a(z), the total small and grown dust surface areas are

$$\begin{eqnarray}&&{\widetilde{\sigma }}_{\mathrm{tot}}^{\mathrm{sm}}=\displaystyle \frac{3{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}}{{\rho }_{{\rm{s}}}\sqrt{{a}_{\min }{a}_{* }}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\widetilde{\sigma }}_{\mathrm{tot}}^{\mathrm{gr}}=\displaystyle \frac{3{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}}{{\rho }_{{\rm{s}}}\sqrt{{a}_{* }{a}_{\max }}}.\end{eqnarray} \tag{ 25 }$$

While the smallest grains in the disk, ${a}_{\min }=0.005$ μm, and the boundary between small and grown dust ensembles, a_* = 1 μm, do not vary along the vertical height, the maximum grain size ${a}_{\max }$ does vary, as dust tends to sediment toward the midplane. Our dust evolution model estimates ${a}_{\max }$ in the disk midplane, where most dust is located. Thus, the value of ${\widetilde{\sigma }}_{\mathrm{tot}}^{\mathrm{gr}}$ can be underestimated by a factor of a few, which is acceptable in the context of our relatively simple dust evolution model.

2.3.2. Photodesorption

In translucent regions, ices can be additionally transferred to the gas phase through photodesorption by stellar or interstellar ultraviolet radiation. In the cold regions of the outer disk and envelope, this effect can be notable. The photodesorption rate coefficient k_pd (s⁻¹) can be estimated as

$$\begin{eqnarray}&&{k}_{\mathrm{pd}}={{YF}}_{\mathrm{UV}}{\sigma }_{{\rm{m}}},\end{eqnarray} \tag{ 26 }$$

where Y (mol photon⁻¹) is the photodesorption yield, F_UV (photons cm⁻² s⁻¹) is the UV photon flux, and σ_m (cm²) is the molecule UV cross section.

The corresponding mass rate of photodesorption per disk unit area η_pd (g cm⁻² s⁻¹) is

$$\begin{eqnarray}&&{\eta }_{\mathrm{pd}}=\chi {n}_{\mathrm{act}}{\widetilde{\sigma }}_{\mathrm{tot}}{N}_{\mathrm{ss}}{m}_{\mathrm{sp}}{k}_{\mathrm{pd}};\end{eqnarray} \tag{ 27 }$$

again assuming χ n_act = 1 and approximating ${\sigma }_{{\rm{m}}}\approx {N}_{\mathrm{ss}}^{-1}$, one gets

$$\begin{eqnarray}&&{\eta }_{\mathrm{pd}}={\widetilde{\sigma }}_{\mathrm{tot}}{m}_{\mathrm{sp}}{{YF}}_{\mathrm{UV}}.\end{eqnarray} \tag{ 28 }$$

We assume that the photodesorption yield is equal to one of water, $Y=3.5\times {10}^{-3}+0.13\exp \left(-336/{T}_{\mathrm{mp}}\right)$ (Westley et al. 1995; note the misprint in the Y₀ value, caption to their Figure 3; see also Walmsley et al. 1999), which is probably the upper boundary for experimentally measured values (Murga et al. 2020).

It is convenient to measure the radiation flux F_UV in the units of a standard UV field, ${F}_{\mathrm{UV}}={F}_{0}^{\mathrm{UV}}G$, where the dimensionless parameter G characterizes the strength of the field. The standard interstellar radiation intensity I₀(E) (photons cm⁻² s⁻¹ sr⁻¹ eV⁻¹) is approximated in the following form (Draine 1978, Equation (11)):

$$\begin{eqnarray}\begin{array}{rcl}{I}_{0}(E) & = & 1.658\times {10}^{6}\left(\displaystyle \frac{E}{\mathrm{eV}}\right)-2.152\times {10}^{5}{\left(\displaystyle \frac{E}{\mathrm{eV}}\right)}^{2}\\ & & +6.919\times {10}^{3}{\left(\displaystyle \frac{E}{\mathrm{eV}}\right)}^{3}.\end{array}\end{eqnarray} \tag{ 29 }$$

The corresponding integrated UV intensity ${I}_{0}^{\mathrm{UV}}$ (photons cm⁻² s⁻¹ sr⁻¹) in the energy range 6.2–13.6 eV (912–2000 Å) is equal to

$$\begin{eqnarray}&&{I}_{0}^{{\rm{U}}{\rm{V}}}={\int }_{6.2\,{\rm{e}}{\rm{V}}}^{13.6\,{\rm{e}}{\rm{V}}}{I}_{0}(E)dE=1.474\times {10}^{7}\,{\rm{p}}{\rm{h}}{\rm{o}}{\rm{t}}{\rm{o}}{\rm{n}}{\rm{s}}\,{{\rm{c}}{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}\,{{\rm{s}}{\rm{r}}}^{-1}.\end{eqnarray} \tag{ 30 }$$

The flux through the unit area from one hemisphere in case of isotropic radiation is equal to

$$\begin{eqnarray}&&{F}_{0}^{{\rm{U}}{\rm{V}}}=\pi {I}_{0}^{{\rm{U}}{\rm{V}}}=4.63\times {10}^{7}\,{\rm{p}}{\rm{h}}{\rm{o}}{\rm{t}}{\rm{o}}{\rm{n}}{\rm{s}}\,{{\rm{c}}{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 31 }$$

In this work, we only consider the interstellar irradiation penetrating the disk and envelope in the vertical direction and neglect the possible input from the central star, as photodesorption is only important at the outer disk boundary, which is shadowed from the star by the disk itself.

We assume a slightly elevated unattenuated interstellar UV field G_env = 5.5G₀ as disks are born in star-forming regions. For the disk midplane, which is illuminated from above and below, this field scales with the UV optical depth τ_UV as

$$\begin{eqnarray}&&{G}_{\mathrm{UV}}=0.5{G}_{\mathrm{env}}{e}^{-{\tau }_{\mathrm{UV}}}.\end{eqnarray} \tag{ 32 }$$

The optical depth in UV toward the disk midplane can be calculated from the surface densities of small and grown dust as

$$\begin{eqnarray}&&{\tau }_{\mathrm{UV}}=0.5\left({\varkappa }_{\mathrm{sm}}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}+{\varkappa }_{\mathrm{gr}}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}\right),\end{eqnarray} \tag{ 33 }$$

where ϰ_sm = 10⁴ and ϰ_gr = 2 × 10² cm² g⁻¹ are typical values for small and grown dust absorption coefficients in the UV (Pavlyuchenkov et al. 2019, Figure 1). The factor 0.5 in Equation (32) reflects our assumption on dominant radiation transfer in the vertical direction (i.e., optical depth in radial direction ≫ optical depth in vertical direction).

2.3.3. Adsorption

The adsorption rate λ(s⁻¹) of a species on a monodisperse noncharged dust grain ensemble with number density n is (Brown & Charnley 1990)

$$\begin{eqnarray}&&\lambda =\pi {a}^{2}n\sqrt{\displaystyle \frac{8{k}_{{\rm{B}}}T}{\pi {m}_{\mathrm{sp}}}}.\end{eqnarray} \tag{ 34 }$$

For a multidispersed dust ensemble, the term π a² n in Equation (34) becomes the total cross section of dust grains in the unit volume, which is four times smaller than ${\widetilde{\sigma }}_{\mathrm{tot}}$ in Equations (24) and (25):

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{\sigma }_{\mathrm{tot}}}{4}\sqrt{\displaystyle \frac{8{k}_{{\rm{B}}}T}{\pi {m}_{\mathrm{sp}}}}.\end{eqnarray} \tag{ 35 }$$

To demonstrate the effect of advection of volatiles on their snowlines, we need to compare the snowline positions obtained in our model with their equilibrium positions. Solutions for Equations (19)–(21) at Δt → ∞ (see Appendix A) for given dust and gas distributions at some disk evolution moment would represent the equilibrium distribution of volatiles, as if timescales of adsorption and desorption were negligibly short compared to local dynamic timescales. We define an equilibrium snowline of a species as a location where its equilibrium gas-phase abundance is equal to the sum of the equilibrium abundances in the ices.

In the presented model, the equilibrium between the gas and ice phases is not necessarily sustained, as we solve time-dependent equations for adsorption and desorption (see Section 2.2). Thus, the obtained distribution of matter is not necessarily consistent with the equilibrium snowline positions. We will use equilibrium snowlines to relate them to the distributions of volatiles and assess the importance of volatile dynamics.

Dust evolution processes lead to small dust conversion to grown dust, and vice versa. Icy mantles should be transported from one dust population to the other as well. We redistribute ices proportionally to dust redistribution; if a certain mass fraction of small dust turns into grown dust, then the same fraction of every ice species on small dust turns into ice on grown dust.

Surface densities of species on dust, ${{\rm{\Sigma }}}_{s,n+1}^{\mathrm{sm}}$ and ${{\rm{\Sigma }}}_{s,n+1}^{\mathrm{gr}}$ at the time step n+1, can be found from their surface densities at time step n,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{s,n+1}^{\mathrm{sm}}={{\rm{\Sigma }}}_{s,n}^{\mathrm{sm}}+{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}{A}_{s},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{s,n+1}^{\mathrm{gr}}={{\rm{\Sigma }}}_{s,n}^{\mathrm{gr}}-{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}{A}_{s},\end{eqnarray} \tag{ 37 }$$

where ΔΣ_d,sm is defined by Equation (9) and A_s is the mass fraction of a given ice on the relevant dust component, which is to be transported:

$$\begin{eqnarray}{A}_{s}=\left\{\begin{array}{cc}{{\rm{\Sigma }}}_{s,n}^{\mathrm{sm}}/{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}^{n}, & \mathrm{if}\ \,{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}\lt 0;\\ {{\rm{\Sigma }}}_{s,n}^{\mathrm{gr}}/{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}^{n}, & \mathrm{if}\ \,{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}\geqslant 0.\end{array}\right.\end{eqnarray} \tag{ 38 }$$

To assess the effect of icy mantles on dust evolution, we use the fragmentation velocity v_frag as a parameter depending on the presence of icy mantles on the surface of grown dust grains. At each time step, after the species surface densities are calculated, we update the values of fragmentation velocity v_frag. To do that, we check whether the total amount of ice at a given location is sufficient to cover all of the grown dust grains present there with at least one monolayer of icy mantle. If so, the dust is treated as "sticky," with high fragmentation velocity; otherwise, it is "bare" and fragile, with lower fragmentation velocity. Although there are laboratory studies that show that water ice might not improve the stickiness of the grains (Musiolik & Wurm 2019), other works show that the presence of icy mantles helps to overcome the fragmentation barrier (Wada et al. 2009; Gundlach & Blum 2015; Okuzumi & Tazaki 2019). In this work, we choose to assume that the icy mantles in general, regardless of their composition, increase v_frag.

The total surface density of the four considered ices on grown dust is

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{ice}}=\displaystyle \sum _{s=1}^{4}{{\rm{\Sigma }}}_{s}^{\mathrm{gr}}.\end{eqnarray} \tag{ 39 }$$

We assign a specific value to v_frag, comparing the ratio between Σ_ice and Σ_d,gr to the threshold value K. To define the value of the threshold K, we use ${{\rm{\Sigma }}}_{\mathrm{ice}}^{\min }$, which is a minimum surface density of ices on grown dust that is sufficient to cover their surfaces with one monolayer. The average size of grown dust can be expressed as $\overline{{a}_{\mathrm{gr}}}=\sqrt{{a}_{* }{a}_{\max }}$. For the monolayer thickness estimated as the size of the water molecule a_ml = 3 × 10⁻⁸ cm and volume densities of ice and dust ρ_ice = 1 and ρ_s = 3 g cm⁻³, the threshold value is

$$\begin{eqnarray}\begin{array}{rcl}K & = & \displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{ice}}^{\min }}{{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}}=\displaystyle \frac{4\pi {\overline{{a}_{\mathrm{gr}}}}^{2}{a}_{\mathrm{ml}}{\rho }_{\mathrm{ice}}}{4/3\pi {\overline{{a}_{\mathrm{gr}}}}^{3}{\rho }_{{\rm{s}}}}=\displaystyle \frac{3{a}_{\mathrm{ml}}{\rho }_{\mathrm{ice}}}{\overline{{a}_{\mathrm{gr}}}{\rho }_{{\rm{s}}}}\\ & = & \displaystyle \frac{3\times {10}^{-8}\,\mathrm{cm}}{\sqrt{{a}_{* }{a}_{\max }}}.\end{array}\end{eqnarray} \tag{ 40 }$$

For example, for the grown dust with a_* = 10⁻⁴ and ${a}_{\max }={10}^{-2}$ cm, we find K = 3 × 10⁻⁵. So, for a reasonable size of grown dust grains, a very small amount of ice is needed to cover them with one molecular layer. Then we use K to define the local fragmentation velocity:

$$\begin{eqnarray}{v}_{\mathrm{frag}}=\left\{\begin{array}{cc}15\,{\rm{m}}\,{{\rm{s}}}^{-1}, & \mathrm{if}\ \,{{\rm{\Sigma }}}_{\mathrm{ice}}/{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}\gt K;\\ 1.5\,{\rm{m}}\,{{\rm{s}}}^{-1}, & \mathrm{if}\ \,{{\rm{\Sigma }}}_{\mathrm{ice}}/{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}\leqslant K.\end{array}\right.\end{eqnarray} \tag{ 41 }$$

The values of the fragmentation velocities for bare and icy dust grains are chosen based on experimental results and typical for icy and bare grains (Wada et al. 2009; Gundlach & Blum 2015).

The evolution of dust and volatiles starts simultaneously with the collapse of the prestellar core. In the predisk phase, the dust growth is slow due to low densities and rare collisions between grains. After t = 55 kyr, the disk is formed, providing a favorable environment for efficient dust growth and initiating the decoupled dynamics of grown dust grains.

Figure 2 shows the distribution of gas and dust and the main dust properties in models 1 (α = 10⁻²) and 2 (α = 10⁻⁴) at three selected time moments. When calculating the dust-to-gas mass ratio, we do not consider ices as part of the dust. Hence, this value in fact represents the rock-to-gas mass ratio (ice masses are not shown). The three time moments shown represent the main stages in disk evolution. At the early time, 152 kyr, the disk is relatively compact and gravitationally unstable. It is highly asymmetric and exhibits a distinct spiral pattern. Later, at an age of 300 kyr, the disk becomes more stable, although weaker spiral arms are still present. At 475 kyr, the disk becomes more uniform and nearly axisymmetric in model 1 and shows only traces of spiral arms in model 2. The gas disk gradually grows in size due to mass loading from the infalling envelope (at the early stages) and spreads out due to viscous spreading (at the later stages). While throughout the disk evolution, some part of the gas is lost in the accretion onto the star and with jets, the dust mass loss is even more profound due to its additional radial drift. This results in a dust-to-gas mass ratio <10⁻² in most of the disk volume at later evolutionary stages.

Zoom In Zoom Out Reset image size

Dust grains grow much faster in the disk than in the envelope. Most of the dust in the disk is found in the form of grown grains, and small dust is depleted. The exception is the very inner 5–20 au of the disk where the fragmentation barrier is low due to the absence of icy mantles. In this region, the dust density can be quite high due to the formation of dust rings (model 2; see Section 3.4). Generally, the dust size increases toward the center of the disk. However, in model 1, there is a sharp drop in dust size in the inner iceless region. In model 2, the dust size reaches its maximum inside the inner dust rings. Grown dust tends to drift toward the star, decreasing the dust-to-gas ratio throughout most of the disk extent and elevating it in the inner parts (see Section 3.2 for the details of dust evolution). Large-scale (tens and hundreds of astronomical units) ringlike structures form at the latest stages, appearing especially bright in the dust-to-gas ratio maps. A more detailed description of the dust ring formation mechanisms can be found in Section 3.4.

A notable difference between the two models is the rate at which the disk spreads out. A more viscous disk in model 1 (α = 10⁻²) tends to spread faster and is larger in size. It also tends to evolve faster and becomes completely axially symmetric at 500 kyr. In model 2 (α = 10⁻⁴), the disk still possesses spiral substructures in both gas and dust at the end of the simulation. The thermal structure of the disks is also affected by the α-value, as viscous heating is one of the main heating mechanisms in the inner 10 au of the disk. As a result, the disk in model 1 is somewhat hotter.

Another important process affected by the α-value is the formation of gravity-bound and pressure-supported clumps, which emerge as a result of gravitational fragmentation in massive disks. Both considered models produce massive gravitationally unstable disks with masses around 0.2 M_⊙; however, in our simulations, long-lived clumps form only at the early stages of disk evolution in model 1. The 152 kyr snapshot for model 2 in Figure 2 also shows signs of clump formation, but the resulting object does not survive for more than 1000 yr. Apart from the gravitational instability described by the Toomre criterion (Toomre 1964), the formation of long-lived clumps requires the cooling rate to be greater than the characteristic growth rate of the gravitational instability (Johnson & Gammie 2003). This condition is usually fulfilled only at large distances from the star, ≳100 au (Rafikov 2005). In the evolutionary stage when the disk is most massive, its size in model 2 is more compact due to lower turbulent viscosity than in model 1. As a result, the regions with appropriate conditions for clump formation are hardly present in model 2.

The ability to compute the long-term evolution of a protoplanetary disk is one of the major advantages of the FEOSAD code. The numerical simulations start with a cloud core collapse and the formation of a gravitationally unstable disk and continue to a 0.5 Myr old axially symmetric disk, which represents a Class II young stellar object. In this section, we consider the azimuthally averaged distribution of gas, dust, and volatiles throughout the disk early history and reveal the long-term effects of icy mantles on the dust evolution.

Figures 3 and 4 present the position–time diagrams showing the azimuthally averaged quantities that describe the evolution of gas and dust for models 1 and 2 (α = 10⁻² and 10⁻⁴, respectively). The radial distributions are shown for all time moments up to 500 kyr. Both long- and short-term trends are observed in the models. The disk forms at around 50 kyr after the onset of gravitational collapse, and during the subsequent ≈450 kyr, the gas disk grows in size and spreads out. Concurrently, the disk loses material through accretion to the star and gradually cools down. During the initial 200 kyr, the disk experiences accretion bursts responsible for short-term disk heating events with a duration of several hundred years. Such bursts are most prominent in model 1. Diluted regions beyond ∼200 au represent the infalling envelope and supply the disk with material during the initial stages of evolution. Dust does not grow above a minimum value of 10⁻⁴ cm in the envelope because of rare collisions and a low fragmentation barrier. The Stokes numbers increase in time in this region, mostly due to decreasing gas densities in the heavily depleted envelope.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the disk where the gas density is high and collisions between grains are frequent, dust grows efficiently up to centimeters in size (more than a decimeter in model 2) and drifts toward the star (Weidenschilling 1977). This leads to a decrease in the dust-to-gas ratio to 10⁻³ and lower in the outer disk after 300 kyr, when the envelope is depleted and the supply of small dust is terminated. However, in the inner ≈10 au of the disk, both small and grown dust accumulate, and the dust-to-gas ratio is elevated. The mechanisms of dust accumulation in the two models are different. In model 1, one of the main factors is a decrease in the Stokes numbers caused by a drop in ${a}_{\max }$. The maximum size of the dust grains drops because the fragmentation velocity is altered by the desorption of ices in this inner disk region. In model 2, dust rings develop in the inner disk due to the formation of a dead zone described in more detail in Section 3.4.

In model 1, the effect of the water snowline on the dust properties is evident. Indicated by the cyan dashed line in Figure 3 is the classic snowline associated with thermal desorption and the temperature gradient in the disk (see Section 3.3.1 for more details on the snowlines). Right inside this snowline, the dust-to-gas ratio is elevated due to the drift of larger dust from outside the snowline. Besides, the content of small dust there is comparable to that of grown dust, and the maximum size of the grains is around 10–100 μm, contrasting with centimeter grains outside the snowline.

This is the effect of fragmentation velocity, which drops to v_frag = 1.5 m s⁻¹ in this region, as dust grains are relieved of their icy mantles in the warm inner disk regions. A factor of 100 contrast in the dust size across the snowline is directly dictated by a factor of 10 change in v_frag (see Equations (18) and (41)). Efficient fragmentation of bare grains lowers the dust size, thus slowing down its radial drift. It also produces more small dust, which is not affected by radial drift. Therefore, both factors, a decrease of ${a}_{\max }$ and increase of the small dust fraction, favor dust accumulation in this region.

Another notable, albeit not very prominent, feature is a thin ring in gas and grown dust, coincident with the position of the water snowline. A mechanism of this ring formation is described in Section 3.4 along with a more thorough illustration of dust and gas distribution around this ring.

The effect of ice mantles on dust in model 1 is prominent due to the principal role of collisional fragmentation in restricting the dust growth. As panel (g) in Figure 3 suggests, in the model with α = 10⁻², the maximum dust size in most of the disk is very close to the dust size limited by the fragmentation barrier a_frag. The solid black lines outline the (orange shaded) disk areas dominated mostly by fragmentation.

It should be noted that at the snowline itself, the dust size does not reach the fragmentation barrier. This effect arises from a sharp change in a_frag combined with azimuthal variations in the gas and dust radial velocities. Although dust, on average, drifts radially inward, it can nevertheless experience cyclic motions around a given radial position during the orbital period. At some azimuths, dust grains with lower ${a}_{\max }$ move from inside the snowline to the region with higher a_frag. Outside the snowline, these grains could have coagulated and grown up to a_frag, if they stayed there long enough. However, they stay there only for a fraction of their Keplerian period (then moving back inside the snowline), which is comparable to the coagulation timescale in this region (∼10 yr), so that the fragmentation limit is mostly not reached. This is a purely multidimensional effect and cannot be observed in 1D disk simulations. More details about azimuthal variations of gas and dust velocities are presented in Appendix B.

In model 2, however, the disk is less turbulent, and a_frag is 2 orders of magnitude higher, as suggested by Equation (18). As a result, collisional fragmentation limits the dust size only inside the water snowline, where v_frag is lowered, while in the bulk of the disk, the radial drift imposes stronger constraints (see panel (g) in Figure 4). Similarly, Stammler et al. (2017) did not see the effect of the CO snowline on dust growth in their modeling due to the dominance of radial drift in this region. There is no sharp change in the dust size at the water snowline, as the border between the regions with drift- and fragmentation-governed dust size is several astronomical units interior to it, particularly at earlier times. Although the transition from drift- to fragmentation-governed regions also results in the change of dust properties, in this case, it is less sharp and does not occur at the snowline (compare the black contours in panel (g) with the inner dust rings in panels (b), (c), and (e) in Figure 4). Besides, the snowline affects the Stokes numbers (panel (h) in Figure 4) and outlines the inner rings at later stages.

Several features of the gas–dust distribution are worth noting. For instance, episodic variations in the gas and dust surface densities are characteristic of model 1. They are best seen in the dust-to-gas ratio distribution (panel (e) in Figure 3). At earlier times (≤200 kyr), they are associated with transient episodes of gravitational instability, which is induced in the inner disk regions by gas–dust clumps that migrate inward. The clump arriving in the inner disk is disintegrated via the action of tidal torques, which increases the density and perturbs the inner disk regions, resulting in the formation of a transient spiral pattern and causing additional inward transport of matter associated with it. The clump also brings dust-rich material from the outer disk. The rise in the protostellar accretion rate caused by the falling clumps creates a luminosity outburst, which heats the disk. Such events are seen as spikes around 150–200 kyr in the gas temperature distribution (see panel (f) in Figure 3) and were first reported by Vorobyov & Basu (2005). The variations in the dust surface density at later stages in model 1 (≥350 kyr) are caused by a "trapping" of drifting dust grains, arising from radial variations in the gas surface density due to the presence of tight spirals in the outer disk. We defer a detailed study of the dust trapping to a follow-up study. In model 2, the prominent feature is a system of dust rings in the inner disk, discussed in more detail in Section 3.4.

3.3.1. Evolution of Volatiles in the Gas and Ices

In this section, we consider the distribution of volatiles in different phases, the positions of their snowlines, and the effects arising at these snowlines due to the collective dynamics of gas and dust. All four species in the gas and on the surface of grown and small dust grains are presented in Figures 5 and 6. The surface densities of the species are shown with respect to the surface density of the corresponding carrier: the species in the gas phase are normalized to the gas surface density, while the species in the ice phase are normalized to the surface density of the dust population on which they reside. Dashed cyan lines show the positions of equilibrium snowlines, which are found using the azimuthally averaged quantities for t → ∞, as described in Section 2.4.

Zoom In Zoom Out Reset image size

The equilibrium snowlines in Figures 5 and 6 reflect sharp changes in the distributions of gas- and ice-phase species. For each species we can indicate two snowlines: the one caused by thermal desorption and the one caused by photodesorption in the outer regions at hundreds of astronomical units from the star. These snowlines are further referred to as thermal and photosnowlines, respectively. The water snowline previously shown in Figures 3 and 4 is the thermal one. While thermal snowlines move closer to the star as the disk cools down due to both accretion luminosity and viscous heating waning in a disk with decreasing density (see, e.g., Makalkin & Dorofeeva 2009), photosnowlines stay at approximately the same distance of several hundred astronomical units. Their positions shift during the first 300 kyr, when the dust surface densities determining UV illumination vary. After the supply of material from the envelope is drained and the dust distributions stabilize, the photosnowlines become steady.

Zoom In Zoom Out Reset image size

Although the volatiles in the model are not necessarily at equilibrium, their resulting distribution is close to that suggested by the equilibrium snowlines. The agreement is quite good for H₂O, CO₂, and CH₄. Only for CO is the mismatch between the equilibrium snowlines and the actual distribution of the species considerable. Here CO is the most volatile of the considered molecules, so its snowlines are virtually absent at the early disk evolution, and they appear only in the later evolution in the outermost disk regions. In the tenuous conditions typical of the outer disk and envelope, all timescales are longer; thus, CO is far from equilibrium.

Although the equilibrium snowlines generally divide the ice- and gas-dominated regions of the disk quite sharply, at certain time moments, the same phase can be present on both sides of the equilibrium snowlines. This is the effect of volatile dynamics, notable for H₂O and CO₂ in Figures 5 and 6. At the early stages, up to 200 kyr, the gas-phase species spread notably through the snowlines into the ice-dominated region. For the thermal snowline, this effect is explained by the presence of warm spiral arms that release water and carbon dioxide vapors in their wakes as they pass through the disk. The positions of equilibrium snowlines are derived using the azimuthally averaged quantities, and the effect of spiral arms is washed out. The shape of snowlines in 2D and the effect of spiral arms on ice distribution are described in more detail in Section 3.3.2.

In the outer regions, interstellar UV radiation responsible for the photodesorption of molecules from the dust surface can also be a source of photodissociation of the gas-phase species, destroying the molecules. This effect is not considered in our model. The gas-phase abundance of H₂O, CO₂, and CH₄ beyond the photosnowline must be much smaller than shown in the left columns of Figures 5 and 6. The only exception is CO, for which the self-shielding may ensure large gas-phase abundances in disk outer regions beyond the photosnowline. To quantify this effect, we show with the white dotted line the locations where the shielding factor θ from van Dishoeck & Black (1988, Table 5) is equal to 0.1. A CO molecule inside this line is likely to survive the photodissociation. Photodissociation could affect ice-phase species as well, but the corresponding reaction rates should probably be lower than those in the gas phase (Murga et al. 2020). The details of this process are still unclear, so we chose to ignore it in our model.

The short-term spikes in the positions of thermal snowlines (most visible in CO₂) are induced by protostellar accretion bursts, which occur in gravitationally unstable disks at the early stages of their evolution (see, e.g., Vorobyov & Basu 2015). In our models, the mass accretion rate is a few × 10⁻⁵ M_⊙ yr⁻¹ during these bursts, corresponding to an increase in the accretion luminosity by a factor of several and the peak accretion luminosities reaching 10–20 L_⊙. While in the quiescent disk, ices are mostly found on grown dust; after the luminosity bursts, the volatiles resettle primarily on small grains, as they dominate in the total surface area, and these ices return to grown dust via coagulation only after tens of kiloyears (see also Section 3.3.2). A more detailed study of the effect of accretion bursts on the evolution of dust and volatiles will be presented in a follow-up study.

Model 1 exhibits variations in the water and carbon dioxide gas abundances in the inner disk. They correlate with the variations in the gas-to-dust ratio (panel (e) in Figure 3). At the same time, the distributions of ices are quite uniform and do not show signs of such variations. These variations are explained by the dynamics of gravitationally unstable disks. Pieces of dust-rich material, such as clumps or rings, migrate toward the star. The amount of dust in them is elevated, but the amount of ice relative to dust does not change, so the ice distribution is unaffected. When the dust-rich material crosses the snowline, the ice evaporates, enriching the gas with the corresponding species. For CO and CH₄, no effect is seen, as the clumps probably originate interior to the methane thermal snowline. In model 2, clump formation is suppressed, as was explained in Section 3.1, so that no variations in the gas-phase volatiles take place.

Another interesting feature of model 2 is multiple thermal snowlines of water and carbon dioxide. Variations in the surface density and temperature arising due to the presence of dust rings inside 10 au after 300 kyr create multiple locations where the equilibrium abundances of gas and ice are equal. Dust is severely depleted between the rings, while the volatiles are still present there thanks to radial transport with the gas. This results in thick icy mantles with the ice surface density reaching and even exceeding the surface density of dust, in this case for both small and grown populations. The formation of the water snowline inside 1 au should be taken with caution. This snowline is formed because of the temperature drop near the sink–disk interface. This drop can be caused by a notable decrease in the disk optical depth thanks to dust depletion in this region. It may, however, be a boundary effect, and further investigation is needed to clarify the nature of the innermost water snowline.

In the vicinity of snowlines, the abundances of volatiles show local peaks, which are characterized by a variety of widths and amplitudes and occur in both the gas and ice phases. There are different mechanisms responsible for these effects, depending on the snowline and the model, including the drift of mantled grains, azimuthally asymmetric radial velocities, and the 2D shape of snowlines (see Section 3.3.2). The effect of ice accumulation right beyond the snowline was first described by Stevenson & Lunine (1988) as a "cold finger" effect (see also Cuzzi & Zahnle 2004); it can also lead to accumulation of refractory material (Cuzzi et al. 2003) in the gas phase.

Accumulation of the gas-phase CO just interior to the snowline is caused by drifting grown dust, which brings solid CO to the inner disk regions, following the evaporation of icy mantles when crossing the corresponding snowline. This effect was also investigated by Stammler et al. (2017) and Krijt et al. (2018), who showed that stronger turbulence leads to weaker enhancement of volatiles because of more efficient smearing of the excess peak by diffusion that scales with α. In our model, the effect of turbulent α is the opposite; namely, the gas-phase abundance peaks at the snowlines of CO and CH₄ are more prominent in model 1 than in model 2 (see Figures 5 and 6).

In the presented simulations, diffusion is not included in the model, so the enhancement in the gas-phase CO and CH₄ is determined by the inward drift of the ice-mantled grains. Around the corresponding snowlines, the radial velocity of dust grains is not directly dependent on α, because the grain size in these regions is determined by radial drift and not by fragmentation. In our model, grown dust drifts faster in the model with α = 10⁻², and more volatiles are accumulated. This is because in model 1 dust size and consequently Stokes number in this region are higher than in model 2.

Despite ignoring turbulent diffusion, our model still displays a diffusion-like redistribution of volatiles, which results in the accumulation of volatiles in the ice phase. This effect originates from azimuthal variations of the radial velocities of gas and dust. Even when the mass is transported inward and azimuthally averaged radial velocities for both dust and gas are negative, the matter can still move locally outward within part of its quasi-Keplerian orbit. These oscillations lead to the effective diffusion of volatiles across the snowline, and the amplitude of velocity variations is typically many times as high as the azimuthally averaged radial velocity (see Appendix B). The effect seems to be stronger for lower α, which could also contribute to the difference in the gas-phase enhancement. A separate study of the accumulation of volatiles at the snowlines is needed to provide a comprehensive description of this effect.

3.3.2. Volatiles in Gravitationally Unstable Disks

Gravitational instability in young massive disks leads to the formation of a distinct spiral structure. At the early evolutionary stages, around 100–200 kyr, the spiral pattern is prominent both in the gas and in dust distributions. It notably affects the disk thermal structure, which is crucial for phase transitions of volatile molecules, and substantially alters the shape of the snowlines. Weaker spirals do not disappear completely even at the late stages and keep affecting gas and dust velocities.

Figure 7 shows the distributions of the two least volatile molecules under consideration, H₂O and CO₂, in the gas and on the surfaces of small and grown dust grains. The inner disk area encompassing the thermal snowlines is shown for a time moment of 100 kyr. Clearly, the equilibrium snowline is not a perfect circle centered on the star. The snowline has a complicated shape that is determined and affected by the spiral pattern. We note that this effect diminishes once the disk gradually stabilizes and becomes more axisymmetric at later evolutionary stages. The complex shape of the snowlines was also noted by Ilee et al. (2017) within a 3D hydrodynamic modeling with chemical evolution. They obtained the deviations of snowline shape from axial symmetry caused by clumps, spirals, and shocks in a gravitationally unstable disk.

Zoom In Zoom Out Reset image size

There is a minor mismatch between the spatial distribution of volatiles and the position of equilibrium snowlines, which are slightly exterior to the sharp changes in the gas and ice abundances of considered species. Thus, there is a contribution of time-dependent processes of adsorption and desorption in the positions of snowlines, albeit their effect is not strong.

In the bulk of the disk, ices are delivered to the surface of grown dust mainly through dust coagulation and growth, as small grains turn into grown ones, carrying along the icy species. Figure 7 demonstrates that in the regions affected by the spirals, the amount of ice on small dust increases, while the corresponding amount on grown dust decreases. As the surface density of small dust in these regions is 2 orders of magnitude lower due to dust growth, it means that most of the considered species are now on the surface of small dust instead of grown dust, while the total amount of ice stays approximately the same.

This peculiar pattern is created by the impact of the spiral density waves that warm up the gas as they run through the disk and sublimate the ices. As the spiral wave retreats, the gas cools down, and the sublimated species freeze back onto dust. The freeze-out occurs predominantly onto the small dust population, which has a greater total surface area despite much lower surface density. Indeed, the ratio between the total surface areas of small and grown dust grains can be expressed from Equations (24) and (25) as $({{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}/{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}})\sqrt{{a}_{\max }/{a}_{\min }}$, which is ≈10 for the conditions between the thermal water snowline and the outer disk boundary. We note that the regions affected by spiral arms are characterized by ${{\rm{\Sigma }}}_{s}^{\mathrm{sm}}/{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{sm}}\sim 1$. Possible complications that can be caused by massive icy mantles are discussed in Section 3.5.

The balance between the adsorption of volatiles to small and grown dust can be altered in favor of grown particles by the Kelvin curvature effect (Ros & Johansen 2013). The saturated vapor pressure is higher for nanometer-sized grains with their curved surface, leading to a lower adsorption rate to small grains. The surface curvature effect is most important when adsorption and desorption are near balance. Ros et al. (2019) argued that around the water snowline, 0.1 mm grains will lose their mantles 100 times as slowly as micron-sized grains. This effect is also responsible for the sintering of ice in dust aggregates (Sirono 2011; Okuzumi et al. 2016). In our model, small particles range from 5 to 1000 nm, so the freeze-out rate could be affected by the Kelvin curvature effect. We do not consider it in our model, but we note that it might lead to diminution of the ice mantles of small grains.

The effect of spiral arms is more pronounced for CO₂ than for H₂O. This is because the spiral arms are strongest in the intermediate disk regions (tens of astronomical units) where the Toomre Q-parameter is smallest. In the innermost disk regions, the Keplerian shear and the temperatures are too high to promote a strong gravitational instability. In the outer disk regions, the gas surface density quickly drops, also reducing the strength of the gravitational instability. As a result, the spiral arms efficiently heat the intermediate disk regions to the level that causes CO₂ desorption, but their effect is reduced in the inner disk, and water is less affected. In model 2 (α = 10⁻⁴), the disk is systematically colder because of smaller viscous heating, and the snowlines move closer to the star. As a consequence, the effect of heating by spiral arms is less pronounced.

In both models, we obtain multiple dust rings in the simulation results. In model 1, there is a thin and faint ring in dust and gas, situated precisely at the thermal water snowline (see panels (a) and (b) in Figure 3). There are also multiple dust rings descending to the star after 400 kyr in this model (see panels (b), (c), and (e) in Figure 3). In model 2, a system of prominent dust rings is created inside 10 au that persists throughout the disk evolution (see panels (b), (c), and (e) in Figure 4). All of these rings are created by various mechanisms that we briefly discuss in this section.

The thin ring in grown dust and gas, coinciding with the position of the thermal snowline of water, is connected to the volatiles and their effect on the fragmentation velocity. An example of the radial cross section of this ring is shown in Figure 8. At this snowline, the fragmentation velocity changes and the dust maximum size decreases, resulting in a sharp drop in the Stokes number. With decreasing Stokes, the inward drift of grown dust grains slows down and creates a factor of ∼1.5 increase in the grown dust surface density. Further accumulation is restricted by collisional fragmentation and transformation of grown dust into small dust. The formation of dust rings at the snowlines, which may assist in planetesimal formation, is described in many works (Cuzzi & Zahnle 2004; Brauer et al. 2008; Drążkowska & Alibert 2017). In our modeling, the ring due to this effect is only produced in the high-viscosity model, and the increase in dust density is modest.

Zoom In Zoom Out Reset image size

Large-scale rings that gradually migrate to the star at the late stages of disk evolution in model 1 are likely formed thanks to the effect known as a "traffic jam." Figure 9 illustrates the effect by showing the azimuthally averaged quantities at t = 475 kyr. The steplike radial distribution of the grown dust surface density represents the dust rings (which are best seen in the dust-to-gas ratio in panel (e) of Figure 3). We note that the gas surface density does not show such an expressed steplike behavior, although the gradient of Σ_g does show radial variations. As a result, the gas pressure gradient that controls the dust drift velocity also shows local maxima and minima, which correlate with the radial velocity of grown dust relative to that of gas (u_r − v_r ). These variations in the relative dust velocity create traffic jams, which result in the formation of the steplike distribution of Σ_d,gr. The maximum sizes of the dust grains and the Stokes numbers are consistent with the expectations of the efficient radial drift of dust grains (e.g., Birnstiel et al. 2016). The radial variations in the gas surface density are caused by weak spiral density waves. It is not yet clear why similar ring structures do not form in model 2. This low-α model has a stronger spiral pattern, which may smear out radial variations. We plan to explore this mechanism in more detail in follow-up studies.

Zoom In Zoom Out Reset image size

Prominent rings in the inner disk regions in model 2 are caused by the formation of a dead zone. The matter that is transported from the disk outer regions mainly by the action of gravitational torques in this low-α model hits a bottleneck in the innermost disk regions where the strength of the spiral pattern diminishes because of a rising disk temperature. We demonstrate this effect by calculating the α-parameter caused by gravitational instability. Following Kratter & Lodato (2016), we express this quantity as

$$\begin{eqnarray}&&{\alpha }_{\mathrm{GI}}={\left|\displaystyle \frac{d\mathrm{ln}{{\rm{\Omega }}}_{{\rm{K}}}}{d\mathrm{ln}r}\right|}^{-1}\displaystyle \frac{1}{4\pi G}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial r}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial \phi },\end{eqnarray} \tag{ 42 }$$

where G is the gravitational constant, Φ is the gravitational potential, and Ω_K is the Keplerian angular velocity.

Panel (i) in Figures 3 and 4 presents the spacetime diagrams of the grown dust surface density and total effective α-parameter (defined as the sum α_eff = α + α_GI) for both considered models. Model 2 is characterized by an α_eff that is a strong function of radial distance with a deep minimum in the inner 5–10 au. This is the region where prominent rings form in the grown dust density distribution. In model 1, the effective α-value is mostly set by the spatially and temporally constant viscous α = 10⁻². As a result, this model lacks a deep minimum in α_eff in the inner disk regions, and the dead zone does not form. Other mechanisms lead to the formation of weaker and more dynamic dust rings in this case, as discussed above.

Ringlike structures in dust continuum emission have been observed with high resolution in many protoplanetary disks (ALMA Partnership et al. 2015; Huang et al. 2018; Long et al. 2018; Cieza et al. 2021), although the apparent width and relative brightness of the rings would also be affected by the dust opacity variations (Akimkin et al. 2020). Nevertheless, among the observed millimeter continuum images of protoplanetary disks, we can identify some structures analogous to the weak-contrast dust rings obtained in our models (see Figure 2) based on the geometry of the rings and assuming that the apparent contrast in brightness is proportional to the surface density of grown dust. There are many analogs of the large-scale, evenly spaced multiple rings seen in model 1, for example, HL Tau and HD 143006, the inner rings in RU Lup and AS 209, and the outer rings in DoAr 25, DL Tau, and GO Tau. A single wide ring at the outer disk edge with a shallow broad gap around a central bright region, similar to the one we obtained in model 2, is observed in Elias 20, Sz 129, and FT Tau.

The evolution of refractory grains and volatiles is closely linked. Icy mantles control the fragmentation velocity and thus the grain size and dynamics. On the other hand, grains of sufficiently large size (and Stokes number) drift radially and redistribute volatiles within the disk. An important parameter for the volatile transport is the ice-to-rock mass ratio in a grain. It can vary in both space and time, as well as depend on the grain size. In this subsection, we focus on water ice, specifically for the temperatures close to the water freeze-out temperature (≈150 K). We also consider the abundances of other molecules in the icy mantles and the evolution of ice composition.

Ice deposition on a grown grain can occur either via direct freeze out or in the coagulation process. In the first case, the mass of the icy mantle should be proportional to the grain surface area, so the ice-to-rock mass ratio is inversely proportional to the grain size. If ice is deposited onto dust grains solely by the coagulation with other grains, the ice-to-rock ratio will be constant for any grain size and equal to the initial one. In real disk conditions, both processes are in action, so the self-consistent treatment of dust and volatile evolution is required to properly assess the ice-to-rock ratio. This is especially important near the condensation fronts, where bidirectional transport of volatiles and grains may significantly redistribute the ice between grains (Stevenson & Lunine 1988; Cuzzi & Zahnle 2004; Drążkowska & Alibert 2017).

Figure 10 shows the ice-to-rock mass ratio for water ice as a function of a_max for models 1 (upper panel) and 2 (lower panel) at 300 kyr. The ratios for small (${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}{\rm{O}}}^{\mathrm{sm}}$/Σ_d,sm) and grown (${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}{\rm{O}}}^{\mathrm{gr}}$/Σ_d,gr) dust are included. Each dot on the plots represents some spatial location within the disk, while the dot color encodes the local disk temperature. The locations with low ice surface densities (<10⁻⁴ g cm⁻²) are marked by pale dots. In both models, there is a significant population of grown dust grains that are relatively cold (blue dots; T < 100 K) and have mantles with the initial ice-to-rock mass ratio ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}{\rm{O}}}^{\mathrm{gr}}/{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{gr}}=3.9\times {10}^{-2}$ (except for the regions where the dust content is very low, marked by pale blue dots). This indicates that the dominant process of ice deposition is coagulation rather than condensation on already-grown grains. However, there is a prominent fraction of grains with ice-to-rock ratios significantly higher or lower than the initial one for both small and grown populations. Their spatial location corresponds to either the water condensation front (model 1) or prominent rings in dust distribution (model 2). The dashed line in both panels shows a threshold corresponding to one monolayer of ice. Grains below the dashed line are not entirely covered with ice and thus have a lower fragmentation barrier.

Zoom In Zoom Out Reset image size

In model 1 with a higher α-parameter, the grown grain size at the water snowline is governed by collisional fragmentation, so the change of fragmentation velocity due to the presence of ice is of relevance. The grains just outside the snowline are hot and relatively small and have a very low water ice-to-rock ratio. The typical ice-to-rock ratio in the model 1 disk corresponds to the cold grains well beyond the snowline and is equal to the initial one. In model 2, with a smaller α-parameter, the dust size at the water snowline is limited by the radial drift rather than by fragmentation. The largest icy grains in the model 2 disk are located just outside the snowline and thus are among the warmest ones. There is an additional population of hot millimeter-sized grains (red dots at ≈10⁻¹–10¹ cm) that belong to the highly dynamical region around the hot dust ring at 2 au. Both populations may have an ice-to-rock ratio up to an order of magnitude higher that the initial one. Note also that the typical temperatures of icy grains can exceed 170 K; due to high density, the characteristic freeze-out temperature in the corresponding disk regions is higher.

Small icy grains demonstrate an even wider range of ice-to-rock ratios. The excess of water in their mantles results from the recondensation of water that is brought from inside the snowline. In our modeling, we assume that the icy mantle does not affect the grain cross section. This assumption holds in most of the disk but appears to break down in the vicinity of snowlines where the ice-to-rock ratio ≫1 due to the accumulation of ices on small dust. While in model 1, the water ice-to-rock ratio does not exceed 1, in model 2, it reaches ∼300 for small grains. For an ice-to-rock ratio of ∼1 and an ice density three times lower than that of rock, the presence of a mantle would increase the grain size only by a factor of $\sqrt[3]{4}\approx 1.6$, which is still of the order of unity. For all grown dust, the relative amount of ice does not exceed 3, so the condition of the mantles not contributing to the dust size is not severely violated. For small grains with ice-to-rock ratios of 20 and 300, the increase in size would be ≈4 and ≈10, respectively. Such an increase is harder to ignore. However, the mass fraction of small dust is notably lower than the one of grown dust, so most collisions occur with grown grains rather than small grains.

The presence of thick icy mantles could affect the dynamics and growth rate of dust, which should be further investigated. Figure 10 only presents ice-to-rock ratios for water, but other volatiles may have noticeable variations in the amount of ice, too. It is especially true for carbon dioxide, for which the effect of accumulation of ices on small dust is also very strong. The regions where different ices dominate over rock would generally not overlap with one another, as the effect works in a limited space near the snowline, so the thick mantles are nearly exclusively composed of one molecule.

Other studies also reported ices mostly residing on smaller grains around the snowline. Stammler et al. (2017) found that outside the CO snowline, ice mantles become thicker, and small grains contain a larger fraction of ice (ice exceeding rock by 2–3 orders of magnitude) than grown grains (ice-to-rock ratio of 0.5–10, depending on the dust size). They also showed that the accumulation region is wider for small grains. Similarly, Krijt et al. (2018) found that when accumulating outside the snowline, CO ice resides on small dust grains, rather than grown pebbles.

The composition of icy mantles changes throughout the disk and evolves significantly compared to the initial one. The ratio of ice surface density to that of gas for different species is shown in Figure 11. Despite the accumulation of ices on small dust, ice surface densities appear to be only a small fraction of the matter in the disk when normalized to gas surface density and azimuthally averaged. A typical ice-to-gas mass ratio is 0.02%–0.08%, only reaching ≈1% in a region near the inner rings (see the out-of-scale peaks in model 2 at 500 kyr in Figure 11). The grains with thick mantles only present a very small subpopulation of the dust, but a more sophisticated model accounting for the presence of these thick mantles is needed to understand their impact on dust evolution.

Zoom In Zoom Out Reset image size

At the earlier stage (100 kyr), ices on small and grown dust are notoriously separated in Figure 11, with ices on grown dust dominating inside ≈100 au and ices on small dust being abundant in the exterior disk parts. Initially deposited on the surfaces of small grains, ices are transferred to grown dust mostly through coagulation, which is not efficient outside 100 au at early evolutionary stages. With time, the disk spreads out and dust growth takes over in most of the ice-dominated space, so at later stages, the ices predominantly reside on the grown grains. At 300 and 500 kyr, only a fraction of ice remains on small dust in the regions beyond 100 au.

The total amount of ice in the disk decreases with time, as drifting grains bring the volatiles to the inner disk, where they turn to the gas and then accrete to the star. The disk global ice-to-gas mass ratio decreases by ≈2 times during the first 0.5 Myr. At later times, most of the ice is on grown grains, as most of the rock material has evolved into grown dust. On both small and grown dust, the amount of ice varies with the distance from the star.

The presence of (thermal) snowlines is reflected in the composition of ices; in the inner disk, the mantles are made of water, then as we move out CO₂ is added, and after that CH₄ and CO are added in the outer regions. The ratios between ices evolve, too. If, at early times, they are close to the initial abundances from Table 2, later, the fraction of CO₂ and CO grows and can exceed the one of water. For example, at 500 kyr at 100–200 au, half of all ice on grown dust is CO ice, and the mantles on small grains consist almost entirely of this species. Carbon dioxide dominates in the composition of mantles between 10 and 30 au in model 2. The amount of CO ice increases with disk age, with most of the CO staying in the gas phase during the first ≈200–300 kyr of disk evolution, as also found in observations (van't Hoff et al. 2020).

There are studies considering the chemical evolution of ices in disks with dust dynamics and evolution. Detailed chemical modeling of quasi-stationary disks shows that at megayear timescales, CO can be chemically depleted from both gas and ice phases by efficient reprocessing to CO₂, CH₄, and CH₃OH (Bosman et al. 2018; Eistrup et al. 2018) or under significant ionization by cosmic-ray and UV radiation (Schwarz et al. 2018, 2019). Booth & Ilee (2019) showed that in viscous disks (α > 10⁻³), pebble drift dominates the chemical depletion of CO, CH₄, and CO₂, although the high cosmic-ray ionization rate (ζ_CR = 10⁻¹⁷) is able to make chemical processing more efficient (as also pointed out by Bosman et al. 2018). The combination of chemical modeling with dust transport suggests CO depletion from the gas by 2 orders of magnitude (Krijt et al. 2020). We follow the disk evolution only for several hundred kiloyears, but it is possible that chemical reprocessing is already in action at these early stages. Adding the chemical evolution of volatiles would be a natural development of the presented model.