Dust Density Distribution and Imaging Analysis of Different Ice Lines in Protoplanetary Disks

Our models follow the radial evolution of a dust density distribution and calculate the growth, fragmentation, and erosion of dust particles, covering objects from 1 μm to 2 m in size. We do not take into account bouncing of particles. The sticking probability depends on the relative velocities before collision. For this velocity, we take into account Brownian motion, vertical settling, turbulent diffusion, and radial and azimuthal drift. For larger particles, the relative velocities increase (e.g., Windmark et al. 2012). The motion of the particles is determined by the interaction with the gas, and it is calculated according to their size (Birnstiel et al. 2010).

We take parameters of a disk around a Herbig star, specifically ${T}_{\star }=9300\,{\rm{K}}$, ${R}_{\star }=2\,{R}_{\odot }$, and ${M}_{\star }=2.5\,{M}_{\odot }$. We assume that the gas surface density remains constant with time, which is assumed to be a power law with an exponential cutoff (e.g., Lynden-Bell & Pringle 1974), that is,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{g}(r)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-\gamma }\times \exp \left[-{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{2-\gamma }\right].\end{eqnarray} \tag{ 1 }$$

The cutoff radius r_c is taken to be 120 au, and $\gamma =1$. The total disk mass is ${M}_{\mathrm{disk}}=0.08\,{M}_{\odot }$. We assume a simple parameterization for the disk temperature that depends on the stellar parameters and is given by (Kenyon & Hartmann 1987)

$$\begin{eqnarray}T(r) & = & {T}_{\star }{\left(\displaystyle \frac{{R}_{\star }}{r}\right)}^{1/2}{\phi }_{\mathrm{inc}}^{1/4}\\ \mathrm{or}\qquad T(r) & \simeq & 426\,{\rm{K}}{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-1/2},\end{eqnarray} \tag{ 2 }$$

where ${\phi }_{\mathrm{inc}}$ is the angle between the incident radiation and the local disk surface, taken to be ${\phi }_{\mathrm{inc}}=0.05$. We assume a logarithmically spaced radial grid from 1 to 500 au with 500 steps. The initial gas-to-dust ratio is 100, and all the dust particles are assumed to be 1 μm in size at the initial time.

We use simplified prescriptions of the fragmentation velocity threshold, which is assumed to change radially at the location of the ice lines corresponding to H₂O, NH₃, and CO₂. We assume that for grains whose mantle composition is dominated by nonpolar molecules (or with very low dipole moment), such as CO, CO₂, or silicates, the fragmentation velocity is of the order of 1 m s⁻¹, in agreement with results from numerical simulations and laboratory experiments (Blum & Wurm 2000; Poppe et al. 2000; Paszun & Dominik 2009; Gundlach & Blum 2015; Musiolik et al. 2016a, 2016b). Under the assumptions of our models, assuming the CO or CO₂ ice line does not make any difference, except that the location of the CO ice line is farther out, and ${v}_{\mathrm{frag}}$ would be 1 m s⁻¹ for both ice lines. In this paper, we assume the CO₂ ice line, but results would be similar if we were to assume the CO ice line.

In the cases where the grain mantle is composed by molecules with permanent electrical dipoles, such as H₂O and NH₃, the fragmentation velocity is assumed to have a higher value. For H₂O, it is assumed to be 10 m s⁻¹ as suggested by numerical and laboratory experiments (Wada et al. 2009, 2011; Gundlach & Blum 2015). For NH₃, we assume a fragmentation velocity of ∼7 m s⁻¹. This is an estimate, as there are no data on NH₃ fragmentation, but it is in correspondence with the contribution of the dispersion force to the total Van der Waals force between dust grains, which is 24% and 57% for H₂O and NH₃, respectively (Israelachvili 1992). For simplicity we assume the particles to be layered with a silicate core and mantles of water, ammonia, and carbon dioxide depending on their location relative to the ice lines.

At the ice lines of H₂O, NH₃, and CO₂, we assume the fragmentation velocity to change as a smooth step function. We used a smooth step function, given by

$$\begin{eqnarray}H(x)=\left\{\begin{array}{ll}\displaystyle \frac{1}{2}\exp \left(\displaystyle \frac{x}{{\rm{\Delta }}x}\right) & \ \mathrm{for}\ x\leqslant 0\\ 1-\displaystyle \frac{1}{2}\exp \left(-\displaystyle \frac{x}{{\rm{\Delta }}x}\right) & \ \mathrm{for}\ x\gt 0\end{array}\right.,\end{eqnarray} \tag{ 3 }$$

where $x=r-{r}_{\mathrm{ice}}$, with ${r}_{\mathrm{ice}}$ being the radial position of a given ice line. These positions are assumed at the location where the disk temperature has the values of the average freezing temperatures of H₂O, NH₃, and CO₂, in agreement with the values reported in Zhang et al. (2015) (Table 1). Notice that for CO₂, we take a freezing temperature of 44 K, which is an average value for CO₂ and CO. Our main motivation is to have this ice line at a distance of ∼90 au, as observed in the disk around HD 163296 (one of the targets that we use to compare our results with observations in Section 4.3). The factor ${\rm{\Delta }}x$ in Equation (3) is a smoothing parameter for the radial change of the fragmentation velocity at a given position and is taken to be ${\rm{\Delta }}x=0.5\,\mathrm{au}$. Depending on how many ice lines are assumed, the fragmentation velocity is taken to be a composition of different smooth step functions (Equation (3)). For instance, for the simplest case, where only the water ice line is considered (as in Birnstiel et al. 2010), the fragmentation velocity is given by

$$\begin{eqnarray}&&{v}_{\mathrm{frag}}=(100\ast {10}^{H(x)})\displaystyle \frac{\mathrm{cm}}{{\rm{s}}}\quad \mathrm{with}\quad x=r-{r}_{\mathrm{ice},{{\rm{H}}}_{2}{\rm{O}}}.\end{eqnarray} \tag{ 4 }$$

Table 1.  Assumed Freezing Temperatures for H₂O, NH₃, and CO₂.

Parameter	H2O	NH3	CO2
${T}_{\mathrm{cond}}$ (K)	150	80	44
${r}_{\mathrm{ice}}$ (au)	8.1	28.5	92.0

Note. According to the averaged values reported in Zhang et al. (2015) and based on Mumma & Charnley (2011) and Martín-Doménech et al. (2014).

Download table as:  ASCIITypeset image

In this work, we assume three different cases: the fragmentation velocity only changes at H₂O ice line (model I), when the fragmentation velocity changes at the H₂O and CO₂ ice lines, being 1 m s⁻¹ inside the H₂O ice line and beyond the CO₂ ice line (model II), and when H₂O, NH₃, and CO₂ ice lines are assumed for the changes of the fragmentation velocity (model III). Figure 1 shows the profiles of the fragmentation velocity for these three cases. With the disk and stellar properties of our models, the first case recreates the results already presented in Banzatti et al. (2015). In this paper, we add to this case the proper radiative transfer calculations to predict images at scattered-light and millimeter wavelengths.

Zoom In Zoom Out Reset image size

In our models, the ice line locations and gas surface density remain constant with time. It is, however, expected that dust dynamics (in particular radial drift and vertical settling) can change the gas surface density and the location of ice lines (Piso et al. 2015; Cleeves 2016; Krijt et al. 2016; Powell et al. 2017; Stammler et al. 2017). In addition, the disk temperature can vary by different effects, such as disk dispersal, which can also change the ice line locations (Panić & Min 2017). Our models are a simplification assuming that the viscous evolution timescales are longer than the dust growth timescales, and that any material that evaporates from the grains does not contribute significantly to the gas surface density, which is a reasonable assumption for the low dust-to-gas ratios that we have at any time of evolution in our models.

For the collisional outcome, the fragmentation probability (${P}_{\mathrm{frag}}$) is assumed to be unity when the relative velocity (${v}_{\mathrm{rel}}$) is above the fragmentation velocity (${v}_{\mathrm{frag}}$) and zero when ${v}_{\mathrm{rel}}$ is between 0 and 0.8${v}_{\mathrm{frag}}$. For intermediate values (between 0.8 and $1\,{v}_{\mathrm{frag}}$), we assume a linear transition of the fragmentation probability, such that ${P}_{\mathrm{stick}}$ = 1 – P_frag (Birnstiel et al. 2011). When the collision velocity is only determined by turbulence, the maximum grain size is given by (Birnstiel et al. 2012)

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

which is usually known as the fragmentation barrier. In Equation (5), ${\rho }_{s}$ is the volume density of dust grains, which is taken to be $1.2\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, and α is a dimensionless quantity that it is typically used to parameterize the disk viscosity ν as (Shakura & Sunyaev 1973)

$$\begin{eqnarray}&&\nu =\alpha {c}_{s}h\quad {\rm{with}}\quad h=\displaystyle \frac{{c}_{s}}{{\rm{\Omega }}},\end{eqnarray} \tag{ 6 }$$

where ${\rm{\Omega }}=\sqrt{{{GM}}_{\star }{r}^{-3}}$ (G is the gravitation constant) and c_s is the sound speed. Particles can reach large sizes drifting toward the star in timescales shorter than the collision timescales (mainly in the other parts of the disk). This barrier to further growth is the drift barrier, and it is given by

$$\begin{eqnarray}&&{a}_{\mathrm{drift}}=\displaystyle \frac{2{{\rm{\Sigma }}}_{d}}{\pi {\rho }_{s}}\displaystyle \frac{{v}_{K}^{2}}{{c}_{s}^{2}}{\left|\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}r}\right|}^{-1},\end{eqnarray} \tag{ 7 }$$

where ${{\rm{\Sigma }}}_{d}$ is the dust density distribution, v_K is the Keplerian angular velocity (i.e., ${v}_{K}=r{\rm{\Omega }}$), and P is the disk pressure $P=\rho {c}_{s}^{2}$, with ρ being the total gas density.

Okuzumi et al. (2016) modeled the dependency of the fragmentation velocity of dust particles on the sintering. We do not take this particle fusion process into account, which can also affect the fragmentation velocity and create gaps and rings in the dust surface density.

For the radiative transfer calculations, we follow the same procedure presented in Pohl et al. (2016). However, for the near-infrared predictions, we do not calculate images in the full Stokes vector as in Pohl et al. (2016), but concentrate on the total intensity. In addition, we assume the temperature as in the dust evolution models. The models are 3D, but we treat the azimuth with only one grid cell since our dust evolution models are azimuthally symmetric. The total dust density is calculated assuming the dust density distribution from the dust evolution models, such that

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}(R,\varphi ,z)=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}(R)}{\sqrt{2\,\pi }\,{H}_{{\rm{d}}}(R)}\,\exp \left(-\displaystyle \frac{{z}^{2}}{2\,{H}_{{\rm{d}}}^{2}(R)}\right),\end{eqnarray} \tag{ 8 }$$

where R and z are spherical coordinates, that is, $R=r\,\sin (\theta )$ and $z=r\,\cos (\theta )$, with θ being the polar angle. In order to include the effect of settling on the models, the dust scale height ${H}_{{\rm{d}}}(R)$ depends on the grain size and the disk viscosity (Youdin & Lithwick 2007; Birnstiel et al. 2010), and it is given by

$$\begin{eqnarray}&&{H}_{{\rm{d}}}=h\times \min \left(1,\sqrt{\displaystyle \frac{\alpha }{\min (\mathrm{St},1/2)(1+{\mathrm{St}}^{2})}}\,\right),\end{eqnarray} \tag{ 9 }$$

where h is the gas scale height (Equation (6)) and St is the Stokes number, which is a parameter that quantifies the coupling of the particles on to the gas; at the midplane St is given by

$$\begin{eqnarray}&&{\rm{St}}=\displaystyle \frac{a{\rho }_{s}}{{{\rm{\Sigma }}}_{g}}\displaystyle \frac{\pi }{2}.\end{eqnarray} \tag{ 10 }$$

For the synthetic images, we assume a distance to the disk of 140 pc and a face-on disk. For the convolved synthetic images, we assume a Gaussian for the point-spread function (PSF) with an FWHM of 004 for the total intensity at 1.6 μm, as a realistic value for observations with VLT/SPHERE or GPI. For images at the continuum millimeter emission, we consider two wavelengths, 0.87 and 3.0 mm, and we convolve the images with a Gaussian beam of 004, in order to have the same resolution as at short wavelength and to mimic ALMA observations with high resolution.

For the dust optical properties, we assume a dust composition as in Ricci et al. (2010), which is a mixture between silicate, carbonaceous, and water ice. In Banzatti et al. (2015), we investigated how different fractions of water ice at the snow line could affect the resulting dust density distributions and hence the observational predictions. We found a weak effect on the results, and for this reason in this paper we keep the dust composition fixed for all the models.

Figure 2 shows the vertically integrated dust density distribution at 1 Myr of evolution when one, two, or three ice lines are considered (see Figure 1, for reference). Results are shown for three different values of α viscosity, specifically $\alpha =[{10}^{-4},{10}^{-3},{10}^{-2}]$. Figure 3 shows the dust density distribution for small ($a\in [1\mbox{--}10]$ μm) versus large grains ($a\in [1\mbox{--}10]\,\mathrm{mm}$) for the same cases as Figure 2. In this section, we describe the radial variations of the dust density distribution for each value of α viscosity. In the Appendix, the vertical dust density distribution (Equation (8)) is shown for each case. These distributions are assumed for the radiative transfer calculations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.1. Case of α = 10⁻²

3.1.2. Case of α = 10⁻³

3.1.3. Case of α = 10⁻⁴

3.2. Radial Intensity Profiles

After the radiative transfer calculations described in Section 2.2, we obtain images at near-infrared and (sub)millimeter emission, specifically at 1.6 μm, 0.87 mm, and 3 mm. Figure 4 shows the radial profiles of the intensity from the synthetic images. Figure 5 shows the same profiles after convolving with a circular PSF with an FHWM of 004 (which corresponds to 5.6 au for the distance that we assume of 140 pc), in order to compare the multiwavelength profiles at the same spatial resolution. The intensity profiles at 1.6 μm are multiplied by r², to compensate for the ${r}^{-2}$ dependency of the stellar illumination.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The intensity profiles reveal different kinds of structures depending on the assumed value for α viscosity and the number of ice lines that are assumed in the models. In our models the radial structures at scattered-light and millimeter emission originate from variations (radial and vertical; see figures in the Appendix) of the total dust density distribution, and not from temperature variations. The temperature that we assume is a simple power law (Equation (2)), as in the dust evolution models. As a test, we checked whether the obtained intensity profiles change when the temperature is also calculated in the radiative transfer models, and we do not find significant changes. This is because the temperature values are not very different from the power-law assumption, and the only significant variations are in the very inner disk (∼1 au), where the maximum temperature is higher from the radiative transfer calculations.

When only the water ice line is assumed (model I), the general trend is a dip of emission near the location of the water ice line at all different wavelengths. The dip is deeper at 1.6 μm than at (sub)millimeter emission, and for lower viscosity the depression of the dip is higher. These results are similar to those of Birnstiel et al. (2015), who found dips of emission in the outer parts of protoplanetary disks at scattered light due to inefficient fragmentation of dust particles. The main difference with the present models is that the dip is located closer in, near the water ice line, where the fragmentation velocity is expected to change, and not in the outer disk. Although the resulting dust surface density distribution in model I shows a depletion of millimeter grains inside the water ice line and a ring-like structure just beyond, the intensity at (sub)millimeter emission does not show such ring-like structure, and instead it also shows a dip of emission. This is because inside the water ice line there is an enhancement of all grain sizes below the maximum grain size, which is determined by fragmentation ${a}_{\mathrm{frag}}$. All these grains contribute to the emission at (sub)millimeter wavelengths, and when fragmentation becomes less efficient beyond the water ice line, there is a depression in the total intensity because grain growth and drift are more efficient. This is opposite to the case of particle trapping by a planet or at the outer edge of a dead zone, where only the dust grains inside the trap emit at (sub)millimeter wavelengths.

For models II and III, when CO₂ and NH₃ ice lines are included in addition to the water ice line, there are changes in the radial profile of the intensity. The main difference is for the intensity at 1.6 μm, where a clear gap of emission is formed, becoming deeper when α has values of 10⁻³ and 10⁻⁴. At the outer edge of that gap, a ring-like structure is also formed, and in the case of three ice lines and low viscosity, there are multiple ring- and gap-like structures. For the case of $\alpha ={10}^{-2}$, the gap in micron-sized particles is shallower than in the other two cases because when α decreases, the growth is more efficient between the water and the other ice lines. As a result, there is a high depletion of micron-sized particles, creating a deep gap in the intensity profile, whose depletion factor is around four orders of magnitude when comparing to the inner part of the disk. The (sub)millimeter emission in the case of multiple ice lines remains similar to that in the case when only the water ice line is included. This is because the millimeter/centimeter-sized particles are not affected by the changes of the fragmentation velocity when more than one ice line is included in the outer disk. In model III and when $\alpha ={10}^{-2}$, there is a sharp decrease of the intensity at 1.6 μm after the gap and ring located near the water ice line. This is the result of a more flared inner disk with high density of small grains that can block the light from the central star, preventing the light from being scattered in the outer surface layers.

In summary, the general trend in all models is just a dip of emission (or shallow gaps) at (sub)millimeter emission near the ice lines. The width of the gap or the separation between the ring-like shapes is smaller at longer wavelengths. The depression of such a dip or gap at (sub)millimeter is less than one order of magnitude and thus much shallower than at 1.6 μm. It is important to note that at longer wavelengths, that is, 3 mm, the outer disk radius is smaller than at shorter wavelengths as expected from radial drift.

In the convolved images the structures do not show a significant change compared to the theoretical images, and they prove that the structures and depletions that are expected from ice lines can be detected with the current capabilities of different telescopes, such as ALMA, VLT/SPHERE, and GPI.

The assumed value for the disk viscosity (α) plays an important role in the final dust distributions and hence in the potential structures. The value of α depends on how angular momentum is transported within the disk, which can have different origins such as MRI and magnetohydrodynamical (MHD) winds (e.g., Balbus & Hawley 1991; Suzuki & Inutsuka 2009; Bai et al. 2016). MRI requires disk ionization for the disk gas to be coupled to the magnetic field, and therefore if MRI drives accretion and turbulence, the value of α depends on the ionization environment of the disk (e.g., Dolginov & Stepinski 1994; Flock et al. 2012; Desch & Turner 2015).

Observationally, measuring the level of turbulence in protoplanetary disks is quite challenging and remains very uncertain, although recent efforts with ALMA observations have allowed us to measure the turbulent velocity dispersion for a couple of protoplanetary disks (Flaherty et al. 2015, 2017; Teague et al. 2016), revealing values predicted by MRI, but still with large uncertainty ($\alpha \sim {10}^{-4}\mbox{--}{10}^{-2}$; e.g., Simon et al. 2015).

The value of α is of great importance in the context of dust evolution models because it determines the turbulent velocities of the dust particles and hence the fragmentation limit (Equation (5)). After combining our results from the dust evolution models with radiative transfer calculations, we demonstrate that the shape (in particular the depth) of the dips or gaps that formed due to variations of fragmentation dust properties near the ice line depends on α. The depletion factor of the gaps or dips becomes higher for lower values of α, although the main differences are between the results of 10⁻² and the other two values considered (10⁻⁴ and 10⁻³). This is because for $\alpha ={10}^{-2}$, the maximum grain size is determined by fragmentation in the entire disk, while for the other two values it is a mix between fragmentation and drift. Therefore, observational insights about the turbulence in disks can provide better constraints on whether or not the origin of gaps is due to ice lines.

In the context of our simulations, the number of rings or gaps depends on how many changes of the fragmentation velocity are assumed. In this paper, we assume up to three variations corresponding to three main volatiles. It is important to notice that multiple gaps/rings are only obtained at the near-infrared emission, while at the millimeter emission there is only a dip or a single gap. Multiple ring structures at scattered light were only obtained in model III with three radial variations of the fragmentation velocity and low values of α-viscosity (Figures 4 and 5).

There are several possibilities for the origin of rings and gaps in protoplanetary disks, including radial variations of the disk viscosity and planet–disk interaction. However, different observational diagnostics can give insights to distinguish between these scenarios.

Models of planet–disk interaction predict that when the planet is massive enough to open a gap in the gas surface density, there is trapping of millimeter-sized particles at the outer edge of the gap, where the density and hence the pressure increase outward (Rice et al. 2006; Zhu et al. 2011; Pinilla et al. 2012a). Depending on the planet mass (which changes the pressure gradient) and disk viscosity, there is a critical grain size that can be trapped. In general, the micron-sized particles are not effectively trapped, and they can be distributed throughout the disk, showing a smooth distribution or a shallower and smaller gap at short wavelengths than at millimeter emission (de Juan Ovelar et al. 2013; Dong et al. 2015). This is an important difference between our current models of ice lines, where a shallower gap is obtained at the millimeter emission. For instance, there is a dip of emission near the water ice line that is much more depleted at $1.6\,\mu {\rm{m}}$ than at 0.87 or 3.0 mm. In addition, in the case of ice lines, the width of the gap or the separation between the rings is smaller at longer wavelengths.

For the gas distribution, in the case of massive planets, a deep gap is also expected in the total gas surface density, for which depth and width again depend on the planet mass disk viscosity. In the case of ice lines, it is not expected that the total gas surface density has strong variations (e.g., Ciesla & Cuzzi 2006; Stammler et al. 2017). The CO molecular line and its isotopologues are usually used to infer the gas surface density distribution, but these can be strongly affected by the CO ice line (e.g., Schoonenberg & Ormel 2017).

In the case where the planet is not massive enough to open a gap in the gas surface density, the dust can still be shaped in rings and gaps, due to changes in the gas velocities near the planet position, which can slow down the radial dust motion (Paardekooper & Mellema 2004, 2006; Dipierro et al. 2016; Rosotti et al. 2016). As in the case of ice lines, the total gas surface density should not show any strong depletions. Nevertheless, the gaps in the case of low-mass planets are expected to be shallower and narrower for smaller grains than for large grains. As a consequence, measuring the depth and the width of an observed gap at near-infrared emission and millimeter emission is one of the keys to discerning between ice lines and low-mass planetary origin.

Alternatives for gap-opening processes are variations of the disk viscosity that can lead to different transport of angular momentum, creating bumps in the gas surface density that can lead to trapping of particles. The trapping mechanism is more effective for large millimeter grains, while the small grains are expected to be diffused and smoothly disturbed in the disk (Pinilla et al. 2012b), opposite to our current results for the ice lines.

TW Hya and HD 163296.—We discuss our results in the context of current observations of these two disks because they are very well studied in the literature, where the CO ice line has been claimed to be imaged by ALMA. In addition, there are constraints in the level of turbulence in these disks, making them excellent laboratories for discussing our current results.

In TW Hya, the ionization environment was investigated by Cleeves et al. (2015), suggesting a dead zone extending up to a distance of 60 au. The CO ice line was suggested to be at 30 au from observations of N₂H⁺, which might be a chemical tracer of CO ice (e.g., Qi et al. 2013). However, van't Hoff et al. (2017) demonstrated that the amount of CO in the gas phase can affect the abundance of N₂H⁺, such that the N₂H⁺ column density peaks further outside the CO snow line. Assuming disk parameters of the TW Hya disk, van't Hoff et al. (2017) found that the CO ice line should be located at 20 au, in agreement with recent ALMA observations (Zhang et al. 2017).

Near a distance of 20 au there is a distinct gap at scattered light follow by a ring-like emission in the TW Hya disk (see Figure 6; Rapson et al. 2015; van Boekel et al. 2017), and ALMA observations also reveal a gap at 0.87 mm that is narrower than the one observed at scattered light (see Figure 6; Andrews et al. 2016). This gap at 20 au follows the trend found for the gap shape of our models when the disk viscosity is low ($\alpha ={10}^{-3}\mbox{--}{10}^{-4}$; see also Figure 10 in Zhang et al. 2017), in agreement with the possibility that the CO ice line is inside an MRI-dead region. It is important to notice that our models are for a disk around a Herbig star and the CO ice line is much farther out (close to the outer edge) than in TW Hya. However, the observed shapes of the gap near 20 au are similar to the results of our models at ∼8 au, where the fragmentation velocity increases by one order of magnitude. The observed gap, in particular in the millimeter emission, is narrower than our prediction. In our models, the width is determined by where the fragmentation velocity changes and this determines whether the maximum grain size is dictated by fragmentation or drift. Our models may be neglecting variations of the fragmentation velocity due to, for instance, other main volatiles that change the stickiness efficiency near the CO ice line that can make the gap smaller.

Zoom In Zoom Out Reset image size

Additional observational insights can be obtained if we can infer the total gas surface density. If a deep gap is also detected in the total gas surface density at this location, this would suggest that the CO ice line is not responsible of the observed gap.

In HD 163296, the level of turbulence has been constrained to be low (Flaherty et al. 2015), and the CO ice line is located at ∼90 au (Qi et al. 2015). Guidi et al. (2016) suggested a gap followed by a ring feature in the millimeter emission just beyond the CO ice line, where there is also an excess of emission at polarized scattered light (Garufi et al. 2014; Monnier et al. 2017). The most recent ALMA observations from Isella et al. (2016) show a strong ring-like signature near ∼90 au as well (Figure 6). The corresponding gap is deeper at scattered-light emission than at millimeter emission. In this case, the trend of a shallower and narrower gap at millimeter emission than at scattered-light emission is also observed (Figure 6), following our current predictions for the ice lines.

Because these comparisons were done with polarized light intensity, we run as a test a 3D single scattering model for one of the cases (model III with $\alpha ={10}^{-3}$), keeping the same disk parameters to calculate the polarized intensity profile at 1.6 μm. The resulting shape of the radial profile of the polarized intensity is similar to the profile of the total intensity.

Both TW Hya and HD 163296 have multiple rings and gaps, and each might have a different origin. Current observations support that the structures observed near the CO ice line are in good agreement with variations of the dust fragmentation efficiency. The other rings and gaps can originate, for instance, from embedded planets or viscosity variations.

Future observational perspectives.—Detecting the water ice line directly from water lines is challenging because these lines are weak and are at high energy levels (Notsu et al. 2016, 2017). An indirect method to detect ice lines is to observe the effect that these locations have on the grain growth (as shown in Birnstiel et al. 2010), which can lead to strong variations of the emission at millimeter wavelengths and possible variations of the millimeter spectral index (Banzatti et al. 2015).

Cieza et al. (2016) observed with ALMA the protoplanetary disk around V883 Ori, which shows a break in the intensity profile and an optical depth discontinuity at ∼42 au, and they suggested that this break originates owing to the water ice line changing the dust properties and evolution. Because this protostar is experiencing an outburst in luminosity, this increases the disk accretion and temperature, pushing the location of the water ice line farther out than in a typical disk around a T Tauri or Herbig star. Future observations of the intensity at different wavelengths, covering different spectral types, are required to test the robustness of this theoretical prediction.

Observations that provide information about the total gas surface density are complementary because near the ice lines strong depletions of the gas are not expected as in the case of massive planets or variations of the disk viscosity (e.g., at the outer edge of a dead zone). Moreover, the combination of optically thick emission that gives constraints on the disk temperate, together with observations at scattered-light and optically thin (sub)millimeter wavelengths, is an additional and crucial key to recognizing ice lines as potential origins of an observed ring or gap at a given location.