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

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

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{gas},\mathrm{init}}(R)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{-1}\exp \left(-\displaystyle \frac{R}{{R}_{c}}\right),\end{eqnarray} \tag{ 1 }$$

where we choose R_c = 40 au and Σ_c = 2.7 g cm⁻² so that the total disk gas mass within the simulation domain is 0.003 M_⊙, 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

$$\begin{eqnarray}&&\bar{T}(R)=44\,{\rm{K}}{\left(\displaystyle \frac{R}{22\mathrm{au}}\right)}^{-0.24}.\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&H/R=0.067{\left(\displaystyle \frac{R}{22\mathrm{au}}\right)}^{0.38}.\end{eqnarray} \tag{ 3 }$$

We fix PDS 70b's mass to 5 M_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_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⁴ yr, while fixing their orbits. After 10⁴ 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_accΩΔ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 τ_acc = (f_accΩ)⁻¹. We choose f_acc = 0.01, with which the planets accrete at 10⁻⁸–10⁻⁷ M_Jup yr⁻¹, 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).

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π 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 (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 Σ_gas,damp over the local viscous timescale following

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{gas},\mathrm{damp}}={{\rm{\Sigma }}}_{\mathrm{gas},\mathrm{init}}\exp (-t/{t}_{\mathrm{vis}}).\end{eqnarray} \tag{ 4 }$$

Here, t_vis = R²/ν is the viscous timescale and $\nu =\alpha {c}_{s}^{2}/{\rm{\Omega }}$, where α is the viscosity parameter, c_s is the sound speed, and Ω is the orbital frequency. A uniform disk viscosity α = 10⁻³ 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_⊕. We assume a dust bulk density of 1.26 g cm⁻³, 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 μ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 μ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

$$\begin{eqnarray}&&\mathrm{St}\simeq 0.07\left(\displaystyle \frac{s}{1\,\mathrm{mm}}\right)\left(\displaystyle \frac{R}{40\,\mathrm{au}}\right)\exp \left(\displaystyle \frac{R}{40\,\mathrm{au}}\right).\end{eqnarray} \tag{ 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 α = 10⁻³. 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 μm to 1 mm, to have the mass per interval in $\mathrm{log}(s)$ be proportional to s^0.5 (this corresponds to a dust size distribution of n(s) ∝ s^−3.5).

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 ∼ 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 ± 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 ≃ 32 au. As 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 significantly reduced. On the other hand, wave modes from the eccentric component of the planet's potential launch at R_LR,ecc = [(m ± 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).

Zoom In Zoom Out Reset image size

Muley et al. (2019) recently showed that accreting planets can have an abrupt increase in its orbital eccentricity to e ∼ 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 ∼2.5 M_Jup/3.5 Myr ≃ 7 × 10⁻⁷ M_Jup yr⁻¹, 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 ∼ 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.

We plot orbital elements of the two planets in Figures 2 and 3. As can be seen from the resonant angle, defined as ψ_cb = 2λ_b − λ_c − ϖ_c, where λ_b and λ_c are mean longitudes of the planets, and ϖ_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).

Zoom In Zoom Out Reset image size

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⁻¹, 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⁻¹. 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⁻¹, respectively.

Zoom In Zoom Out Reset image size

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 ∼ 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 ∼ 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α observations (10⁻⁸ − 10⁻⁷ M_Jup yr⁻¹; Wagner et al. 2018; Haffert et al. 2019; Thanathibodee et al. 2019) is unlikely to destabilize the mean motion resonance.

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/α 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 ≳ α are efficiently filtered at the gap edge, whereas small particles with St ≲ α 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 ≲10 μm have a broad radial distribution, whereas large particles with sizes ≳10 μ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-μm grains¹⁴ are well coupled with gas so leak gradually over time, flowing into the inner disk. It is also interesting to note that some sub-μ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.

Zoom In Zoom Out Reset image size

Figures 5(a)–(c) show the observed continuum image at 855 μ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 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.

Zoom In Zoom Out Reset image size

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 ∼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 ∼1.9:1). If confirmed, HD 100546 would offer another convincing example of giant planets in resonance opening inner cavity in transitional disk.

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}}_{\mathrm{acc}}$ ≃ 10⁻⁸–10⁻⁷ M_Jup yr⁻¹ (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⁻³ (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 τ_vis ≃ T_p/(5πα_CPD), where T_p is the planet's orbital period and α_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{eqnarray}\begin{array}{rcl}{M}_{\mathrm{dust},\mathrm{CPD}} & \simeq & {\dot{M}}_{\mathrm{acc}}{\tau }_{\mathrm{vis}}\\ & = & 2\times {10}^{-5}{M}_{\oplus }\left(\displaystyle \frac{{\dot{M}}_{\mathrm{acc}}}{{10}^{-8}\,{M}_{\mathrm{Jup}}\,{\mathrm{yr}}^{-1}}\right)\left(\displaystyle \frac{{f}_{\mathrm{dtg}}}{{10}^{-3}}\right)\\ & & \times {\left(\displaystyle \frac{{\alpha }_{\mathrm{CPD}}}{{10}^{-3}}\right)}^{-1}\left(\displaystyle \frac{{T}_{p}}{100\,\mathrm{yr}}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

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 × 10⁻⁵ M_⊕ and 4.5 × 10⁻⁵ M_⊕.

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

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\nu } & = & \displaystyle \frac{{\kappa }_{\nu }{B}_{\nu }({T}_{\mathrm{dust},\mathrm{CPD}})}{{d}^{2}}{M}_{\mathrm{dust},\mathrm{CPD}}\\ & = & 0.15\,\mu \mathrm{Jy}\left(\displaystyle \frac{{\kappa }_{\nu }}{2\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}}\right)\left(\displaystyle \frac{{T}_{\mathrm{dust},\mathrm{CPD}}}{20\,{\rm{K}}}\right)\left(\displaystyle \frac{{M}_{\mathrm{dust},\mathrm{CPD}}}{2\times {10}^{-5}\,{M}_{\oplus }}\right).\end{array}\end{eqnarray} \tag{ 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 α_CPD ≲ 10⁻⁵ 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}}_{\mathrm{acc}}\sim {10}^{-7}$ M_Jup² yr⁻¹ and α_CPD = 10⁻³, the CPD continuum flux can reach 100 μJy only when the dust-to-gas mass ratio is f_dtg ≳ 0.01. If dust is depleted by a significant fraction at the gap edge and does not accumulate in CPDs, CPDs with α_CPD ∼ 10⁻³ 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⁻³ M_⊕. Consistent with the estimates from Equation (7), that is 20–30 μJy assuming sub-μm grain opacity of κ_ν = 2 cm² g⁻¹, we find that the CPDs are too faint and are not detected at the rms noise level comparable to the observation.¹⁵ Note, however, that the observed submillimeter flux could be reproduced if sub-μm grains in the CPDs are warmer (≳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 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 wavelengths will be crucial if we were to constrain PDS 70c's orbit.

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/carbonaceous 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).