Evidence for a circumplanetary disc around protoplanet PDS 70 b
We present the first observational evidence for a circumplanetary disc around the protoplanet PDS 70 b, based on a new spectrum in the band acquired with VLT/SINFONI. We tested three hypotheses to explain the spectrum: Atmospheric emission from the planet with either (1) a single value of extinction or (2) variable extinction, and (3) a combined atmospheric and circumplanetary disc model. Goodness-of-fit indicators favour the third option, suggesting circumplanetary material contributing excess thermal emission — most prominent at m. Inferred accretion rates (– yr) are compatible with observational constraints based on the H and Br lines. For the planet, we derive an effective temperature of 1500–1600 K, surface gravity , radius , mass , and possible thick clouds. Models with variable extinction lead to slightly worse fits. However, the amplitude (mag) and timescale of variation ( years) required for the extinction would also suggest circumplanetary material.
0000-0002-0101-8814]Valentin Christiaens \move@AU\move@AF\@affiliationSchool of Physics and Astronomy, Monash University, VIC 3800, Australia \collaboration
0000-0002-3968-3780]Faustine Cantalloube \move@AU\move@AF\@affiliationMPIA, Heidelberg, Germany \collaboration
0000-0002-0433-9840]Simon Casassus \move@AU\move@AF\@affiliationDepartamento de Astronomia, Universidad de Chile, Santiago, Chile \collaboration
0000-0002-4716-4235]Daniel J. Price \move@AU\move@AF\@affiliationSchool of Physics and Astronomy, Monash University, VIC 3800, Australia \collaboration
0000-0002-4006-6237]Olivier Absil \@AF@joinF.R.S.-FNRS Research Associate \move@AU\move@AF\@affiliationSTAR institute, Université de Liège, Liège, Belgium \collaboration
0000-0001-5907-5179]Christophe Pinte \move@AU\move@AF\@affiliationSchool of Physics and Astronomy, Monash University, VIC 3800, Australia \collaboration
0000-0001-8627-0404]Julien Girard \move@AU\move@AF\@affiliationSpace Telescope Science Institute, Baltimore, USA \collaboration
0000-0003-0363-1492]Matias Montesinos \move@AU\move@AF\@affiliationInstituto de Física y Astronomía, Universidad de Valparaíso, Chile \move@AU\move@AF\@affiliationCASSACA, National Astronomical Observatories, Beijing, China \move@AU\move@AF\@affiliationNúcleo Milenio de Formación Planetaria, Chile \collaboration
Madhusudhan2011M11 \defcitealiasMuller2018M18 \defcitealiasChristiaens2019aC19
Discovery of the Galilean moons of Jupiter (Galilei1610) led to the overthrow of Geocentric cosmology, and Galileo’s subsequent denouncement, trial and house arrest. Their hypothesised origin is a circumplanetary disc (CPD) of gas and dust (e.g. Lunine1982; Lissauer1993). Analytic work (Pollack1979; Canup2002; Papaloizou2005) and increasingly sophisticated numerical simulations (e.g. Lubow1999; Ayliffe2009; Gressel2013; Szulagyi2017) have consolidated this hypothesis.
Observational evidence has so far remained elusive. Searches using mm-continuum have produced only non-detections (e.g. Isella2014; Wu2017; Wolff2017). High-contrast infrared (IR) observations have been suggested instead to capture the thermal excess from the CPD (e.g. Zhu2015b; Eisner2015; Montesinos2015). The power of high-contrast IR spectroscopy to characterize young substellar companions has been demonstrated over the last few years (e.g. Allers2013; Bonnefoy2014; Delorme2017). Observed spectra are fitted to either synthetic atmospheric models or observed template spectra which 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. Madhusudhan2011, hereafter Madhusudhan2011; see Madhusudhan2019 for a recent review).
We present evidence for a circumplanetary disc around the recently discovered protoplanet PDS 70 b. PDS 70 is a young K7-type star surrounded by a pre-transitional disc with a large annular gap (Hashimoto2012; Dong2012). The protoplanet was first detected in this gap using VLT/SPHERE (Keppler2018; Muller2018, hereafter Muller2018). Spectral characterisation of the planet was presented in Muller2018, suggesting an effective temperature of 1000-1500 K, surface gravity dex and mass . We present a new spectral characterization of PDS 70 b, including both the measurements presented in Muller2018 and a new spectrum in band obtained with VLT/SINFONI (Christiaens2019a, hereafter Christiaens2019a). The spectrum shows excess emission at m inconsistent with naked model atmospheres of young planets, indicating the presence of circumplanetary material.
Muller2018 presented a spectrum and multi-epoch broadband photometric measurements, acquired with SPHERE/IFS (), SPHERE/IRDIS (/ and /), NICI () and NACO (). We gathered all measurements quoted in Muller2018, but considered updated values for their -band flux estimates. Muller2018 fitted the SED of the star to estimate its flux, without including possible excess disc emission. That value was then multiplied by the contrast of the companion to infer its flux. However, most of excess IR emission compared to the star likely arises from hot dust in the inner disc that would be unresolved from the star at band, and should hence be included. By contrast, the contribution from resolved scattered light from the disc at band is negligible (see images in Keppler2018). Given the importance of the thermal IR flux to the hypothesis of a CPD around PDS 70 b, we re-estimated the flux of the companion considering (1) the flux of the star + unresolved inner disc by interpolating the photometric measurements in the W1 (m) and W2 (m) filters of WISE (Wright2010), and (2) the same values for the contrast of the companion as in Muller2018. The new estimates of the companion are: W m m and W mm for the NICI and NACO data, respectively.
In a recent paper (Christiaens2019a), we inferred the contrast of the protoplanet with respect to the star as a function of wavelength, , in the band at unprecedented spectral resolution ( after spectral binning) using VLT/SINFONI. We employed two different methods for extracting , both leading to flux estimates consistent with each other at all wavelengths. Here, we consider only the spectrum inferred with ANDROMEDA (Cantalloube2015, Cantalloube et al. 2019 in prep.) since it has smaller uncertainties and higher spectral resolution, and avoids the risk of contamination by extended (resolved) disc signals.
Despite the higher quality of the ANDROMEDA , 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 curve before binning it by a factor of 20 (Savitzky1964). We obtained the final -band spectrum of the protoplanet by multiplying with the calibrated spectrum of the star measured with the SpeX spectrograph (Long2018), after resampling the latter at the spectral resolution of the binned SINFONI spectral channels.
In total, our SED has 86 data points; 49 from Muller2018 and 37 obtained with SINFONI (Christiaens2019a). Measurements span 6 years, extending from 2012/03 to 2018/02. Flux estimates at overlapping wavelengths are all consistent with each other except one. Namely, our SINFONI flux estimates (2014/05 epoch) are slightly higher than the 2018/02 epoch SPHERE measurement in the filter (2.25m).
3 Spectral analysis
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 and SINFONI -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 circumplanetary disc with a variable accretion rate (Type III models; Section 3.2). In order to minimise 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/05), SPHERE (2016/05) and NICI (2012/03) epochs, and the other for all other epochs. This division was chosen because (1) the SPHERE IFS ()+IRDIS (/) were obtained simultaneously on 2018/02, and the IRDIS / points of 2015/05 are consistent with the IFS spectrum, and (2) the SINFONI 2014/05 data are brighter than the IRDIS point of 2018/02, while the IRDIS point of 2016/05 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 minimised the following goodness-of-fit indicator:
where is the uncertainty in the flux measurement at wavelength , and weights are defined for photometric and spectroscopic observations following a similar strategy as in Ballering2013 and Olofsson2016. 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 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 (Allard2012; Baraffe2015), and the grid of atmospheric models presented in Madhusudhan2011. 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, , explored between 3.0 and 5.0 (steps of 0.5 dex). We assumed Solar metallicity.
Madhusudhan2011 considered a wide grid of cloud models, with different geometrical and optical thickness, particle size and metallicities. The Madhusudhan2011 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 taking possible values of , or cm s, and the sedimentation parameter (defined as in Ackerman2001). Madhusudhan2011 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 Madhusudhan2011 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 in steps of 0.1 . We also considered dust extinction as an additional free parameter, with allowed values between and 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 (Draine1989). 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. Marocco2014). We assumed a distance of 113 pc (Gaia2018).
Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefBest fit parameters for PDS 70 b
|I. Planet alone||II. Planet alone (variable extinction)||III. Planet and circumplanetary disc|
|Parameter||Range||Best fit||Parameter||Range||Best fit||Parameter||Range||Best fit|
|BT-SETTL atmospheric models|
|M11 atmospheric models|
|is inferred from the best-fit and values, and is inferred from and the best-fit .|
|See Section 3.1 and Madhusudhan2011. A60 refers to model cloud A (thickest) with a modal particle size of 60 m.|
Figures 3.1 and 3.1 compare our best-fit models from the BT-SETTL and Madhusudhan2011 grids, respectively, with the SED of PDS 70 b. Table 3.1 gives the explored parameter ranges, best-fit parameter values and corresponding , for each model. Best-fit Type I BT-SETTL and Madhusudhan2011 models reproduce most of the observed SED but are a poor match to the red end of the SINFONI spectrum, with a discrepancy for most data points at wavelengths m. They lead to reduced goodness-of-fit indicators and 1.20.
Allowing for variable extinction (Type II) yields best-fit models (shaded cyan areas in Figures 3.1 and 3.1) in better agreement with flux estimates longward of m ( and 0.70). Upper and lower edges of the shaded areas correspond to the minimum and maximum extinction values ( and , in Table 3.1) of our best-fit model, with better accounting for the 2014/05 SINFONI, 2016/05 SPHERE and 2012/03 NICI data points and accounting for data points at all other epochs. The difference in extinction is mag for both BT-SETTL and Madhusudhan2011 models.
For both Type I and II models, the best-fit effective temperatures (1100–1500 K), surface gravity (–), planet radius (–) and hence mass (1.7–42.0 ) are in approximate agreement with the previous estimates made in Muller2018. Our mass estimates are uncertain because of the large steps in (0.5 dex) in our model grids, and hence do not rule out a brown dwarf. Our best-fit Madhusudhan2011 models correspond to the thickest cloud models (labelled A); extending to the top of the atmosphere.
3.2 Combined atmospheric+CPD models
For our Type III models (combined atmosphere+CPD emission), we considered the same two grids of atmospheric models as in Section 3.1, coupled with the CPD models presented in Eisner2015. The latter add a single free parameter, the mass accretion rate, which sets the brightness and shape of the CPD spectrum. We explored values of mass accretion rates ranging from -7.0 to -6.0, in steps of 0.1 dex. We did not consider accretion rates smaller than yr because the corresponding models do not contribute significantly at NIR wavelengths. We assume a fixed inner truncation radius of in our CPD models, as in Eisner2015. We thus truncated our grid of planetary radii to for consistency. Other parameters were explored on the same grids as for Type I and II models.
Figures 3.1 and 3.1 show the best-fit combined planet+CPD models (shaded yellow areas). Dashed lines show the contribution from the atmosphere alone. Upper and lower edges of the shaded areas correspond to maximum and minimum mass accretion rates of our best fit model ( and , respectively, in Table 3.1), accounting for the 2014/05 SINFONI and 2016/05 SPHERE data points and, respectively, data points at all other epochs. The planet+CPD best-fit models reproduce better the observed spectrum than pure atmospheric models (with or without variable accretion), with reduced goodness-of-fit indicators and 0.44 using BT-SETTL and Madhusudhan2011 models, respectively. Interestingly, the best-fit parameters for both the CPD and the planet are similar using either BT-SETTL or Madhusudhan2011 models: mass accretion rates ranging between and yr, effective temperature of 1500–1600 K, surface gravity , radius of , and mass of . The Madhusudhan2011 best-fit model also suggests the thickest cloud geometry, with a modal particle size of 60 m.
The estimate of 10 is larger than that inferred from the planet-only BT-SETTL models because the estimated is significantly larger, while the inferred is slightly smaller. For Madhusudhan2011 models the opposite is true because is similar but is smaller. This suggests an older planet when considering a CPD in the model. Interestingly, the estimated planet parameters (, , , ) agree with the BT-SETTL models for a mass of and age 9–11 Myr (Baraffe2015)333Available at https://phoenix.ens-lyon.fr/Grids/BT-Settl/. The inferred age is consistent with estimates for the star in Pecaut2016, but not with the newer estimate of Myr (Muller2018). In contrast, parameters in planet-only models are inconsistent with Baraffe2015 evolutionary models for any combination of mass and age.
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/02 SPHERE measurement and the points at m in our 2014/05 SINFONI spectrum. This is further supported by the disagreement (albeit slight) between the SPHERE and SINFONI measurements at m. Since the variability of classical T-Tauri stars is thought to be related to either variable amounts of extinction from intervening circumstellar dust (e.g. Bozhinova2016) or irregular accretion (e.g. Bouvier2004; Rigon2017), 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. Gressel2013), 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:
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 Madhusudhan2019, 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 (Madhusudhan2011; Bonnefoy2013; Currie2014). Given that our best-fit type I (purely atmospheric) BT-SETTL and Madhusudhan2011 models are similar to the reddest atmospheric models (around m) presented in Muller2018, we do not expect our conclusion of an excess at m to change using different grids of atmospheric models.
The fit is not perfect. Although the best-fit atmosphere+CPD model best reproduces the excess at the end of the band, it does not perfectly reproduce the observed slope around 2.3m. For both the BT-SETTL and Madhusudhan2011 type III models, most photometric points at wavelengths 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).
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 2. As shown in Zhu2015b, different inner truncation radii can lead to different predicted CPD spectra for a given mass accretion rate. Furthermore, neither the models in Zhu2015b nor in Eisner2015 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 . Montesinos2015 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:
Our best-fit accretion rates agree with observations. Assuming a similar relationship between H luminosity and mass accretion rate as T-Tauri stars, Wagner2018a estimated the PDS 70 b accretion rate as yr 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 yr for mag. Assuming a similar relationship for Br emission, Christiaens2019a also constrained the accretion rate to be yr, considering negligible extinction at band. Our best-fit models suggest strong extinction towards the protoplanet (mag), consistent with the presence of a CPD. Our estimates of mass (), accretion rates ( yr), extinction ( mag) and radii (), appear roughly compatible with both observational constraints. Monitoring of the H luminosity would confirm the variability of the accretion rate.
The photometric variability is only observed at relatively long NIR wavelength (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 (au), both the amplitude ( mag) and timescale of extinction variability (less than several years) suggests it would also be caused by circumplanetary dust.
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 (; first transitions at 2.294, 2.323 and 2.352 m) is an indicator of disc presence (Geballe1987; Davis2011). CPD models in Ayliffe2009 and Szulagyi2017a predict temperatures up to several thousand K, which might also produce CO bandhead emission.
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 (Christiaens2019a).
Finally, while this work was in the final stages of review, Szulagyi2019 provided new IR flux predictions for CPDs, based on hydro-dynamical simulations followed by radiative transfer, which are in qualitative agreement with the CPD contribution we inferred, in the case of a planet.
The SED of PDS 70 b is best fit by models that include a circumplanetary disc. Atmospheric models alone are not able to account for the observed flux at wavelengths m. We infer an accretion rate of – yr for a planet with significant extinction, consistent with prior observations. Simultaneous follow-up observations of PDS 70 b with wide spectral coverage in NIR including the H line should confirm the scenario of variable accretion through a circumplanetary disc.
We acknowledge funding from the Australian Research Council via DP180104235, FT130100034 and FT170100040. VC thanks Andre Müller for sharing the SPHERE spectrum of PDS 70 b. MM acknowledges financial support from the Chinese Academy of Sciences (CAS) through CASSACA-CONICYT Postdoctoral Fellowship (Chile).
Facilities: VLT(SINFONI,SPHERE), MASSIVE (www.massive.org.au)
Software: ANDROMEDA (Cantalloube2015), VIP (GomezGonzalez2017), SpeX (http://www.browndwarfs.org/spexprism)