EVIDENCE OF FAST PEBBLE GROWTH NEAR CONDENSATION FRONTS IN THE HL TAU PROTOPLANETARY DISK

The frost lines of abundant ices have long been thought to be important for planet formation, in part because the extra material provided by condensation can significantly enhance the local mass surface density of solids in the disk (Hayashi 1981; Stevenson & Lunine 1988; Desch 2007). This in turn significantly shortens the timescale needed to form protoplanets in core-accretion models (Pollack et al. 1996). Furthermore, the bulk composition of terrestrial planets and the cores of gas or ice giants is determined by that of the planetesimals from which they are built. Because volatility is highly species specific, it is the condensation fronts of the principal solid-vapor reservoirs in the mid-plane of protoplanetary disks that set the rough boundaries of planetesimals with different compositions (Öberg et al. 2011).

Recent experiments and numerical simulations of dust coagulation and settling suggest that condensation fronts may play a crucial role in planetesimal formation via two aspects: first, the growth efficiency of dust grains is composition dependent. Lab experiments show that icy aggregates are significantly more "sticky," and resistant to compaction than silicate aggregates (Güttler et al. 2010; Seizinger & Kley 2013; Kelling et al. 2014). The second, perhaps more important, aspect is that the enhanced local surface mass density near snow lines may produce pressure bumps, high dust-to-gas mass ratios, or viscosity gradients, thus triggering instabilities that form planetesimals (Kretke & Lin 2007; Brauer et al. 2008; Ros & Johansen 2013; Bitsch et al. 2014; Dra̧żkowska & Dullemond 2014; Flock et al. 2015).

Despite the wealth of theoretical predictions concerning condensation fronts in disks, no direct observations of fast dust growth near such fronts have been carried out. The small angular separation and emitting areas of such zones, even in the nearest protoplanetary disks, make such observations challenging. ALMA, with its unprecedented spatial resolution and sensitivity, is changing this situation dramatically. Young, actively accreting objects are good initial targets because their luminosities and disk temperature profiles move the condensation fronts out to larger, resolvable, distances (Kennedy & Kenyon 2008; Blake & Bergin 2015).

HL Tau, a young star in the Taurus molecular cloud (∼10⁵ years), was recently observed by ALMA at 2.9, 1.3, and 0.87 mm with a spatial resolution as good as 3 AU (ALMA Partnership et al. 2015, APR henceforth). Here we propose that the dark concentric rings observed in the (sub)millimeter interferometry of HL Tau are due to fast pebble growth near condensation fronts. We first characterize the location and width of surface brightness dips in the ALMA images before comparing their location with the expected condensation fronts in the mid-plane of the disk. We then examine the dust properties in the dips to determine if they are consistent with pebble growth.

The ALMA Band 3, 6, and 7 (or wavelengths of 2.9, 1.3, and 0.87 mm, respectively) HL Tau observations were carried out in 2014, as part of its long-baseline commissioning and science verification (SV) program. With longest baselines of 15.2 km, the HL Tau observations achieved exceptional spatial resolution, with beam sizes of 11.9 × 8.6, 4.9 × 3.1, and 4.2 × 2.7 AU at 2.9, 1.3, and 0.87 mm (for an adopted distance of 140 pc). More detailed descriptions of the data and calibration can be found in APR, who demonstrated the HL Tau continuum images are ∼axisymmetric and well fit with a series of elliptic rings. Therefore, a first order approximation of the HL Tau surface brightness is a single parameter function of the distance from the star.

For our purposes, it is much easier to work on deprojected, circular images. We thus start with the calibrated measurements and use CASA v4.3 (McMullin et al. 2007) to deproject the visibility data and generate synthesized images for all three bands. Assuming the peak continuum emission represents the location of the central star, we first shift the phase center of each band to that of HL Tau (by ${{{\Delta }}_{\alpha }}=-14\;{\rm mas},{{{\Delta }}_{\delta }}=196\;{\rm mas}$). We then deproject the visibility data using the best-fit inclination angle i = 4672 and PA = 13802 from APR.

From the deprojected visibility data in bands 6 and 7 we have also generated a spectral index (α, where ${{I}_{\nu }}\propto {{\nu }^{\alpha }}$) map. This result is computed using the CLEAN task with nterms = 2, based on the multi-frequency deconvolution algorithm developed by Rau & Cornwell (2011).

We measure the surface brightness versus radial distance in steps of 0.5 AU and present the azimuthally averaged radial brightness distributions in Figure 1. The most consistent features in all of the bands are three emission depressions, or dips, around ∼13, 32, and 63 AU. We note that there are many subtle bumps in the surface brightness distribution, and we refer interested readers to a more detailed description by APR. Here we focus on the three most prominent, ∼radially symmetric dips.

Zoom In Zoom Out Reset image size

In order to constrain the dip widths and depths, a smooth fit to the surface brightness distribution is needed. We use a simple vertically isothermal model,

$$\begin{eqnarray}&&{{I}_{\nu }}(r)={\rm cos} i\times {{B}_{\nu }}({{T}_{r}})(1-{{e}^{-{\Sigma }(r){{\kappa }_{\nu }}\,{\rm sec} i}})\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{T}_{r}}={{T}_{0}}{{(r/[1{\rm AU}])}^{q}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\Sigma }(r)={{{\Sigma }}_{c}}{{(r/{{r}_{c}})}^{-\gamma }}{\rm exp} [-{{(r/{{r}_{c}})}^{2-\gamma }}]\end{eqnarray} \tag{ 3 }$$

where ${{I}_{\nu }}$ is the surface brightness (Jy arcsec⁻²), i the inclination angle, ${{\kappa }_{\nu }}$ the continuum opacity in cm² g⁻¹, and ${\Sigma }(r)$ the dust+gas mass surface density distribution for a steady accretion disk (Pringle 1981; Andrews et al. 2009). The cos i factor in Equation (1) accounts for the effects of the deprojection on the synthesized beam. Since our goal is to produce a smooth fit to the surface brightness distribution rather than a realistic physical model, we use ${{\kappa }_{\nu }}$ = 0.01 cm² g⁻¹$(\nu /230[{\rm GHz}])$ throughout the disk (Ossenkopf & Henning 1994).

The results of the best-fit models are plotted as dashed lines in Figure 1. We then divide the observed ${{I}_{\nu }}(r)$ with the best-fit model in each band to constrain the shape of the dips; see Figure 1 (b). The locations are consistent with those of the dark rings reported in APR. The 13 and 32 AU intensity dips appear to be highly symmetric. Gaussian fits to the normalized radial brightness profile are listed in Table 1. The 63 AU dip is asymmetric, so we only list the radius of the minimum and the dip FWHM. The FWHM of the 13 and 32 AU dips are 5–7 AU, about 2× the minor axis of the synthesized beams at 1.3 and 0.87 mm, indicating that the dips are spatially resolved. The 63 AU dip is resolved at all bands, with an FWHM of 17–18 AU. The minimum surface intensity of the dips range from 37% to 78% that of the smooth continuum from the best-fit models.

To check if the dark rings/dips in the ALMA HL Tau images are correlated with condensation fronts, we calculate the condensation temperatures of major ice species under realistic pressure ranges and compare them with the expected disk mid-plane temperature distribution.

Here we consider the major condensible carriers of these elements based on the composition of comets, which are believed to be the best representatives of primordial icy materials in the Solar Nebula (Mumma & Charnley 2011). We include N₂, one of the major carriers of N as predicted by chemical models of disks (Schwarz & Bergin 2014). As a homonuclear diatomic, N₂ cannot be directly measured in cometary comae.

For "pure"(single component) ices we define the condensation temperature of species i as that where the thermal desorption and accretion rates are equal, i.e.,

$$\begin{eqnarray}&&n_{{\rm ice}}^{i}\times k_{{\rm desorp}}^{i}=n_{{\rm gas}}^{i}\times k_{{\rm accr}}^{i},\end{eqnarray} \tag{ 4 }$$

where $k_{{\rm desorp}}^{i}$ and $k_{{\rm accr}}^{i}$ are calculated using the treatment outlined in Woitke et al. (2009) and Walsh et al. (2010).

Besides pure ices, we also consider clathrate hydrates, special forms of crystalline water ice in which gaseous molecules can be trapped inside lattice cages of the ice. At typical mid-plane disk pressures, abundant volatiles such as methane, H₂S, CO and N₂ are expected to be trapped in the form of clathrate hydrates before they condense as pure ice species, given a sufficiently large water ice surface area to which the vapor has access (Lunine & Stevenson 1985). This interesting feature of clathrate hydrates has been employed to explain the low N/O elemental ratio observed in comets as compared to that of the Sun's photosphere (Iro et al. 2003) and the existence of methane in Titan's atmosphere (Lewis 1971).

A summary of the condensation temperatures of pure ice species from such calculations is shown in Table 2, along with the condensation temperatures of clathrate hydrates from Lunine & Stevenson (1985) and Iro et al. (2003). Our condensation temperatures are consistent with the results of Pollack et al. (1991, pp. 469–512) over the same pressure range. Table 2 shows that most of the clathrate hydrates can be created at higher temperatures than their pure ice condensates. The exceptions are NH₃, which has a very similar temperature for its hydrate and pure ice forms, and CO₂, which condenses as a pure frost at higher temperatures than for which the clathrate is stable.

Table 2.  Condensation Temperatures of the Major Volatiles in Disks

Species	T$_{{\rm cond}}^{a}$	E$_{b}$	Cometary Abundance	References
	(K)	(K)	% of H2O
H2O	128–155	5165	100	1, 5
CO	23–28	890	0.4–30	1, 5
CO2	60–72	2605	2–30	1, 5
CH4	26–32	1000	0.4–1.6	2, 5
CH3OH	94–110	4355	0.2–7	1, 5
N2	12–15	520	...	2, 5
NH3	74–86	2965	0.2–1.4	1, 5
HCN	100–120	4170	0.1–0.6	3, 5
H2S	45–52	1800	0.1–0.6	4, 5
NH3·H2O	78–81	...	...	6
H2S$^{\star }$	77–80	...	...	6
CH$_{4}^{\star }$	55–56 (69–72)	...	...	6, 7
CO$^{\star }$	45–46 (58–61)	...	...	6, 7
N$_{2}^{\star }$	41–43 (55–57)	...	...	6

Note. (a). Condensation temperature ranges for ices corresponding to gas number densities of 10¹⁰–10¹³ cm⁻³, suitable for disk mid-planes. * Condensation temperatures of clathrate hydrate formed from hexagonal ice or hydrate under gas number densities of 10¹²–10¹³ cm⁻³. The values in parentheses are for clathrates formed from amorphous ice.

References. (1) Martín-Doménech et al. (2014), (2) Luna et al. (2014), (3) Sandford & Allamandola (1993), (4) Hasegawa & Herbst (1993), (5) Mumma & Charnley (2011), (6) Iro et al. (2003), (7) Lunine & Stevenson (1985).

Download table as:  ASCIITypeset image

In Figure 2 we compare the dip radii with the expected disk condensation front locations. We also plot the observed brightness profile, which provides a rigorous lower bound to the mid-plane dust temperature. This profile would yield a water ice line at <4–5 AU, unresolved by the ALMA SV data, with the NH₃/ammonia-hydrate condensation front near the 12 AU dip.

Zoom In Zoom Out Reset image size

The actual physical temperature in this embedded disk should be larger, however, especially in the case of optically thin dust emission or significant grain coagulation. As a more detailed estimate of the HL Tau dust temperature structure we follow Men'shchikov et al. (1999, MFH henceforth), who used two-dimensional radiative transfer models to quantitatively match the available spectral energy distribution, intensity, and linear polarization data on HL Tau from near IR to (sub)millimeter wavelengths. Their best-fit model results in T$_{{\rm mid}}$ = 665 ${{(r/{\rm AU})}^{-0.6}}$ K for the mid-plane temperature distribution. This is greater than that estimated from the observed brightness distribution in the ALMA image, at all radii. Interestingly, as Figure 2 shows, the dips overlap nicely with the condensation fronts of the most abundant volatiles using the MFH profile: pure water condenses around the 13 AU dip, while pure NH₃ or ammonia (and hydrogen sulfide) hydrates condense around the 32 AU dip. At farther distances, the condensation front of pure CO₂ and the onset of the CO and N₂ clathrate hydrate stability fields, from an amorphous water ice seed occurring near the 63 AU dip.

We have shown above that the most prominent dips in the HL Tau images are remarkably close to the expected condensation fronts of abundant volatiles. We now investigate the dust properties across these radii.

At (sub)millimeter wavelengths, the spectral index α of the surface brightness (where ${{I}_{\nu }}\propto {{\nu }^{\alpha }}$) is a widely used observable to characterize dust properties in protoplanetary disks (Natta et al. 2007, p. 767). For optically thin dust emission, α can be expressed as α = 2 + β, where ${{\kappa }_{\nu }}\propto {{(\nu /{{\nu }_{0}})}^{\beta }}$. The value of β is sensitive to the maximum dust size, ${{a}_{{\rm max} }}$, for a given size distribution, $n(a)\propto {{a}^{q}}$. If the emission is optically thick, however, α = 2 (the blackbody limit) and information on the dust size is lost.

APR provided a cross-cut of α along the major axis of the as-observed image, but it is not clear if this cross-cut is representative of the full disk. In Figure 3, we show the radial distribution of α, measured between 1.3 and 0.87 mm, in cross-cuts at θ = 0°, 45°, 90° and 135° (starting from west) along with the azimuthally averaged $\alpha (r)$. Despite large azimuthal variations, there is a clear general pattern: inside of 48 AU, $\alpha \sim$ 2, except for increases to values of ∼2.6 and ∼2.4 at the 13 AU and 32 AU dips, respectively. Outside 48 AU, α gradually increases from 2 to 2.7 and then becomes relatively flat, except for two shallow decreases around 70 and 80 AU near local maxima in the surface brightness.

Zoom In Zoom Out Reset image size

One way to explain $\alpha \;\sim$ 2 invokes optically thick dust emission in the inner disk, while that in first two dips is optically thin with α = β + 2. This explanation faces two difficulties. First, a massive disk is needed to create optically thick mm-continuum with modest grain growth. Given a typical continuum opacity of ${{\kappa }_{\nu }}$ = 0.01 cm² g⁻¹$(\nu /[230{\rm GHz}])$, a surface density of ${{{\Sigma }}_{{\rm gas}}}$ = 100 g cm⁻² at 48 AU is needed to reach ${{\tau }_{1.3\,{\rm mm}}}$ = 1, a value 20× larger than that for the MMSN (Hayashi 1981). The second issue is the large discrepancy between the temperature derived from spectral energy distribution fits (the MFH model) and the observed brightness temperature (Figure 2). These should be nearly equal if the millimeter-dust emission is optically thick.

An alternative explanation for $\alpha \;\sim$ 2 that mitigates these difficulties is a scenario where the dust emission for $R\;\lesssim$ 48 AU is effectively optically thin, but with $\beta \sim 0$. This can occur when dust grains grow to sizes (a) larger than the observational wavelength (λ). Under such conditions, the absorption cross section simply equals the particle geometric cross section—and is thus independent of wavelength. Such significant dust growth is likely in the dense inner disk, considering that grains beyond 48 AU have grown to $\geqslant$ millimeter sizes (α = 2.77 ± 0.13 (APR)). This dust growth scenario, however, also faces a problem: for a single dust size-distribution, $\alpha \gt 2$ within the 13/32 AU dips demands smaller dust (and reduced mass surface densities) at these radii. It is difficult to understand why the disk inside of and, particularly beyond, the 13 and 32 AU dips should have experienced more extensive dust growth than the dense inner regions of the disk.

We suggest that the apparent contradiction in the dust growth scenario can be solved by a model with two dust populations. This model differs from the commonly invoked power law size distribution but is commonly encountered in numerical simulations of dust growth (e.g., Dullemond & Dominik 2005; Dra̧żkowska & Dullemond 2014). Our first dust population has intermediate growth, while the second population has grown to an average size for which $a\gg \lambda$. Thus,

$$\begin{eqnarray}&&{{\kappa }_{{\rm dust}1,\nu }}={{\kappa }_{1}}{{(\frac{\nu }{{{\nu }_{0}}})}^{{{\beta }_{1}}}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\kappa }_{{\rm dust}2,\nu }}=\frac{\pi {{{\bar{a}}}^{2}}}{4/3\pi \rho {{{\bar{a}}}^{2}}}=\frac{3}{4\rho }\frac{1}{{\bar{a}}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\kappa }_{{\rm dust},\nu }}=(1-f){{\kappa }_{{\rm dust}1,\nu }}+f{{\kappa }_{{\rm dust}2,\nu }}\end{eqnarray} \tag{ 7 }$$

where $\bar{a}$ and f are the average dust size and mass fraction of the second dust population, and ρ is the physical density of the dust.

A plausible order of magnitude estimate can be found if we assume ${{\kappa }_{{\rm dust}1,\nu }}$ = 1 cm² g⁻¹$(\nu /[230{\rm GHz}])$, or dust grain growth to millimeter size. For the second population, we assume ρ = 1 g cm⁻³, that is, ice dominated grains, and ${{\kappa }_{{\rm dust}2,\nu }}$ = 0.75 ${{(\bar{a}/[{\rm cm}])}^{-1}}$ cm² g⁻¹. We further adopt f = 0.9, i.e., 90% of the dust mass in the second dust population. For modest optical depths in the first population,the overall value of α between 1.3 and 0.87 mm changes with $\bar{a}$ (Figure 3(f)). For $\bar{a}$ near 1 cm, the centimeter-sized dust dominates the opacity and leads to $\alpha \sim 2$. When $\bar{a}$ grows to decimeter-size, the mm-sized grains (population 1) dominate the mass specific dust opacity, even though they contain a small fraction of the dust+ice mass, and α becomes significantly greater than 2.

Thus, with significant dust aggregation in the second population, the flux and α behavior in the dips can be explained. With modest mass surface densities, the fractional area of a given column of the disk covered by the dust becomes <unity, and the brightness temperature drops below the physical temperature even though individual dust aggregates are optically thick. Further (sub)millimeter to centimeter high resolution continuum images can distinguish between these two dust/disk structure scenarios through more accurate constraints on the $\alpha (r)$ distribution, as can direct measurements of the vertical and radial gas temperature distributions via molecular emission.

The analysis presented above demonstrates that the location of the three most prominent dips in the HL Tau (sub)millimeter interferometric images are coincident with the expected condensation fronts of the main volatiles, using the mid-plane temperature distribution derived from observations. Further, the spectral index variation inside the dips can be explained by a bimodal dust size distribution model without the need to invoke significant surface mass density depletions.

From our simple order of magnitude model, we infer that most of the dust mass needs to reside in a population that has grown to decimeter-size scales inside the 13 and 32 AU dips. Ros & Johansen (2013) showed that near the water condensation front dust growth from millimeter to at least decimeter-sized pebbles is possible on a timescale of only 1000 years. This rapid growth is consistent with young age of HL Tau, between 0.1 and 1 Myr (Beckwith et al. 1990; Robitaille et al. 2007).

To date, the numerical simulations for dust growth around condensation fronts have considered only water. It is unclear if pebble growth is sufficiently rapid around other condensation fronts to explain the HL Tau results. The extra solid material delivered by clathrate hydrates cannot significantly enhance the solid mass surface density as does H₂O at the water snowline, since each host molecule needs ∼ six water molecules to form the cage structure. The condensation fronts of pure CO₂ and CO are potentially more important for enhancing mass surface density since they account for ∼30% of water abundance in comets. Nevertheless, the model used by Ros & Johansen (2013) can in principle be applied to other condensation fronts. Clathrate hydrate formation, perhaps driven by the transient warming of amorphous ice (Blake et al. 1991) via accretion bursts, can alter the ice rheology. Changes to the sticking efficiency, porosity and compaction of dust+ice aggregates may thus be central to triggering rapid pebble growth in the outer disk.

Once a critical decimeter-sized pebble population is formed, streaming instabilities can drive the creation of >kilometer-sized planetesimals (Johansen et al. 2014, pp. 547–570). If fast pebble growth does preferentially occur around the condensation fronts of abundant volatiles, this would suggest that snow lines regulate the formation and chemical composition of 1–100 km planetesimals, and ultimately the formation and bulk composition of planets.

As the disk evolves the location of condensation fronts will shift inward, perhaps countered by the effects of episodic accretion. How this might affect dust evolution is uncertain (Hubbard & Ebel 2014). Nonetheless, the possibility that disks such as that encircling HL Tau might be seeded with pebbles suggests that planetesimal formation might occur during early evolutionary stages. This is consistent with the cosmochemical record in our solar system which shows that large and differentiated bodies had already begun to form ≲1 Myr after the condensation of calcium–aluminum-rich inclusions, the oldest minerals in the solar system (Qin et al. 2008; Kleine & Touboul 2009; Kruijer et al. 2014).

This paper makes use of the following ALMA data sets: ADS/JAO.ALMA#2011.0.00015.SV. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The authors gratefully acknowledge support provided by the NSF Astronomy & Astrophysics, NSF INSPIRE (AST-1344133), and NASA Origins of Solar Systems grant programs.