---
abstract: 'We present the first observational evidence for a circumplanetary disk around the protoplanet PDS 70 b, based on a new spectrum in the $K$ 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 disk model. Goodness-of-fit indicators favour the third option, suggesting circumplanetary material contributing excess thermal emission — most prominent at $\lambda \gtrsim 2.3 \mu$m. Inferred accretion rates ($\sim 10^{-7.8}$–$10^{-7.3} M_J$ yr$^{-1}$) are compatible with observational constraints based on the H$\alpha$ and Br$\gamma$ lines. For the planet, we derive an effective temperature of 1500–1600 K, surface gravity $\log(g)\sim 4.0$, radius $\sim 1.6 R_J$, mass $\sim 10 M_J$, and possible thick clouds. Models with variable extinction lead to slightly worse fits. However, the amplitude ($\Delta A_V \gtrsim 3$mag) and timescale of variation ($\lesssim$ years) required for the extinction would also suggest circumplanetary material.'
author:
- Valentin Christiaens
- Faustine Cantalloube
- Simon Casassus
- 'Daniel J. Price'
- Olivier Absil
- Christophe Pinte
- Julien Girard
- Matias Montesinos
bibliography:
- 'PDS70.bib'
title: Evidence for a circumplanetary disk around protoplanet PDS 70 b
---

Introduction
============

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 disk (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 ; see @Madhusudhan2019 for a recent review).

We present evidence for a circumplanetary disk 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 [@Hashimoto2012; @Dong2012]. The protoplanet was first detected in this gap using VLT/SPHERE [@Keppler2018; @Muller2018 hereafter ]. Subsequent search for sub-mm emission from a CPD using ALMA was inconclusive [@Keppler2019]. Spectral characterisation of the planet was presented in , suggesting an effective temperature of $\sim$1000-1500 K, surface gravity $\lesssim 3.5$dex and mass $\leq 17 M_{\rm J}$. We present a new spectral characterization of PDS 70 b, including both the measurements presented in and a new spectrum in $K$ band obtained with VLT/SINFONI [@Christiaens2019a hereafter ]. The spectrum shows excess emission at $\lambda \gtrsim 2.3 \mu$m inconsistent with naked model atmospheres of young planets, indicating the presence of circumplanetary material.

Observations {#Obs+DataRed}
============

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 , but considered updated values for their $L'$-band flux estimates. fitted the 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 compared to the star likely arises from hot dust in the inner disk that would be unresolved from the star at $L'$ band, and should hence be included. By contrast, the contribution from resolved scattered light from the disk at $L'$ 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 $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 \mu$m) and W2 ($4.60 \mu$m) filters of WISE [@Wright2010], and (2) the same values for the contrast of the companion as in . The new $L'$ estimates of the companion are: $(8.55 \pm 3.38) \times 10^{-17}$ W m$^{-2}$ $\mu$m$^{-1}$ and $(6.70 \pm 4.05) \times 10^{-17}$ W m$^{-2} \mu$m for the NICI and NACO data, respectively.

In a recent paper , we inferred the contrast of the protoplanet with respect to the star as a function of wavelength, $c(\lambda)$, in the $K$ band at unprecedented spectral resolution ($R\sim100$ after spectral binning) using VLT/SINFONI. We employed two different methods for extracting $c(\lambda)$, 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) disk signals.

Despite the higher quality of the ANDROMEDA $c(\lambda)$, some spectral channels contain outliers. In order to minimize the risk of bias, we first removed spectral channels with a detection below 3$\sigma$ and lying in strong telluric lines, then used a Savitzky-Golay filter of order 3 with a 81-channel window to smooth the $c(\lambda)$ curve before binning it by a factor of 20 [@Savitzky1964]. We obtained the final $K$-band spectrum of the protoplanet by multiplying $c(\lambda)$ 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 and 37 obtained with SINFONI . 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 $K2$ filter (2.25$\mu$m).

Spectral analysis {#Kspectrum}
=================

We first attempted to fit the observed SED with synthetic spectra modeling pure atmospheric emission (Section \[PhotosphericModels\]). 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 circumplanetary disk with a variable accretion rate (*Type III models*; Section \[CPDModels\]). 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 ($YJH$)+IRDIS ($K1$/$K2$) were obtained simultaneously on 2018/02, and the IRDIS $H1$/$H2$ points of 2015/05 are consistent with the IFS spectrum, and (2) the SINFONI 2014/05 data are brighter than the IRDIS $K2$ point of 2018/02, while the IRDIS $K2$ 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: $$\begin{aligned}
\label{Eq:GoodnessOfFit}
\chi^2 = \sum_i \omega_i \big[ \frac{F_{\rm obs}(\lambda_i) - F_{\rm model} (\lambda_i)}{\sigma_i} \big]^2\end{aligned}$$ where $\sigma_i$ is the uncertainty in the flux measurement $F_{\rm obs}(\lambda_i)$ at wavelength $\lambda_i$, and weights $\omega_i$ 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 $\chi_r^2$ as $\chi^2$ divided by the respective number of degrees of freedom for each type of model.

Atmospheric models {#PhotosphericModels}
------------------

![image](Fig1_v2.png){width="\textwidth"}

![image](Fig2_v2.png){width="\textwidth"}

We considered two grids of synthetic spectra: BT-SETTL models [@Allard2012; @Baraffe2015], and the grid of atmospheric models presented in . 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, $\log(g)$, explored between 3.0 and 5.0 (steps of 0.5 dex). We assumed Solar metallicity.

considered a wide grid of cloud models, with different geometrical and optical thickness, particle size and metallicities. The 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 $\mu$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^2$, $10^4$ or $10^6$ cm$^2$ s$^{-1}$, and the sedimentation parameter $f_{\rm sed}$ [defined as in @Ackerman2001]. 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 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_{\rm J}$ in steps of 0.1 $R_{\rm 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 [@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].

[lcclcclcc]{} & &\
Parameter & Range & Best fit & Parameter & Range & Best fit & Parameter & Range & Best fit\
\
$T_{\mathrm{eff}}$ \[K\] & 1200–1900 & 1500 & $T_{\mathrm{eff}}$ \[K\] & 1200–1900 & 1500 & $T_{\mathrm{eff}}$ \[K\] & 1200–1900 & 1500\
$\log(g)$ & 3.0–5.0 & 3.0 & $\log(g)$ & 3.0–5.0 & 3.0 & $\log(g)$ & 3.0–5.0 & 4.0\
$R_{\rm b}$ \[$R_{\rm J}$\] & 0.1–5.0 & 2.2 & $R_{\rm b}$ \[$R_{\rm J}$\] & 0.1–5.0 & 2.1 & $R_{\rm b}$ \[$R_{\rm J}$\] & 0.1–2.0 & 1.6\
$A_V$ \[mag\] & 0.0–10.0 & 4.34 & $A_{V\mathrm{,1}}$ \[mag\] & 0.0–10.0 & 0.60 & $A_V$ \[mag\] & 0.0–10.0 & 6.40\
& & & $A_{V\mathrm{,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}$\
$M_{\rm b}^{\dagger}$ \[$M_{\rm J}$\] & – & 1.9 & $M_{\rm b}^{\dagger}$ \[$M_{\rm J}$\] & – & 1.7 & $M_{\rm b}^{\dagger}$ \[$M_{\rm J}$\] & – & 9.9\
& & & & & & $\dot{M}_{\rm b}^{\dagger}$ \[$M_{\rm J}$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$\
\
$T_{\mathrm{eff}}$ \[K\] & 600–1800 & 1100 & $T_{\mathrm{eff}}$ \[K\] & 600–1800 & 1200 & $T_{\mathrm{eff}}$ \[K\] & 600–1800 & 1600\
$\log(g)$ & 3.5–5.5 & 4.0 & $\log(g)$ & 3.5–5.5 & 4.0 & $\log(g)$ & 3.5–5.5 & 4.0\
$R_{\rm b}$ \[$R_{\rm J}$\] & 0.1–5.0 & 3.3 & $R_{\rm b}$ \[$R_{\rm J}$\] & 0.1–5.0 & 2.8 & $R_{\rm b}$ \[$R_{\rm J}$\] & 0.1–2.0 & 1.6\
$c$ & \[NC,E,A,...\]$^{\ddagger}$ & A & $c$ & \[NC,E,A,...\]$^{\ddagger}$ & A & $c$ & \[NC,E,A,...\]$^{\ddagger}$ & A60$^{\ddagger}$\
$f_{\rm sed}$ & \[eq.,0,1,2\] & eq. & $f_{\rm sed}$ & \[eq.,0,1,2\] & eq. & $f_{\rm sed}$ & \[eq.,0,1,2\] & eq.\
$K$ \[cm$^2$s$^{-1}$\] & \[eq.,$10^2$–$10^6$\] & eq. & $K$ \[cm$^2$s$^{-1}$\] & \[eq.,$10^2$–$10^6$\] & eq. & $K$ \[cm$^2$s$^{-1}$\] & \[eq.,$10^2$–$10^6$\] & eq.\
$A_V$ \[mag\] & 0.0–10.0& 3.00 & $A_{V\mathrm{,1}}$ \[mag\] & 0.0–10.0 & 0.80 & $A_V$ \[mag\] & 0.0–10.0 & 8.72\
& & & $A_{V\mathrm{,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}$ \[$M_{\rm J}^2$yr$^{-1}$\] & $10^{-7}$–$10^{-6}$ & $10^{-6.7}$\
$M_{\rm b}^{\dagger}$ \[$M_{\rm J}$\] & – & 42.0 & $M_{\rm b}^{\dagger}$ \[$M_{\rm J}$\] & – & 30.2 & $M_{\rm b}^{\dagger}$ \[$M_{\rm J}$\] & – & 9.9\
& & & & & & $\dot{M}_{\rm b}^{\dagger}$ \[$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$\
\
\

Figures \[FinalSpectrum\_BTSETTL\] and \[FinalSpectrum\_M11\] compare our best-fit models from the BT-SETTL and grids, respectively, with the SED of PDS 70 b. Table \[tab:best\_fit\_params\] gives the explored parameter ranges, best-fit parameter values and corresponding $\chi_r^2$, for each model. Best-fit Type I BT-SETTL and models reproduce most of the observed SED but are a poor match to the red end of the SINFONI spectrum, with a $>2\sigma$ discrepancy for most data points at wavelengths $\gtrsim 2.3 \mu$m. They lead to reduced goodness-of-fit indicators $\chi_r^2 \sim 1.01$ and 1.20.

Allowing for variable extinction (Type II) yields best-fit models (shaded cyan areas in Figures \[FinalSpectrum\_BTSETTL\] and \[FinalSpectrum\_M11\]) in better agreement with flux estimates longward of $\sim 2.3 \mu$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 \[tab:best\_fit\_params\]) of our best-fit model, with $A_{V,1}$ better accounting for the 2014/05 SINFONI, 2016/05 SPHERE and 2012/03 NICI data points and $A_{V,2}$ accounting for data points at all other epochs. The difference in extinction is $A_V \gtrsim 3$ mag for both BT-SETTL and models.

For both Type I and II models, the best-fit effective temperatures (1100–1500 K), surface gravity ($\log(g)\sim3.0$–$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 . Our mass estimates are uncertain because of the large steps in $\log(g)$ (0.5 dex) in our model grids, and hence do not rule out a brown dwarf. Our best-fit models correspond to the thickest cloud models (labelled *A*); extending to the top of the atmosphere.

Combined atmospheric+CPD models {#CPDModels}
-------------------------------

For our Type III models (combined atmosphere+CPD emission), we considered the same two grids of atmospheric models as in Section \[PhotosphericModels\], 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 $\log(\dot{M}_b M_b [M_J^2 \textrm{yr}^{-1}])$ ranging from -7.0 to -6.0, in steps of 0.1 dex. We did not consider accretion rates smaller than $10^{-7}M_J^2$ yr$^{-1}$ because the corresponding models do not contribute significantly at NIR wavelengths. We assume a fixed inner truncation radius of $2 R_{\rm J}$ in our CPD models, as in @Eisner2015. We thus truncated our grid of planetary radii to $2 R_{\rm J}$ for consistency. Other parameters were explored on the same grids as for Type I and II models.

Figures \[FinalSpectrum\_BTSETTL\] and \[FinalSpectrum\_M11\] 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 ($M_b\dot{M}_{b,1}$ and $M_b\dot{M}_{b,2}$, respectively, in Table \[tab:best\_fit\_params\]), 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 $\chi_r^2 \sim 0.41$ and 0.44 using BT-SETTL and models, respectively. Interestingly, the best-fit parameters for both the CPD and the planet are similar using either BT-SETTL or models: mass accretion rates ranging between $\sim 10^{-6.4}$ and $\sim 10^{-6.8}$ $M_J^2$ yr$^{-1}$, effective temperature of 1500–1600 K, surface gravity $\log(g)\sim 4.0$, radius of $\sim 1.6 R_J$, and mass of $\sim 10 M_J$. The best-fit model also suggests the thickest cloud geometry, with a modal particle size of 60 $\mu$m.

The estimate of 10 $M_J$ is larger than that inferred from the planet-only BT-SETTL models because the estimated $\log(g)$ is significantly larger, while the inferred $R_b$ is slightly smaller. For models the opposite is true because $\log(g)$ is similar but $R_b$ is smaller. This suggests an older planet when considering a CPD in the model. Interestingly, the estimated planet parameters ($T_{\rm eff}$, $\log(g)$, $R_b$, $M_b$) agree with the BT-SETTL models for a mass of $10 M_J$ and age 9–11 Myr [@Baraffe2015][^1]. The inferred age is consistent with estimates for the star in @Pecaut2016, but not with the newer estimate of $5.4 \pm 1.0$ Myr . In contrast, parameters in planet-only models are inconsistent with @Baraffe2015 evolutionary models for any combination of mass and age.

Discussion {#Discussion}
==========

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 $K2$ measurement and the points at $\gtrsim 2.3 \mu$m in our 2014/05 SINFONI spectrum. This is further supported by the disagreement (albeit slight) between the SPHERE and SINFONI measurements at $\sim 2.25 \mu$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:

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 @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 [; @Bonnefoy2013; @Currie2014]. Given that our best-fit type I (purely atmospheric) BT-SETTL and models are similar to the reddest atmospheric models (around $\sim 2.3 \mu$m) presented in , we do not expect our conclusion of an excess at $\gtrsim 2.3 \mu$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$\mu$m. For both the BT-SETTL and type III models, most photometric points at wavelengths shorter than 2.3 $\mu$m lie below the model (albeit all within 2$\sigma$), while most points longward of 2.3 $\mu$m are slightly brighter than the model (but within 2$\sigma$ 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 2$R_J$. 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 $L_p\sim10^{-4} L_{\odot}$. @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:

1.  *Our best-fit accretion rates agree with observations.* Assuming a similar relationship between H$\alpha$ luminosity and mass accretion rate as T-Tauri stars, @Wagner2018a estimated the PDS 70 b accretion rate as $\sim 10^{-7.8} M_{\rm J}^2 R_{\rm J}^{-1}$ 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}$ yr$^{-1}$ for $A_V \approx 3$mag. Assuming a similar relationship for Br$\gamma$ emission, also constrained the accretion rate to be $< 10^{-6.2} M_{\rm J}^2 R_{\rm J}^{-1}$ yr$^{-1}$, considering negligible extinction at $K$ band. Our best-fit models suggest strong extinction towards the protoplanet ($A_V > 6$mag), consistent with the presence of a CPD. Our estimates of mass ($M_b \sim 10 M_J$), accretion rates ($\sim 10^{-7.5} M_J$ yr$^{-1}$), extinction ($\gtrsim 6$ mag) and radii $R_{\rm b}$ ($\sim 1.6 R_J$), appear roughly compatible with both observational constraints. Monitoring of the H$\alpha$ luminosity would confirm the variability of the accretion rate.

2.  *The photometric variability is only observed at relatively long NIR wavelength ($\sim 2.2 \mu$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 $\mu$m. Even if the variability was due to extinction, given the radial separation of the protoplanet ($\sim 20$au), both the amplitude ($\Delta A_V \gtrsim 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 $\mu$m might suggest CO bandhead emission.* For young stellar objects, CO bandhead emission ($\Delta v=2$; first transitions at 2.294, 2.323 and 2.352 $\mu$m) is an indicator of disk presence [@Geballe1987; @Davis2011]. CPD models in @Ayliffe2009 and @Szulagyi2017a 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 .

Conclusions
===========

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 $\gtrsim~2.3~\mu$m. We infer an accretion rate of $10^{-7.8}$–$10^{-7.3} M_J$ yr$^{-1}$ for a $\sim~10~M_J$ 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$\alpha$ line should confirm the scenario of variable accretion through a circumplanetary disc.

Acknowledgments {#acknowledgments .unnumbered}
===============

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

[^1]: Available at <https://phoenix.ens-lyon.fr/Grids/BT-Settl/>
