Characterization of Ring Substructures in the Protoplanetary Disk of HD 169142 from Multiwavelength Atacama Large Millimeter/submillimeter Array Observations

One of the most important recent discoveries in the planet formation field is the ubiquity of protoplanetary disk substructures (Andrews et al. 2018a; Long et al. 2018). The implications of this discovery are transformational, but still not fully understood. Since the first discoveries of the dominant presence of disk substructures (e.g., Osorio et al. 2014; Atacama Large Millimeter/submillimeter Array (ALMA) Partnership et al. 2015; Andrews et al. 2016; Pérez et al. 2016; Avenhaus et al. 2018), several possible origins have been proposed. The most likely and promising scenario is the dynamic interaction between the disk and one or more young planets (e.g., Zhu & Stone 2014; Bae et al. 2017; Zhang et al. 2018), yet other processes could still play an important role (e.g., Flock et al. 2015; Okuzumi et al. 2016; Pinilla et al. 2017).

An important piece of information to understand the origin and role of disk substructures is their dust content. Various substructure-forming physical processes involve the onset of gas pressure bumps, which can trap and accumulate large dust particles into annular (e.g., Pinilla et al. 2012) or azimuthally asymmetric structures (e.g., Birnstiel et al. 2013; Lyra & Lin 2013). Such non-smooth gas distributions have been proposed to be a key piece of the dust evolution and planet formation process, since they can stop the radial drift of large particles and allow them to grow up to planetesimal sizes (Whipple 1972; Barge & Sommeria 1995; Brauer et al. 2008). Therefore, analyzing the dust particle size distribution in disks and in their substructures is crucial not only to discern between different substructure origins, but also to understand the role that these features play in the planet formation process.

One of the best methods to study the dust size distribution in disks is by analyzing the spectral behavior of the (sub-) millimeter dust continuum emission (e.g., Pérez et al. 2015; Tazzari et al. 2016). The spectral index (α) at these wavelengths depends on the optical depth of the emission (τ_ν) and on the power-law index (β) of the dust opacity (i.e., κ_ν ∝ ν^β; Beckwith et al. 1990). If the dust emission is optically thin and in the Rayleigh–Jeans regime, the spectral index will simply be α = 2 + β. The parameter β depends in turn on the size distribution and maximum grain size of the dust: β ∼ 1.6–1.8 is expected for micron-sized grains, while β ∼ 0−1 for dust populations dominated by sizes larger than millimeter/centimeter (D'Alessio et al. 2001). Therefore, by analyzing the spectral index of optically thin dust emission, one can in principle study the size distribution of dust particles throughout the disk.

Some recent studies using ALMA observations between 0.88 and 2 mm have started to find evidence of significant radial changes in the spectral indices of disks, with higher values in annular gaps and lower values in the ring substructures (Tsukagoshi et al. 2016; Huang et al. 2018a). These trends could hint at spatial changes in the dust size distribution, but they could also be the result of higher optical depths at the position of the rings. In order to discern between these two effects, observations at multiple wavelengths are necessary, including longer wavelengths where the dust emission will be optically thinner (e.g., Carrasco-González et al. 2016; Liu et al. 2017; Macías et al. 2018).

In this paper we study the protoplanetary disk around the nearby (d = 113.6 ± 0.8 pc; Bailer-Jones et al. 2018; Gaia Collaboration et al. 2018) Herbig Ae star HD 169142 (A8 Ve, M_⋆ ≃ 1.65 M_⊙, age ≃ 10 Myr; Carney et al. 2018). This star is surrounded by an almost face-on (i ≃ 13°; Raman et al. 2006; Panic et al. 2008) pre-transitional disk (Meeus et al. 2010; Honda et al. 2012; Espaillat et al. 2014; Osorio et al. 2014) that has been extensively studied at multiple wavelengths. The disk shows a ∼20 au inner cavity and two bright rings of emission, with radii of ∼25 and ∼60 au. This multi-ring morphology was first inferred through near-IR polarimetric observations (Quanz et al. 2013), and later confirmed to be associated with annular rings in the dust surface density of the disk (Osorio et al. 2014). Later studies have analyzed the disk at near-IR (Monnier et al. 2017; Pohl et al. 2017; Bertrang et al. 2018; Ligi et al. 2018) and millimeter wavelengths (Fedele et al. 2017; Macías et al. 2017), showing that the inner cavity and gap are not only depleted from millimeter-sized dust particles, but also from micron-sized grains and gas. Based on these observations, it has been proposed that the disk harbors two or more giant planets that are creating the signature double-ring morphology of HD 169142 (Fedele et al. 2017; Bertrang et al. 2018; Carney et al. 2018). Reggiani et al. (2014) reported the detection of a planet candidate inside the inner cavity through IR imaging, but the confirmation of this source has remained elusive, and it has been suggested that this feature might be instead associated with the inner ring (Biller et al. 2014; Ligi et al. 2018). Multi-epoch Very Large Telescope (VLT)/Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE) observations of HD 169142 have suggested the presence of faint spiral arms in the disk, as well as a planet candidate inside the disk gap (Gratton et al. 2019). More recently, (Pérez et al. 2019) reported high angular resolution ALMA observations that resolved the outer ring into three separate rings. These authors proposed that a mini-Neptune located in the middle of this triple ring could be responsible for the formation of this feature. Overall, HD 169142 shows strong evidence of the presence of multiple young planets, but the origin of the disk substructures has not been fully confirmed yet.

Here we present a multiwavelength analysis of the multi-ring disk of HD 169142 using ALMA observations at 0.89, 1.3, and 3 mm. In Section 2 we report the observations, their calibration, and the imaging details. Section 3 outlines the results from the ALMA cleaned images. In Section 4 we present our modeling approach to the observed visibilities, and estimate the underlying dust size distribution of the disk. Section 5 includes a discussion on the inferred optical depths and dust size distribution, and how this can help us understand the origin of the ring substructures. We summarize and conclude our analysis in Section 6.

We present new ALMA observations obtained at Band 3 (project code: 2016.1.01158.S), as well as archival ALMA data at Band 7 (project code: 2012.1.00799.S), and Band 6 (project codes: 2015.1.00490.S, and 2015.1.01301.S). The Band 7 observations were first reported in Bertrang et al. (2018), but here we present a more detailed reduction of the data that allowed us to obtain a higher sensitivity and angular resolution. Details of the observations and setups are summarized in Table 1. The raw data were calibrated using the reduction scripts provided in the ALMA archive and their corresponding version of CASA (Common Astronomy Software Applications). Further calibration and imaging was performed in CASA version 5.3.0.

After inspecting the calibrated visibilities, each data set was separately self-calibrated with the continuum emission. Phase self-calibration was first performed iteratively decreasing the solution interval in each step, from 120–10 s, or until the peak signal-to-noise ratio (S/N) did not improve from the previous iteration. The hogbom algorithm within the tclean task was used for the imaging during the phase-only self-calibration, together with natural weighting. Afterwards, amplitude self-calibration was performed using the multi-term, multifrequency synthesis algorithm (mtmfs; Rau & Cornwell 2011) in tclean, assuming a linear spectrum (nterms = 2) and point-like components (scales = 0), as well as using natural weighting. Various iterations of amplitude self-calibration were attempted, performing at least one first iteration with a solution interval as long as the scan. In general, self-calibration resulted in substantial improvements of the peak S/N, typically by factors between 5 and 10. The only exception was the Band 3 data, which due to the much lower brightness temperature of the dust emission, resulted in improvements of 20% and 10% for the compact and extended data sets, respectively.

In order to combine and compare the observations, all data sets were corrected for their proper motions and shifted to the first epoch of the Band 6 observations (2016 September 14). To do this, the CASA tasks fixvis and fixplanets were used, in combination with the position and proper motions measured by Gaia (Gaia Collaboration et al. 2018). Images were then obtained at 0.89, 1.3, and 3 mm by combining all the data sets at each band, and using the mtmfs algorithm with nterms = 2 and multiple scales at 0, 1, 3, and 5 times the beam size. Different visibility weightings were tested to reach a compromise between image quality and angular resolution. The final images presented here (Figure 1) were obtained with uniform weighting for Bands 7 and 6, and Briggs weighting (robust = 0.0) for Band 3. The resulting rms sensitivity, angular resolution, and integrated flux density at each band are listed in Table 2.

Zoom In Zoom Out Reset image size

Table 2.  Observation Results

Band	Central Frequency	Wavelength	Beam Shape	rms	Flux Densitya
	(GHz)	(mm)		(mJy beam−1)	(mJy)
7	338.184	0.886	012 × 009, PA = 71°	0.16	554 ± 2
6 + 7	288.755	1.04	011 × 008, PA = 63°	0.13	369 ± 2
6	225.490	1.33	012 × 008, PA = 71°	0.08	166.2 ± 1.3
3	97.493	3.08	022 × 010, PA = −88°	0.017	18.29 ± 0.17

Note.

^aFlux density integrated within a circle 09 in radius centered on HD 169142.

Download table as:  ASCIITypeset image

Finally, the Band 7 and Band 6 data were combined and imaged together in tclean using mtmfs with the same setup as that used for the individual images. By combining these two data sets, the cleaning algorithm is capable of producing a map of the spectral index α. At the same time, the wider frequency extent of the combined data sets provides a larger UV coverage, which results in a higher quality image at an intermediate frequency of 289 GHz. The resulting image is shown in Figure 2. We attempted to obtain a similar image combining the Band 6 and Band 3 data, but the resulting image quality was substantially poorer, probably due to the larger separation in frequency and the significant difference in UV coverage.

Zoom In Zoom Out Reset image size

The multi-ring morphology of the dust emission of HD 169142 is clearly resolved from 0.89–3.1 mm. Two ring-like components are detected peaking at radii ∼022 (∼25 au) and ∼053 (∼60 au), consistent with previous observations (e.g., Osorio et al. 2014; Fedele et al. 2017; Macías et al. 2017; Bertrang et al. 2018; Carney et al. 2018; Pérez et al. 2019). The 0.89 and 1.3 mm images, with a similar spatial resolution, successfully resolve the outer component in the radial direction, but only do it marginally for the inner one. Our 3.1 mm observations, on the other hand, display a very elongated beam that significantly distorts the emission of the inner component.

At the spatial resolution of our observations (Table 2), the disk appears fairly axisymmetric. We note that some small azimuthal asymmetries can be seen in the inner ring at 0.89 and 1.3 mm, which are consistent with previous SPHERE (Bertrang et al. 2018; Gratton et al. 2019), 7 mm VLA (Macías et al. 2017), and very high-resolution 1.3 mm ALMA observations (Pérez et al. 2019). These asymmetries are consistent with tidal interactions with one or more massive planets in the inner and outer gaps of the disk (Bertrang et al. 2018). Nevertheless, they are just marginally resolved in our images, which makes it difficult to analyze them. Therefore, in the following we will focus on the axisymmetric features of HD 169142, and leave the multiwavelength analysis of these asymmetries for future higher spatial resolution observations.

The averaged radial intensity profiles of the three images are shown in Figure 3. These profiles were obtained by averaging the emission within concentric ellipses, projected with the inclination and PA of the disk (i = 13°, PA = 5°; Raman et al. 2006). In order to minimize the effects of the elongated beam at 3.1 mm, the profile at this band is calculated using only two 90° cones to the north and south, along the minor axis of the beam. The profiles clearly show that the outer disk is resolved into a bright ring at ∼053 (∼60 au) followed by some extended emission that could be associated with an additional ring. The 0.89 and 1.3 mm emissions display a small decrease of the emission in between these two components at ∼065 (∼74 au). Recent high spatial resolution observations at 1.3 mm have in fact resolved the outer component of the disk into three separate rings at 0503, 0563, and 0667 (Pérez et al. 2019). The small dip in our profiles is associated with the gap between the last two rings in those observations, but the second and third rings remain unresolved in our images and appear as one single ring at an intermediate radius. VLT/SPHERE polarimetric observations at 0.73 and 1.2 μm also detected an additional annular ring at ∼90 au (Pohl et al. 2017; Bertrang et al. 2018), which may also be seen at 7 mm (Macías et al. 2017). However, our ALMA observations do not detect any emission at this radius.

Zoom In Zoom Out Reset image size

In order to avoid confusion, in the following we will follow the notation by Pérez et al. (2019), and we will refer to the disk substructures as "B" (for bright annular substructures, a.k.a. rings) or "D" (for dark annular substructures, a.k.a. gaps) followed by their rank in radial distance. In this way, we refer to the inner ring as B1, and to the substructures in the outer disk as B2, B3, and B4, which are all separated between each other by D1 (the inner cavity within B1), D2 (the gap between B1 and B2), D3 (between B2 and B3), and D4 (between B3 and B4). Since we do not resolve B2 from B3, we will refer to the combined bright ring at ∼60 au as B2/3.

The image resulting from the combination of Band 6 and 7 data shows a similar morphology to the images in each band separately (see Figure 2). However, this combination now allows us to analyze the variations of the 0.89–1.3 mm spectral index (α) throughout the disk. The bottom-right panel of Figure 2 shows the radial profile of the spectral index map between these two wavelengths. Due to the relatively high spread of values in some zones of the disk, we compute this profile with the median within each projected ellipse. The colored region in this case indicates the 16th and 84th percentiles. We note that we do not include in this figure the systematic error introduced by the absolute flux calibration at each band (∼10% and ∼5% at Bands 7 and 6, respectively). This uncertainty would shift the profile upwards or downwards, but its relative shape would be unaffected. As can be seen, α varies significantly throughout the disk. The α profile is relatively flat around ∼2.5 at the position of the inner ring and the outer disk. Even though the spread of α is higher in the gaps, the profile indicates that its value is higher in these regions.

This behavior has been found in other protoplanetary disks (e.g., Tsukagoshi et al. 2016; Huang et al. 2018a), and it might indicate substantial changes in the size distribution of dust particles in the disk. The opacity law of small dust particles is usually steeper (i.e., with higher β, where κ_ν ∝ ν^β) than that of large grains, so when its emission is optically thin, it will display larger values of α. Dust populations dominated by large dust particles, on the contrary, have much lower values of β (e.g., D'Alessio et al. 2001; Birnstiel et al. 2018), which results in a lower α. Therefore, the spectral index profile of HD 169142 could be explained if small dust particles are filling the gaps, while large dust grains are accumulating in the ring substructures and are dominating their (sub) millimeter emission. However, (sub) millimeter dust emission can also present low values of α if its optical depth (τ) is close to (or larger than) unity. Both effects, low β and high τ, are in principle degenerate and indistinguishable with just two wavelengths. Nevertheless, by analyzing the data at 0.89, 1.3, and 3.1 mm, we should be able to pose strict limits on the radial distribution of τ and β. We note that the observations at these three bands have different UV coverage and spatial resolutions, which could affect the results if the analysis was performed in the image domain. In order to avoid these issues, while also extracting as much information as possible from the different spatial scales probed by the interferometer, we perform our analysis in the visibility domain.

In order to analyze our multiwavelength observations of HD 169142, we aim at fitting an analytical model of the radial intensity profile of the disk to the real part of the observed deprojected visibilities. Following studies of other protoplanetary disks (e.g., Zhang et al. 2016; Macías et al. 2018; Pinilla et al. 2018), we assume that the disk emission is axisymmetric, so that the visibility profile of the disk can be computed as the Hankel transform of its radial intensity profile (Pearson 1999). We can then compare the model and observed visibilities, and explore the parameter space following a Bayesian approach with a Markov chain Monte Carlo (MCMC) algorithm. We note that the small azimuthal asymmetries in the inner ring (B1) are just marginally resolved, so we do not expect them to have important effects at the spatial scales traced by our observations. Therefore, an axisymmetric disk is a good approximation and should give us a good estimate of the average radial structure of the disk.

Instead of fitting a surface brightness profile at each band separately, we fit all the observed visibilities at the same time. To do this, we assume that the disk is geometrically thin and vertically isothermal—a good approximation at (sub) millimeter wavelengths, since the emission is dominated by the large dust grains that are settled onto a thin layer in the disk midplane—and that the intensity at each radius r can be calculated as:

$$\begin{eqnarray}&&{I}_{\nu }(r)={B}_{\nu }({T}_{d}(r))\,(1-{e}^{-{\tau }_{\nu }(r)}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\tau }_{\nu }(r)={\tau }_{0}(r){\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{\beta (r)},\end{eqnarray} \tag{ 2 }$$

where B_ν is the Planck function, T_d is the dust temperature, τ₀ is a reference optical depth at ν₀ = 345 GHz, and β is the power of the dust opacity law (κ_ν ∝ ν^β). In principle, the three free parameters at each radius are therefore T_d, τ₀, and β, for which we will presume a certain functionality. We note that we have made the common assumption that the main source of opacity is the absorption opacity, without considering the scattering component. For a discussion of the possible effects of the latter component, see Sierra et al. (2019).

First, we assume that the dust temperature follows a power-law profile:

$$\begin{eqnarray}&&{T}_{d}(r)={T}_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{q},\end{eqnarray} \tag{ 3 }$$

where T₀ is the dust temperature at the reference radius r₀ = 10 au. In a flared disk in radiative equilibrium, the midplane temperature profile is expected to follow a power law T ∝ r^−0.5. Instead of fixing q to −0.5, we let this parameter vary in our model to account for possible deviations from this simplified view.

For τ₀ and β, we aim at choosing functions that can be flexible enough to reproduce our observations without overfitting them. We do this by assuming that the disk structure is composed by a base extended component and a set of rings where higher densities and/or smaller or larger dust particles may be concentrated. The base component is modeled as power laws for both τ₀ and β:

$$\begin{eqnarray}&&{\tau }_{0}^{\mathrm{base}}(r)={a}^{(\mathrm{base})}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{s},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\beta }^{\mathrm{base}}(r)={b}^{(0)}+{b}^{(\mathrm{base})}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{t},\end{eqnarray} \tag{ 5 }$$

where we take r₀ = 10 au. We note that we include an extra parameter (b⁽⁰⁾) for β^base since β is not expected to reach 0, as opposed to τ₀. This results in five free parameters for the base component. For the rings, we use radially asymmetric Gaussians (i.e., Gaussian rings with independent inner and outer widths):

$$\begin{eqnarray}{\tau }_{0}^{(i)}(r)=\left\{\begin{array}{cc}{a}^{(i)}\,\exp \ \left[-{\left(\displaystyle \frac{r-{x}^{(i)}}{\sqrt{2}{\sigma }_{-}^{(\tau ,i)}}\right)}^{2}\right], & r\leqslant {x}^{(i)}\\ {a}^{(i)}\,\exp \ \left[-{\left(\displaystyle \frac{r-{x}^{(i)}}{\sqrt{2}{\sigma }_{+}^{(\tau ,i)}}\right)}^{2}\right], & r\gt {x}^{(i)}\end{array}\right.\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}{\beta }^{(i)}(r)=\left\{\begin{array}{cc}{b}^{(i)}\,\exp \ \left[-{\left(\displaystyle \frac{r-{x}^{(i)}}{\sqrt{2}{\sigma }_{-}^{(\beta ,i)}}\right)}^{2}\right], & r\leqslant {x}^{(i)}\\ {b}^{(i)}\,\exp \ \left[-{\left(\displaystyle \frac{r-{x}^{(i)}}{\sqrt{2}{\sigma }_{+}^{(\beta ,i)}}\right)}^{2}\right], & r\gt {x}^{(i)}\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where i denotes each ring; σ₋ and σ₊ are the inner and outer widths, respectively, which can be different for τ₀ and β (denoted by their superindex); a and b are the peak of the Gaussians for τ₀ and β, respectively; and x is the radius of the ring. We force the rings in ${\tau }_{0}^{(i)}$ and β⁽ⁱ⁾ to be centered at the same position x⁽ⁱ⁾, but we let the other parameters (scale, inner width, and outer width) vary as free parameters. This results in 3 free parameters per ring for τ₀, 3 per ring for β, plus 1 parameter per ring for the position of its peak, for a total of 7 free parameters per ring. Based on the shape of the radial profiles obtained from the cleaned images (Figure 3), we use three rings in our model to reproduce B1, B2/3, and B4. Motivated by the recent high spatial resolution ALMA observations that resolved the outer disk into three separate rings (Pérez et al. 2019), we ran an additional model using four Gaussian rings. However, due to the lower spatial resolution of our data, our model with four rings does not reproduce the triple-ring morphology of the outer disk, instead predicting a similar morphology to the case with three rings. Therefore, in the following we use the results of our model with three Gaussian rings.

Finally, we add the base and ring components of τ₀ and β and include them in Equation (2), which is then used in Equation (1) together with the temperature power law (Equation (3)) to predict the intensity radial profile at each frequency. The combination of the two types of components of our model (power law plus three Gaussian rings) is flexible enough to reproduce the observed emission while exploring different origins for the rings. If the rings are not associated with spatial variations of the dust size distribution, we expect the model to reproduce them with substantial increases in τ₀ (i.e., a⁽ⁱ⁾ ≫ 0), but no changes in β (i.e., the radial profile in β being dominated by the base component, with b⁽ⁱ⁾ ∼ 0). On the contrary, if the rings are in fact produced by gas pressure bumps that are trapping and accumulating large dust grains, we expect to have increases in τ₀ coupled with decreases in β at the rings (i.e., b⁽ⁱ⁾ ≪ 0). We note that by using τ₀ and β we avoid making any initial assumptions about the dust composition, which would have been necessary if we aimed at directly using the dust opacity (κ_ν) in our model.

In addition to the radial profiles of T_d, τ₀, and β—which determine the radial intensity profile—we include in our model the orientation of the disk (inclination, i, and position angle of major axis, PA), and potential offsets from the expected position of the star for each observation (ΔR.A. and Δdecl.). Finally, we take into account the potential effects of the systematic flux calibration uncertainty by including the systematic errors in our model. These systematics are, in principle, unconstrainable, but by including them in our MCMC as nuisance parameters, we can explore their effects and include the systematic uncertainty in the posterior of the rest of parameters.

In summary, our model (with three rings) has 51 free parameters (28 parameters for the T_d, τ₀, and β profiles; 23 for disk orientation, phase center shifts, and systematic flux corrections). We separately fit the upper sideband and lower sideband of the Band 7 and Band 3 observations (the Band 6 data only had continuum spectral windows in one of the sidebands), so we ended up fitting data at five different frequencies. In order to reduce the computational time of the model and ensure that each observation has a similar weight during the model fitting, we bin the deprojected visibilities into 2000 bins. We use the MCMC ensemble samplers with affine invariance (Goodman & Weare 2010) in the EMCEE package (Foreman-Mackey et al. 2013). We set Gaussian priors for the systematic errors in the flux calibration based on the nominal errors for ALMA at each band (i.e., 10% at Band 7, and 5% at Bands 6 and 3). We also use Gaussian priors for the position offsets based on the expected astrometry uncertainty for ALMA (∼5% of the resolution), and for the inclination and PA of the disk using previous estimates from line observations (i = 13° ± 1° and PA = 5° ± 5°; Raman et al. 2006). For the rest of the parameters we set flat priors between reasonable values. For the width of the Gaussian rings, we chose the lower limit of the flat priors to be ∼001 (∼1/10 of the beam size), since our observations should not be able to probe smaller spatial scales. This resulted in some chains in the model being limited by the priors, indicating that those structures are partially unresolved. Finally, in order to avoid obtaining non-physical values of β, we set a flat prior that limits β to be between 0.2 and 3.1 at all radii. These limits are based on the values of β obtained using the dust opacities from Birnstiel et al. (2018), and exploring different power-law size distributions with its power ranging from 1.5–4.5, and the maximum grain size from 0.1 μm–1 m. We used 200 walkers and ran 100,000 iterations, which was enough to ensure that, despite the relatively large number of free parameters, the chains of the MCMC had converged.

4.1. Model Results

The radial profiles of τ₀, β, and T_d obtained from our model are shown in Figure 4, together with the real part of the deprojected visibilities and the predicted intensity profiles at 0.89, 1.3, and 3.1 mm. The median and percentiles of the model parameters of these profiles are listed in Table 3. Images of our model and the residuals at these three wavelengths are shown in Figure 5. As can be seen, our model successfully reproduces the observations in great detail, with the residuals only appearing mostly due to the small departures from axisymmetry of the inner ring.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Model Parameters

Parameter	Median, 16%, and 84% Percentiles
Temperature profile
T0 (K)	${77}_{-2}^{+2}$
q	$-{0.501}_{-0.011}^{+0.011}$
Radial profile of τ0
Base component
a(base)	${0.00424}_{-0.00020}^{+0.00016}$
s	$-{0.003}_{-0.006}^{+0.003}$
	Ring 1	Ring 2	Ring 3
a(i)	${0.63}_{-0.02}^{+0.02}$	${0.271}_{-0.014}^{+0.016}$	${0.153}_{-0.007}^{+0.006}$
${\sigma }_{-}^{(\tau ,i)}$	$0\buildrel{\prime\prime}\over{.} {0317}_{-0\buildrel{\prime\prime}\over{.} 0005}^{+0\buildrel{\prime\prime}\over{.} 0010}$	$0\buildrel{\prime\prime}\over{.} {0619}_{-0\buildrel{\prime\prime}\over{.} 008}^{+0\buildrel{\prime\prime}\over{.} 006}$	$0\buildrel{\prime\prime}\over{.} {070}_{-0\buildrel{\prime\prime}\over{.} 0005}^{+0\buildrel{\prime\prime}\over{.} 0007}$
${\sigma }_{+}^{(\tau ,i)}$	$0\buildrel{\prime\prime}\over{.} {0258}_{-0\buildrel{\prime\prime}\over{.} 0017}^{+0\buildrel{\prime\prime}\over{.} 0005}$	$0\buildrel{\prime\prime}\over{.} {0218}_{-0\buildrel{\prime\prime}\over{.} 0017}^{+0\buildrel{\prime\prime}\over{.} 0014}$	$0\buildrel{\prime\prime}\over{.} {067}_{-0\buildrel{\prime\prime}\over{.} 0014}^{+0\buildrel{\prime\prime}\over{.} 0011}$
x(i)	$0\buildrel{\prime\prime}\over{.} {2288}_{-0\buildrel{\prime\prime}\over{.} 0007}^{+0\buildrel{\prime\prime}\over{.} 0021}$	$0\buildrel{\prime\prime}\over{.} {5474}_{-0\buildrel{\prime\prime}\over{.} 0009}^{+0\buildrel{\prime\prime}\over{.} 0014}$	$0\buildrel{\prime\prime}\over{.} {631}_{-0\buildrel{\prime\prime}\over{.} 002}^{+0\buildrel{\prime\prime}\over{.} 003}$
Radial profile of βa
Base component
b(0)	${0.06}_{-0.05}^{+0.30}$
b(base)	${1.00}_{-0.2}^{+0.07}$
t	${0.32}_{-0.12}^{+0.04}$
	Ring 1	Ring 2	Ring 3
b(i)	$-{1.22}_{-0.04}^{+0.06}$	$-{0.92}_{-0.10}^{+0.33}$	$-{0.70}_{-0.13}^{+0.43}$
${\sigma }_{-}^{(\beta ,i)}$	$0\buildrel{\prime\prime}\over{.} {0151}_{-0\buildrel{\prime\prime}\over{.} 0011}^{+0\buildrel{\prime\prime}\over{.} 0013}$	$0\buildrel{\prime\prime}\over{.} {113}_{-0\buildrel{\prime\prime}\over{.} 006}^{+0\buildrel{\prime\prime}\over{.} 005}$	$0\buildrel{\prime\prime}\over{.} {025}_{-0\buildrel{\prime\prime}\over{.} 004}^{+0\buildrel{\prime\prime}\over{.} 004}$
${\sigma }_{+}^{(\beta ,i)}$	$0\buildrel{\prime\prime}\over{.} {0235}_{-0\buildrel{\prime\prime}\over{.} 0042}^{+0\buildrel{\prime\prime}\over{.} 0014}$	$0\buildrel{\prime\prime}\over{.} {040}_{-0\buildrel{\prime\prime}\over{.} 004}^{+0\buildrel{\prime\prime}\over{.} 006}$	$0\buildrel{\prime\prime}\over{.} {95}_{-0\buildrel{\prime\prime}\over{.} 14}^{+0\buildrel{\prime\prime}\over{.} 06}$

Note.

^aThe rings in the β profile have the same radius as the rings in the τ₀ profile.

Download table as:  ASCIITypeset image

We find systematic flux calibration corrections that are within the expected uncertainties. For the phase center shifts we obtain values ≲±20 mas, also consistent with the expected astrometric uncertainties of ALMA. The inclination and PA of the model with maximum likelihood are ∼123 and ∼48, respectively, consistent with the values estimated from line observations (13° ± 1° and 5° ± 5°; Raman et al. 2006).

Our model predicts a temperature profile with a reference temperature at 10 au of ${77}_{-2}^{+2}$ K (median, 16th, and 84th percentiles) and a power $-{0.501}_{-0.011}^{+0.011}$. This profile is consistent with the expected slope for a flared irradiated disk, indicating that the temperature profile of the disk is not significantly affected by the presence of the deep gaps in the disk. We can further compare our temperature profile to the expected profile for such a disk (Chiang & Goldreich 1997; Dullemond et al. 2001):

$$\begin{eqnarray}&&T(r)={\left(\displaystyle \frac{\varphi {L}_{\star }}{8\pi {r}^{2}{\sigma }_{\mathrm{SB}}}\right)}^{0.25},\end{eqnarray} \tag{ 8 }$$

where φ is the flaring angle of the disk (as defined in Chiang & Goldreich 1997), L_⋆ is the stellar luminosity, and σ_SB is the Steffan–Boltzmann constant. Assuming L_⋆ = 10 L_⊙ (Fedele et al. 2017), our temperature profile would imply φ = 0.031, or a height of the disk photosphere H_ph ∼ 0.11r at r = 100 au, consistent with the standard range of values found for Class II objects.

The base component of τ₀ has almost no contribution, and as a consequence the inner cavity (D1) and the gap (D2) are almost completely devoid of emission. Furthermore, the general trend of β is to increase with radius (see Figure 4), consistent with the expected effects of radial migration. The values of β at the innermost radii tend to decrease, probably as a consequence of the emission from the inner disk, which was recently detected at millimeter wavelengths in high-resolution observations at r ≤ 1 au (Pérez et al. 2019). In order to properly account for this component our model would need to include an extra Gaussian centered at r = 0 au, but our observations would not have the necessary resolution to constrain this component. This is in fact shown by the high uncertainty in β at r ≲ 10 au. Therefore, we note that the exact behavior of β in this region should be viewed with caution.

On the other hand, our model finds that the multiple rings in HD 169142 can only be explained with increases in τ coupled with substantial decreases in β. The inner ring is centered at $\sim 0\buildrel{\prime\prime}\over{.} {2287}_{-0\buildrel{\prime\prime}\over{.} 0007}^{+0\buildrel{\prime\prime}\over{.} 0021}$ and displays a sharp and narrow decrease in β down to values ∼0.2. This low β is coupled with an increase in τ₀. The ring is wider in the τ₀ profile than in the β one (see Figure 4). As a result, the predicted emission from the inner ring is significantly narrower at longer wavelengths. For the outer disk, the two rings show increases in τ₀ coupled with decreases in β that are not as deep as in the inner ring. Our model predicts its second ring (B2/3) at $\sim 0\buildrel{\prime\prime}\over{.} {5474}_{-0\buildrel{\prime\prime}\over{.} 0009}^{+0\buildrel{\prime\prime}\over{.} 0014}$, while the outermost ring (B4) is centered at $\sim 0\buildrel{\prime\prime}\over{.} {631}_{-0\buildrel{\prime\prime}\over{.} 002}^{+0\buildrel{\prime\prime}\over{.} 003}$. The latter only shows a very slight decrease in β, which overall indicates that strength of the accumulations of large particles in HD 169142 decreases with radius.

We note that given the finite resolution of our observations, our model results on τ and β are smeared out. This effect is minimized by fitting the model in the visibility domain (see Section 5.1), but our resolution still limits the spatial scales of the substructures that we are able to characterize. In particular, the angular resolution of our data (01 and 02) is larger than the pressure scale height (∼003 at 50 au, from our temperature profile), which is the expected size of disk substructures associated with gas perturbations (e.g., Dong et al. 2017).

4.1.1. Dust Surface Density and Particle Size Distribution

We note that our multiwavelength analysis is based on τ₀ and β rather than on the dust surface density (Σ_d) and the dust opacity (κ_ν). In this way, we obtain a more direct and model-independent estimate of the disk properties. However, an increase of τ₀ in our model could be the result of higher Σ_d and/or higher κ₀. At the same time, lower values of β can be obtained with larger maximum grain sizes, or with flatter dust particle size distributions (D'Alessio et al. 2001; Birnstiel et al. 2018).

We can further explore the implications of our results in terms of the dust size distribution by assuming a particular dust composition. We use the dust composition employed in the Disk Substructures at High Angular Resolution Project ALMA Large Program (Birnstiel et al. 2018) and assume the typical power law for the particle size distribution n(a) ∝ a^−p. Using the absorption opacities and Python tools provided by Birnstiel et al. (2018)⁸ we estimate κ₀ and β_{0.89–3 mm} while varying the maximum grain size (a_max) and power of the size distribution (p). Thereafter we can try to reproduce our predicted profiles of τ₀ and β, using Σ_d, a_max, and p as the three free parameters at each radius. This approach results in a model that is in principle unconstrained, with three free parameters for two data points per radius, but that can be used to discard some combinations of parameters (i.e., low β and high τ₀ cannot be obtained with micron-sized particles). We explore the parameter space using an MCMC. For each radius we run an MCMC with 50 walkers for 10,000 iterations, with a burn-in phase of 8000 iterations. The median, 16th, and 84th percentiles of the posteriors at each radius are shown in Figure 6. Overall, Σ_d appears to be most sensitive to τ₀, p to β, and a_max to both τ₀ and β.

Zoom In Zoom Out Reset image size

Despite the high uncertainties, we are able to set some constraints on the three parameters, especially at the position of the inner ring and of the outer disk. In D1, D2, and at r > 100 au, Σ_d is the parameter that is more tightly constrained, showing values close to 0. In these same regions of the disk, p is consistent with the interstellar medium (ISM) value of 3.5, although lower values would also be possible. On the other hand, a_max remains completely unconstrained at these radii. Close to r = 0 au our analysis predicts a possible slight increase in Σ_d and a slight decrease in p. These features are associated with the contribution from the inner disk that is likely not well constrained (see Section 4.1), so we advise caution when interpreting the results in these region.

The extremely low values of β found at the peak of B1 can only be reproduced with a narrow range of values in the parameter space: $p={1.509}_{-0.006}^{+0.014}$, ${a}_{\max }={4.85}_{-0.15}^{+0.15}$ mm, and ${{\rm{\Sigma }}}_{d}\sim {0.89}_{-0.04}^{+0.04}$ g cm⁻². In particular, the low value of p clearly points toward the presence of a remarkably strong accumulation of large pebbles. In the outer disk, Σ_d shows the double peak morphology displayed by B2/3 and B4 in the τ₀ profile, with the two peaks in the range 0.16–1.7 g cm⁻² and 0.063–0.35 g cm⁻². At these same radii, a_max displays a plateau-like shape around 1 cm, with possible values ranging from ∼2 mm to ∼20 cm, and without showing any significant changes between B2/3 and B4. At the position of B2/3, p decreases slightly down to ∼3.2, but it appears to go up to ∼3.5 in B4. Overall, this suggests that the rings B2/3 and B4 are produced by accumulations of large dust grains, but at a much lesser degree than in B1. It is important to note that B2/3 is in fact composed of two unresolved rings (Pérez et al. 2019), so it is possible that our estimates of β are underestimated and that narrower and stronger accumulations of large particles are present at B2 and B3 separately.

Finally, by integrating the dust surface density profile over the disk area we estimate a total dust mass of ${160}_{-90}^{+250}$ M_⊕. Despite the underlying uncertainty introduced by assuming a certain dust composition, our dust mass estimate is significantly accurate, since this method accounts not only for radial changes in the optical depth, but also for possible variations in dust opacity that are not usually considered. Based on modeling of ALMA ¹²CO, ¹³CO, and C¹⁸O line emission, Fedele et al. (2017) estimated a gas mass of 19 M_J. Therefore, our estimated dust mass would imply a dust-to-gas mass ratio of ${0.03}_{-0.02}^{+0.04}$, which is consistent with the ISM value of 0.01, but suggests that it could be substantially higher. A dust-to-gas mass ratio higher than 0.01 would be reasonable for a ∼10 Myr old disk such as HD 169142, since central photoevaporation should have had plenty of time to disperse a considerable amount of gas from the disk (Owen et al. 2013). Nevertheless, we note that this high dust-to-gas mass ratio does not include an uncertainty in the gas mass, and the depletion of CO in the disk could result in a significant underestimation of the total gas mass (e.g., Miotello et al. 2017).

5.1. Optical Depth

In order to reproduce the multiwavelength observations of HD 169142, our model predicts a radial profile of β and τ at 345 GHz, with which we can obtain the predicted τ at the frequencies of our observations. These optical depth profiles are shown in Figure 7.

Zoom In Zoom Out Reset image size

We find optical depths that are lower than unity at all radii and bands, but they show values >0.1 in the inner ring and outer disk even at 3 mm. In particular, the optical depth in B1 is close to 1 at 0.89 and 1.3 mm. This result is consistent with the optical depths found in ring substructures of other protoplanetary disks (e.g., Carrasco-González et al. 2016; Dullemond et al. 2018; Huang et al. 2018a), and it has one important implication: the usual assumption that the (sub) millimeter emission of protoplanetary disks is optically thin could be significantly inaccurate. This assumption has also been called into question in other recent studies (e.g., Tripathi et al. 2017; Andrews et al. 2018b), and it could have important effects. The dust emission is usually assumed to be optically thin in order to obtain dust masses from (sub) millimeter surveys. At the same time, the (sub) millimeter spectral index of disks has been used to estimate β and hence the level of grain growth in disks. If a significant fraction of the disk is optically thick, these mass estimates could be substantially underestimated, and apparent signs of grain growth (i.e., α ≲ 2.5) could be incorrectly identified in massive disks (see also Ricci et al. 2012).

We note that even when resolved multiwavelength observations are available, angular resolution can play an important role in determining accurate optical depths. If the optically thick regions of the disk are compact—as expected for dust traps in radial pressure bumps (see Section 5.2)—observations without sufficient spatial resolution would smear the emission, resulting in apparent lower optical depths.

In our case, we minimize this effect by directly modeling the observed visibilities instead of the cleaned images. In this way, we use the information of all the spatial scales probed by our observations. Figure 7 shows a comparison of the optical depth profiles of our models (solid lines) and the profiles obtained from the cleaned images (dashed lines). As shown in the figure, using the radial profiles from the cleaned images would have resulted in a significant underestimation of the peak optical depths by a factor ∼3. For this same reason, the peak optical depths estimated by our model are in principle lower limits, and even larger optical depths could be present at smaller spatial scales than the ones probed by our observations and modeling.

5.2. Dust Grain Size Distribution

5.3. Origin of Ring Substructures

Several physical mechanisms have been proposed to explain the presence of ring substructures in (sub) millimeter observations of disks. The most common ones can be roughly classified into three groups: planet–disk interactions (e.g., Papaloizou & Lin 1984; Zhu & Stone 2014; Bae et al. 2017), variations in the disk and/or dust properties at snowlines of volatiles (e.g., Kretke & Lin 2007; Ros & Johansen 2013; Okuzumi et al. 2016; Pinilla et al. 2017), and changes in the gas dynamics/viscosity associated with the magnetic field (e.g., Johansen et al. 2009; Bai & Stone 2014; Flock et al. 2015; Ruge et al. 2016). Using the results of our multiwavelength analysis, we can try to constrain the physical mechanisms responsible for the formation of the ring substructures in HD 169142.

5.3.1. Snowlines

The freezeout or sublimation of gas volatiles on the surface of dust grains can significantly change the fragmentation and sticking properties of the particles (Güttler et al. 2010). Additionally, the freezeout of volatiles near their snowlines can also result in substantial growth of the dust particles, as the amount of solid material increases (Ros & Johansen 2013). As a consequence, some studies predict the formation of annular accumulations of large particles near the snowlines (e.g., Pinilla et al. 2017). On the other hand, other studies have predicted that the sintering of dust grains near snowlines should increase their fragmentation rate, decreasing their size, reducing their radial migration velocity, and hence creating annular buildups of small particles that are able to reach τ ∼ 1 (Okuzumi et al. 2016). As discussed above, our results indicate that the ring substructures in HD 169142 are associated with accumulations of large dust grains, implying that, if snowlines are playing a role in HD 169142, they should be increasing the growth efficiency of dust particles.

We can further explore the snowline scenario by directly comparing the temperature profile obtained in our analysis to the expected position of the most important snowlines in the disk. We take the ranges of freezing temperature of CO₂, CO, N₂, and H₂O from Zhang et al. (2015). Their expected locations are plotted over the τ₀ and β profiles on Figure 8. We note that the N₂ snowline would fall at radii >120 au, whereas the H₂O snowline would be located between ∼3 and ∼4.4 au. Thus, these snowlines are not expected to have any effect on the observed substructures and are not plotted.

Zoom In Zoom Out Reset image size

The CO₂ snowline falls within the inner cavity, while the CO snowline should be located between ∼75 and ∼110 au, beyond the B2/3 ring. The position that we predict for the CO snowline is also consistent with estimates based on the DCO⁺ emission (Macías et al. 2017; Carney et al. 2018). The snowline scenario was also tested by Pohl et al. (2017), who modeled VLT/SPHERE scattered light observations and estimated that the H₂O and CO₂ snowlines could be located close to B1 and B2/3. We find a slightly colder temperature profile than the one obtained by these authors, which moves the snowlines of these two volatiles to closer radii, well within the inner cavity (see Figure 8).

Interestingly, given the uncertainties on the exact freezing temperatures, the CO snowline might be consistent with the position of B4, which only showed evidence of a slight increase in the abundance of large dust particles. If confirmed, this could suggest that the outermost ring substructure in HD 169142 is formed through the growth of the small particles filling the outer disk as they move inward, get close to the snowline, and have their growth efficiency enhanced by the condensation of CO on their surface (Ros & Johansen 2013). As dust particles grow, they should also suffer greater drag forces that will result in a faster radial migration. However, some studies suggest that the enhanced surface density of solids at the snowlines could change the disk viscosity and trigger the formation of a gas pressure bump, which could then trap the large dust particles (Kretke & Lin 2007; Bitsch et al. 2014). On the other hand, another shallow ring substructure has been recently detected in near-IR polarimetric observations at ∼90 au, and it has been suggested to be associated with the CO snowline (Bertrang et al. 2018; Carney et al. 2018). These studies probe the small dust particles in the disk atmosphere, so it is possible that, as the disk temperature increases with height, the effects of the CO snowline are seen at larger radii. More observations will be needed to confirm whether B4 and/or the near-IR ring at 90 au are associated with the CO snowline.

In any case, our multiwavelength analysis indicates that at least the most prominent substructures in HD 169142 (B1 and B2/3) are not associated with snowlines. These results support recent studies of larger samples of disks, where no correlation was found between the position of the substructures and the expected position of volatile snowlines (Huang et al. 2018b; Long et al. 2018; van der Marel et al. 2019). We note that these previous studies were based on observations at a single wavelength, so they were forced to assume a certain temperature profile. By combining observations at multiple wavelengths we are able to estimate the temperature profile, obtaining a more stringent constraint on the relationship between the disk substructures and the position of the snowlines in the disk.

5.3.2. Magnetohydrodynamic (MHD) Effects

The interaction between the magnetic field and the disk can significantly affect its gas dynamics (Bai & Stone 2017). These MHD effects can result in the onset of radial pressure bumps that can trap large dust particles and form annular ring substructures. In general, it is difficult to distinguish between MHD effects and other mechanisms associated with pressure bumps such as planet–disk interactions (Ruge et al. 2016), but a few key characteristics of MHD effects can be identified.

The onset of zonal flows due to the magneto rotational instability turbulence (Johansen et al. 2009; Bai & Stone 2014) can produce radial changes in the gas density, but with low amplitudes (∼10%–20%) that are unlikely to induce such strong substructures as the ones in HD 169142 (Simon & Armitage 2014). On the other hand, the transition in disk ionization at the edge of the dead zone, a region of the disk midplane with a lower ionization fraction, can result in a significant change in disk viscosity, and hence, in the formation of an annular pressure bump (Flock et al. 2015; Ruge et al. 2016). Models usually predict this ring to form between ∼50 and ∼80 au, which could be consistent with the position of the outer rings. However, the observed triple-ring morphology (Pérez et al. 2019) implies that at least two of the rings have a different origin. In particular, the narrow and compact rings B2 and B3 would be hard to reconcile with the dead-zone mechanism. More importantly, the predicted dust trap in this scenario is expected to produce azimuthal asymmetries in the form of vortices (Ruge et al. 2016), which is inconsistent with the axisymmetry displayed by the rings in the outer disk.

Overall, MHD effects appear to be unable to explain the ring substructures in HD 169142, similarly to what has been found in other studies (Huang et al. 2018b). In fact, these mechanisms are predicted to be more effective in younger disks, since the magnetic field in a ∼10 Myr source such as HD 169142 should have been mostly removed (Bai & Stone 2017). However, there are still several uncertainties associated with the MHD effects, such as the magnitude and evolution of the magnetic field in disks (e.g., Bai & Stone 2017), and/or the amplitude of some of theses perturbations (Simon & Armitage 2014).

5.3.3. Planet–Disk Interactions

The most commonly proposed scenario to explain the ring substructures in HD 169142 is the gravitational interaction between the disk and two or more giant planets (e.g., Osorio et al. 2014; Reggiani et al. 2014). Based on ∼025 resolution 1.3 mm and CO(2–1) ALMA observations, Fedele et al. (2017) suggested the presence of a ≲1 M_J planet in D1 and a ∼1–10 M_J planet in D2, similar to what other studies found based on the gap profiles (Kanagawa et al. 2015; Dong et al. 2017). Using dust evolution models, Pohl et al. (2017) reproduced the ring positions in their VLT/SPHERE-Infra-Red Dual-beam Imager and Spectrograph, and previous 1.3 mm images with a 3.5 M_J planet in D1, and a second 0.7 M_J planet in D2, but these authors were unable to fit the shallow depth of the D2 gap in scattered light. On the other hand, Bertrang et al. (2018) performed hydrodynamical simulations and successfully reproduced new SPHERE Zurich Imaging Polarimeter observations and ALMA archival 0.89 mm data with a 10 M_J planet in D1 and two 1 M_J planets in D2. One of the latter planets (at 35 au) is consistent with a blob recently revealed in SPHERE Integral Field Spectrograph observations (Gratton et al. 2019). These observations were not sensitive at the radial position of the second planet in D2.

All these studies were based on the assumption that the (sub-) millimeter rings were produced by radial dust traps on pressure bumps, but there was no robust evidence to support this. Our multiwavelength analysis indicates that B1 and B2/3 are associated with annular accumulations of large dust particles, strongly supporting that the disk of HD 169142 harbors multiple giant planets that are disrupting the disk.

Furthermore, Pérez et al. (2019) recently reported the triple-ring morphology of the outer disk of HD 169142 and proposed that a 10 M_⊕ planet at the position of B3 could be responsible for the formation of B2, B3, and B4. This mini-Neptune would create three annular pressure bumps (one at the radius of its orbit, one at shorter radii, and one at longer radii; Dong et al. 2017, 2018) that would in turn trap the large dust particles of the disk forming three rings at (sub) millimeter wavelengths (Pérez et al. 2019). Even though we are unable to resolve the rings B2 and B3, this scenario could be consistent with our analysis, since we see evidence of dust trapping in B2/3 as well as in B4. However, according to the model by Pérez et al. (2019), B4 should be associated with a tight and prominent dust trap, which appears to be inconsistent with our results. In fact, this dust trap also overestimates the emission of B4 in the 1.3 mm observations presented by these authors. Instead, our analysis suggests that B4 is associated with a faint accumulation of particles that might even have an origin unrelated to planets. Overall, higher spatial resolution observations at 2–3 mm will be needed to accurately analyze the dust traps in the disk and discern the architecture of the planet(s) possibly responsible for this triple ring.

We have presented a multiwavelength analysis of the multi-ring protoplanetary disk of HD 169142. We have reported new ALMA observations at 3 mm, which we have analyzed together with archival ALMA data at 0.89 and 1.3 mm. The observations at the three bands clearly resolve the characteristic double-ring morphology of HD 169142, as well as some signs of the small-scale substructure recently revealed in the outer disk (r ≳ 60 au) by high angular resolution observations. The spectral index map between 0.89 and 1.3 mm shows higher values in the gaps, while lower values are found in the inner ring and the outer disk. This behavior has two possible origins that could take place at the same time: either larger particles are being accumulated in the ring substructures, and/or these substructures have optical depths close to 1.

In order to understand the origin of the changes in the spectral index in the disk, we have modeled the observed visibilities using a simple axisymmetric analytical model, which yields a radial profile for the dust temperature (T_d), the reference optical depth at 345 GHz (τ₀), and for the power of the dust opacity law (β; κ_ν ∝ ν^β). From these results we have then estimated the dust surface density and particle size distribution in the disk. The results of our analysis strongly indicate that the ring substructures in HD 169142 are the result of buildups or accumulations of large dust particles. This represents the first unambiguous evidence of the association of these ring substructures with such accumulations in HD 169142. Furthermore, we find evidence of a particularly strong and narrow buildup of large particles in the inner ring of the disk at ∼26 au (B1), where the conditions could be suitable enough to trigger (or have already triggered) the streaming instability. We estimate a total dust mass in the disk of ${0.5}_{-0.3}^{+0.8}$ M_J, which would represent a dust-to-gas mass ratio of ${0.03}_{-0.02}^{+0.04}$ if we assume the gas mass of 0.19 M_J estimated by Fedele et al. (2017), hinting at higher ratios than the usual assumption of 0.01.

We explore different origins for the formation of the annular substructures in HD 169142. Using the results from our model we can discard dust sintering as having an important effect, since it is incompatible with the buildup of large dust particles that we find in the rings. Other mechanisms linked to the snowlines of volatiles might be associated with the outermost ring (B4), which may be located close to the CO snowline. However, other important snowlines do not appear to coincide with any substructure, pointing toward other origins associated with dust trapping at gas pressure bumps. Even though we cannot completely discard MHD effects as the origin of the dust traps in the outer disk, the age of HD 169142, as well as the narrow width and high contrast of the rings, make this scenario unlikely. Overall, our results strongly support the planet origin scenario, in agreement with other recent studies on HD 169142.

Multiwavelength studies represent the most powerful tool to analyze the size distribution of dust particles in the disk and its association with disk substructures. Extending this study to a larger sample of objects will allow us to confirm whether disk substructures can trap large dust particles, and provide a suitable environment for planetesimal formation.

E.M., and C.C.E. acknowledge support from the National Science Foundation under CAREER grant No. AST-1455042 and the Sloan Foundation. M.O., G.A., J.M.T., and J.F.G. acknowledge financial support from the State Agency for Research of the Spanish MCIU through the AYA2017-84390-C2-1-R grant (co-funded by FEDER). M.O., G.A., and J.F.G. acknowledge support from the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). M.F. and G.H.M.B. acknowledge support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement n° 757957).

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00799.S, JAO.ALMA#2015.1.00490.S, JAO.ALMA#2015.1.01301.S, and JAO.ALMA#2016.1.01158.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST 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. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Facilities: VLA - Very Large Array, ALMA. -

Software: CASA (v5.3.0; (McMullin et al. 2007)), EMCEE (Foreman-Mackey et al. 2013), Numpy (Van Der Walt et al. 2011), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration et al. 2018).