# Dust Unveils the Formation of a Mini-Neptune Planet in a Protoplanetary Ring

Sebastián Pérez, Simon Casassus, Clément Baruteau, Ruobing Dong, Antonio Hales, Lucas Cieza

► To cite this version:

Sebastián Pérez, Simon Casassus, Clément Baruteau, Ruobing Dong, Antonio Hales, et al., Dust Unveils the Formation of a Mini-Neptune Planet in a Protoplanetary Ring. The Astronomical Journal, 2019, 158 (1), pp.15. 10.3847/1538-3881/ab1f88. [hal-02369629](https://hal.science/hal-02369629v1)

HAL Id: hal-02369629

https://hal.science/hal-02369629v1

Submitted on 20 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Dust Unveils the Formation of a Mini-Neptune Planet in a Protoplanetary Ring

Sebastián Pérez<sup>1,2,3 ID</sup>, Simon Casassus<sup>3 ID</sup>, Clément Baruteau<sup>4 ID</sup>, Ruobing Dong<sup>5 ID</sup>, Antonio Hales<sup>6,7 ID</sup>, and Lucas Cieza<sup>8 ID</sup>

<sup>1</sup>Universidad de Santiago de Chile, Av. Libertador Bernardo O'Higgins 3363, Estación Central, Santiago, Chile; sebastian.astrophysics@gmail.com

<sup>2</sup>Departamento de Física, Universidad de Santiago de Chile, Av. Ecuador 3493, Estación Central, Santiago, Chile

<sup>3</sup>Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

<sup>4</sup>IRAP, Université de Toulouse, CNRS, UPS, Toulouse, France

<sup>5</sup>Department of Physics & Astronomy, University of Victoria, Victoria, BC, V8P 1A1, Canada

<sup>6</sup>National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA

<sup>7</sup>Atacama Large Millimeter/Submillimeter Array, Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura 763-0355, Santiago, Chile

<sup>8</sup>Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile

Received 2019 February 13; revised 2019 April 25; accepted 2019 April 28; published 2019 June 14

## Abstract

Rings and radial gaps are ubiquitous in protoplanetary disks, yet their possible connection to planet formation is currently subject to intense debates. In principle, giant planet formation leads to wide gaps that separate the gas and dust mass reservoir in the outer disk, while lower mass planets lead to shallow gaps that are manifested mainly on the dust component. We used the Atacama Large Millimeter/submillimeter Array (ALMA) to observe the star HD 169142, host to a prominent disk with deep, wide gaps that sever the disk into inner and outer regions. The new ALMA high-resolution images allow for the outer ring to be resolved as three narrow rings. The HD 169142 disk thus hosts both the wide gap trait of transition disks and a narrow-ring system similar to those observed in the TW Hya and HL Tau systems. The mass reservoir beyond a deep gap can thus host ring systems. The observed rings are narrow in radial extent (width/radius of 1.5/57.3, 1.8/64.2, and 3.4/76.0 in au) and have asymmetric mutual separations: the first and middle ring are separated by 7 au while the middle and outermost ring are distanced by ~12 au. Using hydrodynamical modeling we found that a simple explanation, involving a single migrating low-mass planet (10 $M_\odot$), entirely accounts for such an apparently complex phenomenon. Inward migration of the planet naturally explains the ring's asymmetric mutual separation. The isolation of HD 169142's outer rings thus allows a proof of concept to interpret the detailed architecture of the outer region of protoplanetary disks with low-mass planet formation of mini-Neptune's size, i.e., as in the protosolar nebula.

**Key words:** planet–disk interactions – planets and satellites: formation – protoplanetary disks – submillimeter: planetary systems

**Supporting material:** data behind figure

## 1. Introduction

Planet formation theories describe two main pathways for the formation of protoplanet embryos: gravitational instability followed by fragmentation, which forms gas giant planets, or (sub-)stellar mass companions (e.g., Boss 1997; Kratter et al. 2010); and bottom-up growth by core accretion, which forms terrestrial and icy planets, or the cores of giant planets (e.g., Pollack et al. 1996). However, the timescale of core accretion increases dramatically with stellocentric distance due to the longer dynamical timescale and low surface densities at large radii (e.g., Goldreich et al. 2004), which poses a problem for the formation of planets before the protoplanetary disk dissipates. Mechanisms such as pebble accretion (Ormel & Klahr 2010) are thought to alleviate the tension required to form the cores of Jupiter, Saturn, and the icy giants in the solar system (all within 30 au), and the lifetime of the protosolar disk (Lambrechts & Johansen 2012).

The formation and evolution of pebbles and planetesimals, and hence planetary cores, relies on radial pressure bumps that curb the catastrophic drift of solids toward the star due to aerodynamic drag (e.g., Pinilla et al. 2012). These radial dust

traps have been associated with the ring systems observed in disks in the thermal emission from cold dust grains at radio wavelengths. Forming planets are thought to open gaps in the dust distribution, and the masses of such putative protoplanets are commonly inferred from the gap's depth and width (Rosotti et al. 2016; Dong & Fung 2017). It is frequent to associate one planet per gap (Dipierro et al. 2015; Clarke et al. 2018; Mentiplay et al. 2018), multiple planets to one (wide) gap (Dong et al. 2015), and one planet to multiple (narrow) gaps (e.g., a 0.1 $M_{Jup}$ planet can explain the multiple gaps in AS 209; Zhang et al. 2018).

In the so-called gapped disks, sometimes called transitional disks, dust is evacuated from central cavities or deep radial gaps, resulting in one or two bright rings. Accreting giant planets of masses comparable to or larger than that of Jupiter are able to open deep gaps in the gas and accumulate dust into an outer ring (Crida et al. 2006; Pinilla et al. 2012; Dipierro et al. 2016). By contrast, in disks without a deep gap, dust accumulation in mild pressure maxima produces sequences of several *shallow* rings (e.g., the TW Hya disk; Andrews et al. 2016). Clear signs of low-mass planet formation, as hinted by these fine shallow rings (Dong et al. 2017, 2018a), have been absent in the outer regions of gapped disks.

The 1.7 $M_\odot$ star HD 169142 (estimated age $6^{+6}_{-3}$ Myr; Grady et al. 2007) hosts a prominent disk with deep, wide gaps imaged in the near-infrared (Quanz et al. 2013; Momose et al. 2015; Monnier et al. 2017; Pohl et al. 2017; Bertrang et al. 2018;

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

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Ligi et al. 2018), mid-infrared (Honda et al. 2012), 1.3 mm (Fedele et al. 2017), and centimeter (Osorio et al. 2014; Macías et al. 2018) wavelengths. Circumplanetary features related to giant protoplanets have tentatively been reported in HD 169142 from radio/infrared(IR) imaging (Biller et al. 2014; Osorio et al. 2014; Reggiani et al. 2014; Gratton et al. 2019), but the nature of these signals remains to be confirmed (Ligi et al. 2018).

In this paper, we present new 1.3 mm observations of the dust structures around HD 169142 at ~2 au resolution (Section 2). In these new observations, the outer disk shows an intricate system of fine rings and narrow gaps (described in Section 2.2). Using hydrodynamical modeling (Section 3), a simple connection between the fine ring structure and a single, *migrating*, mini-Neptune protoplanet can be made, which we present and discuss in Section 3.3. Implications are discussed in Section 4.

## 2. Observations and Data Reduction

We obtained 1.3 mm observations by combining the Atacama Large Millimeter/submillimeter Array (ALMA) 12 m array extended (C40-9) and more compact (C40-6) configurations, resulting in baselines ranging from 19 m to up to 13.9 km and a total of 40–46 antennas. The combined observations are sensitive to spatial scales of up to 2′′0, with a spatial resolution of ~20 mas. The long baseline observations were acquired between 2017 September and November (Cycle 4) in four different blocks of ~40 minutes each. Precipitable water vapor ranged between 0.7 and 1.8 mm. Observations of a phase calibrator (J1826–2924) were alternated with the science to calibrate the time-dependent variation of the complex gains. The cycling time for phase calibration was set to 8 minutes and 54 s for the compact and extended configurations, respectively. The ALMA correlator was configured in the frequency division mode (FDM). Two spectral windows with 1.875 GHz bandwidth were set up for detecting the dust continuum, centered at 232.0 GHz and 218.0 GHz, respectively.

All data were calibrated by the ALMA staff using the ALMA Pipeline version 40896 in the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007), including offline water vapor radiometer (WVR) calibration, system temperature correction, as well as bandpass, phase, and amplitude calibrations. The short baseline and long baseline data sets were calibrated independently. Self-calibration of the data was performed to improve coherence (a single iteration of phase-only), after which they were combined using the CASA task CONCAT. A positional offset of (−20, −40) mas between the short and long baseline resulting images was measured and corrected prior to combining the data sets. This offset is consistent with the (−2.8, −38.0) mas yr⁻¹ proper motion listed for this star (SIMBAD) and with the expected astrometric accuracy of ALMA (typically a tenth of the synthesized beam).

### 2.1. ALMA Continuum Imaging

Image reconstruction was initially performed using the CLEAN algorithm (CASA version 5.2, task tclean). As the source is relatively bright, we super-resolve the continuum data using our maximum entropy method (MEM) package UVMEM, which is part of the family of algorithms based on maximum entropy regularization. Here we used the publicly available graphics processing unit (GPU) adaptation GPUVMEM<fnref n="9" /> (Cárcamo et al. 2018), and regularized by minimizing the

Laplacian of the model image with an objective function. The MEM reconstruction provides the highest resolution without compromising on sensitivity. We adopted the MEM image for our analysis. The final image has a peak of 67 mJy beam⁻¹ and an rms (1σ) of 14.7 µJy beam⁻¹ for a synthesized beam of 27 × 20 mas. The total flux density over the whole image is 178 ± 18 mJy (10% flux calibration uncertainty). The flux density of the central source at 1.3 mm is 194 ± 20 µJy. The bright ring at ~0′′22 has an integrated flux of 59 ± 6 mJy, and the flux density of the outer region is 126 ± 13 mJy. A side-by-side comparison of the MEM versus CLEAN image reconstructions is shown in Figure 5 (Appendix A).

### 2.2. A Triple Ring System in the Disk’s Outer Region

The ALMA image (Figure 1(a)) reveals a prominent bright ring (B1, see labeling in Figure 1(b)) followed by substructure in the outer region in the form of a system of three fine rings. In previous observations the three rings (B2, B3, and B4) were seen as a single structure (we refer to this structure as the outer region), with evidence of ripples at the last scattering surface of micron-sized grains (Pohl et al. 2017; Gratton et al. 2019). In these new observations, we see that the outer region is a packed system of three rings, with mutual separations of ~7 au between B2 and B3, and ~12 au between B3 and B4, for a distance to HD 169142 of 113.9 ± 0.8 pc (Gaia Collaboration et al. 2018).

Rings B2 and B3 are detected at a 5–20σ level (1σ = 15 µJy beam⁻¹), while emission from B4 reaches 4–5σ over the image. For a quantitative characterization of the rings, we deproject the image by the disk inclination angle $i = 12.5 ± 0.5$° and transform to polar coordinates centered on the compact millimeter emission detected at the location of the star (see Appendix B). The polar dejections zoomed on each ring are shown in Figure 2. The disk azimuthally averaged surface brightness profile is shown in Figure 1(b).

The locations and widths of the fine rings are determined from Gaussian fits to the average radial profile. Three Gaussian components were fit simultaneously plus a common low-order polynomial to account for low-level background emission. Thus, the observed locations of B2, B3, and B4 are 0.′′503 ± 0.′′005 (57.3 au), 0.′′563 ± 0.′′005 (64.2 au), and 0.′′667 ± 0.′′009 (76.0 au), respectively. Uncertainties only represent the error in fitting the Gaussian centroids. Since the rings are not circular (see below), the radii of each ring depend on azimuth, varying within ±0.′′01 around the mean locations quoted here.

The deconvolved width of each ring is determined by subtracting the beamwidth $\sigma_{\text{beam}} = \text{FWHM}/\sqrt{8\log2}$ from the best-fit Gaussian width in quadrature. Rings B2 and B3 have deconvolved widths of only 1.5 and 1.8 au, respectively, while B4 is more radially extended with a deconvolved width of 3.4 au. Rings B2 and B3 in HD 169142 are likely some of the narrowest structures in protoplanetary disks reported so far (see the Disk Substructures at High Angular Resolution Project, DSHARP, ring characterization in Dullemond et al. 2018, where the narrowest ring is ~3.4 au in deconvolved width).

Figure 2 shows evidence for azimuthal dips along ring B2 spanning azimuthal angles ~50°–100° and ~180°–270°, measured counterclockwise from the disk P.A. Interestingly, all three outer rings have a finite eccentricity (0.09 ± 0.03) and share a common focus. This is demonstrated in Figure 2, which shows ellipses fit to the rings’ radii as a function of azimuthal angle (see Appendix B for details).

<footnote n="9">⁹ The GPUVMEM code is publicly available and can be found here: https://github.com/miguelcarcamov/gpuvmem.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Triple ring system in the outer region of the gapped disk around HD 169142. The 1.3 mm ALMA image is shown with a linear color scale in flux density units of mJy per beam, with a synthesize beam of 27 × 20 mas (shown as an inclined ellipse in the bottom left corner). (b) Azimuthally averaged surface brightness profile. The shaded area around the profile corresponds to ±3σ dispersion around the mean. Rings and gaps are labeled. Unresolved continuum emission is detected at the expected location of the star (which reflects the synthesized beam or resolution). The vertical dashed lines show the locations of the fine outer rings. Shaded rectangles centered on B2, B3, and B4 show the FWHM of each rings.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: polar deprojection zoomed on rings B1, B2, and B3. The structure in B1 cannot be accounted for by neither a circular nor elliptical function as it is most likely dynamically perturbed by a giant planet. The ellipses that best fit the radial centroids of B2 and B3 have both *e* = 0.09 and are shown as red solid lines. The zero azimuth angle corresponds to a position angle equal to 5° in Figure 1, and increases east of north. The color scale is linear in mJy per beam, normalized to the peak of each ring for better visualization. As ring B4 is not always detected above a 5σ level, it is not shown here. Right: intrinsic eccentricities (top) and position angles (P.A.; bottom) of rings B2 and B3 as a function of the disk inclination angle. The shaded areas show ±2σ uncertainties. In the bottom panel, the origin of the position angles is set to the disk P.A. (of 5° on the sky). Assuming B2 and B3 are coplanar, our fit constrains the inclination angle of the disk to 12.5 ± 0.5 (1σ uncertainty). This inclination value is used in the deprojection shown in the left panel. At high and low inclinations, the ellipse P.A. approaches the disk P.A., as expected.</figcaption></figure>

## 2.3. A Structured Ring and Two Deep, Wide Gaps in the Disk's Inner Region

The brightest feature at 1.3 mm is the previously known inner ring B1, whose remarkable radial and azimuthal structure is further emphasized here at finer angular resolutions (Figure 2). Previous hydrodynamical modeling have suggested

that this structure is consistent with dynamical interactions with forming giant planets located inside the gaps D1 and D2 (Bertrang et al. 2018). The narrowness of ring B1 is explained by efficient radial trapping of dust in a pressure maximum, which puts an abrupt halt to the inward radial drift of dust particles $\ge 1$ mm (Pinilla et al. 2012). Interestingly, as seen in Figures 1(a) and 2, the structure in B1 may well be interpreted
as a multiple ring. B1 is not split in the same way as the outer ring, rather it shows a very irregular and faint inner ring.

In addition, compact continuum emission is detected at the star's location. We interpret this emission as coming from thermal emission from an inner disk of ≲1 au in radius, compatible with its near-infrared excess (Chen et al. 2018) as free-free emission is not expected to be significant in HD 169142 (Osorio et al. 2014).

The massive protoplanets required to explain the perturbed morphology of B1 are unlikely to split the outer dust disk into the observed triple ring structure (Bae et al. 2018). Instead, their dynamical interaction with the disk should lead to the clearing of deep and wide gaps such as D1 and D2, which are both ~20 au in width and heavily depleted of dust and gas (Fedele et al. 2017; Pohl et al. 2017; Bertrang et al. 2018; Ligi et al. 2018). Yet, with the fine angular resolution provided by these new ALMA observations, we find that the outer region of HD 169142 is reminiscent of the ring systems in TW Hya (Andrews et al. 2016) and HL Tau (ALMA Partnership et al. 2015), both of which lack the presence of a deep, wide gap in the dust continuum.

# 3. Hydrodynamical Simulations and Radiative Transfer Calculations

## 3.1. Theoretical Background

Several models have been proposed to explain disks with rings. Enhanced dust growth at condensation fronts can produce fine rings (Banzatti et al. 2015; Zhang et al. 2015; Okuzumi et al. 2016), however, the location of rings and gaps are uncorrelated with the expected locations of snowlines (Huang et al. 2018; Long et al. 2018; van der Marel et al. 2019), and these rings are unlikely to be eccentric. Sharp variations in the gas viscosity associated with dead zones (Flock et al. 2015; Miranda et al. 2017) and radially variable magnetic disk winds (Suriano et al. 2018) also yield ring-like structures, but they are driven by magnetohydrodynamic (MHD) instabilities that occur mostly within the central ~20 au.

Theory of disk-planet interactions predicts that a low-mass planet can produce multiple dust rings by the propagation and dissipation of spiral density waves (Goodman & Rafikov 2001; Dong et al. 2017, 2018a). The dissipation of waves progressively moves dust away from the planet's orbit, with the consequence that two narrow dust gaps form, one on each side of the planet's orbit. The expelled dust thus forms two rings exterior to the gaps. Also, dust can remain at the planets (co-)orbital radius, forming a third ring between the two dust gaps (Dong et al. 2017). Previous hydrodynamic simulations have explored this scenario and found that, in a low-viscosity disk (α viscosity of $10^{-4}$ or lower), a single low-mass planet indeed produces multiple rings in the dust distribution separated by annular gaps (Dong et al. 2017, 2018a; Bae et al. 2018), observable at ALMA resolutions. The resulting ring system depends mainly on disk viscosity, disk temperature (aspect ratio), planet mass, and the elapsed time, which will dictate the number of additional spiral arms and hence the potential number of gaps (Bae et al. 2018).

Given the sharpness and small widths of the rings in HD 169142 (the deconvolved width of ring B2 is only ~1.5 au), the α turbulent viscosity should be fairly small allowing a single low-mass planet to produce multiple rings

and gaps. Assuming standard radial profiles for the disk structure, a relation that connects the distance between the two dust gaps, relative to the orbit of the planet, and the planet's mass and disk scale-height can be found. Informed by hydrodynamic calculations (assuming $\alpha \lesssim 10^{-4}$; Dong et al. 2018a), the relation reads

$$ \frac{r_{\text{OG}} - r_{\text{IG}}}{r_{\text{p}}} \approx 2.9 \left( \frac{\gamma + 1}{12/5} \frac{M_{\text{p}}}{M_{\text{th}}} \right)^{-2/5} \left( \frac{h}{r} \right), $$

where $r_{\text{IG}}$ and $r_{\text{OG}}$ are the locations of the inner and outer gaps around the planet, respectively, and $r_{\text{p}}$ is the planet's orbit. $M_{\text{th}} = M_{\star}(h/r)^3$ is the disk thermal mass, $\gamma$ is the polytropic index (equal to 1 for isothermal gas), and $h/r$ is the disk aspect ratio. Figure 1(b) shows that the narrow gaps D3 and D4 are separated by ~9 au. In the outer region, the observed flux density of the rings translates into a brightness temperature of 8–12 K, consistent with a disk aspect ratio of 3%. Thus, for HD 169142's stellar mass ($M_{\star} \approx 1.7 M_{\odot}$), disk aspect ratio ($h/r \approx 0.03$), and approximating $r_{\text{p}}$ by B3's location at ~64 au, a planet of $\lesssim 10 M_{\oplus}$ could produce gaps separated by $r_{\text{D4}} - r_{\text{D3}} \approx 9$ au.

## 3.2. Physical Model and Numerical Setup

We perform hydrodynamic simulations of an outer disk composed of gas and a distribution of dust particles to test whether the three bright rings (B2, B3, B4) unveiled by our ALMA observations can be explained by a single forming planet. The main observables are the distance between the gaps, the mutual separation between the rings, and the observed fluxes. In this study, we do not attempt to model the previously known inner features (D1 and B1) but rather focus on the new substructure of the outer region.

### 3.2.1. Hydrodynamical Simulations

Our 2D hydrodynamical simulation of the gas and dust in the HD 169142 disk has been carried out with the code Dusty FARGO-ADSG. It is an extended version of the publicly available code FARGO-ADSG, which solves the gas hydrodynamics equations on a polar mesh (Masset 2000; Baruteau & Masset 2008a, 2008b), and which models dust as Lagrangian test particles<fnref n="10" /> (Baruteau & Zhu 2016; Fuente et al. 2017). The simulation has been designed to test whether the three bright rings (B2, B3, B4) unveiled by our ALMA observations can be explained by a single forming planet. In our simulation the star's mass is assumed to be $2 M_{\odot}$. We perform a set of simulations sampling the main parameter space (planet mass, disk scale-height, disk alpha viscosity), although no exhaustive exploration is intended as finding a perfect fit is not in the scope of this work.

The gas momentum and continuity equations are solved on a polar mesh with 900 cells logarithmically spaced between 35 and 105 au in radius, and 1200 cells evenly spaced between 0 and $2\pi$ in azimuth. Nonreflecting boundary conditions are used to avoid reflections of the planet wakes at the radial edges of the computational grid.

<footnote n="10">¹⁰ The FARGO-ADSG code is publicly available at http://fargo.in2p3.fr/-FARGO-ADSG-. The version including a Lagrangian treatment of the dust particles, Dusty FARGO-ADSG, can be made available upon request to co-author C. Baruteau.</footnote>
The disk has an initial surface density profile decreasing with radius $r$ as $r^{-1}$ with an exponential cutoff beyond 100 au. The initial surface density is $2.9 \text{ g cm}^{-2}$ at 100 au, very close to the radiative transfer (RT) model which fits the CO istopologue and continuum observations at 1.3 mm (Fedele et al. 2017). This large surface density means that Type I migration cannot be discarded, especially for the low planet masses we are exploring. A locally isothermal equation of state is used with the disk temperature decreasing with radius as $r^{-1/2}$, and equal to 8 K at 70 au (motivated by the observed brightness temperature), slightly lower than previous models. This low temperature, which is required to match the small mutual separations between the bright rings, translates into a disk aspect ratio profile in $r^{1/4}$ and equal to 0.033 at 69 au. Gas self-gravity is included since the disk’s Toomre Q-parameter varies from ~6 to ~9 throughout the computational grid. The $\alpha$ turbulent viscosity is set to $10^{-5}$. This low viscosity is motivated by the sharpness of the observed rings. Accretion onto the planet is not included.

The code solves the equations of motion for 200,000 dust particles with radii between 10 µm and 3 mm. Dust particles feel the gravity of the star, the planet, the disk gas (since gas self-gravity is included), and gas drag. Dust turbulent diffusion is also included as stochastic kicks on the particles position vector (Charnoz et al. 2011). However, dust self-gravity, dust drag, growth, and fragmentation are not taken into account. Dust particles (assumed to be compact spheres of $2 \text{ g cm}^{-3}$ internal density) are inserted at the beginning of the simulation with a number density profile decreasing as $r^{-2}$ between 59 and 90 au. This corresponds to an initial dust surface density decreasing as $r^{-4}$. Although quite steep, the surface density evolves into a profile that increases with $r$ between the rings after 150 orbits. See the perturbed gas density as well as the distribution of dust particles shown in Figure 3. Dust feedback onto the gas is discarded as the dust’s surface density along the B2, B3, and B4 rings remains comfortably smaller than the gas surface density in our model (see Figure 3) for the gas-to-dust mass ratio assumed in the RT calculation (see Section 3.2.2).

A low-mass planet of a mini-Neptune mass (motivated by Section 3.1) is inserted at 69 au at the beginning of the simulation. The planet-to-star mass ratio is set to $1.7 \times 10^{-5}$, which for a $2 M_{\odot}$ star translates into $11 M_{\oplus}$.<fnref n="11" /> The planet migrates due to disk-planet interactions and reaches ~64 au after 150 orbits, which is very close to the location of the B3 ring. Higher mass planets (here we considered planet-to-star mass ratio of $6 \times 10^{-5}$, i.e., a ~33 $M_{\oplus}$ planet) produce significantly wider gaps even at higher viscosities (a simulation with $\alpha = 10^{-4}$ was performed, not shown here).

The time that best reproduces our observations of the outer rings is approximately at 150 orbits of the mini-Neptune, i.e., ~0.05 Myr. We note that the rings are present since orbit ~100 onward (see Figure 3). The model is stable for at least 1000 planet orbits (~0.4 Myr). After this, a global $m = 1$ mode grows in the disk, leading both the planet and the disk to develop some eccentricities. After this, the co-rotating ring gradually diffuses with time as the gas surface density decreases.

### 3.2.2. RT Calculation

We compute the continuum emission from our dust simulations using the public RT code RADMC3D (version 0.41). Twenty logarithmically spaced bins are used to allocate the 200,000 dust particles from the simulation. The dust sizes ($s$) follow a power-law distribution of $n(s) \propto s^{-3.5}$, with minimum and maximum sizes of 1 µm and 1 mm, respectively. Assuming a gas-to-dust ratio of 33 yields a total dust mass in the outer region of ~100 $M_{\oplus}$.

Each 2D distribution of dust particles is turned into a surface density and expanded in the vertical direction assuming hydrostatic equilibrium. The dust’s scale height $h_{i,d}$ is size-dependent and follows

$$h_{i,d} = h \times \sqrt{\frac{D_z}{D_z + \text{St}_i}}, \quad (1)$$

where $h$ is the initial gas pressure scale-height, $St_i$ is the averaged Stokes number for the $i$th dust size bin, and $D_z$ is a turbulent diffusion coefficient in the vertical direction that depends on the level of turbulent activity across the vertical extent of the disk (Yang et al. 2018). We assume $D_z$ to be proportional to the $\alpha$ turbulent viscosity in our 2D hydrodynamic simulation (see Baruteau et al. 2019 for more details). A $D_z$ value equal to $10 \alpha$ yields dust temperatures that match the observed brightness temperature of rings B2 and B3. In the RT calculation, we use 24 cells logarithmically spaced to sample three scale heights in colatitude (with finer sampling toward the midplane).

We expand the hydrodynamic simulation’s grid to include the inner regions and then artificially add the inner disk and the bright ring B1. Their parameterization follows previous RT modeling, which fits the 1.3 mm data (Fedele et al. 2017). The presence of these structures that block stellar radiation are needed to recover the low temperatures of the outer regions. A 7320 K star with a radius of $1.5 R_{\odot}$ (Gaia Collaboration et al. 2018) is placed at the center of the grid and $10^8$ photon packages are propagated to compute the temperature of each dust bin size via the Monte-Carlo RADMC3D mctherm task. The temperatures range between ~13 and ~20 K in the outer rings, consistent with previous multiwavelength modeling (Fedele et al. 2017). These values are two times higher than the locally isothermal prescription used in the hydrodynamical simulation, which is representative of the gas scale-height, rather than a physical temperature of dust grains. Ray tracing is performed to solve the RT equation for continuum emission assuming the dust is a mix of 30% amorphous carbons and 70% silicates. The resulting synthetic image is convolved with the MEM resolution beam and it is shown in Figure 3. The high fidelity of the ALMA image allows direct comparison with the data in the image plane. The hydrodynamical model reproduces the locations of B2, B3, and B4 as well as the level of flux along rings B2 and B3, and to a lesser extend ring B4’s.

### 3.3. Results

A hydrodynamical model including a $10 M_{\oplus}$ planet embedded in a disk with an aspect ratio of $h = 0.033$ reproduces the observations (Figure 4). The initial condition is such that the planet is introduced at 69 au. After about 50,000 years, the planet migrates inward to a ~64 au orbit,

<footnote n="11">¹¹ The stellar mass of HD 169142 is closer to $M_{\star} = 1.7 M_{\odot}$, which yields a planet mass of $9.6 M_{\oplus}$ in our simulation. Hence we round the mini-Neptune’s mass off to ~10 $M_{\oplus}$ throughout the paper.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Gas and dust spatial distribution of a hydrodynamic simulation with a migrating $10 M_\odot$ (mini-Neptune size) planet. The perturbed gas surface densities relative to the initial profile are shown after 100 (a), 150 (b), 200 (c) and 300 (d) orbits of the mini-Neptune. Dust particles overlaid using colored dots (color varies with particle size, see color bar on the rhs of the panel). The planet is marked by a white circle and can be spot through its spiral density waves. Three dust rings are produced after 100 orbits: one at the planet's orbital radius (~64 au), and two on each side of the planet's gap. The two dust rings on both sides of the planet are not symmetrical w.r.t the planet's orbital radius. The asymmetry increases with time as the planet moves inward (the outermost ring keeps a fairly constant location).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Side-by-side comparison between the ALMA 1.3 mm image of HD 169142 (left; same as in Figure 1) and the predicted image based on the hydrodynamic simulation of a migrating mini-Neptune planet of $10 M_⊕$ (right). After 150 orbits, the triple ring system emerges in the outer disk. The inner disk and ring B1 are not included in the hydrodynamic simulation but added a posteriori at the radiative transfer stage.</figcaption></figure>

shepherding ring B3 and pushing rings B2 and B4, thus explaining the asymmetry in the rings mutual separations. This can be appreciated in Figure 3. A simulation with a planet on a fixed orbit (not included) produces rings that are too equidistant for the disk and planet parameters assumed here. If the $10 M_⊕$ core formed closer to the star and then migrated outward, it would not be consistent with the rings locations: the observed asymmetry reveals the direction of migration. This migration could have distinct observational signatures in multiwavelength observations, for example, by probing the spatial segregation of larger grains with longer wavelength observations. Such predictions have recently been reported by Meru et al. (2019) for a planet migrating in a high-viscosity disk. Despite the little knowledge we have about the initial conditions of the dust density distribution, the model reproduces the observed fluxes along B2 and B3, albeit with a slightly brighter ring B4. The initial amount of large particles near B4's location is likely overestimated in our simulations.

As the planet continues to interact dynamically with the dust reservoir, it also starts clearing an azimuthal opening in its vicinity. This produces a dip in azimuth along the middle ring at the planet's location. Whether this dip is observable at 1.3 mm depends on the planet's mass, the maximum size in the dust distribution, and the elapsed time. The more massive the planet, the faster and deeper is the opening of the azimuthal dip. On the other hand, a large maximum particle size produces a more pronounced dip. Dust growth and fragmentation, which are not taken into account in our simulations, probably also have an effect on the observability of this azimuthal opening. The lack of a clear azimuthal dip along ring B3 indicates that the bulk of particles we observe have sizes <1 mm, the planet is rather small ($\lesssim 10 M_⊕$), or that it may have formed less than a few hundred orbits ago ($\sim 50$ kyr), or most likely a combination

of these options. Note that the signal-to-noise ratio at the rings location might not be sufficient to probe the azimuthal structure along B3. Longer wavelength observations might be able to trace the location of the mini-Neptune as the larger grains pile-up behind the planet (note the millimeter particles in Figure 3).

## 3.4. Discussion

### 3.4.1. Could Rings B2, B3, and B4 Be Shaped by Two Low-mass Planets?

A scenario where two low-mass planets (one for each narrow gap) results in the triple ringed structure is unlikely for several reasons. The small separation between the gaps (only 9 au) requires scale heights at the location of the planets which would be too small for a realistic disk temperature; these translate into <5 K at 64 au (following Equation 16 in Dong & Fung 2017). In other words, if you put two planets into the two gaps, they will likely open one common gap instead of two. Moreover, the similarity between D3 and D4 suggests that if they are produced by two planets, these ought to have formed at similar times and with similar masses. The dynamical interactions between the closely orbiting planets might also make them unstable.

### 3.4.2. How Does an Inner Giant Planet Impact the Structure of Rings B2, B3, and B4?

In recent years, it has been shown that planet perturbations can extend far beyond the Hill radius if the disk viscosity is low (Bae et al. 2018). Since the model presented here assumes a very low viscosity, any inner giant planet (in D2, for example) could have an impact on the outer triple ring system. In a debris disk or a protoplanetary disk with low gas surface density, giant
planets could potentially produce fine structures by mean motion resonances with the dust. However, HD 169142 is gas rich and the dust in the outer region is still coupled to the gas distribution. To further test the possibility that the fine structure is somehow related to inner giant planets, we included a giant planet in our simulation, of roughly a Jupiter mass located in the middle of D2 at 38 au (see Appendix C). The giant indeed opens a deep gap in the gas, but fails to generate any additional narrow rings in the outer disk. However, the spiral wakes excited by the giant and the vortices at the gap edge (produced by the Rossby wave instability) induce eccentricity on the outer triple rings associated with the mini-Neptune (Figure 6). Recent scattered-light imaging indeed shows faint spiral features at the disk surface, which can be associated to putative giant planets in the wide gaps (Gratton et al. 2019). The dust particles around the mini-Neptune can acquire noncircular trajectories either (i) by direct gravitational interaction with the mini-Neptune, whose eccentricity varies on account of the inner giant's wakes depositing energy and angular momentum in the planet's horseshoe region (Figure 7), and/or (ii) by being deflected by the (shock) wakes upon crossing them.

The addition of a giant planet adds significant complexities as the final configuration of rings would strongly depend on the initial conditions for the dust distribution and the precise timing of the formation of the mini-Neptune and the giant planet. As protoplanetary disk observations grow in sensitivity and resolution, more and new ingredients need to be considered in the modeling efforts. This paper is meant as a proof of concept to show that a simple one-planet simulation can account for seemingly complex substructure in the outer ring of a transition disk. Finding a perfect fit to the data is not in the scope of this work. Indeed, some important aspects have been left aside such as: dust growth (e.g., Bae et al. 2018) and fragmentation, dust back-reaction onto the gas (e.g., Gonzalez et al. 2017; Dipierro et al. 2018), dust torque due to scattered pebbles (Benítez-Llambay & Pessah 2018), 3D and MHD effects (e.g., Flock et al. 2015; Miranda et al. 2017), and possibly many more. The inclusion of planet migration was needed here to explain the asymmetry in the rings' mutual separations, a step forward toward understanding how low-mass planet migration can be studied from observations of dust radial structures. However, a proper account of dust dynamics, planet accretion, and thermodynamics, in a self-consistent way, is necessary to build the migration history of these low-mass protoplanets (Benítez-Llambay & Pessah 2018).

### 3.4.3. What Can These Observations Tell Us about the Planet Formation Process?

The low mass of the putative protoplanet in HD 169142, in addition to the lack of clear evidence for structures associated with gravitational instabilities in (abundant) multiwavelength observations of this source, suggests that a bottom-up process such as core accretion could be responsible for its formation. This implies that core accretion can potentially operate at ~65 au within the age of the system (~6 Myr; Grady et al. 2007), possibly assisted by pebble accretion as the outer region banks more than 100 $M_⊕$ worth of pebble-sized solids (subject to uncertainties of the dust opacities). The modeling presented here also shows that the mini-Neptune should have formed well after the giant planets carved gaps D1 and D2. The presence of a giant planet in D2 would have produced a pressure maxima in

the outer region which could have enhanced core accretion via a dust trap (Pinilla et al. 2012).

## 4. Concluding Remarks

The new ALMA observations presented here show that a transition disk with wide, deep gaps can also host narrow-ring structures in its outer region, similar to those observed in the HL Tau and the TW Hya systems.

The HD 169142 observation allows us to link the architecture of protoplanetary disks with low-mass planets. The interpretation via hydrodynamics is a proof of concept that links the structure of closely packed double gaps and triple rings with a single and migrating low-mass planet. Planetary migration naturally explains the distinct mutual separations between the narrow rings. The connection was made possible in HD 169142 thanks to the isolation of its outer region. In the absence of a clear gap that separates an outer ring, the superposition of multiple rings due to several planets hampers simple and clear explanations such as that found for HD 169142.

In HD 169142, we have thus found evidence that suggests that low-mass planet formation can occur in the outer regions of disks bearing evidence for giant protoplanets. The planet formation mechanism, likely core accretion or any bottom-up process, can thus produce planet embryos at ~65 au and outside the orbit of inner giants, at least in certain disks.

We thank Ed Fomalont, Anya Yermakova, and Philipp Weber for useful discussions, as well as the anonymous referee for the constructive comments. Financial support was provided by the government of Chile grants Millennium Scientific Initiative RC130007, CONICYT-Gemini 32130007, and CONICYT-FONDECYT grant No. 1171624, 1171246, and 1191934. S.P. acknowledges support from the Joint Committee of ESO and the Government of Chile. The data analysis and some of the simulations were carried out in the Brelka cluster, hosted at DAS/U. de Chile (Fondequip EQM140101). Numerical simulations were carried out on the CalMip machine of the Centre Interuniversitaire de Calcul de Toulouse, which is gratefully acknowledged. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00344.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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work has made use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Facility: ALMA Observatory.

Software: GPUVMEM (Cárcamo et al. 2018), CASA (McMullin et al. 2007), FARGO-ADSG (Baruteau & Masset 2008a, 2008b), RADMC-3D (Dullemond et al. 2012).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** HD 169142 image synthesis. CASA tclean with Briggs weighting, with a robust parameter of 0.0, yielding a resolution of 39 × 23 mas (a), uniform weight with robust −2.0 and resolution 32 × 18 mas (b), and the GPUVMEM deconvolution (c). The effective resolution of the MEM deconvolved image is 27 × 20 mas and is measured directly from a 2D Gaussian fit to the stellar component in the center of the field. The data used to create this figure are available.</figcaption></figure>

# Appendix A
## Comparison between CLEAN and MEM Algorithm

A side-by-side comparison of HD 169142 image synthesis using CLEAN and MEM algorithms is shown in Figure 5. The triple ring structure in the outer region is robust in all image reconstructions. The uniform weighting CLEAN image has a comparable resolution to MEM but at the cost of decreasing sensitivity.

# Appendix B
## Polar Deprojection and Ellipse Fitting

To study the eccentricity of the rings, we extract the radial profile of each ring as a function of azimuth, $R(\theta)$, by fitting Gaussian functions in radius (using the PYTHON iminuit optimization library). Only signal above $5\sigma$ is included in the process. Ring B4 is not considered as it only has a small number of points with high signal-to-noise. B1 is also left aside as its perturbed morphology cannot be reproduced by a simple elliptical curve. The optimization yields the values for $R(\theta)$ and their error bars, $\Delta R(\theta)$, for each ring. The ellipse equation is then fitted to $R(\theta)$ for a given ring using a direct least-square procedure (Fitzgibbon et al. 1996). The fitting procedure yields the ellipse center, its angle of rotation, and its eccentricity ($e = \sqrt{1 - (b/a)^2}$, where *a* and *b* are the semimajor and semiminor axes, respectively). The best-fit ellipse for each ring are plotted on the polar projection in Figure 2. The uncertainties in the eccentricity values are calculated using a Monte-Carlo approach where data points are modified by random numbers of the order of $3\Delta\theta$ and then fitted following the aforementioned procedure. This is repeated a thousand times and the final errors are drawn from the standard deviation. Correlated noise is accounted for by multiplying the centroids' error by a factor, $\sqrt{N_{pix}}$, where $N_{pix} = 5$ is the number of pixels along the synthesized beam's major axis.

The uncertainty of the eccentricity is dominated by our imprecise knowledge of the inclination angle of the disk. In order to determine the significance of the fitted values, we explore a range of inclination angles ranging between 9° and 16°, repeating the fitting procedure described above for each

inclination value. The results are shown in Figure 2(b). Rings B2 and B3 show eccentricities which indeed vary with inclination angle, with a minimum eccentricity of $e \approx 0.085$. Interestingly, the rings have different eccentricities for a given inclination, except at $i \approx 12^\circ.5$ where both rings have the same eccentricity. This inclination angle is close to the standard value of $i \sim 13^\circ$ used in previous multiwavelength modeling (e.g., Fedele et al. 2017; Bertrang et al. 2018). The intrinsic position angles of the ellipses also coincide at $i \approx 12^\circ.5$. This fitting procedure thus suggests a common inclination for both B2 and B3 of $i \approx 12^\circ.5 \pm 0.5$ deg (1$\sigma$). At this inclination, the eccentricity of B2 and B3 share the common value of $0.09 \pm 0.02$ (2$\sigma$ uncertainty; same as the shaded area in Figure 2). A disk inclination angle different from $12^\circ.5$ yields different eccentricity values for B2 and B3, which would suggest the occurrence of warping in the outer regions.

# Appendix C
## Addition of a Giant Planet Associated with Gap D2

The impact that an inner giant could have on the outer narrow rings is explored by including a 0.5 $M_{Jup}$ planet fixed at 38 au, in the middle of the observed gap D2, in addition to the mini-Neptune. The general setup inherits the same parameters as the single mini-Neptune simulation, with some changes. The initial radius of the orbit of the mini-Neptune planet is decreased to 68 au, this is to account for a slightly lower migration rate in the presence of the inner giant. The grid resolution is lowered to 600 × 900 cells in radius and azimuth, respectively. The grid's inner edge is extended down to 20 au as to include gap D2 in the radial domain. The smaller inner edge translates into a significant computational cost. The number of particles is decreased by half (i.e., only 100,000 particles are used) so as to make the computation less expensive. Dust particles are distributed over the same region as in the previous simulation. Dust around the giant planet is not included as this would require increasing the number of particles, making the simulation too expensive. As mentioned in the main text, fitting the already known features B1 and D1 is beyond the scope of this work. Modeling all features in the HD 169142 system requires an extensive search of a large
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Effect of an inner giant planet on the triple ring structure over 100 kyr. The figure shows the perturbed gas density and the distribution of dust particles after 100 (left), 150 (middle), and 300 (right) orbits of the mini-Neptune. The onset of the ring structure around the mini-Neptune happens near 100 orbits, while at 300 orbits (100 kyr) the ring B3 develops an azimuthal asymmetry, similar as in the case without the giant planet. In the presence of the giant, the outer rings associated with the mini-Neptune interact with the giant's wakes, acquiring some eccentricities. The gas distribution in the outer region also becomes eccentric.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Orbital radius (a) and eccentricity (b) of the migrating mini-Neptune in the presence of an inner giant planet. The amplitude of the oscillations in eccentricity and orbital radius reflects the impact of the giant planet's wakes on the mini-Neptune's orbit. The planet's initial eccentricity arises because of the disk's gravity being felt by the planet in addition to that of the star. The planet's eccentricity is calculated assuming a two-body problem with only the star and the planet, not the disk. The value of the planets eccentricity thus reflects the disk-to-star mass ratio.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison between polar deprojection of the ALMA observation (a) and synthetic predictions from hydrodynamic models after 150 orbits of the mini-Neptune. Panel (b) shows the result of the hydrodynamic simulations for a single mini-Neptune (same as Figure 4), while panel (c) shows the mini-Neptune's outer rings under the influence of an inner giant planet at 38 au. All panels have linear color stretch between 0.0 and 0.25 mJy/beam.</figcaption></figure>
parameter space, but the main difficulty lies on the little
knowledge available of the initial conditions and the relative
age of each planet. In our simple two-planet setup, both
companions are introduced simultaneously at the beginning of
the simulation.

Figure 6 shows that the presence of a giant planet at 38 au carves a wide gap in the gas at the location of the D2 dust gap. At the same time, the mini-Neptune shapes the outer region into the triple ring system, consistent with the single-planet simulation. The presence of the giant planet makes the morphology of the outer rings appear eccentric. The asymmetry in the outer rings mutual separation is also reproduced. The migration of the mini-Neptune is slightly affected by the gas giant's wakes (see Figure 7). Vortices develop at the outer edge of the giant's gap but these start decaying shortly after a couple of hundred orbits of the mini-Neptune (see the right panel in Figure 6). The lack of observational evidence for these vortices may suggest that the mini-Neptune formed long after the giant planet carved the gap D2, at a time where the vortices had already decayed. A side-by-side comparison of the observations with the giant planet simulation, including the simulation for a single mini-Neptune, is shown in Figure 8.

ORCID iDs

Sebastián Pérez https://orcid.org/0000-0003-2953-755X

Simon Casassus https://orcid.org/0000-0002-0433-9840

Clément Baruteau https://orcid.org/0000-0002-2672-3456

Ruobing Dong https://orcid.org/0000-0001-9290-7846

Antonio Hales https://orcid.org/0000-0001-5073-2849

Lucas Cieza https://orcid.org/0000-0002-2828-1153

References

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

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

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

Bae, J., Pinilla, P., & Birnstiel, T. 2018, ApJL, 864, L26

Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201

Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJL, 815, L15

Baruteau, C., Barraza, M., Pérez, S., et al. 2019, MNRAS, 486, 304

Baruteau, C., & Masset, F. 2008a, ApJ, 672, 1054

Baruteau, C., & Masset, F. 2008b, ApJ, 678, 483

Baruteau, C., & Zhu, Z. 2016, MNRAS, 458, 3927

Benítez-Llambay, P., & Pessah, M. E. 2018, ApJL, 855, L28

Bertrang, G. H.-M., Avenhaus, H., Casassus, S., et al. 2018, MNRAS, 474, 5105

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

Boss, A. P. 1997, Sci, 276, 1836

Cárcamo, M., Román, P. E., Casassus, S., et al. 2018, A&C, 22, 16

Charnoz, S., Fouchet, L., Aleon, J., et al. 2011, ApJ, 737, 33

Chen, L., Kóspál, Á., Ábrahám, P., et al. 2018, A&A, 609, A45

Clarke, C. J., Tazzari, M., Juhasz, A., et al. 2018, ApJL, 866, L6

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

Dipierro, G., Laibe, G., Alexander, R., et al. 2018, MNRAS, 479, 4187

Dipierro, G., Laibe, G., Price, D. J., et al. 2016, MNRAS, 459, L1

Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73

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

Dong, R., Li, S., Chiang, E., et al. 2017, ApJ, 843, 127

Dong, R., Li, S., Chiang, E., et al. 2018a, ApJ, 866, 110

Dong, R., Liu, S.-y., Eisner, J., et al. 2018b, ApJ, 860, 124

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

Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJL, 869, L46

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A Multi-purpose Radiative Transfer Tool, ascl:1202.015

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

Fitzgibbon, A. W., Pilu, M., & Fisher, R. B. 1996, ICPR, 1, 253

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

Fuente, A., Baruteau, C., Neri, R., et al. 2017, ApJL, 846, L3

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Goldreich, P., Lithwick, Y., & Sari, R. 2004, ARA&A, 42, 549

Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425

Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Gonzalez, J.-F., Pinte, C., Maddison, S. T., et al. 2012, A&A, 547, A58

Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793

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

Gratton, R., Ligi, R., Sissa, E., et al. 2019, A&A, 623, A140

Honda, M., Maaskant, K., Okamoto, Y. K., et al. 2012, ApJ, 752, 143

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42

Kley, W., & Dirksen, G. 2006, A&A, 447, 369

Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

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

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Macías, E., Espaillat, C. C., Ribas, Á., et al. 2018, ApJ, 865, 37

Masset, F. 2000, A&AS, 141, 165

McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA), 127

Mentitplay, D., Price, D. J., & Pinte, C. 2018, MNRAS, 484, L130

Meru, F., Rosotti, G. P., Booth, R. A., et al. 2019, MNRAS, 482, 3678

Miranda, R., Li, H., Li, S., et al. 2017, ApJ, 835, 118

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

Okuzumi, S., Momose, M., Sirono, S.-i., et al. 2016, ApJ, 821, 82

Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43

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

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

Pohl, A., Benisty, M., Pinilla, P., et al. 2017, ApJ, 850, 52

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

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

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

Rosotti, G. P., Juhasz, A., Booth, R. A., et al. 2016, MNRAS, 459, 2790

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., et al. 2018, MNRAS, 477, 1239

van der Marel, N., Dong, R., di Francesco, J., et al. 2019, ApJ, 872, 112

Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27

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

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47