© ESO 2021

1 Introduction

In recent years high-contrast imaging observations have detected a handful of young giant planets on wide orbits (β Pic b, 51 Eri b, HIP 65426 b, HR 8799 bcde, HD 95086 b; Lagrange et al. 2009; Macintosh et al. 2015; Chauvin et al. 2017; Marois et al. 2008; Rameau et al. 2013). Follow-up observations enabled their characterization and continuously provide new insights into their present atmospheric structures and compositions (e.g., De Rosa et al. 2016; Samland et al. 2017; Cheetham et al. 2019; Gravity Collaboration 2019; Stolker et al. 2020b). At the same time, a great deal of effort has been invested in the search for young forming planets still embedded in their protoplanetary disks, with several proposed candidates (e.g., LkCa15 bcd, HD 169142 b, HD 100546 bc, MWC 758 b, Kraus & Ireland 2012; Quanz et al. 2013, 2015; Reggiani et al. 2014, 2018; Biller et al. 2014; Brittain et al. 2014; Sallum et al. 2015; Currie et al. 2017). However, most recent follow-up studies questioned their existence and were not able to confirm the detections (e.g., Follette et al. 2017; Rameau et al. 2017; Ligi et al. 2018; Cugno et al. 2019a).

With two confirmed companions, the PDS70 system (Keppler et al. 2018; Müller et al. 2018; Haffert et al. 2019) offers a unique opportunity to study and characterize protoplanets interacting with their natal environments, sculpting disk gaps (Keppler et al. 2019), and accreting material from their surroundings (Haffert et al. 2019; Wagner et al. 2018; Zhou et al. 2021). Follow-up observations resulted in the first attempts to constrain atmospheric parameters and to describe the planets’ surroundings. Wang et al. (2020b) fitted atmospheric models to determine the radius and mass of the companion PDS70 b (R_p ~ 2−3R_J, M_p ~ 2−4M_J). Stolker et al. (2020a) obtained comparable results fitting a blackbody to the protoplanet spectral energy distribution (SED), indicating that we are actually seeing radiation reprocessed by the planet environment, which is obscuring our view of strong molecular features. More recently, Wang et al. (2021) detected molecular absorption features at the blue and red ends of a K-band spectrum obtained with VLTI/GRAVITY and found that atmospheric models better fit the planet SED. These latest findings support the tentative detection of H₂O feature at 1.4 μm seen in SPHERE/IFS data (Müller et al. 2018). Moreover, they constrained the mass of PDS70 b to be < 10 M_J in order for the system to be dynamically stable.

Facchini et al. (2021) provided the chemical inventory of the PDS70 disk, revealing molecular emission from 12 molecular species with radial variation indicating dynamic chemistry. The comparison of their data with thermo-chemical models of protoplanetary disks strongly suggests a disk-average C/O ratio >1. In addition, they detected emission from CO in the gap hosting the two planets, indicating that CO is likely currently accreting onto the protoplanets’ atmospheres.

It has been suggested that determining the chemical composition of a planet’s atmosphere might play a key role in understanding where, with respect to snowlines, the planet formed (Öberg et al. 2011) and to reconstruct the migration history (Madhusudhan et al. 2014; Mordasini et al. 2016). Planet formation history strongly impacts the chemical enrichment of gas giant atmosphere, which depends on the chemical properties of the disk gas from which the planet envelope is accreted. This, in turn, is affected by disk evolution processes like grain growth, radial drift, and ionization (Helling et al. 2014; Öberg & Bergin 2016; Eistrup et al. 2018; Cridland et al. 2019). Hence, deriving elemental abundance ratios of PDS70 b, still embedded in its natal environment, offers a unique possibility to verify current chemical modeling of atmosphere enrichment at the very early stages.

Spectral cross-correlation (e.g., Sparks & Ford 2002) is a widely used technique in exoplanetary science. It was used by Snellen et al. (2014) to measure the rotational velocity of β Pic b with the high-resolution spectrograph VLT/CRIRES, and by Ruffio et al. (2019) to estimate the orbital radial velocity of the HR8799 bc planets. It is regularly used in the analysis of transit spectroscopy datasets (e.g., Hoeijmakers et al. 2015; Giacobbe et al. 2021). A similar approach has been used on directly imaged exoplanets to infer the presence of molecules in the atmosphere. Using Keck/OSIRIS data, Konopacky et al. (2013) and Barman et al. (2015) identified spectral features from CO, CH₄, and H₂O in the spectrum of HR 8799 b and c, and used cross-correlation of the planet spectra with molecular absorption templates to successfully confirm the presence of H₂O and CO in the planets atmospheres.

Since the spectral signal emitted by the planet is fundamentally different than that of the star, mainly because of the presence of molecules in the atmosphere, medium-resolution spectroscopic data can be used to disentangle the planet spectrum from the stellar speckles in the image. Hoeijmakers et al. (2018) adapted the cross-correlation analysis technique to be used not only on measured planet spectra, but also to recognize weak molecular signatures in residuals of medium-resolution spectroscopy data, obtained after removal of the stellar spectrum, instrument systematics, and planet pseudo-continuum. They successfully applied molecular mapping to K-band β Pic data from VLT/SINFONI, detecting H₂O and CO at the planet location with high confidence (Hoeijmakers et al. 2018), demonstrating the efficacy of this approach. Later, Petit dit de la Roche et al. (2018) applied molecular mapping on the HR 8799 b planet and Petrus et al. (2021) detected H₂O and CO in the atmosphere of HIP 65426 b. Petrus et al. (2021) also compared the atmospheric parameters obtained with the molecular mapping approach with those obtained via matching of synthetic spectra to the data using Bayesian inference. They found that most of the parameter estimates from the two methods agree, confirming that cross-correlation is a powerful tool to investigate giant planet atmospheres.

In this paper we analyze for the first time medium-resolution spectroscopy data of the PDS70 system, with the goal of searching for signatures of molecular absorption, estimating molecular abundances and elemental abundance ratios to better describe the formation history of the planetary-mass object PDS70 b. In Sect. 2 we describe the data we use in this work. In Sect. 3 we give details of our data reduction and stellar PSF subtraction. In Sect. 4 we report on our search for molecular signals and the detection limits, which we discuss in Sect. 5. Finally, we summarize our conclusions in Sect. 6.

2 Observations

We observed the PDS70 system with the SINFONI integral field spectrograph (IFS, Eisenhauer et al. 2003; Bonnet et al. 2004) mounted at the European Southern Observatory’s Very Large Telescope (VLT) facility in Chile. The observations were carried out in two separate blocks on 2019 May 24 and 2019 June 03. Each observing block consisted of 32 exposures, each with DIT[s] × NDIT = 20 × 8, where DIT is the detector integration time and NDIT the number of DIT repetitions per exposure. Thus, the total time on target each night was ~85 min. In addition to the standard calibration, we took sky frames every two science exposures to sample the background. The data were taken in the K band in the highest spatial resolution mode with a plate scale of 12.5 × 25 mas and a field of view of 0.′′8 × 0.′′8. The spectralresolution was R ≈ 5075, measured from the width of OH line spread function. The images were taken in pupil stabilized mode to enable angular differential imaging (ADI, Marois et al. 2008). The achieved field rotation was 62.4° and 78.0° for the first and the second night, respectively. We dithered the position of the star across the detector to reduce the impact of bad pixels during post-processing. For both nights the standard star HIP 81214 was observed for flux and telluric calibration. The median Differential Image Motion Monitor (DIMM) seeing was similar during the two nights, 0.′′65 ± 0.′′07 and 0.′′70 ± 0.′′10 respectively, while the estimated coherence time τ₀ indicates that the second night had on average less turbulent conditions, with τ₀ = 3.1 ± 0.5 ms and τ₀ = 6.5 ± 1.1 ms, respectively. Despite the mean value, we note that both seeing and coherence time indicate that the second night suffered from more variable conditions. Finally, the second observing night had ~ 5 times lower relative humidity.

3 Data reduction and quality assessment

3.1 Basic reduction

The two datasets were calibrated using the SINFONI EsoReflex pipeline (version 3.1.1; Abuter et al. 2006), which included dark subtraction, flat fielding, bad pixel removal, wavelength calibration using neon and argon lamps, detector linearity correction, and extraction of spectral cubes from the raw frames. The EsoReflex pipeline outputs cubes consistent of 2216 frames, covering a spectral range from 1.929 to 2.472 μm.

Since images at the edges of the waveband consisted mostly of NaNs and the shorter wavelength data were dominated by telluric features, we first restricted the spectral range to include frames between 2.088 and 2.440 μm and removed one row or one column at each image side. At this point, only a few NaNs remained in the frames at the location of known bad pixels, which were removed by replacing them with the median of the eight surrounding spatial pixels. Across the different wavelengths, the stellar photocenter shifted by up to one pixel within each cube, especially atlonger wavelengths. Therefore, a 2D Gaussian function was fitted to each individual frame and to the wavelength-combined (median of the cube along the wavelength axis) frames. Then images were shifted to a common photocenter given by the wavelength-collapsed images so that the position of the star is the same throughout the whole cube.

Fig. 1 PDS70 system; the b planet is shown within the solid circle. Six principal components were removed from the first dataset to detect the planet with S∕N  ~ 12. The dashed circle shows the expected position of PDS70 c (Mesa et al. 2019).

3.2 Stellar PSF subtraction

3.2.1 PCA-ADI approach

After the basic reduction described above, each cube was median-combined along the wavelength axis and we obtained two separate datasets for the observing nights, each of which contained 32 images. The images were cropped to a size of 0.′′55.

The stellar PSF was removed using PynPoint, which uses an algorithm based on Principal Component Analysis (PCA) to reconstruct and subtract the principal modes of the PSF (Amara & Quanz 2012; Stolker et al. 2019). We subtracted 1–15 components and used a central mask of 0.′′11 in radius to remove most of the central part of the image without affecting the flux of PDS70 b. Afterwards, the images were derotated to the same orientation according to the parallactic angle and median-combined. We show the final image obtained for the first dataset in Fig. 1, while planet b is hardly detectable in the second dataset. We investigated this difference in Appendix A, concluding that because of more stable weather conditions during the observations the first dataset provides the best chance to study the companion. Thus, in the following we only use the dataset taken on 2019 May 24. Nevertheless, we verified that the results presented in Sect. 4.2 apply to the second dataset and to their combination as well. In Appendix B we treat the wavelength-collapsed cube as an ADI sequence, and we investigate the contrast and position of PDS70 b. Finally, we note that according to the detection limits it is not possible to detect PDS70 c in our data. In the images, a point-like emission is somewhat visible at the planet location when selecting the best combination of number of principal components (PCs), image and central mask sizes. However, the signal-to-noise ratio (S/N) remains low (≲ 5) and it is not clear whether the emission comes from PDS70 c or scattered light from the inner disk rim. Therefore, the following analysis focuses on PDS70 b.

3.2.2 High-resolution spectral differential imaging (HRSDI)

Compared to Sect. 3.2.1, where we focused mainly on our search for the continuum emission, here we explain how we prepared the data cubes for molecular mapping, meaning that we wanted to remove the contribution of stellar spectrum from each spatial pixel, while preserving other spectral features like those produced by molecules in a planet atmosphere. Here, the planet continuum is not needed and it can be removed together with the stellar PSF. As in Hoeijmakers et al. (2018), we used the SkyCalc tool from ESO (Noll et al. 2012; Jones et al. 2013) to verify whether the wavelength solution found by the SINFONI pipeline matches the telluric lines of the synthetic spectrum. We found a mismatch of ~ 60 ± 9.8 km s⁻¹ (~ 2 pix) between the wavelength solution delivered by the pipeline and the sky simulations provided by SkyCalc. We corrected this effect by shifting the wavelength of each frame to the rest frame given by the SkyCalc spectrum.

After the basic reduction explained in Sect. 3.1, the PSF is removed using a spectral differential technique based on Hoeijmakers et al. (2018) and Haffert et al. (2019). Each spatial pixel is divided by its total flux after 12σ clipping to remove remaining hot pixels, and all the resulting spectra are median-combined to generate a reference spectrum, which does not include any planet signal because of the use of the median. The spectrum of each spatial pixel is then divided by this reference spectrum. The residuals are low-pass filtered using a Savitky–Golay filter (Savitzky & Golay 1964) of order 3 with a window of 31 channels (~ 7.6 nm) to remove spectral features that are as broad or broader than the filter size, leaving only high-frequency features. The filtered spectrum is then multiplied again by the reference spectrum. Modeling λ-dependencies over the PSF shape, this process delivers a continuum spectrum of the stellar PSF at each spatial pixel that does not include spectral features coming from the planet and the star. The final residuals are given by the subtraction of the initial spectrum from the local reference filtered spectrum. The spatial pixels of each individual cube are then used for building a reference library through PCA to remove correlated noise within the cube itself. This step is performed cube by cube and each time the first 20 components were removed.

At this point, correlated (instrumental) noise between cubes still dominates the central part of the images and high-frequency residuals that remained in the data could make the detection of molecular signatures from planets still difficult. In particular, telluric residuals could mimic the detection of molecular lines (especially H₂O) at any location in the image, and H₂O and CO emission lines from the inner unresolved circumstellar disk could be detected with cross-correlation (see Sect. 4.2). Those features could have counterbalanced absorption features from the planet atmospheres of the same molecules, impeding their identification. We found that applying PCA a second time is the best method for removing these last noise sources. This time we used the entire dataset (32 cubes) to build the reference library, and we removed two principal components.

After removing the stellar spectrum from each pixel, the star was centered, the cubes were derotated to a common field orientation using their parallactic angles, and the frames corresponding to the same wavelength were median combined, resulting in a residuals cube made of 1479 frames. Finally, a mask of radius 0.′′11 was applied at the center of each frame.

4 Analysis and results

In this section we describe how we exploited the R ≈ 5075 resolution of the SINFONI spectrograph in order to detect and, if possible, quantify the presence of molecules in the atmosphere of PDS70 b. To this end, we cross-correlated the residuals of the data with templates containing characteristic absorption features from specific molecules detected in the atmosphere of low-mass objects.

4.1 Spectral templates

In order to use cross-correlation and identify signatures from individual molecules, we generated templates with petitRADTRANS (Mollière et al. 2019), a Python package for calculating transmission and emission spectra of exoplanets. We used the thermal structure grids calculated with petitCODE (Mollière et al. 2015, 2017) provided in Samland et al. (2017) for objects with an effective temperature of ≈ 1200 K (Stolker et al. 2020a), solar metallicity (as assumed for PDS70A by Müller et al. 2018), and log (g) = 3.0. The surface gravity was chosen to be 3.0 as it was the closest available value assuming M_pl ~ 1.5 M_J and R_pl = 3.0 from Stolker et al. (2020a). Using the vertical thermal profile we created high-resolution emission spectra of planets, with an atmosphere where only one molecular species (H₂O, CO, or CH₄) is present. The mass fraction of the complex molecule was arbitrarily chosen based on two criteria: it must create well-defined molecular absorption lines without saturation, and it must be compatible with mass fractions previously estimated in low-mass objects (e.g., Lavie et al. 2017). These were chosen to be log (XH 2O) = −2.7, log (X_CO) = −2.5, and log (XCH 4) = −4.0, where X represents the mass fraction of the molecule. The mass fractions of H₂ and He wereset to 75 and 25% of the remaining atmosphere mass (Line et al. 2012). We focused on these molecules as chemical equilibrium models predict them to be the most abundant molecules with significant features in the observed wavelength range at this temperature (see Fig. 11 of Stolker et al. 2020b). The thermal structure was not adapted to the molecular opacities as detections mostly depend on the position of the absorption lines, which, unlike line depths, does not change when altering the p–T profile of the atmosphere. Relative abundances were kept constant with altitude (e.g., Todorov et al. 2016; Line et al. 2017). In the calculation H₂ –H₂ and H₂–He continuum opacities were included, while clouds were not considered. To the first order, the presence of clouds would only lower the amplitude of the spectral features, therefore scaling the cross-correlation signal without impacting the S/N of potential detections.

The resulting simulated spectra were resampled to the spectral resolution of the data, and a Savitzky–Golay filter (Savitzky & Golay 1964) with a window of 31 channels (Δλ ≈ 7.6 nm) and order 3 was applied to approximate the continuum emission. The latter was subtracted from the original spectrum in order to obtain the spectral signature of the molecules without low-frequencies contributions from the continuum. The molecular templates, together with the residuals that were cross-correlated with the data, are shown in Fig. 2.

4.2 Cross-correlation

To search for signals due to molecular absorption from H₂O, CO, and CH₄, each spatial pixel in the residual cube is cross-correlated with each of the molecular templates over a range of radial velocities between − 2500 and 2500 km s⁻¹ in steps of 10 km s⁻¹. For each of the templates this process resulted in a new cube made of 501 images, where each image represents the cross-correlation function (CCF) at a different radial velocity (Hoeijmakers et al. 2018). The more similar the model template is to the real emission spectrum of the planet, the larger the value of the CCF at the location and radial velocity of the planet will be.

To verify that our methods are effective in detecting molecules, we downloaded archival data of β Pic taken in the same mode as our PDS70 data and reduced them looking for H₂O and CO signals, as presented in Hoeijmakers et al. (2018). We successfully identified the presence of both molecules (see Appendix C).

Once the pipeline was verified, we searched for peaks in the CCF at the location of PDS70 b found in Appendix B at radial velocities of RV_pl = RV_orbit + RV_bary = 4.3−9.9 ≈−5.6 km s⁻¹, where the first term is the radial velocity component of the Keplerian velocity with respect to the system velocity (Haffert et al. 2019) and the second term represents the barycentric velocity during the observation calculated with the PyAstronomy Python package¹. Given the instrument resolution, each channel corresponds to a shift in radial velocity of ~ 30 km s⁻¹, indicating that a higher precision in the determination of the planet radial velocity is not necessary. No signal was detected for the three molecules over the full RV range. The panels in the top row of Fig. 3 show the molecular maps at RV_pl for H₂O, CO, and CH₄, while the CCFs at the planet position are shown in orange in the bottom row of the same figure for the three molecules.

In addition to the pixel-by-pixel approach explained above, we summed the residual flux at each channel in a circular aperture of 0.′′0275 in radius placed at the planet location to boost the S/N of molecular features. No significant peak in the CCF was detected.

4.3 Detection limits

To assess the quality of the data and the scientific relevance of the non-detection of molecules in the photosphere of PDS70 b, we estimated the limits on molecular abundances by inserting artificial companions in the data before subtracting the PSF and looking for the lowest molecular abundance that could be detected with our method. To do that, we fixed the thermal structure of the atmosphere to self-consistently calculated p–T profiles as we lack empirical constraints on the protoplanet’s thermal structure. Under this assumption, determining precise abundance limits from emission spectra is a difficult task, as changing abundances of the main absorbers theoretically involves major changes in the p–T structure and the emitted SED that could not be taken into account.

Planet spectra were generated with petitRADTRANS, with the thermal structure for T_eff = 1200 K, log (g) = 3.0 from Samland et al. (2017). We varied X_mol, the mass fraction of the molecule we are looking for, while keeping the amount of H₂ and He constant in the rest of the atmosphere (75% H₂, 25% He). The planet was placed at a distance of 113.4 pc (Gaia Collaboration 2018) and the planet radius was assumed to be 3.0 R_J (Wang et al. 2020b; Stolker et al. 2020a). Once generated, the planet emission spectrum was Doppler-shifted by − 5 km s⁻¹ and down-sampled to the spectral resolution of the data. Then the contrast between the injected signal and the stellar flux (estimated in Appendix D) was calculated. Finally, in each frame a copy of the stellar PSF rescaled by the contrast wasinserted in the pre-PSF-subtracted cubes at the position of PDS70 b from Appendix B. The new data with the inserted planets were reduced with the same pipeline detailed in Sect. 3.2.2, and the same templates from Sect. 4.2 were cross-correlated with the residuals to find molecular absorption features.

We estimated the significance of the detection calculating its S/N as the ratio of the peak value of theCCF to the standard deviation of its remaining values at least 20 RV steps (i.e., 200 km s⁻¹) apart from the peak. To correct for autocorrelation effects on the noise measurement, particularly strong in the case of CO, we calculated the autocorrelation function of the templates, rescaled it so that its peak is at the same level as the data-template CCF peak, and then subtracted it from the original CCF for |RV| > 200 km s⁻¹. This way the autocorrelation signal is not counted in estimating the noise. Finally, we adjusted X_mol to obtain S∕N ≈ 5.

Within our assumptions on the atmospheric parameters, the VLT/SINFONI data exclude log (XH 2O) > −4.0, log (X_CO) > −4.1 and log (XCH 4) > −4.9 because molecular mass fractions larger than these would have been detected by our methods with S∕N > 5. The middle row of Fig. 3 shows the molecular maps of PDS70 in H₂O, CO, and CH₄ at the RV of the artificial planet for the limit mass fractions provided above. The signal of the three molecules is detected with S∕N ≳ 5 as shown in the CCF of the bottom row panels (blue lines). The non-detections of H₂O and the low abundance limits we calculated within the framework of our assumptions appear to be in disagreement with the GRAVITY spectrum (Wang et al. 2021), where the presence of H₂O in the atmosphere was inferred thanks to the spectral slope at the red end of the spectrum. We address this possible inconsistency in Sect. 5.5.

In order to account for systematic wavelength calibration effects from the EsoReflex pipeline we tested the impact of a range of errors on the cross-correlation signal. Using the known spectral position of telluric lines in the data we estimated the mean error (σ = 8.4 km s⁻¹) of the wavelength calibration across the field of view. We proceeded to distort the wavelength coordinates of a model spectrum by adding a sine function as a systematic error. The same non-distorted spectrum was used as the template and was cross-correlated with the distorted spectrum. We varied both the frequency and amplitude of the sine function, constrained by the size of the spectrum array (the largest frequency being a single wavelength bin and the smallest the size of the array) and the measured wavelength calibration error, respectively. The H₂O and CO templates were tested individually since wavelength calibration errors might have different impacts on each molecule. We found a maximum loss in S/N of 20% in the extreme case. Thus, the molecular non-detection cannot be explained as instrumental calibration effects

Fig. 2 Molecular templates after removal of the continua (H2O in blue, CO in red, and CH4 in green) and residuals of the data at the position of the planet (gray). No broad feature is visible and only individual absorption lines at the instrument resolution are cross-correlated. The small black dash in the top left corner of the first three panelsrepresents the filter size for continuum removal. Features smaller than this are kept in the residuals, while larger structures are removed as part of the continuum.

Fig. 3 Molecular maps of the PDS70 system and CCFs at the position of the b planet. Top row: H2O, CO, and CH4 molecular maps at RV ~ − 5 km s−1. The central star represents the position of PDS70 A, while the yellow circle indicates where a signal from PDS70 b should be found. Middle row: same as in the top row, but a planet with the molecular abundances indicated at the top of each panel has been inserted in the data and detected with S∕N ~ 5. Bottom row: CCF as a function of the radial velocity at the location of the planet. The orange lines represent for each molecule the CCF of the data, while the blue lines are the CCFs of the inserted artificial planets.

5 Discussion

As chemical equilibrium models predict very low abundances of methane for objects with T_eff ≈ 1200 K (e.g., Zahnle & Marley 2014), we focus the upcoming discussion on H₂O and CO only, as they are the most common molecules expected in self-luminous objects according to chemical models (Zahnle & Marley 2014) and existing observations (Hoeijmakers et al. 2018; Petit dit de la Roche et al. 2018; Petrus et al. 2021; Miles et al. 2020).

5.1 Temperature dependence

The uncertain physics of PDS70 b prevent a clear determination of the thermal structure of the planet. As atmospheric chemistry is highly dependent on the temperature (e.g., Zahnle & Marley 2014), we evaluated the impact of varying the planet’s assumed temperature on the upper limits of the molecular abundances. We estimated the S/N of molecular detections for effective temperatures in the range T_{eff, pT} = 700–1500 K (in steps of 50 K), with log(X_mol) between − 4.5 and − 1.0 dex (in steps of 0.1 dex). Here we denote by T_{eff, pT} the effective temperature used to calculate a given p–T profile (keeping log(g) = 3.0 fixed). The generated spectra have a different T_eff because they are calculated with different abundances. However, for the T_{eff, pT} = 1200 K case, the K-band spectrum of the inserted planet agrees with the flux measured by SPHERE/IRDIS for PDS70 b (Müller et al. 2018), while T_{eff, pT} = 1500 K implies a ~50% brighter planet. An overall brighter planet would induce a stronger signal in the CCF and a given spectral signature would be easier to detect. Hence, we conclude that the flux inconsistencies due to the choice of the thermal structure may lead to a minor underestimate of the molecular abundance limits of the order of 0.3 dex at most, decreasing when nearing T_{eff, pT} = 1200 K. This value was calculated by varying the planet radius in order to match the SPHERE/IRDIS measured K-band flux.

Furthermore, realistic emission spectra are likely the result of absorption features produced by multiple molecular species at different altitudes, thus polluting the signal of each individual molecule. Therefore, when generating the planet spectrum with varying X_mol, we also include the amount of H₂O/CO found by Wang et al. (2020a) in HR8799 c (T_eff ≈ 1200 K, similar to PDS70 b).

In the same way described in Sect. 4.3, artificial planets were inserted in the data and we searched for molecular signatures calculating the S/N. The resulting charts of S/N as a function of X_mol are shown in Fig. 4 for H₂O and CO. The solid black lines show where S∕N = 5. Our data are more sensitive for hotter planets, where in the case of a clear atmosphere we expect to be able to detect mass fractions down to ~10^−3.6. The difference between this value and the one presented in Sect. 4.3 and Fig. 3 is due to the polluting effect of the other molecule in the emission spectrum. At lower temperatures, molecular mapping is less effective for two main reasons: cooler planets are fainter and produce weaker spectral signals, and absorption lines are less prominent as the difference in temperature between atmospheric layers is smaller and the thermal structure tends to be more isothermal. Below 850 K, our method is not sensitive to these molecules. This is not necessarily the case for fully self-consistent temperature structures, as changing X_CO and XH 2O results in the p–T profile being adapted.

For comparison, Fig. 4 also reports the estimated abundances of the directly imaged companions with available measurements in the evaluated temperature range. Abundances for HR 8799 b, HR 8799 e and β Pic b planets are taken from Lavie et al. (2017), Mollière et al. (2020) and Gravity Collaboration (2020). For HR 8799 c we used Wang et al. (2020a), who used atmospheric retrievals on Subaru/CHARIS, Gemini/GPI, and Keck/OSIRIS data to constrain atmosphericproperties of the planet. For 51 Eri b, we used the retrieval analysis performed by Whiteford et al. (2020). We note that our method would not have been able to find molecules in the atmosphere of 51 Eri b, the main reason being the much lower effective temperature (Whiteford et al. 2020).

We note that Fig. 4 strongly depends on the temperature profile of the planetary atmosphere. A profile with a smaller temperature gradient between the inner and outer layers reduces the molecular absorption signature, and therefore our ability to detect such molecules. Pushing this consideration to its limit, an isothermal profile causes the atmosphere to be featureless, resulting in a blackbody emission. Under the assumption of a clear atmosphere,Fig. 4 suggests that the chemistry of PDS70 b is different from that of other directly imaged companions with similar effective temperatures. Furthermore, the non-detection of CO in the planet atmosphere conflicts with the detection of CO gas in the disk gap harboring the protoplanet, as the latter is likely directly accreted onto the young planet (Facchini et al. 2021). This inconsistency led us to consider other sources of higher opacity.

Fig. 4 Detection maps of H2O (left) and CO (right) for clear atmospheres as a function of mass fraction Xmol and effective temperature Teff, pT of the planet. The symbols represent measurements for other directly imaged companions (see left panel inset). The thick solid lines separate the region with S∕N > 5 from regions with S∕N < 5, i.e., they represent the molecular detection limit as a function of planet temperature. The dashed, dotted, and thin solid lines represent the molecular detection limits for different cloud configurations based on Mollière et al. (2020), as described in the main text (Pb = 10−3 bar is not shown). The temperature of PDS70 b ranges between 1200 and 1400 K (Stolker et al. 2020a; Wang et al. 2021). The effectivetemperature of β Pic b is well constrained as Teff = 1742 ± 16 K (Gravity Collaboration 2020) and the arrow in its symbol indicates that it is higher than the maximum Teff, pT considered here and should not be interpreted as a lower limit.

5.2 Clouds

Clouds have been found to be key elements to describe the atmosphere of many gas giant planets (e.g., Madhusudhan et al. 2014; Mollière et al. 2020; Gao et al. 2021). Their presence causes the molecular features to be damped, attenuating their signature in the planet spectra and shifting the SED to the red. Therefore, the abundance limits calculated in Sects. 4.3 and 5.1 are not valid when clouds are included. Although their relevance in describing gas giant atmospheres is undisputed, the details of cloud physics remain difficult to describe and parameterize, and uncertainties and assumptions still dominate their treatment. Hence, estimating upper limits on the abundances in the presence of clouds would cause the parameter space to be too large to be explored. Instead, we used a qualitative approach in which we present the impact that four concrete cloud configurations have on the molecular mass fraction upper limits. The cloud parameterization is based on the findings of Mollière et al. (2020) for the planet HR 8799 e, which has T_eff = 1154 ± 50 K, a value similar to that of PDS70 b (Stolker et al. 2020b). Briefly, in Mollière et al. (2020) five cloud parameters were free to explore the parameter space: mass fraction of Fe and MgSiO₃ at the cloud base compared to their equilibrium chemistry value, settling parameter f_sed, eddy diffusion coefficient K_zz, and width of the log-normal size distribution of the particles σ_g. For each of these parameters we used the value reported in Fig. 4 of Mollière et al. (2020), and we used Eq. (3) from the same paper to describe cloud mass fractions as a function of the altitude. The cloud base pressure P_b was manually set to 1, 0.1, or 10⁻³ bar for both cloud species to obtain three of the four investigated scenarios. The altitude of the clouds in the planet atmosphere strongly depends on the temperature structure. Exploring the parameter space for P_b therefore allows us to consider the impact of different clouds independent from the assumed thermal structure. Finally, the fourth cloud scenario is obtained multiplying the Fe and MgSiO₃ mass fractions at P_b = 0.1 bar with a factor of two (i.e., doubling the cloud’s optical thickness).

We repeated the calculation from Sect. 5.1 and we report the lines at S∕N = 5 in Fig. 4. The case with P_b = 10⁻³ bar is not shown as clouds at high altitudes present low particle densities that are unable to affect the molecular spectral features. Lines for this scenario would overlap the solid lines for clear atmospheres. As expected, as a consequence of clouds, the upper limits we derived are higher, meaning that only larger molecular abundances could have been detected. These values are only valid for this specific realization and they would be different for different parameterizations. Nevertheless, they give a sense of the impact that clouds might have on the detectability of molecules in planet atmospheres. Again, molecular mapping is very sensitive to the presence of H₂O and CO for hot objects. Under these assumptions our data strongly suggests that much higher optical depths than those used above are necessary to hide the spectral features produced by molecules, and even clouds may not be sufficient to explain the non-detections in the SINFONI data, assuming that the chemistry of PDS70 b is similar to that of the other directly imaged objects.

Fig. 5 Signal-to-noise ratio of the detection of H2O and CO as a function of extinction affecting the emission spectrum of the planet. Assuming PDS70 b has the same chemical composition as HR 8799 c, an additional extinction with AV ≈ 16−17 mag (AK ≈ 1.2 mag, using the extinction law from Wang & Chen 2019) is required in order not to detect molecules with S∕N > 5. The vertical dashed line shows the maximum value obtained by SED fitting from Wang et al. (2021). This extinction would have revealed molecules in the planet atmosphere with S∕N ~ 7–9.

5.3 Extinction

Assuming a composition for the planet similar to what has been measured for other gas giants, we can estimate the minimum extinction necessary toimpede the detection of the absorption spectrum with molecular mapping. Following Wang et al. (2021), we used the near-infrared (1.0−3.3 μm) extinction law derived by Wang & Chen (2019) to attenuate the emitted spectrum, with A_λ ∝ λ^−2.07. We applied an extinction factor 10^{−A_λ∕2.5} to the generated planet spectrum, where A_λ is the extinction to the radiation emitted by a planetary atmosphere at wavelength λ. Applying such an absorbing element causes the line features to be less prominent in the residuals (although the line contrast relative to the planet continuum remains the same), and therefore the molecules are harder to detect. We benchmarked again our data againstthe molecular abundances from Wang et al. (2020a) for HR 8799 c. We applied different levels of visual extinction A_V to its emission spectrum, from 0 mag to 20 mag in steps of 0.2 mag, using Wang & Chen (2019) to calculate the extinction in the K band. For each of the assumed values we estimated the S/N of the detection of H₂O and CO. These are shown in Fig. 5. Assuming the HR 8799 c chemical composition, our data indicates that an extinction of A_V ≈ 16−17 mag (A_K ≈ 1.2 mag) is required to attenuate molecular signatures ofH₂O in the object’s emission spectrum, while a lower extinction (A_V ~ 10 mag) is necessaryto explain the non-detection of CO. These values, particularly the one for H₂O, are much higher than the highest value from Wang et al. (2021), which is A_V = 8.3 mag when coupled with Exo-REM atmospheric models (Charnay et al. 2018). Assuming only one component for the flux, our analysis suggests that the intrinsic planet emission would need to be a factor ~ 3 larger than is measured in the K band in Sect. 3.2.1 and in previous observations (Müller et al. 2018; Christiaens et al. 2019a; Wang et al. 2021). A larger radius could be invoked to explain the larger intrinsic flux, but this would in turn increase again the amplitude of molecular features in the spectrum, which is inconsistent with non-detection via molecular mapping. To attenuate the absorption signature while maintaining the flux at the observed level, additional radiation reprocessed by the planet dusty environment, likely coming from both the planet interior and the accretion shock, needs to be invoked (Stolker et al. 2020a; Wang et al. 2021). Following Stolker et al. (2020a), we can use Ṁ = 5 × 10⁻⁷M_J yr⁻¹, M_p = 1M_J, and R_p = 2.0 R_J to estimate the accretion luminosity L_acc = G M_p Ṁ∕R_p ≃ 7.0 × 10⁻⁵ L_⊙. If we consider the accretion process to be fully thermalized, we expect its reprocessed emission to appear in form of a blackbody. To the first order we can use the photometric radius R_phot = 3.0 R_J from Stolker et al. (2020a) to constrain its effective temperature (~950 K) and its contribution to the K-band flux (~30% of the flux measured by GRAVITY, Wang et al. 2021). We note that Stolker et al. (2020a) used the lower limit of the mass accretion rate measured by Hashimoto et al. (2020), and that using higher mass accretion rates increases L_acc and in turn the K-band contribution. A more sophisticated radiative transfer modeling will be crucial to investigate this scenario.

At this point, it is not possible to conclude where the dust responsible for such a high extinction is located, meaning that it could be (i) in the form of very thick atmospheric clouds, as we addressed in Sect. 5.2; (ii) external to the planet in the form of a dusty shell enshrouding it; or (iii) due to circumstellar disk material. A combination of these elements is also possible. van Holstein et al. (2021) estimated an upper limit on the degree of polarization in K_s -band of 5% which may exclude the presence of a highly inclined circumplanetary disk (that partially obscures the planet atmosphere) as origin of the extinction. The presence of circumstellar disk material depends on the mass of PDS70 b. More massive objects are expected to carve larger gaps and clear them from circumstellar material, therefore lowering extinction effects due to the protoplanetary disk (see, e.g., Sanchis et al. 2020; Szulágyi & Ercolano 2020). A mass in the 1− 4 M_J range was inferred for PDS70 b (Wang et al. 2020b, 2021; Stolker et al. 2020a), which might not be high enough to completely clear the protoplanetary disk gap. However, the presence of a second companion in the same cavity and at larger separations, together with the dust being trapped in the outer ring (Keppler et al. 2019), could prevent replenishment of new material and keep the gap clear.

Conversely, when adopting the estimated A_V from Wang et al. (2021) (A_V ~ 5−8 mag), we conclude that the molecular abundances in the atmosphere of PDS70 b are lower than what we observed in more mature gas giant planets of similar temperature. Precise measurements on a sufficiently large number of planetary mass objects is required in order to further investigate this possibility.

The absence of absorption features together with the possibility of a very high extinction as indicated by our analysis could suggest a scenario similar to Fomalhaut b (Kalas et al. 2008), where scattered light is detected from a point-like source located at ~18 AU. As Fomalhaut b emission was only detected at optical wavelengths, and not in thermal emission, several ideas have been proposed to explain the observations, some involving a planetary object. Examples are a low-mass planet surrounded by a ring system (Kalas et al. 2008) or by a swarm of collisional satellites (Kennedy & Wyatt 2011), but also a transient dispersing dust cloud (e.g., Janson et al. 2012; Gaspar & Rieke 2020) or a background object (Neuhäuser et al. 2015). However, despite the similarities, several empirical observations of PDS70 b seem to disprove this hypothesis. First, a scattered light scenario cannot explain the dip at the end of the GRAVITY spectrum (Wang et al. 2021). Second, the Hα line emission measured in Haffert et al. (2019) appears to be shifted with respect to the expected radial velocity if it were scattered light from material moving with Keplerian motion, and its width is narrower than that from PDS70A. Third, it is possible that PDS70 b is located in the shadow of the inner disk resolved by Wang et al. (2021), thus it might be only weakly irradiated from the star. Finally, in reflected light the protoplanet’s SED is expected to follow the stellar SED.However, the measured flux at short wavelengths decreases rapidly below the detection limit in the Y band (Müller et al. 2018). In summary, although scattered light could contribute to its emission spectrum, it seems unlikely that this is the dominant source of radiation for PDS70 b.

Fig. 6 Comparison of the detectability of molecules in the SINFONI data with measured abundances from substellar objects. The blue squares represent directly imaged companions HR 8799 bce, β Pic b, and 51 Eri b (Lavie et al. 2017; Wang et al. 2020a; Mollière et al. 2020; Gravity Collaboration 2020; Whiteford et al. 2020); the blue lines show the measured range for H2O in κ And b (Todorov et al. 2016). The blue arrow indicates that the measured value is higher than the range of abundances shown in the plot. The violet circles represent a sample of T dwarfs analyzed in Line et al. (2017), and here the arrows refer to upper limits. The solid black line represents the region where no molecule would have been detected in our SINFONI data because the mass fractions were too low. The dashed black line encloses the region where, conversely, both molecules would have been detected.

Fig. 7 Effect of the extinction on the spectral measurement of the pseudo-continuum and on molecular mapping. Top: extincted best-fit BT-Settl models resampled at the SINFONI (gray) and at the GRAVITY (colored circles) spectral resolution. Even with high extinction, the broadband absorption features at the edges of the K-band spectrum can be identified. Bottom: spectral features at the SINFONI resolution after the removal of the continuum. Higher extinctions decrease the spectral signatures molecular mapping is looking for, impeding the detection of molecules.

5.4 Chemical comparison with other objects

We compare PDS70 b chemical abundance limits with measured values from the literature to put its chemistry in the context of other observed substellar objects. We mainly focus on the chemical abundances, leaving aside the temperature differences between the objects. We injected artificial planets in the data (T_eff = 1200 K, log (g) = 3.0) with different mass fractions of H₂O and CO andthen tried to detect both individually. The color map in Fig. 6 gives the value of the highest S/N detection for the two molecules. The solid contour line encloses the region where no molecule was detected with S∕N > 5, while the dashed line highlights the region where both molecules could be identified in the data (i.e., both molecules have S∕N > 5).

Previous works have exploited the capabilities of atmospheric retrievals (Benneke & Seager 2012; Lee et al. 2012; Line et al. 2012, 2013; Barstow et al. 2013, 2014; Waldmann et al. 2015b,a; Heng & Lyons 2016) to determine molecular abundances in planet atmospheres, mainly focusing on transiting hot Jupiters and low-mass field objects. Only a few of the directly imaged planets have been studied with atmospheric retrievals (e.g., HR 8799 bcde and 51 Eri b, Lee et al. 2013; Lavie et al. 2017; Wang et al. 2020a; Mollière et al. 2020; Whiteford et al. 2020), and their abundances are shown with blue squares in Fig 6. Furthermore, Todorov et al. (2016) applied atmospheric retrieval on the spectrum of the substellar companion κ And b (T_eff ~ 2000 K), concluding that clouds played a minor role in the atmosphere. They put constraints on the water abundance (blue dashed line), and a limit on the presence of methane, but were not able to say anything regarding CO.

Line et al. (2017) applied atmospheric retrieval analysis to 11 late T dwarfs (T_eff ~ 600−700 K) looking for physical and/or compositional trends. Late T dwarfs usually show deep spectral features due to H₂O absorption, which, combined with cloud-free atmospheres, allow us to probe a wide range of pressures and altitudes. They found constraints on H₂O, CH₄, and NH₃, and only upper limits on CO, CO₂, and H₂S consistent with non-detections (violet circles in Fig. 6). Given their water abundances (log (XH 2O) ~ −2.6 to − 2.3), they show chemical properties different from those of PDS70 b.

Regarding measured abundances in transiting hot Jupiters, the large uncertainties (often 2–3 orders of magnitude for H₂O and 5–6 for CO; e.g., Line et al. 2014) prevent us from making an informed comparison, as the chemical constraints are too loose. For this reason we do not report them in Fig. 6. Future missions like the James Webb Space Telescope (JWST, Gardner et al. 2006) and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL, Tinetti et al. 2016) will improve measurements of chemical composition of hot Jupiters and our understanding of their atmospheres.

Assuming that we are directly observing the atmosphere of the protoplanet, our analysis suggests that the chemical composition of PDS70 b is substantially different from those found for other substellar objects whose molecular composition is known, particularly under the assumption of a clear atmosphere.

Fig. 8 Cross-correlation function for injected planets emitting as the BT-Settl+ISM extinction (top) and DRIFT-PHOENIX (bottom) best-fit models proposed in Wang et al. (2021). CO is detected in both cases, while H2O only in the BT-Settl model. The S/N of each detection is shown (top right corner). The insets in the top left of each panel show the detection maps for H2O (blue) and CO (red).

5.5 Comparison with previous studies on PDS70 b

Wang et al. (2021) interpreted the stronger support of the data for atmospheric models than a blackbody spectrum as proof that we are indeed looking at the atmosphere of the planet. Our VLT/SINFONI data does not support this hypothesis, but we cannot discard it either. Dust extinction has a different impact on the detectability of absorption lines than on broad atmospheric features, as shown in Fig. 7: with increasing extinction the S/N of absorption lines is expected to get smaller, making a molecular mapping detection more difficult. This is shown in the bottom panel, where for A_V = 0, 5, 15 mag we show the absorption signature of the BT-Settl model reported in Wang et al. (2021) as the best-fit model for PDS70 b (T_eff = 1392 K, R_p = 1.96 R_J, log (g) = 3.83). Conversely, dust extinction may also modify the spectral slope at wavelengths where broad features are present, but the triangular K-band shape associated with H₂O bands remains visible even with A_V = 15 mag. The effect can be seen in the three top panels of Fig. 7, which display in gray the extincted BT-Settl model at the SINFONI resolution, while the colored circles show the spectrum down-sampled to the resolution of VLTI/GRAVITY. It is therefore plausible that dust extinction disproportionately affects line absorption features and broadband features, making molecular mapping in our SINFONI data unsuccessful in detecting molecules in the planet atmosphere. Hence, to better constrain the properties of the dust enshrouding PDS70 b, higher spectral resolution data able to either detect absorption lines or put tighter constraints on the amount of dust are required.

In addition to detecting molecular absorption at the blue and red end of the K band, Wang et al. (2021) tested several atmosphericmodels, augmenting them with dust extinction and possible circumplanetary disk contribution. Here we considered the two models that provided the strongest Bayesian evidence (BT-Settl and DRIFT-PHOENIX; Allard et al. 2012; Helling et al. 2008), and for these models we considered the model augmentation following the same criteria (ISM extinction and no augmentation, respectively). In the first phase we used the models (without ISM extinction for BT-Settl) to search for cross-correlation signals. Combining the signature of different molecules present in self-consistent models is expected to boost the S/N and reveal the detection of the atmosphere of the planet. This approach provided successful results in previous works (Hoeijmakers et al. 2018; Petrus et al. 2021). Nonetheless, using the pre-computed model grids indicated from the best fit of the SED did not allow the detection of molecular signals via cross-correlation. In the second step we inserted the models in the data as explained in Sect. 4.3, and then tried to detect H₂O and CO. In both cases, CO was expected to be detected, while H₂O only in the BT-Settl + ISM extinction model, even though with a low significance (S∕N ~ 4.5). The CCFs are presented in Fig. 8, while the insets show the molecular maps. Given the confidence of the expected detection (S/N reported in the legend of Fig. 8), which is inconsistent with the non-detections reported in Sect. 4, we can exclude that either of the models is able to properly describe the emission coming from PDS70 b.

However, even though the pure DRIFT-PHOENIX model maximizes the likelihood, provides high Bayesian evidence, and delivers a good fit to the spectral and photometric data points, the combination of log (g) and R_p reported in Table 4 of Wang et al. (2021) seems to indicate that PDS70 b has a mass of ~ 230 M_J (≈ 10–500 M_J considering the 95% interval). As pointed out in the same work, such a high mass is in clear contrast with the 10 M_J upper limit on the planet mass based on orbital fitting enforcing stability (Wang et al. 2021). Hence, although broad H₂O features were detected in the VLTI/GRAVITY spectrum and overall dusty models seem to explain the SED data points better, none of the presented atmospheric models is able to satisfy all the constraints arising from different analysis approaches (orbital fitting, molecular mapping, SED modeling). These self-consistent model grids have proven to be effective in describing more evolved substellar companion atmospheres, and are widely used by the community. Nevertheless, we argue that their efficacy in providing a full description of companion emission spectra does not apply to young protoplanets that are still forming like PDS70 b.

6 Summary and conclusions

Medium-resolution (R ≈ 5075) data from PDS70 b taken with VLT/SINFONI confirmed the complexity of constraining the chemical composition of forming giant planets. We applied molecular mapping to detect absorption features from the planet atmosphere in the K band, but we did not find any molecular signal. The non-detection of molecules in the protoplanet atmosphere, despite very deep data in terms of time on target, confirms that a dusty environment likely enshrouds the forming planet. The constraints are not tight enough to yield information about the physical properties of the dust. We showed that every self-consistent atmospheric model used in previous works to describe the SED contains spectral features that should have been detected in the SINFONI data. This inconsistency underlines how molecular mapping can provide crucial information to SED fitting with atmospheric models when studying young forming protoplanets, and both techniques need to be explored in more detail to understand the contradictory results. An approach involving multiple techniques, including molecular mapping, is probably required to combine the differentpieces composing planet formation and further constrain the physics and chemistry governing this process.

Higher resolution data with VLT/ERIS (Davies et al. 2018), VLT/CRIRES+ (Follert et al. 2014), or even with the proposed HiRISE (Vigan et al. 2018) might help shed light on the presence of molecules and spectral features of the object in the near-infrared, while observations at longer wavelengths, possibly with JWST/MIRI could help in better constraining the emission from the circumplanetary environment.

Acknowledgements

We would like to thank the anonymous referee, whose careful and constructive comments improved the quality of this manuscript. G.C. and S.P.Q. thank the Swiss National Science Foundation for financial support under grant number 200021_169131. 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. A.B. and S.P.Q. acknowledge the financial support of the SNSF. 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. I.S. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement No 694513. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1 and MA 9185/1-1) and from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. T.S. acknowledges the support from the ETH Zurich Postdoctoral Fellowship Program and the Netherlands Organisation for Scientific Research (NWO) through grant VI.Veni.202.230. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A Data selection

To verify the data quality we calculated S/N and detection limits for the individual observations, and for the combination of the two datasets. The procedure for calculating the S/N follows the prescriptions from Mawet et al. (2014), and the detection limits were estimated using PynPoint as explained in Stolker et al. (2020b) after the planet was removed from the data (see below). In short, after measuring at a given separation the azimuthal noise and self-subtraction effects on a high S/N injected planet, we estimated the contrast for a false positive fraction of 2.9 × 10⁻⁷, corresponding to 5σ confidence in the Gaussian limit, and including the penalty factor due to small number statistics (Mawet et al. 2014). In the top panel of Fig. A.1 we show the S/N for the first dataset, for the second dataset, and for the combination of the two as a function of the removed PCs. The highest S/N is always obtained when looking at the first dataset only (S∕N ~ 7.5−12), and the S/N estimated for the second dataset never reaches values >6. This dataset actually degrades the overall quality when combined with that from May (green curve). As the next step we generated contrast curves to estimate the detection limits for planets detectable with confidence level > 5 σ. The resulting detection limits, shown in the bottom panel of Fig. A.1, indicate that with the second dataset it is impossible to detect PDS70 b (shown as a red circle; see Appendix B for contrast estimation). When combined with the first dataset, the achieved detection limits are lower than when the first dataset alone is considered. Hence, the first of the two datasets was taken under generally better conditions, leading to a more robust detection of PDS70 b.

Fig. A.1 Evaluation of the two VLT/SINFONI datasets. Top: S/N of PDS70 b for the first dataset (orange dashed line), the second dataset (blue dotted line), and the combination of the two (green solid line) as a function of the number of removed principal components. Bottom: contrast curves obtained for the three considered datasets. The colors and linestyles are the same as in the top panel. PDS70 b separation and contrast estimates (see Appendix B) is shown as a red square. PDS70 c position and contrast are from Haffert et al. (2019).

Fig. A.2 Results after using our pipeline on an archival β Pic b dataset. Top: molecular maps for H2O and CO obtained from a dataset of β Pic published inHoeijmakers et al. (2018). North is pointing to the top, east to the left. Bottom: cross-correlation function extracted atthe planet position. The peaks representing the detection of water and carbon monoxide are clearly visible.

Appendix B Wavelength-collapsed cubes

A commonly used method to characterize the astrometry and photometry of companions in direct imaging is the injection of negative artificial planets (e.g., Stolker et al. 2019; Cugno et al. 2019b). This forward modeling technique allows us to measure therelative brightness and position of a companion despite self-subtraction, which can undermine direct measurements from the final PSF-subtracted images. Our approach was to iteratively minimize the residuals at the position of the companion after inserting a negative, rescaled copy of the PSF at its location and performing PSF-subtraction with PCA. In our case the PSF template was taken directly from each image, since the images were never saturated.

In order to sample posterior distribution for separation, position angle (PA), and contrast, we used the Markov chain Monte Carlo (MCMC) algorithm from PynPoint presented in Stolker et al. (2019). We let 200 walkers sample the parameter space with chains of 500 steps. For each walker, the first 100 steps were not considered (burn-in phase). The residuals were minimized in a circular aperture of diameter 2.0 FWHM (0.′′11) centered on the companion position removing six PCs at each step, as this is the setting that delivered the highest S/N (see top panel of Fig. A.1). Since uncertainties delivered by the MCMC algorithm were slightly asymmetric, we conservatively considered the largest of the two values. The results of the chain were examined, meaning that a negative signal with the retrieved contrast was injected in the images and after PSF subtraction the images were visually inspected to verify that no residuals from the companion are visible. The residual minimization was effective for all considered cases.

In order to quantify any potential bias in the measurements due to speckles and stellar PSF residuals, we injected artificial point sources at the same separation and with the same contrast as PDS70 b, at 360 equi-distant PAs (Stolker et al. 2020b). The same minimization algorithm was used to retrieve the parameter values and determine possible biases. The resulting offset between the injected and the retrieved values was calculated as the median of the distribution for each parameter. The uncertainties are conservatively estimated to be the maximum between the 16th and 84th percentiles of the offset distribution.

The final values for each parameter was obtained by adding the offset to the value estimated using the MCMC chain. Final uncertainties result from the sum in quadrature of the uncertainties of the MCMC chain and of the bias estimation. In addition, uncertainties on the plate scale and position angle obtained by Meshkat et al. (2015) were added in quadrature, although they were estimated for a companion at larger separation (Christiaens et al. 2019b). These are 0.4 mas in the separation and 0.5° in the PA. All values (and uncertainties) in terms of position and contrast for PDS70 b are listed in Table B.1.

The calculated contrast is in agreement with the contrast from IRDIS K1+K2 2016 data from Keppler et al. (2018), but our value is smaller (meaning the planet is measured to be brighter) than was estimated by Müller et al. (2018) and Stolker et al. (2020a) for a similar dataset taken in 2018. Finally, our contrast matches the value estimated by Christiaens et al. (2019b) from SINFONI H + K data. The position that we found is consistent with the past detections of PDS70 b (Keppler et al. 2018; Müller et al. 2018; Haffert et al. 2019; Wang et al. 2020b; Stolker et al. 2020a).

Appendix C Molecular mapping on β Pictoris b

In this section our aim is to prove that our data reduction and analysis is successful when molecules in the planet photosphere are abundant enough and their features are not damped by dust. We do this by analyzing a VLT/SINFONI dataset of β Pic already presented in Hoeijmakers et al. (2018). Compared to their work, we only used data from one of the two nights (2014 November 09, 50 min of time on target) when the star was in the field of view. The observation taken on 2014 October 09 (96 min) was carried out with β Pic A outside of the 0.′′8 × 0.′′8 field of view, and was therefore systematically different from our PDS70 data. We applied the same pipeline used on the PDS70 data (Sect. 3) and cross-correlated the residuals with templates containing H₂O and CO (Sect. 4.1). The resulting maps are shown in the top panels of Fig. A.2, with the molecular signature from the planet atmosphere correctly identified southwest from the star. The corresponding cross-correlation function is shown in the bottom panel of Fig. A.2 for both molecules. We obtained S/N${}_{{\text{H}}_2\text{O}}$ of ~ 10.4, and S∕N_CO ~ 9.2, lower than the values reported in Hoeijmakers et al. (2018), as we only consider one observing night. The analysis on the β Pic dataset shows that our pipeline is effective in identifying signals from molecular absorption in planetary atmospheres.

Appendix D Flux calibration

To estimate the PDS70 A spectrum we used the observation of the standard star HIP 81214 that was executed as part of the calibration. To obtain the spectrum of HIP 81214 we used the same pipeline applied to the PDS70 data, obtaining a cube covering the spectral range between 2.088 and 2.44 μm. The photometry in each channel was measured in an aperture of radius 2λ/D for both the HIP 81214 cube and the median-combination of the centered PDS70 cubes, and we calculated the ratio of these spectra in each channel. We used BT-NextGen models (Allard et al. 1997) to fit HIP 81214 photometric points from Gaia, HIPPARCOS, and 2MASS (Gaia Collaboration 2018; ESA 1997; Cutri et al. 2003) with species (Stolker et al. 2020b). The temperature and radius obtained from our fit agreed with the values from the Gaia catalog², while from photometric datapoints it was not possible to constrain metallicity and surface gravity. We assumed a constant 100% transmission as the instrumental response is taken into account in the calculation of the spectra from the SINFONI data, and any effect cancels out when calculating the ratio with the stellar spectra. This synthetic spectrum is multiplied by the flux ratio between theHIP 81214 and the PDS70 A spectra in order to obtain PDS70A flux density. We note that this does not correspond to the pure stellar spectrum, as substantial contribution from an inner disk are to be expected (Wang et al. 2021). Since the night was not photometric and aperture photometry on PDS70 A shows flux variations of up to 50%, we scaled the total spectrum of PDS70 A to the one calculated in Stolker et al. (2020a) in order to have continua that match each other. Their spectrum was obtained fitting a power-law function to the 2MASS JHK and WISE W1 and W2 magnitudes of the system. This step makes the stellar spectrum consistent with 2MASS K-band photometry.

References

All Tables

All Figures

	Fig. 1 PDS70 system; the b planet is shown within the solid circle. Six principal components were removed from the first dataset to detect the planet with S∕N  ~ 12. The dashed circle shows the expected position of PDS70 c (Mesa et al. 2019).
In the text

Fig. 2 Molecular templates after removal of the continua (H2O in blue, CO in red, and CH4 in green) and residuals of the data at the position of the planet (gray). No broad feature is visible and only individual absorption lines at the instrument resolution are cross-correlated. The small black dash in the top left corner of the first three panelsrepresents the filter size for continuum removal. Features smaller than this are kept in the residuals, while larger structures are removed as part of the continuum.

In the text

Fig. 3 Molecular maps of the PDS70 system and CCFs at the position of the b planet. Top row: H2O, CO, and CH4 molecular maps at RV ~ − 5 km s−1. The central star represents the position of PDS70 A, while the yellow circle indicates where a signal from PDS70 b should be found. Middle row: same as in the top row, but a planet with the molecular abundances indicated at the top of each panel has been inserted in the data and detected with S∕N ~ 5. Bottom row: CCF as a function of the radial velocity at the location of the planet. The orange lines represent for each molecule the CCF of the data, while the blue lines are the CCFs of the inserted artificial planets.

In the text

Fig. 4 Detection maps of H2O (left) and CO (right) for clear atmospheres as a function of mass fraction Xmol and effective temperature Teff, pT of the planet. The symbols represent measurements for other directly imaged companions (see left panel inset). The thick solid lines separate the region with S∕N > 5 from regions with S∕N < 5, i.e., they represent the molecular detection limit as a function of planet temperature. The dashed, dotted, and thin solid lines represent the molecular detection limits for different cloud configurations based on Mollière et al. (2020), as described in the main text (Pb = 10−3 bar is not shown). The temperature of PDS70 b ranges between 1200 and 1400 K (Stolker et al. 2020a; Wang et al. 2021). The effectivetemperature of β Pic b is well constrained as Teff = 1742 ± 16 K (Gravity Collaboration 2020) and the arrow in its symbol indicates that it is higher than the maximum Teff, pT considered here and should not be interpreted as a lower limit.

In the text

Fig. 5 Signal-to-noise ratio of the detection of H2O and CO as a function of extinction affecting the emission spectrum of the planet. Assuming PDS70 b has the same chemical composition as HR 8799 c, an additional extinction with AV ≈ 16−17 mag (AK ≈ 1.2 mag, using the extinction law from Wang & Chen 2019) is required in order not to detect molecules with S∕N > 5. The vertical dashed line shows the maximum value obtained by SED fitting from Wang et al. (2021). This extinction would have revealed molecules in the planet atmosphere with S∕N ~ 7–9.

In the text

Fig. 6 Comparison of the detectability of molecules in the SINFONI data with measured abundances from substellar objects. The blue squares represent directly imaged companions HR 8799 bce, β Pic b, and 51 Eri b (Lavie et al. 2017; Wang et al. 2020a; Mollière et al. 2020; Gravity Collaboration 2020; Whiteford et al. 2020); the blue lines show the measured range for H2O in κ And b (Todorov et al. 2016). The blue arrow indicates that the measured value is higher than the range of abundances shown in the plot. The violet circles represent a sample of T dwarfs analyzed in Line et al. (2017), and here the arrows refer to upper limits. The solid black line represents the region where no molecule would have been detected in our SINFONI data because the mass fractions were too low. The dashed black line encloses the region where, conversely, both molecules would have been detected.

In the text

Fig. 7 Effect of the extinction on the spectral measurement of the pseudo-continuum and on molecular mapping. Top: extincted best-fit BT-Settl models resampled at the SINFONI (gray) and at the GRAVITY (colored circles) spectral resolution. Even with high extinction, the broadband absorption features at the edges of the K-band spectrum can be identified. Bottom: spectral features at the SINFONI resolution after the removal of the continuum. Higher extinctions decrease the spectral signatures molecular mapping is looking for, impeding the detection of molecules.

In the text

	Fig. 8 Cross-correlation function for injected planets emitting as the BT-Settl+ISM extinction (top) and DRIFT-PHOENIX (bottom) best-fit models proposed in Wang et al. (2021). CO is detected in both cases, while H2O only in the BT-Settl model. The S/N of each detection is shown (top right corner). The insets in the top left of each panel show the detection maps for H2O (blue) and CO (red).
In the text

Fig. A.1 Evaluation of the two VLT/SINFONI datasets. Top: S/N of PDS70 b for the first dataset (orange dashed line), the second dataset (blue dotted line), and the combination of the two (green solid line) as a function of the number of removed principal components. Bottom: contrast curves obtained for the three considered datasets. The colors and linestyles are the same as in the top panel. PDS70 b separation and contrast estimates (see Appendix B) is shown as a red square. PDS70 c position and contrast are from Haffert et al. (2019).

In the text

	Fig. A.2 Results after using our pipeline on an archival β Pic b dataset. Top: molecular maps for H2O and CO obtained from a dataset of β Pic published inHoeijmakers et al. (2018). North is pointing to the top, east to the left. Bottom: cross-correlation function extracted atthe planet position. The peaks representing the detection of water and carbon monoxide are clearly visible.
In the text