Testing planet formation from the ultraviolet to the millimeter

The protoplanetary disc around AS 209 has five annular gaps in dust continuum emission at orbital radii $r\approx 9,\,24,\,35,\,61,$ and $99$ au (Guzmán et al., 2018; Huang et al., 2018). Zhang et al. (2018) demonstrated that the four gaps from $24-99$ au can be produced by a single planet situated within the outermost gap at 99 au (see their figure 19). Favre et al. (2019) and Zhang et al. (2021) observed a broad gap in C¹⁸O (2-1) emission extending from $\sim$60 to 100 au. Favre et al. (2019) modeled the emission using the DALI thermo-chemical code (Bruderer, 2013) and took the gas density profile to include two gaps coincident with the dust gaps at 61 and 99 au. They showed these two gas gaps could be smeared by the ALMA beam into a single wide gap like that observed in molecular emission. At the bottom of the 99 au gap where the planet is believed to reside, they fitted $\Sigma_{\rm gas}\approx 0.04-0.08\,\rm g/cm^{2}$ (see the blue and green curves in their figure 6).

Zhang et al. (2021) modeled the C¹⁸O emission using the thermo-chemical code RAC2D (Du & Bergin, 2014), including only a single gap. They initialized their models with an ISM-like $X_{\rm CO}$ which decreased from CO freeze-out and chemistry over the simulation duration of 1 Myr. They found $\Sigma_{\rm gas}\approx 0.08\,\rm g/cm^{2}$ at 99 au (their figure 16), similar to Favre et al. (2019). Zhang et al. (2021) also considered the possibility that there is no gap in H₂ and that the observed gap in C¹⁸O is due to $X_{\rm CO}$ being somehow locally depleted beyond the predictions of RAC2D. Assuming a total gas-to-dust ratio of $10$ by mass (their table 2), they estimated in this super-depleted $X_{\rm CO}$ scenario $\Sigma_{\rm gas}\approx 0.4\,\rm g/cm^{2}$ at 99 au (their figure 5; see also Alarcón et al. 2021 who advocated on separate grounds for a depleted CO:H₂ ratio in AS 209).

We combine the results of Favre et al. (2019) and Zhang et al. (2021) to consider $\Sigma_{\rm gas}=0.04-0.4\,\rm g/cm^{2}$ for the gap AS 209-G2 at $r=99$ au (Table 1).

Dust continuum observations of HD 163296 reveal gaps centred at 10, 48, and 86 au (Isella et al., 2016; Huang et al., 2018), each opened by a separate planet according to Zhang et al. (2018). Isella et al. (2016) fitted the emission at $r\gtrsim 25$ au from three CO isotopologues assuming an ISM-like value of $X_{\rm CO}$ and found $\Sigma_{\rm gas}\approx 2.5-10\,\rm g/cm^{2}$ at 48 au and $\Sigma_{\rm gas}\approx 0.1-2\,\rm g/cm^{2}$ at 86 au (their figure 2, blue curve; note the gap radii in their figure are about 20% larger than the values we use because they used a pre-GAIA source distance). The thermo-chemical models of Zhang et al. (2021) predict roughly consistent values of $\Sigma_{\rm gas}\approx 1\,\rm g/cm^{2}$ and $\Sigma_{\rm gas}\approx 0.6\,\rm g/cm^{2}$ respectively (see their figure 16). If instead the CO:H₂ abundance is much lower than predicted by their models and the gas-to-dust ratio is $60$ (their table 2), then $\Sigma_{\rm gas}\approx 40\,\rm g/cm^{2}$ and $\Sigma_{\rm gas}\approx\,20\,\rm g/cm^{2}$ (their figure 5).

The above estimates of $\Sigma_{\rm gas}$ are combined to yield the full ranges quoted in Table 1, $\Sigma_{\rm gas}=1-40\,\rm g/cm^{2}$ for HD 163296-G2 at 48 au, and $\Sigma_{\rm gas}=0.1-20\,\rm g/cm^{2}$ for HD 163296-G3 at 86 au.

IM Lup presents a shallow gap at $r\approx 117$ au in millimeter continuum images (Huang et al., 2018) which Zhang et al. (2018) attributed to a 10-30 $M_{\oplus}$ planet. Zhang et al. (2021) found $\Sigma_{\rm gas}\approx 0.1\,\rm g/cm^{2}$ in their RAC2D modeling which includes CO freeze-out and chemistry (their figure 16). Alternatively, for a super-depleted CO:H₂ scenario and gas-to-dust ratio of 100 (their table 2), $\Sigma_{\rm gas}\approx 10\,\rm g/cm^{2}$ (their figure 5). We allow for both possibilities and take $\Sigma_{\rm gas}=0.1-10\,\rm g/cm^{2}$.

HD 169142 features an annular gap at 50 au in near-infrared scattered light (Quanz et al., 2013) and mm-wave continuum emission (Fedele et al., 2017). In C¹⁸O emission, Fedele et al. (2017) found evidence for a cavity inside $r\approx 50$ au and used DALI to infer $\Sigma_{\rm gas}\approx 0.1-0.2\,\rm g/cm^{2}$ near the cavity outer edge (their figure 5, top left panel, blue curve). Midplane gas temperatures according to the same thermo-chemical model (figure 5, bottom center panel) are $\sim$$30-40$ K, higher than the CO freeze-out temperature of $\sim$20 K. CO may still be depleted for other reasons (Schwarz et al. 2018) but without a model we cannot quantify the uncertainty.

Near-infrared scattered light images reveal annular gaps at $r\approx 21$ au (TW Hya-G1) and 85 au (TW Hya-G2) (van Boekel et al., 2017). The inner gap may also have been detected in the mm continuum by Andrews et al. (2016), albeit slightly outside the gap in scattered light (their figure 2). Nomura et al. (2021) studied TW Hya in the $J=3-2$ transition of C¹⁸O. Near the inner gap they measured a wavelength-integrated line flux of 5 mJy beam^-1 km s^-1 (their figure 2, middle panel, orange curve). We extrapolate their radial intensity profile, which cuts off at 80 au, to the outer gap at 85 au and estimate there a flux of 1 mJy beam^-1 km s^-1. Using temperatures of 60 K and 20 K for the inner and outer gaps respectively (section 2.3.3), we convert these fluxes to C¹⁸O column densities (e.g. Mangum & Shirley, 2015, their equation B2). Then assuming $X_{\rm CO}=10^{-4}$ and a ¹⁸O$:^{16}$O gas-phase abundance ratio of $2\times 10^{-3}$ (e.g. Qi et al., 2011), we calculate $\Sigma_{\rm gas}=0.04\,\rm g/cm^{2}$ and $0.008\,\rm g/cm^{2}$ for the two gaps. These surface densities could be underestimated by factors of 70 and 30 (assuming a gas-to-dust ratio of 100), respectively, if CO is somehow super-depleted (Zhang et al. 2019, their figure 8, second panel, blue curve). Thus we have $\Sigma_{\rm gas}\approx 0.04-2.8\,\rm g/cm^{2}$ for TW Hya-G1 and $\Sigma_{\rm gas}\approx 0.008-0.24\,\rm g/cm^{2}$ for TW Hya-G2.

Facchini et al. (2021) studied PDS 70 in the $J=2-1$ transition of $\rm C^{18}O$, but the beam size at this wavelength was too large to resolve the cavity within which the two planets reside. We therefore resort to using the C¹⁸O emission just outside the cavity to place an upper limit on the gas surface density near the planets. Following the same procedure as for TW Hya, we convert the observed C¹⁸O flux outside the cavity into an H₂ surface density. We use 30 K as the gas excitation temperature at 80 au, close to the H₂CO excitation temperature measured by Facchini et al. (2021, their figure 7). The gas surface density so derived is $\Sigma_{\rm gas}=8\times 10^{-2}\,\rm g/cm^{2}$. We pair this upper limit with a lower limit derived from optically thick ¹²CO in the $J=3-2$ transition. The smaller beam associated with the higher frequency of this transition allowed the cavity to be resolved by Facchini et al. (2021). Near the positions of the planets at $r\sim 22-34$ au, the observed ¹²CO line flux is 37 K km s^-1, which for an assumed temperature of $T\sim 60$ K (scaled from 30 K at 80 au using $T\propto r^{-3/7}$ for a passively irradiated disk; Chiang & Goldreich 1997) yields a lower bound on $\Sigma_{\rm gas}$ of $8\times 10^{-4}\,\rm g/cm^{2}$. The temperatures quoted above are high enough to avoid CO freeze-out, though we cannot rule out that gas-phase CO is depleted from other effects (Schwarz et al., 2018).

For all nineteen of the DSHARP planets, we estimate $\Sigma_{\rm dust}$ at gap center using equation (1) and the optical depth profiles in figure 6 of Huang et al. (2018).

Three of the DSHARP gaps, AS 209-G2 ($r=99$ au), Elias 24 ($57$ au), and HD 143006-G1 ($22$ au), have bottoms that fall below the 2$\sigma$ noise floor at $\tau=10^{-2}$, although apparently only marginally so. For Elias 24 and HD 143006-G1, we set $\tau=10^{-2}$. For AS 209-G2, while portions of the gap lie below the noise floor, at gap center there is actually a local maximum (ring) of emission that rises above the noise. We use the peak optical depth of this ring, $\tau\approx 2\times 10^{-2}$, to evaluate $\Sigma_{\rm dust}$. The “W"-shaped gap profile for AS 209-G2 resembles transient structures found in two-fluid simulations of low-mass planets embedded in low-viscosity discs (Dong et al., 2017; Zhang et al., 2018).

In their re-analysis at higher spatial resolution of the DSHARP data, Jennings et al. (2021, their figure 12) found evidence for deeper dust gaps than originally measured by Huang et al. (2018) for GW Lup, HD 163296-G2, SR 4, AS 209-G1, and HD 143006-G1. The re-analyzed gaps in the mm continuum are deeper by at least a factor of 10 and in most cases much more. We make a note of these sources in Table 1, but continue to use the Huang et al. (2018) results to compute $\Sigma_{\rm dust}$ for consistency with the rest of our analysis.

For TW Hya-G1 at 21 au we calculate $\Sigma_{\rm dust}=0.06-0.4\,\rm g/cm^{2}$ using equations (1)-(2) with $I_{\nu}=10$ mJy beam^-1 (figure 2, middle panel, purple curve of Nomura et al. 2021) and $T=60$ K (section 2.3.3). For TW Hya-G2 at 85 au we estimate $\Sigma_{\rm dust}=0.002-0.02\,\rm g/cm^{2}$ using $I_{\nu}\sim 10^{-1}\,$ mJy beam^-1 (this is an extrapolation of the intensity profile measured by Nomura et al. 2021 which does not extend beyond $\sim$70 au) and $T=20$ K. Our estimates of $\Sigma_{\rm dust}$ are consistent with those plotted in figure 3 of Nomura et al. (2021).

For the gap in HD 169142 and central clearing in PDS 70, the continuum emission falls below the 3$\sigma$ noise floor of the observations (0.21 mJy beam^-1 and 0.0264 mJy beam^-1, respectively; Fedele et al. 2017; Benisty et al. 2021). For assumed temperatures of $T=40$ K in HD 169142 and $T=60$ K in PDS 70, the corresponding upper limits on the dust column are $\Sigma_{\rm dust}\leq 10^{-2}\,\rm g/cm^{2}$ and $\Sigma_{\rm dust}\leq 2\times 10^{-2}\,\rm g/cm^{2}$, respectively.

For most of the DSHARP sample, we adopt the $h/r$ estimates listed in table 3 of Zhang et al. (2018), computed assuming a passive disc irradiated by its host star. For AS 209-G2, all three gaps in HD 163296, and IM Lup, we use more sophisticated estimates of $h/r$ from table 2 of Zhang et al. (2021). In AS 209, Zhang et al. (2021) used the observed ¹³CO emission surface to fit for $h/r=0.06(r/100\,\rm au)^{0.25}$ (since their CO isotopologue data do not extend down to AS 209-G1 at 9 au, we only apply this scale height relation to AS 209-G2 at 99 au). Note that Zhang et al. (2018) found a similar value, $h/r=0.05$ at 99 au, reproduced the multi-gap structure of AS 209. For HD 163296 and IM Lup, Zhang et al. (2021) obtained excellent fits to the spectral energy distributions (their figure 4) using $h/r=0.084(r/100\,\rm au)^{0.08}$ and $h/r=0.1(r/100\,\rm au)^{0.17}$, respectively.³³3 The values of $h$ obtained by Zhang et al. (2021), if interpreted as a standard isothermal hydrostatic disc gas scale height, imply gas temperatures that differ from the gas temperatures obtained from thermo-chemical and radiative transfer modeling by these same authors (their figure 7). This is because the latter temperatures are solved assuming a fixed density field (specified by $h$), and no provision is made to iterate the density field to ensure simultaneous hydrostatic and radiative equilibrium. There is no rigorous justification for using either set of temperatures to compute the gas scale height; we use $h$ as tabulated by Zhang et al. (2021) for convenience (their parameters $H_{100}$ and $\psi$ in table 2, inserted into their equation 3 with $\chi=1$).

Fedele et al. (2017) simultaneously fitted the surface density and temperature profiles in HD 169142. We use their aspect ratio of $h/r=0.07$ at the gap position; this is similar to the value of $h/r=0.08$ found by Dong & Fung (2017) in fitting the gap width to their planet-disk simulations. The aspect ratio relates to the disc midplane temperature following $h/r=c_{s}/(\Omega r)$ where $c_{s}=\sqrt{kT/\mu m_{\rm p}}$ is the sound speed, $\mu=2.4$ is the mean molecular weight, $k$ is Boltzmann’s constant, $m_{\rm p}$ is the proton mass, $\Omega=\sqrt{GM_{\star}/r^{3}}$ is the Keplerian orbital frequency, and $G$ is the gravitational constant. The implied temperature at the gap position is 40 K.

Calahan et al. (2021) obtained a good fit to the spectral energy of TW Hya using $h/r=0.11\,(r/400\,\rm au)^{0.1}$ (their figure 5). From this scaling we calculate $h/r=0.08$ and $h/r=0.09$ at the inner ($r=21$ au) and outer ($r=85$ au) gap positions, respectively. These aspect ratios correspond to temperatures of 58 K and 20 K.

Absent direct estimates of the temperature in the PDS 70 cavity, we estimate the disc scale height using arguments similar to those in section 2.1.6. Facchini et al. (2021) derived an H₂CO excitation temperature near the disc midplane of 33 K at 80 au (see their figure 7). Setting this equal to the gas kinetic temperature there, we calculate $h=c_{s}/\Omega$ at 80 au, and then scale to the location of the planets assuming $h/r\propto r^{2/7}$ (Chiang & Goldreich, 1997). We find $h/r=0.07$ and 0.08 at the locations of PDS 70b and c. Keppler et al. (2018) found a similar temperature profile in their radiative transfer models of the PDS 70 disc.

Zhang et al. (2018) used their grid of planet-disk hydrodynamical simulations to estimate the planet masses that can reproduce the radial widths and in some cases the radial positions of DSHARP gaps. Their table 3 gives $M_{\rm p}$ for different dust grain size distributions and the Shakura & Sunyaev (1973) viscosity parameter $\alpha$. Here we consider $\alpha=10^{-5}-10^{-3}$. Values of $\alpha\gtrsim 10^{-2}$ appear incompatible with various mm-wave observations (Pinte et al., 2016; Teague et al., 2016; Flaherty et al., 2017) including the multi-gap structure of some discs (Dong et al., 2017; Zhang et al., 2018). Planet masses inferred by Zhang et al. (2018) increase with viscosity as $M_{\rm p}\propto\alpha^{1/3}$; more massive planets are required to counteract faster diffusion of gas to open the same width gap. Thus to obtain an upper limit on $M_{\rm p}$ we use the entries corresponding to $\alpha=10^{-3}$ in table 3 of Zhang et al. (2018), and for a lower limit we take their $M_{\rm p}$ values for $\alpha=10^{-4}$ and divide by two in a crude extrapolation down to $\alpha=10^{-5}$. This procedure gives a range on $M_{\rm p}$ that typically spans a factor of 8, and that encompasses the variation in $M_{\rm p}$ due to different dust size distributions.

For HD 163296-G3, we increase the upper bound on $M_{\rm p}$ to $1M_{\rm J}$ based on the fit to non-Keplerian kinematical data from Teague et al. (2018) and Teague et al. (2021). Note that their mass estimate is based on planet-disc simulations assuming $\alpha=10^{-3}$, and that lowering $\alpha$ lowers the planet mass required to reproduce the same non-Keplerian velocity signature ($M_{\rm p}\propto\alpha^{0.41}$ according to equation 16 of Zhang et al. 2018). For Elias 24, the $M_{\rm p}$ range in Table 1 overlaps with mass estimates inferred from infrared photometry of an unconfirmed point source at $r\approx 55$ au (Jorquera et al. 2021). The candidate point source in HD 163296 is located at $r\approx 67$ au (Guidi et al. 2018) and as such cannot be identified with the planets inferred at $r\approx 48$ and 86 au (Zhang et al. 2018) of interest here.

For most of the DSHARP systems we use stellar masses derived from stellar evolution tracks and listed in table 1 of Andrews et al. (2018). For IM Lup, AS 209, and HD 163296 we use instead $M_{\star}=1.2,\,1.1,$ and $2.0\,M_{\odot}$ respectively, calculated from gas disc rotation curves by Teague et al. (2021) as part of the MAPS survey (Oberg et al., 2021).

By comparing the observed gaps in HD 169142 and TW Hya with their grid of planet-disk simulations, Dong & Fung (2017) estimated planet masses of $M_{\rm p}=0.2-2\,M_{\rm J}$ in HD 169142, $0.05-0.5\,M_{\rm J}$ in TW Hya-G1, and $0.03-0.3\,\,M_{\rm J}$ in TW Hya-G2. These values were derived for viscosities $\alpha=10^{-4}-10^{-2}$. For consistency with our adopted range of $\alpha=10^{-5}-10^{-3}$ (see section 2.4.1), we apply the same $M_{\rm p}\propto\alpha^{1/3}$ scaling from Zhang et al. (2018) and shift all limits from Dong & Fung (2017) down by a factor of two.

The stellar masses of HD 169142 and TW Hya are $M_{\star}=1.65\,M_{\odot}$ and $M_{\star}=0.8\,M_{\odot}$, respectively (Blondel & Djie, 2006; Ducourant et al., 2014).

Masses of PDS 70b and c derived from near-infrared spectral energy distributions and cooling models vary widely (Müller et al., 2018; Keppler et al., 2018; Christiaens et al., 2019; Mesa et al., 2019; Wang et al., 2020; Stolker et al., 2020). For simplicity we take $M_{\rm p}=1-10\,M_{\rm J}$ for both planets, which brackets most literature estimates relying on cooling models.⁴⁴4In later sections, and in particular the discussion surrounding Figures 5 and 6, we will consider the possibility that the observed near-infrared emission does not arise from a planet passively cooling into empty space (as assumed by the cooling models) but one which sits below a hot accretion shock boundary layer. In this scenario our adoption of planetary masses for PDS 70b and c drawn from the passive cooling models will lead to an inconsistency. Our hope is that our estimated mass range of $1-10\,M_{\rm J}$ is not seriously impacted; it should not be as long as the posited accretion luminosities are not much larger than the observed near-infrared luminosities ($10^{-4}L_{\odot}$; Wang et al. 2020). It may also help to have the accretion columns be localized to “hot spots” on the planet, like they are for magnetized T Tauri stars (e.g. Gregory et al. 2006), leaving the rest of the planet passively cooling. This range is consistent with masses for PDS 70b derived by Wang et al. (2021) based on dynamical stability constraints. They calculated a 95% upper limit of 10 $M_{\rm J}$, with roughly equal probability per log mass from 1-10 $M_{\rm J}$ (see their figure 5).

The mass of the star is $M_{\star}=0.88\,M_{\odot}$ (Keppler et al., 2019).