Constraining the Nature of the PDS 70 Protoplanets with VLTI/GRAVITY^∗

The observations were carried out using the GRAVITY instrument (Gravity Collaboration et al. 2017) on the VLTI using the four Unit Telescopes (UT). The log of the observations are given in Table 1. Atmospheric conditions ranged from very good (atmospheric coherence time τ₀ = 20 ms) to average (τ₀ ≈ 2 ms). The first observation, in 2018, is a classical interferometric observation of the star (PI M. Benisty, ID 0101.C-0281). The 2019 observations were obtained as a backup of the AGN large program (PI E. Sturm, ID 1103.B-0626, for PDS 70 b), and from director discretionary time (PI A. Vigan, ID 2103.C-5018, for PDS 70 c). Last, the 2020 observations were obtained as part of the ExoGRAVITY large program (PI S. Lacour, ID 1104.C-0651).

The 2018 observations were carried out using the single-field on-axis mode. On-axis means the beam splitter was used (Pfuhl et al. 2014), therefore sending 50% of the flux to the fringe tracker (Lacour et al. 2019) and 50% on the science channel. Single-field means the fringe tracker and the science fibers observed the same object: in this case, the star PDS 70 A. The observations were followed by observations of the calibrator HD 124058. The reduction of this data set was standard using the ESO GRAVITY pipeline ²⁷ (Lapeyrere et al. 2014).

The observations of the exoplanets used the dual-field on-axis mode. Dual-field means the fringe tracker and the science fibers observed different objects. Thanks to the splitter, the fringe tracker observed the star for phase referencing, and the science fibers observe the planet. Non-common path phase aberrations were calibrated by interleaving the observation of the planet with single-field on-axis observations.

The coherent flux was extracted following a standard procedure with the ESO GRAVITY pipeline. From this first step, we obtain V_onplanet(b, t, λ) and V_onstar(b, t, λ), the coherent flux observed on the star and the planet as a function of baseline b, time t, and wavelength λ.

The removal of stellar contamination was performed during a second step. The computation is described in detail in Appendix A of Gravity Collaboration et al. (2020). The code is available as a Python library developed by our team. ²⁸ The main objective of the algorithm is to calculate R(λ, b, t), the ratio of the uncontaminated coherent flux between the star and planet:

$$\begin{eqnarray}&&R(b,t,\lambda ,)=\displaystyle \frac{{\gamma }_{\mathrm{planet}}}{{\gamma }_{\mathrm{star}}}\ \displaystyle \frac{{V}_{\mathrm{planet}}(b,t,\lambda )}{{V}_{\mathrm{star}}(b,t,\lambda )}\end{eqnarray} \tag{ 1 }$$

where V_star(b, t, λ) and V_planet(b, t, λ) are the coherent flux of both objects in the absence of stellar speckle. The additional term γ comes from the fact that the science fiber is not exactly positioned at the location of the target in the focal plane (see Appendix A). The astrometry was obtained from the argument of R, the ratio of the coherent flux:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{R}(b,t,\lambda )=-\displaystyle \frac{2\pi }{\lambda }({\rm{\Delta }}{\rm{R}}.{\rm{A}}.u+{\rm{\Delta }}\mathrm{Decl}.v)\end{eqnarray} \tag{ 2 }$$

in which (u,v) are the coordinates in the frequency domain and (ΔR. A.,ΔDecl. ) are the sky-coordinates of the planet relative to the star.

The values, obtained by χ² minimization, are given in Table 2. The error bars given here correspond to the precision of the measurement. They are estimated from the scatter of the astrometric values (between each file, or by splitting the files into independent measurements). The typical precision is 100 μas. The systematic errors, which are ultimately limiting the accuracy of the astrometry, are theoretically smaller. They were estimated to be 16.5 μ (Lacour et al. 2014).

The coherent fluxes V_star(b, t, λ) and V_planet(b, t, λ) are not normalized and therefore include the shape of the spectrum. To extract the spectrum F_planet(λ), we assumed the planet to be unresolved. The amplitude of coherent flux of the planet is then equal to the planetary flux:

$$\begin{eqnarray}&&{V}_{\mathrm{planet}}(b,t,\lambda )={F}_{\mathrm{planet}}(\lambda )\exp (i{\rm{\Phi }}(b,t,\lambda )),\end{eqnarray} \tag{ 3 }$$

with the phase derived from the astrometry obtained as described in the previous section.

The star itself has an angular size much smaller than 0.1 mas and therefore was not resolved by our observations. However, the visibilities are still below 1 because the interferometer did partially resolve the hot inner disk:

$$\begin{eqnarray}&&{V}_{\mathrm{star}}(b,t,\lambda )={F}_{\mathrm{star}}(\lambda ){J}_{\mathrm{star}}(b,t,\lambda )\end{eqnarray} \tag{ 4 }$$

where J_star(b, t, λ) is the visibility drop due to resolving the central source (star and inner disk). Hence:

$$\begin{eqnarray}&&{F}_{\mathrm{planet}}(\lambda )=\displaystyle \frac{{\gamma }_{\mathrm{star}}{F}_{\mathrm{star}}(\lambda )}{{\gamma }_{\mathrm{planet}}}\unicode{x02039}\displaystyle \frac{R(\lambda ,b,t){J}_{\mathrm{star}}(b,t,\lambda )}{\exp (i{\rm{\Phi }}(b,t,\lambda ))}{\unicode{x0203A}}_{b,t}.\end{eqnarray} \tag{ 5 }$$

The notation corresponds to the mean notation: the coherent flux ratio is averaged over b and t. The injection efficiencies γ_star and γ_planet are real numbers between 0 and 1, where 1 is the theoretical maximum. The value of γ_planet was nearly 1 for most epochs except for the first epoch of PDS 70 c when we did not know the precise position of the planet. During this observation, the fiber was pointing 16.5 mas away from the planet, which gave an injection efficiency of ρ = 0.84. The calculation of this term is given in Appendix A.

To obtain J_star(b, t, λ), we used the observation of PDS 70 A, calibrated by the star HD 124058 (assuming a diameter for the calibrator of 0.136 ± 0.002 mas). The obtained visibilities are shown in Figure 1. They are mostly consistent with a single constant:

$$\begin{eqnarray}&&{J}_{\mathrm{star}}(b,t,\lambda )=\mathrm{cst}.\end{eqnarray} \tag{ 6 }$$

The value of the constant does depend on the normalization of the stellar flux. We found that 94% ± 1% of the flux on-axis of the star is unresolved, independent of baseline. That is, 6% of the flux came in excess emission from the hot inner circumstellar disk that is resolved with GRAVITY. An inner circumstellar disk has been predicted from SED analysis (Dong et al. 2012; Hashimoto et al. 2012; Long et al. 2018), detected in scattered light (Keppler et al. 2018; Mesa et al. 2019), and resolved in the millimeter wavelength (Francis & van der Marel 2020).

Zoom In Zoom Out Reset image size

Because the inner disk is resolved, we cannot simply use the 2MASS K-band photometry of the system (Cutri et al. 2003) to normalize the stellar model of the star to calibrate the planetary spectra, since it would include excess emission from the circumstellar disk, which is not part of F_star(λ). Fortunately, the combined scattered light and thermal emission from the inner disk at shorter wavelengths contributes ≲5% of the total flux from the system (Dong et al. 2012), comparable to the 1σ errors on the stellar photometry. Therefore, we used a BT-NextGen stellar atmosphere (Allard et al. 2012) determined from a joint evolutionary-atmospheric model fit to literature optical and near-infrared photometry using the procedure described in Wang et al. (2020). We did not impose a prior on the effective temperature as in Wang et al. (2020), but all other aspects of the fit were the same. From the resulting posterior distributions of the fitted and derived parameters, we measured an age of 8 ± 1 Myr, a mass of 0.88 ± 0.02 M_⊙, an effective temperature of ${4109}_{-30}^{+36}$ K (higher than the spectroscopically derived value of 3972 ± 36 K; Pecaut & Mamajek 2016), and a surface gravity of $\mathrm{log}g=4.23\pm 0.02$ [dex]. The mass estimated from SED fitting is significantly higher than previous literature estimates from SED fitting (0.76 ± 0.02 M_⊙; Müller et al. 2018), but in better agreement with the dynamical mass estimates from circumstellar disk modeling (0.875 ± 0.03 M_⊙; Keppler et al. 2019) and from the orbit fits presented in Section 3. A synthetic stellar spectrum was computed by randomly drawing 200 samples from the Markov Chain Monte Carlo (MCMC) chain, taking the median and standard deviation of the flux at each wavelength as the adopted spectrum and corresponding uncertainties, respectively. We used this spectrum as the spectrum of the star (F_star) to calibrate our planet spectra. The statistical uncertainties in the model stellar spectrum are much lower than the uncertainties on the planetary spectra, but do not include any estimate of systematic errors in the stellar models.

Comparing our stellar model to the 2MASS K-band photometry, we also measured a significant K-band excess of 14% ± 3%, caused by emission from circumstellar material close to the star. This is much greater than the 6% excess emission resolved by GRAVITY. As the stellar variability is dominated by the rotation modulation of the star (Thanathibodee et al. 2019), it is unlikely to be responsible for this disagreement, as the photometric measurements used in the stellar SED fit should average over this variability. Rather, the GRAVITY observations probe spatial scales of 0.2–0.5 au, which are right in the middle of the spatial extent of the inner disk based on models from Dong et al. (2012) and Long et al. (2018). The remaining ∼8% could be emitted closer in to the star where it is not resolved by GRAVITY, or further out at larger spatial scales that are outside the field of view of GRAVITY (∼50 mas). We did not see significant change in the visibilities over the baselines observed, indicating the emission must be significantly closer in or significantly further out. Inner disk models have the inner edge of the disk at 0.05 au so there could be emission that is a factor of ∼3 closer in (Dong et al. 2012; Long et al. 2018). On the other hand, Keppler et al. (2018) detected the inner disk in polarized light with VLT/SPHERE and found the inner disk could be extended out to 20 au, and Francis & van der Marel (2020) found an outer cutoff of 10 au at millimeter wavelengths with ALMA. This also agrees with disk-modeling analysis that found the inner disk to end at 15 au (Long et al. 2018). Any emission outside of 6 au would not have been seen by GRAVITY (would not couple into the single mode fibers). The location of the emission affects our photometric calibration as emission at larger separations do not need to be included in F_star whereas emission unresolved by GRAVITY needs to be included in F_star, which ultimately changes the absolute brightness of the planets. Without further information at the moment, we assumed that half of the remaining 8% excess emission comes from the star. For simplicity, we assumed that it has the same spectrum as the star, so we essentially multiplied our model stellar spectrum by 1.04. A 4% change is already much smaller than the 1σ uncertainties on the planet spectra, so the exact scale factor and spectral shape of the excess dust emission should negligibly affect our results.

Since the star is young and accreting with clear Hα emission (Thanathibodee et al. 2019), line emission from the star could affect the calibration of the planetary spectra. However, the stellar Paβ line has not been detected (Long et al. 2018). Given that the Brγ line is expected to be even weaker, and given that it should be unresolved, the impact of the Brγ line on a single spectral element of our planetary spectra should be negligible given the relatively large error bars of a single spectral channel (∼15% of the total flux). Similarly, Long et al. (2018) did not find appreciable CO line emission in the K-band, which too should be unresolved in our data, and thus have a negligible impact on our spectrum. Thus, we concluded that our omission of emission lines in our model stellar spectrum is a reasonable approximation.

With the coherent flux ratio R(λ, b, t) and this model spectrum of the star F_star(λ) scaled by 1.04, Equation (5) gave us the spectra for both planets. The resulting spectrum as well as the uncertainties are plotted in Figure 2.

Zoom In Zoom Out Reset image size

We also reanalyzed the SPHERE IFS data on PDS 70 c published in Mesa et al. (2019). Specifically, we reanalyzed the data from 2018-02-24, which were obtained in exquisite conditions (0.40'' seeing). The data were acquired with SPHERE (Beuzit et al. 2019) in its IRDIFS-EXT mode where IFS (Claudi et al. 2008) and IRDIS (Dohlen et al. 2008) observe in parallel, with IFS covering the YJH bands and IRDIS in K1 and K2 band (Vigan et al. 2010). The data were collected with the apodized pupil Lyot coronagraph (Carbillet et al. 2011; Guerri et al. 2011) in its N_ALC_YJH_S configuration optimized for the H band. The raw data were preprocessed using the vlt-sphere ²⁹ open-source pipeline (Vigan 2020) to produce calibrated (x, y, λ) data cubes of coronagraphic images and off-axis reference PSFs.

Stellar PSF subtraction and spectral extraction of PDS 70 c was done using pyKLIP version 2.1 (Wang et al. 2015). We used angular differential imaging (ADI; Liu 2004; Marois et al. 2006) and spectral differential imaging (SDI; Sparks & Ford 2002) to build up a model of the stellar PSF. We used any frame where PDS 70 c moved by one pixel due to ADI and SDI to calculate the principal components to model the star. We used the first 10 principal components, as this gave us the best signal-to-noise on the planet. We used the forward modeling framework described in Pueyo (2016) and Greenbaum et al. (2018) to measure the spectrum of PDS 70 c. We injected eight simulated planets at the same separation but at different azimuthal positions as PDS 70 c, measured their spectra in the same way, and used the scatter in their measured spectra to estimate the uncertainties on the spectrum of PDS 70 c.

Due to the fact the planet is adjacent to the edge of the circumstellar disk, there is concern that the spectral extraction is biased by the disk even with forward modeling. To assess this, we injected five simulated planets at similar separations as PDS 70 c but at other azimuthal positions in the image where the simulated planets would be adjacent to the disk edge and computed the average bias in the flux after spectral extraction. We found and corrected for biases that were at most the size of the 1σ uncertainty of each spectral channel. We also verified that the scatter in flux between these five planets is consistent at the 20% level with the uncertainty we estimated for PDS 70 c in the previous paragraph. In both cases, the errors in the spectral measurements are dominated by the low signal-to-noise of the planet, and not by systematics due to the presence of disk signal.

We found that PDS 70 c is only detected in the last 12 spectral channels of the IFS data (i.e., H-band). Unlike Mesa et al. (2019), the Y- and J-band spectra are consistent with non-detections, and the scatter we measured at those wavelengths is due to noise. We note that the PDS 70 c spectrum from Mesa et al. (2019) only significantly deviates from their median spectrum of the circumstellar disk in the H-band (see their Figure 5c), which may indicate their PDS 70 c spectrum at the Y and J bands are contaminated by the disk. Our measured spectrum appears to be less affected by the disk likely because we used a more aggressive stellar PSF-subtraction routine that filtered out more disk signal. We also found larger uncertainties per spectral channel on average. Our extracted spectrum of PDS 70 c is plotted in Figure 3. Both reductions have indications of correlated noise, as neighboring spectral channels have less scatter than the error bars would imply for uncorrelated spectral channels.

Zoom In Zoom Out Reset image size

With only two epochs of GRAVITY measurements for each planet, we were able to constrain the positions and velocities of PDS 70 b and c with 100 μas precision (see Table 2, but the accelerations and ultimately the orbital elements of the planet are still limited by the precision of astrometry from single-dish telescopes. We supplemented the GRAVITY astrometry with imaging astrometry from Müller et al. (2018), Mesa et al. (2019), Haffert et al. (2019), and Wang et al. (2020) (see Table 8 in Appendix B).

To define the orbit, we used the following orbital parameters: semimajor axis (a), eccentricity (e), inclination (i), argument of periastron (ω), longitude of the ascending node (Ω), epoch of periastron in units of fractional orbital period (τ), system parallax, and the component masses of each body (Blunt et al. 2020). We defined the orbital elements for each planet in Jacobi coordinates as they vary less when accounting for the effects of multiple planets (see next paragraph). The reference epoch for τ for both planets was MJD 55,000 (2009-06-18). To keep the orbits realistic, we used the same priors as Wang et al. (2020) to impose PDS 70 b and c have non-crossing orbits as well as near coplanarity of the planets and circumstellar disk. We rejected all orbits for which periastron of PDS 70 c is inside of apastron of PDS 70 b. We also applied a Gaussian prior centered at 0 with a standard deviation of 10° on the coplanarity of planet b with the disk, planet c with the disk, and planet b with planet c. We fixed the disk plane to i = 1283 (equivalent to i = 517 but for clockwise orbits) and PA = 1567 based on Keppler et al. (2019) measurements of the outer disk. The only prior we changed is the one on the stellar mass: we used a Gaussian prior centered at 0.88 M_⊙ with a standard deviation of 0.09 M_⊙ based on our stellar SED fit in Section 2.3 but with 10% errors to account for model systematics. All our priors are listed in Table 3.

Table 3. Orbital Parameters for PDS 70 b and c

Orbital Element	Prior	Near-Coplanar	Dynamically Stable
ab (au)	LogUniform(1, 100) a	${20.1}_{-1.0\,(-1.9)}^{+0.9\,(+1.9)}$	${20.8}_{-0.7\,(-1.1)}^{+0.6\,(+1.3)}$
eb	Uniform(0, 1) a	${0.22}_{-0.07\,(-0.13)}^{+0.07\,(+0.14)}$	${0.17}_{-0.06\,(-0.10)}^{+0.06\,(+0.11)}$
ib (°)	$\sin (i)$ b	${132.8}_{-3.1\,(-5.7)}^{+3.6\,(+7.7)}$	${131.0}_{-2.6\,(-4.8)}^{+2.9\,(+6.4)}$
ωb (°)	Uniform(0, 2π)	${164}_{-13\,(-25)}^{+12\,(+22)}$	${161}_{-12\,(-23)}^{+12\,(+23)}$
Ωb (°)	Uniform(0, 2π) b	${170.5}_{-4.4\,(-9.0)}^{+4.7\,(+10.4)}$	${169.7}_{-3.6\,(-6.9)}^{+4.1\,(+8.8)}$
τb	Uniform(0, 1)	${0.403}_{-0.034\,(-0.068)}^{+0.035\,(+0.074)}$	${0.402}_{-0.034\,(-0.065)}^{+0.037\,(+0.075)}$
ac (au)	LogUniform(1, 100) a	${33.2}_{-2.3\,(-4.0)}^{+2.5\,(+4.6)}$	${34.3}_{-1.8\,(-3.0)}^{+2.2\,(+4.6)}$
ec	Uniform(0, 1) a	${0.051}_{-0.035\,(-0.049)}^{+0.052\,(+0.104)}$	${0.037}_{-0.025\,(-0.035)}^{+0.041\,(+0.088)}$
ic (°)	$\sin (i)$ b	${131.6}_{-2.9\,(-5.2)}^{+3.0\,(+6.1)}$	${130.5}_{-2.4\,(-5.2)}^{+2.5\,(+4.9)}$
ωc (°)	Uniform(0, 2π)	${56}_{-24\,(-45)}^{+27\,(+60)}$	${53}_{-23\,(-43)}^{+29\,(+65)}$
Ωc (°)	Uniform(0, 2π) b	${161.0}_{-4.9\,(-10.4)}^{+4.2\,(+8.4)}$	${161.7}_{-4.3\,(-8.8)}^{+4.1\,(+8.1)}$
τc	Uniform(0, 1)	${0.51}_{-0.08\,(-0.14)}^{+0.09\,(+0.19)}$	${0.49}_{-0.07\,(-0.13)}^{+0.09\,(+0.20)}$
Parallax (mas)	${ \mathcal N }$(8.8159, 0.0405)	${8.821}_{-0.042\,(-0.081)}^{+0.041\,(+0.080)}$	${8.821}_{-0.040\,(-0.079)}^{+0.041\,(+0.084)}$
Mb (MJup)	Uniform(1, 15)	${7.9}_{-4.7\,(-6.5)}^{+4.9\,(+6.8)}$	${3.2}_{-1.6\,(-2.1)}^{+3.3\,(+8.4)}$
Mc (MJup)	Uniform(1, 15)	${7.8}_{-4.7\,(-6.4)}^{+5.0\,(+6.9)}$	${7.5}_{-4.2\,(-6.1)}^{+4.7\,(+7.0)}$
M* (M⊙)	${ \mathcal N }$(0.88, 0.09)	${0.967}_{-0.065\,(-0.126)}^{+0.065\,(+0.131)}$	${0.982}_{-0.066\,(-0.130)}^{+0.066\,(+0.128)}$

Notes. The orbital parameters used here are defined in Blunt et al. (2020) and are in Jacobi coordinates. For each parameter, the median value of the posterior is listed, with superscript and subscript describing the 68% credible interval (95% credible interval in parentheses).

^a Additional prior on periastron of c is larger than apastron of b. ^b Additional Gaussian prior on the coplanarity of b, c, and the disk.

Download table as:  ASCIITypeset image

We are sensitive to perturbations on the visual orbit of one of the planets around the star due to the other planet. Essentially, our visual orbits are defined by the relative separation between the planet and the star. A second planet, to first order, perturbs the star's position and causes the measured distance between the first planet and the star to change, creating epicycles in the visual orbit (note that this effect is different from direct planet–planet gravitational interactions, which we will not consider and are much smaller in amplitude). We followed the prescription defined in Brandt et al. (2020) where only the perturbations of inner planets are accounted for. Thus, the visual orbit of PDS 70 c relative to the star is sensitive to the orbit and mass of PDS 70 b, but not the other way around. For a Jupiter-mass planet at 20 au, the peak-to-valley amplitude of this perturbation is 400 μas. To properly model the GRAVITY astrometry, we needed to account for this effect. However, we note that with only two GRAVITY epochs per planet, we did not constrain the masses. Simply, the orbital elements we inferred would have been different if we assumed the planets were massless rather than Jovian mass. In this work, we added uniform priors on planet mass between 1 and 15 Jupiter masses for each planet. Even though recent work (e.g., Stolker et al. 2020; Wang et al. 2020) inferred masses closer to 1 M_Jup than 15 M_Jup, we purposely extended the prior range to higher masses to assess if we could rule out high-mass solutions via dynamical stability arguments.

The orbital parameters were inferred by Bayesian parameter estimation using an unreleased version of the orbitize! package with commit id 83356d9 (Blunt et al. 2020), which uses the parallel-tempered affine-invariant sampler ptemcee (Foreman-Mackey et al. 2013; Vousden et al. 2016). This version of orbitize! automatically handles the covariances of the uncertainties in R.A. and decl. that result due to the u–v coverage of the observations. We ran the sampler using 20 temperatures, 1000 walkers per temperature, and 100,000 steps per walker. Convergence was assessed by visual inspection of the walker chains, and by checking that we ran the sampler for at least 100 autocorrelation times. We also accounted for the perturbations on the visual orbits of each planet due to the other planet in the system as described in the previous paragraph. The posterior was formed using the last 40,000 steps from each walker at the lowest temperature. The visual orbit for PDS 70 b is plotted in Figure 4 and the posterior credible intervals (CIs) are listed in Table 3.

Zoom In Zoom Out Reset image size

With the new GRAVITY data, we constrained the semimajor axis of PDS 70 b (a_b ) to ±2 au (95% credible interval), but a_c to only ±5 au for the same credible interval. Perfectly circular orbits for PDS 70 b are disfavored by the current data whereas PDS 70 c is consistent with circular orbits. The period ratio of PDS 70 c to PDS 70 b is ${2.13}_{-0.24}^{+0.27}$ (68% credible interval; ${2.13}_{-0.45}^{+0.56}$ for the 95% credible interval), putting it near the 2:1 mean-motion resonance (MMR) as has been proposed by Haffert et al. (2019) and Bae et al. (2019). For the first time, we are able to strongly disfavor all other first-order mean-motion resonances such as the 3:2 and 4:3 MMR. Given these planets are thought to be Jovian mass (Stolker et al. 2020; Wang et al. 2020), if they are locked in MMR, it would likely require the strength of a first-order resonance (e.g., André & Papaloizou 2016). Thus, the 2:1 MMR is the single likely candidate for orbital resonance.

The period ratio, slight eccentricity, and masses of the planets bear strong resemblance to the HR 8799 system where at least the innermost two planets (HR 8799 d and e) harbor eccentricities near 0.1, are likely locked in a 2:1 MMR, and have a period ratio slightly larger than 2. Wang et al. (2018) proposed that HR 8799 d and e arrived at this orbital configuration due to resonant migration in the protoplanetary disk: radial migration in MMR excites the planets' eccentricity while eccentricity damping due to the viscous circumstellar disk repelled the planets to period ratios greater than 2.

We investigated whether PDS 70 b and c are in a similar dynamical scenario by searching for dynamically stable orbits. Following the procedure in Wang et al. (2018), we performed rejection sampling on our posterior of orbital parameters by imposing a stability prior that the orbital configuration is stable for the system's age of 8 Myr, noting our results are not extremely sensitive to the exact choice of age. For each orbit in the posterior, we used the REBOUND N-body package with the IAS15 integrator to advance the system backward in time for 8 Myr (Rein & Liu 2012; Rein & Spiegel 2015). We assigned component masses for all three bodies based on the mass posteriors from our orbit fit. Note that only the mass of the star was constrained by our orbit fit, and that the posterior masses from the two protoplanets were dictated by the prior. We considered a system unstable if the two planets pass within one mutual hill radius of each other (∼3.3 au) or if any planet is ejected to 500 au. During the simulations, we logged the 2:1 resonance angle between the two protoplanets that we define as

$$\begin{eqnarray}&&{\theta }_{c:b}={\lambda }_{b}-2{\lambda }_{c}+{\varpi }_{b},\end{eqnarray} \tag{ 7 }$$

where ϖ = Ω + ω is the longitude of periastron, and λ = ϖ + M is the mean longitude (M is the mean anomaly). We used the same algorithm as Wang et al. (2018) to identify where in time series of θ_c:b is it librating or circulating (see their Figure 7), and saved the fraction of time the angle was librating over 8 Myr in each simulation.

In these simulations that just account for the gravitational interactions of the three bodies, we found that 41% of the orbital posterior is dynamically stable and only 3% of the stable orbits have the two planets in resonance lock (where θ_c:b is librating >95% of the time). In fact, the majority of stable orbits had θ_c:b circulating the entire time, indicating that the planets were not in resonance for any significant period of time in the simulations. The relatively small fraction of stable orbits in resonance lock is likely due to the significant uncertainties in the orbital parameters, with many combinations of orbital parameters lying well outside of any region with MMR can occur. This difficulty in finding resonant orbits has also been seen in HR 8799 (Wang et al. 2018).

However, we note that we have not included planet-disk interactions and gas drag, which could affect the planets' orbits. It also does not rule out that these planets could in the future migrate into orbital resonance. Thus, we will avoid investigating the detailed dynamical interactions of the system with our simulations alone. Encouragingly, Bae et al. (2019) accounted for these effects and showed that having planets in the approximate orbital configuration of PDS 70 b and c migrate into resonance while accreting from the circumstellar disk would create a circumstellar disk and gap that is consistent with the millimeter wavelength observations and imply mass accretion rates consistent with the Hα luminosities. Furthermore, they predicted that such a migration into resonance would pump the eccentricity of PDS 70 b to ∼0.1–0.2, which agrees very well with our inferred eccentricity of 0.17 ± 0.06 for PDS 70 b assuming dynamical stability. The current observations are thus consistent with planets being in 2:1 resonance. However, we cannot reject other scenarios that do not require the planets to be in resonance at this time.

We used the dynamical stability prior to place upper limits on the masses of the two planets, and plot the 1D marginalized posterior distribution of their masses in Figure 5. We found nearly no constraint on the mass of PDS 70 c from enforcing stability, but we found a 95% upper limit of 10 M_Jup for PDS 70 b, consistent with the masses predicted by Wang et al. (2020).

Zoom In Zoom Out Reset image size

The mass of the host star is also constrained by the orbital motion of the planets. Given that the masses of young stars are more difficult to constrain from photometry or spectroscopy alone compared to main-sequence stars, the dynamical mass constraints on the host star is another piece of useful information from the orbit fit. We found a stellar mass of 0.982 ± 0.066 M_⊙ in our dynamically stable solutions. Compared to the dynamical mass estimate of 0.875 ± 0.03 M_⊙ from fitting velocity maps of the circumstellar gas (Keppler et al. 2019) and the model-dependent mass estimate of 0.88 ± 0.02 M_⊙ from the stellar SED fit described in Section 2.3, our dynamical mass estimate from the planetary orbits is systematically high, although it is consistent at the 1.6σ level. Due to the short orbital arc, there remains ∼10% uncertainties in our dynamical mass estimate from the orbital motion of the planets, since orbital period, semimajor axis, and stellar mass are degenerate. Because of this, using our dynamical mass posterior as a prior in our stellar SED fits results in no change in the derived stellar spectrum. If we instead fix the stellar mass to 1 M_⊙ in our SED fits, the K-band excess predicted by SED fits drops to 10%, but the J and H band photometry become 2σ discrepant with the model. Extending the orbital coverage with more astrometric monitoring will improve our stellar mass estimate.

We used the following literature measurements in the fits along with our GRAVITY K-band spectrum: VLT/SPHERE YJH spectrum at R ∼ 30 (Müller et al. 2018), VLT/SPHERE photometry at H and K bands (Müller et al. 2018), and 3–5 μm photometry from Keck/NIRC2, Gemini/NICI, and VLT/NACO (Müller et al. 2018; Stolker et al. 2020; Wang et al. 2020). All of the literature photometry we used are listed in Table 7 in Appendix B. We excluded fitting the VLT/SINFONI spectrum from Christiaens et al. (2019a) as it disagrees with both the GRAVITY spectrum and the SPHERE photometry at the same wavelengths by being ∼30% brighter. Whether this is astrophysical variability (the SINFONI data was taken ∼4 yr earlier) or instrumental systematics is uncertain at this point, so we did not consider it here for simplicity.

With only a tentative water absorption band between the J and H bands, Wang et al. (2020) found that a single blackbody was the most justified model. However, our GRAVITY spectrum shows a dip at the blue end of the K band that is consistent with the water-absorption band seen in substellar atmospheres. To perform this test quantitatively, we performed Bayesian model comparison between the different fits. We fit each model using the same Bayesian framework as Wang et al. (2020). We used a Gaussian process with the same square exponential kernel to empirically estimate the correlated noise in the SPHERE YJH spectrum when fitting the atmospheric models to the data. The GRAVITY spectrum has its covariance estimated as part of the data reduction and we used this covariance matrix when accounting for its correlated noise in the likelihood.

For the priors, we picked uniform priors in effective temperature (T_eff) between 1000 K and 1500 K and uniform priors in effective radius (R) between 0.5 and 5 Jupiter radii. For grids with surface gravity (log(g)), metalicity ([M/H]), and carbon to oxygen ratio (C/O), we used uniform priors with the bounds spanned by the edges of the model grids: for BT-SETTL, $3.5\lt \mathrm{log}(g)\lt 5.5;$ for DRIFT-PHOENIX, $3.0\lt \mathrm{log}(g)\lt 5.5$ and −0.3 < [M/H] < 0.3; for Exo-REM, $3.0\lt \mathrm{log}(g)\lt 4.5$, −0.5 < [M/H] < 0.5, and 0.3 < C/O < 0.75. We used pymultinest to sample the posterior distribution and numerically compute the evidence of each model (Buchner et al. 2014). The median and 95% credible intervals of each parameter are listed in Table 4 in the "Plain Models" section. The evidence allows us to compute the Bayes factor B to test the relative probability of two models:

$$\begin{eqnarray}&&B\equiv \displaystyle \frac{P({M}_{1}| D)}{P({M}_{2}| D)}=\displaystyle \frac{P(D| {M}_{1})P({M}_{1})}{P(D| {M}_{2})P({M}_{2})}.\end{eqnarray} \tag{ 8 }$$

Table 4. Model Fits to SED of PDS 70 b

Model	Teff (K)	R (RJup)	$\mathrm{log}(g)$	[M/H]	C/O	AV (mag)	${a}_{\max }$ (μm)	βdust	B
Plain Models
Blackbody	${1244}_{-41}^{+50}$	${2.47}_{-0.24}^{+0.24}$	⋯	⋯	⋯	⋯	⋯	⋯	1
BT-SETTL	${1252}_{-19}^{+22}$	${1.97}_{-0.03}^{+0.05}$	${3.51}_{-0.01}^{+0.04}$ a	⋯	⋯	⋯	⋯	⋯	5.5 × 10−6
DRIFT-PHOENIX	${1384}_{-48}^{+47}$	${1.85}_{-0.16}^{+0.24}$	${5.25}_{-1.43}^{+0.24}$	$-{0.15}_{-0.15}^{+0.36}$	⋯	⋯	⋯	⋯	96
Exo-REM	${1051}_{-12}^{+17}$	${3.62}_{-0.44}^{+0.15}$	${3.50}_{-0.00}^{+0.01}$ a	${0.04}_{-0.05}^{+0.46}$	${0.70}_{-0.05}^{+0.01}$	⋯	⋯	⋯	5.3 × 10−9
ISM Extinction
Blackbody	${1352}_{-109}^{+137}$	${2.16}_{-0.30}^{+0.37}$	⋯	⋯	⋯	${2.9}_{-2.7}^{+3.7}$	⋯	⋯	0.69
BT-SETTL	${1418}_{-106}^{+69}$	${1.98}_{-0.14}^{+0.19}$	${3.87}_{-0.36}^{+0.77}$	⋯	⋯	${6.1}_{-2.2}^{+3.3}$	⋯	⋯	84
DRIFT-PHOENIX	${1442}_{-71}^{+52}$	${1.85}_{-0.18}^{+0.25}$	${4.62}_{-1.21}^{+0.84}$	$-{0.12}_{-0.17}^{+0.36}$	⋯	${3.0}_{-2.7}^{+4.7}$	⋯	⋯	59
Exo-REM	${1234}_{-78}^{+95}$	${2.88}_{-0.35}^{+0.44}$	${3.89}_{-0.35}^{+0.42}$	${0.01}_{-0.37}^{+0.44}$	${0.61}_{-0.21}^{+0.12}$	${8.3}_{-2.2}^{+1.5}$	⋯	⋯	46
Power-Law Dust Extinction
Blackbody	${1273}_{-55}^{+71}$	${2.34}_{-0.29}^{+0.27}$	⋯	⋯	⋯	${2.9}_{-2.8}^{+6.4}$	${0.102}_{-0.078}^{+0.592}$	$-{2.4}_{-2.8}^{+2.3}$	0.14
BT-SETTL	${1402}_{-89}^{+71}$	${2.53}_{-0.80}^{+0.74}$	${3.89}_{-0.35}^{+0.53}$	⋯	⋯	${3.5}_{-2.6}^{+4.3}$	${2.41}_{-2.13}^{+6.26}$	$-{3.9}_{-1.1}^{+3.1}$	3.8
DRIFT-PHOENIX	${1413}_{-59}^{+54}$	${1.80}_{-0.18}^{+0.27}$	${5.04}_{-1.58}^{+0.43}$	$-{0.07}_{-0.21}^{+0.33}$	⋯	${3.1}_{-2.9}^{+6.3}$	${0.097}_{-0.073}^{+0.369}$	$-{2.4}_{-2.7}^{+2.2}$	20
Exo-REM	${1223}_{-79}^{+105}$	${2.81}_{-0.45}^{+0.52}$	${3.84}_{-0.31}^{+0.38}$	${0.00}_{-0.39}^{+0.43}$	${0.61}_{-0.17}^{+0.12}$	${4.7}_{-2.9}^{+4.1}$	${0.70}_{-0.28}^{+6.49}$	$-{4.1}_{-1.1}^{+3.5}$	1.1
Second Blackbody							T2 (K)	R2 (RJup)
Blackbody	${1252}_{-40}^{+48}$	${2.42}_{-0.23}^{+0.24}$	⋯	⋯	⋯	⋯	${351}_{-225}^{+190}$	${24.5}_{-19.7}^{+23.3}$	0.53
BT-SETTL	${1266}_{-81}^{+68}$	${1.34}_{-0.42}^{+0.36}$	${3.88}_{-0.35}^{+0.60}$	⋯	⋯	⋯	${1037}_{-182}^{+107}$	${3.1}_{-0.7}^{+1.3}$	15
DRIFT-PHOENIX	${1393}_{-45}^{+54}$	${1.81}_{-0.19}^{+0.21}$	${5.30}_{-1.28}^{+0.19}$	$-{0.11}_{-0.18}^{+0.38}$	⋯	⋯	${301}_{-189}^{+147}$	${28.4}_{-24.2}^{+20.1}$	38
Exo-REM	${1121}_{-59}^{+80}$	${1.44}_{-0.39}^{+0.59}$	${4.16}_{-0.50}^{+0.31}$	${0.14}_{-0.47}^{+0.32}$	${0.54}_{-0.22}^{+0.16}$	⋯	${1078}_{-83}^{+59}$	${3.1}_{-0.4}^{+0.6}$	0.44
ISM Extinction + Second Blackbody
Blackbody	${1367}_{-119}^{+123}$	${2.11}_{-0.27}^{+0.38}$	⋯	⋯	⋯	${3.1}_{-2.9}^{+3.2}$	${369}_{-241}^{+158}$	${26.5}_{-21.4}^{+21.7}$	0.34
BT-SETTL	${1392}_{-82}^{+82}$	${1.96}_{-0.17}^{+0.20}$	${3.83}_{-0.32}^{+0.69}$	⋯	⋯	${5.4}_{-1.9}^{+3.4}$	${407}_{-128}^{+265}$	${26.6}_{-21.6}^{+21.5}$	147
DRIFT-PHOENIX	${1449}_{-54}^{+47}$	${1.84}_{-0.18}^{+0.26}$	${5.11}_{-1.67}^{+0.38}$	$-{0.11}_{-0.18}^{+0.33}$	⋯	${3.0}_{-2.2}^{+4.5}$	${384}_{-263}^{+323}$	${3.8}_{-1.9}^{+4.8}$	2.2
Exo-REM	${1235}_{-76}^{+91}$	${2.86}_{-0.33}^{+0.34}$	${3.89}_{-0.36}^{+0.36}$	${0.03}_{-0.32}^{+0.43}$	${0.60}_{-0.21}^{+0.13}$	${8.0}_{-2.2}^{+1.7}$	${318}_{-202}^{+207}$	${24.7}_{-20.2}^{+23.4}$	14
Accreting CPD							${M}_{p}\dot{M}$ (${M}_{\mathrm{Jup}}^{2}$/yr)	Rin (RJup)
Blackbody	${1257}_{-43}^{+50}$	${2.38}_{-0.26}^{+0.24}$	⋯	⋯	⋯	⋯	${3.3}_{-2.1}^{+6.4}\times {10}^{-7}$	${2.7}_{-1.2}^{+1.1}$	0.24
BT-SETTL	${1259}_{-26}^{+126}$	${1.68}_{-0.25}^{+0.12}$	${3.55}_{-0.05}^{+0.18}$ a	⋯	⋯	⋯	${2.9}_{-0.8}^{+11.5}\times {10}^{-7}$	${1.1}_{-0.1}^{+2.2}$	0.55
DRIFT-PHOENIX	${1409}_{-47}^{+45}$	${1.76}_{-0.17}^{+0.23}$	${5.26}_{-1.33}^{+0.22}$	${0.02}_{-0.30}^{+0.26}$	⋯	⋯	${3.2}_{-2.1}^{+6.0}\times {10}^{-7}$	${2.9}_{-1.2}^{+1.1}$	14
Exo-REM	${1076}_{-43}^{+77}$	${2.59}_{-0.42}^{+0.38}$	${3.52}_{-0.02}^{+0.11}$ a	${0.26}_{-0.29}^{+0.22}$	${0.66}_{-0.16}^{+0.05}$	⋯	${2.8}_{-0.7}^{+2.7}\times {10}^{-7}$	${1.1}_{-0.1}^{+0.2}$	2.8 × 10−5

Note. For each parameter, a 95% credible interval centered about the median is reported. The superscript and subscript denote the upper and lower bounds of that range.

^a Mode of posterior reached edge of model grid.

Download table as:  ASCIITypeset image

In this equation, P is the probability of a quantity, M₁ and M₂ are the two models that are being compared, and D is the data. The left side is the relative probability of M₁ compared to M₂ given the current data. On the right side, P(D∣M) is the evidence of a given model, and P(M) is the prior probability of a given model. Assuming equal weight for all models, as we do not think one model is better justified than any of the others, the Bayes factor of two models is equal to the ratio of evidences. We benchmarked all of the models we considered against the simple blackbody model (i.e., we set it as M₂) given it has been the preferred model in previous work (Stolker et al. 2020; Wang et al. 2020). We list the values of B for each model relative to the plain blackbody model in the rightmost column of Table 4.

Given the accreting nature of PDS 70 b, it is possible that the planet is extincted by accreting materials or by circumplanetary and circumstellar dust. Indeed previous atmosphere modeling indicated that both planets should be shrouded by its dust (Stolker et al. 2020; Wang et al. 2020). The emission we observed could be a planetary atmosphere attenuated by obscuring dust. Although it is not entirely accurate, we first considered a simple interstellar medium (ISM) extinction law. An ISM extinction law has been shown to be an adequate approximation of stars shrouded by their circumstellar disks, so it is not unreasonable (Looper et al. 2010a, 2010b). We used an extinction law derived for stars attenuated by the interstellar medium in the near-infrared by Wang & Chen (2019). This extinction law follows the form of

$$\begin{eqnarray}&&{A}_{\lambda }={A}_{V}{\left(\displaystyle \frac{\lambda }{{\lambda }_{V}}\right)}^{-\beta },\end{eqnarray} \tag{ 9 }$$

where A_λ is the magnitudes of extinction at a wavelength λ, A_V is the extinction in the V band with center wavelength λ_V = 0.55 μm, and β is the power-law index of 2.07 derived by Wang & Chen (2019). As we fixed the power-law index, A_V is the only new free variable introduced. We placed a uniform prior on A_V between 0 and 10 mags. We repeated the fit of the four model grids, but now with the model flux attenuated by

$$\begin{eqnarray}&&{F}_{\mathrm{obs}}={10}^{-{A}_{\lambda }/2.5}{F}_{\mathrm{emit}},\end{eqnarray} \tag{ 10 }$$

where F_obs is the observed flux and F_emit is the original flux from the model grids. We recorded the best-fit parameters and B relative to the plain blackbody model with no extinction in Table 4 under "ISM Extinction."

Given that the ISM extinction law may not be fully representative of accreting dust in a circumplanetary environment where we expect grain growth (e.g., Birnstiel et al. 2012; Kataoka et al. 2013; Piso et al. 2015), a more general case where the dust particles follow a variable power law in grain sizes may better describe the data. Such dust extinction prescriptions have been shown to fit dusty free-floating brown dwarfs (Marocco et al. 2014; Hiranaka et al. 2016) and directly imaged companions (Bonnefoy et al. 2016; Delorme et al. 2017). Thus, we considered replacing the ISM extinction law with a power-law dust extinction prescription. We assumed MgSiO₃ dust with particle size distribution $n\propto {a}^{{\beta }_{\mathrm{dust}}}$ with a being the radius of the dust, and β_dust being the power-law exponent. We set a minimum dust radius of 1 nm, and vary the maximum dust radius ${a}_{\max }$, as the minimum grain size does not significantly impact the spectrum. For a given ${a}_{\max }$ and β_dust, we computed the extinction cross section (σ_dust) of the dust as a function of wavelength with PyMieScatt (Sumlin et al. 2018) by using the refractive indices from Scott & Duley (1996) and Jaeger et al. (1998). To relate the cross-sectional area to an amount of attenuated flux per wavelength, we used the relation

$$\begin{eqnarray}&&{F}_{\mathrm{obs}}={e}^{-{\tau }_{\mathrm{dust}}}{F}_{\mathrm{emit}}={10}^{-{A}_{V}/2.5}{e}^{-\tfrac{{\sigma }_{\mathrm{dust}}}{{\sigma }_{\mathrm{dust},{\rm{V}}}}}{F}_{\mathrm{emit}},\end{eqnarray} \tag{ 11 }$$

where σ_dust,V is the cross-sectional absorption area averaged across the V band, and τ_dust is the optical depth of the dust. Our uniform prior on ${a}_{\max }$ was between 0.01 and 10 μm and and our uniform prior on β_dust was between −10 and 0. Our prior on A_V remained between 0 and 10 mags.

We also considered augmenting the forward models with circumplanetary disk (CPD) models. We first used a simple blackbody component to model the CPD as has been considered in the past work (Stolker et al. 2020; Wang et al. 2020). CPD models have indicated that the bulk of the thermal emission from a circumplanetary disk would come from the inner edge of the disk (Zhu 2015; Szulágyi et al. 2019). For the second blackbody, we adopted priors for the temperature of the second blackbody component (T₂ and R₂, respectively) that are motivated by these modeling studies. The T₂ prior was a uniform prior between 100 K and T_eff, the effective temperature of the first component. The R₂ prior was a uniform prior between R, the effective radius of the first component, and 50 R_Jup. These priors are not uniform in T₂ and R₂, but if we marginalized over T_eff and R, we get priors that only weakly favor lower temperatures and larger radii.

Since extinction of the planetary atmosphere may play a significant role, we considered the case of an extincted atmosphere model plus a second blackbody component, similar to what was done in Christiaens et al. (2019a). This case attempts to model circumplanetary dust absorbing the light from the protoplanet and reradiating it away at longer wavelengths. We used the simple ISM extinction law, as it has fewer free parameters, even though we note that the slope may not be perfectly accurate for circumplanetary dust. We only applied the extinction to the atmospheric model and not the second blackbody component. We used the same priors on A_V , T₂, and R₂ as previously.

Last, we considered using the more sophisticated accreting CPD model from Zhu (2015) that model the emission from a CPD with density and temperature gradients and account for molecular and atomic opacities. The resulting spectra are parameterized by the product of the planet's mass and its mass-accretion rate (${M}_{p}\dot{M}$) and the inner edge of the CPD (R_in). The spectra are only weakly sensitive to the outer disk edge. Zhu (2015) produced models with the outer radius being 50 and 1000 R_in, and we marginalized over the two outer radii in our SED fits, as there was no statistically significant difference.

In all, we tried six different modifications to the four forward models, resulting in 24 models. We plotted the best-fit model for the model modification with the highest Bayes factor for each of the four atmospheric forward models in Figure 6 and listed all the results in Table 4. We discussed the model selection and implications further in Section 4.3.

Zoom In Zoom Out Reset image size

In addition to the GRAVITY K-band spectrum and the re-extracted SPHERE IFS YJH spectrum, we also used the K-band photometry from Mesa et al. (2019) and L-band photometry from Wang et al. (2020). We listed the exact numbers for the literature photometry that we used in Table 7 in Appendix B. In this systematic exploration of models, we did not fit the 855 μm continuum emission coming from the location of PDS 70 c (Isella et al. 2019), as the Exo-REM and accreting CPD grids did not extend to those wavelengths, and primarily focused on the 1–5 μm SED where the bulk of the planetary emission should be (see the end of the section for more fits with this data point).

We used the same four base forward models. We modified the prior for T_eff of the blackbody models to be between 700 to 1200 K instead, as the previous prior range for PDS 70 b was too high. We did not modify the T_eff for the other forward models because they did not go below 1000 K. For all four models, the range of T_eff remained 500 K, so the impact of a different T_eff prior on the evidence of the blackbody models should be negligible.

We also used the same extinction and CPD modifications as for PDS 70 b. The only change was changing the prior limits for A_V to be between 0 and 20 mags instead of 0 and 10 mags, as preliminary analysis indicated the extinction could be greater than 10 mags. Increasing the prior range on A_V may decrease the evidence of the models with extinction slightly, but we accepted this in order to have a more flexible extinction prescription.

The 95% CI centered about the median of each parameter of each model fit along with the Bayes factor of each model relative to the single blackbody model with no modifications are listed in Table 5. The best-fit spectrum of the model modification with the highest Bayes factor for each of the four base forward models are plotted in Figure 7.

Zoom In Zoom Out Reset image size

Table 5. Model Fits to SED of PDS 70 c

Model	Teff (K)	R (RJup)	$\mathrm{log}(g)$	[M/H]	C/O	AV (mag)	${a}_{\max }$ (μm)	βdust	B
Plain Models
Blackbody	${1005}_{-58}^{+65}$	${2.39}_{-0.42}^{+0.53}$	⋯	⋯	⋯	⋯	⋯	⋯	1
BT-SETTL	${1155}_{-8}^{+37}$	${1.27}_{-0.12}^{+0.05}$	${3.50}_{-0.00}^{+0.01}$ a	⋯	⋯	⋯	⋯	⋯	3.4 × 10−11
DRIFT-PHOENIX	${1054}_{-48}^{+60}$	${1.98}_{-0.31}^{+0.39}$	${4.54}_{-0.85}^{+0.61}$	$-{0.22}_{-0.07}^{+0.46}$	⋯	⋯	⋯	⋯	6.2 × 107
Exo-REM	${1024}_{-22}^{+30}$	${1.89}_{-0.21}^{+0.26}$	${3.50}_{-0.00}^{+0.02}$ a	${0.47}_{-0.15}^{+0.03}$ a	${0.65}_{-0.02}^{+0.02}$	⋯	⋯	⋯	3.1 × 10−17
ISM Extinction
Blackbody	${1401}_{-289}^{+92}$	${1.40}_{-0.16}^{+0.65}$	⋯	⋯	⋯	${14.1}_{-8.7}^{+3.8}$	⋯	⋯	19
BT-SETTL	${1289}_{-104}^{+99}$	${1.60}_{-0.19}^{+0.30}$	${4.10}_{-0.54}^{+0.54}$	⋯	⋯	${18.1}_{-4.0}^{+1.8}$	⋯	⋯	3.6 × 106
DRIFT-PHOENIX	${1086}_{-68}^{+128}$	${1.97}_{-0.49}^{+0.44}$	${5.05}_{-0.81}^{+0.41}$	$-{0.11}_{-0.17}^{+0.37}$	⋯	${5.3}_{-4.4}^{+5.1}$	⋯	⋯	5.5 × 107
Exo-REM	${1112}_{-98}^{+133}$	${2.42}_{-0.54}^{+0.66}$	${3.87}_{-0.30}^{+0.48}$	${0.01}_{-0.42}^{+0.44}$	${0.50}_{-0.19}^{+0.20}$	${18.6}_{-3.2}^{+1.3}$	⋯	⋯	2.2 × 106
Power-Law Dust Extinction
Blackbody	${1338}_{-279}^{+150}$	${1.23}_{-0.24}^{+0.88}$	⋯	⋯	⋯	${16.0}_{-8.5}^{+3.8}$	${0.31}_{-0.24}^{+0.24}$	$-{2.8}_{-2.3}^{+2.4}$	13
BT-SETTL	${1327}_{-115}^{+148}$	${1.24}_{-0.17}^{+0.51}$	${4.38}_{-0.83}^{+0.84}$	⋯	⋯	${15.2}_{-8.5}^{+4.5}$	${0.41}_{-0.12}^{+2.79}$	$-{2.5}_{-2.6}^{+2.3}$	4.3 × 105
DRIFT-PHOENIX	${1074}_{-64}^{+78}$	${1.92}_{-0.42}^{+0.50}$	${4.73}_{-0.72}^{+0.66}$	$-{0.14}_{-0.15}^{+0.39}$	⋯	${11.3}_{-10.3}^{+8.0}$	${0.089}_{-0.070}^{+0.595}$	$-{2.5}_{-3.2}^{+2.3}$	2.1 × 107
Exo-REM	${1118}_{-104}^{+130}$	${1.96}_{-0.50}^{+0.55}$	${3.92}_{-0.35}^{+0.50}$	$-{0.15}_{-0.33}^{+0.57}$	${0.47}_{-0.15}^{+0.22}$	${16.3}_{-8.1}^{+3.4}$	${0.37}_{-0.08}^{+0.42}$	$-{2.0}_{-2.8}^{+1.9}$	2.9 × 105
Second Blackbody							T2 (K)	R2 (RJup)
Blackbody	${1024}_{-22}^{+51}$	${2.25}_{-0.31}^{+0.18}$	⋯	⋯	⋯	⋯	${268}_{-153}^{+190}$	${21.9}_{-18.5}^{+25.5}$	0.12
BT-SETTL	${1152}_{-33}^{+47}$	${0.80}_{-0.18}^{+0.19}$	${3.52}_{-0.02}^{+0.09}$ a	⋯	⋯	⋯	${811}_{-104}^{+72}$	${4.0}_{-0.9}^{+2.4}$	4.1 × 10−3
DRIFT-PHOENIX	${1054}_{-48}^{+63}$	${1.96}_{-0.30}^{+0.30}$	${4.46}_{-0.76}^{+0.48}$	$-{0.19}_{-0.10}^{+0.41}$	⋯	⋯	${269}_{-154}^{+238}$	${22.4}_{-18.9}^{+24.5}$	1.6 × 107
Exo-REM	${1064}_{-54}^{+68}$	${0.62}_{-0.11}^{+0.24}$	${3.94}_{-0.37}^{+0.38}$	${0.22}_{-0.49}^{+0.25}$	${0.53}_{-0.21}^{+0.19}$	⋯	${878}_{-66}^{+54}$	${3.5}_{-0.7}^{+1.1}$	3.1 × 10−5
ISM Extinction + Second Blackbody
Blackbody	${1186}_{-155}^{+122}$	${1.80}_{-0.38}^{+0.59}$	⋯	⋯	⋯	${7.8}_{-6.2}^{+2.1}$	${265}_{-152}^{+204}$	${23.4}_{-20.2}^{+24.1}$	1.3
BT-SETTL	${1291}_{-91}^{+101}$	${1.59}_{-0.20}^{+0.23}$	${4.11}_{-0.54}^{+0.53}$	⋯	⋯	${17.8}_{-4.3}^{+2.0}$	${278}_{-167}^{+238}$	${21.4}_{-18.3}^{+26.2}$	8.3 × 105
DRIFT-PHOENIX	${1088}_{-65}^{+121}$	${1.90}_{-0.42}^{+0.37}$	${4.91}_{-0.73}^{+0.46}$	$-{0.02}_{-0.25}^{+0.29}$	⋯	${5.5}_{-4.8}^{+4.7}$	${252}_{-136}^{+207}$	${24.8}_{-20.8}^{+23.1}$	9.6 × 106
Exo-REM	${1123}_{-104}^{+117}$	${2.36}_{-0.51}^{+0.67}$	${3.92}_{-0.34}^{+0.44}$	${0.02}_{-0.34}^{+0.43}$	${0.49}_{-0.17}^{+0.20}$	${18.5}_{-3.2}^{+1.4}$	${261}_{-149}^{+205}$	${21.6}_{-17.0}^{+26.6}$	5.4 × 105
Accreting CPD							${M}_{p}\dot{M}$ (${M}_{\mathrm{Jup}}^{2}$/yr)	Rin (RJup)
Blackbody	${1026}_{-24}^{+56}$	${2.22}_{-0.35}^{+0.20}$	⋯	⋯	⋯	⋯	${2.2}_{-1.1}^{+6.6}\times {10}^{-7}$	${3.1}_{-1.4}^{+0.8}$	0.038
BT-SETTL	${1153}_{-27}^{+50}$	${0.81}_{-0.22}^{+0.17}$	${3.52}_{-0.02}^{+0.09}$ a	⋯	⋯	⋯	${3.1}_{-1.1}^{+3.4}\times {10}^{-7}$	${1.3}_{-0.15}^{+0.13}$	1.8 × 10−4
DRIFT-PHOENIX	${1059}_{-53}^{+57}$	${1.89}_{-0.27}^{+0.28}$	${4.40}_{-0.83}^{+0.53}$	$-{0.18}_{-0.12}^{+0.42}$	⋯	⋯	${2.6}_{-1.5}^{+6.6}\times {10}^{-7}$	${3.0}_{-1.3}^{+0.9}$	3.4 × 106
Exo-REM	${1051}_{-47}^{+81}$	${0.63}_{-0.12}^{+0.25}$	${3.79}_{-0.27}^{+0.37}$	${0.29}_{-0.-54}^{+0.20}$	${0.55}_{-0.23}^{+0.18}$	⋯	${2.7}_{-0.7}^{+11.2}\times {10}^{-7}$	${1.2}_{-0.2}^{+2.2}$	3.2 × 10−4

Note. For each parameter, a 95% credible interval centered about the median is reported. The superscript and subscript denote the upper and lower bounds of that range.

^a Mode of posterior reached edge of model grid.

Download table as:  ASCIITypeset image

Given that Isella et al. (2019) used the 855 μm detection to demonstrate the existence of a CPD disk around PDS 70 c, we ran a few fits including this photometric point to verify this conclusion and characterize the CPD. As baseline models, we repeated the single and two blackbody fits with this longer wavelength measurement. From the fits above to the 1–5 μm data, we found that the model with the most support from the data was the plain DRIFT-PHOENIX, and that augmenting it with a cooler blackbody component was acceptable (see Table 5 and Section 4.3). We thus also refit the plain DRIFT-PHOENIX model and the DRIFT-PHOENIX model supplemented with a cooler blackbody component. In the fits, we extended the upper limit on the prior for R₂ to 5000 R_Jup (2.4 au) and the lower limit of T₂ to 10 K based on the results from (Isella et al. 2019) and Wang et al. (2020) that point to a very large and cold CPD. Since the data used in the fit changed, we avoided direct model comparisons between these fits and the fits to only the 1–5 μm data, and only compared these four models among themselves. We define B₈₅₅ as the Bayes factor between one of these models and the plain blackbody fit that includes the 855 μm data point. We list the results of the model fits in Table 6.

Table 6. Model Fits to SED of PDS 70 c Including 855 μm Photometry

Model	Teff (K)	R (RJup)	$\mathrm{log}(g)$	[M/H]	T2 (K)	R2 (RJup)	B855
Plain Models
Blackbody	${1005}_{-57}^{+55}$	${2.39}_{-0.37}^{+0.52}$	⋯	⋯	⋯	⋯	1
DRIFT-PHOENIX	${1051}_{-44}^{+62}$	${1.99}_{-0.32}^{+0.33}$	${4.53}_{-0.91}^{+0.57}$	$-{0.21}_{-0.08}^{+0.40}$	⋯	⋯	7.0 × 107
2nd Blackbody
Blackbody	${1021}_{-20}^{+49}$	${2.27}_{-0.31}^{+0.17}$	⋯	⋯	${119}_{-105}^{+150}$	${430}_{-168}^{+1198}$	1.7 × 104
DRIFT-PHOENIX	${1055}_{-48}^{+58}$	${1.96}_{-0.29}^{+0.28}$	${4.46}_{-0.79}^{+0.45}$	$-{0.19}_{-0.10}^{+0.41}$	${125}_{-94}^{+132}$	${419}_{-151}^{+518}$	1.7 × 1012

Note. For each parameter, a 95% credible interval centered about the median is reported. The superscript and subscript denote the upper and lower bounds of that range.

Download table as:  ASCIITypeset image

For the purposes of model selection, we denoted any model within a Bayes factor of 100 of the best fitting model (i.e., highest Bayes factor) to be "adequate." The relative probability of adequate models are >1% compared to the best-fitting model, which we considered good enough to not be excluded. First, we discuss the fits that only consider the 1–5 μm data. For PDS 70 b, we found that the BT-SETTL model modified with both extinction and a second blackbody component has the most support from the data. For PDS 70 c, the plain DRIFT-PHOENIX model has the highest Bayes factor by being able to fit the data the best without unnecessary free parameters. Thus we considered models with B > 1.5 and B > 6 × 10⁵ to be adequate for PDS 70 b and c, respectively.

Based on the Bayes factors, the new GRAVITY K-band spectra are able to reject the pure blackbody model for the photosphere of both protoplanets in favor of the three planetary atmosphere models. In particular, the falling slopes in both the short- and long-wavelength ends of the K band are incompatible with blackbody predictions (see inset of Figures 6 and 7), and require opacity sources such as water, molecular hydrogen, and carbon monoxide absorption to create the observed slopes in the GRAVITY spectra. For PDS 70 b, this corroborates the tentative 1.4 μm water absorption feature seen in the SPHERE IFS data (Müller et al. 2018). The difference in Bayes factors is far steeper for PDS 70 c. This appears to be due to the fact that the slope of the GRAVITY K-band spectrum for PDS 70 c is in much starker disagreement with the predictions made by the blackbody model. Thus, we will mainly focus on the three planetary atmosphere models, as all three are adequate fits given the appropriate modifications.

The plain DRIFT-PHOENIX model, in addition to being the most favored model for PDS 70 c, is the model with the second highest support for PDS 70 b. Adding modifications to the DRIFT-PHOENIX model did not improve the fit, resulting in lower Bayes factors. This can also be seen in the range of A_V , T₂, and R₂ parameters derived in the fits with modifications. The ranges of these parameters are typically consistent with the lower bounds of the priors for these parameters, implying they are minimally altering the DRIFT-PHOENIX spectrum.

The BT-SETTL and Exo-REM models, on the other hand, are poor fits to the data without modifications, with a Bayes factor orders of magnitude worse than both the plain DRIFT-PHOENIX and blackbody models. However, adding some sort of extinction to change the overall 1–4 μm slope drastically improved their fit, pulling their Bayes factor to within a factor of 100 of the best-fitting model. The ISM extinction amplitude of 3.9 < A_V < 9.4 mag for BT-SETTL fits to PDS 70 b corresponds to an 0.23 < A_K < 0.54 mag, which is similar to the extinction values found for dusty brown dwarfs (Marocco et al. 2014; Delorme et al. 2017).

Switching from ISM extinction with one free parameter to a variable power-law dust extinction with three free parameters caused drops in the Bayes factor in nearly all cases, except for the Exo-REM model of PDS 70 c (although this model's B was too low to be considered adequate). We do not think this implies that we are seeing extinction from ISM-like grains, but rather that the current data are insufficient to constrain more-flexible extinction models. In all cases, we ruled out extreme size distributions with β_dust < −5 that are dominated solely by small particles. While the maximum dust size (${a}_{\max }$) is relatively unconstrained for PDS 70 b, most of our fits ruled out dust particles larger than about 1 μm for PDS 70 c. However, we note that only the DRIFT-PHOENIX model with power-law dust extinction has an adequate B for PDS 70 c, so it is unclear how robust this conclusion is. Rather, we are worried that the free parameters in the model are compensating for other model deficiencies. Overall, the lack of improvement in the B indicates the current data is unable to characterize the properties of any obscuring dust.

The DRIFT-PHOENIX and Exo-REM models have free parameters to describe the composition of the atmosphere ([M/H] and C/O). In all the adequate fits to the data, these parameters are essentially unconstrained (e.g., [M/H] spans the whole prior range for acceptable DRIFT-PHOENIX models of PDS 70 b). There are a few edge cases that are excluded (e.g., C/O < 0.4 is excluded for adequate Exo-REM models of PDS 70 b), but we take such constraints with caution as atmosphere models can spuriously constrain C/O when there are other inaccuracies in the model (e.g., the plain Exo-REM fits to both planets have the smallest uncertainties on C/O, but the lowest B of all models).

For all three planetary atmosphere models, the implied masses based on the retrieved $\mathrm{log}(g)$ and radii generally favor masses > 10 M_Jup. However, our priors are biased to high masses as the model grids generally do not go down to a sufficiently low surface gravity for the ∼2 R_Jup effective radii we measured: a 1 M_Jup and 2 R_Jup planet has $\mathrm{log}(g)=2.8$, which is below the bounds of all our model grids. If we instead use the mass posterior for PDS 70 b from Section 3.2 as a prior on $\mathrm{log}(g)$, we obtained surface gravity values for PDS 70 b near the lower bound of all of the model grids, but none of the other atmospheric parameters changed significantly. As spectroscopic masses from surface gravity and radius have been shown to be unreliable for brown dwarf atmospheres of comparable temperatures (e.g., Zhang et al. 2020), we avoid overinterpreting the results on these protoplanetary photospheres.

It appears that the 1–5 μm data alone does not provide significant evidence CPD emission. Evidence for a second blackbody component by itself is marginal in our fits. For both PDS 70 b and c, the models with a second blackbody component for both the blackbody and DRIFT-PHOENIX models have smaller B than the plain models. Adding a second blackbody does improve the Bayes factor from the plain models for the BT-SETTL and Exo-REM models for both planets, but only the BT-SETTL model for PDS 70 b with a hot compact second component has an adequate Bayes factor.

The addition of extinction to the BT-SETTL and Exo-REM atmospheric models combined with the second blackbody component generally improved the Bayes factor significantly more than the addition of the second blackbody component alone. The BT-SETTL model with extinction and a second blackbody component has the highest Bayes factor for all of the models considered for PDS 70 b. However, all other extincted models with a second blackbody result in a lower Bayes factor than those with the addition of just ISM extinction alone.

Switching from the pure blackbodies to accreting CPD models from Zhu (2015) only decreases B, so there is no evidence that these models are better. The ${M}_{p}\dot{M}$ we derived are consistent with mass and mass accretion values from evolutionary models (Wang et al. 2020). At these low mass-accretion rates, the CPD SEDs look similar to blackbody emission (Zhu 2015), but may be less flexible than the blackbody model (e.g., the blackbody model prior range is flexible enough that it can negligibly alter the planetary SED in the observed spectral ranges if needed), resulting in a worse fit.

Evaluation of CPD emission would not be complete without considering emission at 855 μm from PDS 70 c, which argues for circumplanetary dust emission from PDS 70 c (Isella et al. 2019). This data point has a significant impact on the evidence for CPD emission given its large spectral lever arm. Unlike in the previous case of considering just 1–5 μm data, the DRIFT-PHOENIX model augmented with a cooler blackbody has the highest evidence by a factor of 10⁴, strongly ruling out models without a CPD (as seen in Table 6). Reassuringly, the parameters of the atmospheric model remain unchanged from Table 5, demonstrating that this 855 μm data point only probes CPD emission. Thus, our previous conclusions regarding the atmospheric properties of these protoplanets using soley the 1–5 μm data should hold. With this single data point constraining the CPD properties, the radius and temperature of the second blackbody component are dengenerate, but are consistent with the values found by Isella et al. (2019).

Thus, we found that there is some evidence for a second blackbody component for both planets when only considering the 1–5 μm data. The inclusion of the ALMA 855 μm detection for PDS 70 c definitively rejects models without a second blackbody component, demonstrating the need for observations at longer wavelengths to characterize the circumplanetary dust. These findings are consistent with those of Stolker et al. (2020), who found that the sole driver of the second blackbody model for PDS 70 b was their M-band photometry point, and Isella et al. (2019), who originally presented to detection of the CPD around PDS 70 c at 855 μm.

We interpret the favoring of planetary atmosphere models over featureless blackbody models to indicate that we indeed are seeing into the atmospheres of these protoplanets and that the accreting dust is not completely blocking all molecular signatures as was proposed by Wang et al. (2020). The effective radii of PDS 70 b from the best-fitting BT-SETTL model with extinction and an added blackbody component is between 1.8 and 2.2 R_Jup. Similarly for PDS 70 c, the plain DRIFT-PHOENIX model inferred radii between 1.7 and 2.3 R_Jup. Regardless, these effective radii are much smaller than has been found from previous works (Müller et al. 2018; Wang et al. 2020), and are even starting to be consistent with hot-start evolutionary models (Baraffe et al. 2003). Using the protoplanetary evolutionary models from Ginzburg & Chiang (2019a), for this age and luminosity, these effective radii are consistent with the lowest mean opacities for the atmospheres of the protoplanets ( < 2 × 10⁻² cm²/g using the values for PDS 70 b). Such low opacities have been predicted to occur due to the accretion of dust grains that have undergone grain growth or by coagulation of grains after accretion onto the planet, resulting in a distribution of grain sizes that favors more larger-sized grains than typical ISM distributions (Mordasini 2014; Piso et al. 2015).

The plain DRIFT-PHOENIX model without any modifications has some of the highest Bayes factors for both planets. Given that these models were originally designed to fit dusty, but older, brown dwarfs, we investigated why they have the most support from the data we have obtained so far. We found that this is due to the fact the DRIFT-PHOENIX models do not reproduce the L–T transition by having the clouds clear up at lower effective temperatures, but rather produce thicker clouds at T_eff < 1600 K that create a redder near-infrared spectrum (Witte et al. 2011). Given that our retrieved T_eff are all lower than 1600 K, all of the best-fit plain DRIFT-PHOENIX models should appear more dusty than typical substellar atmospheres due to the model creating a large dust cloud purely from atmospheric physics. However, the PDS 70 planets are known to be accreting (Haffert et al. 2019), with the accreting dust expected to shroud the atmosphere (Wang et al. 2020). Thus, we caution against interpreting these model-fitting results as indicating that PDS 70 b and c are just extremely cloudy substellar objects. Instead, it may be that this known deficiency in the DRIFT-PHOENIX models of producing extremely cloudy planets for T_eff < 1600 K may be emulating extinction from accreting dust to current measurement precision.

The plain DRIFT-PHOENIX models may not be so different from the extincted BT-SETTL and Exo-REM models that also have significant support from the data. These extincted models also seek to redden the planetary atmosphere to better match the overall 1–5 μm SED of both planets. The similar extinction amplitudes with dusty brown dwarfs that are not actively accreting may be a coincidence due to large uncertainties in the extinction characteristics. The sedimentation timescale of the dust in these protoplanet atmospheres should be on the order of 10 yr (Wang et al. 2020), so it is unlikely that we are seeing lingering dust from accretion in the field brown dwarfs. The formation of aerosols in the upper atmosphere through some undetermined process has been proposed to explain the dusty brown dwarf population instead (Hiranaka et al. 2016).

Overall, we interpret the fact that DRIFT-PHOENIX models and extincted BT-SETTL and Exo-REM models having the most support form the data to indicate that the planetary atmospheres are indeed significantly extincted by dust from the the planet formation process. The current spectral data does not well constrain the dust properties, so it is difficult to say how consistent the dust properties are compared to the 10⁻² cm²/g that is implied by evolutionary models. However, the A_V required for the extincted BT-SETTL and Exo-REM models are consistent with the constraint that A_Hα > 2 mag to be consistent with the non-detection of either planet with H β spectroscopy (Hashimoto et al. 2020). The drastically higher extinction of A_V ∼ 10 for PDS 70 c compared to PDS 70 b implies significantly more dust shrouding PDS 70 c. This could be inherent to the planet or circumplanetary environment since only PDS 70 c has a significant millimeter wavelength signal at its position, which implies a larger CPD than PDS 70 b (Isella et al. 2019). On the other hand, the inferred accretion rates are lower for PDS 70 c, which would imply that it harbors less dust around it (Haffert et al. 2019; Wang et al. 2020). Alternatively, the extinction could be enhanced due to additional extinction by the flared edge of the circumstellar disk (Keppler et al. 2018), which has been ignored in this and previous studies of PDS 70 c (Mesa et al. 2019; Wang et al. 2020). Better characterization of the 3D vertical structure of the circumstellar disk can help pinpoint if the extinction is due to the circumplanetary or circumstellar environment.

Although both protoplanets have been seen to emit in Hα (Wagner et al. 2018; Haffert et al. 2019), no discernible Brγ emission has been seen in previous observations (Christiaens et al. 2019b) or in our GRAVITY spectra. Here, we quantified some upper limits on the Brγ luminosity and its constraints on the accretion rates.

We can decompose the planetary flux density that we measure into continuum and line emission:

$$\begin{eqnarray}&&{F}_{\mathrm{planet}}(\lambda )={F}_{\mathrm{planet},\mathrm{cont}}(\lambda )+{F}_{\mathrm{planet},\mathrm{Br}\gamma }(\lambda ).\end{eqnarray} \tag{ 12 }$$

The continuum emission F_planet,cont is the broad K-band spectral shape that we have measured in our GRAVITY spectrum and analyzed in the previous SED fitting sections. Since the Hα emission does not appear to be resolved at 10 times higher spectral resolution than our GRAVITY observations, we expect that any Brγ emission would be unresolved with a spectral shape that is simply the line spread function of the GRAVITY instrument centered at the Brγ line. Assuming a Gaussian line spread function LSF(λ) with standard deviation σ_LSF , the flux density of the planet's line emission F_planet,Brγ can be written as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{planet},\mathrm{Br}\gamma }(\lambda ) & = & \displaystyle \frac{{f}_{\mathrm{planet},\mathrm{Br}\gamma }}{\sqrt{2\pi }{\sigma }_{{LSF}}}\exp \left(-\displaystyle \frac{{\left(\lambda -{\lambda }_{{\rm{Br}}\gamma }\right)}^{2}}{2{\sigma }_{{LSF}}^{2}}\right)\\ & \equiv & \ {f}_{\mathrm{planet},\mathrm{Br}\gamma }{{LSF}}_{{\rm{Br}}\gamma }.\end{array}\end{eqnarray} \tag{ 13 }$$

Here, f_planet,Brγ is the flux of the Brγ line integrated over all wavelengths, λ_Brγ = 2.166 μm, and σ_LSF = 0.0018 μm.

We estimated the integrated flux of the Brγ line by performing a matched filter of the continuum subtracted GRAVITY spectra with LSF_Brγ , the expected spectral shape of Brγ emission line. The matched filter across the GRAVITY spectral channels is written as:

$$\begin{eqnarray}&&{f}_{\mathrm{planet},\mathrm{Br}\gamma }=\displaystyle \frac{{\left({F}_{\mathrm{planet}}-{F}_{\mathrm{planet},\mathrm{cont}}\right)}^{T}{C}^{-1}({{LSF}}_{{\rm{Br}}\gamma })}{{\left({{LSF}}_{{\rm{Br}}\gamma }\right)}^{T}{C}^{-1}({{LSF}}_{{\rm{Br}}\gamma })},\end{eqnarray} \tag{ 14 }$$

where C is the covariance matrix of F_planet(λ) that we computed as part of our spectral extraction in Section 2.3. The error on f_planet,Brγ is then defined as

$$\begin{eqnarray}&&{\sigma }_{f,\mathrm{Br}\gamma }=\displaystyle \frac{1}{\sqrt{{\left({{LSF}}_{{\rm{Br}}\gamma }\right)}^{T}{C}^{-1}({{LSF}}_{{\rm{Br}}\gamma })}}.\end{eqnarray} \tag{ 15 }$$

We approximated the continuum emission F_planet,cont by applying a 41-channel median filter on the original planetary spectrum, F_planet. Subtracting this continuum emission from the original spectrum of each planet gave us an estimated wavelength-integrated Brγ flux of 0.1 ± 1.7 × 10⁻²⁰ W/m² for PDS 70 b, and 0.3 ± 1.3 × 10⁻²⁰ W/m² for PDS 70 c. Both are fully consistent with non-detections. These correspond to 3σ upper limits of 5.1 × 10⁻²⁰ W/m² for PDS 70 b and 4.0 × 10⁻²⁰ W/m² for PDS 70 c. The PDS 70 b upper limit is consistent with the 5σ upper limit of 8.3 × 10⁻²⁰ W/m² derived by Christiaens et al. (2019b).

We also followed the accretion luminosity and mass accretion analysis from Christiaens et al. (2019b), acknowledging that the relations are not calibrated to planetary mass companions forming in the disk, so unknown biases exist. Using the Calvet et al. (2004) relation between Brγ and total accretion luminosity for T Tauri stars, we find upper limits on the total accretion luminosity of <9.6 × 10⁻⁵ L_⊙ and <7.7 × 10⁻⁵ L_⊙ for PDS 70 b and PDS 70 c, respectively. Then, we can relate total accretion luminosity (L_acc) to mass accretion rate $\dot{M}$ using the equation $\dot{M}=1.25{L}_{\mathrm{acc}}{R}_{\mathrm{planet}}/{{GM}}_{\mathrm{planet}}$ that Christiaens et al. (2019b) used from Gullbring et al. (1998). For PDS 70 b, assuming a mass of 3 M_Jup and a radius of 2 R_Jup based on evolutionary model fits (Wang et al. 2020), we found an upper limit on the mass accretion rate of <2.9 × 10⁻⁷ M_Jup/ yr for PDS 70 b. For PDS 70 c, assuming a mass of 2 M_Jup and a radius of 2 R_Jup from the same evolutionary models, we found an upper limit on the mass accretion rate of <3.4 × 10⁻⁷ M_Jup/ yr. Both rates are consistent with most literature results (Wagner et al. 2018; Christiaens et al. 2019b; Haffert et al. 2019; Wang et al. 2020). This PDS 70 b upper limit from Brγ is incompatible with the lower limit derived by Hashimoto et al. (2020) based on the detection of Hα but non-detection of Hβ, and only marginally consistent with the range of mean accretion rates derived by Wang et al. (2020) from evolutionary models. The disagreements are not surprising given that we used relations that were calibrated to stars and not planets. Models of planetary accretion are needed to translate our Brγ upper limits to more realistic mass accretion limits.

With baselines up to 130 m, the VLTI at K-band wavelengths can achieve an angular resolution (λ/2B) of ∼2 mas (and can marginally resolve features at smaller angular scales given a sufficient signal-to-noise ratio in the data). The PDS 70 planetary system has a parallax of 8.8 mas (Gaia Collaboration et al. 2018), meaning that GRAVITY is able to spatially resolve scales down to at least 0.2 au (400 R_Jup) in projected separation. Here, we present our attempt to resolve the circumplanetary environments of the protoplanets.

Qualitatively, the marker of a resolved circumplanetary environment would be a drop in the coherent flux coming from the planet and its circumplanetary environment. We can look at two indicators to assess whether any emission was spatially resolved. The first indicator is the total coherent flux measured by GRAVITY, which probes spatial scales of a few mas, compared to the total flux within ∼40 mas measured by single-dish telescopes in the K band. If GRAVITY is spatially resolving the source, the coherent flux should be lower than the incoherent flux. With uncertainties of ∼10%–20% from SPHERE photometry, we did not see a significant drop for either PDS 70 b or c.

The second indicator would be a drop of the coherent flux as a function of increasing baseline. It would indicate that the longest baselines are spatially resolving structure that appears point like to the shorter baselines. Again, we did not see any significant drop of at least ∼10% in the flux going to longer baselines (see Figure 8). However, we used this technique to quantify upper limits on the spatial extent of the circumplanetary region. This is done by fitting synthetic visibility models to the measured visibilities.

Zoom In Zoom Out Reset image size

In this section, we considered the case where all the K-band emission is coming from a uniform circular disk. This could be true early on in the planet's formation, when the planet has just finished the runaway accretion phase and could have radii of ∼1000 R_Jup or ∼0.5 au (Ginzburg & Chiang 2019a). This is unlikely for these planets, given that the photometrically derived radii (see Section 4) are orders of magnitude smaller, and that the mass accretion rates measured with Hα and through evolutionary models have indicated that these planets are near the end of their formation process when their radii have contracted to near their final radii (Wagner et al. 2018; Haffert et al. 2019; Wang et al. 2020). Alternatively, this could approximate the case where the emission is dominated by the CPD (Zhu 2015; Szulágyi et al. 2019), although in Section 4.3, we favored emission from the planetary atmosphere. Regardless, we considered this a simple and limiting case for our ability to resolve the planet or its circumplanetary environment.

In the previous sections, we assumed that each planet was point like. This assumption resulted in Equation (3), where the contrast is equal to the planetary flux. We can remove this assumption by including a normalized visibility term J_planet(b, t, λ) proportional to the ratio between coherent flux, V_planet(b, t, λ), and total flux from the protoplanetary object, F_planet(λ):

$$\begin{eqnarray}&&{V}_{\mathrm{planet}}(b,t,\lambda )={J}_{\mathrm{planet}}(b,t,\lambda ){F}_{\mathrm{planet}}(\lambda ).\end{eqnarray} \tag{ 16 }$$

The normalized visibility can therefore be obtained from the ratio of the coherent flux between the planet and star, R(λ, b, t), using the equation:

$$\begin{eqnarray}&&{J}_{\mathrm{planet}}(b,t,\lambda )=R(\lambda ,b,t)\displaystyle \frac{{J}_{\mathrm{star}}(b,t,\lambda ){F}_{\mathrm{star}}(\lambda )}{\rho \ {F}_{\mathrm{planet}}(\lambda )},\end{eqnarray} \tag{ 17 }$$

where J_star(b, t, λ), F_star(λ), and ρ are established in Section 2.3. Additionally, it requires the knowledge of F_planet(λ). We assumed in this section that F_planet(λ) is a BT-SETTL model of temperature 1300 K. The exact model used is not important given our measurement precision. What plays a major role is the absolute scaling of the model: the higher it is, the lower the normalized visibility of the planet will be. In this work, we scaled the model to agree with the K-band magnitude of 16.5 ± 0.1 for PDS 70 b (Müller et al. 2018) and 17.7 ± 0.2 for PDS 70 c (Mesa et al. 2019) that were obtained through imaging.

The visibilities as a function of projected baseline length are plotted in Figure 8. To increase the signal-to-noise ratio, we averaged the visibility J_planet(b, t, λ) over all wavelengths and time: only one visibility value is obtained for each epoch and baseline. We used a uniform disk model, a Bessel function of first order (J₁) with R_planet as the exoplanet radius and $\sqrt{{u}^{2}+{v}^{2}}$ the projected length of the baseline in units of λ:

$$\begin{eqnarray}&&{J}_{\mathrm{planet}}(b,t,\lambda )=\displaystyle \frac{{{\rm{J}}}_{1}\left(2\pi {R}_{\mathrm{planet}}\sqrt{{u}^{2}+{v}^{2}}\right)}{\pi {R}_{\mathrm{planet}}\sqrt{{u}^{2}+{v}^{2}}}.\end{eqnarray} \tag{ 18 }$$

Fitting these models to both PDS 70 b and c, we derived a 1σ upper limit on the radii of the exoplanets of, respectively, 0.2 and 0.4 mas. Using a 3σ upper limit, we found a maximum planetary radius of 143 R_Jup for PDS 70 b and 285 R_Jup for PDS 70 c. These upper limits are consistent with our photometrically derived radii from Section 4 and protoplanet evolutionary models that are much smaller (Ginzburg & Chiang 2019a).

In reality, we expect an exoplanet of size ≈ 1–2 R_Jup with a circumplanetary disk (CPD) that is more extended. Thus, we constructed a two-component model where a small uniform disk represents the planet, and a larger uniform disk represents the CPD. For simplicity, we fix the radius of the planet to be 2 R_Jup, which is consistent with our photometrically derived radii from Section 4. We also assumed that the CPD contributed 10% of the K-band flux. This is near the upper bound of what is predicted by our models that include a second blackbody component in the SED fits. Alternatively, the CPD could be bright due to scattered starlight, since scattered starlight dominates the K-band emission of the outer circumstellar disk. Regardless, a CPD brightness much lower than this would be completely undetectable given our measurement errors, so our analysis here is only suitable for constraining a bright CPD. Thus, we created a model that included the contribution of both the planet and the CPD:

$$\begin{eqnarray}&&{J}_{\mathrm{planet}}(b)=0.9\displaystyle \frac{2{{\rm{J}}}_{1}({x}_{\mathrm{planet}})}{{x}_{\mathrm{planet}}}+0.1\displaystyle \frac{2{{\rm{J}}}_{1}({x}_{\mathrm{CPD}})}{{x}_{\mathrm{CPD}}}\end{eqnarray} \tag{ 19 }$$

where the only free parameter is the radius of the CPD (R_CPD):

$$\begin{eqnarray}&&{x}_{\mathrm{planet}}=2\pi {R}_{\mathrm{planet}}\sqrt{{u}^{2}+{v}^{2}}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{x}_{\mathrm{CPD}}=2\pi {R}_{\mathrm{CPD}}\sqrt{u{{\prime} }^{2}+v{{\prime} }^{2}}\end{eqnarray} \tag{ 21 }$$

Note that the projected baseline length, $\sqrt{u{{\prime} }^{2}+v{{\prime} }^{2}}$, which corresponds to a position in the frequency plane, is modified to account for the inclination of the CPD. The CPD will look shorter in the direction perpendicular to its position angle due to viewing geometry. We assumed the CPD has the same orientation and inclination as the circumstellar disk. We took values derived from ALMA millimeter-wavelength measurements: the inclination of the disk is i = 517 and the position angle is α = 1567 (Keppler et al. 2019). Hence, the u–v plane is shifted to the new coordinates:

$$\begin{eqnarray}&&u^{\prime} =\cos i\ (u\cos \alpha -v\sin \alpha )\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&v^{\prime} =u\sin \alpha +v\cos \alpha .\end{eqnarray} \tag{ 23 }$$

For PDS 70 b, we found a 1σ upper limit of R_CPD <1 mas, which we translated to a 3σ upper limit of 0.3 au. Although we did not detect any signature of a CPD, we have demonstrated the first observations that can resolve scales of <1 au around a protoplanet.

For PDS 70 c, the fit is dominated by the error caused by the normalization of the visibilities (the visibility of 1 is assumed to correspond to a coherent flux of Kmag = 17.7 ± 0.2). If we neglected this uncertainty, the fit converged to a large CPD of radius R_CPD > 10 mas. Because the extent of such CPD seems unrealistic (see below), we are inclined to think that the drop in visibility is purely instrumental, due to the uncertainty in the SPHERE K-band magnitude.

The size of a CPD is governed by either the Hill or Bondi radius, whichever is smaller. The Hill radius (R_H) for an eccentric planet is defined as

$$\begin{eqnarray}&&{R}_{{\rm{H}}}=a(1-e)\sqrt[3]{{M}_{{\rm{p}}}/(3{M}_{* })},\end{eqnarray} \tag{ 24 }$$

where a is the semimajor axis, e is the eccentricity, M_p is the planet mass, and M_* is the stellar mass. For PDS 70 b, if we use a planet mass of 3 M_Jup based on evolutionary models (Wang et al. 2020), a stellar mass of 0.9 M_⊙, a semimajor axis of 20 au, and an eccentricity of 0.2 from our orbit fit (Section 3), then we derive R_H = 1.6 au.

The Bondi radius for an isothermal ideal gas is

$$\begin{eqnarray}&&{R}_{{\rm{B}}}=\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{c}_{s}^{2}}=\displaystyle \frac{{{GM}}_{{\rm{p}}}\mu }{{k}_{{\rm{B}}}T},\end{eqnarray} \tag{ 25 }$$

where μ is the mean molecular weight, G is the gravitational constant, k_B is the Boltzmann constant, M_p is the planet mass, T is the gas temperature, and c_s is the sound speed. Assuming the gas is a blackbody in thermal equilibrium with the star, the gas temperature is 50 K. Assuming the same planet mass as before and μ = 2 amu for molecular hydrogen, we found a Bondi radius of 8 au. Our upper limit of R_CPD <0.3 au thus corresponds to 0.2 R_H or 0.04 R_B.

With R_H <R_B, we are likely looking at a planet that is slightly above the thermal mass, the mass above which accretion starts transitioning into a 2D process from a 3D one (Ginzburg & Chiang 2019b; Rosenthal et al. 2020). We compared our upper limit on PDS 70 b's CPD to the CPD models of Szulágyi et al. (2017) and Fung et al. (2019) that explore planet masses above the thermal mass regime. Szulágyi et al. (2017) reported results in R_H whereas Fung et al. (2019) provided disk sizes in R_B. For planets with comparable Hill radii, Szulágyi et al. (2017) found CPDs with spatial scales of up to 0.1 R_H, which agrees well with our upper limit of 0.2 R_H. Fung et al. (2019) found that the CPD are about 0.1 R_B for planets just above the thermal mass. We found a upper limit in units of Bondi radii that is approximately three times smaller. We note though that we assumed a very bright CPD (10% of the total K-band flux), and that a fainter CPD would escape detection even if it was more extended than 0.04 R_B. Our upper limit of 0.3 au is also consistent with the upper limit of 0.1 au for a CPD around PDS 70 b derived using SED fitting by Stolker et al. (2020), accounting for a non-detection at 855 μm.

For PDS 70 c, using a planet mass of 2 M_Jup (Wang et al. 2020), a semimajor axis of 30 au, and an eccentricity of 0, we found R_H = 2.7 au. Assuming a gas temperature of 40 K at 30 au, we find R_B = 10 au. The limits from Szulágyi et al. (2017) and Fung et al. (2019) imply CPD sizes of 0.5 au (5 mas) and 1 au (9 mas) respectively. Both are smaller than the 10 mas mentioned in the previous suction that would explain the visibility normalization offset, arguing for photometric calibration to be responsible for the offset.