<figure><img src="image_1.png" /></figure>

# The Circumstellar Disk HD 169142: Gas, Dust, and Planets Acting in Concert?\*

A. Pohl¹,², M. Benisty³,⁴ ␐, P. Pinilla⁵ ␐, C. Ginski⁶,⁷, J. de Boer⁶ ␐, H. Avenhaus⁸, Th. Henning¹, A. Zurlo⁹,¹⁰,¹¹ ␐, A. Boccaletti¹², J.-C. Augereau⁴, T. Birnstiel¹³ ␐, C. Dominik⁷ ␐, S. Facchini¹⁴ ␐, D. Fedele¹⁵, M. Janson¹,¹⁶ ␐, M. Keppler¹, Q. Kral¹²,¹⁷, M. Langlois⁹,¹⁸, R. Ligi⁹, A.-L. Maire¹ ␐, F. Ménard⁴ ␐, M. Meyer⁸ ␐, C. Pinte⁴ ␐, S. P. Quanz⁸ ␐, J.-F. Sauvage¹⁹, É. Sezestre⁴, T. Stolker⁷ ␐, J. Szulágyi⁸, R. van Boekel¹, G. van der Plas⁴ ␐, M. Villenave⁴, A. Baruffolo²⁰ ␐, P. Baudoz¹², D. Le Mignant⁹, D. Maurel⁴, J. Ramos¹, and L. Weber²¹

¹Max-Planck-Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany; pohl@mpia.de

²Heidelberg University, Institute of Theoretical Astrophysics, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany

³Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UUMI 3386 and Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

⁴Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁵Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

⁶Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands

⁷Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

⁸Institute for Astronomy, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland

⁹Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388, Marseille, France

¹⁰Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile

¹¹Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

¹² LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France

¹³University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 München, Germany

¹⁴Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstrasse 1, D-85748 Garching, Germany

¹⁵INAF—Osservatorio Astrofisico di Arcetri, L.go E. Fermi 5, I-50125 Firenze, Italy

¹⁶Department of Astronomy, Stockholm University, AlbaNova University Center, SE-10691 Stockholm, Sweden

¹⁷Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

¹⁸CRAL, UMR 5574, CNRS, Université Lyon 1, 9 avenue Charles André, F-69561 Saint Genis Laval Cedex, France

¹⁹ONERA, Optics Department, BP 72, F-92322 Chatillon, France

²⁰INAF—Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy

²¹Geneva Observatory, University of Geneva, Chemin des Maillettes 51, 1290 Versoix, Switzerland

Received 2017 July 27; revised 2017 October 10; accepted 2017 October 17; published 2017 November 16

## Abstract

HD 169142 is an excellent target for investigating signs of planet-disk interaction due to previous evidence of gap structures. We perform *J*-band (~1.2 µm) polarized intensity imaging of HD 169142 with VLT/SPHERE. We observe polarized scattered light down to 0.′′16 (~19 au) and find an inner gap with a significantly reduced scattered-light flux. We confirm the previously detected double-ring structure peaking at 0.′′18 (~21 au) and 0.′′56 (~66 au) and marginally detect a faint third gap at 0.′′70–0.′′73 (~82–85 au). We explore dust evolution models in a disk perturbed by two giant planets, as well as models with a parameterized dust size distribution. The dust evolution model is able to reproduce the ring locations and gap widths in polarized intensity but fails to reproduce their depths. However, it gives a good match with the ALMA dust continuum image at 1.3 mm. Models with a parameterized dust size distribution better reproduce the gap depth in scattered light, suggesting that dust filtration at the outer edges of the gaps is less effective. The pileup of millimeter grains in a dust trap and the continuous distribution of small grains throughout the gap likely require more efficient dust fragmentation and dust diffusion in the dust trap. Alternatively, turbulence or charging effects might lead to a reservoir of small grains at the surface layer that is not affected by the dust growth and fragmentation cycle dominating the dense disk midplane. The exploration of models shows that extracting planet properties such as mass from observed gap profiles is highly degenerate.

**Key words:** planet-disk interactions – protoplanetary disks – radiative transfer – scattering – techniques: polarimetric

## 1. Introduction

About two decades ago, our own solar system was the only available laboratory to test models of planet formation. Today, we know that planetary systems are common around other stars and that their architectures are very diverse. The initial conditions and evolution of protoplanetary disks, where planets form, must have a direct influence on most

fundamental properties of their planetary systems (Mordasini et al. 2012, 2016). It is therefore essential to improve our knowledge of the structure of protoplanetary disks by observing and studying them at high spatial scales and with various tracers that enable us to characterize different disk regions. This, indirectly, can constrain the physical processes that influence the disk evolution (e.g., gap opening by a planet, dust growth and settling, photo-evaporation). Even with the advent of a new generation of extreme adaptive optics instruments, the detection of forming planets within their host disks is still challenging, but one can look for indirect

<footnote>* Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO program 095.C-0273.</footnote>
signatures of planet formation, such as the imprints that it
leaves on the disk.

In recent years, high-resolution images of protoplanetary
disks have been published, both in scattered light that traces the
(sub-)micron-sized dust particles in the upper disk layers and in
the (sub)millimeter ((sub-)mm) regime that traces larger dust
grains (mm and centimeter-sized (cm) grains), while the bulk
mass is not directly observable. While for a long time
protoplanetary disks were thought to be smooth and contin-
uous, a variety of small-scale features are now frequently
detected in these images and seem rather common, if not
ubiquitous. Large cavities (~few tens of au) are detected in a
number of objects (e.g., Williams & Cieza 2011), the transition
disks, that often have spectral energy distributions (SEDs) with
a clear dip in the mid-infrared (MIR), indicating a lack of dust
in the inner regions (Strom et al. 1989). Smaller cavities and
gaps in the inner au are also present (Menu et al. 2015) but
cannot be easily directly imaged, nor do they leave a clear
imprint on the SED. Multiple-arm spiral features are observed,
mostly in scattered light (Muto et al. 2012; Garufi et al. 2013;
Grady et al. 2013; Avenhaus et al. 2014; Benisty et al. 2015;
Stolker et al. 2016; Benisty et al. 2017), but also more recently
in the sub-mm wavelength range (CO gas lines: Christiaens
et al. 2014; Tang et al. 2017; continuum emission: Pérez
et al. 2016). These spiral arms have large opening angles, and
their origin is still not fully understood. The presence of one or
more rings and gaps in disks seems to be quite a common
feature; they are found in both young (e.g., HL Tau; ALMA
Partnership et al. 2015; Carrasco-González et al. 2016) and
rather old systems (e.g., TW Hya; Rapson et al. 2015; Andrews
et al. 2016; Tsukagoshi et al. 2016; van Boekel et al. 2017) and
around stars of very different spectral types (e.g., de Boer
et al. 2016; Ginski et al. 2016; van der Plas et al. 2017).
Various mechanisms have been proposed in the literature that
can be assigned to three main categories: structures caused by
fluid dynamics, dust evolution effects, and planet-disk
perturbations. More precisely, these possibilities include zonal
flows from magneto-rotational instability (e.g., Simon &
Armitage 2014; Béthune et al. 2016), gap/bump structures in
the surface density close to the dead-zone outer edge (e.g.,
Flock et al. 2015; Pinilla et al. 2016; Ruge et al. 2016), efficient
particle growth at condensation fronts near ice lines or a
depletion of solid material between ice lines (Zhang et al. 2015;
Pinilla et al. 2017; Stammler et al. 2017), aggregate sintering
zones (Okuzumi et al. 2016), secular gravitational instabilities
(Youdin 2011; Takahashi & Inutsuka 2014), or planet-disk
interactions (e.g., Zhu et al. 2011, 2012; Dong
et al. 2015, 2016; Rosotti et al. 2016). Finally, dips or dark
regions can be interpreted as shadows by inner disk material
(e.g., Marino et al. 2015; Pinilla et al. 2015b; Stolker
et al. 2016; Canovas et al. 2017).

We focus on a multiple-ring system in this study: specifically,
the ~6<sup>+6</sup><sub>-3</sub> Myr old Herbig A5/A8 star HD 169142 (Dunkin
et al. 1997; Grady et al. 2007), located at a distance<sup>22</sup> of
d = 117 ± 4 pc (Gaia Collaboration et al. 2016). With this new
distance of 117 pc, the star is intrinsically less luminous by a
factor of ~0.65. The age estimate by Grady et al. (2007) is based
on Hertzsprung-Russell (HR) placement of the companion
2MASS 18242929-2946559. Moving this star down in the HR
diagram (Figure 9 in Grady et al. 2007) leads to a revised age

estimate of ~10 Myr. Its SED shows a strong infrared excess
indicating a young gas-rich disk with many emission-line features
(Kama et al. 2016; Riviere-Marichalar et al. 2016; Seok &
Li 2017) and a clear dip of emission in the infrared regime
(Grady et al. 2007; Meeus et al. 2010), qualifying it as a transition
disk. The near-infrared (NIR) flux indicates the presence of hot
dust close to the sublimation radius, resolved by NIR interfero-
metric observations (Lazareff et al. 2017). HD 169142 still
experiences gas accretion onto the star, with estimates of the mass
accretion rate varying between 0.7 and 2.7 × 10<sup>-9</sup> M<sub>☉</sub> yr<sup>-1</sup>
(Grady et al. 2007; Wagner et al. 2015). Garufi et al. (2017) noted
that HD 169142 has a reduced NIR excess compared to
continuous Herbig disks or those hosting spirals. The NIR and
MIR fluxes were also found to vary by up to ~45% over a
temporal baseline of 10 yr, indicating strong variability in the
innermost regions (Wagner et al. 2015).

The outer disk has a low inclination (i = 13°, PA = 5°), as
derived by CO mm observations (Raman et al. 2006; Panić
et al. 2008) and confirmed with high-contrast imaging in the
NIR (Quanz et al. 2013; Momose et al. 2015; Monnier
et al. 2017). These images show, from small to larger
separations from the star, a wide inner cavity, a bright
(unresolved) ring, a second wide gap, and an outer disk that
extends up to 1"/7. The inner cavity appears devoid of small
dust grains, while the second gap is not. Observations with the
Atacama Large Millimeter Array (ALMA) at 1.3 mm, obtained
with a resolution of 0.′′28 × 0.′′18, also show two rings
(0.′′17–0.′′28 and 0.′′48–0.′′64) and a gap between them (Fedele
et al. 2017). The mm continuum extends up to 0.′′64, while the
gas extends up to twice as far. The channel maps of the 2–1 line
transition of the three CO isotopologues reveal the presence of
gas inside the dust gaps. Model fitting provides a drop in the
gas surface density by a factor of 30–40. The two rings are also
detected in Very Large Array (VLA) observations at longer
wavelengths (7 mm; Osorio et al. 2014; Macías et al. 2017),
and the azimuthally averaged radial intensity profiles indicate
the marginal detection of a new gap at ~0.′′7, very close to the
CO ice line (Macías et al. 2017). In addition to the disk
features, a candidate massive companion was proposed,
slightly inside the inner ring, at a separation of ~0.′′11 and
~0.′′16, respectively (Biller et al. 2014; Reggiani et al. 2014).
Osorio et al. (2014) reported the detection of a compact 7 mm
emission source with VLA external to the inner ring. The
detection of pointlike structures in the context of potential
planetary companions is discussed further in Ligi et al. (2017).

In this paper, we report new polarized differential images of
HD 169142 obtained in the J-band with the Spectro-Polari-
meter High-contrast Exoplanet REsearch instrument (SPHERE;
Beuzit et al. 2008) at the Very Large Telescope (VLT),
complemented with ALMA continuum data from Fedele et al.
(2017). We investigate whether the observed rings can be
explained by the trapping of dust particles as a consequence of
the presence of two planets. This study on HD 169142 serves
as a prototype in which it is demonstrated that multiwavelength
observations are needed to constrain the dust size distribution
and physical mechanisms at work in the disk. The paper is
organized as follows. Section 2 describes the observations and
data processing. Section 3 reports on the detected disk features,
Section 4 provides a physical disk model, and, in Section 5, we
discuss our findings.

<footnote n="22">22 Note that we are using the revised value by Gaia, while most of the papers in the literature use d = 145 pc.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left: J-band azimuthally polarized intensity image $Q_{\phi}$ in logarithmic scale for better visualization. Right: $Q_{\phi} \times r^2$ in linear scale with annotations for the gap and ring structures. Each image pixel is multiplied by the square of its distance to the star, $r^2$, to compensate for the stellar illumination drop-off with radius. All flux scales are normalized to half of the brightest pixel along the inner ring. The region masked by the coronagraph is indicated by the gray circle. North is up, and east is toward the left.</figcaption></figure>

## 2. Observations and Data Reduction

The observations were obtained at the VLT at Cerro Paranal, Chile, on 2015 May 2 with the SPHERE instrument. SPHERE is equipped with an extreme adaptive optics (AO) system (Fusco et al. 2006; Petit et al. 2014; Fusco et al. 2014) that feeds three science channels allowing for high-angular-resolution and high-contrast imaging, spectroimaging, and/or polarimetry at visible and NIR wavelengths. The observations were obtained through the Guaranteed Time program. HD 169142 was observed in the J-band filter ($\lambda_0 = 1.258$, $\Delta\lambda = 0.197$ µm) using the polarimetric imaging mode of the infrared dual-band imager and spectrograph (IRDIS; Dohlen et al. 2008; Langlois et al. 2014), with a plate scale of 12.25 mas pixel$^{-1}$ (Maire et al. 2016) and a 145 mas diameter coronagraphic focal mask (N_ALC_YJ_S, inner working angle (IWA) of 0.08; Boccaletti et al. 2008). HD 169142 was observed for ~53 minutes on-source under moderate AO conditions (seeing of 0.9). The analysis of the reference point-spread function (PSF) that is estimated from a noncoronagraphic total intensity measurement shows that the observations reach a 33.8 mas × 40.8 mas resolution (FWHM along the x and y directions) and a Strehl ratio of 56%.

We observed HD 169142 using the polarimetric differential imaging technique (PDI; e.g., Kuhn et al. 2001; Apai et al. 2004) that measures the linear polarization of the light scattered by dust grains in the disk and enables one to efficiently remove the unpolarized contribution, including that from the star. This allows us to image, with high contrast, the polarized signal from the disk. In this mode, the instrument splits the beam into two orthogonal polarization states. The control of the polarization orientation was performed with a half-wave plate that was set to four positions shifted by 22.5° in order to construct a set of linear Stokes images. We reduce the data according to the double-difference method (Kuhn et al. 2001) and obtain the Stokes parameters $Q$ and $U$. If we assume that there is only one scattering event for each photon, the scattered light from a circumstellar disk seen at a low inclination angle is expected to be linearly polarized in the azimuthal direction. It is therefore convenient to describe the

polarization vector field in polar coordinates (Schmid et al. 2006; Avenhaus et al. 2014). We therefore define the polar-coordinate Stokes parameters $Q_{\phi}$, $U_{\phi}$ as

$$Q_{\phi} = +Q \cos(2\phi) + U \sin(2\phi), \quad (1)$$

$$U_{\phi} = -Q \sin(2\phi) + U \cos(2\phi), \quad (2)$$

where $\phi$ is the position angle of the location of interest (x, y) with respect to the star location. In this coordinate system, the azimuthally polarized flux appears as a positive signal in the $Q_{\phi}$ image, whereas the $U_{\phi}$ image remains free of disk signal and can be used as an estimate of the residual noise in the $Q_{\phi}$ image (Schmid et al. 2006). This is only valid for disks with face-on geometry, since multiple scattering effects in inclined disks can cause a considerable physical signal in $U_{\phi}$ (e.g., TCha; Pohl et al. 2017). The correction for instrumental polarization is done using a $U_{\phi}$ minimization by subtracting scaled versions of the total intensity frame from the Stokes $Q$ and $U$ frames. The final data images are corrected for the true north (by rotating them by 1.775° in the counterclockwise direction; Maire et al. 2016). We do not attempt to perform an absolute flux calibration of our images due to the inherent problems with measuring flux in PDI images.

## 3. Polarized Intensity Images

Figures 1 and 10 (Appendix A) show the polarized scattered-light images $Q_{\phi}$ and $U_{\phi}$, respectively, obtained in the J band. The $U_{\phi}$ image contains very low signal, suggesting that the assumption of single scattering is valid (cf. Canovas et al. 2015). Figure 1 is similar to previously published scattered-light images of HD 169142, in particular those of Momose et al. (2015) and Monnier et al. (2017), but it brings the highest signal-to-noise ratio view of the inner ring. It shows a number of features. We detect, from outside in:

(a) A faint gap (Gap #1) at ~0.70–0.73 (81–85 au). Beyond this radius, the image shows diffuse scattered light up to ~1.5 (~176 au). The marginal detection of this gap can be seen in the normalized, azimuthally averaged radial profile
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Normalized radial (left) and azimuthal intensity profiles (right) obtained after deprojection of the $r^2$-scaled J-band $Q_\phi$ image. The radial cut is obtained after azimuthally averaging and normalized to the maximum brightness of the inner ring. The red vertical dotted line is the limit of our IWA. The normalized azimuthal cuts are obtained after averaging radially between $0''_{\mathrm{4}} 14$ and $0''_{\mathrm{4}} 22$ (inner ring, purple squares) and between $0''_{\mathrm{65}}$ and $0''_{\mathrm{4}} 65$ (outer ring, green diamonds). The green curve is shifted vertically for clarity. The plotted error bars are the standard deviation in each bin in the $U_\phi$ image on a pixel basis.</figcaption></figure>

of small dust grains, the ALMA image shows no continuum detection in both the inner (Gap #3) and wide (Gap #2) gaps. This indicates that dust particles, independently of their sizes, are filtered out in the inner gap but that the filtering mechanism at play in the outer gap affects small and large dust grains differently. The inner ring (Ring #2) is also well detected in VLA 7 and 9 mm observations (Osorio et al. 2014; Macías et al. 2017), although at a slightly larger radius (~25 au, compared to ~21 au in scattered light). Furthermore, Macías et al. (2017) also reported on the detection of a third gap at ~85 au, consistent with the marginal detection in the SPHERE polarized intensity data.

## 4. Disk Modeling

We start our models by introducing planet-induced depressions in a uniform disk gas density profile to mimic the position and shape of the observed gaps. We present physical simulations, including dust evolution and trapping processes, to constrain the disk’s dust distribution and investigate whether planet–disk interactions are responsible for the detected substructures. We take the approach of fixing as many parameter values as possible and do not attempt a best-fitting procedure. Because of the high parameter degeneracy when physical processes related to dust evolution are involved, we do not explore a large grid of these models. Our concept is complementary to the one presented by Monnier et al. (2017), who showed a parameterized model without connecting the gap and ring structures to a planetary origin or dust evolution. In their model, the scale height of an inner and outer disk region and the density scaling factor for the outer gap are determined via a fitting process.

### 4.1. Model Set-ups

In our models, we consider two spatially separated planets that are massive enough to open a gap in the gas surface density. The planet cores are assumed to be at fixed orbits and are not allowed to migrate. We note that we constrain the total number of planets to two, although multiple low-mass planets in close-by resonances could exist to cause the second, wide gap

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** SPHERE/IRDIS $r^2$-scaled J-band $Q_\phi$ image overlaid with contours of the ALMA 1.3 mm continuum image (from Fedele et al. 2017). The white ellipse in the bottom left corner shows the ALMA beam with a size of $0''_{\mathrm{28}} \times 0''_{\mathrm{18}}$. The $Q_\phi$ image is normalized in the same way as in the right panel of Figure 1, but the color scale is chosen such that the structures in the outer disk are enhanced and the inner ring is saturated.</figcaption></figure>

(Gap #2). The perturbed gas surface density profiles $\Sigma_g$ depend on the planet masses and disk viscosity. To derive $\Sigma_g$, we consider the analytical solution of Crida et al. (2006), in which the gravitational and pressure torques are assumed to be zero very close to the planet. For this reason, we implement a correction for the depth of the gap using the empirical relation from Fung et al. (2014). The resulting gas surface density distributions are used as inputs to model the dust evolution considering the dust dynamics, including the processes of coagulation, fragmentation, and erosion of dust particles (Birnstiel et al. 2010; Pinilla et al. 2015a). For the background surface density profile, we use an exponentially tapered power law with a power index of 1 and a tapered radius of six times the location of the inner planet. As the inner gap appears
relatively devoid of scattering material and free from larger grains, we consider a planet-to-stellar mass ratio of $2 \times 10^{-3}$ (3.5 $M_{Jup}$) for the inner planet, such that the gap is deep enough to lead to a filtration of particles of all sizes. For the second gap, filled with small particles, we consider the mass of the outer planet close to the mass estimate obtained in Osorio et al. (2014) in this region, and we choose 0.7 and 0.3 $M_{Jup}$, the latter being the minimum mass needed to open a gap in the gas surface density (and hence to have a pressure trap at the outer edge of the gap) under our assumptions. The locations of the planets are chosen according to the current SPHERE observations and are $r_1 = 14$ au and $r_2 = 53$ au, such that the pressure maxima are close to the observed peaks of the mm emission. The companion masses considered in our simulations are compatible with the detection limits obtained in total intensity with IRDIS and the Integral Field Spectrograph (IFS) of SPHERE. These data will be presented in detail in the F100 SPHERE High-Contrast Imaging Survey for Exoplanets (SHINE) data analysis and detection performance paper (M. Langlois et al., in preparation). The disk temperature profile is a power law ($\sim r^{-0.5}$; cf. Equation (25) in Birnstiel et al. 2010) such that at 1 au the temperature is $\sim 230$ K. We assume an $\alpha$ viscosity of $10^{-3}$ throughout the disk and note that this choice also influences the planet masses assumed as described above. Furthermore, we consider a disk mass of $5 \times 10^{-3} M_{\odot}$, which is consistent with the value range reported by Panić et al. (2008), and a disk radial extension from 1 to 300 au. The initial gas-to-dust ratio is 100, and particles are initially 1 $\mu$m in size. The model follows the evolution of 180 grain sizes (from 1 $\mu$m to 2 m) and calculates the dust density distribution at each radius for timescales from $10^4$ to $10^7$ yr. We do not consider the effect of ice lines on the dust dynamics.

To compute synthetic images, we consider the resulting dust distribution as input to the radiative transfer code RADMC-3D (Dullemond et al. 2012). From the vertically integrated dust density distribution, we derive the dust density for each grain size $\Sigma_d(r, a)$. From the temperature profile $T(r)$ used in the dust evolution, the pressure scale height $H_p(r)$ is determined. We take the approach in Pohl et al. (2016) and calculate the dust scale height for each grain size $a$ following Birnstiel et al. (2010) as

$$H_d(r, a) = H_p(r) \times \min \left( 1, \sqrt{\frac{\alpha}{\min(\text{St}, 1/2)(1 + \text{St}^2)}} \right), \quad (3)$$

where $\alpha$ is the turbulent viscosity and St is the Stokes number, a dimensionless parameter that indicates how well a dust grain is coupled to the gas. In the Epstein regime, valid for most regions of protoplanetary disks and where the molecular hydrogen mean free path is larger than 4/9 times the grain size, the Stokes number at the midplane can be written as

$$\text{St} = \frac{\rho_s a \pi}{\Sigma_g 2}, \quad (4)$$

where $\rho_s$ is the volume density of the dust grain, typically $\sim 1.2 \text{ g cm}^{-3}$ according to the averaged values of the volume density for silicates. Dust grains with sizes corresponding to $\text{St} \sim 1$ are subject to the strongest gas drag and move fast to the regions of pressure maxima (Brauer et al. 2008). From the dust surface density and scale height, we compute the volume

density profile for each grain size as

$$\rho(R, \phi, z, a) = \frac{\Sigma_d(R, a)}{\sqrt{2\pi} H_d(R, a)} \exp\left(-\frac{z^2}{2H_d(R, a)^2}\right), \quad (5)$$

where $R = r \sin(\theta)$ and $z = r \cos(\theta)$ are cylindrical coordinates and $\theta$ is the polar angle. The opacity calculation of each grain size bin takes into account porous spheres with a dust mixture composed of astronomical silicates (Draine 2003), carbonaceous material (Zubko et al. 1996), and water ice (Warren & Brandt 2008). The fractional abundances of 7%, 21%, and 42% (amount of vacuum is 30%) are adopted from Ricci et al. (2010). The temperature structure of each dust grain size is determined with a Monte Carlo radiative transfer simulation, and synthetic scattered-light images are computed including the full treatment of polarization. Mie theory is used to compute the Mueller matrix elements. These images are convolved by an elliptical Gaussian PSF ($0.034 \times 0.041$) chosen to mimic the angular resolution of our SPHERE observations, and each pixel is multiplied by the square of its distance to the star to compensate for the stellar illumination drop-off with distance. For the synthetic mm observations, we consider a beam size of $0.3 \times 0.2$ (Fedele et al. 2017).

For comparison reasons, we additionally perform simplified models by neglecting the self-consistent dust evolution—that is, the dust growth—and its dynamics. However, there are mm grains in these models, so significant evolution has taken place here as well. In this second approach, we consider the same initial gas density profile perturbed by the two giant planets and assume a fixed gas-to-dust ratio exploring the range from 50 to 100, a simple approach typically used to compare scattered-light images with hydrodynamical simulations of planet-disk interaction (e.g., Dong & Fung 2017). This approach is valid as long as the micron-sized particles are well coupled to the gas and no self-consistent dust settling is included. However, these simplified models are expected to differ from dust evolution models because several processes, such as growth and fragmentation, can change the dust distribution in the disk, which at the same time changes the dynamics, in particular when pressure maxima are present (cf. Figure 5). In these simplified models, an average opacity is used considering a power-law distribution for the grain size, where the number density follows $n(a) \propto a^{-3.5}$ with $a_{\min} = 0.01 \mu\text{m}$ and $a_{\max} = 1 \text{ mm}$. For all models, we consider the stellar parameters $T_{\text{eff},*} = 8400$ K, $M_* = 1.65 M_{\odot}$, and $R_* = 1.5 R_{\odot}$, hence $10 L_{\odot}$ (Dunkin et al. 1997; Blondel & Djie 2006; Fedele et al. 2017). The stellar luminosity adopted in Fedele et al. (2017) is based on the new distance estimate from Gaia ($d = 117$ pc). For the stellar spectrum, a Kurucz spectrum of a star with metallicity $[\text{Fe/H}] = 0$ and a surface gravity of $\log g = 4.5$ is taken into account (cf. Folsom et al. 2012).

In addition, as we do not know the shape of the innermost disk (masked by the coronagraph), and because the presence of a tiny amount of dust could alter the brightness signal close to the inner peak, we set the dust density to a floor value within 10 au and apply a smooth Gaussian taper to create a rounded inner rim for the inner ring. This step is especially needed in the simplified models, because in the models with dust evolution included, most of the dust particles are filtered out and trapped at the outer edge of the gap opened by the innermost planet;
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Vertically integrated dust density distribution after 5 Myr of evolution, when two massive planets (left: 3.5 and 0.7 $M_{Jup}$; right: 3.5 and 0.3 $M_{Jup}$) are embedded in the disk at 14 and 53 au, respectively.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Vertical disk dust density structure assumed for the radiative transfer calculations following Equation (5) for an outer planet mass of 0.7 $M_{Jup}$. The cumulative density distribution for small dust grains only from 1 to 5 µm (left), all dust grains (middle), and the simplified parametric approach (right) are shown. Note that the radial scale is logarithmic for a better visualization.</figcaption></figure>

hence, in these models, the inner disk is almost empty of
grains. However, we note that a tiny inner disk exists as
presented in Lazareff et al. (2017).

## 4.2. Results

Figure 5 (left panel) shows the dust density distribution after 5 Myr of evolution for planet masses of 3.5 and 0.7 $M_{Jup}$ located at 14 and 53 au, respectively. A pressure bump is formed at the outer edge of each gap that acts as a particle trap and helps to reduce the radial drift. The higher the mass of the planet, the more efficient the trapping and higher the mm flux (Pinilla et al. 2012, 2015b). This trend is also seen in Figure 5 (right panel), where the second planet has a lower mass (0.3 $M_{Jup}$), leading to less efficient trapping there. Although the trapping of mm grains is effective, the values for the planet mass (at the considered disk turbulence) are chosen such that the small grains are not fully filtered out. Note that there is a degeneracy between the choice of disk mass, temperature, $\alpha$ turbulence, and planet mass; thus, we do not claim to infer mass limits for potential planets. The density contrast between the two rings also depends on whether the planets formed simultaneously or sequentially (Pinilla et al. 2015a) and whether they migrate. Applying Equation (5) to the dust density distributions results in the vertical density structure illustrated in Figures 6 (left and middle panels) and 12 (Appendix C). While the small grains are distributed radially over the disk extension and all the way up to the disk surface layers according to their dust scale

height, the large grains are concentrated at the pressure bump
regions close to the midplane. For comparison, Figure 6 (right
panel) shows the vertical density structure of our simplified
approach, where no grain growth model is involved and a fixed
gas-to-dust ratio of 100 is considered. In this case, a larger amount
of dust is still present within the two gap regions and in the outer
disk because the dust radial drift is neglected.

4.2.1. *Scattered Light*

Figure 7 shows the radial surface brightness profiles of the $Q_{\phi} \times r^{-2}$ model images compared to the observational radial profile (also shown in Figure 3). All profiles are normalized to the peak flux of the inner ring (Ring #2). The brightness of this ring in scattered light is highest because of the geometry of the disk scattering surface. The incidence angle of stellar radiation is steepest here so that the disk receives and scatters the most light per unit surface area. This is due to density effects given the large amount of dust there and the temperature profile. The curves in the left panel are based on our dust evolution modeling approach. The solid green curve representing the model with a higher outer planet mass well reproduces the two main ring locations observed with SPHERE, but there is a discrepancy for the brightness contrast between the rings. Although the overall width of Gap #2 matches the observations, a gap much deeper and with sharper edges than observed is produced by our model. Reducing the mass of the outer
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of surface brightness radial profiles of the $Q_\phi \times r^2$ model image at the J band. The values are scaled by the square of the distance from the central star in order to compensate for the falloff of the stellar irradiation. Models with self-consistent dust evolution are shown in the left panel, and those with a simplified parametric grain size distribution are shown in the right panel. The vertical dashed lines mark the positions of the brightness peaks (left) and the location of Gap #2 (right).</figcaption></figure>

planet helps to create a slightly shallower outer gap edge (dash-dotted blue curve). However, in this case, the amount of dust that is trapped in the outer disk is lower, and, consequently, the brightness of Ring #1 decreases. In addition, the peak of this outer ring moves slightly inward when the planet mass is reduced, differing from the observations. In both cases, the outer disk in our models appears too faint in scattered light, as all dust grains originally located in the outer disk have already grown and drifted inward, even after an evolutionary time of 0.5 Myr (see Appendix D, Figure 13 and dashed olive curve in Figure 14).

The orange line in the right panel illustrates a good match for our simplified approach that ignores dust growth and fragmentation processes and assumes a power law for the dust size distribution. The gap location between the two main rings can be reproduced well when moving the second planet position further in from 53 to 42 au. The reason why the planet position needs to be different between the dust evolution models and the simplified approach is as follows. The gap in the small grains is similar in shape to that of the gas surface density, since the small grains are well coupled to the gas. Contrarily, in the dust evolution models, there is a dominating peak of emission at the pressure maximum, which is further out from the outer edge of the gas gap and where small grains are continuously reproduced by fragmentation due to turbulent motions (Pinilla et al. 2012). The shallowness of the gap is a better match compared to the dust evolution approach. However, a simultaneous fit of gap depth and ring positions is not possible either. This trend is also seen in the HD 169142 model-fitting results by Monnier et al. (2017), where either the gap depth or the outer ring position is off compared to the GPI J- and H-band profiles.

**Mixed midplane-surface dust models.** At this point of the analysis, it seems that the model without dust evolution does a significantly better job reproducing the SPHERE scattered-light observations, which is, however, not the case for the mm dust continuum, as demonstrated later in Section 4.2.2. One has to keep in mind that dust evolution assumptions are developed for the disk midplane, where dust growth is quite efficient due to the high densities. Since the dust evolution models are only 1D, the vertical disk structure chosen influences the situation at larger height. More precisely, the coagulation equation itself is

not only calculated in the midplane, but it averages the processes with presumed weights over the vertical structure. Then, it assumes that the size distribution at a given location develops as a whole, followed by a redistribution of the grains. There is also a reservoir of small grains produced that are going through the growth and fragmentation cycle. It might be that the vertical exchange in the dust evolution is not working properly and that there are small grains at the disk surface that do not grow quickly at high altitude where the densities are lower. If this population at the top layer is indeed isolated, its coagulation compared to the midplane situation will also be on a different timescale. Moreover, charging effects could play a role for dust evolution processes at the disk surface, which would keep the particles very small.

Hence, as a test, we introduce a new population of small grains (0.01–0.5 µm) that follow the initial gas density distribution of the dust evolution model with a mass fraction of $8 \times 10^{-4} \Sigma_{\text{disk,gas}}$. The mass in the other size bins is reduced correspondingly to keep the same dust mass as that for the original dust evolution simulation. This model is displayed in Figure 7 with the dash-dotted dark green line. It helps to decrease the gap depth and increase the scattered light in the outer disk.

We note that our models that do not contain any dust inside 10 au provide a good match to the SED for wavelengths longer than 10 µm that trace the outer disk. However, the addition of a small inner belt between ~0.05 and 0.09 au allows us to reproduce the NIR excess as well. To not be seen in our scattered-light model predictions, any dusty material in the inner disk must be confined within ~0.09 au.

### 4.2.2. Millimeter Dust Continuum Emission

Figure 8 shows synthetic mm continuum images at 1.3 mm for our two representative models (solid lines in Figure 7) alongside the ALMA data from Fedele et al. (2017). The left panel considers our self-consistent dust growth model. This results in an inner dust cavity, an inner ring between ~15 and 35 au (0.′′13 and 0.′′3), and an outer ring between ~55 and 80 au (0.′′47 and 0.′′68), with a gap in between. Both the inner cavity and the gap are depleted in mm-sized dust particles. Our model is qualitatively consistent with the ALMA dust continuum image showing rings at ~20–35 au and
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Simulated intensity images of HD 169142 for ALMA Band 6 (1.3 mm) based on the radiative transfer models. The images are smoothed to the same angular resolution as in Fedele et al. (2017). The intensity units are mJy beam⁻¹.</figcaption></figure>

~56–83 au (0"48–0"64). This agreement supports the view that the efficient dust-trapping scenario by means of the two giant planets may be at work in HD 169142. The relatively sharp outer edge of the continuum map gives further evidence of large dust grains radially drifting inward (cf. Birnstiel & Andrews 2014; Facchini et al. 2017). For completeness, the synthetic image for our simplified fixed gas-to-dust ratio model without dust evolution treatment is also shown. The clear depletion of dust particles within the gap region is not seen in this case. Furthermore, the outer ring is more extended, and both the inner and outer edges are less well defined, leading to a fuzzier overall disk structure. We note that the flux is underpredicted in both model scenarios compared to the actual ALMA measurement.

# 5. Discussion

## 5.1. Fragmentation and Trapping Efficiency

The main reason for the discrepancy between the dust evolution models and the simplified approach is the fact that, in the dust evolution models, micron-sized particles efficiently grow to larger sizes already at early times of evolution. The growth changes their coupling to the gas (i.e., their Stokes number) and hence their dynamics (e.g., their drift velocities increase when they grow). In these models, small grains are continuously reproduced thanks to destructive collisions that can occur because of turbulent motions and radial drift. In the particular case of two planets embedded in the disk presented in this study, the radial drift is reduced at the pressure bumps, and fragmentation occurs due to turbulence that replenishes these regions with small grains. These small grains are more affected by turbulent motions and more difficult to trap, and thus they can be dragged along with the gas. For this reason, a small amount of micron-sized particles can still flow through the gap (cf. Figure 5). This amount of small grains is, however, not enough to reproduce the observed surface brightness profile inside the gap (as it is in the case of a constant gas-to-dust ratio). For sub-micron grains, this amount would be significantly higher, and these smaller grains scatter more efficiently in our direction, too. A possible solution for this discrepancy is to make fragmentation more efficient, for example, by increasing the turbulent motions of the grains, that is, increasing

the α viscosity. However, when α increases, a more massive planet is needed to open a gap (Crida et al. 2006), which can lead to a new discrepancy with the gap width. Moreover, the higher the turbulent motions, the more difficult it is to trap the mm-sized particles, because of the high dust diffusion that allows particles to escape from pressure bumps (de Juan Ovelar et al. 2016). As an alternative, in order to have less growth and more fragmentation, the maximum fragmentation velocity threshold of the particles can be decreased, which mainly depends on the grain composition and its structure. Nonetheless, while having more fragmentation might help to increase the surface brightness inside the gap, this can also lead to less dust trapping, which can create differences with the current mm observations. Note that decreasing the initial minimum grain size in our current dust evolution models would not help us have a better match with observations, because small grains quickly grow regardless of their initial size.

## 5.2. Dust Evolution as a Function of z

Following the analysis in Section 4.2.1, it becomes clear that the assumptions in current 1D dust evolution models are tuned for coagulation processes happening in the disk midplane. The evolution of gas and dust is modeled in a vertically integrated way assuming a steady-state disk model, although there might actually be a strong dependency on the vertical disk height. Hence, our treatment of vertical exchange and vertical settling might also be inaccurate. For instance, if there is indeed a population of small grains isolated at the top layer, this would suggest very weak turbulence. This provokes quite efficient settling, even for small grains. What might work is a population of small, charged grains that is kept from settling by, for example, magnetic fields. Hence, these small grains would be unaffected by efficient dust growth and could be permanently present at the disk surface. As a consequence, scattered-light detections at optical and NIR wavelengths would not be affected by significant dust growth. Contrarily, the surface layers might have higher turbulence, which is expected because they are hotter and highly ionized (e.g., Dzyurkevich et al. 2013). However, with only higher turbulence, grains are also mixed downward and get into contact with the lower turbulence regions deeper in the disk, where they settle and take part in the coagulation there. Thus, a locally higher turbulence is not a way to isolate grains; it is a way to move
them faster. Consequently, it could be that the gas velocities at very high altitude are such that fragmentation also works in a thin surface layer to locally replenish the reservoir of even sub-micron grains.

### 5.3. Mass of Gap-opening Planets

As shown in Section 4.2, planets with masses of 3.5 and 0.7 $M_{\text{Jup}}$ located at 14 and 53 au are needed in our dust evolution models in order to create effective pressure bumps that trap particles at the location of the rings seen in scattered light. Note that these values are compatible with the mass detection limits derived from contrast curves in total intensity SPHERE IRDIS and IFS data. The minimum planet mass limit in our model is chosen such that the planet perturbs the gas profile and efficient trapping can be generated ($M_p \gtrsim 0.3 M_{\text{Jup}}$). While a planet mass of 0.3 $M_{\text{Jup}}$ is too low to clear the full extent of Gap #2, the 0.7 $M_{\text{Jup}}$ planet is able to reproduce the gap width. The presence of multiple planets below this mass whose gaps overlap is an alternative possibility (e.g., Dodson-Robinson & Salyk 2011). Numerical studies have shown that less massive planets do not open a gap in the gas but effectively open a gap in the dust (Paardekooper & Mellema 2004, 2006; Picogna & Kley 2015; Dipierro et al. 2016; Rosotti et al. 2016; Dipierro & Laibe 2017). The gas azimuthal velocities can be perturbed such that the drift velocities of the particles are reduced, leading to a traffic-jam effect without creating local pressure maxima. In addition to the gas viscous forces, Dipierro et al. (2016) and Dipierro & Laibe (2017) also included the contribution from the tides of an embedded planet and showed that a low-mass planet can open a gap in the dust only if the tidal torque exceeds the drag torque outside the planetary orbit. In this scenario, a shallow gap can be carved out, but it is rather unlikely that this effect can create the strong rings in the distribution of small and large grains in HD 169142. It should be tested whether a combination of pressure bumps, self-consistent dust evolution, and the consideration of disk-planet tidal interactions can lead to a coherent picture for the gap and ring appearances.

Pioneering studies from Kanagawa et al. (2015), Rosotti et al. (2016), and Dong & Fung (2017) look at the inverse problem, meaning to derive planet masses from observed gap profiles. For this method, a number of assumptions about the disk structure and dynamics are made when simulating the gap shape and deriving the correlation with planet mass. The mass of the putative second planet in our HD 169142 model is consistent with the numerical analysis by Dong & Fung (2017), who derived disk and planet properties based on the morphology of gaps in NIR scattered-light images. They estimated a mass between 0.2 and $2.1 M_{\text{Jup}}$ for an $\alpha$ viscosity varying from $10^{-4}$ to $10^{-2}$. Kanagawa et al. (2015) suggested a mass $\gtrsim 0.4 M_{\text{Jup}}$ by measuring the gap depth in VLA 7 mm data. Although this is principally consistent with the other estimates, a measurement based on mm data only is complicated due to dust/gas coupling effects (Rosotti et al. 2016). This makes an exact definition of the gap width difficult, and its value depends on the disk lifetime. As discussed in Rosotti et al. (2016), a more robust indicator of the planet mass from (sub-)mm images is the location of the bright ring tracing the gas pressure maximum. This is the reason why we intend to reproduce the ring positions rather than the gap locations with the modeling approach in this paper. The inclusion of dust growth and fragmentation processes would certainly change the conclusions from Rosotti et al. (2016) and Dong & Fung (2017) as dust evolution dynamics

affects the gap depth, the slope of the gap edges, and the position of the rings and their contrast.

### 5.4. Dust Evolution Timescale

All of our results based on dust evolution modeling consider a dust evolutionary timescale of 5 Myr and that the giant planets embedded were formed simultaneously. The disk and planet age can affect the appearance of the radial profiles in polarized intensity at the NIR (see Figure 14 in the Appendix D) and in total intensity at mm wavelengths. On the one hand, the outer ring (Ring #1) becomes narrower at longer times of evolution, which produces a rather sharp outer disk edge and shifts it toward smaller radii. This in turns lowers the brightness signal in the outer disk. While this is consistent with the mm data, the amount of small dust particles decreases with time, and the NIR observations cannot be reproduced. If longer times of evolution are taken (~10 Myr, which is consistent with the revised age of the system), there would be a higher discrepancy between the dust evolution models and the NIR observations unless additional trapping mechanisms play a role all across the disk. In contrast, at very early timescales of 0.1–0.5 Myr, when all grains are not yet at the pressure maxima, the wide gap (Gap #2) remains shallower. On the other hand, analogous to the sequential planet formation scenario presented in Pinilla et al. (2015b), it could be possible that the outer planet forms earlier than the inner planet (or vice versa). This can affect the amount of dust in both traps and adjust the contrast between the two rings. However, we do not have any constraint on whether the two planets have been forming at the same time or consecutively. Together with the uncertainty when the putative planets have been forming at all, this means that the dust evolution after 5 Myr could still be a good proxy for the situation in the HD 169142 system.

### 5.5. Gaps and Rings in the Context of Ice Lines

Although observational signposts of embedded planets are the most widely used explanation to interpret ring structures in disks, the relation to ice lines of various materials is another possible scenario (Zhang et al. 2015; Okuzumi et al. 2016). Ice lines of different volatile species can significantly affect the dynamics of dust evolution processes including growth and fragmentation, which in turn has an effect on the observational appearance of rings and gaps at different wavelengths (Pinilla et al. 2017). The freeze-out temperatures of the main volatiles, such as water ($H_2O$), ammonia ($NH_3$), carbon dioxide ($CO_2$), and carbon monoxide ($CO$), are estimated to have average values of ~142, ~80, ~66, and ~26 K, respectively (Zhang et al. 2015).

Our radiative transfer models show that the midplane temperatures at the inner ring position are such that the $H_2O$ and $NH_3/CO_2$ ice lines are located close to the inner and outer edges of this ring, respectively, at mm emission. It is recognizable that the gap at scattered light lies between the ice lines of $H_2O$ and $CO_2$ when comparing the surface layer temperatures with the volatile freeze-out temperatures. Thus, the $H_2O$, $NH_3$, and $CO_2$ ice lines nearly coincide with the scattered-light ring positions. Furthermore, the CO ice line at the midplane is located at ~110 au, which is close to the outermost gap at ~85 au, consistent with $\text{DCO}^+(3-2)$ and $\text{C}^{18}\text{O}(2-1)$ ALMA observations presented in Macías et al. (2017). In Figure 9, there is an uncertainty for the specific location of these ice lines that depends on the freezing temperatures that we assume, and, for the CO ice line, it could be between ~95 and 145 au. The current observations suggest that
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Temperature map as a function of radius and polar angle in spherical coordinates for the simplified model. The ice lines for H₂O, NH₃, CO₂, and CO are indicated with yellow contours. The vertical dashed lines illustrate the edges of the two bright rings in mm dust emission.</figcaption></figure>

the accumulation of large dust grains close to the CO ice line is a
possible mechanism to explain the origin of this outermost gap.

5.6. Shadowing Effects and Time Variability

We note that because of the modeling procedure (analytical
gas profile coupled with 1D dust evolution), the observational
signatures presented in this paper are always azimuthally
symmetric. This is, for example, not necessarily true for a
massive enough planet for which an eccentric gap and vortex
formation at its edge are expected (e.g., Ataiee et al. 2013). As
mentioned in Section 2, there are significant asymmetries along
both the inner and outer ring regions in polarized intensity. It is
noticeable that the maximum polarization of the outer ring is
along the minor axis. This is opposite to several other disks
showing a brighter polarized intensity along the major axis,
expected due to the polarization efficiency being highest for
90° scattering in the Rayleigh and Mie scattering regime. A
significant scattering angle effect is also not expected at the low
inclination of HD 169142. Momose et al. (2015) invoked
corrugations of the scattering surface in the outer region as a
possible origin. Alternatively, asymmetries in the outer disk
emission can be caused by shadowing of the inner disk region.
For HD 169142, there is a slightly inclined inner disk at sub-au
distance (i = 21°, PA = 100°–130°; Lazareff et al. 2017) that
is known to be variable and might contain an extended dust
envelope, as suggested by Wagner et al. (2015). Azimuthal
brightness variations in the inner ring of the scattered light
could be caused, for example, by perturbations by a
protoplanet, optical depth variations through the suggested
dust envelope, accretion flows, or turbulence in the inner disk.
The local brightness enhancements along the innermost
scattered-light ring at different time epochs are discussed in
Ligi et al. (2017). The azimuthal inhomogeneities in the inner
ring (Ring #2) could in turn cause radial shadowing on
Ring #1 and the remaining outer disk. The pace of variations
in the illumination pattern of the outer disk depends on the
precession timescale of the inner disk material. Given that there
is no apparent difference in the rings' brightness asymmetries
in the three observational data sets in polarized intensity
(Subaru/HiCIAO: Momose et al. 2015; Gemini South/GPI:
Monnier et al. 2017; VLT/SPHERE: this paper) that span a
time period of 3 yr, the shadowing scenario for the outer disk
seems rather unlikely. However, it cannot be ruled out either, as
the precession timescale for the inner disk, for example, in the

context of a hypothetical star–companion system, can be rather
long (several hundred to thousand years).

**6. Conclusions**

In this paper, we present scattered-light observations of the
protoplanetary disk around the Herbig Ae star HD 169142
obtained with the VLT/SPHERE at the J band and compare
our results with recent ALMA data of this target. Together with
TW Hya, HD 163296, and HD 97048, it is one of a handful of
disks around young stars that have been observed at very high
angular resolution at NIR and mm wavelengths, where in each
case both images show similar substructures even if their scales
differ but also show different morphologies. For HD 169142,
we confirm the previous detection of two ringlike features
separated by a wide gap and an additional inner gap and report
on the marginal detection of a third gap in the outer disk, as
well as azimuthal brightness variations along both rings. We
present azimuthally symmetric radiative transfer models based
on planet–disk interaction processes that account for the main
observational features and discuss the influence of dust
evolution and particle trapping on the gap and ring properties.
We place our findings in the context of planet masses inferred
from the gap-opening process and ice-line chemistry. Our
measurements and modeling results suggest the following:

1. The location and width of the gap, as well as the peak positions in polarized scattered light, of HD 169142 are reproduced with our model based on dust evolution processes when two giant planets of 3.5 and 0.7 *M*<sub>Jup</sub> are embedded in the disk. The observed gap, however, possesses a shallower outer flank than expected for planet–disk interaction signatures. There is also a significant discrepancy in the gap depth, as micron-sized particles rapidly grow in the presence of pressure bumps. Small grains distributed all over the disk wherever there is gas, as in our simplified approach, decreases the gap depth such that there is good agreement with the observed shape. This also helps to increase the scattered-light flux in the outer disk. A more efficient fragmentation by increasing the turbulent motion of dust particles or to adjust the fragmentation velocities could help to overcome this deficit in small grains. Including the contribution from the tides of an embedded planet can lead to a shallower dust gap, in case the planet hypothesis is correct at all.

2. The assumptions in current dust evolution models are tuned for the disk midplane, and the vertical exchange does not work properly; thus, the coagulation timescale might be different at higher disk altitudes. A population of small (sub-)micron-sized grains might exist in the upper surface layers that is unaffected by quick growth due to lower densities and different turbulence there. Thus, the vertical disk structure and its consequences for dust evolution processes also have a significant role in interpreting scattered-light images.

3. In order to obtain a consistent picture with the mm observations, the accumulation of large grains in the dust trap of a pressure bump is needed. This generates the bright emission rings and the sharp outer disk edge as detected in the mm continuum image of HD 169142. A simplified parameterized dust size distribution is not able to reproduce the high dust depletion factor required.

4. A scenario with a grain size-dependent gap opening that still allows a perturbation in the radial pressure gradient is
required. We emphasize that inferring the mass of gap-
opening planets from simplified models is degenerate and
depends on the choice of disk mass, temperature, and α
turbulence. Constraining planet masses becomes even
more uncertain when including more physical processes
that are expected to occur in protoplanetary disks, such as
grain growth, fragmentation, and vertical disk instabilities.

5. Observing the total amount of gas and using different techniques that allow us to get better constraints on the grains sizes in disks, such as mm-wave dust polarization (Kataoka et al. 2015, 2016; Pohl et al. 2016; Yang et al. 2016), may allow us to further explain the origin of the gaps and derive the properties of potential embedded planets.

Eventually, 2D dust evolution models are needed in order to
have a self-consistent treatment of radial transport and vertical
settling and consider turbulence changes across the vertical
direction of the disk. In principle, we can start to use
multiwavelength analyses, such as the one presented in this
paper, to provide feedback on the model assumptions and
calibrate our understanding of microphysical dust processes
(sticking, fragmentation, compact versus fluffy grains, etc.).

We acknowledge the ESO Paranal staff for their support during
the observations. We thank J. Monnier for sharing his data. We
are thankful to L. Pérez, A. Garufi, and G. Bertrang for the useful
discussions about this project. A.P. is supported by the Reimar-
Lüst Fellowship of the Max Planck Society and a member of the
International Max Planck Research School for Astronomy and
Cosmic Physics at Heidelberg University, IMPRS-HD, Germany.
We acknowledge financial support from the Programme National
de Planétologie (PNP) and the Programme National de Physique
Stellaire (PNPS) of CNRS-INSU. This work has also been
supported by a grant from the French Labex OSUG@2020
(Investissements d'avenir—ANR10 LABX56). M.B., F.M., and
G.P. acknowledge funding from ANR of France under contract
number ANR-16-CE31-0013 (Planet Forming Disks). P.P.
acknowledges support by NASA through Hubble Fellowship
grant HST-HF2-51380.001-A awarded by the Space Telescope
Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under
contract NAS 5-26555. H.A., M.M., and S.P.Q. acknowledge
support from the Millennium Science Initiative (Chilean Ministry
of Economy) through grant RC130007 and further financial
support by FONDECYT, grant 3150643, as well as financial
support within the framework of the National Centre for
Competence in Research PlanetS supported by the Swiss National
Science Foundation. T.B. acknowledges funding from the
European Research Council (ERC) under the European Unions
Horizon 2020 research and innovation program under grant
agreement No. 714769. D.F. acknowledges support from the
Italian Ministry of Education, Universities and Research project
SIR (RBSI14ZRHR). Q.K. acknowledges funding from STFC via
the Institute of Astronomy, Cambridge Consolidated Grant.
SPHERE is an instrument designed and built by a consortium
consisting of IPAG (Grenoble, France), MPIA (Heidelberg,
Germany), LAM (Marseille, France), LESIA (Paris, France),
Laboratoire Lagrange (Nice, France), INAF—Osservatorio di
Padova (Italy), Observatoire astronomique de l'Université de
Genève (Switzerland), ETH Zurich (Switzerland), NOVA
(Netherlands), ONERA (France), and ASTRON (Netherlands)
in collaboration with ESO. SPHERE was funded by ESO, with

additional contributions from CNRS (France), MPIA (Germany),
INAF (Italy), FINES (Switzerland), and NOVA (Netherlands).
SPHERE also received funding from the European Commission
Sixth and Seventh Framework Programmes as part of the Optical
Infrared Coordination Network for Astronomy (OPTICON) under
grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant
number 226604 for FP7 (2009–2012), and grant number 312430
for FP7 (2013–2016). This work has made use of data from the
European Space Agency (ESA) mission Gaia (http://www.
cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC; http://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.

Appendix A
$U_\phi$ Image and Polar Mapping

<figure><img src="image_11.png" /><figcaption>**Figure 10.** The *J*-band $U_{\phi} \times r^2$ image in linear scale. Each pixel is multiplied by the square root of its distance to the star, $r^2$, to compensate for the stellar illumination drop-off with radius. The normalization is similar to the $Q_{\phi}$ image, but the dynamical range of the color bar is adjusted. The region masked by the coronagraph is indicated by the gray circle. North is up, and east is toward the left.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Polar map of the $Q_{\phi} \times r^2$ image in linear scale. The vertical dotted lines indicate PA = 100° and 200°.</figcaption></figure>
## Appendix B
### Ellipse Fitting

<table><caption>**Table 1**

Ellipse Parameters Fitted to the Two Rings and Gap #2 in Our Scattered-light Images</caption><thead><tr><th></th><th>Ring #1</th><th>Gap #2</th><th>Ring #2</th></tr></thead><tbody><tr><td>ΔR.A. (mas)</td><td>28.4 ± 5.5</td><td>12.4 ± 5.3</td><td>4.4 ± 2.9</td></tr><tr><td>ΔDecl. (mas)</td><td>18.9 ± 5.7</td><td>33.5 ± 5.4</td><td>5.3 ± 2.8</td></tr><tr><td>Semimajor axis (mas)</td><td>536.4 ± 18.2</td><td>375.2 ± 14.5</td><td>173.8 ± 9.1</td></tr><tr><td>Semimajor axis (au)</td><td>62.8 ± 2.1</td><td>43.9 ± 1.7</td><td>20.3 ± 1.1</td></tr><tr><td>Semiminor axis (mas)</td><td>522.7 ± 1.2</td><td>365.6 ± 1.5</td><td>169.3 ± 2.0</td></tr><tr><td>Semiminor axis (au)</td><td>61.2 ± 0.1</td><td>42.8 ± 0.2</td><td>19.8 ± 0.2</td></tr></tbody></table>

**Note.** The position angle of the disk is fixed to 5° and the inclination to 13°. We give the offset of the ellipses from the star position, as well as the size of the semimajor and semiminor axes for each fitted feature.

## Appendix C
### Modeled 2D Disk Dust Density Distribution

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Vertical disk dust density structure assumed for the radiative transfer calculations following Equation (5) for an outer planet mass of 0.3 $M_{Jup}$. The cumulative density distributions for only small dust grains from 1 to 5 µm (left) and for all dust grains (right) are shown. Note that the radial scale is logarithmic for a better visualization.</figcaption></figure>

## Appendix D
### Effect of Dust Evolution Timescale

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Vertically integrated dust density distribution after 0.5 Myr (left) and 10 Myr (right) of evolution, when two massive planets (3.5 and 0.3 $M_{Jup}$) are embedded in the disk at 14 and 53 au, respectively.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Comparison of surface brightness radial profiles of the $Q_{\phi} \times r^2$ model image at the *J* band for different dust evolution timescales. The values are scaled by the square of the distance from the central star in order to compensate for the falloff of the stellar irradiation. The vertical dashed lines mark the positions of the brightness peaks.</figcaption></figure>

## ORCID iDs

M. Benisty https://orcid.org/0000-0002-7695-7605

P. Pinilla https://orcid.org/0000-0001-8764-1780

J. de Boer https://orcid.org/0000-0002-1028-7522

A. Zurlo https://orcid.org/0000-0002-5903-8316

T. Birnstiel https://orcid.org/0000-0002-1899-8783

C. Dominik https://orcid.org/0000-0002-3393-2459

S. Facchini https://orcid.org/0000-0003-4689-2684

M. Janson https://orcid.org/0000-0001-8345-593X

A.-L. Maire https://orcid.org/0000-0002-2591-4138

F. Ménard https://orcid.org/0000-0002-1637-7393

M. Meyer https://orcid.org/0000-0003-1227-3084

C. Pinte https://orcid.org/0000-0001-5907-5179

S. P. Quanz https://orcid.org/0000-0003-3829-7412

T. Stolker https://orcid.org/0000-0002-5823-3072

G. van der Plas https://orcid.org/0000-0001-5688-187X

A. Baruffolo https://orcid.org/0000-0002-1114-4355

## References

ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJL, 820, L40

Apai, D., Pascucci, I., Brandner, W., et al. 2004, A&A, 415, 671

Ataiee, S., Pinilla, P., Zsom, A., et al. 2013, A&A, 553, L3

Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014, ApJ, 781, 87

Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6

Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42

Béthune, W., Lesur, G., & Ferreira, J. 2016, A&A, 589, A87

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, Proc. SPIE, 7014, 18

Biller, B. A., Males, J., Rodigas, T., et al. 2014, ApJL, 792, L22

Birnstiel, T., & Andrews, S. M. 2014, ApJ, 780, 153

Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79

Blondel, P. F. C., & Djie, H. R. E. T. A. 2006, A&A, 456, 1045

Boccaletti, A., Abe, L., Baudrand, J., et al. 2008, Proc. SPIE, 7015, 70151B

Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859

Canovas, H., Hardy, A., Zurlo, A., et al. 2017, A&A, 598, A43

Canovas, H., Ménard, F., de Boer, J., et al. 2015, A&A, 582, L7

Carrasco-González, C., Henning, T., Chandler, C. J., et al. 2016, ApJL, 821, L16

Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJL, 785, L12

Crida, A., Morbidelli, A., & Masset, F. 2006, Icar, 181, 587

de Boer, J., Salter, G., Benisty, M., et al. 2016, A&A, 595, A114

de Juan Ovelar, M., Pinilla, P., Min, M., Dominik, C., & Birnstiel, T. 2016, MNRAS, 459, L85

Dipierro, G., & Laibe, G. 2017, MNRAS, 469, 1932

Dipierro, G., Laibe, G., Price, D. J., & Lodato, G. 2016, MNRAS, 459, L1

Dodson-Robinson, S. E., & Salyk, C. 2011, ApJ, 738, 131

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, Proc. SPIE, 7014, 70143

Dong, R., & Fung, J. 2017, ApJ, 835, 146

Dong, R., Fung, J., & Chiang, E. 2016, ApJ, 826, 75

Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93

Draine, B. T. 2003, ApJ, 598, 1026

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, ascl:1202.015

Dunkin, S. K., Barlow, M. J., & Ryan, S. G. 1997, MNRAS, 286, 604

Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114

Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16

Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72

Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68

Folsom, C. P., Bagnulo, S., Wade, G. A., et al. 2012, MNRAS, 422, 2072

Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88

Fusco, T., Rousset, G., Sauvage, J.-F., et al. 2006, OExpr, 14, 7515

Fusco, T., Sauvage, J.-F., Petit, C., et al. 2014, Proc. SPIE, 9148, 1

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2

Garufi, A., Meeus, G., Benisty, M., et al. 2017, A&A, 603, A21

Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105

Ginski, C., Stolker, T., Pinilla, P., et al. 2016, A&A, 595, A112

Grady, C. A., Muto, T., Hashimoto, J., et al. 2013, ApJ, 762, 48

Grady, C. A., Schneider, G., Hamaguchi, K., et al. 2007, ApJ, 665, 1391

Kama, M., Bruderer, S., Carney, M., et al. 2016, A&A, 588, A108

Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJL, 806, L15

Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78

Kataoka, A., Tsukagoshi, T., Momose, M., et al. 2016, ApJL, 831, L12

Kuhn, J. R., Potter, D., & Parise, B. 2001, ApJL, 553, L189

Langlois, M., Dohlen, K., Vigan, A., et al. 2014, Proc. SPIE, 9147, 1

Lazareff, B., Berger, J.-P., Kluska, J., et al. 2017, A&A, 599, A85

Ligi, R., Vigan, A., Gratton, R., et al. 2018, MNRAS, 473, 1774

Macías, E., Anglada, G., Osorio, M., et al. 2017, ApJ, 838, 97

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, Proc. SPIE, 9908, 990834

Marino, S., Perez, S., & Casassus, S. 2015, ApJL, 798, L44

Meeus, G., Pinte, C., Woitke, P., et al. 2010, A&A, 518, L124

Menu, J., van Boekel, R., Henning, T., et al. 2015, A&A, 581, A107

Momose, M., Morita, A., Fukagawa, M., et al. 2015, PASJ, 67, 83

Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20

Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012, A&A, 547, A111

Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41

Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJL, 748, L22

Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82

Osorio, M., Anglada, G., Carrasco-González, C., et al. 2014, ApJL, 791, L36

Paardekooper, S.-J., & Mellema, G. 2004, A&A, 425, L9

Paardekooper, S.-J., & Mellema, G. 2006, A&A, 453, 1129

Panić, O., Hogerheijde, M. R., Wilner, D., & Qi, C. 2008, A&A, 491, 219

Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Sci, 353, 1519

Petit, C., Sauvage, J.-F., Fusco, T., et al. 2014, Proc. SPIE, 9148, 0

Picogna, G., & Kley, W. 2015, A&A, 584, A110

Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81

Pinilla, P., Birnstiel, T., & Walsh, C. 2015a, A&A, 580, A105

Pinilla, P., de Juan Ovelar, M., Ataiee, S., et al. 2015b, A&A, 573, A9

Pinilla, P., Flock, M., Ovelar, M. d. J., & Birnstiel, T. 2016, A&A, 596, A81

Pinilla, P., Pohl, A., Stammler, S. M., & Birnstiel, T. 2017, ApJ, 845, 68

Pohl, A., Kataoka, A., Pinilla, P., et al. 2016, A&A, 593, A12

Pohl, A., Sissa, E., Langlois, M., et al. 2017, A&A, 605, A34

Quanz, S. P., Avenhaus, H., Buenzli, E., et al. 2013, ApJL, 766, L2

Raman, A., Lisanti, M., Wilner, D. J., Qi, C., & Hogerheijde, M. 2006, AJ, 131, 2290

Rapson, V. A., Kastner, J. H., Millar-Blanchaer, M. A., & Dong, R. 2015, ApJL, 815, L26

Reggiani, M., Quanz, S. P., Meyer, M. R., et al. 2014, ApJL, 792, L23

Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15

Riviere-Marichalar, P., Merín, B., Kamp, I., Eiroa, C., & Montesinos, B. 2016, A&A, 594, A59

Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790
Ruge, J. P., Flock, M., Wolf, S., et al. 2016, A&A, 590, A17

Schmid, H.-M., Downing, M., Roelfsema, R., et al. 2012, Proc. SPIE, 8446, 84468Y

Schmid, H. M., Joos, F., & Tschan, D. 2006, A&A, 452, 657

Seok, J. Y., & Li, A. 2017, ApJ, 835, 291

Simon, J. B., & Armitage, P. J. 2014, ApJ, 784, 15

Stammler, S. M., Birnstiel, T., Panić, O., Dullemond, C. P., & Dominik, C. 2017, A&A, 600, A140

Stolker, T., Dominik, C., Avenhaus, H., et al. 2016, A&A, 595, A113

Strom, K. M., Strom, S. E., Edwards, S., Cabrit, S., & Skrutskie, M. F. 1989, AJ, 97, 1451

Takahashi, S. Z., & Inutsuka, S.-i. 2014, ApJ, 794, 55

Tang, Y.-W., Guilloteau, S., Dutrey, A., et al. 2017, ApJ, 840, 32

Thalmann, C., Schmid, H. M., Boccaletti, A., et al. 2008, Proc. SPIE, 7014, 70143F

Tsukagoshi, T., Nomura, H., Muto, T., et al. 2016, ApJL, 829, L35

van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132

van der Plas, G., Wright, C. M., Ménard, F., et al. 2017, A&A, 597, A32

Wagner, K. R., Sitko, M. L., Grady, C. A., et al. 2015, ApJ, 798, 94

Warren, S. G., & Brandt, R. E. 2008, JGRD, 113, D14220

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Yang, H., Li, Z.-Y., Looney, L., & Stephens, I. 2016, MNRAS, 456, 2794

Youdin, A. N. 2011, ApJ, 731, 99

Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJL, 806, L7

Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6

Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47

Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321