Characterization of Ring Substructures in the Protoplanetary Disk of HD 169142 from Multi-Wavelength ALMA Observations

TL;DR

This study uses multi-wavelength ALMA observations to analyze the ring structures in the protoplanetary disk of HD 169142, revealing dust traps likely caused by planet-disk interactions, which are crucial for understanding planet formation.

Contribution

It provides the first unambiguous evidence of radial dust traps in HD 169142's rings through detailed multi-wavelength analysis and modeling.

Findings

Multiple rings are caused by accumulations of large particles at gas pressure bumps.

Maximum dust grain size in rings is about 1 cm, with flatter size distributions than the ISM.

Inner ring may host conditions for streaming instability, relevant for planet formation.

Abstract

We present a detailed multi-wavelength characterization of the multi-ring disk of HD 169142. We report new ALMA observations at 3 mm and analyze them together with archival 0.89 and 1.3 mm data. Our observations resolve three out of the four rings in the disk previously seen in high-resolution ALMA data. A simple parametric model is used to estimate the radial profile of the dust optical depth, temperature, density, and particle size distribution. We find that the multiple ring features of the disk are produced by annular accumulations of large particles, probably associated with gas pressure bumps. Our model indicates that the maximum dust grain size in the rings is $\sim1$ cm, with slightly flatter power-law size distributions than the ISM-like size distribution ($p\sim3.5$) found in the gaps. In particular, the inner ring ($\sim26$ au) is associated with a strong and narrow buildup of dust particles that could harbor the necessary conditions to trigger the streaming instability. According to our analysis, the snowlines of the most important volatiles do not coincide with the observed substructures. We explore different ring formation mechanisms and find that planet-disk interactions are the most likely scenario to explain the main features of HD 169142. Overall, our multi-wavelength analysis provides some of the first unambiguous evidence of the presence of radial dust traps in the rings of HD 169142. A similar analysis in a larger sample of disks could provide key insights on the impact that disk substructures have on the dust evolution and planet formation processes.

Tables

Table 1. Table 1: Summary of ALMA Observations

Project	P.I.	Date	On-source	$N_{\mathrm{ant}}$	Baselines	Freq. range	Flux	Bandpass	Phase
Code			time (min)		(m)	(GHz)	cal.	cal.	cal.
Band 7
https://almascience.nrao.edu/aq/?project_code=2012.1.00799.S (catalog 2012.1.00799.S)	M. Honda	2015 July 26	41.4	41	15–1574	331.040–332.915	J1924-2914	J1924-2914	J1826-2924
						343.050–344.925
		2015 July 27	21.2	41	15–1574	331.040–332.915	J1924-2914	J1924-2914	J1826-2924
						343.050–344.925
		2015 08 Aug	41.4	43	35–1574	331.040–332.915	Pallas	J1924-2914	J1826-2924
						343.050–344.925
Band 6
https://almascience.nrao.edu/aq/?project_code=2015.1.00490.S (catalog 2015.1.00490.S)	M. Honda	2016 September 14	49.6	38	15–3200	232.029–233.904	J1733-1304	J1924-2914	J1820-2528
		2016 September 14	49.6	38	15–3200	232.029–233.904	J1924-2914	J1924-2914	J1820-2528
https://almascience.nrao.edu/aq/?project_code=2015.1.01301.S (catalog 2015.1.01301.S)	J. Hashimoto	2016 September 17	29.3	40	15–3100	231.472–233.472	J1924-2914	J1924-2914	J1820-2528
Band 3
https://almascience.nrao.edu/aq/?project_code=2016.1.01158.S (catalog 2016.1.01158.S)	M. Osorio	2017 May 07	7.8	40	15–1124	89.495–93.495	J1924-2914	J1924-2914	J1826-2924
						101.495–105.495
		2017 September 09	25.4	40	41–7552	89.495–93.495	J1924-2914	J1924-2914	J1826-2924
						101.495–105.495

Table 2. Table 2: Observation Results

Band	Central frequency	Wavelength	Beam shape	rms	Flux DensityaaFlux density integrated within a circle $0{\text{.}}^{\prime \prime }9$ in radius centered on HD 169142.
	(GHz)	(mm)		(mJy beam-1)	(mJy)
7	338.184	0.886	$0{\text{.}}^{\prime \prime }12\times 0{\text{.}}^{\prime \prime }09$, PA$={71}^{\circ }$	0.16	$554\pm 2$
6+7	288.755	1.04	$0{\text{.}}^{\prime \prime }11\times 0{\text{.}}^{\prime \prime }08$, PA$={63}^{\circ }$	0.13	$369\pm 2$
6	225.490	1.33	$0{\text{.}}^{\prime \prime }12\times 0{\text{.}}^{\prime \prime }08$, PA$={71}^{\circ }$	0.08	$166.2\pm 1.3$
3	97.493	3.08	$0{\text{.}}^{\prime \prime }22\times 0{\text{.}}^{\prime \prime }10$, PA$=-{88}^{\circ }$	0.017	$18.29\pm 0.17$

Table 3. Table 3: Model parameters.

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

Equations

$$I_{\nu }(r)=B_{\nu }(T_d(r))(1-e^{-{\tau }_{\nu }(r)}),$$

$${\tau }_{\nu }(r)={\tau }_0(r){\left(\frac{\nu }{{\nu }_0}\right)}^{\beta (r)},$$

$$T_d(r)=T_0{\left(\frac{r}{r_0}\right)}^q,$$

$${\tau }_0^{base}(r)=a^{(base)}{\left(\frac{r}{r_0}\right)}^s,$$

$${\beta }^{base}(r)=b^{(0)}+b^{(base)}{\left(\frac{r}{r_0}\right)}^t,$$

$${\tau }_0^{(i)}(r)=\{\begin{array}{ll}{a}^{(i)}\text{exp}\left[-{\left(\frac{r-x^{(i)}}{\sqrt{2}{\sigma }_{-}^{(\tau ,i)}}\right)}^2\right],& r⩽x^{(i)}\\ a^{(i)}\text{exp}\left[-{\left(\frac{r-x^{(i)}}{\sqrt{2}{\sigma }_{+}^{(\tau ,i)}}\right)}^2\right],& r>x^{(i)}\end{array}$$

$${\beta }^{(i)}(r)=\{\begin{array}{ll}{b}^{(i)}\text{exp}\left[-{\left(\frac{r-x^{(i)}}{\sqrt{2}{\sigma }_{-}^{(\beta ,i)}}\right)}^2\right],& r⩽x^{(i)}\\ b^{(i)}\text{exp}\left[-{\left(\frac{r-x^{(i)}}{\sqrt{2}{\sigma }_{+}^{(\beta ,i)}}\right)}^2\right],& r>x^{(i)}\end{array}$$

$$T(r)={\left(\frac{\phi L_{\star }}{8\pi r^2{\sigma }_{SB}}\right)}^{0.25},$$

Full text

protoplanetary disks — planet-disk interactions — stars: individual (HD 169142) — stars: pre-main sequence — techniques: interferometric

††facilities: VLA, ALMA††software: CASA (v5.3.0; (McMullin et al., 2007)), EMCEE (Foreman-Mackey et al., 2013), Numpy (van der Walt, Colbert, Varoquaux, 2011), Matplotlib (Hunter, 2007), Astropy (Astropy Collaboration et al., 2018)

\AuthorCallLimit

=20

1 Introduction

One of the most important recent discoveries in the planet formation field is the ubiquity of protoplanetary disk substructures (Long et al., 2018; Andrews et al., 2018a). 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; 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, Zhu, & Hartmann, 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, Dullemond, & Pinilla, 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, Dullemond, & Henning, 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-)mm dust continuum emission (e.g., Pérez et al., 2015; Tazzari et al., 2016). The spectral index ($\alpha$) at these wavelengths depends on the optical depth of the emission ($\tau_{\nu}$) and on the power-law index ($\beta$) of the dust opacity (i.e., $\kappa_{\nu}\propto\nu^{\beta}$; Beckwith et al., 1990). If the dust emission is optically thin and in the Rayleigh Jeans regime, the spectral index will simply be $\alpha=2+\beta$. The parameter $\beta$ depends in turn on the size distribution and maximum grain size of the dust: $\beta\sim 1.6-1.8$ is expected for micron-sized grains, while $\beta\sim 0-1$ for dust populations dominated by sizes larger than mm/cm (D’Alessio, Calvet & Hartmann, 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 mm 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\pm 0.8$ pc; Gaia Collaboration et al., 2018; Bailer-Jones et al., 2018) Herbig Ae star HD 169142 (A8 Ve, $M_{\star}\simeq 1.65~{}M_{\odot}$, age$\simeq 10$ Myr; Carney et al., 2018). This star is surrounded by an almost face-on ($i\simeq 13^{\circ}$; 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 $\sim$20 au inner cavity and two bright rings of emission, with radii $\sim 25$ and $\sim 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 (Ligi et al., 2018; Monnier et al., 2017; Pohl et al., 2017; Bertrang et al., 2018) and mm wavelengths (Fedele et al., 2017; Macías et al., 2017), showing that the inner cavity and gap are not only depleted from mm-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 VLT/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 the 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 multi-wavelength analysis of the multi-ring disk of HD 169142 using ALMA observations at 0.89 mm, 1.3 mm, 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 of 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.

2 Observations

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.

\floattable

After inspecting the calibrated visibilities, each dataset was separately self-calibrated with the continuum emission. Phase self-calibration was first performed iteratively decreasing the solution interval in each step, from 120 s to 10 s, or until the peak signal-to-noise ratio (SNR) 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, multi-frequency 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 SNR, 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 datasets, respectively.

In order to combine and compare the observations, all datasets were corrected for their proper motions and shifted to the first epoch of the Band 6 observations (2016 Sep 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 mm, 1.3 mm, and 3 mm by combining all the datasets 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 (Fig. 1) were obtained with uniform weighting for Band 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.

\floattable

Finally, the Band 7 and Band 6 data were combined and imaged together in tclean using mtmfs with the same setup as for the individual images. By combining these two datasets, the cleaning algorithm is capable of producing a map of the spectral index $\alpha$. At the same time, the wider frequency extent of the combined datasets 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 Fig. 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 coverages.

3 Observational Results

The multi-ring morphology of the dust emission of HD 169142 is clearly resolved from 0.89 mm to 3.1 mm. Two ring-like components are detected peaking at radii $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}22$ ($\sim 25$ au) and $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}53$ ($\sim 60$ au), consistent with previous observations (e.g. Osorio et al., 2014; Fedele et al., 2017; Macías et al., 2017; Carney et al., 2018; Bertrang et al., 2018; Pérez et al., 2019). The 0.89 mm 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 mm 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 multi-wavelength analysis of these asymmetries for future higher spatial resolution observations.

The averaged radial intensity profiles of the three images are shown in Fig. 3. These profiles were obtained by averaging the emission within concentric ellipses, projected with the inclination and PA of the disk ($i=13^{\circ}$, PA$=5^{\circ}$; 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^{\circ}$ 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 $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}53$ ($\sim 60$ au) followed by some extended emission that could be associated with an additional ring. The 0.89 mm and 1.3 mm emission display a small decrease of the emission in between these two components at $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}65$ ($\sim 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 $0\hbox to0.0pt{.\hss}^{\prime\prime}503$, $0\hbox to0.0pt{.\hss}^{\prime\prime}563$, and $0\hbox to0.0pt{.\hss}^{\prime\prime}667$ (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 $\mu$m and 1.2 $\mu$m also detected an additional annular ring at $\sim 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.

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 $\sim 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 Fig. 2). However, this combination now allows us to analyze the variations of the 0.89-1.3 mm spectral index ($\alpha$) throughout the disk. The bottom right panel of Fig. 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 ($\sim 10\%$ and $\sim 5\%$ at Band 7 and 6, respectively). This uncertainty would shift the profile upwards or downwards, but its relative shape would be unaffected. As can be seen, $\alpha$ varies significantly throughout the disk. The $\alpha$ profile is relatively flat around $\sim 2.5$ at the position of the inner ring and the outer disk. Even though the spread of $\alpha$ 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 $\beta$, where $\kappa_{\nu}\propto\nu^{\beta}$) than that of large grains, so when its emission is optically thin, it will display larger values of $\alpha$. Dust populations dominated by large dust particles, on the contrary, have much lower values of $\beta$ (e.g D’Alessio, Calvet & Hartmann, 2001; Birnstiel et al., 2018), which results in a lower $\alpha$. 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)mm emission. However, (sub)mm dust emission can also present low values of $\alpha$ if its optical depth ($\tau$) is close to (or larger than) unity. Both effects, low $\beta$ and high $\tau$, are in principle degenerate and indistinguishable with just two wavelengths. Nevertheless, by analyzing the data at 0.89 mm, 1.3 mm, and 3.1 mm, we should be able to pose strict limits on the radial distribution of $\tau$ and $\beta$. We note that the observations at these three bands have different uv coverages 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.

4 Visibility Modeling

In order to analyze our multi-wavelength 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; Pinilla et al., 2018; Macías 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)mm wavelengths, since the emission is dominated by the large dust grains that are settled onto a thin layer in the disk mid-plane –, and that the intensity at each radius $r$ can be calculated as:

[TABLE]

[TABLE]

where $B_{\nu}$ is the Planck function, $T_{d}$ is the dust temperature, $\tau_{0}$ is a reference optical depth at $\nu_{0}=345$ GHz, and $\beta$ is the power of the dust opacity law ($\kappa_{\nu}\propto\nu^{\beta}$). In principle, the three free parameters at each radius are therefore $T_{d}$, $\tau_{0}$, and $\beta$, 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).

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

[TABLE]

where $T_{0}$ is the dust temperature at the reference radius $r_{0}=10$ au. In a flared disk in radiative equilibrium, the mid-plane temperature profile is expected to follow a power-law $T\propto 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 $\tau_{0}$ and $\beta$, 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 $\tau_{0}$ and $\beta$:

[TABLE]

[TABLE]

where we take $r_{0}=10$ au. We note that we include an extra parameter ($b^{(0)}$) for $\beta^{base}$ since $\beta$ is not expected to reach 0, as opposed to $\tau_{0}$. This results in 5 free parameters for the base component. For the rings, we use radially asymmetric Gaussians (i.e., Gaussian rings with independent inner and outer widths):

[TABLE]

[TABLE]

where $i$ denotes each ring; $\sigma_{-}$ and $\sigma_{+}$ are the inner and outer widths, respectively, which can be different for $\tau_{0}$ and $\beta$ (denoted by their superindex); $a$ and $b$ are the peak of the Gaussians for $\tau_{0}$ and $\beta$, respectively; and $x$ is the radius of the ring. We force the rings in $\tau_{0}^{(i)}$ and $\beta^{(i)}$ to be centered at the same position $x^{(i)}$, 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 $\tau_{0}$, 3 per ring for $\beta$, 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 (Fig. 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 $\tau_{0}$ and $\beta$ and include them in eq. (2), which is then used in eq. (1) together with the temperature power-law (eq. (3)) to predict the intensity radial profile at each frequency. The combination of the two type 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 $\tau_{0}$ (i.e., $a^{(i)}>>0$), but no changes in $\beta$ (i.e., the radial profile in $\beta$ being dominated by the base component, with $b^{(i)}\sim 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 $\tau_{0}$ coupled with decreases in $\beta$ at the rings (i.e., $b^{(i)}<<0$). We note that by using $\tau_{0}$ and $\beta$ we avoid making any initial assumptions about the dust composition, which would have been necessary if we aimed at directly using the dust opacity ($\kappa_{\nu}$) in our model.

In addition to the radial profiles of $T_{d}$, $\tau_{0}$, and $\beta$ –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 ($\Delta$RA and $\Delta$DEC). 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}$, $\tau_{0}$, and $\beta$ 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 end up fitting data at 5 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 Band 6 and 3). We also use Gaussian priors for the position offsets based on the expected astrometry uncertainty for ALMA ($\sim 5\%$ of the resolution), and for the inclination and PA of the disk using previous estimates from line observations ($i=13\pm 1^{\circ}$ and PA$=5\pm 5^{\circ}$; 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 $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}01$ ($\sim 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 $\beta$, we set a flat prior that limits $\beta$ to be between 0.2 and 3.1 at all radii. These limits are based on the values of $\beta$ obtained using the dust opacities from Birnstiel et al. (2018), and exploring different power-law size distributions with its power ranging from 1.5 to 4.5, and the maximum grain size from 0.1 $\mu$m to 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 $\tau_{0}$, $\beta$, and $T_{d}$ obtained from our model are shown in Fig. 4, together with the real part of the deprojected visibilities and the predicted intensity profiles at 0.89 mm, 1.3 mm, 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 Fig. 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.

\floattable

We find systematic flux calibration corrections that are within the expected uncertainties. For the phase center shifts we obtain values $\lesssim\pm 20$ mas, also consistent with the expected astrometric uncertainties of ALMA. The inclination and PA of the model with maximum likelihood are $\sim 12.3^{\circ}$ and $\sim 4.8^{\circ}$, respectively, consistent with the values estimated from line observations ($13^{\circ}\pm 1^{\circ}$ and $5^{\circ}\pm 5^{\circ}$; 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, Dominik, & Natta, 2001):

[TABLE]

where $\varphi$ is the flaring angle of the disk (as defined in Chiang & Goldreich, 1997), $L_{\star}$ is the stellar luminosity, and $\sigma_{SB}$ is the Steffan-Boltzmann constant. Assuming $L_{\star}=10~{}L_{\odot}$ (Fedele et al., 2017), our temperature profile would imply $\varphi=0.031$, or a height of the disk photosphere $H_{ph}\sim 0.11r$ at $r=100$ au, consistent with the standard range of values found for Class II objects.

The base component of $\tau_{0}$ 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 $\beta$ is to increase with radius (see Fig. 4), consistently with the expected effects of radial migration. The values of $\beta$ at the innermost radii tend to decrease, probably as a consequence of the emission from the inner disk, which was recently detected at mm wavelengths in high resolution observations at $r\leqslant 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 $\beta$ at $r\lesssim 10$ au. Therefore, we note that the exact behavior of $\beta$ in this region should be taken with caution.

On the other hand, our model finds that the multiple rings in HD 169142 can only be explained with increases in $\tau$ coupled with substantial decreases in $\beta$. The inner ring is centered at $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}2287^{+0\hbox to0.0pt{.\hss}^{\prime\prime}0021}_{-0\hbox to0.0pt{.\hss}^{\prime\prime}0007}$ and displays a sharp and narrow decrease in $\beta$ down to values $\sim 0.2$. This low $\beta$ is coupled with an increase in $\tau_{0}$. The ring is wider in the $\tau_{0}$ profile than in the $\beta$ one (see Fig. 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 $\tau_{0}$ coupled with decreases in $\beta$ that are not as deep as in the inner ring. Our model predicts its second ring (B2/3) at $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}5474^{+0\hbox to0.0pt{.\hss}^{\prime\prime}0014}_{-0\hbox to0.0pt{.\hss}^{\prime\prime}0009}$, while the outermost ring (B4) is centered at $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}631^{+0\hbox to0.0pt{.\hss}^{\prime\prime}003}_{-0\hbox to0.0pt{.\hss}^{\prime\prime}002}$. The latter only shows a very slight decrease in $\beta$, 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 $\tau$ and $\beta$ are smeared out. This effect is minimized by fitting the model in the visibility domain (see §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 ($0\hbox to0.0pt{.\hss}^{\prime\prime}1$ and $0\hbox to0.0pt{.\hss}^{\prime\prime}2$) is larger than the pressure scale height ($\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}03$ 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 multi-wavelength analysis is based on $\tau_{0}$ and $\beta$ rather than on the dust surface density ($\Sigma_{d}$) and the dust opacity ($\kappa_{\nu}$). In this way, we obtain a more direct and model independent estimate of the disk properties. However, an increase of $\tau_{0}$ in our model could be the result of higher $\Sigma_{d}$ and/or higher $\kappa_{0}$. At the same time, lower values of $\beta$ can be obtained with larger maximum grain sizes, or with flatter dust particle size distributions (D’Alessio, Calvet & Hartmann, 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 DSHARP ALMA Large Program (Birnstiel et al., 2018) and assume the typical power-law for the particle size distribution $n(a)\propto a^{-p}$. Using the absorption opacities and Python tools provided by Birnstiel et al. (2018)111Opacities and Python scripts are availabe at https://github.com/birnstiel/dsharp_opac we estimate $\kappa_{0}$ and $\beta_{0.89-3~{}{\rm 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 $\tau_{0}$ and $\beta$, using $\Sigma_{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 $\beta$ and high $\tau_{0}$ 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 10000 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, $\Sigma_{d}$ appears to be most sensitive to $\tau_{0}$, $p$ to $\beta$, and $a_{max}$ to both $\tau_{0}$ and $\beta$.

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, $\Sigma_{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 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 $\Sigma_{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 §4.1), so we advise caution when interpreting the results in these region.

The extremely low values of $\beta$ found at the peak of B1 can only be reproduced with a narrow range of values in the parameter space: $p=1.509^{+0.014}_{-0.006}$, $a_{max}=4.85^{+0.15}_{-0.15}$ mm, and $\Sigma_{d}\sim 0.89^{+0.04}_{-0.04}$ g cm*-2*. In particular, the low value of $p$ clearly points toward the presence of a remarkably strong accumulation of large pebbles. In the outer disk, $\Sigma_{d}$ shows the double peak morphology displayed by B2/3 and B4 in the $\tau_{0}$ profile, with the two peaks in the range $0.16-1.7$ g cm*-2* and $0.063-0.35$ g cm*-2*. At these same radii, $a_{max}$ displays a plateau-like shape around 1 cm, with possible values ranging from $\sim 2$ mm to $\sim 20$ cm, and without showing any significant changes between B2/3 and B4. At the position of B2/3, $p$ decreases slightly down to $\sim 3.2$, but it appears to go up to $\sim 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 $\beta$ 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^{+250}_{-90}$ 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 12CO, 13CO, and C18O line emission, Fedele et al. (2017) estimated a gas mass of $19$ MJ*. Therefore, our estimated dust mass would imply a dust-to-gas mass ratio of $0.03^{+0.04}_{-0.02}$, 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 $\sim 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 underestimate of the total gas mass (e.g., Miotello et al., 2017).

5 Discussion

5.1 Optical Depth

In order to reproduce the multi-wavelength observations of HD 169142, our model predicts a radial profile of $\beta$ and $\tau$ at 345 GHz, with which we can obtain the predicted $\tau$ at the frequencies of our observations. These optical depth profiles are shown in Figure 7.

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 mm 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; Huang et al., 2018a; Dullemond et al., 2018) and it has one important implication: the usual assumption that the (sub)mm 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)mm surveys. At the same time, the (sub)mm spectral index of disks has been used to estimate $\beta$ 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., $\alpha\lesssim 2.5$) could be incorrectly identified in massive disks (see also Ricci et al., 2012).

We note that even when resolved multi-wavelength 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 §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 $\sim 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

As mentioned in §3, the spectral index profile of the disk of HD 169142 displays higher values in the gaps than in the inner ring and outer disk. This behavior could in principle be explained with two different scenarios: the accumulation of larger grains in the rings while smaller grains fill the gaps, and/or optical depths close to 1 in the rings. Previous studies have found similar trends in the spectral index of other disks (e.g. Tsukagoshi et al., 2016; Huang et al., 2018a), but so far very few studies have successfully found unambiguous evidence of accumulations of large dust particles in ring substructures (e.g. Carrasco-González et al., 2016; Liu et al., 2017; Macías et al., 2018). By analyzing our multi-wavelength observations we are able to disentangle both effects and demonstrate that the inner ring and the outer disk of HD 169142 are produced by an increase in $\tau$ coupled with a decrease in $\beta$. This behavior indicates that large dust particles (indicated by a low $\beta$) are being accumulated in the ring substructures, confirming previous results from the modeling of mm observations at individual wavelengths (Osorio et al., 2014; Fedele et al., 2017).

Furthermore, our estimate of the dust surface density and particle size distribution shows low $\Sigma_{d}$ and ISM-like values of $p$ in the gaps (D1 and D2) and at $r>100$ au. These trends are consistent with these regions being almost completely depleted of large dust particles, as also indicated by recent high resolution 1.3 mm observations (Pérez et al., 2019). However, D1 and D2 still harbor some amount of micron-sized dust particles, since recent polarimetric observations show that the gaps are shallower at near-IR wavelengths (Monnier et al., 2017; Pohl et al., 2017; Bertrang et al., 2018). The ISM-like values of $p$ that we find are consistent with this scenario, but since our observations are not sensitive to micron-sized particles, we are unable to constrain the amount of small particles filling the gaps.

On the other hand, our analysis of the particle size distribution also indicates that the low values of $\beta$, coupled with high values of $\tau$, at the inner ring and outer disk can be reproduced with dust populations that have a flatter size distribution than the ISM ($p<3.5$), and have maximum grain sizes between 2 mm and 20 cm. These results strongly indicate that the ring substructures in HD 169142 are the result of accumulations of large dust grains. Such buildups of solids might be associated with increases in growth efficiency beyond certain volatile snowlines, but the more likely scenarios involve the presence of gas pressure bumps that are able to trap the large dust particles (see §5.3). These dust traps have been hypothesized to be an important ingredient in the planet formation process, since they can stop the radial drift of large particles, concentrate them, and enable them to grow up to planetesimal sizes (Pinilla et al., 2012). The discovery of the likely ubiquity of disk substructures has represented a strong support to this theory (Andrews et al., 2018a; Dullemond et al., 2018; Long et al., 2018; van der Marel et al., 2019), but few studies have been able to find unequivocal evidence of their role as dust traps (e.g., Liu et al., 2017; Macías et al., 2018, this work). If similar evidence is found for a larger sample of objects, it could represent the final confirmation that ring substructures represent the solution to the drift and fragmentation barriers in the planetary formation process.

In the case of HD 169142, the extremely low $p$ at B1 (see Fig. 6) is particularly interesting. Such a flat size distribution could only be formed in an extremely tight accumulation of particles that has significantly enhanced the dust growth efficiency, as fragmentation is expected to move $p$ toward $\sim 3.5$ (Birnstiel, Ormel, & Dullemond, 2011). According to our results, the dust surface density at this position reaches $\sim 1$ g cm*-2*, which would imply dust-to-gas mass ratios $\sim 1$ when compared with the gas surface density estimated by Fedele et al. (2017). If confirmed, this inner ring would represent an ideal location to trigger the streaming instability and hence form the seeds for new young planets (Youdin & Goodman, 2005; Auffinger & Laibe, 2018). It is in fact possible that the streaming instability was triggered in the past at this ring, thus resulting in a run-away growth process responsible for the low values of $p$. Interestingly, the predicted $a_{max}$ at this position is not particularly high ($\sim 5$ mm), despite the low value of $p$ (Abod et al., 2018). Nevertheless, we have assumed a certain dust composition, which could affect the results. In addition, our assumption that the dust size distribution can be described by a power-law with a single power might not be appropriate for such accumulations of particles. Lastly, the streaming instability might have induced the formation of small azimuthal clumps in the ring that are not resolved in our data, and where larger particles might be accumulating. Evidence of these possible asymmetries has in fact been revealed by the VLA (Macías et al., 2017) and ALMA (Pérez et al., 2019). A multi-wavelength analysis at higher spatial resolution will be needed to confirm the flat particle size distribution and remarkably strong dust traps at the inner ring of HD 169142.

Finally, we note that our results are derived purely from the dust continuum emission, without including information about the gas. A more complete description of the dust trapping mechanism can be derived by analyzing both components (gas and dust), but such a complex analysis is out of the scope of this paper. For a more detailed modeling of the effects of dust trapping, taking into account the gas density and viscosity, we refer to Sierra et al. (2019).

5.3 Origin of Ring Substructures

Several physical mechanisms have been proposed to explain the presence of ring substructures in (sub)mm 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, Zhu, & Hartmann, 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, Youdin, & Klahr, 2009; Bai & Stone, 2014; Flock et al., 2015; Ruge et al., 2016). Using the results of our multi-wavelength 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 $\tau\sim 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 CO2, CO, N2, and H2O from Zhang, Blake, & Bergin (2015). Their expected locations are plotted over the $\tau_{0}$ and $\beta$ profiles on Figure 8. We note that the N2 snowline would fall at radii $>120$ au, whereas the H2O snowline would be located between $\sim 3$ and $\sim 4.4$ au. Thus, these snowlines are not expected to have any effect on the observed substructures and are not plotted.

The CO2 snowline falls within the inner cavity, while the CO snowline should be located between $\sim 75$ au and $\sim 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 H2O and CO2 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 Fig. 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 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 $\sim 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 a 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 multi-wavelength 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 (Long et al., 2018; Huang et al., 2018b; 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 Effects

The interaction between the magnetic field and the disk can significantly affect its gas dynamics (Bai & Stone, 2017). These magnetohydrodynamic (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 (MRI) turbulence (Johansen, Youdin, & Klahr, 2009; Bai & Stone, 2014) can produce radial changes in the gas density, but with low amplitudes ($\sim 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 mid-plane 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 $\sim 50$ and $\sim 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 $\sim 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 $\sim 0\hbox to0.0pt{.\hss}^{\prime\prime}25$ resolution 1.3 mm and CO(2–1) ALMA observations, Fedele et al. (2017) suggested the presence of a $\lesssim 1$ MJ planet in D1 and a $\sim 1-10$ MJ 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-IRDIS and previous 1.3 mm images with a 3.5 MJ planet in D1, and a second $0.7$ MJ 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-ZIMPOL observations and ALMA archival 0.89 mm data with a 10 MJ planet in D1 and two 1 MJ planets in D2. One of the latter planets (at 35 au) is consistent with a blob recently revealed in SPHERE-IFS 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-)mm rings were produced by radial dust traps on pressure bumps, but there was no robust evidence to support this. Our multi-wavelength 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)mm 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 not related with 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.

6 Summary and Conclusions

We have presented a multi-wavelength 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 mm 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\gtrsim 60$ au) by high angular resolution observations. The spectral index map between 0.89 mm 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 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 ($\tau_{0}$), and for the power of the dust opacity law ($\beta$; $\kappa_{\nu}\propto\nu^{\beta}$). 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 $\sim 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.8}_{-0.3}$ MJ, which would represent a dust-to-gas mass ratio of $0.03^{+0.04}_{-0.02}$ 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 buildups 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.

Multi-wavelength 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 Number 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 (gran 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.