# An Ideal Testbed for Planet–Disk Interaction: Two Giant Protoplanets in Resonance Shaping the PDS 70 Protoplanetary Disk

Jaehan Bae, Zhaohuan Zhu, Clément Baruteau, Myriam Benisty, Cornelis

Dullemond, Stefano Facchini, Andrea Isella, Miriam Keppler, Laura Perez, Richard Teague

► To cite this version:

Jaehan Bae, Zhaohuan Zhu, Clément Baruteau, Myriam Benisty, Cornelis Dullemond, et al. An Ideal Testbed for Planet–Disk Interaction: Two Giant Protoplanets in Resonance Shaping the PDS 70 Protoplanetary Disk. The Astrophysical Journal Letters, 2019, 884 (2), pp.L41. 10.3847/2041-8213/ab46b0. hal-02369599

HAL Id: hal-02369599

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

Submitted on 1 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

# An Ideal Testbed for Planet–Disk Interaction: Two Giant Protoplanets in Resonance Shaping the PDS 70 Protoplanetary Disk

Jaehan Bae<sup>1,13</sup>①, Zhaohuan Zhu<sup>2</sup>①, Clément Baruteau<sup>3</sup>①, Myriam Benisty<sup>4,5,6</sup>①, Cornelis P. Dullemond<sup>7</sup>①,
Stefano Facchini<sup>8</sup>①, Andrea Isella<sup>9</sup>①, Miriam Keppler<sup>10</sup>①, Laura M. Pérez<sup>4</sup>⑤, and Richard Teague<sup>11,12</sup>⑤

## Abstract

While numerical simulations have been playing a key role in the studies of planet–disk interaction, testing numerical results against observations has been limited so far. With the two directly imaged protoplanets embedded in its circumstellar disk, PDS 70 offers an ideal testbed for planet–disk interaction studies. Using two–dimensional hydrodynamic simulations we show that the observed features can be well explained with the two planets in formation, providing strong evidence that previously proposed theories of planet–disk interaction are in action, including resonant migration, particle trapping, size segregation, and filtration. Our simulations suggest that the two planets are likely in 2:1 mean motion resonance and can remain dynamically stable over million-year timescales. The growth of the planets at $10^{-8}-10^{-7} M_{\text{Jup}} \text{ yr}^{-1}$, rates comparable to the estimates from $H\alpha$ observations, does not destabilize the resonant configuration. Large grains are filtered at the gap edge and only small, (sub-)µm grains can flow to the circumplanetary disks (CPDs) and the inner circumstellar disk. With the submillimeter continuum ring observed outward of the two directly imaged planets, PDS 70 provides the first observational evidence of particle filtration by gap-opening planets. The observed submillimeter continuum emission at the vicinity of the planets can be reproduced when (sub-)µm grains survive over multiple CPD gas viscous timescales and accumulate therein. One such possibility is if (sub-)µm grains grow in size and remain trapped in pressure bumps, similar to what we find happening in circumstellar disks. We discuss potential implications to planet formation in the solar system and mature extrasolar planetary systems.

*Unified Astronomy Thesaurus concepts: Circumstellar disks (235); Planet formation (1241); Hydrodynamical simulations (767)*

## 1. Introduction

Numerical simulations have been playing a key role in the studies of planet–disk interaction and planet formation, by confirming analytic theories, examining nonlinear phenomena, allowing us to explore a broad parameter space, and making observational predictions. Despite the important role, testing results from numerical simulations against observations has been limited because observing planets in formation has been very challenging. The situation is, however, gradually changing. Thanks to increasingly powerful observing facilities and techniques we are now able to peer into the birthplaces of planets, routinely finding disk substructures hinting at ongoing planet formation (e.g., Andrews et al. 2018; Avenhaus et al. 2018; Pinte et al. 2018; Teague et al. 2018), although whether or not planets are indeed the cause of the observed features has to be further investigated.

With two directly detected protoplanets embedded in its circumstellar disk (Keppler et al. 2018; Wagner et al. 2018; Haffert et al. 2019), PDS 70 offers an ideal testbed for planet–disk

interaction studies. The two planets, PDS 70b and c, are located 195 and 234 mas (Keppler et al. 2018; Haffert et al. 2019) from the 5.4 Myr old, K7 pre-main-sequence star PDS 70 (Müller et al. 2018) on the sky. The deprojected distances between the planets and the star (~20 and 35 au) suggest that the two planets are in or near 2:1 mean motion resonance (Haffert et al. 2019). They cleared their vicinity and opened an inner cavity in the circumstellar disk, initially identified in IR observations (Dong et al. 2012; Hashimoto et al. 2012). Submillimeter continuum observations have also revealed an inner cavity with the emission concentrated beyond the two planets’ orbits (Hashimoto et al. 2015; Long et al. 2018; Keppler et al. 2019). Rotation velocity measurements of CO gas revealed that the gas pressure changes over radius in a way suggesting that grains are trapped in a pressure bump (Keppler et al. 2019), although it is pointed out that the cavity is too wide to explain with PDS 70b alone, the only planet known to be in the system at that time. The size of the inner cavity is larger in the submillimeter continuum (~74 au; Keppler et al. 2019) than in IR (~52 au with the Gaia DR2 distance of 113 pc; Dong et al. 2012). IR and submillimeter observations

<footnote n="1">¹ Department of Terrestrial Magnetism, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA; jbae@carnegiescience.edu</footnote>

<footnote n="2">² Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154, USA</footnote>

<footnote n="3">³ IRAP, Université de Toulouse, CNRS, UPS, F-31400 Toulouse, France</footnote>

<footnote n="4">⁴ Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="5">⁵ Unidad Mixta Internacional Franco-Chilena de Astronomía (CNRS, UMI 3386), Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="6">⁶ Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="7">⁷ Zentrum für Astronomie, Heidelberg University, Albert Ueberle Str. 2, D-69120 Heidelberg, Germany</footnote>

<footnote n="8">⁸ European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, Rice University, 6100 Main Street, MS-108, Houston, TX 77005, USA</footnote>

<footnote n="10">¹⁰ Max Planck Institute for Astronomy, Königstuhl 17, D-69117, Heidelberg, Germany</footnote>

<footnote n="11">¹¹ Department of Astronomy, University of Michigan, 311 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA</footnote>

<footnote n="12">¹² Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote>Received 2019 August 20; revised 2019 September 18; accepted 2019 September 20; published 2019 October 15</footnote>

<footnote n="13">¹³ NHFP Sagan Fellow.</footnote>
suggest that both planets are surrounded by a dusty circumplanetary disk (CPD; Christiaens et al. 2019; Isella et al. 2019), while the nature of the submillimeter emission detected near PDS 70b, which is 74 mas offset from the location of PDS 70b inferred from Hα/IR emission, is yet to be explained (Isella et al. 2019).

Interestingly enough, many of the aforementioned observed features have been previously seen in and proposed to happen in protoplanetary disks by numerical simulations. Multiple planets in resonance are proposed as a possible cause of inner cavities seen in (pre)transitional disks (Zhu et al. 2011). In the pressure bump forming beyond a planet's orbit, particles can be efficiently trapped, and we expect segregation of grain sizes such that, under typical protoplanetary disk conditions, smaller grains are distributed over a larger radial extent, resulting in a smaller cavity size (Pinilla et al. 2012; Zhu et al. 2012). The outer edge of the gap can act as a filter so small particles can penetrate the planet-induced gap, whereas large particles are filtered and remain trapped beyond the gap (Rice et al. 2006; Zhu et al. 2012). The existence of CPDs around giant planets is also conceptually and numerically predicted (Quillen & Trilling 1998; Lubow et al. 1999; Ayliffe & Bate 2009; Ward & Canup 2010), as a consequence of angular momentum conservation similar to the formation of circumstellar disks around protostars.

In this Letter, we carry out two-dimensional hydrodynamic calculations adopting physical properties of the PDS 70 circumstellar disk and the two planets embedded therein. By simulating the dynamics of disk gas and dust in response to two growing planets, we show that the signposts of planet-disk interaction predicted by numerical simulations and the observed features of the PDS 70 disk show a good agreement with each other. We believe the resemblance between simulations and observations strongly supports the previously proposed theories, including resonant migration, particle trapping, size segregation, and filtration, are in action.

The Letter is organized as follows. In Section 2, we introduce our numerical model. In Section 3, we present results from numerical simulations, focusing on the evolution of the circumstellar disk and planets' orbits. In Section 4, we generate simulated submillimeter continuum images from our simulations to compare with the observations and discuss how our findings can help better understand planet-disk interaction, CPDs, and planet formation. We present our conclusions in Section 5.

## 2. Model Setup

### 2.1. Initial Gas Disk

We adopt an initial disk gas surface density profile that falls off with an exponential tail against radius R:

$$ \Sigma_{\text{gas,init}}(R) = \Sigma_c \left(\frac{R}{R_c}\right)^{-1} \exp\left(-\frac{R}{R_c}\right), \quad (1) $$

where we choose $R_c = 40$ au and $\Sigma_c = 2.7$ g cm$^{-2}$ so that the total disk gas mass within the simulation domain is 0.003 $M_\odot$, similar to the model used in Keppler et al. (2019).

To set up the initial disk temperature profile, we iterate Monte Carlo radiative transfer (MCRT) calculations using RADMC-3D (Dullemond et al. 2012) until we obtain a converged three-dimensional disk density and temperature profile (see Appendix A). The least-squares power-law fit to the

density-weighted, vertically integrated temperature $\bar{T}$ is

$$ \bar{T}(R) = 44 \text{ K} \left(\frac{R}{22 \text{ au}}\right)^{-0.24}. \quad (2) $$

Adopting a stellar mass of 0.85 $M_\odot$ (Keppler et al. 2019) and a mean molecular weight of 2.4, this temperature profile corresponds to a disk aspect ratio $H/R$ of

$$ H/R = 0.067 \left( \frac{R}{22 \text{ au}} \right)^{0.38}. \quad (3) $$

### 2.2. Planets

We fix PDS 70b's mass to 5 $M_{\text{Jup}}$, consistent with previous estimates (Keppler et al. 2018; Müller et al. 2018). We test three different masses for PDS 70c: 2.5, 5, and 10 $M_{\text{Jup}}$. The two planets are initially placed at 20 and 35 au, in an agreement with the observed deprojected distances to the central star (Haffert et al. 2019; Keppler et al. 2019). With the semimajor axes the orbital period ratio is 2.3:1, so the two planets are located slightly outside of 2:1 mean motion resonance initially. In addition to the three simulations, we carry out a simulation with PDS 70b only.

We run simulations for 2 Myr, which corresponds to about 20,000 orbits at PDS 70b's initial radial location. We linearly increase planet masses over the first $10^4$ yr, while fixing their orbits. After $10^4$ yr, we allow the planets to gravitationally interact with each other as well as with the circumstellar disk. As we will show below, the two planets settle into 2:1 mean motion resonance within about 0.1 Myr in all three cases, and the common gap opened by the planets reaches a quasi-steady-state by 0.2 Myr.

To examine whether increases in planet masses affect the stability of the system, we allow the planets to accrete gas starting at 0.2 Myr. In practice, the gas density in cells within a fraction of planets' Hill radius is reduced by a fraction $f_{\text{acc}}\Omega\Delta t$ and added to the planets each hydrodynamic time step, following the approach presented in Kley (1999) and Dürmann & Kley (2017). This results in the gas depletion timescale of $\tau_{\text{acc}} = (f_{\text{acc}}\Omega)^{-1}$. We choose $f_{\text{acc}} = 0.01$, with which the planets accrete at $10^{-8}-10^{-7} M_{\text{Jup}} \text{ yr}^{-1}$, consistent with the estimates based on Hα observations (Wagner et al. 2018; Haffert et al. 2019). Note that this parameterized planet accretion is adopted to examine the dynamical stability in a way that the total mass and momentum are conserved, rather than to realize actual accretion processes within planets' Hill sphere.

The gravitational potential of planets is softened over 60% of the local gas scale height, in order to mimic the overall magnitude of the torque in three dimensions (Müller et al. 2012).

### 2.3. Hydrodynamic Simulations

We carry out two-dimensional, locally isothermal hydrodynamical simulations using the Dusty FARGO-ADSG code (Baruteau et al. 2019). This is an extended version of the publicly available FARGO-ADSG (Masset 2000; Baruteau & Masset 2008a, 2008b), with Lagrangian test particles implemented (Baruteau & Zhu 2016).

The simulation domain extends from 2.2 to 198 au in the radial direction and covers the entire $2\pi$ in azimuth. We adopt 672 logarithmically spaced grid cells in the radial direction and 936 uniformly spaced grid cells in the azimuthal direction. At the radial boundaries we adopt a wave-damping zone
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (Left) The two-dimensional gas surface density distribution, (middle) particle distribution, and (right) azimuthally averaged radial gas and dust density distributions at *t* = 0.6 Myr. From top to bottom, models without PDS 70c, with *M*<sub>c</sub> = 2.5 *M*<sub>Jup</sub>, *M*<sub>c</sub> = 5 *M*<sub>Jup</sub>, and *M*<sub>c</sub> = 10 *M*<sub>Jup</sub>. Dashed circles in the left and middle panels show the location of the azimuthally averaged gas pressure maximum. In the middle panels, the red circles show the size of the Hill radius of the planets. In the right panel, the red dashed line shows 1:2 outer Lindblad resonance location of PDS 70c (PDS 70b in the model without PDS 70c), while the black dashed line shows the gas pressure maximum. The dotted curve shows the initial gas surface density profile. The color histograms show the azimuthally averaged dust surface density with sizes of (dark blue) 0.1–1 µm, (green) 1–10 µm, (red) 10–100 µm, and (yellow) 0.1–1 mm, adopting 1 au bin size in radius.</figcaption></figure>
(de Val-Borro et al. 2006) to suppress wave reflection, from 2.2 to 2.64 au and from 176 to 198 au. Because we are interested in long-term, Myr timescale evolution, we decrease the surface density in the wave-damping zone $\Sigma_{gas,damp}$ over the local viscous timescale following

$$ \Sigma_{gas,damp} = \Sigma_{gas,init} \exp(-t/t_{vis}). \quad (4) $$

Here, $t_{vis} = R^2/\nu$ is the viscous timescale and $\nu = \alpha c_s^2/\Omega$, where $\alpha$ is the viscosity parameter, $c_s$ is the sound speed, and $\Omega$ is the orbital frequency. A uniform disk viscosity $\alpha = 10^{-3}$ is applied.

In addition to the gas component we insert 100,000 Lagrangian test particles at $t = 0.5$ Myr, well after the overall disk structure reaches a quasi-steady-state. Test particles are inserted between 50 and 100 au, with a uniform dust-to-gas mass ratio of 2.5% across this radial region. This results in a total dust mass of about $10 M_{\oplus}$. We assume a dust bulk density of $1.26 \text{ g cm}^{-3}$, which corresponds to that of aggregates with 30% silicate matrix and 70% water ice. Particle sizes are determined such that we have approximately the same number of test particles per decade of a size between 0.1 $\mu$m and 1 mm. We choose the maximum particle size of 1 mm because the fragmentation-limited maximum grain size (e.g., Birnstiel et al. 2012) between 50 and 100 au is a few hundred $\mu$m with the disk and particle properties we adopt. With the initial gas surface density in Equation (1), the Stokes number of the test particles can be expressed as

$$ St \approx 0.07 \left( \frac{s}{1 \text{ mm}} \right) \left( \frac{R}{40 \text{ au}} \right) \exp \left( \frac{R}{40 \text{ au}} \right). \quad (5) $$

Test particles feel the gravity of the star and the planets. In addition, they interact with the circumstellar disk gas via aerodynamic drag. Turbulent diffusion is included as stochastic kicks on the particles' position following the method presented in Charnoz et al. (2011), adopting $\alpha = 10^{-3}$. Test particles do not provide feedback onto the planets and the disk gas. The size evolution of particles is not included in the simulations.

Since we aim to explain the submillimeter continuum flux associated with the CPDs (Section 4.2), we assign mass to test particles so that we can keep track of the CPD dust mass. In practice, this mass assignment is done such that the dust mass at each radius in the initial disk is distributed over a range of dust size $s$, from 0.1 $\mu$m to 1 mm, to have the mass per interval in $\log(s)$ be proportional to $s^{0.5}$ (this corresponds to a dust size distribution of $n(s) \propto s^{-3.5}$).

## 3. Results

### 3.1. Disk Evolution

We start by discussing the overall circumstellar disk evolution. In Figure 1, we present the two-dimensional gas surface density distribution, test particle distribution, and azimuthally averaged radial distributions of the gas and dust surface density at $t = 0.6$ Myr. When only PDS 70b exists in the disk, the planet's outer gap edge becomes eccentric ($e \sim 0.1$) because of the eccentric Lindblad resonance (Lubow 1991a, 1991b; Papaloizou et al. 2001). The wave modes from the circular component of the planet's potential are excited at $R_{LR,circ} = [(m \pm 1)/m]^{2/3} R_p$, where $m$ denotes the azimuthal wavenumber of the wave modes. The outer Lindblad resonance located farthest away from the planet is the $m = 1$ mode, which occurs at $1.6 R_p \approx 32$ au. As

<figure><img src="image_4.png" /><figcaption>Figure 2. Time evolution of (a) semimajor axes, (b) orbital period ratio, (c) resonant angle, (d) orbital eccentricities, and (e) accretion rates onto the planets. Note that the two planets settle into 2:1 mean motion resonance within the first 0.1 Myr and remain dynamically stable for the duration of the simulation. The first 1 Myr of the simulations are shown for visualization purposes only.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Same as Figure 2, but results with (left panels) $M_c = 5 M_{Jup}$ and (right panels) $M_c = 10 M_{Jup}$.</figcaption></figure>

shown in the radial gas density profile in Figure 1, PDS 70b opens a wide gap with the gas density peaking at about 45 au. The total angular momentum exchange via the circular component of the planet's potential, which is the summation of the individual contribution over the entire azimuthal wavenumbers, is therefore

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

significantly reduced. On the other hand, wave modes from the eccentric component of the planet's potential launch at $R_{LR,ecc} = [(m \pm l)/m]^{2/3} R_p$ ($l$ is an integer greater than 1), which is always beyond the Lindblad resonance of their counterpart circular component $R_{LR,circ}$ in the outer disk. Thus, the overall
angular momentum exchange can be dominated by the eccentric component, making the outer gap edge eccentric. In the PDS 70b- only model particles with sizes 0.1–1 mm (hereafter submillimeter particles), which dominate the submillimeter continuum flux, are trapped in the gas pressure peak at 45 au; this is insufficient to explain the continuum peak at ~74 au in the submillimeter observation (Keppler et al. 2019).

Muley et al. (2019) recently showed that accreting planets can have an abrupt increase in its orbital eccentricity to $e \sim 0.25$ as they grow in mass. These planets can carve a wider gap than they otherwise would, helping explain the large continuum cavity. We, however, do not observe such an increase in orbital eccentricity, at least for the duration of our simulations. We conjecture this could be because the reported eccentricity growth prefers large accretion rates. In the fiducial simulation of Muley et al. (2019), the time-averaged accretion rate until the onset of the eccentricity growth is $\sim 2.5 M_{Jup}/3.5 \text{ Myr} \simeq 7 \times 10^{-7} M_{Jup} \text{ yr}^{-1}$, more than an order of magnitude larger than the accretion rate estimates from observations.

When both PDS 70b and c are inserted the two planets open a common gap. The gas disk and particle ring at the common gap outer edge remain circular when PDS 70c's mass is 2.5 or 5 $M_{Jup}$. In the two models, submillimeter particles are trapped in the pressure peak at ~63 au, which is further out compared with the PDS 70b-only model and is in a better agreement with the continuum ring location in the submillimeter observation (Keppler et al. 2019). The gas surface density around PDS 70c is larger than that around PDS 70b when $M_c = 2.5 M_{Jup}$, while it is comparable to the gas density around PDS 70b when $M_c = 5 M_{Jup}$. This, together with the fact that the CO emission is significantly more depleted around PDS 70b's orbit (Keppler et al. 2019), suggests that PDS 70c likely has a smaller mass than PDS 70b.

When PDS 70c's mass is 10 $M_{Jup}$ the outer gap edge in the gas disk and particle ring become eccentric, with the eccentricity varying between 0.2 and 0.4 over time. Similar to the PDS 70b-only model, the nonzero disk eccentricity arises as the gap is sufficiently wide such that the eccentric component of PDS 70c's potential dominates over its circular component. Because such a large continuum ring eccentricity of $e \sim 0.2–0.4$ can be ruled out by submillimeter observations (Keppler et al. 2019), we conclude that PDS 70c's mass has to be smaller than 10 $M_{Jup}$.

## 3.2. Planets' Orbital Evolution

We plot orbital elements of the two planets in Figures 2 and 3. As can be seen from the resonant angle, defined as $\psi_{cb} = 2\lambda_b - \lambda_c - \varpi_c$, where $\lambda_b$ and $\lambda_c$ are mean longitudes of the planets, and $\varpi_c$ is the outer planet's longitude of perihelion, the two planets migrate toward each other and settle into 2:1 mean motion resonance within the first 0.1 Myr of the simulations for all three cases. We find that the two planets remain dynamically stable for the next 2 Myr. Such a rapid adjustment in their orbits into a resonant configuration suggests that PDS 70b and c are likely in 2:1 mean motion resonance. The exact period ratio is slightly larger than 2:1 because the planets experience repulsion as they interact with each other's spiral arms (Baruteau & Papaloizou 2013).

When PDS 70c's mass is 2.5 $M_{Jup}$ the two planets migrate outward (Figure 2(a)). Between 0.2 and 2 Myr, PDS 70b and c's semimajor axes increase at 0.30 and 0.51 au $Myr^{-1}$, respectively. On the other hand, when PDS 70c's mass is

comparable to or larger than PDS 70b's mass ($M_c = 5 M_{Jup}$ and 10 $M_{Jup}$ models), the planets experience an inward migration (Figure 3). When PDS 70c's mass is 5 $M_{Jup}$, PDS 70b and c's semimajor axes decrease at 0.39 and 0.73 au $Myr^{-1}$. When PDS 70c's mass is 10 $M_{Jup}$, PDS 70b and c's semimajor axes decrease at 0.27 and 0.43 au $Myr^{-1}$, respectively.

In all three models PDS 70c's orbital eccentricity remains relatively small, less than about 0.05, whereas PDS 70b's orbital eccentricity converges to $e \sim 0.05–0.2$ with the exact value dependent upon PDS 70c's mass. Note that a smaller orbital eccentricity for the outer planet is known to be a generic feature of 2:1 mean motion resonance (Baruteau & Papaloizou 2013). Even the largest PDS 70b's orbital eccentricity we find in our simulations ($e \sim 0.2$) cannot be ruled out with the existing observations (Müller et al. 2018).

The accretion rates onto the planets gradually decrease over time. Our simulations suggest that the mass growth at rates consistent with recent $H\alpha$ observations ($10^{-8}–10^{-7} M_{Jup} \text{ yr}^{-1}$; Wagner et al. 2018; Haffert et al. 2019; Thanathibodee et al. 2019) is unlikely to destabilize the mean motion resonance.

# 4. Discussion

## 4.1. Particle Trapping, Filtration, and Size Segregation

Giant planets open a deep gap and trap large grains at the gap edge (e.g., Paardekooper & Mellema 2004; Fouchet et al. 2007; Zhu et al. 2012). The ratio between the Stokes number and the disk viscosity parameter $St/\alpha$ is important in determining the width of particle distribution in a pressure bump and whether or not particles penetrate the gap (Zhu et al. 2012; Dullemond et al. 2018): large particles having Stokes number $St \gtrsim \alpha$ are efficiently filtered at the gap edge, whereas small particles with $St \lesssim \alpha$ are mixed well with gas and follow the gas distribution.

In Figure 4 we show the trajectories of test particles from the $M_c = 2.5 M_{Jup}$ model. As shown, small particles with sizes $\lesssim 10 \mu\text{m}$ have a broad radial distribution, whereas large particles with sizes $\gtrsim 10 \mu\text{m}$ move toward the pressure bump aided by aerodynamic drag and remain trapped within a narrower radial region. While large grains incapable of penetrating the gap establish a quasi-steady-state radial distribution rapidly, note that small, sub-$\mu\text{m}$ grains<fnref n="14" /> are well coupled with gas so leak gradually over time, flowing into the inner disk. It is also interesting to note that some sub-$\mu\text{m}$ grains are captured in CPDs while penetrating the gap (see also Figure 1). The rate at which small particles penetrate the gap appears to be dependent on PDS 70c's mass. Among the three models with two planets, we find that the inner disk is fed with small dust most efficiently when $M_c = 2.5 M_{Jup}$. We conjecture this is because the outward resonant migration facilitates the interaction between PDS 70c and the dust reservoir in the outer disk.

Figures 5(a)–(c) show the observed continuum image at 855 $\mu\text{m}$ and simulated continuum images based on the particle distribution of the PDS 70b-only model and the $M_c = 2.5 M_{Jup}$ model presented in Figure 1 (see Appendix B for details about simulated observations). As shown, the $M_c = 2.5 M_{Jup}$ model reproduces the flux and morphology of the outer continuum ring reasonably well, but the continuum ring in the PDS 70b-only

<footnote n="14">¹⁴ The exact size of grains available to penetrate a gap can differ in other disks depending upon various disk properties, including the gas surface density, the gap depth, and the level of disk turbulence, as well as the grain internal density (Zhu et al. 2012).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** An overlay of the trajectories of test particles in the $M_c = 2.5 M_{Jup}$ model: (top left) 0.1–1 µm, (top right) 1–10 µm, (bottom left) 10–100 µm, and (bottom right) 0.1–1 mm particles. At the gas pressure peak (63 au) particles have Stokes numbers $St \simeq 0.8$ (s/1 mm). Note that only small particles with sizes $\lesssim 1$ µm penetrate the common gap opened by the planets. Note also that some of the small particles penetrating the gap are captured in the circumplanetary disks.</figcaption></figure>

model locates closer to the center of the system compared with the observation. With the submillimeter continuum ring observed outward of the two directly imaged planets, PDS 70 system provides the first observational evidence that gap-opening planets can trap particles beyond their orbits.

Observations of transitional disks at different wavelengths, probing different regions in the disk and/or grains with different sizes, show varying cavity size. In general, the inner cavity seen in molecular gas lines and optical/IR scattered light observations is smaller in size than seen in (sub-)millimeter continuum observations (e.g., van der Marel et al. 2016). This is consistent with what is seen in PDS 70 observations (Dong et al. 2012; Hashimoto et al. 2012; Keppler et al. 2019) and the gas and particle distribution seen in our simulations (Figures 1 and 4).

We note that, similar to PDS 70, two planet candidates are detected within the inner cavity of the transitional disk around HD 100546. Near-IR spectroscopic monitoring of fundamental rovibrational CO emission lines revealed a spectroastrometric evidence of an orbiting companion within the inner cavity, at a deprojected separation of $\sim 12$ au from the central star (Brittain et al. 2013, 2014, 2019). In addition, a point-source submillimeter continuum is detected with the Atacama Large Millimeter/submillimeter Array (ALMA) at a deprojected separation of 7.8 au (Pérez et al. 2019). The deprojected distances suggest that the two planets may be in/ near 2:1 mean motion resonance (orbital period ratio $\sim 1.9:1$). If confirmed, HD 100546 would offer another convincing example of giant planets in resonance opening inner cavity in transitional disk.

## 4.2. Circumplanetary Disks: Can We Explain the Submillimeter Continuum Emission?

Submillimeter observations with ALMA revealed spatially unresolved continuum emission at the vicinity of PDS 70b and c, possibly originating from dusty CPDs (Isella et al. 2019; Figure 5(a)). Here, we discuss if the observed continuum flux can be explained with thermal emission from the dust in the CPDs.

If the CPDs are in a steady state such that the CPD gas density remains constant over time, we can assume that gas is supplied from the circumstellar disk onto CPDs at a rate comparable to the planets’ accretion rate, $\dot{M}_{acc} \simeq 10^{-8}-10^{-7} M_{Jup} \text{ yr}^{-1}$ (Wagner et al. 2018; Haffert et al. 2019; Thanathibodee et al. 2019). Coupled with the gas flow from the circumstellar disk, only small, sub-µm dust is replenished as can be seen from Figures 1 and 4. We use $f_{dtg}$ for the mass ratio between this small dust to circumstellar disk gas at the outer edge of the gap. In our simulations, we find that this is of the order of $10^{-3}$ (Figure 1). Since sub-µm grains are expected to be well coupled to the CPD gas, they would accrete onto the planets over the CPD gas viscous timescale, which can be written as $\tau_{vis} \simeq T_p/(5\pi\alpha_{CPD})$, where $T_p$ is the planet’s orbital period and $\alpha_{CPD}$ is the viscosity parameter of the CPD (Zhu et al. 2011). We thus expect that the CPDs would have a steady-state dust mass of

$$ 
\begin{aligned}
M_{\text{dust,CPD}} &\simeq \dot{M}_{\text{acc}} \tau_{\text{vis}} \\
&= 2 \times 10^{-5} M_{\oplus} \left( \frac{\dot{M}_{\text{acc}}}{10^{-8} M_{\text{Jup}} \text{ yr}^{-1}} \right) \left( \frac{f_{\text{dtg}}}{10^{-3}} \right) \\
&\qquad \times \left( \frac{\alpha_{\text{CPD}}}{10^{-3}} \right)^{-1} \left( \frac{T_p}{100 \text{ yr}} \right).
\end{aligned}
\qquad (6) $$
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) The observed continuum emission at 855 µm, reproduced from Isella et al. (2019). (b) A simulated continuum image with the PDS 70b-only model, using spatial resolution (93 × 74 mas synthesized beam), rms noise level (18 µJy beam⁻¹), disk inclination (51°7), and position angle (156°7) comparable to the observation. (c) Same as the middle left panel, but with the $M_c = 2.5 M_{Jup}$ model. The CPD dust mass is distributed over a range of dust sizes from 0.1 µm to 1 mm (see the text). (d) Same as the middle right panel, but using a 30 × 30 mas synthesized beam and 14 µJy beam⁻¹ rms noise, requiring 5 hr of on-source integration. The synthesized beam is shown in red at the bottom right corner of each panel.</figcaption></figure>

Assuming semimajor axes of 22 and 35 au for PDS 70b and c and a stellar mass of 0.85 $M_⊕$, the orbital periods of the planets are 112 and 225 yr, respectively. This results in $M_{dust,CPD}$ of $2.2 \times 10^{-5} M_⊕$ and $4.5 \times 10^{-5} M_⊕$.

Adopting the distance $d = 113$ pc (Gaia Collaboration et al. 2016, 2018), the optically thin continuum flux at 855 µm (350 GHz) is

$$F_{\nu} = \frac{\kappa_{\nu} B_{\nu} (T_{\text{dust,CPD}})}{d^2} M_{\text{dust,CPD}} \\ = 0.15 \ \mu\text{Jy} \left( \frac{\kappa_{\nu}}{2 \text{ cm}^2 \text{ g}^{-1}} \right) \left( \frac{T_{\text{dust,CPD}}}{20 \text{ K}} \right) \left( \frac{M_{\text{dust,CPD}}}{2 \times 10^{-5} M_{\oplus}} \right). \quad (7)$$

With an opacity of ~2 cm² g⁻¹ for sub-µm grains at 855 µm (Figure 8) and a stellar irradiation-dominated CPD temperature of ~30–35 K (Appendix B), Equation (7) implies that the continuum flux is estimated to be <1 µJy, far too small to explain the observed continuum flux of ~100 µJy (Isella et al. 2019). We may invoke in situ grain growth within CPDs; however, although a larger submillimeter grain opacity of ~10 cm² g⁻¹ helps, the estimated continuum flux is still a few µJy, more than an order of magnitude smaller than the observed flux.

To explain the observed continuum flux we thus need either a small CPD viscosity of $\alpha_{CPD} \lesssim 10^{-5}$ or accumulation of dust in CPDs over multiple viscous timescales so that the total dust mass is much larger than the steady-state dust mass estimated. One possibility for the latter is that sub-µm grains grow in situ in the CPDs to submillimeter sizes and the submillimeter grains remain trapped in pressure bumps, created perhaps by existing moons and/or radially varying mass transport throughout the disk. Since submillimeter grains in CPDs have short radial drift timescales compared with the gas viscous timescale (e.g., Zhu et al. 2018), pressure bumps are necessary if we were to explain the continuum flux with submillimeter grains. Interestingly, this overall picture is very similar to what we infer to commonly happen in protoplanetary disks, where (sub-)µm grains are supplied from the interstellar environment, grow to millimeter sizes (and

beyond), and are trapped in pressure bumps. Also, this picture is similar to some Galilean satellite formation models in the proto-Jovian disk (e.g., Canup & Ward 2002).

We note that the above order-of-magnitude estimates are consistent with results from more detailed calculations of CPD’s submillimeter continuum flux, for instance, that in Zhu et al. (2018). Based on Figure 3 of Zhu et al. (2018) we can infer that, with $M_p \dot{M}_{acc} \sim 10^{-7} M_{Jup} \text{ yr}^{-1}$ and $\alpha_{CPD} = 10^{-3}$, the CPD continuum flux can reach 100 µJy only when the dust-to-gas mass ratio is $f_{dtg} \gtrsim 0.01$. If dust is depleted by a significant fraction at the gap edge and does not accumulate in CPDs, CPDs with $\alpha_{CPD} \sim 10^{-3}$ will be much fainter than 100 µJy.

In our numerical simulations, test particles captured in CPDs are not accreted onto the planets. Particles hence accumulate in the CPDs as if there are pressure bumps. In the $M_c = 2.5 M_{Jup}$ model, the total dust mass in the CPDs at $t = 0.6$ Myr reaches (2–3) × $10^{-3} M_⊕$. Consistent with the estimates from Equation (7), that is 20–30 µJy assuming sub-µm grain opacity of $\kappa_{\nu} = 2 \text{ cm}^2 \text{ g}^{-1}$, we find that the CPDs are too faint and are not detected at the rms noise level comparable to the observation.<fnref n="15" /> Note, however, that the observed submillimeter flux could be reproduced if sub-µm grains in the CPDs are warmer ($\gtrsim 100$ K) than the stellar irradiation-dominated temperature. If CPD’s internal heating and/or planet’s accretion heat play an important role, the observed submillimeter CPD flux could be explained with sub-µm grains only.

As an alternative, one may invoke in situ grain growth within the CPDs. In order to examine such a possibility, we redistribute the total CPD dust mass over a range of dust sizes from 0.1 µm to 1 mm, assuming a power-law size distribution with an exponent of -3.5. With this implementation, most of the CPD dust mass is in submillimeter grains and the observed CPD continuum flux can be reproduced as can be seen in Figure 5(c), thanks to a larger grain opacity. Again, note that

<footnote n="15">¹⁵ We find that reducing the rms noise level to 5 µJy beam⁻¹ would reveal the CPDs at 3σ level.</footnote>
the CPDs need to have pressure bumps to trap those large
grains because they otherwise are subject to rapid radial drift.

Future higher angular resolution observations will be able to
separate PDS 70c's CPD from the outer continuum ring
(Figure 5(d)). Because the upper layers of the near (i.e.,
western) side of the outer circumstellar disk will block optical/
IR emission from PDS 70c as it moves behind the upper layers,
high spatial resolution observations at (sub)millimeter wave-
lengths will be crucial if we were to constrain PDS 70c's orbit.

4.3. *Potential Implications to Planet Formation in the Solar System and Mature Exoplanetary Systems*

In the solar nebula, an increasing number of meteoritic isotope
measurements suggests that the solar protoplanetary disk had two
genetically distinct reservoirs, which coexisted and remained
spatially separated (e.g., Warren 2011; Kruijer et al. 2017). One
possible way to explain this so-called noncarbonaceous/carbonac-
eous meteorite dichotomy is if Jupiter (and possibly Saturn, too)
has grown early to open a gap around its orbit within ~1 Myr,
preventing large particles from flowing into the inner disk (Kruijer
et al. 2017). The submillimeter continuum ring located beyond
the two forming planets in the PDS 70 disk provides strong
observational evidence that such a filtration mechanism also could
have operated in the solar protoplanetary disk.

The *Kepler* mission revealed that more than 30% of nearby
Sun-like, FGK-type stars have at least one close-in super-Earth/
mini-Neptune (Petigura et al. 2013; Zhu et al. 2018), whereas the
solar system does not have such a planet. Similar to what we see in
the PDS 70 disk and our simulations, it is likely that Jupiter and
Saturn reduced the inward solid mass flux in the solar nebula by
trapping large grains beyond their orbits. As a result, a significantly
less amount of solid would have been available in the inner disk, in
which case the growth of terrestrial planets had to be limited (see
also, e.g., Haugbølle et al. 2019; Lambrechts et al. 2019). Future
statistical comparisons of the warm/cold Jupiter occurrence rate
between Earth hosting stars and super-Earth/mini-Neptune hosting
stars may help reveal if giant planets indeed play a role in
determining the final mass of terrestrial planets in the system.

Our simulations suggest that giant planets could settle into
mean motion resonance early while they grow embedded in their
gaseous host disk. In the solar system, it is proposed that Jupiter
and Saturn have captured in mean motion resonance while they
are embedded in the solar nebula (Masset & Snellgrove 2001;
Morbidelli & Crida 2007; Walsh et al. 2011). An early settlement
into mean motion resonance is also consistent with the directly
imaged four giant planets orbiting in the HR 8799 debris disk
that are suggested to be in 8:4:2:1 mean motion resonance
(Konopacky et al. 2016; Wang et al. 2018).

5. Conclusion

PDS 70 offers an ideal testbed for planet–disk interaction
studies. Using two-dimensional hydrodynamic planet–disk
interaction simulations, we show that the signposts of planet–
disk interaction predicted by numerical simulations show an
excellent agreement with observed features in the PDS 70 disk.
This strongly suggests that previously proposed theories of
planet–disk interaction, including resonant migration, particle
trapping, size segregation, and filtration, are indeed in action. In
particular, the submillimeter continuum ring observed outward
of the two directly imaged planets provide the first

observational evidence that gap-opening planets can hold
particles with appropriate sizes beyond their orbits.

By studying planets in formation and the coevolution with their
host disk, we can also infer the formation history of mature
planetary systems. The fact that giant planets filter large grains
suggests gas giants can have an influence on the formation and
characteristics of terrestrial planets in the system, as they can
induce a chemical inhomogeneity and reduce the inward solid
mass flux. It will be interesting to know if it is common for giant
planets to settle in a mean motion resonance while they are
growing embedded in the host disk. Transitional disks having
inner cavities might be good targets for this purpose.

Although we show that two already-grown giant planets placed
close to 2:1 mean motion resonance can quickly settle into an
orbital resonance, our simulations do not address how the two
planets have reached such a configuration in the first place. On one
hand, it is possible that the two planets form and grow near the
commensurability. In this case, the two planets should have grown
at a rate of ≳10⁻⁶ M<sub>Jup</sub> yr⁻¹ on average, potentially punctuated by
even higher accretion rates during episodic accretion events.
Whether or not the planets can remain dynamically stable with
such mass growth rates have to be examined. On the other hand,
the two planets could have formed and grown at distance, but later
migrated toward each other and captured in a mean motion
resonance. In this case, the migration must be sufficiently slow so
that the planets do not cross the 2:1 mean motion resonance but
also they do not trigger a dynamical instability. Future simulations
taking into account both orbital migration and mass growth from
one to a few tens of Earth-mass cores will help infer the full history
of PDS 70b and c's formation and evolution.

We thank the anonymous referee for prompt and helpful reports.
This Letter makes use of the following ALMA data: ADS/JAO.
ALMA #2015.1.00888.S, ADS/JAO.ALMA #2017.A.00006.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. J.B. acknowledges support by NASA
through the NASA Hubble Fellowship grant #HST-HF2-
51427.001-A awarded by the Space Telescope Science Institute,
which is operated by the Association of Universities for Research
in Astronomy, Incorporated, under NASA contract NAS5-26555.
J.B. acknowledges support from NASA grant NNX17AE31G and
computing resources provided by the NASA High-End Computing
(HEC) Program through the NASA Advanced Supercomputing
(NAS) Division at Ames Research Center. Z.Z. acknowledges
support from the National Aeronautics and Space Administration
through the Astrophysics Theory Program with grant No.
NNX17AK40G and support from the Sloan Research Fellowship.
L.P. acknowledges support from CONICYT project Basal AFB-
170002 and from FONDECYT Iniciación project #11181068.
R.T. acknowledges support from the Smithsonian Institution as a
Submillimeter Array (SMA) Fellow.

Software: Dusty FARGO-ADSG (Baruteau et al. 2019),
RADMC-3D (Dullemond et al. 2012).
<figure><img src="image_9.png" /><figcaption>Figure 6. (a) Blue curves show the density-weighted, vertically integrated temperature distribution $\bar{T}$ from the first 10 MCRT iterations. The blue dashed curve shows the power-law fit to the temperature profile obtained from the 10th iteration: $\bar{T}(R) = 44$ K ($R/22$ au)$^{-0.24}$. The black curve shows the stellar irradiation-dominated temperature $T_{irr}$ (Equation (8)). (b) The two-dimensional R-Z temperature distribution after 10 MCRT iterations.</figcaption></figure>

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

## Appendix A
### Initial Disk Temperature

To set up the initial disk temperature profile, we first construct a stellar irradiation-dominated temperature profile $T_{irr}$:

$$T_{irr}(R) = \left( \frac{fL_{*}}{4\pi R^{2}\sigma_{SB}} \right)^{1/4} \quad (8)$$

Here, $f = 0.1$ is introduced to account for the non-normal irradiation at the disk surface, $L_{*} = 0.35 L_{\odot}$ (Pecaut & Mamajek 2016; Keppler et al. 2018) is the stellar luminosity, and $\sigma_{SB}$ is the Stefan-Boltzmann constant. Assuming that the disk has a constant temperature over height at each radius, we construct the three-dimensional gas density structure that satisfies the hydrostatic equilibrium in the vertical direction:

$$-\frac{GM_{*}Z}{(R^2 + Z^2)^{3/2}} - \frac{1}{\rho} \frac{\partial P}{\partial Z} = 0. \quad (9)$$

Here, $G$ is the gravitational constant, $M_{*}$ is the stellar mass, $Z$ is the height, $\rho$ is the three-dimensional gas density, and $P = \rho RT/\mu$ is the gas pressure with $R$ and $\mu$ being the gas constant and mean molecular weight. We adopt the meridional boundary of the three-dimensional MCRT calculation domain at $Z/R = \pm 0.45$ (i.e., $\pm 24^{\circ}$ against the midplane). At the location of the two planets, this covers more than 6 scale heights with the stellar irradiation-dominated disk temperature.

We fix the gas surface density profile, as described in Equation (1), and iterate MCRT calculations to compute the three-dimensional disk temperature profile using RADMC-3D (Dullemond et al. 2012). After each MCRT calculation, the three-dimensional gas density structure is updated with the new temperature. We stop the iteration when the density-weighted, vertically integrated temperature distribution $\bar{T}$, defined as

$$\bar{T}(R) = \frac{\int \rho(R, Z) T(R, Z) dZ}{\int \rho(R, Z) dZ}, \quad (10)$$

does not vary more than 1% at each radius from the previous iteration. We find that $\bar{T}$ converges within the first few iterations (Figure 6) and the final temperature profile is not sensitive to the choice of the input temperature profile.

In each Monte Carlo calculation we use $10^9$ photon packages. We assume a total of $2 \times 10^{-7} M_{\odot}$ of small dust between 0.01

and 0.1 µm with a power-law size distribution, adopting a power-law exponent of $-3.5$. We assume that small grains are perfectly coupled with disk gas. We assume these small grains are compact monomers consisting of 60% silicate and 40% amorphous carbon, having an internal density of $2.7 \text{ g cm}^{-3}$. We adopt optical constants of silicate and amorphous carbon from Draine & Lee (1984) and Li & Greenberg (1997), respectively.

## Appendix B
### Simulated Continuum Observation

To produce the simulated continuum images presented in Figure 5, we take the gas density distribution at $t = 0.6$ Myr and run MCRT iterations as explained in Appendix A. One additional step we had is that, before we expand the two-dimensional surface density from the hydrodynamic simulation in the vertical direction, we remove the material within the Hill sphere of the planets and fill the region with the azimuthally averaged circumstellar disk surface density at the corresponding radial location. This is necessary because the CPDs have their own scale heights that are much smaller than that of the circumstellar disk. The temperatures at the location of the CPDs from the MCRT iterations are therefore stellar irradiation-dominated temperature. If internal (e.g., viscous) and/or planet's accretion heat play a role, the CPD temperature can be higher than the stellar irradiation-dominated temperature.

The temperature profile from MCRT iterations is shown in Figure 7. The temperature is slightly lower than the initial temperature within the gap and is higher beyond the gap because the outer disk is more directly exposed to the stellar photons.

We assume that grains are composed of 30% silicate matrix and 70% water ice, having an internal density of $1.26 \text{ g cm}^{-3}$. The optical constants of water ices and astrosilicates are adopted from the Jena database<fnref n="16" /> and Draine & Lee (1984), respectively. In the left panel of Figure 8, we show the size-averaged absorption opacity as a function of the observing wavelength, assuming a power-law size distribution with a power-law exponent of $-3.5$ and the minimum and maximum grain sizes of 0.1 µm and 1 mm. Presented in the right panel is size-dependent absorption and scattering opacities at $\lambda = 855$ µm. Simulated continuum images are produced considering both absorption and anisotropic scattering using the Henyey-Greenstein approximation.

<footnote n="16">¹⁶ https://www.astro.uni-jena.de/Laboratory/Database/databases.html</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Similar to Figure 6, but with the gas density distribution at $t = 0.6$ Myr in the $M_c = 2.5 M_{\text{Jup}}$ model. In panel (a), the shaded region shows the radial location of the common gap opened by PDS 70b and c. The orange curve shows the disk midplane temperature obtained after 10 MCRT iterations. The black curve shows the initial $\bar{T}(R)$ distribution (Equation (2)).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** (Left) The absorption (solid) and scattering (dashed) opacities $\kappa_{\nu}^{\text{abs}}$ and $\kappa_{\nu}^{\text{sca}}$ as a function of the observing wavelength $\lambda$. (Right) The absorption (solid) and scattering (dashed) opacities at $\lambda = 855$ µm, as a function of dust size.</figcaption></figure>

## ORCID iDs

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Zhaohuan Zhu https://orcid.org/0000-0003-3616-6822

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

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

Cornelis P. Dullemond https://orcid.org/0000-0002-7078-5910

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

Andrea Isella https://orcid.org/0000-0001-8061-2207

Miriam Keppler https://orcid.org/0000-0001-7250-074X

Laura M. Pérez https://orcid.org/0000-0002-1199-9564

Richard Teague https://orcid.org/0000-0003-1534-5186

## References

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

* Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

* Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657

* 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., & Papaloizou, J. C. B. 2013, ApJ, 778, 7

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

* Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

* Brittain, S. D., Carr, J. S., Najita, J. R., Quanz, S. P., & Meyer, M. R. 2014, ApJ, 791, 136

* Brittain, S. D., Najita, J. R., Carr, J. S., et al. 2013, ApJ, 767, 159

* Brittain, S. D., Najita, J. R., & Carr, J. S. 2019, ApJ, 883, 37

* Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404

* Charnoz, S., Fouchet, L., Aleon, J., & Moreira, M. 2011, ApJ, 737, 33

* Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJL, 877, L33

* de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529

* Dong, R., Hashimoto, J., Rafikov, R., et al. 2012, ApJ, 760, 111

* Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

* 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, Astrophysics Source Code Library, ascl:1202.015

* Dürmann, C., & Kley, W. 2017, A&A, 598, A80

* Fouchet, L., Maddison, S. T., Gonzalez, J.-F., & Murray, J. R. 2007, A&A, 474, 1037

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

* Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

* Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, NatAs, 3, 749

* Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJL, 758, L19

* Hashimoto, J., Tsukagoshi, T., Brown, J. M., et al. 2015, ApJ, 799, 43

* Haugbølle, T., Weber, P., Wielandt, D. P., et al. 2019, AJ, 158, 55

* Isella, A., Benisty, M., Teague, R., et al. 2019, ApJL, 879, L25

* Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

* Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118

* Kley, W. 1999, MNRAS, 303, 696

* Konopacky, Q. M., Marois, C., Macintosh, B. A., et al. 2016, AJ, 152, 28

* Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, PNAS, 114, 6712

* Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83

* Li, A., & Greenberg, J. M. 1997, A&A, 323, 566

* Long, Z. C., Akiyama, E., Sitko, M., et al. 2018, ApJ, 858, 112

* Lubow, S. H. 1991a, ApJ, 381, 259
Lubow, S. H. 1991b, ApJ, 381, 268

Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001

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

Masset, F., & Snellgrove, M. 2001, MNRAS, 320, L55

Morbidelli, A., & Crida, A. 2007, Icar, 191, 158

Muley, D., Fung, J., & van der Marel, N. 2019, ApJL, 879, L2

Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2

Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123

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

Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263

Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794

Pérez, S., Casassus, S., Hales, A., et al. 2019, arXiv:1906.06305

Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, PNAS, 110, 19273

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

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJL, 860, L13

Quillen, A. C., & Trilling, D. E. 1998, ApJ, 508, 707

Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJL, 860, L12

Thanathibodee, T., Calvet, N., Bae, J., Muzerolle, J., & Hernández, R. F. 2019, ApJ, in press (arXiv:1909.06450)

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58

Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJL, 863, L8

Walsh, K. J., Morbidelli, A., Raymond, S. N., O'Brien, D. P., & Mandell, A. M. 2011, Natur, 475, 206

Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192

Ward, W. R., & Canup, R. M. 2010, AJ, 140, 1168

Warren, P. H. 2011, E&PSL, 311, 93

Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101

Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850

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