SINTERING-INDUCED DUST RING FORMATION IN PROTOPLANETARY DISKS: APPLICATION TO THE HL TAU DISK

We construct a radial temperature profile T(r) of the HL Tau disk based on the data of the surface brightness profiles provided by ALMA Partnership et al. (2015). We deproject the intensity maps of the disk's dust continuum at ALMA Bands 3, 6, and 7 into circularly symmetric views, assuming the disk inclination angle of $46\_\_AMP\_\_fdg;7$ and the position angle of 138^◦ (ALMA Partnership et al. 2015). We then obtain the radial profiles of the intensities I_ν by azimuthally averaging the deprojected images. The upper panel of Figure 1 shows the derived radial emission profiles. Here, the intensities are expressed in terms of the Planck brightness temperature ${T}_{{\rm{B}}}$. Shown in the lower panel is the spectral index between Bands 6 and 7, ${\alpha }_{{\rm{B}}6-{\rm{B}}7}\equiv \mathrm{ln}({I}_{{\rm{B7}}}/{I}_{{\rm{B6}}})/\mathrm{ln}({\nu }_{{\rm{B7}}}/{\nu }_{{\rm{B6}}})$, where ${\nu }_{{\rm{B6}}}=233.0\;{\rm{GHz}}$ and ${\nu }_{{\rm{B7}}}=343.5\;{\rm{GHz}}$ are the frequencies at Bands 6 and 7, respectively.

As noted by ALMA Partnership et al. (2015), the HL Tau disk has a pronounced central emission peak at $\lesssim 10\;{\rm{au}}$ and three major bright rings at ∼20, 40, and 80 au. The central emission peak and two innermost bright rings have a spectral index of ${\alpha }_{{\rm{B}}6-{\rm{B}}7}\approx 2$. In general, this indicates that the emission at these wavelengths is either optically thick, or optically thin but from dust particles larger than millimeters. The brightness temperature is equal to the gas temperature in the former case and is lower in the latter case. While Zhang et al. (2015) adopted the latter interpretation, we now pursue the former interpretation. Specifically, we assume that the Band 7 emission is optically thick and hence $T(r)={T}_{{\rm{B}}}$ at the center, 20 au, and 40 au. The simplest profile satisfying this assumption is the single power law

$$\begin{eqnarray}&&T(r)=310{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-0.57}\;{\rm{K}},\end{eqnarray} \tag{ 1 }$$

which is shown by the dashed line in the upper panel of Figure 1. We adopt this temperature profile in this study.

The density structure of the HL Tauri gas disk is unknown. Therefore, we simply assume that the gas surface density ${{\rm{\Sigma }}}_{{\rm{d}}}(r)$ obeys a power law with an exponential taper (Hartmann et al. 1998; Kitamura et al. 2002; Andrews et al. 2009),

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}(r)=\displaystyle \frac{(2-\gamma ){M}_{{\rm{disk}}}}{2\pi {r}_{c}^{2}}{\left(\displaystyle \frac{r}{{r}_{{\rm{c}}}}\right)}^{-\gamma }\mathrm{exp}\left[-{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{2-\gamma }\right],\end{eqnarray} \tag{ 2 }$$

where ${r}_{{\rm{c}}}$ and ${M}_{{\rm{disk}}}$ are a characteristic radius and the total mass of the gas disk, respectively, and γ ($0\lt \gamma \lt 2$) is the negative slope of ${{\rm{\Sigma }}}_{{\rm{g}}}$ at $r\ll {r}_{{\rm{c}}}$. We take $\gamma =1$ as the fiducial value but also consider $\gamma =0.5$ and 1.5. The dependence of our simulation results on γ will be studied in Section 6.1. The values of ${M}_{{\rm{disk}}}$ and ${r}_{{\rm{c}}}$ are fixed to $0.2{M}_{\odot }$ and $150\;\mathrm{au}$, respectively. The adopted disk mass is about twice the upper end of the previous mass estimates for the HL Tau disk (Guilloteau et al. 2011; Kwon et al. 2011, 2015). We assume such a massive disk because the dust mass in the disk decreases with time due to the radial drift of dust particles (see Section 5.1). Figure 2 shows the surface density profiles for $\gamma =1$, 0.5, and 1.5.

Zoom In Zoom Out Reset image size

Because the disk is assumed to be vertically isothermal, the vertical distribution of the gas density obeys a Gaussian with the midplane value ${\rho }_{{\rm{g}}}={{\rm{\Sigma }}}_{{\rm{g}}}/(\sqrt{2\pi }{H}_{{\rm{g}}})$, where ${H}_{{\rm{g}}}={c}_{s}/{\rm{\Omega }}$ is the gas scale height, ${c}_{{\rm{s}}}$ is the sound speed, and Ω is the Keplerian frequency. The isothermal sound speed is given by ${c}_{{\rm{s}}}=\sqrt{{k}_{{\rm{B}}}T/{m}_{\mu }}$, where ${k}_{{\rm{B}}}$ is the Boltzmann constant and m_μ is the mean molecular mass of the disk gas assumed to be ${m}_{\mu }=2.3\;{\rm{amu}}$. The Keplerian frequency is given by ${\rm{\Omega }}=\sqrt{{{GM}}_{*}/{r}^{3}}$, where G is the gravitational constant and M_* is the stellar mass. We adopt ${M}_{*}=1{M}_{\odot }$ so that the sum of the stellar and disk masses in the fiducial model, ${M}_{*}+{M}_{{\rm{disk}}}=1.2{M}_{\odot }$, is within the range of the previous estimates for the HL Tau system (Beckwith et al. 1990; Sargent & Beckwith 1991; ALMA Partnership et al. 2015).

We note here that the assumed disk model is marginally gravitationally stable: the Toomre stability parameter $Q\equiv {c}_{{\rm{s}}}{\rm{\Omega }}/(\pi G{{\rm{\Sigma }}}_{{\rm{g}}})$ satisfies $Q\gtrsim 1$ at all radii for all choices of γ. This can be seen in Figure 2, where the surface density corresponding to Q = 1 is shown by the dashed line.

We assume that the volatile composition of the HL Tau disk is similar to that of comets in our solar system. We select six major cometary volatiles in addition to ${{\rm{H}}}_{2}{\rm{O}}$ and take their abundances relative to ${{\rm{H}}}_{2}{\rm{O}}$ to be consistent with cometary values (Mumma & Charnley 2011). The volatiles we select are ammonia (${\mathrm{NH}}_{3}$), carbon dioxide (${\mathrm{CO}}_{2}$), hydrogen sulfide (${{\rm{H}}}_{2}{\rm{S}}$), ethane (${{\rm{C}}}_{2}{{\rm{H}}}_{6}$), methane (${\mathrm{CH}}_{4}$), and carbon monoxide (CO). We neglect another equally abundant species, methanol (${\mathrm{CH}}_{3}\mathrm{OH}$), because the snow line of ${\mathrm{CH}}_{3}\mathrm{OH}$ is very close to that of more abundant ${{\rm{H}}}_{2}{\rm{O}}$ (Sirono 2011b). Table 1 lists the abundances we adopt and the observed ranges of cometary abundances taken from Mumma & Charnley (2011).

Table 1.  Abundances of Major Cometary Volatiles Relative to ${{\rm{H}}}_{2}{\rm{O}}$ (in Percent)

Species	Cometary Valuea	Adopted Value fj
${{\rm{H}}}_{2}{\rm{O}}$	100	100
${\mathrm{NH}}_{3}$	0.2–1.4	1
${\mathrm{CO}}_{2}$	2–30	10
${{\rm{H}}}_{2}{\rm{S}}$	0.12–1.4	1
${{\rm{C}}}_{2}{{\rm{H}}}_{6}$	0.1–2	1
${\mathrm{CH}}_{4}$	0.4–1.6	1
CO	0.4–30	10

Note.

^aMumma & Charnley (2011).

Download table as:  ASCIITypeset image

The equilibrium vapor pressures of volatiles determine the temperatures at which sublimation and sintering occurs. In this study, we approximate the equilibrium vapor pressure for each volatile species j by the Arrhenius form

$$\begin{eqnarray}&&{P}_{{\rm{ev}},j}=\mathrm{exp}\left(-\displaystyle \frac{{L}_{j}}{T}+{A}_{j}\right)\;\mathrm{dyn}\;{\mathrm{cm}}^{-2},\end{eqnarray} \tag{ 3 }$$

where L_j is the heat of sublimation in Kelvin and A_j is a dimensionless constant. Table 2 summarizes the values of L_j and A_j for the seven volatile species considered in this study. For ${\rm{CO}}$, ${\mathrm{CH}}_{4}$, ${\mathrm{CO}}_{2}$, and ${\mathrm{NH}}_{3}$, we follow Sirono (2011b) and take the values from Table 2 of Yamamoto et al. (1983). The values for ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ are derived from the analytic expression of the vapor pressure by Moses et al. (1992, see their Table III; note that we here neglect the small offset of T in their original expression). For ${{\rm{H}}}_{2}{\rm{S}}$, we determined L_j and A_j by fitting Equation (3) to the vapor pressure data provided by (Haynes 2014, page 6–92). Figure 3 shows ${P}_{{\rm{ev}},j}$ of the seven volatile species as a function of r for the temperature distribution given by Equation (1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Vapor Pressure Parameters

Species	Lj (K)	Aj	Ref.	${L}_{j,\mathrm{tuned}}$ (K)
${{\rm{H}}}_{2}{\rm{O}}$	6070	30.86	1	5463
${\mathrm{NH}}_{3}$	3754	30.21	2	3379
${\mathrm{CO}}_{2}$	3148	30.01	2	⋯
${{\rm{H}}}_{2}{\rm{S}}$	2860	27.70	3	⋯
${{\rm{C}}}_{2}{{\rm{H}}}_{6}$	2498	30.24	4	2248
${\mathrm{CH}}_{4}$	1190	24.81	2	⋯
CO	981.8	26.41	2	⋯

References. (1) Bauer et al. (1997), (2) Yamamoto et al. (1983), (3) Haynes (2014), (4) Moses et al. (1992).

Download table as:  ASCIITypeset image

Strictly speaking, the vapor pressure data given in Table 2 only apply to pure ices. In protoplanetary disks, volatiles may be trapped inside the ${{\rm{H}}}_{2}{\rm{O}}$ mantle of dust grains instead of being present as pure ices. If this is the case, the volatiles would sublimate not only at the sublimation temperatures for pure ices but also at higher temperatures where monolayer desorption from the ${{\rm{H}}}_{2}{\rm{O}}$ substrate, phase transition of ${{\rm{H}}}_{2}{\rm{O}}$ ice, or codesorption with the ${{\rm{H}}}_{2}{\rm{O}}$ ice takes place (Collings et al. 2003, 2004; Martín-Doménech et al. 2014). However, all these high-temperature desorption processes are irrelevant to neck formation (sintering) because the desorbed molecules are unable to recondense onto grain surfaces at such high temperatures. By using the vapor pressure data for pure ices, we effectively neglect all these desorption processes.

When estimating the locations of the snow lines, it is important to note that the vapor pressure data in the literature are subject to small but still non-negligible uncertainties. For example, the values of the sublimation energies we use for ${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{NH}}_{3}$, ${\mathrm{CO}}_{2}$, and ${\rm{CO}}$ are 10%–20% higher than those derived from the very recent temperature programmed desorption experiments by Martín-Doménech et al. (2014, see their Table 4). Luna et al. (2014) compiled the sublimation energies of major cometary volatiles from different experimental methods, and showed that the published sublimation energies have a standard deviation of 14%, 8%, 11%, and 8% for ${\mathrm{NH}}_{3}$, ${\mathrm{CO}}_{2}$, ${\mathrm{CH}}_{4}$, and ${\rm{CO}}$, respectively (see their Table 2 and Figures 4 and 5). Such uncertainties might be present in the vapor pressure data for other volatiles species. As we demonstrate in Section 6.4, even a 10% uncertainty in L_j can lead to a 20%–30% uncertainty in the location of its snow line, and a better match between our simulation results and the ALMA observation can be achieved if the sublimation energies of ${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{NH}}_{3}$, and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ are taken to be 10% lower than the fiducial values. We denote these tuned sublimation energies by ${L}_{j,\mathrm{tuned}}$ (see Table 2). For ${{\rm{H}}}_{2}{\rm{O}}$ and ${\mathrm{NH}}_{3}$, the tuned sublimation energies are more consistent with the results of Martín-Doménech et al. (2014).

For each volatile species j, we define the snow line as the orbit inside which the equilibrium pressure ${P}_{{\rm{ev}},j}$ exceeds the partial pressure P_j. Assuming that the disk gas is well mixed in the vertical direction, P_j is related to the surface number density of j-molecules in the gas phase, N_j, as

$$\begin{eqnarray}&&{P}_{j}=\displaystyle \frac{{N}_{j}{k}_{{\rm{B}}}T}{\sqrt{2\pi }{H}_{{\rm{g}}}}.\end{eqnarray} \tag{ 4 }$$

In this study, we do not directly treat the evolution of N_j, but instead estimate them by assuming that the ratio between N_j and the surface number density of ${{\rm{H}}}_{2}{\rm{O}}$ molecules in the solid phase is equal to the cometary abundance f_j given in Table 1. We also assume that the mass fraction of ${{\rm{H}}}_{2}{\rm{O}}$ ice inside the aggregates is 50%. Under these assumptions, N_j can be expressed as

$$\begin{eqnarray}&&{N}_{j}=\displaystyle \frac{0.5{f}_{j}{{\rm{\Sigma }}}_{{\rm{d}}}}{{m}_{\mathrm{mol},{{\rm{H}}}_{2}{\rm{O}}}}.\end{eqnarray} \tag{ 5 }$$

The adopted relative abundances give the relations ${P}_{j}=0.1{P}_{{{\rm{H}}}_{2}{\rm{O}}}$ for $j={\rm{CO}}$ and ${{\rm{CO}}}_{2}$, and ${P}_{j}=0.01{P}_{{{\rm{H}}}_{2}{\rm{O}}}$ for $j={\mathrm{NH}}_{3}$, ${{\rm{H}}}_{2}{\rm{S}}$, ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$, and ${\mathrm{CH}}_{4}$. For the purpose of calculating the locations of the snow lines, the simplification made here is acceptable as a first-order approximation, because the locations of the snow lines are predominantly determined by the strong dependence of ${P}_{{\rm{ev}},j}$ on T(r) and are much less sensitive to a change in P_j.

To see the approximate locations of the snow lines in our HL Tau disk model, we temporarily assume the standard dust-to-gas mass ratio of ${{\rm{\Sigma }}}_{{\rm{d}}}=0.01{{\rm{\Sigma }}}_{{\rm{g}}}$ throughout the disk. The dashed, dotted–dashed, and dotted lines in Figure 3 show the partial pressure curves ${P}_{j}=\{1,0.1,0.01\}{P}_{{{\rm{H}}}_{2}{\rm{O}}}$ for the fiducial disk model ($\gamma =1$). For each volatile species, the location of the snow line is given by the intersection of ${P}_{{\rm{ev}},j}$ and P_j, which is indicated by a filled circle in Figure 3. In this example, the ${{\rm{H}}}_{2}{\rm{O}}$ snow line lies at a radial distance of 3.5 au, which is well interior to the innermost dark ring of the HL Tau disk lying at ∼ 13 au (see Figure 1). Zhang et al. (2015) suggested that the ${{\rm{H}}}_{2}{\rm{O}}$ snow line lies on the 13 au dark ring, assuming a higher gas temperature than ours. The snow lines of ${\mathrm{NH}}_{3}$, ${\mathrm{CO}}_{2}$, ${{\rm{H}}}_{2}{\rm{S}}$, and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ are narrowly distributed over the intermediate region of 10–30 au, and those of ${\mathrm{CH}}_{4}$ and ${\rm{CO}}$ are located in the outermost region of 100–150 au.

Sintering is the process of neck growth, and its timescale is inversely proportional to the rate at which the neck radius increases (e.g., Swinkels & Ashby 1981). The timescale depends on the size of monomers, with larger monomers generally leading to slower sintering. In this study, we simply assume monodispersed monomers and treat their radius a₀ as a free parameter (see Section 6.3 for parameter study). We only consider ${a}_{0}\leqslant 4\;\mu {\rm{m}}$ because sintering is too slow to affect dust evolution in a protoplanetary disk beyond this size range (see below).

When neck growth is driven by vapor transport of volatile j, its timescale is given by (Sirono 2011b)

$$\begin{eqnarray}&&{t}_{{\rm{sint}},j}=4.7\times {10}^{-3}\displaystyle \frac{{a}_{0}^{2}{(2\pi {m}_{{\rm{mol}},j})}^{1/2}{({k}_{{\rm{B}}}T)}^{3/2}}{{P}_{{\rm{ev}},j}(T){V}_{{\rm{mol}},j}^{2}{\gamma }_{j}},\end{eqnarray} \tag{ 6 }$$

where ${m}_{{\rm{mol}},j}$, ${V}_{{\rm{mol}},j}$, and ${\gamma }_{j}$ are the molecular mass, molecular volume, and surface energy of the species, respectively. The small prefactor $4.7\times {10}^{-3}$ comes from the fact that the neck radius is much smaller than a₀. In general, ${t}_{{\rm{sint}},j}$ rapidly decreases with increasing T because ${P}_{{\rm{ev}},j}$ strongly depends on T. For species other than ${{\rm{H}}}_{2}{\rm{S}}$ and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$, we use the same set of ${V}_{{\rm{mol}},j}$ and ${\gamma }_{j}$ adopted by Sirono (2011b). The molecular volumes of ${{\rm{H}}}_{2}{\rm{S}}$ and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ are estimated as ${V}_{\mathrm{mol},{{\rm{H}}}_{2}{\rm{S}}}=5.7\times {10}^{-23}\;{\mathrm{cm}}^{3}$ and ${V}_{\mathrm{mol},{{\rm{C}}}_{2}{{\rm{H}}}_{6}}=7.1\times {10}^{-23}\;{\mathrm{cm}}^{3}$ assuming that the densities of ${{\rm{H}}}_{2}{\rm{S}}$ and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ solids are $1\;{\rm{g}}\;{\mathrm{cm}}^{-3}$ and $0.7\;{\rm{g}}\;{\mathrm{cm}}^{-3}$ (Moses et al. 1992), respectively. We assume that the surface energy of ${{\rm{H}}}_{2}{\rm{S}}$ ice is equal to that of ${{\rm{H}}}_{2}{\rm{S}}$ liquid: $30\;\mathrm{erg}\;{\mathrm{cm}}^{-2}$ (Meyer 1977). For ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$, we use ${\gamma }_{{{\rm{C}}}_{2}{{\rm{H}}}_{6}}=40\;\mathrm{erg}\;{\mathrm{cm}}^{-2}$, which is the value at $\sim 50\;{\rm{K}}$ (Moses et al. 1992). The locations of the sintering lines discussed below are insensitive to the values of ${V}_{{\rm{mol}},j}$ and ${\gamma }_{j}$ because of the strong temperature dependence of ${P}_{{\rm{ev}},j}$.

The necks not only grow but also are destroyed by at least two processes.

• 1.  
  The necks evaporate when the ambient gas temperature exceeds the sublimation temperature of the volatile j that constitutes the necks. This occurs at $r\lt {r}_{{\rm{snow}},j}$, where ${r}_{{\rm{snow}},j}$ is the orbital radius of the snow line of the volatile.
• 2.  
  The necks break when the aggregate is plastically deformed by another aggregate upon collision. Unsintered aggregates are known to experience substantial plastic deformation even if the collision velocity is much below the fragmentation threshold (Dominik & Tielens 1997; Wada et al. 2008). Therefore, fully sintered aggregates form only if they do not collide with each other until the sintering is completed, i.e., only if the sintering timescale is shorter than their collision timescale, which we will denote by ${t}_{{\rm{coll}}}$. Because ${t}_{{\rm{sint}},j}$ falls off rapidly toward the central star, there exists a location $r={r}_{{\rm{sint}},j}$ inside which ${t}_{{\rm{sint}},j}\lt {t}_{{\rm{coll}}}$ for a given volatile species j. In this study, we call these locations the sintering lines.

Taken together, each volatile species j causes aggregate sintering only inside the annulus defined by ${r}_{{\rm{snow}},j}\lt r\lt {r}_{{\rm{sint}},j}$. We call such regions the sintering zones. Our sintering zones are essentially equivalent to the "sintering regions" of Sirono (2011b).

In order to find the location of the sintering zones, one needs to estimate ${t}_{{\rm{coll}}}$. In principle, the collision time can be calculated using the size, number density, and relative speed of aggregates, as we do in our simulations (see Equation (12)). In this subsection, we avoid such detailed calculations and instead use the formula ${t}_{{\rm{coll}}}=100{{\rm{\Omega }}}^{-1}$, which is an approximate expression for the collision timescale of macroscopic aggregates in a turbulent disk with the dust-to-gas mass ratio of 0.01 (Takeuchi & Lin 2005; Brauer et al. 2008). This is a rough estimate (see Sato et al. 2015), but still provides a reasonable estimate for ${r}_{{\rm{sint}},j}$ because ${t}_{{\rm{sint}},j}$ is a steep function of r.

In Figure 4, we plot the sintering timescales of the seven volatile species as a function of r for our HL Tau disk model with $\gamma =1$. We also indicate ${t}_{{\rm{coll}}}=100{{\rm{\Omega }}}^{-1}$ by the dotted lines, the location of the sintering lines (${r}_{{\rm{sint}},j}$) by the open circles, and the locations of the sintering zones (${r}_{{\rm{snow}},j}\lt r\lt {r}_{{\rm{sint}},j}$) by the horizontal bars. We consider two different values of a₀—0.1 and $4\;\mu {\rm{m}}$ (panels (a) and (b), respectively)—to highlight the importance of this parameter in our sintering model. For ${a}_{0}=0.1\;\mu {\rm{m}}$, the width of each sintering zone is 30%–50% of ${r}_{{\rm{sint}},j}$, and the sintering zones of ${\mathrm{NH}}_{3}$, ${\mathrm{CO}}_{2}$, and ${{\rm{H}}}_{2}{\rm{S}}$ significantly overlap with each other. A larger a₀ leads to longer sintering timescales and hence to narrower sintering zones. For ${a}_{0}=4\;\mu {\rm{m}}$, the sintering zones for species other than ${{\rm{H}}}_{2}{\rm{O}}$ almost disappear. For this reason, we restrict ourselves to ${a}_{0}\leqslant 4\;\mu {\rm{m}}$ in the following sections.

We simulate the global evolution of particles in a gas disk using the single-size approximation. We assume that the total solid mass at each orbital radius r is dominated by particles with a mass $m={m}_{*}(r)$. We then follow the evolution of the solid surface density ${{\rm{\Sigma }}}_{{\rm{d}}}(r)$ and "representative" (mass-dominating) particle mass ${m}_{*}(r)$, taking into account aggregate collision and radial drift (see Equations (7) and (8) below). The single-size approach (or mathematically speaking, moment approach) has often been used in the modeling of particle growth in planetary atmospheres (e.g., Ferrier 1994; Ormel 2014) as well as in protoplanetary disks (e.g., Kornet et al. 2001; Garaud 2007; Birnstiel et al. 2012; Estrada et al. 2015; Krijt et al. 2016, see also Appendix of Sato et al. 2015 for the mathematical background of the single-size approximation and a comparison between single-size and full-size simulations). This approach allows us to track the evolution of the mass budget of dust in a disk without using the computationally expensive Smoluchowski's coagulation equation. A drawback of this approach is that one has to assume the aggregate size distribution at each orbit whenever it is needed. In this study, we assume a power-law size distribution when we predict dust emission from the disk (see Section 4.5).

Under the single-size approximation, the evolution of ${{\rm{\Sigma }}}_{{\rm{d}}}$ and m_* is described by (Ormel 2014; Sato et al. 2015)

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}({{rv}}_{{\rm{r}}}{{\rm{\Sigma }}}_{{\rm{d}}})=0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {m}_{*}}{\partial t}+{v}_{{\rm{r}}}\displaystyle \frac{\partial {m}_{*}}{\partial r}=\displaystyle \frac{{\rm{\Delta }}{m}_{*}}{{t}_{{\rm{coll}}}},\end{eqnarray} \tag{ 8 }$$

where ${v}_{{\rm{r}}}$ and ${t}_{{\rm{coll}}}$ are the radial drift velocity and mean collision time of the representative particles, respectively, and ${\rm{\Delta }}{m}_{*}$ is the change of m_* upon a single aggregate collision. Equation (7) expresses the mass conservation for solids, while Equation (8) states that the growth rate of representative particles along their trajectory, ${{Dm}}_{*}/{Dt}\equiv \partial {m}_{*}/\partial t\;+{v}_{{\rm{r}}}\partial {m}_{*}/\partial r$, is equal to ${\rm{\Delta }}{m}_{*}/{t}_{{\rm{coll}}}$. The mass change ${\rm{\Delta }}{m}_{*}$ is equal to m_* when the collision results in pure sticking, while ${\rm{\Delta }}{m}_{*}\lt {m}_{*}$ when fragmentation or erosion occurs. Our Equations (7) and (8) correspond to Equations (3) and (8) of Ormel (2014), respectively, although Equation (8) of Ormel (2014) assumes ${\rm{\Delta }}{m}_{*}={m}_{*}$. The expressions for ${v}_{{\rm{r}}}$, ${t}_{{\rm{coll}}}$, and ${\rm{\Delta }}{m}_{*}$ will be given in Sections 4.2–4.4, respectively.

We solve Equations (7) and (8) by discretizing the radial direction into 300 logarithmically spaced binds spanning from ${r}_{{\rm{in}}}=1\;\mathrm{au}$ to ${r}_{{\rm{out}}}=1000\;\mathrm{au}$. The parameter sets and initial conditions used in the simulations are described in Section 4.7.

The motion of a particle in a gas disk is characterized by the dimensionless Stokes number ${\rm{St}}\equiv {\rm{\Omega }}{t}_{{\rm{s}}}$, where Ω is the Keplerian frequency and ${t}_{{\rm{s}}}$ is the particle's stopping time. When the particle radius is much smaller than the mean free path of the gas molecules in the disk, as is true for the particles in our simulations, the stopping time is given by Epstein's drag law. The Stokes number of a representative aggregate at the midplane can then be written as (Birnstiel et al. 2010)

$$\begin{eqnarray}&&{\rm{St}}=\displaystyle \frac{\pi }{2}\displaystyle \frac{{\rho }_{{\rm{int}}}{a}_{*}}{{{\rm{\Sigma }}}_{{\rm{g}}}},\end{eqnarray} \tag{ 9 }$$

where ${\rho }_{{\rm{int}}}=0.26\;{\rm{g}}\;{\mathrm{cm}}^{-3}$ and ${a}_{*}\equiv {(3{m}_{*}/4\pi {\rho }_{{\rm{int}}})}^{1/3}$ are the internal density and radius of the aggregate, respectively. We use Equation (9) whenever we calculate ${\rm{St}}$.

The drift velocity ${v}_{{\rm{r}}}$ is given by (Adachi et al. 1976; Weidenschilling 1977)

$$\begin{eqnarray}&&{v}_{{\rm{r}}}=-2\eta {v}_{{\rm{K}}}\displaystyle \frac{{\rm{St}}}{1+{\mathrm{St}}^{2}},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&\eta =-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{c}_{{\rm{s}}}}{{v}_{{\rm{K}}}}\right)}^{2}\displaystyle \frac{d\mathrm{ln}{P}_{{\rm{g}}}}{d\mathrm{ln}r}\end{eqnarray} \tag{ 11 }$$

is the parameter characterizing the sub-Keplerian motion of the gas disk, ${v}_{{\rm{K}}}=r{\rm{\Omega }}$ is the Keplerian velocity, and ${P}_{{\rm{g}}}={\rho }_{{\rm{g}}}{c}_{{\rm{s}}}^{2}$ is the midplane gas pressure. In our disk model, ${{dP}}_{{\rm{g}}}/{dr}\lt 0$ and therefore ${v}_{{\rm{r}}}\lt 0$ everywhere. At $r\ll {r}_{{\rm{c}}}$, we approximately have $\eta \approx 1.8$ × ${10}^{-3}{(r/1{\rm{au}})}^{0.43}$ and $\eta {v}_{{\rm{K}}}\approx 52{(r/1{\rm{au}})}^{-0.07}\quad {\rm{m}}\;{{\rm{s}}}^{-1}$.

We evaluate the collision term ${\rm{\Delta }}{m}_{*}/{t}_{{\rm{coll}}}$ assuming that collisions between representative aggregates dominate the evolution of m_*. Under this assumption, the collision time is approximately given by

$$\begin{eqnarray}&&{t}_{{\rm{coll}}}=\displaystyle \frac{1}{4\pi {a}_{*}^{2}{n}_{*}{\rm{\Delta }}v},\end{eqnarray} \tag{ 12 }$$

where $4\pi {a}_{*}^{2}$, n_*, and ${\rm{\Delta }}v$ are the collisional cross section, number density, and collision velocity of the representative aggregates, respectively. We use the midplane values for n_* and ${\rm{\Delta }}v$. We do not consider the erosion of representative aggregates by a number of small grains (Seizinger et al. 2013; Krijt et al. 2015) because there remain uncertainties in the threshold velocity for erosive collisions as a function of the projectile mass (see Section 2.3.2 of Krijt et al. 2015).

To evaluate n_*, we consider disk turbulence and assume that vertical settling balances with turbulent diffusion for the representative aggregates. We parameterize the strength of disk turbulence with the dimensionless parameter ${\alpha }_{{\rm{t}}}=D/{c}_{{\rm{s}}}{H}_{{\rm{g}}}$, where D is the particle diffusion coefficient in the turbulence. For simplicity, ${\alpha }_{{\rm{t}}}$ is assumed to be independent of time and distance from the midplane, but we allow ${\alpha }_{{\rm{t}}}$ to depend on r (see Section 4.7). Under this assumption, n_* at the midplane can be written as

$$\begin{eqnarray}&&{n}_{*}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{\sqrt{2\pi }{H}_{{\rm{d}}}{m}_{*}},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&{H}_{{\rm{d}}}={\left(1+\displaystyle \frac{{\rm{St}}}{{\alpha }_{{\rm{t}}}}\displaystyle \frac{1+2{\rm{St}}}{1+{\rm{St}}}\right)}^{-1/2}{H}_{{\rm{g}}}.\end{eqnarray} \tag{ 14 }$$

is the scale height of the representative aggregates (Dubrulle et al. 1995; Youdin & Lithwick 2007).

The collision velocity ${\rm{\Delta }}v$ is given by the root sum square of the contributions from Brownian motion, gas turbulence (Ormel & Cuzzi 2007), and size-dependent drift relative to the gas disk (Adachi et al. 1976; Weidenschilling 1977). The expressions for these contributions can be found in, for example, Section 2.3.2 of Okuzumi et al. (2012). The contributions from turbulence and drift motion are functions of the Stokes numbers of two colliding aggregates, ${{\rm{St}}}_{1}$ and ${{\rm{St}}}_{2}$. In this study, we set ${{\rm{St}}}_{1}={\rm{St}}$ and ${{\rm{St}}}_{2}=0.5{\rm{St}}$ because we consider collisions between aggregates similar in size. Sato et al. (2015) and Krijt et al. (2016) have recently shown that such a choice best reproduces the results of coagulation simulations that resolve the full size distribution of the aggregates. As long as ${\rm{St}}\ll 1$, the collision velocity is an increasing function of ${\rm{St}}$.

For macroscopic aggregates satisfying ${10}^{-4}\lesssim {\rm{St}}\ll 1$, either turbulence or radial drift mostly dominates their collision velocity. For this range of ${\rm{St}}$, the collision velocities driven by turbulence and radial drift are approximately given by ${\rm{\Delta }}{v}_{{\rm{t}}}\approx \sqrt{2.3{\alpha }_{{\rm{t}}}{\rm{St}}}{c}_{{\rm{s}}}$ and ${\rm{\Delta }}{v}_{{\rm{r}}}\approx 2\eta {v}_{{\rm{K}}}| {{\rm{St}}}_{1}-{{\rm{St}}}_{2}| \approx \eta {v}_{{\rm{K}}}{\rm{St}}$, respectively, where we use ${{\rm{St}}}_{1}={\rm{St}}$ and ${{\rm{St}}}_{2}=0.5{\rm{St}}$ (for the expression of ${\rm{\Delta }}{v}_{{\rm{t}}}$, see Equation (28) of Ormel & Cuzzi 2007).

We denote the change in m_* due to a single collision between two mass-dominating aggregates by ${\rm{\Delta }}{m}_{*}$. Following Okuzumi & Hirose (2012), we model ${\rm{\Delta }}{m}_{*}$ as

$$\begin{eqnarray}&&{\rm{\Delta }}{m}_{*}=\mathrm{min}\left\{1,-\displaystyle \frac{\mathrm{ln}({\rm{\Delta }}v/{\rm{\Delta }}{v}_{{\rm{frag}}})}{\mathrm{ln}5}\right\}{m}_{*},\end{eqnarray} \tag{ 15 }$$

where the fragmentation threshold ${\rm{\Delta }}{v}_{{\rm{frag}}}$ characterizes the sticking efficiency of the colliding aggregates. The mass change is positive for ${\rm{\Delta }}v\lt {\rm{\Delta }}{v}_{{\rm{frag}}}$ and negative for ${\rm{\Delta }}v\gt {\rm{\Delta }}{v}_{{\rm{frag}}}$ as shown in Figure 5(a). Equation (15) is a fit to the data of the collision simulations for unsintered aggregates by Wada et al. (2009, their Figure 11). For sintered aggregates, Equation (15) overestimates the sticking efficiency at low collision velocities where the aggregates bounce rather than stick (Sirono 1999). However, we show in Section 7.1 that the bouncing hardly alters the evolution of m_* in the sintering zone. For simplicity, we use Equation (15) for both unsintered and sintered aggregates.

Zoom In Zoom Out Reset image size

To account for the effect of sintering on the aggregate sticking efficiency, we model the fragmentation threshold ${\rm{\Delta }}{v}_{{\rm{frag}}}$ as

$$\begin{eqnarray}{\rm{\Delta }}{v}_{{\rm{frag}}} & = & {\rm{\Delta }}{v}_{{\rm{frag,NS}}}\mathrm{exp}\left(-\displaystyle \frac{{t}_{{\rm{coll}}}}{{t}_{{\rm{sint,eff}}}}\right)\\ & & +{\rm{\Delta }}{v}_{{\rm{frag,S}}}\left[1-\mathrm{exp}\left(-\displaystyle \frac{{t}_{{\rm{coll}}}}{{t}_{{\rm{sint,eff}}}}\right)\right],\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{t}_{{\rm{sint,eff}}}\equiv {\left(\sum _{j}{t}_{{\rm{sint}},j}^{-1}\right)}^{-1}\end{eqnarray} \tag{ 17 }$$

is the effective sintering timescale (see Equation (6) for the definition of the individual sintering timescale ${t}_{{\rm{sint}},j}$) and ${\rm{\Delta }}{v}_{{\rm{frag,NS}}}$ and ${\rm{\Delta }}{v}_{{\rm{frag,S}}}(\lt {\rm{\Delta }}{v}_{{\rm{frag,NS}}})$ are the thresholds for unsintered and sintered aggregates, respectively. In Equation (17), the summation is taken over all solid-phase volatiles (i.e., volatiles satisfying $r\gt {r}_{{\rm{snow}},j}$). As described in Section 3.3, the snow line location ${r}_{{\rm{snow}},j}$ for each species is calculated from the relation ${P}_{{\rm{ev}},j}({r}_{{\rm{snow}},j})={P}_{j}({r}_{{\rm{snow}},j})$ with Equations (4) and (5). At $r\lt {r}_{\mathrm{snow},{{\rm{H}}}_{2}{\rm{O}}}$, where all volatile ices sublimate, we exceptionally set ${\rm{\Delta }}{v}_{{\rm{frag}}}={\rm{\Delta }}{v}_{{\rm{frag,S}}}$ to mimic the low sticking efficiency of bare silicate grains compared to (unsintered) ice-coated grains (e.g., Chokshi et al. 1993).

Equation (16) is constructed so that ${\rm{\Delta }}{v}_{{\rm{frag}}}\approx {\rm{\Delta }}{v}_{{\rm{frag,NS}}}$ in the non-sintering zones (${t}_{{\rm{coll}}}\ll {t}_{{\rm{sint,eff}}}$) and ${\rm{\Delta }}{v}_{{\rm{frag}}}\approx {\rm{\Delta }}{v}_{{\rm{frag,S}}}$ in the sintering zone (${t}_{{\rm{coll}}}\gg {t}_{{\rm{sint,eff}}}$). Simulations of aggregate collisions suggest that ${\rm{\Delta }}{v}_{{\rm{frag,NS}}}\sim 50\;{\rm{m}}\;{{\rm{s}}}^{-1}$ (Wada et al. 2009) and ${\rm{\Delta }}{v}_{{\rm{frag,S}}}\sim 20\;{\rm{m}}\;{{\rm{s}}}^{-1}$ (Sirono & Ueno 2014) if the colliding aggregates are identical and made of $0.1\;\mu {\rm{m}}$ sized ice monomers. The theory of particle sticking (Johnson et al. 1971), on which the simulations by Wada et al. (2009) are based, indicates that ${\rm{\Delta }}{v}_{{\rm{frag,NS}}}$ scales with the monomer size as ${a}_{0}^{-5/6}$ (Chokshi et al. 1993; Dominik & Tielens 1997). For ${\rm{\Delta }}{v}_{{\rm{frag,S}}}$, the scaling is yet to be studied, so we simply assume the same scaling as for ${\rm{\Delta }}{v}_{{\rm{frag,NS}}}$. We thus model ${\rm{\Delta }}{v}_{{\rm{frag,NS}}}$ and ${\rm{\Delta }}{v}_{{\rm{frag,S}}}$ as

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{frag,NS}}}=50{\left(\displaystyle \frac{{a}_{0}}{0.1\mu {\rm{m}}}\right)}^{-5/6}\;{\rm{m}}\;{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{frag,S}}}=20{\left(\displaystyle \frac{{a}_{0}}{0.1\mu {\rm{m}}}\right)}^{-5/6}\;{\rm{m}}\;{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 19 }$$

Figure 5(b) shows ${\rm{\Delta }}{v}_{{\rm{frag}}}$ versus ${t}_{{\rm{coll}}}/{t}_{{\rm{sint,eff}}}$ for ${a}_{0}=0.1\;\mu {\rm{m}}$.

We calculate the absorption cross section of porous aggregates using the analytic expression by Kataoka et al. (2014, their Equation (18)), which is based on Mie calculations with effective medium theory. Monomers are treated as composite spherical grains made of astronomical silicates, carbonaceous materials, and water ice with the mass abundance ratio of 2.64:3.53:5.55 (Pollack et al. 1994). We calculate the effective refractive index of the monomers using the Bruggeman mixing rule. The optical constants of silicates, carbons, and water ice are taken from Draine (2003), Zubko et al. (1996, data for ACH2 samples), and Warren (1984, data for $T=-60^\circ {\rm{C}}$), respectively. We neglect the contribution of volatiles other than ${{\rm{H}}}_{2}{\rm{O}}$ to the monomer optical properties. The effective refractive index of porous aggregates is computed using the Maxwell–Garnett rule in which the monomers are regarded as inclusions in vacuum.

Because we adopt the single-size approach, we only track the evolution of aggregates dominating the dust surface density. However, these aggregates do not necessarily dominate millimeter dust emission from a disk. While the mass-dominating aggregates are generally the largest aggregates in the population (e.g., Birnstiel et al. 2012; Okuzumi et al. 2012), smaller ones can dominate the millimeter opacity of the population when the largest aggregates are significantly larger than a millimeter in radius. In order to take into account this effect, we only assume a size distribution when we calculate dust opacities. Specifically, we assume a power-law distribution

$$\begin{eqnarray}{N}_{{\rm{d}}}^{\prime }(a)=\left\{\begin{array}{ll}{{Ca}}^{-3.5}, & {a}_{{\rm{min}}}\lt a\lt {a}_{*}\\ 0, & {\rm{otherwise}},\end{array}\right.\end{eqnarray} \tag{ 20 }$$

where ${N}_{{\rm{d}}}^{\prime }(a)$ is the column number density of aggregates per unit aggregate radius a ($\lt {a}_{*}$), ${a}_{{\rm{min}}}$ is the minimum aggregate radius, and C is the normalization constant determined by the condition ${\int }_{0}^{\infty }{{mN}}_{{\rm{d}}}^{\prime }(a){\rm{d}}a={{\rm{\Sigma }}}_{{\rm{d}}}$. We fix ${a}_{{\rm{min}}}$ to be $0.1\;\mu {\rm{m}}$ with the understanding that millimeter opacities are insensitive to the choice of ${a}_{{\rm{min}}}$ as long as ${a}_{{\rm{min}}}\ll 1\;{\rm{mm}}$. The slope of −3.5 is based on the classical theory of fragmentation cascades (Dohnanyi 1969; Tanaka et al. 1996; see Birnstiel et al. 2011 for how the coagulation of the fragments modifies this value). Therefore, Equation (20) would overestimate the amount of fragments when the collisions between the largest aggregates (which are the source of the fragments) do not lead to their catastrophic disruption. Possible effects of this simplification are discussed in Section 7.2.

The upper panel of Figure 6 shows our dust opacities ${\kappa }_{{\rm{d}},\nu }$ at wavelengths $\lambda =0.87$, 1.3, and 2.9 mm (corresponding to ALMA Bands 7, 6, and 3, respectively) as a function of a_*. Note that the opacities are expressed in units of ${\mathrm{cm}}^{2}$ per gram of dust. In the lower panel of Figure 6, we plot the opacity slope measured at $\lambda =0.87$–1.3 mm, ${\beta }_{0.87-1.3\mathrm{mm}}\quad \equiv \mathrm{ln}({\kappa }_{\text{0.87 mm}}/{\kappa }_{\text{1.3 mm}})/\mathrm{ln}({\nu }_{\text{0.87 mm}}/{\nu }_{\text{1.3 mm}})$. Our model gives ${\beta }_{0.87-1.3\mathrm{mm}}\approx 1.7$ at ${a}_{*}\ll 1\;{\rm{cm}}$ and ${\beta }_{0.87-1.3\mathrm{mm}}\approx 0.8$ in the opposite limit.

Zoom In Zoom Out Reset image size

We calculate the intensities I_ν of dust thermal emission at each orbital radius r as

$$\begin{eqnarray}&&{I}_{\nu }(r)=[1-\mathrm{exp}(-{\tau }_{\nu }(r))]{B}_{\nu }(T(r)),\end{eqnarray} \tag{ 21 }$$

where ${B}_{\nu }(T)$ is the Planck function,

$$\begin{eqnarray}&&{\tau }_{\nu }(r)=\displaystyle \frac{{\kappa }_{{\rm{d}},\nu }(r){{\rm{\Sigma }}}_{{\rm{d}}}(r)}{\mathrm{cos}i}\end{eqnarray} \tag{ 22 }$$

is the line of sight optical depth, and i is the disk inclination. We use $i=46\_\_AMP\_\_fdg;7$ for the HL Tau disk (ALMA Partnership et al. 2015). The Planck brightness temperature ${T}_{{\rm{B}}}$ is computed by solving the equation ${I}_{\nu }={B}_{\nu }({T}_{{\rm{B}}})$ for ${T}_{{\rm{B}}}$. Equation (22) assumes that the dust disk is geometrically thin (i.e., the radial distance over which ${\kappa }_{{\rm{d}},\nu }(r)$ and ${{\rm{\Sigma }}}_{{\rm{d}}}(r)$ vary is longer than the dust scale height). This assumption is, however, not always satisfied in our simulations, as we discuss in detail in Section 6.3.

We will also use the flux density

$$\begin{eqnarray}&&{F}_{\nu }=\displaystyle \frac{2\pi \mathrm{cos}i}{{d}^{2}}{\int }_{{r}_{{\rm{in}}}}^{{r}_{{\rm{out}}}}{I}_{\nu }(r){rdr},\end{eqnarray} \tag{ 23 }$$

where d is the distance to HL Tau, ${r}_{{\rm{in}}}=1\;\mathrm{au}$ and ${r}_{{\rm{out}}}=1000\;\mathrm{au}$ are the boundaries of our computational domain (see Section 4.1), and the factor $\mathrm{cos}i$ accounts for the ellipticity of the disk image. In accordance with ALMA Partnership et al. (2015), we set $d=140\;{\rm{pc}}$, the standard mean distance to Taurus.

When comparing our simulation results with the ALMA observation, it is useful to smooth the simulated radial emission profiles at the spatial resolution of ALMA. In this study, we do this in two steps. First, we generate projected images of our simulation snapshots assuming the disk inclination of $46\_\_AMP\_\_fdg;7$. For simplicity, the geometrical thickness of the disks is neglected in this process. Second, we smooth the "raw" images along their major axis using a circular Gaussian with the FWHM angular resolutions of $\sqrt{85.3\times 61.1}$, $\sqrt{35.1\times 21.8}$, and $\sqrt{29.9\times 19.0}$ mas for $\lambda =2.9$, 1.3, and 0.87 mm, in accordance with the ALMA observation of HL Tau at Bands 3, 6, and 7, respectively (ALMA Partnership et al. 2015). Because we assume $d=140\;{\rm{}}\;{\rm{pc}}$, these angular resolutions translate into the spatial resolutions of ≈10, 3.9, and 3.3 au at Bands 3, 6, and 7, respectively.

We conduct 10 simulation runs with different sets of model parameters. Columns 1 through 4 of Table 3 list the run names and parameter choices (γ, a₀, ${\alpha }_{{\rm{t}}}$) for the simulation runs. Run Sa0 is our fiducial model and assumes $\gamma =1$, ${a}_{0}=0.1\;\mu {\rm{m}}$, and ${\alpha }_{{\rm{t}}}=0.03{(r/10\mathrm{au})}^{1/2}$. Model Sa0-NoSint is the same as model Sa0 but neglects sintering. Runs Sa0-Lgam and Sa0-Hgam are designed to study the dependence of the results on the gas surface density slope γ. Runs Sa0-Lalp and Sa0-Halp will be used to study how the sintering-induced ring formation scenario constrains the radial distribution of ${\alpha }_{{\rm{t}}}$ in the HL Tau disk. Runs La0 and LLa0 assume more fragile aggregates (i.e., larger a₀) and weaker turbulence than in the fiducial run. As we will see, the set of a₀ and ${\alpha }_{{\rm{t}}}$ controls the degree of dust sedimentation (i.e., the geometrical thickness of the dust disk). Runs Sa0-tuned and La0-tuned are the same as runs Sa0 and La0, respectively, except that they adopt slightly lower L_j for ${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{NH}}_{3}$, and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ (see ${L}_{j,\mathrm{tuned}}$ in Table 2). These runs will be used to quantify possible uncertainties of our results that might arise from the uncertainties in the vapor pressure data.

Table 3.  List of Simulation Runs

Run	γ	a0	${\alpha }_{{\rm{t}}}$	${t}_{{\rm{snap}}}$	Fν (Jy) at $t={t}_{{\rm{snap}}}$			Section
		(μm)		(Myr)	2.9 mm	1.3 mm	0.87 mm
Sa0	1	0.1	$0.03{(r/10\mathrm{au})}^{1/2}$	0.26	0.070	0.79	2.2	5
Sa0-NoSinta	1	0.1	$0.03{(r/10\mathrm{au})}^{1/2}$	0.12	0.063	0.79	2.3	5.5
Sa0-Lgam	0.5	0.1	$0.03{(r/10\mathrm{au})}^{1/2}$	0.29	0.076	0.76	2.1	6.1
Sa0-Hgam	1.5	0.1	$0.03{(r/10\mathrm{au})}^{1/2}$	0.05	0.072	0.80	2.3	6.1
Sa0-Lalp	1	0.1	0.03	0.18	0.066	0.77	2.2	6.2
Sa0-Halp	1	0.1	0.1	0.29	0.075	0.75	2.1	6.2
La0	1	1	${10}^{-3}{(r/10\mathrm{au})}^{1/2}$	0.41	0.064	0.78	2.3	6.3
LLa0	1	4	${10}^{-4}{(r/10\mathrm{au})}^{1/2}$	0.45	0.062	0.80	2.4	6.3
Sa0-tunedb	1	0.1	$0.03{(r/10\mathrm{au})}^{1/2}$	0.27	0.068	0.77	2.2	6.4
La0-tunedb	1	1	${10}^{-3}{(r/10\mathrm{au})}^{1/2}$	0.41	0.063	0.79	2.3	6.4

Notes.

^aNo sintering. ^bUses ${L}_{j,{\rm{tuned}}}$ instead of L_j for ${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{NH}}_{3}$, and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ (see Table 2).

Download table as:  ASCIITypeset image

The initial conditions are given by ${{\rm{\Sigma }}}_{{\rm{d}}}(t=0,r)\;=0.01{{\rm{\Sigma }}}_{{\rm{g}}}(r)$ and ${a}_{*}(t=0,r)={a}_{0}$, where we have assumed that the dust-to-gas mass ratio of the initial disk is 0.01.

In our simulations, the observational appearance of the disk changes with time because dust particles grow and drift inward. In particular, the millimeter emission of the disk diminishes as the particles drain onto the central star. For this reason, we select from each simulation run one snapshot that best reproduces the millimeter flux densities of the HL Tau disk reported by ALMA Partnership et al. (2015). Specifically, we calculate the relative errors between the simulated and observed flux densities at $\lambda =0.87$, 1.3, and 2.9 mm as a function of time, and search for the time $t={t}_{{\rm{snap}}}$ at which the sum of the relative errors is minimized. For reference, the flux densities reported by the ALMA observation are 0.0743, 0.744, and 2.14 Jy at $\lambda =0.87$, 1.3, and 2.9 mm (Bands 7, 6, and 3), respectively (ALMA Partnership et al. 2015).

Figure 7 illustrates such an analysis for our fiducial simulation run Sa0. This figure shows the simulated time evolution of the flux densities at the three wavelengths as well as the evolution of the total dust mass ${M}_{{\rm{d}}}$ within the computational domain. The dust mass and flux densities decrease on a timescale of $\sim 1\;\mathrm{Myr}$, which reflects the timescale on which the dust near the disk's outer edge (which dominates ${M}_{{\rm{d}}}$) grows into rapidly drifting pebbles (see, e.g., Sato et al. 2015). Comparing the flux densities from the simulation with those from the ALMA observations (shown by the dashed horizontal line segments in the lower panel of Figure 7), we find that the sum of the relative errors in the flux densities is minimized when $t=0.26\;\mathrm{Myr}$. At this time, the flux densities in the simulation are 0.070, 0.79, and 2.2 Jy at $\lambda \;=$ 2.9, 1.3, and 0.87 mm, respectively, which is in agreement with the ALMA measurements to an accuracy of less than 6%.

Zoom In Zoom Out Reset image size

Columns 5 through 8 of Table 3 list the values of ${t}_{{\rm{snap}}}$ and ${F}_{\nu }(t={t}_{{\rm{snap}}})$ for all simulation runs. We find that ${t}_{{\rm{snap}}}$ falls within the range 0.1–$0.5\;\mathrm{Myr}$. Since ${t}_{{\rm{snap}}}$ may be regarded as the time after disk formation, our results are consistent with the idea that HL Tau is younger than 1 Myr. In fact, the age predicted from our simulations depends on the disk mass ${M}_{{\rm{disk}}}$ assumed: a higher ${M}_{{\rm{disk}}}$ leads to a larger ${t}_{{\rm{snap}}}$ because it takes longer for dust emission to decay to the observed level when the initial dust mass is larger. However, a disk mass much in excess of $0.2\;{M}_{\odot }$ seems to be unrealistic because the disk would then be gravitationally unstable at outer radii (see Section 2.2).

Below we fix t to be ${t}_{{\rm{snap}}}=0.26\;\mathrm{Myr}$ and look at the radial distribution of dust in detail. Figures 8(a) and (b) show the radial distribution of the representative aggregate radius a_* and dust surface density ${{\rm{\Sigma }}}_{{\rm{d}}}$ at this time. We also show in Figure 8(c) the Stokes number ${\rm{St}}$ of the representative aggregates, which is more directly related to their dynamics than a_*. In these figures, the vertical stripes indicate the locations of the sintering zones. Here, the sintering zone of each volatile j is defined by the locations where $r\gt {r}_{{\rm{snow}},j}$ and $r\lt {r}_{{\rm{sint}},j}$, with the latter being equivalent to ${t}_{{\rm{coll}}}\lt {t}_{{\rm{sint}},j}$ (see Figure 8(e) for the radial distribution of ${t}_{{\rm{coll}}}$ and ${t}_{{\rm{sint,eff}}}$). The sintering zones of ${\mathrm{NH}}_{3}$, ${\mathrm{CO}}_{2}$, and ${{\rm{H}}}_{2}{\rm{S}}$ partially overlap with each other and form a single sintering zone. The exact locations of the sintering zones are 3–6 au (${{\rm{H}}}_{2}{\rm{O}}$), 11–23 au (${\mathrm{NH}}_{3}$–${\mathrm{CO}}_{2}$–${{\rm{H}}}_{2}{\rm{S}}$), 24–33 au (${{\rm{C}}}_{2}{{\rm{H}}}_{6}$), 80–106 au (${\mathrm{CH}}_{4}$), and 116–160 au (${\rm{CO}}$). Strictly speaking, the locations of the sintering zones are time-dependent because the volatile partial pressures P_j and aggregate collision timescale ${t}_{{\rm{coll}}}$ evolve with ${{\rm{\Sigma }}}_{{\rm{d}}}$ (see Equations (4), (5), and (12)). However, a comparison of Figures 4 and 8 shows that the sintering zones little migrate during this $0.26\;\mathrm{Myr}$. This is because the locations of the sintering zones depend on the radial distribution of the gas temperature T (which is taken to be time-independent) much more strongly than on the distribution of ${{\rm{\Sigma }}}_{{\rm{g}}}$.

Zoom In Zoom Out Reset image size

Figures 8(a) and (b) show that sintering produces a clear pattern in the radial distribution of the dust component. We see that dust aggregates in the sintering zones tend to have a high surface density and a small radius compared to those in the adjacent non-sintering zones. The small aggregate size is a direct consequence of fragmentation induced by sintering. To see this, we plot in Figure 8(f) the collision velocity ${\rm{\Delta }}v$ and fragmentation threshold ${\rm{\Delta }}{v}_{{\rm{frag}}}$ as a function of r. In the non-sintering zones, we find ${\rm{\Delta }}{v}_{{\rm{frag}}}\approx 50\;{\rm{m}}\;{{\rm{s}}}^{-1}$ and ${\rm{\Delta }}v\;\approx 25$–$35\;{\rm{m}}\;{{\rm{s}}}^{-1}$, implying that no disruptive collisions occur for the unsintered aggregates (as we show below, the maximum size of the unsintered aggregates is determined by radial drift rather than fragmentation). In the sintering zones, ${\rm{\Delta }}{v}_{{\rm{frag}}}$ is decreased to $20\;{\rm{m}}\;{{\rm{s}}}^{-1}$ and ${\rm{\Delta }}v$ is suppressed down to the same value. Since ${\rm{\Delta }}v$ is an increasing function of a_* (as long as ${\rm{St}}\ll 1$), this indicates that the sintered aggregates disrupt each other so that ${\rm{\Delta }}v$ never exceeds ${\rm{\Delta }}{v}_{{\rm{frag}}}$. The disrupted aggregates pile up there because the inward drift speed $| {v}_{{\rm{r}}}|$ decreases with decreasing a_*. These pileups provide the high surface densities in the sintering zones.

To understand the radial distribution of ${{\rm{\Sigma }}}_{{\rm{d}}}$ and a_* more quantitatively, we look at the radial inward mass flux of drifting aggregates,

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{d}}}\equiv 2\pi r| {v}_{{\rm{r}}}| {{\rm{\Sigma }}}_{{\rm{d}}}.\end{eqnarray} \tag{ 24 }$$

The radial distribution of ${\dot{M}}_{{\rm{d}}}$ is shown in Figure 8(d). The mass flux is radially constant ($\approx 5\times {10}^{-9}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}$) at $\lesssim 100\;\mathrm{au}$, which indicates that the radial dust flow in this region can be approximated by a steady flow. Such a quasi-steady dust flow is commonly realized when some mechanism like radial drift or fragmentation limits dust growth (see, e.g., Birnstiel et al. 2012; Lambrechts & Johansen 2014). Substituting $| {v}_{{\rm{r}}}| \approx 2\eta {v}_{{\rm{K}}}{\rm{St}}$ (${\rm{St}}\ll 1$) into Equation (24), we obtain the relation between ${{\rm{\Sigma }}}_{{\rm{d}}}$ and ${\rm{St}}$,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}}}=\displaystyle \frac{{\dot{M}}_{{\rm{d}}}}{2\pi r| {v}_{{\rm{r}}}| }\approx \displaystyle \frac{{\dot{M}}_{{\rm{d}}}}{4\pi \eta {{rv}}_{{\rm{K}}}{\rm{St}}}.\end{eqnarray} \tag{ 25 }$$

Because ${\rm{St}}\propto {a}_{*}$, we have ${{\rm{\Sigma }}}_{{\rm{d}}}\propto 1/{a}_{*}$ for constant ${\dot{M}}_{{\rm{d}}}$.

Using Equation (25) with the assumption that either radial drift or fragmentation limits dust growth, one can analytically estimate the radial distribution of ${{\rm{\Sigma }}}_{{\rm{d}}}$ and a_* in the non-sintering and sintering zones. When radial drift limits dust growth, the maximum aggregate size is determined by the condition (Okuzumi et al. 2012)

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{{\rm{coll}}}}{{t}_{{\rm{drift}}}}\approx \displaystyle \frac{1}{30},\end{eqnarray} \tag{ 26 }$$

where ${t}_{{\rm{coll}}}$ is the collision timescale given by Equation (12) and

$$\begin{eqnarray}&&{t}_{{\rm{drift}}}\equiv \displaystyle \frac{r}{| {v}_{{\rm{r}}}| }\approx \displaystyle \frac{1}{2\eta {\rm{\Omega }}{\rm{St}}}\end{eqnarray} \tag{ 27 }$$

is the timescale of radial drift. Substituting ${H}_{{\rm{d}}}\approx {(1+{\rm{St}}/{\alpha }_{{\rm{t}}})}^{-1/2}{H}_{{\rm{g}}}$ (${\rm{St}}\ll 1$) and ${m}_{*}/\pi {a}_{*}^{2}=(4/3){\rho }_{{\rm{int}}}{a}_{*}=(8/3\pi ){{\rm{\Sigma }}}_{{\rm{g}}}{\rm{St}}$ into Equation (13), the collision timescale can be evaluated as

$$\begin{eqnarray}&&{t}_{{\rm{coll}}}\approx 0.53\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{{H}_{{\rm{g}}}{\rm{St}}}{{\rm{\Delta }}v\sqrt{1+{\rm{St}}/{\alpha }_{{\rm{t}}}}}.\end{eqnarray} \tag{ 28 }$$

Furthermore, the collision velocity can be approximated as the root square sum of the turbulence-driven velocity ${\rm{\Delta }}{v}_{{\rm{t}}}$ and differential radial drift velocity ${\rm{\Delta }}{v}_{{\rm{r}}}$,

$$\begin{eqnarray}&&{\rm{\Delta }}v\approx \sqrt{{({\rm{\Delta }}{v}_{{\rm{t}}})}^{2}+{({\rm{\Delta }}{v}_{{\rm{r}}})}^{2}}\approx \sqrt{2.3{\alpha }_{{\rm{t}}}{c}_{{\rm{s}}}^{2}{\rm{St}}+{(\eta {v}_{{\rm{K}}}{\rm{St}})}^{2}},\end{eqnarray} \tag{ 29 }$$

where we have used the approximate expressions for ${\rm{\Delta }}{v}_{{\rm{t}}}$ and ${\rm{\Delta }}{v}_{{\rm{r}}}$ given in Section 4.3. Substituting Equations (25) and (27)–(29) into Equation (26), we obtain the equation for the maximum Stokes number in the drift-limited growth,

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{St}}}_{{\rm{drift}}}^{5/2}}{\sqrt{(2.3{\alpha }_{{\rm{t}}}{c}_{{\rm{s}}}^{2}+{\eta }^{2}{v}_{{\rm{K}}}^{2}{{\rm{St}}}_{{\rm{drift}}})(1+{{\rm{St}}}_{{\rm{drift}}}/{\alpha }_{{\rm{t}}})}}=\displaystyle \frac{0.0025{\dot{M}}_{{\rm{d}}}}{{\eta }^{2}{{rv}}_{{\rm{K}}}{c}_{{\rm{s}}}{{\rm{\Sigma }}}_{{\rm{g}}}},\end{eqnarray} \tag{ 30 }$$

where we labeled ${\rm{St}}$ by the subscript "drift" to emphasize drift-limited growth. In Figure 9(c), we compare ${{\rm{St}}}_{{\rm{drift}}}$ to ${\dot{M}}_{{\rm{d}}}=5\times {10}^{-9}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}$ with the Stokes number directly obtained from run Sa0 at $t={t}_{{\rm{snap}}}$. We find that ${{\rm{St}}}_{{\rm{drift}}}$ reproduces ${\rm{St}}$ in the non-sintering zones, implying that radial drift limits the growth of unsintered aggregates. We are also able to estimate a_* and ${{\rm{\Sigma }}}_{{\rm{d}}}$ in the non-sintering zones by substituting ${\rm{St}}={{\rm{St}}}_{{\rm{drift}}}$ into ${a}_{*}=(2/\pi ){{\rm{\Sigma }}}_{{\rm{g}}}{\rm{St}}/{\rho }_{{\rm{int}}}$ and Equation (25), respectively. These are shown by the dashed lines in Figures 9(a) and (b).

Zoom In Zoom Out Reset image size

If fragmentation limits dust growth, the maximum Stokes number is simply determined by the balance

$$\begin{eqnarray}&&{\rm{\Delta }}v({\rm{St}})={\rm{\Delta }}{v}_{{\rm{frag}}}.\end{eqnarray} \tag{ 31 }$$

We will denote the solution to this equation by ${{\rm{St}}}_{{\rm{frag}}}$. If we approximate ${\rm{\Delta }}v$ by Equation (29), Equation (31) can be rewritten as a quadratic equation for ${\rm{St}}$, and its positive root gives

$$\begin{eqnarray}&&{{\rm{St}}}_{{\rm{frag}}}=\displaystyle \frac{-2.3{\alpha }_{{\rm{t}}}{c}_{{\rm{s}}}^{2}+\sqrt{{(2.3{\alpha }_{{\rm{t}}}{c}_{{\rm{s}}}^{2})}^{2}+4{(\eta {v}_{{\rm{K}}}{\rm{\Delta }}{v}_{{\rm{frag}}})}^{2}}}{2{(\eta {v}_{{\rm{K}}})}^{2}}.\end{eqnarray} \tag{ 32 }$$

The dotted–dashed and dotted lines in Figure 9(c) show ${{\rm{St}}}_{{\rm{frag}}}$ for ${\rm{\Delta }}{v}_{{\rm{frag}}}={\rm{\Delta }}{v}_{{\rm{frag,NS}}}$ and ${\rm{\Delta }}{v}_{{\rm{frag,S}}}$ (denoted by ${{\rm{St}}}_{{\rm{frag,NS}}}$ and ${{\rm{St}}}_{{\rm{frag,S}}}$), respectively. We can see that ${{\rm{St}}}_{{\rm{frag,S}}}$ reproduces ${\rm{St}}$ in the sintering zones, which confirms that fragmentation limits the growth of sintered aggregates. One can estimate the values of a_* and ${{\rm{\Sigma }}}_{{\rm{d}}}$ in the sintering zones by substituting ${\rm{St}}={{\rm{St}}}_{{\rm{frag,S}}}$ into ${a}_{*}=2{{\rm{\Sigma }}}_{{\rm{g}}}{\rm{St}}/(\pi {\rho }_{{\rm{int}}})$ and Equation (25). These estimates are in excellent agreement with the simulation results, as shown by the dotted lines of Figures 9(a) and (b).

It is worth mentioning at this point that the radial pattern of dust shown in Figure 8 fades out as the disk becomes depleted of dust. As the dust-to-gas mass ratio ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ decreases, the collision timescale ${t}_{{\rm{coll}}}$ of the aggregates increases, and consequently the maximum aggregate size set by radial drift (Equation (26)) decreases. The radial pattern disappears when the radial drift barrier dominates the fragmentation barrier at all r because sintering has no effect on radial drift. This is illustrated in Figure 10, where we plot the radial distribution of a_* and ${{\rm{\Sigma }}}_{{\rm{d}}}$ for model Sa0 at different values of t. We find that the radial pattern present at $t={t}_{{\rm{snap}}}=0.26\;\mathrm{Myr}$ has disappeared by $t=1\;\mathrm{Myr}$. In this particular example, ${{\rm{St}}}_{{\rm{drift}}}$ falls below ${{\rm{St}}}_{{\rm{frag,S}}}$ at all r when ${M}_{{\rm{d}}}\lt {10}^{-11}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}$. In models La0 and LLa0, which assume weaker turbulence, dust evolution is slower than in model Sa0 owing to the lower turbulence-driven collision velocity (see ${t}_{{\rm{snap}}}$ in Table 3). However, even in these cases, the sintering-induced ring patterns are found to decay in 2 Myr. We note that the lifetime of the pattern would be longer for radially more extended (${r}_{{\rm{c}}}\gt 150\;\mathrm{au}$) disks, because the lifetime of dust flux in a disk generally scales with the orbital period at the disk's outer edge (Sato et al. 2015).

Zoom In Zoom Out Reset image size

We move on to the observational appearance of the sintering-induced dust rings at millimeter wavelengths. The upper left panel of Figure 11 shows the radial distribution of the line of sight optical depths ${\tau }_{\nu }$ (Equation (22)) from the snapshot of run Sa0 at $t={t}_{{\rm{snap}}}$. We here present the optical depths at three wavelengths $\lambda \quad =$ 0.87, 1.3, and 2.9 mm, which correspond to ALMA Bands 7, 6, and 3, respectively.

Zoom In Zoom Out Reset image size

Overall, the optical depths in the sintering zones are higher than in the non-sintering zones. This mainly reflects the higher dust surface density in the sintering zones (see Figure 8(b)). The radial variation of ${\kappa }_{{\rm{d}},\nu }$ represents only a minor contribution to the radial variation of ${\tau }_{\nu }$, in particular in outer regions where the representative aggregates are smaller than $\sim 1\;{\rm{cm}}$ in radius. At $\lambda \;=$ 0.87 and 1.3 mm, the three inner sintering zones (of ${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{NH}}_{3}$–${\mathrm{CO}}_{2}$–${{\rm{H}}}_{2}{\rm{S}}$, and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$) are optically thick, while the two outer sintering zones (of ${\mathrm{CH}}_{4}$ and ${\rm{CO}}$) are optically thin or marginally thick. The ${\rm{CO}}$ sintering zone is much darker than the other sintering zones because the disk surface density drops at $r\gt {r}_{{\rm{c}}}=150\;\mathrm{au}$. The non-sintering zones are optically thin or marginally thick at all three wavelengths. The opacity index at $\lambda =$ 0.87–1.3 mm, ${\beta }_{0.87-1.3\mathrm{mm}}\equiv \mathrm{ln}({\kappa }_{\text{0.87 mm}}/{\kappa }_{\text{1.3 mm}})$ $/\mathrm{ln}({\nu }_{\text{0.87 mm}}/{\nu }_{\text{1.3 mm}})$ is shown in the lower left panel of Figure 11. We see that ${\beta }_{0.87-1.3\mathrm{mm}}\sim 1$ at $\lesssim 70\;\mathrm{au}$ and approaches the interstellar value ∼1.7 beyond 80 au.

The upper right panel of Figure 11 shows the distribution of the brightness temperatures ${T}_{{\rm{B}}}$ for the same snapshot. For comparison, we also plot the gas temperature T of our disk model given by Equation (1). We find that the three innermost sintering zones are optically thick (${T}_{{\rm{B}}}\approx T$) at $\lambda \quad =$ 0.87 and 1.3 mm.

An interesting observational signature of the sintering-induced rings appears in the radial variation of the spectral slope. In the lower right panel of Figure 11, we plot the spectral index at 0.87–1.3 mm, ${\alpha }_{0.87-1.3\mathrm{mm}}\equiv \mathrm{ln}({I}_{\text{0.87 mm}}/{I}_{\text{1.3 mm}})$ $/\mathrm{ln}({\nu }_{\text{0.87 mm}}/{\nu }_{\text{1.3 mm}})$, as a function of r. In the three innermost sintering zones, we have ${\alpha }_{0.87-1.3\mathrm{mm}}\approx 2$ because these zones are optically thick at these wavelengths. If these regions were optically thin, we would have ${\alpha }_{0.87-1.3\mathrm{mm}}\approx 3$ because ${\beta }_{0.87-1.3\mathrm{mm}}\approx 1$ (see the lower left panel of Figure 11). In the non-sintering zones lying at ∼6–11 au and ∼33–80 au, we obtain ${\alpha }_{0.87-1.3\mathrm{mm}}\approx 2.3$–2.5, which is between the values in the optically thick and thin limits. This reflects the fact that the non-sintering zones are marginally thick at 0.87–1.3 mm.

Now we make more detailed comparisons between the simulation results and the ALMA observation of the HL Tau disk. We smooth the radial profiles of the intensities I_ν from run Sa0 ($t={t}_{{\rm{snap}}}$) at the ALMA resolutions as described in Section 4.6. In the center panels of Figure 12, the solid lines show the radial profiles of the brightness temperatures ${T}_{{\rm{B}}}$ and spectral slope ${\alpha }_{{\rm{B}}6-{\rm{B}}7}$ obtained from the smoothed I_ν. For comparison, ${T}_{{\rm{B}}}$ and ${\alpha }_{{\rm{B}}6-{\rm{B}}7}$ from the raw I_ν are shown by the dotted lines. The left panels show the profiles from the ALMA observations (same as Figure 1). We also show in the right panels the results from the non-sintering run Sa0-NoSint to clarify the features of the sintering-indued structures.

Zoom In Zoom Out Reset image size

After smoothing, the innermost emission dip lying at 6–11 au has been partially smeared at 0.87 mm (Band 7) and 1.3 mm (Band 6). The emission dip in the smoothed images has a lower spectral slope than that in the raw images directly obtained from the simulation. This is a consequence of the frequency-dependent angular resolution: Because Band 6 has a coarser resolution than Band 7, the emission dip seen at Band 6 is more significantly buried than that seen at Band 7, resulting in a decrease in the spectral slope after smoothing. At 2.9 mm (Band 3), the radial structure of ${T}_{{\rm{B}}}$ at $\lesssim 10\;{\rm{au}}$ has been significantly smoothed out.

We find that our simulation reproduces many observational features of the HL Tau disk. First, the simulation predicts a central emission peak that closely resembles the observed one. This central emission peak is associated with the ${{\rm{H}}}_{2}{\rm{O}}$ sintering zone as shown in Figure 11. The radial extent of the simulated central peak is $\approx 6\;\mathrm{au}$, which is comparable to $\approx 8\;\mathrm{au}$ of the observed peak. The simulation perfectly reproduces the emission features of the central region: the magnitudes and radial slopes of ${T}_{{\rm{B}}}$ at all three ALMA Bands, and the millimeter spectral slope of ${\alpha }_{{\rm{B}}6-{\rm{B}}7}\approx 2$. In our simulation, the spectral slope simply reflects the high optical thickness of the central region, and has nothing to do with the optical properties of the aggregates in the region. Sintering plays an essential role in the buildup of the optically thick region; without sintering, the disk would be entirely optically thin at these bands, as shown by run Sa0-NoSint (see the right panels of Figure 12). The fact that the central emission peak has a lower intensity at Band 3 than at Bands 6 and 7 is explained as a consequence of the lower spatial resolution ($\approx 10\;\mathrm{au}$) at Band 3 (i.e., this compact emission is underresolved at this band).

The simulation also predicts two bright rings in the region of ∼10–30 au that might be identified with the two innermost bright rings of HL Tau observed at ∼20 and 40 au. In our simulation, the two emission rings are associated with the sintering zones of ${\mathrm{NH}}_{3}$–${\mathrm{CO}}_{2}$–${{\rm{H}}}_{2}{\rm{S}}$ (11–23 au) and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ (24–33 au). These rings are optically thick at Bands 6 and 7 (and therefore ${\alpha }_{{\rm{B}}6-{\rm{B}}7}\approx 2$) and are optically thin at Band 3. These features are consistent with those of the two innermost bright rings of HL Tau. However, the separation of the predicted rings is much smaller than in the observed rings, as we discuss below.

Farther out in the disk, the simulation predicts an optically thin emission peak at 80 au associated with ${\mathrm{CH}}_{4}$ sintering. Its location coincides with the 80 au bright ring in the ALMA image of HL Tau, and its brightness temperatures agree with those of the observed ring at all three wavelengths to within a factor of two. The simulation also predicts a less pronounced peak at 120 au associated with ${\rm{CO}}$ sintering. As mentioned in Section 5.4, this 120 au peak is much less pronounced than other inner peaks because this location is close to the outer edge of our modeled gas disk. Interestingly, the observed HL Tau disk also has a minor emission peak exterior to the 80 au ring ($r\sim 97\;\mathrm{au};$ see ALMA Partnership et al. 2015).

The innermost dark ring seen at 6–11 au in the simulated image is optically marginally thick and has ${\alpha }_{{\rm{B}}6-{\rm{B}}7}\gt 2$, which is consistent with the observed innermost dark ring at $\sim 13\;\mathrm{au}$. Our simulation also explains why this innermost emission dip is much shallower at Band 3 than at Bands 6 and 7: as mentioned, this is simply because the spatial resolution at Band 3 is no high enough to distinguish the dark ring from the central emission peak.

However, there are some discrepancies between the prediction from the fiducial model and the observation of the HL Tau disk. For example, the second innermost ring predicted by the model, which is associated with the ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ sintering zone, is about 10 au interior to that observed by ALMA. For this reason, the dark ring just inside the ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ sintering zone is much narrower than the second innermost dark ring of HL Tau extending from 30 to 40 au. In addition, as we explain in detail in Section 6.3, the fiducial model predicts that dust particles are vertically well mixed in the gas disk, which seems to be inconsistent with the observations of HL Tau, suggesting that large dust particles settle to the disk midplane (Kwon et al. 2011; Pinte et al. 2016). In the following section, we examine if these discrepancies can be removed or alleviated by tuning the parameters in our model.

Because the gas density distribution of the HL Tau disk is unknown, it is important to quantify how strongly our results depend on the assumption about the gas density profile. We present the results of two simulation runs—Sa0-Lgam and Sa0-Hgam—in which we change γ in Equation (2) to 0.5 and 1.5, respectively. The other parameters are unchanged form the fiducial run.

Figure 13 compares the radial profiles of the brightness temperature and spectral slope at 0.87 and 2.9 mm from these runs with those from run Sa0. The three snapshots are taken at different times but still provide similar flux densities (see Table 3) because ${t}_{{\rm{snap}}}$ is defined as such. We find that the variation of γ within the range 0.5–1.5 little affects the emission properties of the disk, with the variation of ${T}_{{\rm{B}}}$ within a factor of 2 and the variation of ${\alpha }_{0.87-1.3\mathrm{mm}}$ as small as $\lesssim 25\%$ at all r. This is mainly because the steady-state radial distribution of the dust surface density ${{\rm{\Sigma }}}_{{\rm{d}}}$ is insensitive to γ. If we measure the aggregate size by St, the steady-state distribution of ${{\rm{\Sigma }}}_{{\rm{d}}}$ is determined from Equation (25) with either ${\rm{St}}={{\rm{St}}}_{{\rm{drift}}}$ (Equation (30)) or ${\rm{St}}={{\rm{St}}}_{{\rm{frag}}}$ (Equation (32)). The γ dependence of the radial dust flux ${\dot{M}}_{{\rm{d}}}$ is weak as long as we fix ${M}_{{\rm{disk}}}$ and ${r}_{{\rm{c}}}$. ${{\rm{St}}}_{{\rm{frag}}}$ depends on γ only through $\eta \propto \gamma +1.8$ (for $T\propto {r}^{-0.57}$ and $r\ll {r}_{{\rm{c}}}$), and its variation is small as long as we vary γ within the range 0.5–1.5. For ${{\rm{St}}}_{{\rm{drift}}}$, the dependence is less obvious from Equation (30), but it turns out that the weak γ dependences of ${\eta }^{2}{{\rm{\Sigma }}}_{{\rm{d}}}$ and ${\dot{M}}_{{\rm{d}}}$ partly cancel out in this equation.

Zoom In Zoom Out Reset image size

The radial distribution of turbulence strength ${\alpha }_{{\rm{t}}}$ is another important uncertainty in our simulations. Our fiducial model assumes ${\alpha }_{{\rm{t}}}\propto \sqrt{r}$, which gives a turbulence-driven collision velocity nearly independent of r (${\rm{\Delta }}{v}_{{\rm{t}}}\;\propto \sqrt{{\alpha }_{{\rm{t}}}}{c}_{{\rm{s}}}\propto \sqrt{{\alpha }_{{\rm{t}}}T}\propto {r}^{-0.035}$). In fact, a radially constant collision velocity is required in our model to simultaneously reproduce dark rings at small r and bright rings at large r.⁷ To demonstrate this, we consider two models in which ${\alpha }_{{\rm{t}}}$ is fixed to 0.03 and 0.1 at all r (models Sa0-Lalp and Sa0-Halp, respectively). These values correspond to the values of ${\alpha }_{{\rm{t}}}$ in the fiducial Sa0 model at $r\approx 10\;\mathrm{au}$ and $100\;\mathrm{au}$, respectively.

Figure 14 compares the result of fiducial run Sa0 with those of runs Sa0-Lalp and Sa0-Halp. We find that model Sa0-Lalp fails to produce an emission peak at $\approx 80\;\mathrm{au}$, because turbulence is too weak to cause fragmentation of sintered aggregates at that location. By contrast, model Sa0-Halp produces only shallow dips at $r\lesssim 30\;\mathrm{au}$. In this high-${\alpha }_{{\rm{t}}}$ model the maximum aggregate size at $r\lesssim 30\;\mathrm{au}$ is limited by turbulent fragmentation, even in the non-sintering zones. The suppressed dust growth causes a slowdown of radial drift, which acts to fill the density dips in the non-sintering zones. Thus, a model assuming a high ${\alpha }_{{\rm{t}}}$ in the outer regions and a low ${\alpha }_{{\rm{t}}}$ in the inner regions best reproduces the observation of HL Tau.

Zoom In Zoom Out Reset image size

Theoretically, a ${\alpha }_{{\rm{t}}}$ that is lower at smaller r is expected for turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley 1991). Because the origin of the MRI is the coupling between the gas disk and magnetic fields, MRI turbulence tends to be weaker at locations where the ionization degree is lower. In protoplanetary disks, a lower ionization degree corresponds to a higher gas density and hence to a smaller orbital radius, because ionizing cosmic-rays or X-rays are attenuated at large column densities and recombination is faster in denser gas (e.g., Sano et al. 2000; Bai 2011).

As mentioned in Section 4.6, our calculations of dust thermal emission neglect the geometrical thickness of the dust subdisk. In reality, if a dust disk has an inclination of $\sim 45^\circ$ and a finite vertical extent ${H}_{{\rm{d}}}$, then any radial structure of the disk emission is smeared out over this length scale in the direction of the minor axis of the disk image. Pinte et al. (2016) point out that the scale height of large dust particles in the HL Tau disk must be as small as $\sim 1\;{\rm{au}}$ at $r=100\;\mathrm{au}$ in order to be consistent with the well separated morphology of the bright rings observed by ALMA. Such a small dust scale height strongly indicates that dust settling has occurred in the gas disk; without settling, the dust scale height would be $\sim 10\;{\rm{au}}$ at 100 au.

However, in our fiducial model, the settling of representative aggregates is severely prevented by turbulent diffusion. According to Equation (14), significant settling of dust particles (${H}_{{\rm{d}}}\ll {H}_{{\rm{g}}}$) requires ${\rm{St}}\gg {\alpha }_{{\rm{t}}}$. This condition is not satisfied in fiducial run Sa0, because the value of ${\rm{St}}$ of the representative aggregates observed in the simulation is comparable to the value of ${\alpha }_{{\rm{t}}}$ assumed (see Figure 8(c)). In this run, ${\alpha }_{{\rm{t}}}$ is arranged to have a high value so that aggregates disrupt and pile up in the sintering zones. If turbulence were weak, sintered aggregates would not experience disruption as long as we maintain the assumption ${\rm{\Delta }}{v}_{{\rm{frag,S}}}=20\;{\rm{m}}\;{{\rm{s}}}^{-1}$.

One way to reconcile sintering-induced ring formation with dust settling in our model is to assume weaker turbulence and a lower fragmentation velocity. Within our dust model, a lower value of ${\rm{\Delta }}{v}_{{\rm{frag}}}$ corresponds to a larger monomer size a₀ (see Equations (18) and (19)). However, aggregates made of larger monomers tend to be sintered less slowly, as demonstrated in Section 3.4. To examine if there is a range of a₀ where the aggregates can experience settling and sintering simultaneously, we performed two simulations (La0 and LLa0). In run La0, we increase a₀ to $1\;\mu {\rm{m}}$ while lowering ${\alpha }_{{\rm{t}}}$ to ${10}^{-3}{(r/10\mathrm{au})}^{1/2}$. Run LLa0 is a more extreme case of ${a}_{0}=4\;\mu {\rm{m}}$ and ${\alpha }_{{\rm{t}}}={10}^{-4}{(r/10\mathrm{au})}^{1/2}$. The other parameters are the same as in the fiducial run Sa0.

Figure 15 shows the raw radial profiles of ${T}_{{\rm{B}}}$ and the 0.87–1.3 mm spectral index at $t={t}_{{\rm{snap}}}$ obtained from runs La0 and LLa0. The profiles after smoothing at the ALMA resolutions are shown in Figure 16. Note that we neglect the geometric thickness of the dust disk when calculating ${T}_{{\rm{B}}}$. Indeed, the dust disks in models La0 and LLa0 are geometrically thin, unlike in model Sa0 because turbulent diffusion is inefficient in these two models. We plot in Figure 17 the dust scale height ${H}_{{\rm{d}}}$ from Equation (14) as a function of r for the three models. At $r=100\;\mathrm{au}$, we find ${H}_{{\rm{d}}}\approx 3\;\mathrm{au}$ for model La0 and ${H}_{{\rm{d}}}\approx 1\;\mathrm{au}$ for model LLa0. Models Sa0, La0, and LLa0 give similar results for the radial distribution of ${T}_{{\rm{B}}}$, particularly in inner disk regions. However, model LLa0 fails to produce an emission peak at 80 au because the ${\mathrm{CH}}_{4}$ sintering zone completely disappears for ${a}_{0}=4\;\mu {\rm{m}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, we find that dust settling gives a strong constraint on the turbulence strength and monomer grain size in the HL Tau disk. Dust settling and pile up occur simultaneously only if ${10}^{-4}\lt {\alpha }_{{\rm{t}}}\lesssim {10}^{-3}$ and $1\;\mu {\rm{m}}\lesssim {a}_{0}\lt 4\;\mu {\rm{m}}$. If ${a}_{0}\gtrsim 4\;\mu {\rm{m}}$, sintering would be too slow to provide an appreciable emission peak at $\sim 80\;\mathrm{au}$. If ${a}_{0}\ll 1\;\mu {\rm{m}}$, sintered aggregates would be disrupted only when ${\alpha }_{{\rm{t}}}\gg {10}^{-3}$, but such a strong turbulence would inhibit dust settling.

Interestingly, the above results suggest that the grains constituting the aggregates in the HL Tau disk are considerably larger than the interstellar grains ($\lesssim 0.25\;\mu {\rm{m}}$ in radius; e.g., Mathis et al. 1977). They are also larger than the grains constituting interplanetary dust particles of presumably cometary origin (typically 0.1–0.5 μm in diameter; e.g., Rietmeijer 1993). However, such a large grain size is not excluded by previous observations of HL Tau. The near-infrared scattered light images of the envelope of HL Tau are best reproduced by models that assume the maximum particle size of more or less $1\;\mu {\rm{m}}$ (Lucas et al. 2004; Murakawa et al. 2008). Because the scattered light probes the envelope's surface layer where coagulation is inefficient, the observed micron-sized particles might be monomers rather than aggregates.

As mentioned in Section 5.5, the fiducial model does not fully explain the exact locations and widths of all major HL Tau rings. For example, the innermost dark ring predicted by model Sa0 lies somewhat closer to the central star than the observed one. If we define the position of the innermost dark ring as the location where ${T}_{{\rm{B}}}/T$ at Band 6 is maximized, the radius of the predicted innermost dark ring (8 au) is $\approx 40$% smaller than that of the observed one (13 au). Furthermore, the second innermost dark ring in model Sa0 is much narrower than that observed because the ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ snow line lies very close to the ${{\rm{H}}}_{2}{\rm{S}}$ sintering line. Of course, a different temperature profile would provide a different configuration of the sintering zones. However, as long as T(r) is assumed to obey a single power law, it is generally impossible to move some of the sintering zones while leaving the positions of the others unchanged. For example, a temperature profile that is slightly higher than Equation (1) would shift the ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ sintering zone to 40 au, but would also shift the ${\mathrm{CH}}_{4}$ sintering zone to 100 au. Moreover, a different T(r) would make our prediction for the brightness temperatures in the optically thick regions less good.

However, the locations of the sintering zones depend not only on the gas temperature profile but also on the sublimation energies L_j in ${P}_{{\rm{ev}},j}(T)$. As mentioned in Section 3.2, there is typically a 10% uncertainty in the published data of L_j. In general, a 10% uncertainty in the sublimation energy causes a $\sim 10\%$ uncertainty in the sublimation temperature ${T}_{{\rm{subl}},j}$ because ${P}_{{\rm{ev}},j}$ is a function of the ratio ${L}_{j}/T$. Assuming the temperature profile given by Equation (1), we have ${r}_{{\rm{snow}},j}\propto {T}_{{\rm{subl}},j}^{-1.8}$ and hence $| \delta {r}_{{\rm{snow}},j}| /{r}_{{\rm{snow}},j}\approx 1.8| \delta {T}_{{\rm{subl}},j}| /{T}_{{\rm{subl}},j}\sim 2| \delta {L}_{j}| /{L}_{j}$, where $\delta {L}_{j}$, $\delta {T}_{{\rm{subl}},j}$, and $\delta {r}_{{\rm{snow}},j}$ denote the uncertainties of L_j, ${T}_{{\rm{subl}},j}$, and ${r}_{{\rm{snow}},j}$, respectively. Consequently, a 10% uncertainty in L_j leads to a $\sim 20\%$ uncertainty in the snow line location. Such an uncertainty can be significant because the separation of the observed rings is only a fraction of their radii.

To demonstrate the potential importance of uncertainties in the sublimation energies, we performed two simulations (Sa0-tuned and La0-tuned). The parameters adopted in these simulations are the same as those in Sa0 and La0, respectively, except that the sublimation energies of ${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{NH}}_{3}$, and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ are lowered by 10% from our baseline values (see column 5 of Table 2). According to the estimate shown previously, such modifications shift the sintering zones of these volatiles outward by $\sim 20\%$. The resulting radial profiles of ${T}_{{\rm{B}}}$ and ${\alpha }_{0.87-1.3\mathrm{mm}}$ before smoothing are shown in Figure 18. In model Sa0-tuned, the three innermost sintering zones are located at 4–8 au (${{\rm{H}}}_{2}{\rm{O}}$), 13–23 au (${\mathrm{NH}}_{3}$–${\mathrm{CO}}_{2}$–${{\rm{H}}}_{2}{\rm{S}}$), and 29–40 au (${{\rm{C}}}_{2}{{\rm{H}}}_{6}$) instead of 3–6 au, 11–23 au, and 24–33 au as in model Sa0. Figure 19 shows the radial emission profiles after smoothing, which we compare with the ALMA observation. For the sake of comparison, we define the radii of the bright and dark rings by the orbital radii at which ${T}_{{\rm{B}}}/T$ at Band 6 is locally maximized and minimized, respectively. Then, the radii of four innermost bright/dark rings are found to be 13, 23, 32, and 38 au for the observation, 9, 16, 24, and 30 au for model Sa0-NoSint, and 9, 15, 24, and 29 au for model La0-NoSint. Thus, model Sa0-NoSint reproduces the ring radii in the observed image to an accuracy of $\lesssim 30\%$, which is about 10% better than model Sa0. Model La0-NoSint is slightly less accurate, with a maximum error of 35%, but it is still 10% better than the original model La0. Furthermore, the second innermost dark ring in these models is much wider than that in model Sa0. As a consequence, the radial distribution of ${\alpha }_{{\rm{B}}6-{\rm{B}}7}$ now has a clear peak structure at the position of this ring, and the peak value $\approx 2.3$ agrees with the observed value $\approx 2.5$ to within a relative error of 10%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Thus, assuming sublimation energies of ${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{NH}}_{3}$, and ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$ that are only 10% lower than the baseline values significantly improves our predictions for the radial configuration of the dust rings. Interestingly, for ${{\rm{H}}}_{2}{\rm{O}}$ and ${\mathrm{NH}}_{3}$, the tuned values of the sublimation energies are more consistent with the results of recent experiments by Martín-Doménech et al. (2014; 5165 K for ${{\rm{H}}}_{2}{\rm{O}}$ and 2965 K for ${\mathrm{NH}}_{3}$) than our fiducial values. However, not all new data for sublimation energies improve our predictions. Martín-Doménech et al. (2014) also measured the sublimation energies of ${\mathrm{CO}}_{2}$ and ${\rm{CO}}$, and the measured values (2605 K for ${\mathrm{CO}}_{2}$ and 890 K for ${\rm{CO}}$) are also lower than ours by about 10%. A 10% change in the sublimation energy of ${\mathrm{CO}}_{2}$ has little effect on the resulting ring patterns, because the ${\mathrm{CO}}_{2}$ sintering zone partially overlaps with the sintering zones of ${\mathrm{NH}}_{3}$ and ${{\rm{H}}}_{2}{\rm{S}}$. A 10% decrease in the sublimation energy of ${\rm{CO}}$ shifts the ${\rm{CO}}$ sintering zone to 143–197 au, which makes the correspondence between the CO sintering zone and the faint 97 au ring of HL Tau (see Section 5.5) less good.

As mentioned in the introduction, sintered aggregates tend to bounce rather than stick when they collide at a velocity below the fragmentation threshold. For example, sintered aggregates made of 0.1 μm sized icy grains bounce at $1\;{\rm{m}}\;{{\rm{s}}}^{-1}\lesssim {\rm{\Delta }}v\;\lesssim 20\;{\rm{m}}\;{{\rm{s}}}^{-1}$ (Sirono 1999; S. Sirono 2016, in preparation). This effect was neglected in this study by simply applying Equation (15) to both sintered and unsintered aggregates. In principle, such a simplification causes an overestimate of the maximum size of aggregates in the sintering zones. However, this effect is expected to be minor because sintered aggregates grow only moderately even without the bouncing effect. To see this, we go back to the radial profile of the radius a_* of the representative aggregates from run Sa0 shown in Figure 8(a). Because the radial inward flow of the aggregates is nearly stationary (see Figure 8(d)), this figure shows how the size of individual representative aggregates evolve as they move inward. We find that a_* is radially constant in the ${\mathrm{CH}}_{4}$ and ${\rm{CO}}$ sintering zones, meaning that the aggregates do not grow at all when they go through these zones. Appreciable growth of sintered aggregates occurs only in inner regions of the ${{\rm{H}}}_{2}{\rm{O}}$ and ${\mathrm{NH}}_{3}$–${\mathrm{CO}}_{2}$–${{\rm{H}}}_{2}{\rm{S}}$ sintering zones. However, even at these locations, the aggregate size increases only by a factor of less than two until they reach the inner edges of the sintering zones. Therefore, we can conclude that inclusion of bouncing collisions would little change the evolution of representative aggregates in the sintering zones.

Our single-size approach (Section 4.1) relies on the assumption that the mass budget of dust at each orbital radius is dominated by a single population of aggregates having mass $m={m}_{*}(r)$. This assumption might be inadequate at the boundaries of sintering and non-sintering zones, around which two populations of aggregates of different characteristic sizes (i.e., sintered and unsintered aggregates) can coexist. However, this effect would only be important in the close vicinity of the boundaries because the sintering timescale is a steep function of r and the aggregates in our simulations do not drift faster than they collide with each other.

A probably more critical limitation of the single-size approach is that one has to assume the size distribution of fragments produced by the collisions of the mass-dominating aggregates. In this study we avoided detailed modeling of the fragmentation process by assuming a simple power-law fragment size distribution (Equation (20)) independent of the collision velocity of the largest aggregates. The assumed power-law distribution would reasonably approximate the true fragment size distribution when the collisions of the mass-dominating aggregates are highly disruptive. In fact, however, unsintered aggregates in our simulations do not experience catastrophic disruption because their collision velocity is always below the catastrophic disruption threshold ${\rm{\Delta }}{v}_{{\rm{frag}}}$. For example, our fiducial run Sa0 shows that ${\rm{\Delta }}{m}_{*}/{m}_{*}\approx 0.2$–0.4 in the non-sintering zones, which implies that fragments would carry away only a few tens of percent of the total mass of two colliding mass-dominating aggregates in these zones. Equation (20) might also overestimate the amount of fragments from sintered aggregates because they actually bounce off, rather than fragment at low collision velocities (see Section 7.1). The preliminary results of our aggregate collision simulations (S. Sirono 2015, in preparation) show that fragments from two colliding sintered aggregates carry away only a few percent of their total mass even when ${\rm{\Delta }}v\approx {\rm{\Delta }}{v}_{{\rm{frag}}}$.

Therefore, it important to assess how the predictions from our models depend on the amount of fragments assumed. We consider a fragment size distribution that is similar to Equation (20), but assumes a reduced amount of fragments in the size range $a\lt {a}_{*}/10$,

$$\begin{eqnarray}{N}_{{\rm{d}}}^{\prime }(a)=\left\{\begin{array}{ll}{{Ca}}^{-3.5}, & {a}_{*}/10\lt a\lt {a}_{*},\\ {f}_{{\rm{red}}}{{Ca}}^{-3.5}, & {a}_{{\rm{min}}}\lt a\lt {a}_{*}/10,\\ 0, & {\rm{otherwise}},\end{array}\right.\end{eqnarray} \tag{ 33 }$$

where the factor ${f}_{{\rm{red}}}(\lt 1)$ encapsulates the reduction of fragment production, and C is again determined by the condition ${\int }_{0}^{\infty }{{mN}}_{{\rm{d}}}^{\prime }(a){\rm{d}}a={{\rm{\Sigma }}}_{{\rm{d}}}$. We consider two values— ${f}_{{\rm{red}}}=0.3$ and 0.03—based on the estimates for unsintered and sintered aggregates mentioned above.

We now recalculate the radial profiles of ${T}_{{\rm{B}}}$ and ${\alpha }_{{\rm{B}}6-{\rm{B}}7}$ for model Sa0-tuned using Equation (33) instead of Equation (20). The results are shown in the center and right panels of Figure 20. These, together with the result for ${f}_{{\rm{red}}}=1$ (the center panels in Figure 19), show that both ${T}_{{\rm{B}}}$ and ${\alpha }_{{\rm{B}}6-{\rm{B}}7}$ in the two innermost emission dips decrease as ${f}_{{\rm{red}}}$ is decreased. In these inner regions, the radius a_* of the largest (mass-dominating) aggregates is significantly larger than millimeters, and therefore the fragments smaller than millimeters have a non-negligible contribution to the millimeter dust opacity. In this case, the opacity index ${\beta }_{0.87-1.3\mathrm{mm}}$ decreases toward zero, which is the value in the geometric optics limit, as the amount of the fragments decreases (see, e.g., Draine 2006). This explains why a lower value of ${f}_{{\rm{red}}}$ leads to a lower value of the spectral slope in the optically thin inner dips where the relation ${\alpha }_{0.87-1.3\mathrm{mm}}\approx {\beta }_{0.87-1.3\mathrm{mm}}+2$ applies. In the case of ${f}_{{\rm{red}}}=0.03$, the spectral slope in the innermost emission dip falls below two, due to the effect of the frequency-dependent angular resolution already mentioned in Section 5.5. We find that the results for ${f}_{{\rm{red}}}=0.3$ and 0.03 better reproduce the observed depths of the two innermost emission dips than the result for ${f}_{{\rm{red}}}=1$. However, these low-${f}_{{\rm{red}}}$ models yield poorer agreement with the observed value of ${\alpha }_{{\rm{B}}6-{\rm{B}}7}$ in these emission dips. Varying the value of ${f}_{{\rm{red}}}$ within the range 0.03–1 has no effect on the predictions for ${T}_{{\rm{B}}}$ and ${\alpha }_{{\rm{B}}6-{\rm{B}}7}$ outside the two innermost emission dips. Taken together, we cannot judge at this point which value of ${f}_{{\rm{red}}}$ best reproduces the observed appearance of HL Tau. In any case, the effects of assuming the lower values of ${f}_{{\rm{red}}}$ are not significant as long as ${f}_{{\rm{red}}}$ is higher than 0.3, which is the value we expect for unsintered aggregates in the dark rings.

Zoom In Zoom Out Reset image size

We have assumed that vapor transport within an aggregate, which drives sintering, does not alter the size of monomer grains. In reality, this is not true because the growth of necks must be compensated by the shrinkage of the bodies of monomers. Furthermore, if the monomers within an aggregate are not uniform in size, a small number of large monomers may grow—the phenomenon known as Ostwald ripening (see, e.g., Lifshitz & Pitaevskii 1981)—while smaller ones may completely evaporate leaving silicate grains (Kuroiwa & Sirono 2011; Sirono 2011a).⁸ The evolution of the monomer size distribution is negligible in the sintering zones of minor volatiles like ${\mathrm{NH}}_{3}$ and ${\mathrm{CH}}_{4}$ (as noted in the introduction, sintering can still occur because the neck volume is much smaller than the monomer volume). However, this is not the case for water ice because it constitutes more than half of the monomer volume.

As pointed out by Kuroiwa & Sirono (2011), the evolution of the monomer size distribution due to water vapor transport would result in a decrease in the sticking efficiency of the icy aggregates. Kuroiwa & Sirono (2011) showed that water vapor transport within an aggregate produces a small number of large ice-rich monomers and a large number of small bare silicate monomers. Aggregates made of the two populations of monomers would be fragile, like sintered aggregates made of equally ice-coated grains, because the total binding energy (∝number of contacts) of the former would be determined by weak silicate–silicate contacts. These aggregates would experience catastrophic disruption and pile up near the ${{\rm{H}}}_{2}{\rm{O}}$ snow line in a similar way to sintered aggregates.

In addition, the evaporation of small monomers could directly result in the fragmentation of icy aggregates. Sirono (2011a) demonstrated that an aggregate made of initially submicron-sized ${{\rm{H}}}_{2}{\rm{O}}$ monomers breaks up into large detached monomers of radii ∼1–10 $\mu {\rm{m}}$ (see also Kuroiwa & Sirono 2011). This spontaneous breakup would also produce a pileup of dust near the ${{\rm{H}}}_{2}{\rm{O}}$ snow line, but the surface density contrast provided by this effect could be much more significant than that by than sintering-induced fragmentation, because the micron-sized monomers are much more strongly coupled to the disk gas than aggregates. Sirono (2011a) notes that the dust surface density in the ${{\rm{H}}}_{2}{\rm{O}}$ sintering zone can become high enough to trigger the formation of icy planetesimals via the gravitational instability (Goldreich & Ward 1973; Sekiya 1998; Youdin & Shu 2002).

Thus, the evolution of the monomer size distribution due to vapor transport would further promote the fragmentation and accumulation of aggregates in the ${{\rm{H}}}_{2}{\rm{O}}$ sintering zone. This may not affect the observational appearance of the HL Tau disk because the ${{\rm{H}}}_{2}{\rm{O}}$ sintering zone is already optically thick, even without the monomer size evolution. However, this could have a significant effect on planet(esimal) formation near the ${{\rm{H}}}_{2}{\rm{O}}$ snow line.

We neglected the evolution of the porosity of icy aggregates by assuming the fixed aggregate internal density of $0.24\;{\rm{g}}\;{\mathrm{cm}}^{-3}$. In reality, the porosity of aggregates (much smaller than planetesimals) can vary greatly through collisions (e.g., Ormel et al. 2007; Suyama et al. 2008; Okuzumi et al. 2009), compression by gas drag (Kataoka et al. 2013a, 2013b), and condensation and sintering of ${{\rm{H}}}_{2}{\rm{O}}$ ice (again, minor volatile ices occupy only a small fraction of grain volume and are therefore negligible here). In the absence of the compaction due to condensation and sintering, the internal density of icy aggregates in protoplanetary disks can decline to as low as 10⁻⁴–10⁻³ g cm⁻³ (Okuzumi et al. 2012; Kataoka et al. 2013b).

The porosity evolution has two important effects on the evolution of icy aggregates. First, fluffy aggregates tend to grow to the size regime where the gas drag onto the aggregates is determined by Stokes' law. This particularly occurs in the inner regions of protoplanetary disks, where the mean free path of the disk gas molecules is short. In the Stokes drag regime, the growth limit given by Equation (30) no longer applies because that limit assumes Epstein's drag law. Okuzumi et al. (2012) found that fluffy aggregates that are subject to Stokes' drag law grow fast enough to overcome the radial drift barrier. In a massive protoplanetary disk as considered in this study, icy aggregates break through the drift barrier and form icy planetesimals at $r\lesssim 10\;{\rm{au}}$ (see the bottom right panel of Kataoka et al. 2013a). Second, the millimeter absorption opacity of large porous aggregates is similar to that of small compact aggregates because ${\kappa }_{{\rm{d}},\nu }$ can be approximated as a function of the product ${\rho }_{{\rm{int}}}a$ (Kataoka et al. 2014).

In this study, these important effects have been neglected by fixing the aggregate porosity to a relatively low value. Rapid coagulation beyond the Epstein drag regime could strongly reduce the optical depth in the non-sintering zone lying at $r\sim 10\;{\rm{au}}$, where the radial drift barrier determines the maximum size of the relatively compact aggregates that we assumed in this study. This effect should be studied in detail in future work.

On the other hand, the above two effects are less important at $r\gtrsim 10\;{\rm{au}}$, where even the highly fluffy aggregates would grow within the Epstein drag regime. To show this, we first note that neither ${{\rm{St}}}_{{\rm{drift}}}$ (Equation (30)) nor ${{\rm{St}}}_{{\rm{frag}}}$ (Equation (32)) depends on ${\rho }_{{\rm{int}}}$ (note again that Epstein's law is used in the derivation of Equation (30)). This means that porosity evolution has no effect on the maximum attainable Stokes number at each r, and therefore does not change the radial distribution of ${{\rm{\Sigma }}}_{{\rm{d}}}$ because the radial advection velocity v_r is a function of ${\rm{St}}$ (see Equation (10)). Porosity evolution does not change the dust absorption opacity ${\kappa }_{{\rm{d}},\nu }$ either, because the product ${\rho }_{{\rm{int}}}{a}_{*}$ is proportional to ${\rm{St}}$ in the Epstein regime (see Equation (9)). This is a natural consequence of the fact that both the stopping time in the Epstein regime and the absorption opacity are determined by the mass-to-area ratio $m/\pi {a}^{2}(\sim {\rho }_{{\rm{int}}}a)$ of the aggregates. Taken together, the radial distribution of the absorption optical depth ${\tau }_{\nu }\propto {\kappa }_{{\rm{d}},\nu }{{\rm{\Sigma }}}_{{\rm{d}}}$ is independent of the aggregate porosity, as long as the aggregates grow within the Epstein drag regime.

Nevertheless, porosity evolution might have some important observational consequence in the outer disk region. The millimeter continuum emission from the HL Tau disk is known to be polarized in the direction roughly perpendicular to the disk major axis at a level of 0.9% in the average degree of polarization (Stephens et al. 2014). Yang et al. (2015) and Kataoka et al. (2015b) point out that self-scattering of thermal emission by dust rings (Kataoka et al. 2015a) can explain the polarization observation. However, this mechanism produces a degree of polarization of $\approx 1\%$ only when the size of the dust particles falls within a particular range. For example, if the particles are compact and their size distribution obeys our Equation (20), the self-scattering scenario is consistent with the observation only when the maximum particle size ${a}_{{\rm{max}}}$ ($={a}_{*}$ in this study) is $\sim 0.1\;{\rm{mm}}$ (Kataoka et al. 2015b; Yang et al. 2015), which is considerably smaller than the size of mass-dominating aggregates predicted in this study. In the compact particle model, a particle size larger than a millimeter is ruled out because the scattered emission would be too weakly polarized (see also Kataoka et al. 2015a). By contrast, large but fluffy aggregates are known to produce highly polarized scattered light like constituent monomers do (Min et al. 2016), and therefore might be able to explain the polarization observation of HL Tau.

We neglected the condensation growth of ice particles near the snow lines. In fact, a recent study by Ros & Johansen (2013) has shown that condensation growth can be effective at locations slightly outside the ${{\rm{H}}}_{2}{\rm{O}}$ snow line. This effect might change our predictions for dust evolution near ${{\rm{H}}}_{2}{\rm{O}}$ snow line because condensation dominantly takes place on the smallest particles (i.e., monomers). If condensation dominates over sintering-induced fragmentation and bouncing, the regions slightly outside the snow lines would be seen as ${\text{}}{dark}$ rings at millimeter wavelengths, as noted by Zhang et al. (2015). However, it is unclear whether condensation growth becomes important for volatiles less abundant than ${{\rm{H}}}_{2}{\rm{O}}$, like ${\rm{CO}}$ and ${\mathrm{NH}}_{3}$. To address this open question, we will incorporate condensation growth into our global dust evolution simulations in future work.

Perhaps the strongest prediction from our model, as well as from the model of Zhang et al. (2015), is that the multiple dust rings may not be peculiar to the HL Tau disk. In principle, sintering-induced ring formation operates as long as the disk is not too depleted of dust (or not too old), the monomers are sufficiently small to cause rapid sintering, and turbulence is strong enough to cause the disruption of sintered aggregates. If these conditions are satisfied, axisymmetric dust rings emerge slightly outside the snow lines of relatively major volatiles. An intriguing object in this context is TW Hya. For this system, the CO snow line was indirectly detected at an orbital distance of ∼30 au (Qi et al. 2013), and furthermore, near-infrared scattered light images of the disk suggest the presence of two axisymmetric dust gaps at ∼80 au (Debes et al. 2013) and ∼20 au (Akiyama et al. 2015; Rapson et al. 2015). A latest ALMA observation has shown that there is also a gap in millimeter dust emission in the vicinity of the 20 au scattered light gap (Nomura et al. 2015). The gaps at 80 au and 20 au could be associated with non-sintering zones exterior to and interior to the CO sintering zone, respectively. However, it is not obvious whether our ring formation mechanism directly applies to TW Hya because this object is known to be much older than HL Tau (at least 3 Myr; Vacca & Sandell 2011). This question should also be addressed in future work.