Evidence for a Circumplanetary Disk around Protoplanet PDS 70 b

Discovery of the Galilean moons of Jupiter (Galilei 1610) led to the overturning of Geocentric cosmology and Galileo's subsequent denouncement, trial, and house arrest. Their hypothesized origin is a circumplanetary disk (CPD) of gas and dust (e.g., Lunine & Stevenson 1982; Lissauer 1993). Analytic work (Pollack et al. 1979; Canup & Ward 2002; Papaloizou & Nelson 2005) and increasingly sophisticated numerical simulations (e.g., Lubow et al. 1999; Ayliffe & Bate 2009; Gressel et al. 2013; Szulágyi et al. 2017) have consolidated this hypothesis.

Observational evidence has so far remained elusive. Searches using the mm-continuum have produced only nondetections (e.g., Isella et al. 2014; Wolff et al. 2017; Wu et al. 2017). High-contrast infrared (IR) observations have been suggested instead to capture the thermal excess from the CPD (e.g., Eisner 2015; Montesinos et al. 2015; Zhu 2015). The power of high-contrast IR spectroscopy to characterize young substellar companions has been demonstrated over the last few years (e.g., Allers & Liu 2013; Bonnefoy et al. 2014; Delorme et al. 2017). Observed spectra are fitted to either synthetic atmospheric models or observed template spectra that enable the estimation of effective temperature, surface gravity, and radius of the planet, which can then be used to estimate its mass and age. Constraints on clouds or haze in the atmosphere can also be obtained (e.g., Madhusudhan et al. 2011, hereafter M11; see Madhusudhan 2019 for a recent review).

We present evidence for a CPD around the recently discovered protoplanet PDS 70 b. PDS 70 is a young K7-type star surrounded by a pre-transitional disk with a large annular gap (Dong et al. 2012; Hashimoto et al. 2012). The protoplanet was first detected in this gap using the Very Large Telescope (VLT)/SPHERE (Keppler et al. 2018; Müller et al. 2018, hereafter M18). Spectral characterization of the planet was presented in M18, suggesting an effective temperature of ∼1000–1500 K, surface gravity ≲3.5 dex, and mass ≤17M_J. We present a new spectral characterization of PDS 70 b, including both the measurements presented in M18 and a new spectrum in the K-band obtained with VLT/SINFONI (Christiaens et al. 2019, hereafter C19). The spectrum shows excess emission at λ ≳ 2.3 μm that is inconsistent with naked model atmospheres of young planets, indicating the presence of circumplanetary material.

M18 presented a YJH spectrum and multi-epoch broadband photometric measurements, acquired with SPHERE/IFS (YJH), SPHERE/IRDIS (H1/H2 and K1/K2), NICI (L'), and NACO (L'). We gathered all measurements quoted in M18, but considered updated values for their L'-band flux estimates. M18 fitted the spectral energy distribution (SED) of the star to estimate its L' flux, without including possible excess disk emission. That value was then multiplied by the contrast of the companion to infer its flux. However, most of excess IR emission at the L'-band compared to the star likely arises from hot dust in the inner disk that would be unresolved from the star, and should hence be included. By contrast, the contribution from resolved scattered light from the disk at the L'-band is negligible (see images in Keppler et al. 2018). Given the importance of the thermal IR flux to the hypothesis of a CPD around PDS 70 b, we re-estimated the L' flux of the companion considering (1) the L' flux of the star + unresolved inner disk by interpolating the photometric measurements in the W1 (3.35 μm) and W2 (4.60 μm) filters of the Wide-field Infrared Survey Explorer (Wright et al. 2010), and (2) the same values for the contrast of the companion as in M18. The new L' estimates of the companion are: (8.55 ± 3.38) × 10⁻¹⁷ W m⁻² μm⁻¹ and (6.70 ± 4.05) × 10⁻¹⁷ W m⁻² μm for the NICI and NACO data, respectively.

In a recent paper (C19), we inferred the contrast of the protoplanet with respect to the star as a function of wavelength, c(λ), in the K-band at unprecedented spectral resolution (R ∼ 100 after spectral binning) using VLT/SINFONI. We employed two different methods for extracting c(λ), both leading to flux estimates consistent with each other at all wavelengths. Here, we consider only the spectrum inferred with ANDROMEDA (Cantalloube et al. 2015, F. Cantalloube et al. 2019, in preparation) because it has smaller uncertainties and higher spectral resolution, and avoids the risk of contamination by extended (resolved) disk signals.

Despite the higher quality of the ANDROMEDA c(λ), some spectral channels contain outliers. In order to minimize the risk of bias, we first removed spectral channels with a detection below 3σ and lying in strong telluric lines, then used a Savitzky–Golay filter of order 3 with a 81-channel window to smooth the c(λ) curve before binning it by a factor of 20 (Savitzky & Golay 1964). We obtained the final K-band spectrum of the protoplanet by multiplying c(λ) with the calibrated spectrum of the star measured with the SpeX spectrograph (Long et al. 2018), after resampling the latter at the spectral resolution of the binned SINFONI spectral channels.

In total, our SED has 86 data points; 49 from M18 and 37 obtained with SINFONI (C19). Measurements span 6 yr, extending from 2012 March to 2018 February. Flux estimates at overlapping wavelengths are all consistent with each other except one. Namely, our SINFONI flux estimates (2014 May epoch) are slightly higher than the 2018 February epoch SPHERE measurement in the K2 filter (2.25 μm).

We first attempted to fit the observed SED with synthetic spectra modeling pure atmospheric emission (Section 3.1). We considered a single value of extinction for all epochs, treated as a free parameter in the spectral fit (referred to as Type I models hereafter). Given the discrepancy between the SPHERE K2 and SINFONI K-band measurements, we also considered a fit with variable amount of extinction for different epochs (Type II models). Finally, we also examined models consisting of combined emission from an atmosphere and a CPD with a variable accretion rate (Type III models; Section 3.2). In order to minimize the number of free parameters, we considered only two possible values of extinction (for Type II models) and two accretion rates (for Type III models); one value for the SINFONI (2014 May), SPHERE (2016 May), and NICI (2012 March) epochs, and the other for all other epochs. This division was chosen because (1) the SPHERE IFS (YJH)+IRDIS (K1/K2) were obtained simultaneously on 2018 February, and the IRDIS H1/H2 points of 2015 May are consistent with the IFS spectrum, and (2) the SINFONI 2014 May data are brighter than the IRDIS K2 point of 2018 February, while the IRDIS K2 point of 2016 May is in better agreement with the SINFONI data. For the NICI and NACO points, we arbitrarily assigned them to the first and second group, respectively.

For all model types, we minimized the following goodness-of-fit indicator:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}{\omega }_{i}{\left[\displaystyle \frac{{F}_{\mathrm{obs}}({\lambda }_{i})-{F}_{\mathrm{model}}({\lambda }_{i})}{{\sigma }_{i}}\right]}^{2}\end{eqnarray} \tag{ 1 }$$

where σ_i is the uncertainty in the flux measurement F_obs(λ_i) at wavelength λ_i, and weights ω_i are defined for photometric and spectroscopic observations following a similar strategy as in Ballering et al. (2013) and Olofsson et al. (2016). Weights are proportional to the FWHM of the filters used (for broadband photometric measurements), or the spectral resolution (for SPHERE/IFS and SINFONI data). The sum of all weights is normalized to the total number of points. We define a reduced goodness-of-fit indicator ${\chi }_{r}^{2}$ as χ² divided by the respective number of degrees of freedom for each type of model.

3.1. Atmospheric Models

We considered two grids of synthetic spectra: BT-SETTL models (Allard et al. 2012; Baraffe et al. 2015), and the grid of atmospheric models presented in M11. These models treat dust and clouds differently. BT-SETTL models account for dust formation using a parameter-free cloud model. They consider cloud microphysics—the particle size of each species is calculated self-consistently based on condensation and sedimentation mixing timescales. Free parameters are the effective temperature, varied between 1200 and 1900 K (steps of 100 K), and the surface gravity, $\mathrm{log}(g)$, explored between 3.0 and 5.0 (steps of 0.5 dex). We assumed Solar metallicity.

M11 considered a wide grid of cloud models, with different geometrical and optical thickness, particle size, and metallicities. The M11 models do not include microphysics. They consider different cloud spatial structure and particle sizes; labeled A, AE, AEE, or E, based on the rapidity with which clouds are cut off at their upper end. Several modal particle sizes are considered, including 1, 60, and 100 μm. The grid also includes cloud-free models (NC), with both equilibrium and non-equilibrium chemistry. In the latter case, two additional free parameters arise: the eddy diffusion coefficient K_zz taking possible values of 10², 10⁴, or 10⁶ cm² s⁻¹, and the sedimentation parameter f_sed (defined as in Ackerman & Marley 2001). M11 varied these parameters on a grid of effective temperature and surface gravity ranging from 600 to 1800 K (steps of 100 K) and from 3.5 to 5.0 (steps of 0.5 dex), respectively.

For both BT-SETTL and M11 models, we treated the planet radius as a free parameter to scale the total emergent flux. We explored values between 0.1 and 5.0 R_J in steps of 0.1 R_J. We also considered dust extinction as an additional free parameter, with allowed values between A_V = 0 and 10.0 mag (steps of 0.2 mag). For Type II models, we explored both minimum and maximum extinction values within this grid. We considered the extinction curve for interstellar dust (Draine 1989). Considering other dust species (e.g., typical species found in the atmosphere of brown dwarfs) would increase the number of free parameters, and is not expected to give qualitatively different slopes after dereddening (see e.g., Marocco et al. 2014). We assumed a distance of 113 pc (Gaia Collaboration et al. 2018).

Figures 1 and 2 compare our best-fit models from the BT-SETTL and M11 grids, respectively, with the SED of PDS 70 b. Table 1 gives the explored parameter ranges, best-fit parameter values, and corresponding ${\chi }_{r}^{2}$, for each model. Best-fit Type I BT-SETTL and M11 models reproduce most of the observed SED but are a poor match to the red end of the SINFONI spectrum, with a >2σ discrepancy for most data points at wavelengths ≳2.3 μm. They lead to reduced goodness-of-fit indicators ${\chi }_{r}^{2}\sim 1.01$ and 1.20.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Best-fit Parameters for PDS 70 b

I. Planet Alone			II. Planet Alone (Variable Extinction)			III. Planet and CPD
Parameter	Range	Best Fit	Parameter	Range	Best Fit	Parameter	Range	Best Fit
BT-SETTL Atmospheric Models
${T}_{\mathrm{eff}}$ [K]	1200–1900	1500	Teff [K]	1200–1900	1500	Teff [K]	1200–1900	1500
$\mathrm{log}(g)$	3.0–5.0	3.0	$\mathrm{log}(g)$	3.0–5.0	3.0	$\mathrm{log}(g)$	3.0–5.0	4.0
Rb [RJ]	0.1–5.0	2.2	Rb [RJ]	0.1–5.0	2.1	Rb [RJ]	0.1–2.0	1.6
AV [mag]	0.0–10.0	4.34	${A}_{V,1}$ [mag]	0.0–10.0	0.60	AV [mag]	0.0–10.0	6.40
			${A}_{V,2}$ [mag]	0.0–10.0	4.00	${M}_{{\rm{b}}}{\dot{M}}_{{\rm{b}},1}$ $[{M}_{{\rm{J}}}^{2}$ yr−1]	10−7–10−6	10−6.4
						${M}_{{\rm{b}}}{\dot{M}}_{{\rm{b}},2}$ $[{M}_{{\rm{J}}}^{2}$ yr−1]	10−7–10−6	10−6.8
Mba [MJ]	⋯	1.9	Mba [MJ]	⋯	1.7	Mba [MJ]	⋯	9.9
						${\dot{M}}_{{\rm{b}}}$ a[MJ yr−1]	⋯	10−7.8–10−7.4
		${\chi }_{r}^{2}\sim 1.01$			${\chi }_{r}^{2}\sim 0.52$			${\chi }_{r}^{2}\sim 0.41$
M11 Atmospheric Models
Teff [K]	600–1800	1100	Teff [K]	600–1800	1200	Teff [K]	600–1800	1600
$\mathrm{log}(g)$	3.5–5.5	4.0	$\mathrm{log}(g)$	3.5–5.5	4.0	$\mathrm{log}(g)$	3.5–5.5	4.0
Rb [RJ]	0.1–5.0	3.3	Rb [RJ]	0.1–5.0	2.8	Rb [RJ]	0.1–2.0	1.6
c	[NC, E, A, ...]b	A	c	[NC, E, A, ...]b	A	c	[NC, E, A, ...]b	A60b
fsed	[eq., 0, 1, 2]	eq.	fsed	[eq., 0, 1, 2]	eq.	fsed	[eq., 0, 1, 2]	eq.
K [cm2 s−1]	[eq., 102–106]	eq.	K [cm2 s−1]	[eq., 102–106]	eq.	K [cm2 s−1]	[eq., 102–106]	eq.
AV [mag]	0.0–10.0	3.00	${A}_{V,1}$ [mag]	0.0–10.0	0.80	AV [mag]	0.0–10.0	8.72
			${A}_{V,2}$ [mag]	0.0–1.2	4.00	${M}_{{\rm{b}}}{\dot{M}}_{{\rm{b}},1}$ $[{M}_{{\rm{J}}}^{2}$ yr−1]	10−7–10−6	10−6.3
						${M}_{{\rm{b}}}{\dot{M}}_{{\rm{b}},2}$ [MJ2 yr−1]	10−7–10−6	10−6.7
Mba [MJ]	⋯	42.0	Mba [MJ]	⋯	30.2	Mba [MJ]	⋯	9.9
						${\dot{M}}_{{\rm{b}}}$ a $[{M}_{{\rm{J}}}$ yr−1]	⋯	10−7.7–10−7.3
		${\chi }_{r}^{2}\sim 1.20$			${\chi }_{r}^{2}\sim 0.70$			${\chi }_{r}^{2}\sim 0.44$

Notes.

^aM_b is inferred from the best-fit $\mathrm{log}(g)$ and R_b values, and ${\dot{M}}_{{\rm{b}}}$ is inferred from M_b and the best-fit ${M}_{{\rm{b}}}{\dot{M}}_{{\rm{b}}}$. ^bSee Section 3.1 and M11. A60 refers to model cloud A (thickest) with a modal particle size of 60 μm.

Download table as:  ASCIITypeset image

Allowing for variable extinction (Type II) yields best-fit models (shaded cyan areas in Figures 1 and 2) that are in better agreement with flux estimates longward of ∼2.3 μm (${\chi }_{r}^{2}\sim 0.52$ and 0.70). Upper and lower edges of the shaded areas correspond to the minimum and maximum extinction values A_V,1 and A_V,2, in Table 1) of our best-fit model, with A_V,1 better accounting for the 2014 May SINFONI, 2016 May SPHERE and 2012 March NICI data points and A_V,2 accounting for data points at all other epochs. The difference in extinction is A_V ≳ 3 mag for both BT-SETTL and M11 models.

For both Type I and II models, the best-fit effective temperatures (1100–1500 K), surface gravity ($\mathrm{log}(g)\sim 3.0\mbox{--}4.0$), planet radius (2.1–3.3R_J), and hence mass (1.7–42.0 M_J) are in approximate agreement with the previous estimates made in M18. Our mass estimates are uncertain because of the large steps in $\mathrm{log}(g)$ (0.5 dex) in our model grids, and hence do not rule out a brown dwarf. Our best-fit M11 models correspond to the thickest cloud models (labeled A); extending to the top of the atmosphere.

3.2. Combined Atmospheric+CPD Models

We explored the hypothesis of variability for PDS 70 b because of the absence of atmospheric models red enough to account for both the 2018 February SPHERE K2 measurement and the points at ≳2.3 μm in our 2014 May SINFONI spectrum. This is further supported by the disagreement (albeit slight) between the SPHERE and SINFONI measurements at ∼2.25 μm. Because the variability of classical T-Tauri stars is thought to be related to either variable amounts of extinction from intervening circumstellar dust (e.g., Bozhinova et al. 2016; Scholz et al. 2019) or irregular accretion (e.g., Bouvier et al. 2004; Rigon et al. 2017), we tested similar hypotheses in our type II and III models, respectively. Accretion variability has also been predicted in magnetohydrodynamics simulations of forming planets (e.g., Gressel et al. 2013), further justifying type III models.

The best fit to the SED of PDS 70 b is obtained with an atmosphere+CPD model. Nonetheless, our caveats are as follows.

• 1.  
  We considered a limited range of atmospheric models. Atmospheric models have been proposed in recent years with different levels of complexity, including, for example, microphysics, non-equilibrium chemistry, or clouds/hazes (see Madhusudhan 2019, and references therein). We used the two most complete publicly available synthetic atmospheric model grids, which successfully reproduce the spectrum of adolescent giant exoplanets such as Beta Pic b or HR 8799 b, c, d, and e (M11; Bonnefoy et al. 2013; Currie et al. 2014). Given that our best-fit type I (purely atmospheric) BT-SETTL and M11 models are similar to the reddest atmospheric models (around ∼2.3 μm) presented in M18, we do not expect our conclusion of an excess at ≳2.3 μm to change using different grids of atmospheric models.
• 2.  
  The fit is not perfect. Although the best-fit atmosphere+CPD model best reproduces the excess at the end of the K-band, it does not perfectly reproduce the observed slope around 2.3 μm. For both the BT-SETTL and M11 type III models, most photometric points at wavelengths that are shorter than 2.3 μm lie below the model (albeit all within 2σ), while most points longward of 2.3 μm are slightly brighter than the model (but within 2σ also).
• 3.  
  The assumptions behind our CPD models may be incorrect. Incorrect assumptions may explain why our CPD models do not reproduce perfectly the observed steep slope. We fixed the inner truncation radius to 2R_J. As shown in Zhu (2015), different inner truncation radii can lead to different predicted CPD spectra for a given mass accretion rate. Furthermore, neither the models in Zhu (2015) nor in Eisner (2015) take into account radiative feedback from the protoplanet itself. The best-fit effective temperature and radius found for PDS 70 b suggest a protoplanet luminosity L_p ∼ 10⁻⁴ L_⊙. Montesinos et al. (2015) showed that the effect of such bright protoplanet would be to further increase the IR excess of the CPD with a steeper spectral slope, hence possibly improving the match with the K-band spectrum of PDS 70 b. New dedicated simulations are required to verify this hypothesis.

However, several lines of evidence support the hypothesis of circumplanetary material.

• 1.  
  Our best-fit accretion rates agree with observations. Assuming a similar relationship between Hα luminosity and mass accretion rate as T-Tauri stars, Wagner et al. (2018) estimated the PDS 70 b accretion rate to be$\sim {10}^{-7.8}{M}_{{\rm{J}}}^{2}{R}_{{\rm{J}}}^{-1}\,{\mathrm{yr}}^{-1}$ in the absence of extinction. However, the protoplanet is likely to be embedded within the circumplanetary material from which it feeds. The observed accretion rate would then be $\sim {10}^{-6.7}{M}_{{\rm{J}}}^{2}{R}_{{\rm{J}}}^{-1}\,{\mathrm{yr}}^{-1}$ for A_V ≈ 3 mag. Assuming a similar relationship for Brγ emission, C19 also constrained the accretion rate to be $\lt {10}^{-6.2}{M}_{{\rm{J}}}^{2}{R}_{{\rm{J}}}^{-1}\,{\mathrm{yr}}^{-1}$, considering negligible extinction at the K-band. Our best-fit models suggest strong extinction toward the protoplanet (A_V > 6 mag), which is consistent with the presence of a CPD. Our estimates of mass (M_b ∼ 10M_J), accretion rates (∼10^−7.5 M_J yr⁻¹), extinction (≳6 mag), and radii R_b (∼1.6R_J) appear roughly compatible with both observational constraints. Monitoring of the Hα luminosity would confirm the variability of the accretion rate.
• 2.  
  The photometric variability is only observed at relatively long NIR wavelength (∼2.2 μm). The best-fit models involving variable extinction lead to only slightly worse fits to the data than the atmosphere+CPD best-fit models. However, for the former models the amplitude of the variability is larger at short than at long NIR wavelengths, which is not observed despite multiple epoch observations at wavelengths shorter than 1.7 μm. Even if the variability was due to extinction, given the radial separation of the protoplanet (∼20 au), both the amplitude (ΔA_V ≳ 3 mag) and timescale of extinction variability (less than several years) suggests it would also be caused by circumplanetary dust.
• 3.  
  Tentative excess with respect to the atmosphere+CPD models at 2.29–2.35 μm might suggest CO bandhead emission. For young stellar objects, CO bandhead emission (Δv = 2; first transitions at 2.294, 2.323 and 2.352 μm) is an indicator of disk presence (Geballe & Persson 1987; Davis et al. 2011). CPD models in Ayliffe & Bate (2009) and Szulágyi (2017) predict temperatures up to several thousand K, which might also produce CO bandhead emission.
• 4.  
  Presence of a spiral arm. Our conclusion regarding the presence of circumplanetary material around PDS 70 b is consistent with recent images obtained with VLT/SINFONI, suggesting the presence of an outer spiral arm likely feeding the CPD (C19).

The SED of PDS 70 b is best fit by models that include a CPD. Atmospheric models alone are not able to account for the observed flux at wavelengths ≳ 2.3  μm. We infer an accretion rate of 10^−7.8–10^−7.3 M_J yr⁻¹ for a ∼10 M_J planet with significant extinction, which is consistent with prior observations. Simultaneous follow-up observations of PDS 70 b with wide spectral coverage in the NIR, including the Hα line, should confirm the scenario of variable accretion through a CPD.

We acknowledge funding from the Australian Research Council via DP180104235, FT130100034 and FT170100040. V.C. thanks Andre Müller for sharing the SPHERE spectrum of PDS 70 b. M.M. acknowledges financial support from the Chinese Academy of Sciences (CAS) through CASSACA-CONICYT Postdoctoral Fellowship (Chile).

Facilities: VLT(SINFONI - , SPHERE) - .

Software: ANDROMEDA (Cantalloube et al. 2015), VIP (Gomez Gonzalez et al. 2017), SpeX (http://www.browndwarfs.org/spexprism).