Issue
A&A
Volume 644, December 2020
Article Number A13
Number of page(s) 18
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202038878
Published online 26 November 2020

© ESO 2020

1 Introduction

The formation of planets occurs in circumstellar disks (CSDs) around pre-main sequence stars. Spatially resolved observations have revealed a ubiquity of substructures in the gas and dust distribution of those disks, such as gaps and spiral arms (e.g., Andrews et al. 2018; Avenhaus et al. 2018). These features may point to the gravitational interaction of embedded planets with their natal environment (e.g., Pinilla et al. 2012; Dong et al. 2015), but the direct detection of these potential protoplanets remains challenging (e.g., Currie et al. 2019; Cugno et al. 2019), possibly due to their low intrinsic brightness and extinction effects by dust (e.g., Brittain et al. 2020; Sanchis et al. 2020). Nevertheless, direct detections of forming planets are critical to advance our empirical understanding of the physical processes by which planets accumulate gas and dust from, and interact with, their circumstellar environment.

The CSD of PDS 70 is a unique example in which two embedded planets were directly detected with high-resolution instruments. PDS 70 is a weak-line T Tauri, K7-type (Pecaut & Mamajek 2016) star with an estimated age of 5.4 ± 1.0 Myr (Müller et al. 2018); it is surrounded by a gapped accretion disk (Hashimoto et al. 2012) and is located in the Scorpius-Centaurus OB association (Gregorio-Hetem & Hetem 2002; Preibisch & Mamajek 2008). Keppler et al. (2018) discovered PDS 70 b within the gap of the disk with SPHERE and archival L′ band data. This planet is located at a favorable position (close to the major axis of the disk) where projection and extinction effects are minimized. Later, a second planetary companion, PDS 70 c, was discovered by Haffert et al. (2019) in Hα, together with a detection of Hα emission from PDS 70 b (see also Wagner et al. 2018). Such measurements of hydrogen emission lines place constraints on the physics of the accretion flow and shock, extinction, and the mass and accretion rate of the planets (Thanathibodee et al. 2019; Aoyama & Ikoma 2019; Hashimoto et al. 2020).

While PDS 70 b and c have been suggested to be planetary-mass objects and have been confirmed to be comoving with PDS 70 (Keppler et al. 2018; Müller et al. 2018; Haffert et al. 2019; Mesa et al. 2019), their atmospheric and circumplanetary characteristics remain poorly understood. The H and K band fluxes of PDS 70 b are consistent with a mid L-type object, but its near-infrared (NIR) colors are redder than those of field dwarfs and low-gravity objects (Keppler et al. 2018). A comparison with cloudy atmosphere models by Müller et al. (2018) shows that a wide range of temperatures (1000–1600 K) and radii (1.4–3.7 RJ) could describe the spectral energy distribution (SED) from the Y to L′ bands. A detailed SED analysis by Christiaens et al. (2019) has revealed an excess of the K band emission with respectto predictions by atmospheric models. The authors show that the SED is consistent with a combination of emission from a planet atmosphere (1500–1600 K) and a circumplanetary disk (CPD). Most recently, Wang et al. (2020) presented NIRC2 L′ imaging and analyzed the SED with atmospheric models and blackbody spectra. From this, the authors conclude that the data is best described by a blackbody spectrum with Teff=120453+52${T_{\textrm{eff}}}\,{=}\,1204^{+52}_{-53}$ K and R=2.720.34+0.39$R\,{=}\,2.72^{+0.39}_{-0.34}$ RJ.

In this work, we report the first detection of PDS 70 b in the 4–5 μm range as part of the MIRACLES survey (Stolker et al. 2020). The object was observed with NACO at the Very Large Telescope (VLT) in Chile and detected with both the Brα and M′ filters. Mid-infrared (MIR) wavelengths are in particular critical to uncover potential emission from a circumplanetary environment. We will analyze the photometry, colors, and SED of the object to gain insight into its main characteristics.

2 Observations and data reduction

2.1 High-contrast imaging with VLT/NACO

We observed PDS 70 with VLT/NACO (Lenzen et al. 2003; Rousset et al. 2003) in the NB4.05 (Brα; λ0 = 4.05 μm, Δλ = 0.02 μm) and M′ (λ0 = 4.78 μm, Δλ = 0.59 μm) filters (ESO program ID: 0102.C-0649(A)) as part of the MIRACLES survey, which aims at the systematic characterization of directly imaged planets and brown dwarfs at 4–5 μm (Stolker et al. 2020). The data were obtained without coronagraph in pupil-stabilized mode, while dithering the star across the detector to sample the thermal background emission. The total telescope time was 3 and 4 h for the NB4.05 and M′ filters, respectively, split over multiple nights with observing blocks (OBs) of 1 h each. This resulted in a total of 1.75 and 2 h of on-source telescope time for NB4.05 and M′. A detailed description of the observing strategy for the MIRACLES survey is available in Stolker et al. (2020) but a few specifics for the observations of PDS 70 are provided here.

The observations with the NB4.05 filter were executed on UT 2019 February 23 and UT 2019 March 15. A detector integration time (DIT) of 1.0 s and NDIT of 61 or 65 was used, resulting in 1680 (first and second OB) and 1792 (third OB) frames. During the first night, two OBs were executed in good conditions (seeing ≲0.′′8) while during the second night (i.e., the third OB), the conditions were slightly worse (0.′′ 75–0.′′95), resulting in an average angular resolution of 115 mas (1 FWHM). Aperture photometry (2 FWHM in diameter) of the star revealed flux variations of 4.6% across the three datasets, which in particular reflects the variable conditions during the third OB. The total, non-intermittent, and non-overlapping field rotation was 50° but gaps in the parallactic angle range between OBs helped with minimizing the self-subtraction during post-processing.

With asimilar setup, we observed the target with the M′ filter on UT 2019 February 20, 21, and 22, with two OBs executed during the second night. The detector was windowed to a field of view of 256 × 256 pixels to allow for a short integration time of 35 ms without frame loss. With an NDIT of 1500 integrations and 14 exposures (i.e., data cubes) for each of the two dithering positions, this resulted in 42 000 frames per OB. The seeing was approximately stable during three of the observations with average values in the range of 0.′′ 7–0.′′8. During the second OB, the seeing was about 1.′′ 0–1.′′2 with a short increase to 2.′′0. After a frame selection and combining the data from the four OBs, the stellar flux varied by about 6.5% and the FWHM of the PSF was 134 mas. The total, continuous field rotation was 56°, but 82° if the gaps in the parallactic angle coverage are included.

2.2 Data processing and calibration

The data were processed with PynPoint1 which is a generic, end-to-end pipeline for high-contrast imaging data (Amara & Quanz 2012; Stolker et al. 2019). We used the latest release of the package (version 0.8.3) for the pre- and post-processing, and the relative photometric and astrometric calibration. The pre-processing was done for each dataset separately and the frames from the different OBs were combined before the PSF subtraction. We used an implementation of full-frame principal component analysis (PCA; Amara & Quanz 2012; Soummer et al. 2012) to remove the quasi-static structures of the stellar PSF.

In general, we followed the processing and calibration procedure that is described in Stolker et al. (2020). However, in addition to subtracting the mean background (based on the adjacent data cubes in which the star was dithered to a different detector position) we also applied an additional correction with PCA (Hunziker et al. 2018). Specifically, we decomposed the stack of all background images (after subtracting the mean of the stack) at a given dithering position and projected the science data on the first principal component (PC). The central region (8 FWHM in diameter) was masked during the projection but included when subtracting the model. This provided better results on visual inspection compared to a mean background subtraction alone.

After pre-proccessing and combining the OBs, subsets of 8 (NB4.05) and 330 (M′) images weremean-collapsed, resulting in a final stack of 501 and 502 images for NB4.05 and M′, respectively.Then, we extracted the photometry and astrometry of the companions relative to their star in the following way. First, the dependence on the number of PCs was tested (1–5 PCs for NB4.05 and 1–10 PCs for M′), second, we used an MCMC approach to estimate the statistical uncertainty for a fixed number of PCs by removing the planet signal with a negative copy of the PSF, and thirdly, a bias and systematic uncertainty was estimated by injecting and retrieving artificial planets (see Stolker et al. 2020 for details). For the calibration, we used a field of view of 57 pixels, we subtracted three (NB4.05) and five (M′) PCs, and we applied an one-on-one injection of the PSF templates (the stellar flux had remained within the linear regime of the detector). The estimation of a potential bias and systematic error is challenging since the planet is only at 1.5 λ/D in M′, and both disk signal and noise residuals are present at the same separation (see Fig. 1). Therefore, to not introduce a bias, we excluded position angles with relatively bright disk or noise features for the estimation of the systematic error (see Table 1).

From the relative calibration, we determined the apparent magnitudes in the NB4.05 and M′ filters. We first used the species2 toolkit (Stolker et al. 2020)to convert the 2MASS JHK, and WISE W1 and W2 magnitudes of the PDS 70 system into fluxes. We then fitted a power law function to these in log-log space. The stellar magnitudes in NB4.05 and M′ were then computed by integrating the model spectrum across the filter profiles (see Table 1). We note that this approach assumes that the photometry in the considered spectral range is dominated by continuum emission from the star and inner disk. Therefore, potential Brα emission due to accretion onto the star is ignored, but that is a reasonable assumption given the low accretion rate of (0.6–2.2) × 10−10 M yr−1 for PDS 70 (Thanathibodee et al. 2020).

thumbnail Fig. 1

Detection of the PDS 70 planetary system and CSD with the NACO NB4.05 (left panel) and M′ (right panel) filters. The images show the mean-combined residuals of the PSF subtraction on a color scale that has been normalized to the peak of the stellar PSF. The flux in the NB4.05 image has been increased by a factor of 1.8 for clarity. The detected emission from PDS 70 b and c (only marginal in NB4.05) is encircled. North is up and east is left.

Table 1

Photometry and error budget.

2.3 Reanalysis of archival data

In addition to the new NB4.05 and M′ data, we reanalyzed archival NACO L′ data from Keppler et al. (2018) (ESO program ID: 097.C-0206(A)), in line with the systematic 3–5 μm analysis for the MIRACLES survey, and additionally coronagraphic SPHERE H23 and K12 data from Keppler et al. (2018) and Müller et al. (2018) (ESO program IDs: 095.C-0298(A) and 1100.C-0481(D)). Below, we provide a few details on the data quality and processing, but we refer to the respective papers for more information on these datasets. The calibration was done in a similar way as the NB4.05 and M′ data. While the H band is also covered by the NIR spectrum that we adopted from Müller et al. (2018) (see Sect. 3.3.1), the K band flux was in particular critical for estimating the photospheric temperature and radius of PDS 70 b (see Sect. 3.3).

The NACO L′ data were obtained with seeing conditions in the range of 0.′′ 55–0.′′7. We removed 18% of the frames (based on aperture photometry at the position of the star) after which the stellar flux varied by 29% across the dataset. The flux had not saturated the detector so we applied an one-on-one PSF injection during the relative calibration. Therefore, the variation in the stellar flux is not expected to have introduced a bias in the extracted planet flux. A total of 14 464 frames were selected across a parallactic rotation of 85°.

The archival H23 and K12 datasets had been obtained with the IR dual-band imager (IRDIS; Dohlen et al. 2008; Vigan et al. 2010) of SPHERE (Beuzit et al. 2019). We analyzed both H23 epochs fromKeppler et al. (2018) but only use the results from UT 2015 May 04 since the second dataset (UT 2015 June 01) was obtained in poor observing conditions with a seeing larger than 1′′. During the first epoch, the seeing was 0.′′ 35 at the start of the observations, but degraded to >1′′ at 1/3 of the sequence. The stellar halo appeared bright and asymmetric, possibly due to a low-wind effect (~5 m s−1) and/or a wind-driven halo (Cantalloube et al. 2018). We only used 30 frames that were obtained in good conditions, which were selected by measuring the flux of the background star at ~2.′′4 north of PDS 70. Similarly, we only used the off-axis PSF exposures from the start of the observations because these were obtained in conditions that were similar to the selected frames with the star behind the coronagraph. The flux in the unsaturated PSF exposures has been scaled to the coronagraphic data to account for the difference in exposure time and the transmission of the neutral density filter.

There are two archival SPHERE/IRDIS K12 datasets available, which had been obtained on UT 2016 May 15 and 2018 February 25. We analyzed both datasets but only used the results from the second epoch because the assessment of the first epoch revealed large-scale noise residuals after the PSF subtraction, which may have biased the photometry. The second dataset was obtained in good observing conditions but the seeing degraded toward the end of the sequence. Therefore, similar to the H23 data, we selected 24 frames from the start of the sequence, based on the photometry of the background star, and the unsaturated exposures from the start of the observations.

3 Results

3.1 Detection of the PDS 70 system

The mean-combined residuals from the PSF subtraction after subtracting two (NB4.05) and three (M′) PCs are presented in Fig. 1. The choice of the number of PCs for the image is dictated by the brightness of the planets; to characterize them a somewhat larger number of PCs was removed (see Sect. 2.2) to better suppress the residual speckle noise. The images reveal a bright source at the expected position of PDS 70 b.

While planet b is visible in both filters, planet c is only marginally detected with the NB4.05 filter and not visible in the M′ image. Here, the position of planet c relative to the near side of the disk may have prevented a detection in M′ due to the reduced angular resolution compared to NIR wavelengths. In the NB4.05 image, planet c is blended with the disk signal, and therefore the extracted flux is potentially biased. We estimated the bias due to the CSD signal by injecting and retrieving the contrast of an artificial planet at a location with comparable disk flux but somewhat offset from the c planet, yielding an approximate correction of ~0.1 mag (see Table 1).

The results from the photometric extraction of the companions are listed in Table 1, both for the new and archival data. The final contrast is calculated by adding the bias offset and combining the error components in quadrature. The error budget of the planet photometry is dominated by the error from the relative calibration while the error on the stellar magnitude (expected to be a few tens of a magnitude) is negligible. For the coronagraphic SPHERE H23 and K12 data, we have included an additional error component that was derived from the flux of the background star, which varied by about ~1% after the frame selection. The astrometry is available in Table A.1 but these results will not be analyzed.

In addition to the point sources, also scattered light from the near side of the gap edge of the CSD, which is illuminated by the central star, is visible in both datasets. Therefore, the scattering opacity of the dust grains in the disk surface is non-negligible even at these relatively long wavelengths. Interestingly, only the near side of the disk is visible which points to an asymmetry in the scattering phase function of the dust. This finding suggests that the dust grains are comparable to or larger than the observed wavelength (4–5 μm).

3.2 Color and magnitude comparison

3.2.1 Color–magnitude diagram

The absolute brightness of PDS 70 b in the new and archival data is derived from the calibrated magnitudes in Table 1 and the Gaia distance of 113.4 ± 0.5 pc (Gaia Collaboration 2018). We determined absolute magnitudes of 9.40 ± 0.27 mag and 8.52 ± 0.27 mag in the NB4.05 and M′ filters, respectively. The uncertainty on the parallax is negligible in the error budget. With the K1 and M′ magnitudes,we place PDS 70 b in a color–magnitude diagram to show its photometric characteristics with respect to field and low-gravity dwarfs (Dupuy & Liu 2012; Dupuy & Kraus 2013; Liu et al. 2016), other directly imaged planets planets and brown dwarfs (Marois et al. 2010; Bonnefoy et al. 2011, 2014; Ireland et al. 2011; Galicher et al. 2011; Bailey et al. 2013; Daemgen et al. 2017; Chauvin et al. 2017; Delorme et al. 2017; Rajan et al. 2017; Stolker et al. 2019, 2020), predictions by the AMES-Cond and AMES-Dusty evolutionary and atmospheric models (Chabrier et al. 2000; Allard et al. 2001; Baraffe et al. 2003), and blackbody spectra. The color–magnitude diagram was created with the species toolkit (Stolker et al. 2020) and is shown in Fig. 2. We note that the SPHERE K1 magnitude was adopted for PDS 70 b, HIP 65426 b, and HD 206893 B. Since the K1 filter is close to the central wavelength of a typical K band filter, the color between such filters is ≲0.1 mag, which has been quantified by considering all available DRIFT-PHOENIX spectra (Helling et al. 2008) with Teff in the range of 1000–2000 K. Such a color effect is small compared to the uncertainty on the KM′ color of these three objects.

The M′ flux of PDS 70 b is consistent with a mid to late M-type field dwarf and comparable in brightness to ROXs 42 Bb and GSC 06214 B, which are both young, planetary-mass companions. The latter is known to have a circumsubstellar disk (Bowler et al. 2011). Compared to the L-type directly imaged planets β Pic b and HIP 65426 b, PDS 70 b is brighter in M′ by ~1 mag. In addition to the absolute brightness, we derived a K1–M′ color of 2.86 ± 0.27 mag, which is significantly redder than any of the planets and brown dwarfs in the color–magnitude diagram. Specifically, PDS 70 b is about 2 mag redder than the young, planetary-mass companions and 1.3 mag redder than β Pic b. Most comparable in color are HIP 65426 b and HD 206893 B but the difference is still 0.4 mag and these objects are ≳1 mag fainter in M′. Interestingly, these are two of the reddest low-mass companions (Milli et al. 2017; Cheetham et al. 2019), with unusual M′ colors that might be caused by enhanced cloud densities close to their photosphere (Stolker et al. 2020).

The empirical comparison shows that PDS 70 b is brighter and/or redder than any of the other directly imaged planets. In addition, we compare the data with synthetic photometry from the AMES-Cond (cloudless) and AMES-Dusty (efficient mixing of dust grains) models, which have been computed from the isochrone data at an age of 5 Myr. The comparison in Fig. 2 shows that the observed flux in M′ is about 1.6 mag brighter than the AMES-Dusty predictions for an object of the same color, which would have a mass of 3–4 MJ. This flux difference corresponds to a factor of ≈2.1 in radius. In the model spectra, the dust causes a veiling of the molecular features and a shift of the photosphere to higher (cooler) altitudes. Consequently, the IR colors become redder and the M′ flux larger, in particular because of weaker CO absorption at 4.6 μm. While the radius had been calculated self-consistently in these models, the offset with the PDS 70 b magnitude may indicate that either the radius is larger than predicted and/or the atmosphere is even dustier than what is modeled.

A comparison of the photometric characteristics with the synthetic fluxes from a blackbody spectrum shows indeed that PDS 70 b is consistent with a blackbody temperature of ~1000 K and a radius of ≈5 RJ (see Sect. 4 for a more detailed estimation of the blackbody parameters). This is in tension with the predicted radii in the AMES-Dusty and AMES-Cond models, either of which have ≈ 1.4–1.8 RJ for 1–10 MJ at 5 Myr (see isochrones in Fig. 5). As was pointed out some time ago (Fortney et al. 2005; Marley et al. 2007), at these early ages (≲50–100 Myr) the (arbitrary) choice of the starting luminosity or radius in the models still matters a lot; put differently, the planet may have formed (much) later than the star. Whether considering a younger cooling age sufficiently alleviates the tension is discussed in Sect. 4.1.

PDS 70 b is located in the gap of a CSD and is actively accreting from its environment. Therefore, the planet might be partially obscured by (dusty) material in its vicinity, which is expected to attenuate the planet’s spectrum. To understand the impact of the dust on the color and magnitude of the object, we show reddening vectors in Fig. 2 for spherical grains with a homogeneous, crystalline enstatite composition (Scott & Duley 1996; Jaeger et al. 1998). The extinction cross sections were calculated with PyMieScatt (Sumlin et al. 2018) by assuming a log-normal size distribution with a geometric standard deviation of 2 (Ackerman & Marley 2001). For grains with a geometric mean radius of 0.1 μm, the extinction would cause a reddening of the KM′ color, whichwould result in an under- and overestimated blackbody temperature and radius, respectively. For 1 μm grains, the color is close to gray so potential extinction would cause an underestimation of the planet radius. The radius of PDS 70 b will be estimated and discussed in more detail in Sect. 4.1.

thumbnail Fig. 2

Color–magnitude diagram of MM $_{M'}$ versus KM′. The field objects are color-coded by M, L, and T spectral types (see discrete colorbar), the young and low-gravity objects are indicated with a gray square, and the directly imaged companions are labeled individually. PDS 70 b is highlighted with a red star. The blue and orange lines show the synthetic colors computed from the AMES-Cond and AMES-Dusty evolutionary tracks at an age of 5 Myr. Blackbody radiation curves are shown for 8 RJ, 4 RJ, and 2 RJ (black dashed lines, from top to bottom). The black arrows indicate the reddening by MgSiO3 grains with amean radius of 0.1 and 1 μm, and AM $A_{M'}$ of 0.05 and 2 mag, respectively.

thumbnail Fig. 3

Color–color diagram of HM′ versus L′–M′. The field objects are color-coded by M, L, and T spectral types (see discrete colorbar), the young and low-gravity dwarf objects are indicated with a gray square, and the directly imaged companions are labeled individually. PDS 70 b is highlighted with a red star. The blue and orange lines show the synthetic colors computed from the AMES-Cond and AMES-Dusty evolutionary tracks at an age of 5 Myr. The black dashed line shows the synthetic colors of a blackbody spectrum. The black arrows indicate the reddening by MgSiO3 grains with amean radius of 0.1 and 1 μm, and AM $A_{M'}$ of 0.05 and 5 mag, respectively.

3.2.2 Color–color diagram

While color–magnitude diagrams reveal trends related to the intrinsic brightness of an object, color–color diagrams are independentof the distance and radius. Therefore, they are a useful diagnostic for understanding correlations between colors which are related to the atmospheric characteristics. In the case of a forming planet, the interpretation is more complicated because the colors are also affected by the accretion luminosity and the presence of circumplanetary material. This may cause a reddening of the IR fluxes due to reprocessed radiation and extinction of the atmospheric flux.

The data from Fig. 2 are used together with available H (either broad- or narrowband) and L′ photometry of directly imaged companions (Biller et al. 2010; Ireland et al. 2011; Currie et al. 2012, 2013, 2014; Bonnefoy et al. 2014; Milli et al. 2017; Chauvin et al. 2017; Rajan et al. 2017; Keppler et al. 2018). We created a color–color diagram of HM′ versus L′ –M′ with species, which is displayed in Fig. 3. PDS 70 b is positioned in a red part of the diagram with a HM′ color of 4.44 ± 0.27 mag and a L′ –M′ color of 0.88 ± 0.35 mag. For HM′, we computed the synthetic MKO H band photometry (18.24 ± 0.04 mag) from the SPHERE spectrum of Müller et al. (2018), although the difference between the broadband H and narrowband H2 photometry is only ~0.1 mag. Both colorsare consistent with HD 206893B and the L′ –M′ color is also comparable to HIP 65426 b.

The color characteristics of PDS 70 b are clearly distinct from more evolved objects. Specifically, the sequence of field objects and cloudless atmosphere models show approximately gray colors at high temperatures, while toward lower temperatures the L′ –M′ color becomes bluer and then redder because of CO and CH4 absorption, respectively (see e.g., Stolker et al. 2020). Similarly, the increasing strength of H2O absorption in the H band causes a redder HM′ color toward lower temperatures. Interestingly, the H band spectrum of PDS 70 b shows only weak evidence of H2O absorption (Müller et al. 2018) so the origin of the very red HM′ color is presumably different.

Spectra of giant planets and brown dwarfs are usually not well described by blackbody emission due to molecular absorption which causes a strong variation in the photosphere temperature with wavelength. Indeed, the comparison of the colors in Fig. 3 shows that, for a given temperature, the blackbody colors are redder than the colors of M- and L-type field objects, as well as the predictions from the atmospheric models. Several of the directly imaged objects lie close to the blackbody curve but the uncertainties (on the M′ flux in particular) are large. The spectrum of a low-gravity atmosphere may indeed approach a blackbody spectrum if the quasi-continuum cross-sections of the dust grains dominate the atmospheric opacity.

3.3 The spectral energy distribution from 1 to 5 μm

3.3.1 Modeling approach of the SED

The obtained NB4.05 and M′ fluxes enable us to extend the SED of PDS 70 b into the 4–5 μm regime. To construct the SED, we adopted the Y to H band spectrum from Müller et al. (2018), which had been obtained with the integral field spectrograph (IFS) of SPHERE (Claudi et al. 2008), and also the NIRC2 L′ photometry from Wang et al. (2020). These data were combined with the new NB4.05 and M′ fluxes from this work, and the reanalyzed photometry of the NACO L′ and the SPHERE/IRDIS H23 and K12 data. For consistency in the SED analysis, we recalibrated the NIRC2 L′ magnitude with the stellar spectrum from Sect. 2.2 to 14.59 ± 0.18. In the K band, we only considered the reanalyzed SPHERE photometry while the SINFONI spectrum from Christiaens et al. (2019) was excluded due to a discrepancy in the calibrated fluxes between these datasets (see top panel in Fig. 4). Finally, we adopted a root mean square (rms) noise at 855 μm of 18 μJy beam−1 (i.e., 7.4 × 10−23 W m−2 μm−1) from the ALMA continuum imagery by Isella et al. (2019) as the approximate “forced photometry” (see discussion by Samland et al. 2017)at the position of PDS 70 b.

The SED is shown in the top panel of Fig. 4 across the 1–5 μm range. Apart from the potential broad, H2O absorption feature around 1.4 μm (Müller et al. 2018), we could not identify any obvious other molecular features (e.g., H2O, CH4, or CO) in the SED on visual inspection. Such absorption features might to be expected given the constraints on the temperature of the atmosphere, which is comparable to the HR 8799 planets (cf. Bonnefoy et al. 2016; Greenbaum et al. 2018; Mollière et al. 2020). We note that some of the smaller fluctuations in the SPHERE and SINFONI spectra may possibly be attributed to correlated noise. With this in mind, we attempt a simplified fitting approach by describing the spectrum with one or two blackbody components (see also Wang et al. 2020). A spectrum based on a single blackbody temperature may naturally describe a photosphere in which the dust opacity dominates over line absorption, with the temperature set either by the internal luminosity of the planet or by the accretion luminosity (see discussion in Sect. 4.1.3). Later on, a second temperature component is included to account for excess emission at thermal wavelengths (≳3 μm), for example due to to reprocessed radiation in a CPD.

The fit of the photometric and spectroscopic data was done with species. The toolkit uses the nested sampling implementation of MultiNest (Feroz & Hobson 2008) through the Python interface of PyMultiNest (Buchner et al. 2014). For the parameter estimation, we used a Gaussian log-likelihood function (see Greco & Brandt 2016), lnL(D|M)=12[(SIFSF)TC1(SIFSF)+i=19(dimi)2σi2],\begin{equation*}\ln \mathcal{L} (D|M)\,{=}\,{-}\frac{1}{2}\left[(\mathbf{S}_{\mathrm{IFS}}-\mathbf{F})^{T} \mathbf{C}^{-1} (\mathbf{S}_{\mathrm{IFS}}-\mathbf{F}) + \sum_{i\,{=}\,1}^{9} \frac{(d_i - m_i)^2}{\sigma_i^2}\right], \end{equation*}(1)

where D is the data, M the model, SIFS the IFS spectrum, F the model spectrum, C the (modeled) covariances for the IFS spectrum (see Eq. (2)), di the photometric flux for filter i, mi the synthetic flux from the blackbody spectrum, and σi the uncertainty on the flux di. The second term of Eq. (1) contains the sum over the nine photometric fluxes that were included in fit.

Spectra from integral field units are known to be affected by correlated noise (Greco & Brandt 2016). We therefore follow the approach byWang et al. (2020) and model the covariances of the SPHERE spectrum as a Gaussian process with a squared exponential kernel (Czekala et al. 2015; Wang et al. 2020), Cij=f2σiσjexp((λiλj)22l2)+(1f2)σiσjδij,\begin{equation*}C_{ij}\,{=}\,f^2 \sigma_i \sigma_j \exp{\left(-\frac{(\lambda_i - \lambda_j)^2}{2\ell^2}\right)} + (1 - f^2) \sigma_i \sigma_j \delta_{ij}, \end{equation*}(2)

where Cij is the covariance between wavelengths λi and λj, σi the total uncertainty on the flux of wavelength λi, f the relative amplitude of correlated noise with respect to the total uncertainty, and the correlation length. The correlation length and amplitude were fitted while adjusting the covariance matrix in the log-likelihood function (see Eq. (1)). Finally, each model spectrum was smoothed with a Gaussian filter to match the spectral resolution of the IFS data (R = 30) and resampled to the IFS wavelengths with SpectRes (Carnall 2017).

3.3.2 Parameter estimation and model evidence

The posterior distributions of the temperature, radius, and calibration parameters were sampled with 5000 live points and using uniform priors for all parameters except the correlation length. For the latter, we used a log-uniform sampling of the prior space. The marginalized distributions are shown in Figs. B.1 and B.2 for the cases of fitting one and two blackbody components, respectively. A comparison of the best-fit solution, randomly drawn spectra from the posterior, and the data are shown in Fig. 4.

When fitting one blackbody component, we constrained the temperature and radius of the photospheric region to 1193 ± 20 K and 3.0 ± 0.2 RJ, and we derived from this a luminosity of log(LL) = −3.79 ± 0.02. The overall spectral morphology appears well described by blackbody emission except for the deviation between the J and H bands. Also the 3–5 μm fluxes match reasonably well with the blackbody emission, thereby confirming the findings by Wang et al. (2020). Specifically, the NB4.05 and M′ fluxes deviate from the best-fit spectrum by 1σ. For the covariance model that describes the correlated noise in the IFS spectrum, we determined a length scale of ≈0.04 μm and a fractional amplitude of 0.54 ± 0.19. While the upper limit on the ALMA band 7 flux was included in the fit, its impact on the retrieved parameters is negligible because all the single-blackbody model spectra are below the rms noise at band 7.

Although the M′ flux only deviates by 1σ from the best-fit model, we also attempted a fit with two blackbody components to test if such a spectrum provides a better match at wavelengths ≳4 μm. A second blackbody component could for example describe the excess emission from a CPD, which will be discussed in Sect. 4.2. Here, we restricted the temperature and radius of the second component to values that are smaller and larger, respectively, than the first component by rejecting samples that did not met this condition. We also restricted the temperature prior of the second component to 0–600 K and the radius to 1–350 RJ, that is, extending up to ~0.1 times the Hill radius for a 1 MJ planet at 22 au (Tanigawa et al. 2012).

When fitting two blackbody components, the retrieved temperature (T1 = 1194 ± 20 K) and radius (R1 = 3.0 ± 0.2 RJ) of the first component are very similar to those from fitting a single blackbody component. For the second component, we constrained the temperature to T2 ≲ 256 K and the radius to R2 ≲ 245 RJ. The sparse wavelength coverage and large uncertainties at wavelengths longer than 4 μm leave a degeneracy between the temperature and radius of the second component (see Fig. B.2). Specifically, a large fraction of the samples is only driven by the upper limit at 855 μm while not fitting the M′ flux, since it is only a 1σ deviation from the first blackbody component. The posterior of T2 peaks toward 0 K, which is fully degenerate with the radius, R2, going to large values. Therefore, in the bottom panel of Fig. 4, we selected random samples with T2 > 100 K since the surface layers of a CPD are expected to be heated by accretion (e.g., Aoyama et al. 2018). When considering all posterior samples, we derived a luminosity ratio of log(L1/L2)=0.71.0+1.8$\log(L_1/L_2)\,{=}\,0.7^{+1.8}_{-1.0}$ for the two components (see Fig. B.2). Thus the luminosity of the second component would be about an order of magnitude smaller than the first component.

In addition to the parameter estimation, nested sampling has the advantage of providing the marginalized likelihood (i.e., the model evidence), which enables pair-wise model comparisons. The Bayes factor is used to quantify the evidence of favoring a certain model, and is given by the ratio of the evidence of two models in case the prior probability is the same for both models, B=Z(D|M0)Z(D|M1),\begin{equation*}B\,{=}\,\frac{\mathcal{Z}(D|M_0)}{\mathcal{Z}(D|M_1)}, \end{equation*}(3)

where Z(D|Mi)$\mathcal{Z}(D|M_i)$ is the evidence of data D given model Mi. In our case, the Bayes factor is calculated from the evidence ratio of fitting the SED with one or two blackbody components. We obtained a Bayes factor of 2.3, whichindicates weak evidence for favoring a model with one blackbody component when considering the Jeffreys’ scale (e.g., Trotta 2008).

thumbnail Fig. 4

Spectral energy distribution of PDS 70 b. Top and bottom panels: results from fitting one and two blackbody components, respectively (the flux units are different between the two panels). The photometric and spectroscopic data (colored markers) are shown in comparison with the best-fit blackbody spectrum (black line), and randomly drawn samples from the posterior distribution (gray lines). The residuals are shown relative to the data uncertainties. The Hα and Hβ (upper limit) fluxes (Hashimoto et al. 2020) are shown for reference but were not used in the fit.

4 Discussion

4.1 Implications from the luminosity and photospheric radius

Summarizing the fits to one or two blackbody component(s), the blackbody emission radius of the component peaking at smaller wavelength (i.e., L in Fig. B.1 and L1 in Fig. B.2) is Rphot ≈ 3.0 ± 0.2 RJ, while the corresponding luminosity is log(LSEDL) = −3.79 ± 0.02. Both numbers are comparable to the results of Wang et al. (2020) for one blackbody, two blackbodies (taking the luminosity only of the first), and even, as an extreme, their fit to the BT-Settl models. Thus, our Rphot and LSED seem robust. Here, we want to analyze systematically what they imply for the physical properties of PDS 70 b.

Since PDS 70 b is presumably still forming, one should think carefully about the evolutionary track models used for the analysis. An important aspect is the time evolution of the models. Cooling tracks need to assume an initial entropy for a given mass and thus, equivalently, an initial radius and luminosity (e.g., Arras & Bildsten 2006; Marleau & Cumming 2014). By definition, this state of the planet is “initial” with respect to the phase of pure cooling. It is set by the formation process and thus also referred to as the “post-formation” state (Marleau & Cumming 2014). As pointed out by Fortney et al. (2005) and Marley et al. (2007) and discussed for instance by Mordasini et al. (2012a, their Sect. 8.1), different formation scenarios will lead to different post-formation entropies. Therefore, the time label used in cooling track models is not guaranteed to be meaningful at early times. This holds in particular for a planet in the middle of formation, as PDS 70 b could conceivably be, but also if the current accretion rate is negligible such that the planet is evolving at essentially constant mass.

Predictions for the entropy of forming and “newborn” planets do exist (Mordasini et al. 2017) but here we take a more general approach. We ignore any time information and consider a grid of hydrostatic gas giant models labeled only by mass, Mp, and radius, Rp. This is possible because for a given atmospheric model, a non-irradiated gas giant planet has only two independent parameters, as discussedin Arras & Bildsten (2006), and also Marleau & Cumming (2014). In that latter work, these were Mp and the entropy, s, with Rp or the luminosity seen as functions of Mp and s, while here we consider Mp and Rp to be the two independent parameters. This allows us to drop the time label, thus circumventing the uncertainties about cold-, warm-,or hot-start conditions that are linked to the relevant physical processes (e.g., Mordasini 2013; Berardo et al. 2017; Marleau et al. 2017, 2019a). Therefore, our approach is independent of the entropy during and at the end of formation.

For the interpretation of the results, we assume that the observable bolometric luminosity of the one blackbody or both is in general the sum of three components: Ltot=Lint(Mp,Rp)+Lacc(Mp,Rp,M˙)+LCPD,\begin{equation*}{L_{\textrm{tot}}}\,{=}\,{L_{\textrm{int}}}({M_{\textrm{p}}},{R_{\textrm{p}}}) + {L_{\textrm{acc}}}({M_{\textrm{p}}},{R_{\textrm{p}}},{\dot{M}}) + {L_{\textrm{CPD}}}, \end{equation*}(4)

where Lint is the luminosity from the planet’s interior, possibly including some compression luminosity below the surface in the case that there is an accretion shock (Berardo et al. 2017); Lacc=ηGMpM˙(1/Rp1/Racc)${L_{\textrm{acc}}}\,{=}\,\eta G{M_{\textrm{p}}}{\dot{M}}\left(1/{R_{\textrm{p}}}-1/{R_{\textrm{acc}}}\right)$ is the luminosity from accretion at the surface of the planet; and LCPD is the sum of any thermal (e.g., Zhu 2015; Eisner 2015) and shock (e.g., Aoyama et al. 2018) emission from a possible CPD. This assumes that any shock luminosity from the planet surface or CPD is reprocessed and thermalized, with only a negligible fraction escaping at least as Hα; Aoyama et al. (2018) report that for a planetary shock, only a small fraction of Lacc goes into Hα. In the classical, highly simplified picture of material going directly from the CSD to the planet, the accretion radius is Racc ~ RHill (Bodenheimer et al. 2000), so that the 1∕Racc term is negligible compared to 1∕Rp. If however the gas releases part of its potential energy between the CSD and the CPD, the effective Racc would be closer to Rp but still possibly somewhat larger. Finally, we assume complete local radiative efficiency at the shock, η ≈ 1, following the results of Marleau et al. (2017, 2019a).

Therefore, in the following, we analyze what requiring Ltot = LSED implies. Here, we explore the case in which there is no CPD present or the emission from the CPD is negligible in the total luminosity budget, that is, LCPD = 0, motivated by the lack of evidence for a second blackbody component in the SED (see Sect. 3.3.2). We also assume that RaccRp. Alternative scenarios in which LCPD contributes to the bolometric luminosity will be discussed in Sect. 4.2. Finally, we note that Eq. (4) is valid under the assumption of isotropic radiation while, in particular for the shock (Lacc), this may not be accurate. We will deal with this in a crude manner below by considering also the case Lacc = 0.

We will use the BEX-Cond models (Bern EXoplanet cooling tracks; Marleau et al. 2019b), which graft the AMES-Cond atmospheres (Allard et al. 2001; Baraffe et al. 2003) onto the standard Bern planet structure code completo21(Mordasini et al. 2012a,b; Linder et al. 2019). The precise choice of atmospheric model, such as AMES-Cond, AMES-Dusty, or that of Burrows et al. (1997), does not influence much the mapping from Mp and Rp to Lint. In fact, AMES-Cond and AMES-Dusty both use exactly the same Rp(t) and Lint(t) tracks (these models differ only in the photometric fluxes), and apart for systematic shifts the results would be very similar for another set of atmospheres. The most important and generic feature of such models is that Lint increases with both Mp and Rp. In general, the functional form of this dependency Lint(Mp, Rp) is different than that of Lacc(Mp, Rp) ∝ MpRp.

4.1.1 Constraints from the luminosity on the planetary radius and mass

We explore first what only the derived luminosity LSED implies for the physical radius of PDS 70 b, defined as the (very nearly) hydrostatic structure terminating in general at the photosphere or at the shock location. We will return to Rphot only in Sects. 4.1.2 and 4.1.3. As a limiting case, we consider at first Lacc = 0 in Eq. (4), that is, we assume that by some geometric effects (for instance magnetospheric accretion at the planet’s poles, away from the observer) most of the (reprocessed) accretion luminosity is not reaching an observer on Earth, and therefore that the observed luminosity is coming only from the photosphere while LCPD is assumed to be zero. The left panel of Fig. 5 shows the Mp and Rp combinations consistent with this extreme assumption, that is, the models that have Ltot = Lint equal to LSED. At any mass there is an Rp such that Ltot = LSED, with Rp ≈ 2.8 RJ at most, down to 1.5 RJ at 10 MJ. Formally, there is no upper mass limit.

The right panel of Fig. 5 assumes instead = min = 5 × 10−7 MJ yr−1 in Eq. (4), corresponding to the lower limit on the accretion rate derived by Hashimoto et al. (2020). Now, only small masses Mp ≲1.5 MJ are allowed since Lacc raises everywhere Ltot significantly. For example, at Mp = 5 MJ, LSED would need to be higher than derived by at least ≈0.5 dex (≈25σLSED$\sigma_{{L_{\textrm{SED}}}}$) for there to be a matching luminosity. The discrepancy is larger for larger Mp values. Where Ltot can be matched, however, the possible radii3 range from Rp ≈ 1.1 to ≈ 2.8 RJ. This is thus the physical radius of PDS 70 b implied by LSED alone (i.e., ignoring Rphot), assuming = 5 × 10−7 MJ yr−1, and the corresponding mass is Mp ≈ 0.5–1.5 MJ.

As a rough check, we inspected the Bern population synthesis4, both from Generation Ib (Mordasini et al. 2012a,b, 2017) and from the newest, Generation III (Emsenhuber et al. 2020a,b) to see whether this combinationof (, Mp, Rp) is met. We find that, not considering the time at which this happens in the population synthesis, planets accreting at min have Rp ≈ 1.3–1.7 RJ for Mp ≈ 0.5–2 MJ, reaching up to Rp ≈ 2.5 RJ down to Mp ≈ 0.3 MJ. Since this range of Rp is within our allowed range Rp ≈ 1.1–2.8 RJ, the would be consistent with these formation models.

As mentioned, the value from Hashimoto et al. (2020) is a lower limit. Already at ≈ 3min, the implied mass from the structure models (dotted white line in the right panel of Fig. 5) is5 Mp≲0.6 MJ. This mass might seem small but at least in the Bern population synthesis, there are planets in the corresponding region of (, Mp, Rp), again not taking time into account. Thus such low-mass solutions might be possible. On the other hand, if the lower limit on is overestimated, then the derived Mp is underestimated (see the = 0.3min case in Fig. 5) because LtotLaccMp. Hence, to keep the luminosity constant (i.e., equal to LSED), a smaller is compensated by an increase in Mp. Nevertheless, we need to see whether the derived Mp range matches other constraints.

One constraint is the presence of a gap. From Kanagawa et al. (2016) a suitable combination of disk parameters (scale height and viscosity) could lead to a gap even at low masses. For example, with an aspect ratio of Hprp = 0.067 (Bae et al. 2019) at the separation of the planet (rp ≈ 22 au), an estimated gap width of 20 au, and a turbulence parameter of α = 10−4, given the evidence for weak turbulence in protoplanetary disks (e.g., Flaherty et al. 2020), the relation from Kanagawa et al. (2016) implies a planet mass Mp ≈ 1 MJ. For this combination of parameters, the derived mass is likely an upper limit since the gap in the PDS 70 disk is opened by the combined effect of two planets. In any case, in a first approximation the low Mp value inferred from the luminosity seems compatible with the presence of a gap. More detailed, radiation-hydrodynamical modeling of the disk would clearly be warranted.

A second aspect concerns the orbital stability of the system. Bae et al. (2019) studied the dynamics of the PDS 70 system by fixing Mp = 5 MJ for planet b and varying the mass of c from 2.5 to 10 MJ. They concluded that the orbits are likely in a 2:1 mean-motion resonance (as had been suggested by Haffert et al. 2019) and can remain dynamically stable over millions of years. It would be interesting to repeat their simulations with a lower mass Mp ≈ 1 MJ for planet b, possibly also considering a lower mass for PDS 70 c. A low Mp value would be in agreement with the N-body simulations by Mesa et al. (2019), who showed that the two-planet system would be stable with masses of 2 MJ, whereas dynamical perturbations occurred in their simulations with higher-mass planets.

Coming back to the luminosity constraint, we compare LSED to the AMES (hot-start) isochrones, discussing the validity of this approach afterward. Figure 5 shows that LSED intersects the 5-Myr isochrone (roughly the age of the star) at Rp = 1.6 RJ when not considering a contribution of Lacc in Ltot (left panel), and at Rp ≈ 1.4–1.5 RJ for ≈ (0.3–3)min (right panel). Taken at face value, the corresponding masses are Mp ≈ 6.5 MJ for Lacc = 0 and again Mp ≈ 0.5–1.5 MJ for ≈ (3–1)min.

However, several caveats apply. One is that the AMES models were made for isolated planets, whereas during formation there can be a spread of at the very least 0.3 dex in Lint at a given mass at 5 Myr for what are effectively hot starts (see Figs. 2 and 4 of Mordasini et al. 2017). This will affect the derived radius. Another concern is that there might be a formation delay of perhaps a few Myr, which would be significant at this age (Fortney et al. 2005). In the case of HIP 65426, which has a mass of 2 M (Chauvin et al. 2017), Marleau et al. (2019b) estimated roughly a formation time near 2 Myr and argued that this should increase with lower stellar mass (PDS 70 has a mass of 0.8 M; Keppler et al. 2018). Finally, the true shape of the physical isochrone could be different than in the AMES track, which was not guided by a formation model. In particular the post-formation radius as a function of mass could conceivably be non-monotonic, allowing for several solutions to Ltot(t = 5 Myr) = LSED.

We show as an extreme comparison the 1-Myr hot-start isochrone in the left panel of Fig. 5. This would imply a somewhat larger radius Rp ≈ 1.8 RJ. For Lacc = 0, the mass would be clearly smaller, with Mp ≈ 3 MJ, whereas for = min the mass would be similar to the 5 Myr isochrone (not shown explicitly in the plot). We note that the AMES isochrone at the maximal age (the system age) does provide an upper mass limit within the hot-start assumption; younger ages necessarily imply smaller masses, as Fig. 5 makes clear. However, the initial radius could be smaller. While extreme cold starts at the Marley et al. (2007) are disfavored (e.g., Berardo et al. 2017; Mordasini et al. 2017; Snellen & Brown 2018; Wang et al. 2018; Marleau et al. 2019a), warm starts seem a realistic possibility. In this case, a 5-Myr isochrone would match the luminosity at a smaller radius and larger mass—thus the mass upper limit from the hot start is in fact a “lower upper limit,” meaning it is not informative.

Interestingly, the mass that we derived from the luminosity LSED and the adopted lower limit on the accretion rate (i.e., the right panel in Fig. 5) appears lower then what has been inferred in previous studies. Keppler et al. (2018) compared the H and L′ band magnitudes with predictions by evolutionary models and estimated a mass of 5–9 MJ and 12–14 MJ with a hot-start and warm-start formation, respectively. Similarly, the estimates by Müller et al. (2018) also assumed that the NIR fluxes trace directly the planet atmosphere. The authors determined a mass of 2–17 MJ by fitting the SED with atmospheric model spectra. One difference in our analysis is that it is based on the bolometric luminosity, but the key point is that it takes the accretion luminosity into account. Therefore, it is not surprising that we derive a different (i.e., lower) planet mass. We note that our quoted mass error bars are smaller both in relative and absolute terms with respect to previous studies, but our error bars do not reflect the (large) uncertainty on the accretion rate.

More recently, Wang et al. (2020) estimated a mass for PDS 70 b of 2–4 MJ from the bolometric luminosity by using the evolutionary models of Ginzburg & Chiang (2019) and assuming a system age of 5.4 Myr. Our mass constraint (Mp ≈ 0.5–1.5 MJ with = min) is somewhat comparable with the findings by Wang et al. (2020), indicating a relatively low mass compared to earlier estimates. However, we want to stress again that the low planet mass that we estimated hinges on the adopted accretion rate. If the accretion rate is overestimated then the mass is underestimated, as can be seen from the × 0.3 example in Fig. 5. The mass of planet b that was derived by Hashimoto et al. (2020) from the Hα flux is significantly larger (~12 MJ) than our value. However, the authors noted that the line profile was not resolved, hence only an upper limit on the free-fall velocity could be determined (v0 = 144 km s−1). Since the planet mass scales quadratically with the velocity (see Eq. (3) in Hashimoto et al. 2020), it would require a factor of ≈ 3 smaller velocity to lower the estimated mass from 12 to 1.5 MJ. We note that a velocity of v0 = 48 km s−1 would still be twice as large as the minimum velocity that is required to produce Hα emission (i.e., v0,min ≈ 25 km s−1; see Fig. 6 by Aoyama et al. 2018).

thumbnail Fig. 5

Total luminosity from standard models of isolated planets as a function of mass and radius, but withalso the accretion luminosity from the planet surface shock (Eq. (4) with LCPD = 0). The thick black contour and gray shade highlight the measured luminosity and 5σ uncertainty (see main text for details), log(LSEDL) = −3.79 ± 0.02, and the thin black contours are steps of 1.0 dex relative to LSED. The accretion rate was set to = 0 in the left panel and = 5 × 10−7 MJ yr−1 (Hashimoto et al. 2020) in the right panel. The inferred photospheric radius and 1σ uncertainty, Rphot ± σ, are shown with a horizontally dashed line and gray shaded area. The light yellow curves show as reference the mass–radius predictions by the AMES (Cond/Dusty) structure model at 5 Myr (solid) and 1 Myr (dashed). The dotted white lines in the right panel show Ltot = LSED with scaled by a factor of 0.3 and 3. The crosshatched regions indicate parts of the parameter space for which no predictions are available from the structure model: because of electron degeneracy pressure at small radii, and because no stable hydrostatic structure exists at large radii.

4.1.2 Comparing the planetary and photospheric radii

How does the planet radius discussed so far compare to the derived photospheric radius derived above, Rphot ≈ 3.0 RJ? In both panels of Fig. 5, and in particular for min, the models with Ltot = LSED all have Rp < Rphot, with a substantial difference between the two. Put differently, there is no mass predicted by the structure model for which the radius is equal to Rphot and the total luminosity is equal to LSED simultaneously. A non-zero contribution from an accretion luminosity only exacerbates the tension.

For a range of masses between 1 and 10 MJ, the discrepancy ΔRRphotRp is typically at least 0.7–1.5 RJ, which is 3.5σRphot$3.5\sigma_{{R_{\textrm{phot}}}}$ at Mp = 1 MJ and 7.5σRphot$7.5\sigma_{{R_{\textrm{phot}}}}$ at Mp = 10 MJ. It does decrease toward low masses for both Lacc = 0 and ≠0, such that formally there is a narrow match within the 1–2σ regions of Rphot and LSED. However, this is at the maximum radius possible for a convective hydrostatic planet and thus seems unlikely, especially given that PDS 70 b probably has evolved at least for a short time, even if not the full age of the system (5.4 Myr). In short, there is in fact no satisfying solution within one or two σRphot$\sigma_{{R_{\textrm{phot}}}}$.

If one takes the AMES isochrone at 5 Myr, the discrepancy between the physical and photospheric radii is at the very least (taking Lacc = 0) ΔR = 1.4 RJ, or 7σRphot$\,{\approx}\,7\sigma_{{R_{\textrm{phot}}}}$. With the extreme case of a 4.4-Myr formation delay, and thus the 1-Myr isochrone, the difference is still 6σRphot$6\sigma_{{R_{\textrm{phot}}}}$. In any case, we argued that the AMES cooling models are possibly not directly appropriate for (maybe still forming) young planets.

There are four non-mutually exclusive possible implications from this discrepancy between the inferred physical and photospheric radii:

  • (i)

    LSED is underestimated. This could be the case if Lacc dominates LSED (after reprocessing) and is emitted anisotropically, which is conceivable. A dominating Lacc in turn seems plausible if is higher than min from Hashimoto et al. (2020). Alternatively, or in addition, extinction in the system may lead not only to radiation (Lint and/or Lacc) being shifted to longer wavelengths, but also to it being re-emitted away from the observer. This could possibly come from non-isotropic scattering by dust grains in the CSD (through the upper layers of which we are observing the PDS 70 b region) or in a CPD. In any case, such geometry effects would let LSED represent only a fraction of Ltot, allowing in principle for mass–radius solutions given classical planet structure models.

  • (ii)

    Rphot is overestimated. Assuming that the data constrain the shape and thus the approximate Teff of the spectrum, this is equivalent to (i). The derived Teff and Rphot from fitting synthetic spectra are sometimes correlated, therefore a different model spectrum may give a larger Teff and a smaller Rphot, withoutchanging LSED much. For example, a decrease in the radius by 50% would correspond to an increase in the temperature by ≈ 40%, such that the luminosity remains constant. Alternatively, extinction by small dust grains could have altered the SED, possibly mimicking a larger radius and smaller temperature (see Fig. 2).

  • (iii)

    The structure models (classical, non-accreting gas giants that turn out to be fully convective) do not apply and Lint should be smaller at a given (Mp, Rp), or equivalently Rp should be larger at a given (Mp, Lint), assuming Lint still increases with both Mp and Rp. Recent modeling work by Berardo et al. (2017) suggests that this might hold, at least qualitatively, but the effect might not be large enough.

  • (iv)

    Rp is the physical radius of the planet, implying there is Rosseland-mean optically thick material between Rp and Rphot.

The last possibility is a particularly interesting one that we consider in more detail in the next subsection.

4.1.3 Constraints on the vicinity of the planet?

We will assume that there is optically thick material between Rp and Rphot. This material could be flowing onto the planet or be in some layer of the CPD. However, given that CPDs are thought to be a fraction of the size of the Hill sphere (e.g., Lubow et al. 1999) and that RphotRHill ~ 3500 RJ for a ~1 MJ planet at 22 au, this is most likely material flowing onto the planet (see Sect. 4.2 for a more detailed discussion on the presence/absence of a CPD).

In principle, the extinction could be due to gas or to dust opacity. However, the absence of strong molecular features (as argued in Sect. 3.3.1), contrary to what would be expected from gas at T ~ 1000 K, suggests that the opacity is grayer and dust-dominated (Wang et al. 2018). One can estimate whether the(possibly high) temperatures near the planet would allow for dust to exist within Rphot ≈ 3.0 RJ. From Isella & Natta (2005), dust is destroyed at Tdest ≈ 1280–1340 K at ρ ~ 10−10–10−9 g cm−3 (see below). Assuming that the luminosity is approximately constant in the accretion flow onto the planet (Marleau et al. 2017, 2019a), Teff = 1193 K at Rphot = 3.0 RJ implies a local Teff, loc = 1307 K at a radius 2.5 RJ (Teff = 1687 K at 1.5 RJ). The temperature of the gas and dust, in turn, is given by solving the implicit approximate equation TTeff,loc/41/4×(1+1.5κRρRp)1/4$T\,{\approx}\,{T_{\textrm{eff,\,loc}}}/4^{1/4}\,{\times}\,(1&#x002B;1.5{\kappa_{\textrm{R}}}\rho{R_{\textrm{p}}})^{1/4}$, where κR (T) is the Rosseland mean opacity (see Eq. (32) of Marleau et al. 2019a). Given the numerical values, we estimate that at both positions the dust could be partially destroyed. This suggests that the dust destruction, which is strongly sensitive to the temperature(Bell & Lin 1994; Semenov et al. 2003), could occur over a non-negligible spatial scale comparable to ΔR. This effect is seen as the temperature plateau in Fig. 9 of Marleau et al. (2019a). Geometrical effects in the accretion flow will affect the details but in an average sense, the region between the planet radius and the photosphere could well be partially filled with dusty material, with full abundance near Rphot decreasing smoothly toward Rp.

For the radiation to be reprocessed between Rp and Rphot, there must be at least a few Rosseland-mean optical depths between the two radii, that is, ΔτR ~ ⟨κRρ⟩(RphotRp) > 1 (i.e., the infrared extinction must be AIR ≳ 1 mag). For a filling factor ffill, the typical gas preshock density is ρ=M˙/(4πRp2ffillvff)~$\rho\,{=}\,{\dot{M}}/(4\pi{R_{\textrm{p}}}^2{f_{\textrm{fill}}}{v_{\textrm{ff}}})\,{\sim}$ 10−10 g cm−3, with the preshock free-fall velocity given by vff22GMp/Rp${v_{\textrm{ff}}}^2\,{\approx}\, 2G{M_{\textrm{p}}}/{R_{\textrm{p}}}$ (e.g., Zhu 2015). With ΔR ~ 1 RJ, the requirement ΔτR > 1 implies κR ≳ 1 cm2 g−1 as a conservative lower limit. This is the opacity per unit gas mass, that is, the opacity per unit dust mass κR, • times the dust-to-gas ratio fd/g. At these temperatures near or below 2000 K, the gas opacity is κR ≲ 10−2 cm2 g−1 (Malygin et al. 2014), which implies that dust dominates the total opacity budget, so that the requirement would be κR, •≳ 100∕(fd/g∕0.01) cm2 g−1. For the canonical fd/g = 0.01, this needed opacity is in line with the calculation of Semenov et al. (2003) and seems in general reasonable given the uncertainties about the exact dust composition, porosity, non-sphericity, material properties, etc. If the accreting gas comes from high latitudes in the CSD (Tanigawa et al. 2012; Morbidelli et al. 2014; Teague et al. 2019), the settling of the dust to the midplane could imply a lower fd/g in the accretion flow (Uyama et al. 2017), perhaps by as much as a few orders of magnitude. Even in this case, however, the required κR, • seems consistent with predictions from Woitke et al. (2016) for an appropriate size distribution and material properties for the dust. Finally, this rough estimate assumes a spherically symmetric accretion flow; if ffill < 1 in the accretion flow and this concentration of matter is along the line of sight, the required minimum opacity would be lower, proportionally to ffill, and thus easier to reach. Therefore, altogether, a photospheric radius that is larger than the planet radius could be explained by dusty material in the vicinity of PDS 70 b.

4.2 Mid-infrared excess from a circumplanetary disk?

The formation of a giant planet is characterized by several distinct phases of growth. At an early stage, the accretion flow from the CSD feeds directly the atmosphere of the object, through spherical accretion of gas and solids entering the planet’s Hill sphere (e.g., Pollack et al. 1996; Cimerman et al. 2017). As the planet grows further, the gaseous envelope may collapse, thereby triggering a runaway accretion and the potential formation of a CPD (e.g., Canup & Ward 2002; D’Angelo et al. 2003). Hydrodynamical simulations have indeed shown that a CPD can remain, spanning a fraction of the planet’s Hill radius (e.g., Lubow et al. 1999; Ayliffe & Bate 2012; Tanigawa et al. 2012; Szulágyi et al. 2016). The disk will act as important mediator for channeling the infalling gas and dust toward the planet (e.g., Tanigawa et al. 2012) and the accretion onto the planet–disk system may leave a strong imprint on the bolometric luminosity (Papaloizou & Nelson 2005).

In Sect. 4.1, we assumed that the SED luminosity reflected the planet’s interior and accretion luminosity, while ignoring a contribution from a CPD. The analysis in Sect. 3.3 revealed indeed weak evidence that the SED is better described by one blackbody component instead of two. This is mainly because the second component is only constrained by the 1σ deviation of the M′ flux and the non-detection with ALMA at the expected position of PDS 70 b. Therefore, an alternative interpretation based on the deviation of the M′ flux will be very speculative. Nonetheless, we will briefly discuss our findings in the context of a CPD that could be present.

If there is no excess emission at MIR wavelengths, the L′, NB4.05, and M′ fluxes trace the same photospheric region as the NIR part of the SED. In that case, the SED is described by a single temperature and radius, which can be characterized by an extended dusty environment, as discussed in Sect. 4.1. A non-detection of a CPD in M′ and with ALMA at 855 μm may indicate that the CPD is either very faint (e.g., low in temperature and/or mass), that the physical conditions near the planet do not allow (yet) the formation of a CPD, or that the CPD may have already been dispersed. This finding, combined with the constraint on the mass of PDS 70 b from Sect. 4.1.1 (~1 MJ) may guide the calibration of CPD models.

Alternatively, we speculate that the slight excess emission in M′ could trace a second component from a cooler (≲256 K) and more extended region (≲245 RJ) that is associated with a CPD. This could either be thermal emission coming directly from the disk or reprocessed emission from the accretion shock on the surface of the disk, as given for example by Aoyama et al. (2018). From the retrieved temperatures and radii of the two blackbody components, we derived that the luminosity of the cooler component, L2, is approximately an order of magnitude smaller than that of the first component, L1 (although they could be comparable within 1σL1/L2$\sigma_{L_1/L_2}$ since log(L1/L2)=0.71.0+1.8$\log(L_1/L_2)\,{=}\,0.7^{&#x002B;1.8}_{-1.0}$; see Fig. B.2). If we assume that the CPD is heated by the luminosity of the planet, this may indicate that ≲10% of the planet flux is reprocessed by the CPD. Here, the percentage of reprocessed emission is an upper limit since the CPD is also expected to be heated by accretion from the CSD and/or viscous heating. Such a process may in fact dominate the luminosity budget of the CPD.

Previously, Christiaens et al. (2019) suggested that part of the K band flux originates from a CPD, since the considered atmosphere models could not explain the absolute flux and slope of the SINFONI spectrum. Although there is a discrepancy between SINFONI and SPHERE K band fluxes, Fig. 4 shows that the SPHERE photometry is consistent with a blackbody spectrum (see also Wang et al. 2020), therefore possibly not requiring excess flux from a CPD at these wavelengths. However, this needs to be confirmed. Instead, we identified very marginal excess emission in the M′ band, but we stress that the result is not significant. More precise photometry at 4–5 μm is required to constrain the circumplanetary characteristics of PDS 70 b, for example with the aperture masking interferometry (AMI) mode of the NIRISS instrument (Artigau et al. 2014) on board the James Webb Space Telescope (JWST).

The spectral appearance of a forming planet and the disk surrounding it will deviate from that of an isolated object with atmospheric emission alone. Predictions by Zhu (2015), based on a simplified, steady-state disk model, showed thatan accreting CPD can be brighter at near- and mid-IR wavelengths than the planet itself if Mp is sufficiently large. Therefore, the peak in the observed SED may solely trace the hottest region of the CPD instead of the planet atmosphere.This would imply that the actual planet is not visible and we mainly detect the (reprocessed) luminosity from the disk, that is, LSEDLCPD.

Considering a 1 MJ planet and = 5 × 10−7 MJ yr−1 (Hashimoto et al. 2020), the predictions in Fig. 1 by Zhu (2015) show that the NIR fluxes are dominated by the emission from the planet atmosphere instead of the viscous heating in the CPD (i.e., when Mp ~ 10−7106 MJ2yr1$10^{-6}~M_{\textrm{J}}^2\,\textrm{yr}^{-1}$ and Teff ~ 1000 K), unless the inner radius of the CPD is very small (Rin = 1 RJ) and/or MpM˙106 MJ2yr1${M_{\textrm{p}}}{\dot{M}} \gtrsim 10^{-6}~M_{\textrm{J}}^2\,\textrm{yr}^{-1}$. In the case the SED mostly traces emission from the CPD, the spectral slope is expected to be less steep at longer wavelengths due to the radial temperature gradient in the disk. This contrasts with the observed SED, which is consistent with a single blackbody (see Sect. 3.3), therefore pointing to a photospheric region that is characterized by a single temperature.

Since the adopted accretion rate from Hashimoto et al. (2020) is a lower limit, we applied the fitting procedure from Sect. 3.3.1 on the predicted magnitudes by Zhu (2015) to test what blackbody temperature and radius would be retrieved if the accretion rate is larger and the SED only traces CPD emission. For this, we considered the full-disk case with Rin = 2 RJ and MpM˙=105 MJ2yr1${M_{\textrm{p}}}{\dot{M}}\,{=}\,10^{-5}~M_{\textrm{J}}^2\,\textrm{yr}^{-1}$. We adopted the J- to M′-band magnitudes and added arbitrary error bars of 0.1 mag. When fitting a single blackbody, we retrieved Teff = 994 ± 15 K and R = 11.6 ± 0.6 RJ. The flux density peaks at ~3 μm, which is similar to the SED of PDS 70 b (see top panel in Fig. 4), and the M′ flux from the CPD model shows a 10–20% excess with respect to the best-fit blackbody spectrum (due to the temperature gradient in the disk). Interestingly, while the temperature is somewhat comparable to the photospheric temperature of PDS 70 b (Teff = 1193 ± 20 K), the retrieved radius from the predicted CPD fluxes is clearly larger than the photospheric radius of PDS 70 b (R = 3.0 ± 0.2 RJ). This brief assessment may suggest that both the accretion rate and photospheric radius of PDS 70 b are too small to interpret the SED as LSEDLCPD, so the photosphere traces presumably a more compact, dusty environment around the planet instead of a CPD. However, a more detailed analysis would be required to confirm this.

Apart from a luminosity contribution by a viscously heated CPD, the accretion flow and shock (on the planet surface and/or disk) may further alter the energy distribution. For example magnetospheric accretion from the disk onto the planet could also heat the photosphere of the planet, thereby enhancing the flux at shorter wavelengths (Zhu 2015). The importance of such accretion processes remain poorly constrained and can additionally be variable and subject to outbursts (Lubow & Martin 2012; Brittain et al. 2020).

5 Summary and conclusions

We have reported on the first detection of PDS 70 b at 4–5 μm. We used high-resolution observations with NACO at the VLT to image the forming planet with the NB4.05 (Brα) and M′ filters. PDS 70 c is tentatively recovered in NB4.05 and the near side of the gap edge of the CSD is detected in scattered light. We have also reanalyzed the photometry of PDS 70 b from archival SPHERE H23 and K12, and NACO L′ imaging data.

The absolute M′ flux of PDS 70 b is compatible with a late M-type dwarf, and the young, planetary-mass objects ROXs 42 Bb and GSC 06214 B. The NIR–M′ colors, on the other hand, are redder than any of the known directly imaged planets and most comparable to the dusty, L-type companions HD 206893 B and HIP 65426 b. While the M′ magnitude and related colors are unusual compared to other directly imaged planets, they are consistent with blackbody emission from an extended region that is several times the radius of Jupiter.

With the new NB4.05 and M′ photometry, we modeled the available SED data (including a SPHERE/IFS spectrum) by assuming a blackbody and derived a photospheric temperature of Teff = 1193 ± 20 K and radius of Rphot = 3.0 ± 0.2 RJ, which is consistent with the blackbody analysis of the 1–4 μm SED by Wang et al. (2020). Apart from small-scale deviations (partially due to expected correlated noise in the NIR spectra) and the tentative H2O feature at 1.4 μm, the photometric and spectroscopic data appear to be well described by a single blackbody temperature and radius. From the sampled posterior distributions, we derived a bolometric luminosity of log (LL) = −3.79 ± 0.02.

The derived luminosity and photospheric radius enabled us to place constraints on the planetary radius and mass of PDS 70 b. We used standard models for isolated gas giant planets to infer the mass–radius solutions corresponding to the measured luminosity, while taking into account the accretion luminosity. The time-independent approach of the analysis makes it unaffected by the uncertain cooling time of the object. Here we summarize the main findings and conclusions from this analysis:

  • (i)

    In the limiting case that Lacc = 0 (e.g., due to a geometric effect), there are solutions of the radius for all considered masses (up to 10 MJ), but always smaller than Rphot.

  • (ii)

    When including Lacc in the luminosity budget (based on the Hashimoto et al. 2020 estimate of the accretion rate), only masses up to 1.5 MJ have solutions for which the observed luminosity is equal to the combination of the intrinsic and accretion luminosity.

  • (iii)

    Considering these two cases, we constrain the mass of PDS 70 b to Mp ≈ 0.5–1.5 MJ and the physical radius to Rp ≈ 1–2.5 RJ. This is consistent with predictions from population synthesis models of forming planets and an approximate estimate based on the gap width.

  • (iv)

    The discrepancy between the photospheric and planetary radius could imply that the planet is enshrouded by a dusty, extended environment, which is consistent with the approximate blackbody spectrum and the dearth of strong molecular features.

  • (v)

    The derived photospheric radius is orders of magnitude smaller than the planet’s Hill radius. In the case of a dusty envelope, this indicates that the extended region is actively replenished by dust that is coupled to the gaseous accretion flow from the CSD (see also Wang et al. 2020).

  • (vi)

    Alternatively, the discrepancy may indicate that the actual luminosity is larger than the observed luminosity, for example due to anisotropic emission or scattering, extinction, or that the structure models may not apply because PDS 70 b is still forming.

The M′ flux shows a slight deviation from the best-fit results when considering a single blackbody temperature. We modeled the MIR excess with a second blackbody component and obtained an approximate upper limit on the temperature and radius of potential emission from a CPD, Teff ≲ 256 K and R ≲245 RJ, but the Bayes factor indicates weak evidence that the data is better described by a model with a single blackbody component. Higher precision photometry at MIR wavelengths is required to place stronger constraints on potential emission from a CPD, for example with the improved 4–5 μm imaging capabilities of VLT/ERIS, the AMI mode of NIRISS instrument on JWST, and in the further future M′ and N band photometry with ELT/METIS.

Acknowledgements

We would like to thank André Müller, Dino Mesa, and Valentin Christiaens for kindly sharing their SPHERE and SINFONI spectra and we thank Yuhiko Aoyama, Timothy Gebhard, Greta Guidi, and Jason Wang and the ExoGRAVITY team for clarifying discussions. We also thank the referee whose constructive comments improved the quality of this manuscript. T.S. acknowledges the support from the ETH Zurich Postdoctoral Fellowship Program. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1) and from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. T.S., G.C., and S.P.Q. thank the Swiss National Science Foundation for financial support under grant number 200021_169131. P.M. acknowledges support from the European Research Council under the European Union’s Horizon 2020 research and innovation program under grant agreement No. 832428. K.O.T acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 679633; Exo-Atmos). Part of this work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. S.P.Q. acknowledges the financial support of the SNSF.

Appendix A Astrometric calibration

In this appendix, we provide an overview of the calibrated astrometry. The final separation is calculated by adding the bias offset and combining the two error components in quadrature. The final position angle is calculated by adding the bias and true north offset and combining the three error components in quadrature. For NACO, we adopted a plate scale of 27.2 mas pixel−1 and true north of − 0. °44 ± 0. °1 from Cheetham et al. (2019). For SPHERE/IRDIS, we adopted from Maire et al. (2016) a plate scale of 12.25 and 12.26 mas pixel−1 for the H23 and K12 filters, respectively, and true north of − 1. °75 ± 0. °08.

Table A.1

Astrometry and error budget.

Appendix B Posterior distributions

Figures B.1 and B.2 show the 1D and 2D projections of the posterior samples from fitting the photometricand spectroscopic data of PDS 70 b with a model spectrum consisting of one and two blackbody components, respectively.Throughout this work, we have used the median of each parameter as the best-fit value, and the 16th and 84th percentiles as the 1σ uncertainties. For the second blackbody component, we have quoted the 84th percentile as the upper limit on the temperatureand radius. For a single blackbody component, the fitted (photospheric) temperature and radius are Teff and R, while for two blackbody components,these are given as T1 and R1, T2 and R2 for the first andsecond component. For the SPHERE/IFS spectrum, we have fitted the logarithm of the correlation length, log SPHERE and fractional amplitude of the correlated noise, fSPHERE (see Sect. 3.3.1 for details).

thumbnail Fig. B.1

Posterior distributions from fitting a single Planck function to the SED of PDS 70 b. The 1D marginalized distributions are shown in the diagonal panels and the 2D parameter projections in the off-axis panels. The bolometric luminosity, log LL, has been calculated from the posterior samples of Teff and R.

thumbnail Fig. B.2

Posterior distributions from fitting a combination of two Planck functions to the SED of PDS 70 b. Further details are provided in the caption of Fig. B.1.

References

  1. Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872 [NASA ADS] [CrossRef] [Google Scholar]
  2. Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357 [NASA ADS] [CrossRef] [Google Scholar]
  3. Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [NASA ADS] [CrossRef] [Google Scholar]
  4. Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aoyama, Y., & Ikoma, M. 2019, ApJ, 885, L29 [NASA ADS] [CrossRef] [Google Scholar]
  6. Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84 [NASA ADS] [CrossRef] [Google Scholar]
  7. Arras, P., & Bildsten, L. 2006, ApJ, 650, 394 [NASA ADS] [CrossRef] [Google Scholar]
  8. Artigau, É., Sivaramakrishnan, A., Greenbaum, A. Z., et al. 2014, SPIE Conf. Ser., 9143, 914340 [Google Scholar]
  9. Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44 [NASA ADS] [CrossRef] [Google Scholar]
  10. Ayliffe, B. A., & Bate, M. R. 2012, MNRAS, 427, 2597 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJ, 884, L41 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bailey, V., Hinz, P. M., Currie, T., et al. 2013, ApJ, 767, 31 [NASA ADS] [CrossRef] [Google Scholar]
  13. Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
  15. Berardo, D., Cumming, A., & Marleau, G.-D. 2017, ApJ, 834, 149 [NASA ADS] [CrossRef] [Google Scholar]
  16. Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Biller, B. A., Liu, M. C., Wahhaj, Z., et al. 2010, ApJ, 720, L82 [NASA ADS] [CrossRef] [Google Scholar]
  18. Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 2 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bonnefoy, M., Lagrange, A. M., Boccaletti, A., et al. 2011, A&A, 528, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Bonnefoy, M., Currie, T., Marleau, G. D., et al. 2014, A&A, 562, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Bowler, B. P., Liu, M. C., Kraus, A. L., Mann, A. W., & Ireland, M. J. 2011, ApJ, 743, 148 [NASA ADS] [CrossRef] [Google Scholar]
  23. Brittain, S. D., Najita, J. R., Dong, R., & Zhu, Z. 2020, ApJ, 895, 48 [CrossRef] [Google Scholar]
  24. Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856 [NASA ADS] [CrossRef] [Google Scholar]
  26. Cantalloube, F., Por, E. H., Dohlen, K., et al. 2018, A&A, 620, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404 [NASA ADS] [CrossRef] [Google Scholar]
  28. Carnall, A. C. 2017, ArXiv e-prints [arXiv:1705.05165] [Google Scholar]
  29. Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464 [NASA ADS] [CrossRef] [Google Scholar]
  30. Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, A&A, 605, L9 [CrossRef] [EDP Sciences] [Google Scholar]
  31. Cheetham, A. C., Samland, M., Brems, S. S., et al. 2019, A&A, 622, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJ, 877, L33 [Google Scholar]
  33. Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662 [NASA ADS] [CrossRef] [Google Scholar]
  34. Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, SPIE Conf. Ser., 7014, 70143E [Google Scholar]
  35. Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A&A, 622, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Currie, T., Fukagawa, M., Thalmann, C., Matsumura, S., & Plavchan, P. 2012, ApJ, 755, L34 [NASA ADS] [CrossRef] [Google Scholar]
  37. Currie, T., Burrows, A., Madhusudhan, N., et al. 2013, ApJ, 776, 15 [NASA ADS] [CrossRef] [Google Scholar]
  38. Currie, T., Burrows, A., Girard, J. H., et al. 2014, ApJ, 795, 133 [Google Scholar]
  39. Currie, T., Marois, C., Cieza, L., et al. 2019, ApJ, 877, L3 [Google Scholar]
  40. Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., & Green, G. M. 2015, ApJ, 812, 128 [NASA ADS] [CrossRef] [Google Scholar]
  41. Daemgen, S., Todorov, K., Silva, J., et al. 2017, A&A, 601, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. D’Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 [NASA ADS] [CrossRef] [Google Scholar]
  43. Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Dohlen, K., Langlois, M., Saisse, M., et al. 2008, SPIE Conf. Ser., 7014, 70143L [Google Scholar]
  45. Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5 [NASA ADS] [CrossRef] [Google Scholar]
  46. Dupuy, T. J., & Kraus, A. L. 2013, Science, 341, 1492 [NASA ADS] [CrossRef] [Google Scholar]
  47. Dupuy, T. J., & Liu, M. C. 2012, ApJS, 201, 19 [NASA ADS] [CrossRef] [Google Scholar]
  48. Eisner, J. A. 2015, ApJ, 803, L4 [NASA ADS] [CrossRef] [Google Scholar]
  49. Emsenhuber, A., Mordasini, C., Burn, R., et al. 2020a, A&A, submitted [arXiv:2007.05561] [Google Scholar]
  50. Emsenhuber, A., Mordasini, C., Burn, R., et al. 2020b, A&A, submitted [arXiv:2007.05562] [Google Scholar]
  51. Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449 [NASA ADS] [CrossRef] [Google Scholar]
  52. Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109 [CrossRef] [Google Scholar]
  53. Fortney, J. J., Marley, M. S., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Astron. Nachr., 326, 925 [NASA ADS] [CrossRef] [Google Scholar]
  54. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011, ApJ, 739, L41 [Google Scholar]
  56. Ginzburg, S., & Chiang, E. 2019, MNRAS, 490, 4334 [CrossRef] [Google Scholar]
  57. Greco, J. P., & Brandt, T. D. 2016, ApJ, 833, 134 [NASA ADS] [CrossRef] [Google Scholar]
  58. Greenbaum, A. Z., Pueyo, L., Ruffio, J.-B., et al. 2018, AJ, 155, 226 [Google Scholar]
  59. Gregorio-Hetem, J., & Hetem, A. 2002, MNRAS, 336, 197 [NASA ADS] [CrossRef] [Google Scholar]
  60. Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
  61. Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJ, 758, L19 [NASA ADS] [CrossRef] [Google Scholar]
  62. Hashimoto, J., Aoyama, Y., Konishi, M., et al. 2020, AJ, 159, 222 [NASA ADS] [CrossRef] [Google Scholar]
  63. Helling, C., Dehn, M., Woitke, P., & Hauschildt, P. H. 2008, ApJ, 675, L105 [NASA ADS] [CrossRef] [Google Scholar]
  64. Hunziker, S., Quanz, S. P., Amara, A., & Meyer, M. R. 2018, A&A, 611, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Ireland, M. J., Kraus, A., Martinache, F., Law, N., & Hillenbrand, L. A. 2011, ApJ, 726, 113 [Google Scholar]
  66. Isella, A., & Natta, A. 2005, A&A, 438, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25 [Google Scholar]
  68. Jaeger, C., Molster, F. J., Dorschner, J., et al. 1998, A&A, 339, 904 [Google Scholar]
  69. Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43 [NASA ADS] [Google Scholar]
  70. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Lenzen, R., Hartung, M., Brandner, W., et al. 2003, Proc. SPIE, 4841, 944 [NASA ADS] [CrossRef] [Google Scholar]
  72. Linder, E. F., Mordasini, C., Mollière, P., et al. 2019, A&A, 623, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ, 833, 96 [NASA ADS] [CrossRef] [Google Scholar]
  74. Lubow, S. H., & Martin, R. G. 2012, ApJ, 749, L37 [CrossRef] [Google Scholar]
  75. Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 [NASA ADS] [CrossRef] [Google Scholar]
  76. Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, SPIE Conf. Ser., 9908, 990834 [Google Scholar]
  77. Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Marleau, G.-D., & Cumming, A. 2014, MNRAS, 437, 1378 [NASA ADS] [CrossRef] [Google Scholar]
  79. Marleau, G.-D., Klahr, H., Kuiper, R., & Mordasini, C. 2017, ApJ, 836, 221 [Google Scholar]
  80. Marleau, G.-D., Mordasini, C., & Kuiper, R. 2019a, ApJ, 881, 144 [Google Scholar]
  81. Marleau, G.-D., Coleman, G. A. L., Leleu, A., & Mordasini, C. 2019b, A&A, 624, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541 [NASA ADS] [CrossRef] [Google Scholar]
  83. Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Nature, 468, 1080 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  84. Mesa, D., Keppler, M., Cantalloube, F., et al. 2019, A&A, 632, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Milli, J., Hibon, P., Christiaens, V., et al. 2017, A&A, 597, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131 [CrossRef] [EDP Sciences] [Google Scholar]
  87. Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 [NASA ADS] [CrossRef] [Google Scholar]
  88. Mordasini, C. 2013, A&A, 558, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012a, A&A, 547, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Mordasini, C., Alibert, Y., Georgy, C., et al. 2012b, A&A, 547, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Mordasini, C., Marleau, G. D., & Mollière, P. 2017, A&A, 608, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Papaloizou, J. C. B., & Nelson, R. P. 2005, A&A, 433, 247 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Pecaut, M. J.,& Mamajek, E. E. 2016, MNRAS, 461, 794 [NASA ADS] [CrossRef] [Google Scholar]
  95. Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
  97. Preibisch, T., & Mamajek, E. 2008, The Nearest OB Association: Scorpius-Centaurus (Sco OB2) (San Francisco: Astronomical Society of the Pacific Monograph Publications), 5, 235 [Google Scholar]
  98. Rajan, A., Rameau, J., De Rosa, R. J., et al. 2017, AJ, 154, 10 [NASA ADS] [CrossRef] [Google Scholar]
  99. Rousset, G., Lacombe, F., Puget, P., et al. 2003, Proc. SPIE, 4839, 140 [NASA ADS] [CrossRef] [Google Scholar]
  100. Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Sanchis, E., Picogna, G., Ercolano, B., Testi, L., & Rosotti, G. 2020, MNRAS, 492, 3440 [CrossRef] [Google Scholar]
  102. Scott, A., & Duley, W. W. 1996, ApJS, 105, 401 [NASA ADS] [CrossRef] [Google Scholar]
  103. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Snellen, I. A. G., & Brown, A. G. A. 2018, Nat. Astron., 2, 883 [NASA ADS] [CrossRef] [Google Scholar]
  105. Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [NASA ADS] [CrossRef] [Google Scholar]
  106. Stolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, A&A, 621, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182 [EDP Sciences] [Google Scholar]
  108. Sumlin, B. J., Heinson, W. R., & Chakrabarty, R. K. 2018, J. Quant. Spectr. Rad. Transf., 205, 127 [CrossRef] [Google Scholar]
  109. Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853 [NASA ADS] [CrossRef] [Google Scholar]
  110. Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47 [NASA ADS] [CrossRef] [Google Scholar]
  111. Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378 [NASA ADS] [CrossRef] [Google Scholar]
  112. Thanathibodee, T., Calvet, N., Bae, J., Muzerolle, J., & Hernández, R. F. 2019, ApJ, 885, 94 [NASA ADS] [CrossRef] [Google Scholar]
  113. Thanathibodee, T., Molina, B., Calvet, N., et al. 2020, ApJ, 892, 81 [CrossRef] [Google Scholar]
  114. Trotta, R. 2008, Contemp. Phys., 49, 71 [Google Scholar]
  115. Uyama, T., Tanigawa, T., Hashimoto, J., et al. 2017, AJ, 154, 90 [NASA ADS] [CrossRef] [Google Scholar]
  116. Vigan, A., Moutou, C., Langlois, M., et al. 2010, MNRAS, 407, 71 [NASA ADS] [CrossRef] [Google Scholar]
  117. Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJ, 863, L8 [Google Scholar]
  118. Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192 [Google Scholar]
  119. Wang, J. J., Ginzburg, S., Ren, B., et al. 2020, AJ, 159, 263 [CrossRef] [Google Scholar]
  120. Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Zhu, Z. 2015, ApJ, 799, 16 [NASA ADS] [CrossRef] [Google Scholar]

3

For a narrow mass around Mp ≈ 1–1.5 MJ, there are two solutions: a small- and a large-Rp solution, with, respectively, a small (large) Lint and large (small) Lacc, summing up to Ltot = LSED.

4

The data can be visualized at and downloaded from the Data Analysis Centre for Exoplanets (DACE) platform at https://dace.unige.ch

5

At these low masses, the structure models become sensitive to other parameters such as metallicity or core mass, so that this value is to be taken with a grain of metal. However, that no high-mass models are possible here is robust.

All Tables

Table 1

Photometry and error budget.

Table A.1

Astrometry and error budget.

All Figures

thumbnail Fig. 1

Detection of the PDS 70 planetary system and CSD with the NACO NB4.05 (left panel) and M′ (right panel) filters. The images show the mean-combined residuals of the PSF subtraction on a color scale that has been normalized to the peak of the stellar PSF. The flux in the NB4.05 image has been increased by a factor of 1.8 for clarity. The detected emission from PDS 70 b and c (only marginal in NB4.05) is encircled. North is up and east is left.

In the text
thumbnail Fig. 2

Color–magnitude diagram of MM $_{M&#x0027;}$ versus KM′. The field objects are color-coded by M, L, and T spectral types (see discrete colorbar), the young and low-gravity objects are indicated with a gray square, and the directly imaged companions are labeled individually. PDS 70 b is highlighted with a red star. The blue and orange lines show the synthetic colors computed from the AMES-Cond and AMES-Dusty evolutionary tracks at an age of 5 Myr. Blackbody radiation curves are shown for 8 RJ, 4 RJ, and 2 RJ (black dashed lines, from top to bottom). The black arrows indicate the reddening by MgSiO3 grains with amean radius of 0.1 and 1 μm, and AM $A_{M&#x0027;}$ of 0.05 and 2 mag, respectively.

In the text
thumbnail Fig. 3

Color–color diagram of HM′ versus L′–M′. The field objects are color-coded by M, L, and T spectral types (see discrete colorbar), the young and low-gravity dwarf objects are indicated with a gray square, and the directly imaged companions are labeled individually. PDS 70 b is highlighted with a red star. The blue and orange lines show the synthetic colors computed from the AMES-Cond and AMES-Dusty evolutionary tracks at an age of 5 Myr. The black dashed line shows the synthetic colors of a blackbody spectrum. The black arrows indicate the reddening by MgSiO3 grains with amean radius of 0.1 and 1 μm, and AM $A_{M&#x0027;}$ of 0.05 and 5 mag, respectively.

In the text
thumbnail Fig. 4

Spectral energy distribution of PDS 70 b. Top and bottom panels: results from fitting one and two blackbody components, respectively (the flux units are different between the two panels). The photometric and spectroscopic data (colored markers) are shown in comparison with the best-fit blackbody spectrum (black line), and randomly drawn samples from the posterior distribution (gray lines). The residuals are shown relative to the data uncertainties. The Hα and Hβ (upper limit) fluxes (Hashimoto et al. 2020) are shown for reference but were not used in the fit.

In the text
thumbnail Fig. 5

Total luminosity from standard models of isolated planets as a function of mass and radius, but withalso the accretion luminosity from the planet surface shock (Eq. (4) with LCPD = 0). The thick black contour and gray shade highlight the measured luminosity and 5σ uncertainty (see main text for details), log(LSEDL) = −3.79 ± 0.02, and the thin black contours are steps of 1.0 dex relative to LSED. The accretion rate was set to = 0 in the left panel and = 5 × 10−7 MJ yr−1 (Hashimoto et al. 2020) in the right panel. The inferred photospheric radius and 1σ uncertainty, Rphot ± σ, are shown with a horizontally dashed line and gray shaded area. The light yellow curves show as reference the mass–radius predictions by the AMES (Cond/Dusty) structure model at 5 Myr (solid) and 1 Myr (dashed). The dotted white lines in the right panel show Ltot = LSED with scaled by a factor of 0.3 and 3. The crosshatched regions indicate parts of the parameter space for which no predictions are available from the structure model: because of electron degeneracy pressure at small radii, and because no stable hydrostatic structure exists at large radii.

In the text
thumbnail Fig. B.1

Posterior distributions from fitting a single Planck function to the SED of PDS 70 b. The 1D marginalized distributions are shown in the diagonal panels and the 2D parameter projections in the off-axis panels. The bolometric luminosity, log LL, has been calculated from the posterior samples of Teff and R.

In the text
thumbnail Fig. B.2

Posterior distributions from fitting a combination of two Planck functions to the SED of PDS 70 b. Further details are provided in the caption of Fig. B.1.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.