© The Authors 2023

1 Introduction

The development of adaptive optics (AO) in recent years has allowed ground-based instruments such as the Gemini Planet Imager (GPI; Macintosh 2014) and the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE; Beuzit et al. 2019) to detect tens of substellar companions (e.g. Bowler et al. 2017; Janson et al. 2019; Bohn et al. 2020). Direct imaging allows spectroscopic observations of the companions in order to probe the properties and composition of their atmospheres. High-contrast imaging is sensitive to the thermal near-infrared (NIR) emission of recently formed giant planets and brown dwarfs. The young age of these objects makes them excellent targets for testing planet-formation theories (e.g. Spiegel & Burrows 2012; Mordasini et al. 2009a,b), by comparing their luminosity with evolutionary track predictions for models of differing initial entropy, such as hot- or cold-start models (e.g. Allard et al. 2012). In addition, advancements in the treatment of clouds in atmospheric models allow us to better characterise the detected companions (e.g. Baudino et al. 2015; Mollière et al. 2020; Carrión-González et al. 2020).

51 Eridani b is the first discovered planet by the GPI exoplanet survey (Macintosh et al. 2015). It was first characterised using both J- and H-band spectra from GPI, and L_P-band photometry from Keck/NIRC2. This young giant planet shows strong methane spectral signatures, an unusual feature in most directly imaged exoplanets. The planet orbits 51 Eridani A, a young F0IV star member of the β Pictoris moving group (Zuckerman et al. 2001; Bell et al. 2015). The latest estimate for the isochronal age of the system from Gaia EDR3 is ~10 Myr (Lee et al. 2022), a much younger age than the commonly adopted ~20 Myr (see Table 1 for more physical parameters of the star). Using photometry from the Transiting Exoplanet Survey Satellite (TESS), Sepulveda et al. (2022) recently determined that 51 Eri is a γ Doradus pulsator with a core rotation period of days. The star is part of a hierarchical triple system along with the M-dwarf binary GJ 3305AB, with the two separated by ~2000au (Feigelson et al. 2006; Kasper et al. 2007). 51 Eri is located at 29.90±0.06 pc, as derived from the precise parallax measurement by the Gaia mission (Gaia Collaboration 2021). From 24 μm Spitzer observations, 51 Eri is known to have an infrared (IR) excess (Rebull et al. 2008), and a debris disk was detected using Herschel observations at 70 and 100 μm with a very low IR fractional luminosity of L_ir/L_* = 2.3 × 10⁻⁶ and a lower limit on the inner radius of 82 au (Riviere-Marichalar et al. 2014). From WISE observations, Patel et al. (2014) report a warm disk (T ~ 180 K) at a radius of 5.5 au assuming blackbody radiation. Therefore, it is likely that the architecture of 51 Eri could resemble that of our Solar System with a two-belt debris disk.

The planet 51 Eri b was confirmed to be bound to the system in a follow-up paper by De Rosa et al. (2015). In De Rosa et al. (2020), the authors presented a revised version of the orbital parameters of the planet (inclination deg, semi-major axis au, and orbital period yr assuming a mass of 1.75 M_⊙ for the host star) from GPI observations, which they found to be consistent with the parameters derived by Maire et al. (2019) from SPHERE observations. The inclination found in these studies suggests that the planet is not coplanar with the binary GJ 3305AB (i = 92.1 ± 0.2, Montet et al. 2015). The discovery paper used the luminosity and age of the system to derive a mass estimate for 51 Eri b of 2-12 M_Jup for a cold-start formation scenario (Macintosh et al. 2015). Later estimates reported by Samland et al. (2017; hereafter SAM17) ranged from 2.4 to 5 M_Jup for hot-start, and up to 12 M_Jup for warm-start models. De Rosa et al. (2020) showed that a dynamical mass measurement of the planet with Gaia may be possible if the mass of the planet is M ≳ 4 M_Jup. Recently, Müller & Helled (2021) reported mass and metallicity estimates for 51 Eri b derived from synthetic cooling tracks and the planet's luminosity. Assuming an age range of 17–23 Myr and a hot-start formation scenario, they obtain a mass of M = 2.3 M_Jup and a metallicity of [Fe/H] = 0.11. Another study by Dupuy et al. (2022) presented an upper limit for the mass of the planet of M < 11 M_Jup at 2σ using cross-calibration of HIPPARCOS and Gaia EDR3 astrometry. They also revise the luminosity of 51 Eri b using a photometric approach and find log(L_bol/L_⊙) = −5.5 ± 0.2 dex. Additionally, they derived a lower limit on the initial specific entropy of the planet, which rules out cold-start formation scenarios.

In addition to the discovery paper, two atmospheric analyses have been carried out to characterise 51 Eri b. By combining SPHERE/IFS YJ and YH, and GPI H spectra, along with photometry from SPHERE (broad-band H, H23 and K12) and Keck/NIRC2 (L_P), SAM17 found the atmosphere to be cloudy. They report T_eff = 760 ± 20 K, , log g = 4.26 ± 0.25 dex, [Fe/H] = 1.0 ± 0.1 dex, and for their best-fit model. On the other hand, Rajan et al. (2017) used updated GPI J and H spectra and updated Keck/NIRC2 L_P photometry from the discovery paper, along with new Keck/NIRC2 M_S photometry and determined the atmosphere to be partially cloudy. Their best-fit model yielded T_eff = 605−737 K, [Fe/H] = 1.0, and log g = 3.5–4.0 dex. The three studies made use of self-consistent models and differ mainly in the degree of cloudiness of the atmosphere.

In this work, we present new NIR spectro-photometric observations of 51 Eri b obtained using VlT/SPHERE as part of the SHINE survey (Desidera et al. 2021; Langlois et al. 2021; Vigan et al. 2021). These observations were carried out as a follow up to the ones presented in Samland et al. (2017) and have the highest signal-to-noise ratio (S/N) achieved so far (S/N ~23 for K1). We use the radiative transfer code petitRADTRANS to model the atmospheric spectrum of the planet. Additionally, we attempt to reproduce the results in SAM17 using a retrieval approach.

The paper is structured as follows: in Sect. 2 we describe our spectro-photometric observations. In Sect. 3, we detail our data-reduction and spectrum-extraction procedures, as well as the derived detection limits in planet contrast and mass. The description of the atmospheric retrieval runs is detailed in Sect. 4. We present a detailed analysis of certain parameters of the planet and a discussion of the results in Sect. 5, and finally, we summarise our results in Sect. 6.

2 Observations

New data for 51 Eri were obtained with the VLT/SPHERE high-contrast instrument (Beuzit et al. 2019) within the SHINE (SpHere INfrared survey for Exoplanets, Chauvin et al. 2017) Guaranteed Time Observations (GTO) program on the night of September 28, 2017. These observations were carried out in the IRDIS_EXT mode, which combines IRDIS in dual-band imaging mode (DBI; Vigan et al. 2010) in the K1K2 (K1 = 2.110 μm and K2 = 2.251 μm) filters with IFS (Claudi et al. 2008) in the YH spectral bands (between 0.95 and 1.65 μm with spectral resolution R ~ 33). An apodised pupil Lyot coronagraph with a focal mask diameter of 185 marcsec, was used for the observations (Carbillet et al. 2011). In order to reduce residual speckle noise, the observations were carried out close to meridian passage using the pupil-stabilised mode, which allows the use of ADI post-processing (Marois et al. 2006).

To calibrate the flux and centre of the images, unsaturated non-coronagraphic images (hereafter referred to as the point spread function (PSF)), as well as coronagraphic images with the deformable mirror (DM) waffle mode (see Langlois et al. 2013, for more details on this mode) were acquired at the beginning and at the end of the observing sequence. The waffle mode generates four artificial replicas of the star in a 'cross' pattern, commonly known as satellite spots. These spots are used to measure the star's position at the centre of the pattern. In order to minimise the uncertainties in the frame centring and the astro-metric error, and monitor the photometric stability throughout the sequence, the science frames were also obtained using this mode. Finally, night-time sky images were acquired to estimate the background level in the science frames. The pixel scale and the true north (or north angle offset) were obtained using astrometric calibrators, as described in Maire et al. (2016). The astrometric calibration of the data set was done using an observation of an astrometric field performed during the observing run. The usual calibration images (i.e. flat-field, bias, and spectral calibration) were obtained by the internal calibration hardware of the instrument during daytime after the observations. The observations were conducted in favourable overall conditions (see Table 2), except for the presence of clouds near the end of the observing sequence.

3 Data reduction and spectrum extraction

3.1 Data reduction

The data were reduced with the SPHERE Data Center pipeline (Delorme et al. 2017) using the SPHERE Data Reduction Handling (DRH) software (version 15.0; Pavlov et al. 2008). This basic reduction consists in performing the usual sky-background subtraction, flat fielding, bad-pixel identification and interpolation, star-centring corrections and, for IFS, the calibration of the wavelengths and of the cross-talk between spectral channels. We then removed poor-quality frames where a significant drop in stellar flux is detected from the photometry of satellite spots. Because of interfering cirrus clouds near the end of the observing sequence, only the first 140 of the total of 154 IFS frames were used for the post-processing analysis. We also discarded the second PSF for the flux normalisation and only used the first PSF frames. The conditions were very stable from the start of observations for the first PSF frames until the onset of cirrus clouds. Additionally, we tested different SPHERE data-reduction recipes and pipelines (see Appendix A for a detailed discussion on how they compare).

For the post-processing data analysis, we used the ANgular Differential OptiMal Exoplanet Detection Algorithm (ANDROMEDA; Cantalloube et al. 2015), which utilises angular differential imaging (ADI) and an inverse problem approach based on a maximum-likelihood estimator (Mugnier et al. 2009) to search for companion candidates. It carries out a pair-wise subtraction of frames at different rotation angles and performs a cross-correlation of the signature that a point source would leave in the residual image (see Cantalloube et al. 2015, for more details). The main outputs from ANDROMEDA are 2D maps of: the estimated contrast of detected point sources at every location in the image; the corresponding standard deviation on this contrast; and the S/N ratio obtained dividing the contrast by the standard deviation. Figure 1 shows the resulting S/N detection map from ANDROMEDA for our IFS data. Other than 51 Eri b, we detect no additional point sources.

Additionally, we performed the reduction with TRAP, a temporal, non-local systematic-uncertainty-modelling algorithm, in order to look for point sources at small separations (see Samland et al. 2021). Regarding the S/N, with ANDROMEDA we get S/N = 23.31 for K1 and S/N = 2.82 for K2, while with TRAP we get S/N = 18.03 for K1 and S/N = 3.52 for K2. In both cases, there is an improvement in S/N from SAM17 (7.46 and 1.26 for K1 and K2, respectively using ANDROMEDA). On the other hand, the contrast limits are improved with the TRAP reduction as discussed below in Sect. 3.3. Due to the higher S/N achieved in K1 , and for consistency throughout the paper, we decided to use the results of ANDROMEDA for the following steps.

Fig. 1 Median combined S/N detection map from ANDROMEDA based on SPHERE/IFS data. The circle indicates the position of 51 Eri b. The azimuthal wings around the planet signal is the characteristic planet signature that ANDROMEDA fits for ADI data.

3.2 Spectrum calibration

To construct the spectrum of 51 Eri b, we multiplied the planet contrast at each wavelength from ANDROMEDA by a template spectrum of the host star. This template spectrum was obtained as follows: we used a model stellar spectrum from the BT-NextGen library (Allard et al. 2012) with T_eff = 7200 K, log g = 4.0dex, [Fe/H] = 0.0 dex, and no alpha enhancement (overabundance of He with respect to metallicity, [α/Fe]), as these parameters are the closest to the ones determined from highresolution spectra for 51 Eri (T_eff = 7256 K, log g = 4.13 dex, and [Fe/H] = 0.0 dex; Prugniel et al. 2007). We then fit this model spectrum to the spectral energy distribution (SED) of the star using the χ² minimisation in the Virtual Observatory SED Analyzer (Bayo et al. 2008) to obtain the flux scaling factor to account for the distance of 51 Eri. The SED was built with photometry from Tycho B_T, V_T (Hoeg et al. 1997), WISE W3 (Cutri et al. 2013), Johnson U, V, B (Mermilliod 2006), and IRAS 12 μm (Helou & Walker 1988). Finally, we scaled the model spectrum to the resolution of our IFS data using the Python function SpectRes¹.

A similar procedure along with the respective transmission curves for the filters was used to obtain updated IRDIS photometry for K1 and K2. As the planet is not detected significantly in the K2 filter, we applied forced photometry with ANDROMEDA. This consists in performing a photometric measurement in the K2 images at the position of the planet in the K1 frames. To estimate the uncertainties in the flux of both the spectrum and the photometry, we propagated the standard deviation obtained with ANDROMEDA for the IFS and IRDIS data, respectively. Our results are shown in Table 3.

Figure 2 shows the spectrum of 51 Eri b using our IFS data, along with the K1K2 photometric points. For comparison, we have overplotted the IFS/YH spectrum and the K1K2 photometry presented in SAM17, as well as the updated GPI J and H spectra and the K1 and K2 spectra from Rajan et al. (2017). From the latter, we see that the GPI K1 and K2 spectra are consistent with our derived K1K2 SPHERE photometry. It is also worth noticing that even with the revised version of the GPI J and H spectra, discrepancies with the SPHERE data still persist. Furthermore, we now see a slight discrepancy in the H-band spectrum as well, although mostly within the error bars. As explained in detail in SAM17, the discrepancies can most likely be attributed to systematic uncertainties. To further highlight the differences between the two IFS/YH spectra, in Fig. 3 we present the relative error as a function of wavelength. For both figures, we removed the spectral channels that were not used in SAM17 (1.14 and 1.41 μm) to allow for a better comparison. Overall, our data exhibit a lower relative error – except for the telluric H₂O absorption band regions around 1.1 and 1.35 − 1.4 μm (see Appendix B) –, which results in a higher S/N ratio than the data in SAM17.

Fig. 2 Newly obtained 51 Eri b YH spectrum and K1K2 photometry. Overplotted for comparison are the YH spectrum and K1K2 photometry from Samland et al. (2017), as well as the updated GPI J and H spectra and the K1 and K2 spectra from Rajan et al. (2017).

3.3 Detection limits

We used both ANDROMEDA and TRAP to derive 5σ contrast curves for our IRDIS and IFS data. For the IRDIS/K1–K2 bands, the analysis setup was SDI+ANDROMEDA assuming no planet flux in K2. For the IFS–YH, we also used SDI+ANDROMEDA assuming a T5 spectral template for putative planets. To convert contrast to mass limits we used the evolutionary tracks of Baraffe et al. (2003) along with the atmosphere model of Baraffe et al. (2015). For the star, we used the 2MASS JHK magnitudes (Cutri et al. 2003), the L_P magnitude from Macintosh et al. (2015), the newly calculated distance from Gaia Collaboration (2021), and the latest age estimate from Lee et al. (2022). The resulting detection limits are shown in Fig. 4, we show the curves from ANDROMEDA and overplotted in dashed lines the curves obtained with TRAP for the IFS data. The mass limit for the IFS data is cut to the lowest mass computed by the model grid for both pipelines. The Κ1 mass curve reaches smaller values than the K2 mass curve because the K2 filter matches a methane absorption band that strongly dims the flux of cold giant planets. Our detection limits are corrected for the coronagraphic transmission (Boccaletti et al. 2018) and for small sample statistics (Mawet et al. 2014). Our results from ANDROMEDA tell us that the IFS data are sensitive to planets more massive than ~1 M_Jup beyond ~6 au, and more massive than 2 M_Jup beyond 4.5 au. While with TRAP we see an additional improvement in the sensitivity of planets down to 2 M_Jup at 3 au. We therefore achieve an improvement of about 2 M_Jup in sensitivity over previous studies (e.g. SAM 17).

Fig. 3 Flux uncertainty over wavelength as a function of wavelength for the YH spectra of 51 Eri b from this work compared to the one presented in Samland et al. (2017).

4 Retrieval analysis

We used the ID radiative transfer code petitRADTRANS² (pRT; Mollière et al. 2019) in its scattering implementation (Mollière et al. 2020) in combination with PyMultiNest³ (Buchner 2014) to derive the posterior distributions of the thermal structure, chemical composition, and cloud properties of 51 Eri b. The code takes as an input the spectra of the planet along with prior distributions for the metallicity, C/O ratio, log g, radius, a list of molecules to be included, quench pressure, and cloud parameters such as f_sed, K_zz, and log τ_cloud.

Fig. 4 5σ detection limits. Planet contrast (top) with respect to the star as a function of separation to the star, and planet mass detection limits (bottom) as a function of the separation from the star for the SPHERE/IRDIS_EXT data of 51 Eri b. The mass limit for the IFS data is cut to the lowest mass computed by the model grid (1 MJup) for both ANDROMEDA and TRAP.

4.1 Modelling setup

4.1.1 petitRADTRANS setup

The atmospheric retrieval modelling setup of pRT is described in detail in Mollière et al. (2020). In the following, we describe the parameters used in our retrievals:

Retrieved parameters

The following parameters are of prime interest in the retrieval, and we assign a flat prior to each; see Table 4.

Additional parameters

An important additional parameter is the effective temperature of the planet, T_eff. This is not an explicit input (and thus retrievable) parameter of pRT, but instead has to be calculated by generating a second spectrum for a given set of parameters that covers a wide spectral range in order to estimate the bolometric flux of the planet. Due to the required large wavelength coverage, this can be quite time consuming and is usually only carried out on a subset of the equal-weighted posterior distribution.

In addition to these parameters, there are a number of internal 'nuisance' parameters to which prior ranges are also assigned. These are a set of connecting temperatures (t₁, t₂, t₃), an internal temperature t_int, and the two parameters log δ and α for the optical depth model τ = δΡ^α. These parameters are later used to determine the atmospheric pressure-temperature (P–T) profile. This P–T model is described in Mollière et al. (2020).

Clouds

From our first retrieval runs on the new data, we observed that the best-fit solutions tended to be non-cloudy. According to previously reported results, the photosphere of 51 Eri b is thought to be at least partially cloudy (e.g. Samland et al. 2017; Rajan et al. 2017). We decided to implement the parameter log τ_cloud to enforce clouds in the retrieval, which represents the logarithm of the cloud optical depth at the photospheric region of the clear atmosphere (τ = 1 ). These are wavelength-averaged optical depths estimated over the wavelength range of the retrievals. By changing the range of the prior of log τ_cloud we were able to enforce clouds in the photospheric region. In order to verify whether cloud cover is actually present, we split our models into 'nominal' and 'enforced clouds' models. These are the same models, differing only in the prior range of log τ_cloud, [−7, 3] for the 'nominal', and [−1, 3] for 'enforced clouds' (cf. Sect. 4).

For the cloud species, we use Na₂S and KCl, which according to Morley et al. (2012) are the most important species at the previously estimated effective temperature of 51 Eri b (700–750 K).

Molecular species

The species contributing to the line opacities in our model are CO, H₂O, CH₄, NH₃, CO₂, H₂S, Na, K, PH₃, VO, TiO, and FeH. In addition, we include H₂ and He as species contributing to both Rayleigh scattering and collision-induced absorption. The species are retrieved under equilibrium chemistry assumptions and including quenching pressure.

Fig. 5 petitRADTRANS best-fit spectrum of 51 Eri b for the 'nominal' retrieval run (shown in black) on our new SPHERE spectro-photometric data (teal circles and purple squares) along with the photometric points included in SAM17 and (Rajan et al. 2017; shown as squares). The photometric points describe the average flux in the respective filter, and the x-error bar represents the filter widths. The 34 randomly drawn samples from the posterior probability distribution are shown in grey in order to show the spread of model parameter combinations to fit the data. Residuals in multiples of 1 σ uncertainties of the data for the best-fit model are shown below.

4.1.2 MultiNest setup

To fit model spectra to the data by sampling the posterior probability, we used the nested sampling library PyMultiNest (Buchner 2014), which is in turn based on MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2019). Nested sampling (Skilling 2004, 2006) is a powerful method, which in contrast to MCMC methods is better at exploring the parameter space and is less prone to falling into a local minimum.

Our derived model parameters are chosen to be the median of the marginalised, equal-weighted posterior distribution, and the uncertainties quoted refer to the 16th and 84th percentiles of these posteriors. We note that the parameters used to generate actual best-fit spectra are generally different from the medians mentioned above and correspond to the highest log-likelihood. We use importance nested sampling with flat priors (see Sect. 4.1.1) and 4000 initial live points to sufficiently cover the parameter space. To ensure a high sample-acceptance fraction, we use MultiNest's constant efficiency mode and a sampling efficiency of 0.05.

4.1.3 Input data

Data set

Our YH-spectrum, as all SPHERE/IFS spectra, comprises 39 channels, all of which were fed into the retrieval process. In addition, we use the IRDIS H2/H3 photometric filters and broad-band H from SAM17, and K1/K2 from this work. Also included are the L_P and M_S Keck/NIRC2 data points from Rajan et al. (2017). We do not include GPI spectra in these retrievals due to the difficulty in quantifying systematic differences in the calibration of different instruments. A study combining data from SPHERE and K-band GRAVITY to expand the wavelength coverage using the retrieval approach is reserved for future work (Samland et al., in prep.).

As described in SAM17 and originally in Greco & Brandt (2016), the spectral covariance of the residual speckle noise should be taken into account when computing the likelihood of a model that matches the data from an IFS-type instrument. We computed the correlation matrix Ψ for our spectrum in the same way as SAM17. Consequently, in the log-likelihood computation in the retrieval code, we used

(1)

instead of the simple

(2)

In both equations, S represents the observed spectrum, and F the model spectrum. See Appendix C for more information. The SPHERE YH spectrum and the covariance matrix is available in Vizier/CDS.

4.2 Results and comparison with older data

4.2.1 Retrieval results

Figure 5 shows the best-fit spectra for the nominal model with the new data. Overall, the best-fit model is able to reproduce our YH spectrum, as well as the H2H3, K1K2, and broad-band H photometric points within the uncertainties. However, we note that the J-band flux from the best-fit model is lower than in our data by as much as 3–4 σ, which also occurs in the 'enforced clouds' retrieval (see Fig. D.1). This discrepancy appears to be caused by systematic uncertainties either in the opacities provided to the model, or, more likely, in the wavelength calibration strategy applied by the reduction pipeline that we are not able to remove. In either case, by using the correlation matrix, these data points do not have a meaningful impact on the spectrum fit. Additionally, the M_S photometric point is not fitted in either the 'nominal' or the 'enforced clouds' retrievals. Given the large photometric uncertainties in the M_S data, the best-fit photometry model lies within 2σ of the data. As further discussed in Sect. 5, future mid-IR instruments might help to better constrain this photometric band.

Figure 6 shows the corresponding corner plot for the 'nominal' model with the posterior PDFs of selected parameters. The best-fit spectra along with the corner plots for the 'enforced clouds' retrieval with our new data are shown in Appendix D. The median of the posterior probability distribution, and the uncertainties representing a 1σ uncertainty range for selected parameters (chosen due to their physical relevance) of our four sets of retrievals are shown in Table 5.

Fig. 6 Corner plot of the posterior PDFs of the 'nominal' retrieval run on the new data set.

4.2.2 Atmospheric retrieval on previous data

In order to determine the extent to which the use of a retrieval method alone impacts the outcome in terms of the derived atmospheric parameters with respect to the SAM17 grid of self-consistent models, we ran an additional set of pRT retrievals (nominal and enforcing clouds) on the original data used by SAM17. These attempts can also tell us the extent to which the outcomes are a result of the better quality of the SPHERE 2017 data. We included all the data cited in the work of Samland et al., that is, the YJ and the YH spectra, and the H2H3, K1K2, L_P, and broad-band H photometric points, all resulting from 2015 and 2016 observations, as well as the GPI spectrum published in the discovery paper (Macintosh et al. 2015). We note that SAM17 used an SDI+ANDROMEDA data reduction, which resulted in the reference channels of their spectra being unusable in the model comparison. In addition, they masked out a number of channels due to low S/Ns. We took these considerations into account in our reproduction attempts. We do not include the revised version of the J and H spectra, or the K1K2 spectra from Rajan et al. (2017), for consistency in our attempt to reproduce the results of SAM17. The results of the nominal and enforced cloud retrievals are shown in Appendix D.

When comparing the parameter values for the 'SAM17 nominal' retrieval in Table 5 to the ones in Table 5 of SAM17 (top row, 'PTC-C', shown in the last row of Table 5), we notice that we find a significantly lower metallicity, and very little indication of the presence of clouds. Other parameters differ as well, but are mostly within the 16th and 84th percentile limits. We discuss this issue, and in particular our attempt to determine whether the object is cloudy or not, in depth in Sect. 5.6.

5 Discussion

In the following, we present a discussion of the results obtained for certain parameters of interest of 51 Eri b using the retrieval approach. Next we discuss the differences between using atmospheric retrievals and self-consistent models.

5.1 C/O ratio

The atmospheric carbon-to-oxygen ratio has been linked to the formation scenario of exoplanets (Öberg et al. 2011). The different condensation temperatures of water (H₂O), carbon oxide (CO), and carbon dioxide (CO₂) govern the location of their 'ice lines' in different parts of the protoplanetary disk, which results in different values of C/O throughout the disk. A difference in the C/O ratio of a gaseous giant planet compared to the C/O ratio of its host star can provide information about the formation of the planet, depending on whether their C/O ratio is super- or substellar; however, see Mollière et al. (2022) on how challenging it is, in practice, to go from C/O to formation.

We find the C/O ratio of 51 Eri b to be consistent along retrieval runs for the same data set (C/O ~ 0.4 ± 0.08 for the new data and C/O ~ 0.8 ± 0.17 for the data in SAM17), respectively. Compared to the solar C/O ratio assumed by SAM17 (i.e. C/O = 0.55), all our retrieved values within error bars differ by ~0.1. However, no definitive conclusions can be drawn because the value for the C/O ratio of the star remains undetermined. The lower flux from the M_S photometric point hints to the presence of carbon monoxide (CO), but further observations would be needed to use it to constrain the C/O ratio of 51 Eri b.

5.2 [Fe/H]

In our retrievals using the data in SAM17, we obtain a metallicity [Fe/H] in the range of −0.09 to 0.30, including uncertainties, which differs from the results of SAM17, who find [Fe/H] = 1.0±0.1. However, when comparing with retrieval results for two benchmark brown dwarfs, these latter authors find that they tend to fall in the lower end of the estimated metallicity range for the host star. In our case, we observe a similar behaviour, even for the new value we obtain of [Fe/H] ~ 0.26, which is slightly super-solar. As discussed at length in SAM17, the derived metallicity strongly depends on the K-band flux, and indeed our retrieval tends to slightly underfit the K1-flux, whereas it was overfitted in SAM17. When artificially enhancing the importance of the K1-flux point by lowering its uncertainty by a factor of ten, the resulting fit for the SAM17 data shows an increased metallicity , and a higher . These values are closer to those of SAM17, but still do not agree. The initial and the remaining difference indicate that parameters derived from retrievals can differ significantly from self-consistent models.

5.3 Clouds and log τ_clouds

Our initial retrievals always resulted in the best-fit models tending towards non-cloudy solutions. Based on previous results, which suggested that the atmosphere of 51 Eri b is at least partially cloudy, we included an additional prior to enforce clouds and check the robustness of the retrieval: log τ_clouds. However, even for the 'enforced clouds' retrievals, the value for log τ_clouds always tended towards the lower limit of the prior (i.e. towards cloud-free solutions; see Table 5). Both f_sed and K_zz describe the cloud properties (Ackerman & Marley 2001). Our best-fit values for fsed in the 'enforced clouds' retrievals are in agreement with previously reported f_sed values for brown dwarfs in a similar temperature range as 51 Eri b (e.g. GJ 758 B, GJ 504 b; Vigan et al. 2016; Skemer et al. 2016, respectively). We also note that our values are within the ranges of f_sed found by SAM17. A higher f_sed corresponds to vertically thinner clouds with larger particle sizes. Our derived values for K_zz are within the assumed values in both Macintosh et al. (2015) and SAM17.

In Fig. 7, we show the P–T profiles resulting from the retrieval runs along with the 1σ, 2σ, and 3σ confidence intervals for our four different cases. Overplotted are the corresponding self-consistent P–T profiles obtained from petitCODE when feeding in the best-fit parameters of the retrievals. We note that the self-consistent P–T profiles are less isothermal than the retrieval ones in all cases, following the characteristic atmospheric temperature gradient for models in radiative-convective equilibrium. In this scenario, the only way to reproduce the observed low NIR fluxes is to add clouds, which contradicts our results with the retrieval models. This discrepancy has been observed in other studies (see Sect. 5.6 for an in-depth discussion), and has yet to be resolved, which means we cannot draw conclusions about the cloudiness of exoplanet atmospheres.

Fig. 7 Retrieved pressure-temperature profiles in white, with confidence intervals (1σ, 2σ, and 3σ) for our four different retrievals shown in Table 5. Overplotted are the resulting self-consistent P–T structures in pink.

5.4 Radius and surface gravity

The values we obtain for the radius and surface gravity are in agreement with previous results for the data in SAM17 within uncertainties. As seen in Figs. 6 and D.2, the retrieval finds the best-fitting models towards the lower R_p prior boundary for the new data, while this is not the case for the retrievals using the data in SAM17 (Figs. D.4 and D.6). We attribute these variations to the different input data. Nevertheless, all values are consistent with the radius of Jupiter within uncertainties, which according to planetary evolution models are slightly smaller than expected for the possible age of the system (i.e. 10–20 Myr, ~ 1–1.2 R_Jup; Mordasini et al. 2012).

5.5 Mass

To derive an estimate of the mass of 51 Eri b, we used the posterior distribution for the surface gravity and radius of our 'nominal' model, and the relation M = g/g_Jup × (R/R_Jup)², where g_Jup = 24.79 m s⁻², and R_Jup = 6.99 × 10⁷ m are the surface gravity and the volumetric mean radius of Jupiter, respectively. We obtain a mass of M = 3.9 ± 0.4 M_Jup. Additionally, we used the photometric values for the K1K2 bands along with the evolutionary models for extrasolar giant planets from Baraffe et al. (2003). We used both estimates of the age of the system and we obtain a value of M = 2.4 M_Jup at 10 Myr, and M = 2.6 M_Jup at 20 Myr. All mass values of the planet are smaller than the value found by SAM17 , but the formation scenario can strongly impact the mass (and the surface gravity) of the planet, and small masses are possible within the 'hot-' and 'warm-start' scenarios.

5.6 Atmospheric retrievals versus self-consistent models

Atmospheric retrievals are now a commonly used tool for fitting the spectra of exoplanets (for a recent review, see Madhusudhan 2019). The general idea of retrievals is that an atmospheric forward model calculates planetary spectra based on a parameterised atmospheric structure, which is then compared to the data. This parameterisation is key, because finding the atmospheric temperature, composition, and cloud structure in a physically self-consistent way is a numerically expensive step. Self-consistent models typically assume that the atmosphere is in radiative-convective equilibrium, and couple this assumption to a model solving for the chemical composition and cloud structure of the atmosphere in an iterative fashion (see e.g. Hubeny 2017, for a review). Furthermore, many processes, especially those connected to clouds, are not sufficiently well understood. If the underlying assumptions in the self-consistent model are incorrect, this may lead to unfounded conclusions about the properties and parameters of the atmospheric model. At the same time, self-consistent models allow us to compare our complete physical understanding against what the data show. Conversely, retrievals attempt to constrain the atmospheric structure mostly from the data alone (if uninformative priors are chosen), circumventing these issues. However, this requires data of sufficiently high S/N and spectral coverage.

Retrievals have proven challenging so far for cloudy, directly imaged planets. This is connected to a degeneracy, or at least a correlation: a cloud generally reddens the spectrum by hiding the deep hot regions of the atmosphere from view. If the cloud was not present, too much flux would escape from the atmosphere, especially in the opacity windows in the Y, J, and H bands. Moreover, an atmospheric model in radiative-convective equilibrium generally results in a large atmospheric temperature gradient, such that the only way to reproduce the red spectral energy distributions (SEDs) of cloudy planets and brown dwarfs is to add clouds. In Tremblin et al. (2015, 2016, 2017, 2019), atmospheric instabilities that decrease the atmospheric temperature gradient when compared to the equilibrium solution have been suggested to mimic the effect of clouds in the NIR. This can be easily understood: if the deep atmosphere is less hot, no clouds are required to lead to a reduced flux in the YJH-bands. For retrievals with a flexible atmospheric temperature and cloud parameterisation, this may result in atmospheric structures that are cloud-free and too isothermal when compared to classical self-consistent models. Due to the ease with which red exoplanet spectra can be fit with overly isothermal temperature profiles, it is not surprising that a retrieval can fall into this 'trap': different temperature structure realisations are explored during a retrieval, and if the cloud model is not the 'perfect' one, and leads to a slightly worse fit, there is no reason for the retrieval to add a cloud of appreciable opacity to the atmosphere.

This cloud-free retrieval problem appears to be emerging in recent studies (Mollière et al. 2020; Zhang et al. 2021; Kammerer et al. 2021), and we also see it occurring here. Options to enforce a more cloudy solution may include making the atmospheric temperature parameterisation less flexible, such that a strong temperature gradient is always present, which then needs to be corrected for by adding a cloud (which was identified as a workaround in Kammerer et al. 2021 when using the atmospheric retrieval code ATMO; Tremblin et al. 2015, 2016). The danger is that such temperature profiles may not be sufficiently complex to capture the atmospheric temperature structure even for a cloud-free planet, leading to potentially biased results for the atmospheric properties. Enforcing a minimum cloudiness in the atmosphere via a cloud optical depth prior, as we attempted to do here, may be another option, but we remind the reader that we tried to make the retrievals as cloud-free and isothermal as allowed by the retrieval. A truly satisfactory method for approaching this cloud-P-T correlation has yet to be found.

It is conceivable that these atmospheres are truly mostly cloud-free, and that this is why the retrievals tend towards these solutions, as also argued in the Tremblin et al. papers. However, even synthetic cloudy spectra are retrieved as cloud-free if the cloud model is modified between the synthetic observation and retrieval step (Mollière et al. 2020), such that the absence of clouds in the retrievals cannot be regarded as proof of the absence of clouds in the atmospheres of real planets. A promising avenue is presented in Burningham et al. (2021). The authors found that adding mid-IR data tracing silicate cloud absorption features at 10 μm leads to definitely cloudy solutions, as well as temperature structures that are less isothermal compared to the retrievals that neglected the mid-IR data (Burningham et al. 2017). The James Webb Space Telescope's (JWST) observations of cloudy exoplanets and brown dwarfs with its mid-IR instrument MIRI (Wright et al. 2004) may therefore hold great potential to resolve – at least partially – the cloud-temperature gradient correlation. A remaining challenge is that the 10 μm region probes lower pressures than the NIR (YJH bands), and therefore probing the silicate feature could merely help to constrain the P–T profile in the upper atmosphere. Indeed, Burningham et al. (2021) need to include a second deep cloud to produce the NIR reddening. However, we note that the above could be a general solution for planets but could not be applied to 51 Eri b, because the planet is too cold to have silicate clouds, and too close to the host star to be observed with JWST/MIRI. A possible solution for cool, closer-in giant planets like 51 Eri b could be to use the Mid-infrared ELT Imager and Spectrograph (METIS, Quanz et al. 2015) to expand the wavelength coverage (up to 13 μm) and search for the silicate feature. Finally, we advocate that results from atmospheric retrievals should never be discussed in isolation, but instead along with the results from self-consistent models. This can be done in different steps:

1. For example, via comparisons to self-consistent atmospheric structures obtained from using the retrieval's best-fit parameters for the atmospheric composition, gravity, effective temperature, and cloud parameters. This includes comparing pressure-temperature structures, as described in Sect. 5.3, but also abundance structures with the expected values from chemical equilibrium (e.g. Line et al. 2015, 2017; Gandhi & Madhusudhan 2018; Zalesky et al. 2019). We note that in our case, chemical equilibrium is already imposed within the pRT retrievals.

2. Or by comparing the retrieved cloud location and abundances with the expectations from physics and chemistry (e.g. compare the intersection of the saturation vapour pressure curve and the P–T profile given the expected cloud base, with the retrieved location of the cloud; see Fig. 3 of Burningham et al. 2021). Once more, in our case, this is already included in the pRT retrievals.

3. Or by analysing the robustness of the results when changing the model setup. Regarding the cloud model (i): our model is very strongly physically and chemically motivated, but another option would be a power-law opacity implementation to mimic the clouds (see Cloud Model 2 in Mollière et al. 2020); or various other prescriptions (see e.g. Burningham et al. 2017; Barstow 2020). Regarding the P–T parameterisation (see e.g. the discussion in Mollière et al. 2020) (ii): the P–T parameterisation could also effectively be changed by imposing certain priors during the retrieval; for example by limiting the second derivative of the temperature structure such that solutions that would give a smooth P–T profile are favoured, as proposed by Line et al. (2015), their Sect. 2.4.2. Our current retrievals do not have this implemented, but it could be added in our future work. In general, we note that the robustness of retrieval results when changing the parameterisation does not necessarily indicate correctness.

6 Summary and conclusions

In this work, we present VLT/SPHERE spectro-photometric observations of 51 Eridani b. The new YH spectrum and K1K2 observations show improved S/N compared to previously reported data, allowing us to revise the published flux measurements. We used the radiative transfer code petitRADTRANS, which uses a retrieval approach to fit the atmospheric parameters. In addition, we attempted to reproduce previous results (obtained with self-consistent models) using this approach and compared the outcomes of retrievals to self-consistent models. Our results can be summarised as follows:

• We extracted the spectrum of 51 Eri b using the ANDROMEDA algorithm (Fig. 2). We obtained new photometric measurements for the K1K2 filters (M_K1 = 15.11 ± 0.04 mag, M_k2 = 17.11 ± 0.38 mag; Table 3).

• The detection limits derived from our data show an increased sensitivity and rule out the presence ofplanets more massive than 2 M_Jup at 3 au, and 1 M_Jup beyond 4.5 au (Fig. 4).

• Our initial retrieval runs tended towards clear atmospheres. To verify the robustness of these results, we introduced an additional fit parameter (log τ_clouds) to enforce clouds. We report the results of four different cases in Table 5: a 'nominal' and an 'enforced clouds' version for our new data in conjunction with the photometry in SAM17 and the MS and L_P provided by Rajan et al. (2017); and for the same data used in SAM17.

• We are able to obtain a good fit to the observations with pRT (e.g. Fig. 5), with the exception of the MS photometric point, which can be explained by the large uncertainty of the data. Further mid-IR observations in this band could improve the fit and help constrain the C/O ratio of the planet. We observe that even the 'enforced clouds' retrieval runs tend towards non-cloudiness (log τ_clouds = -0.85 ± 0.16), which differs from previous results obtained using self-consistent models. This discrepancy may be due to the larger and more flexible parameter space that can be explored with retrievals as opposed to self-consistent models. In particular, the isothermal P–T profiles may imitate the effect of clouds.

• Overall, our results (C/O = 0.38 ± 0.09, [Fe/H] = 0.26 ± 0.30 dex, T_eff = 807±45 K and log g = 4.05 ± 0.37) are in agreement with previously reported parameters of the planet. One of the major disagreements is the metallicity, which we find to be close to solar with the new data. Once more, this can be explained by the different methods used by atmospheric retrievals and self-consistent models to fit the data. We estimate the mass of the planet to be between 2 and 4 M_Jup, which is consistent with both 'hot-' and 'warm-start' formation scenarios.

• As an additional test, we used the best-fit parameters from the retrievals to obtain the pressure-temperature structure using a self-consistent model (Fig. 7). The results show a larger temperature gradient for the self-consistent models, which suggests that the T-gradient-cloud correlation plays a role.

Our results highlight the challenges that are still to be overcome when modelling exoplanet atmospheres, as well as the importance of observations at longer wavelengths to determine the presence or absence of clouds. Observations with future instruments that allow the study of additional cloud absorption features, such as ELT/METIS, would be required to allow final conclusions to be made about the cloud-temperature gradient degeneracy.

Acknowledgements

SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zürich (Switzerland), NOVA (Netherlands), ON ERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon (OSUL/CRAL). A.L.M. acknowledges financial support from the Agence Nationale de la Recherche, the European Research Council under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 819155), and the F.R.S.-FNRS.

Appendix A Wavelength calibration and spectral differential imaging

SPHERE/IFS data reduction pipelines remove instrumental signatures, calibrate and compute wavelength solutions, extract spectra of individual lenslets from the 2D detector, and reassemble the data into a 3D data cube with one spectral and two spatial dimensions. The wavelength calibration relies on a range of monochromatic lasers projected on the detector by the calibration unit.

The default of ESO's EsoRex pipeline (Freudling et al. 2013) and the Data Reduction Handling software used by the SPHERE consortium (Pavlov et al. 2008) is to determine the wavelength solution by fitting a second-order polynomial to the spectral calibration data. In addition to using spectral lines for the absolute wavelength calibration, the vlt-sphere Python package (Vigan 2020)⁴ aims at a more refined calibration of the dispersion solution by tracing the radial separation between diagonally opposite satellite spots for each spectral plane in the 3D data cube.

Figure A.1 (top) visualises the respective dispersion solutions inherent to three different IFS pipelines. Data cubes reduced by EsoRex version 0.42.0 include the median of the shortest wavelengths, and the median dispersion of a linear fit to the wavelength solution as keywords in the FITS header. SPHERE DRH and vlt-sphere (version 1.4.3, with a wavelength calibration issue fixed) provide a separate FITS file with wavelengths corresponding to each spectral plane of the 3D data cube. In the case of SPHERE DRH, this is based on the 2nd order polynomial fit.

As the choice of dispersion solution determines the spectral band-width of individual spectral channels, it also influences the recovered spectral energy distribution of the detected astro-physical sources. This has to be considered, e.g. when applying retrieval techniques to the observational spectra. For Spectral Differential Imaging (SDI) data sets, the dispersion solution serves a second purpose by providing the radial λ/D scaling of the speckles.

We notice that the strongest gradients in the vlt-sphere dispersion solution coincide with edges of telluric H₂O absorption bands (Fig. A.1, bottom). While the star itself could be considered as a flat continuum source between neighbouring spectral channels, the edges of telluric absorption bands result in strong gradients in the number of photons recorded as a function of wavelength. As a consequence, at the blue-ward edge of a telluric absorption band a channel records more shorter wavelength ('bluer') than longer wavelength ('redder') photons. The opposite happens at the red-ward edge of an absorption band. The centroids of satellite spots at the blue edge of an absorption band are thus weighted in favour of short wavelength photons, resulting in a smaller radial separation of opposing spots on the detector. The peak of satellite spots at the red edge of an absorption band are slightly further apart. Thus the vlt-sphere 'dispersion solution' is not representative of the intrinsic (smooth) response of the IFS AMICI prism to a 'flat spectrum' source, but representative of the response to a source with the spectral characteristics of the satellite spots.

To correct for the above mentioned effect, one could mask the channels around the water absorption bands and use a cubic relation to fit the position of the satellite spots, which would correct the quadratic dispersion computed by fitting the three (four) diode lasers observed in the wavelength calibration for the YJ (YH) IFS modes.

Fig. A.1 Comparison of different data reduction pipelines for the SPHERE/IFS data. Top: SPHERE/IFS dispersion solutions of the reduced 3D data cubes resulting from the standard ESO pipeline EsoRex, the Data Reduction, and Handling (DRH) of the SPHERE consortium, and the vlt-sphere Python package tracing the separations of satellite spots. Bottom: Satellite spots showing strong spectral gradients in the count rates at the edges of telluric absorption bands.

Appendix B Telluric monitoring and correction

Continuous satellite spots not only provide means for a continuous spatial registration of the star, but also offer a simultaneous monitoring of the (grey) atmospheric transmission, and of variations in the strength of telluric H₂O absorption bands. Figure B.1 visualises the variations in atmospheric conditions over the duration of the full sequence of 154 frames of the 2017-09-28 IFS data set. In the top panel we plot the inverse of the normalised IFS satellite spot count rates observed in the 1.4 μm H₂O absorption band. For better comparison with the Paranal atmospheric monitoring data as made available by the ESO archive⁵, we smoothed the SPHERE/IFS data to the same coarse time sampling of ≈120s. The most noticeable feature is the sharp rise in absorption near the end of the sequence due to an incoming cloud layer. In the middle panel of Fig. B.1 we show the contemporaneous sky infrared temperature as measured by the Low Humidity and Temperature Profiling (LHATPRO) instrument (Querel & Kerber 2014). The rise in the telluric absorption seen in the IFS data approximately coincides with the increase in the sky infrared temperature due to clouds. In the bottom panel of Fig. B.1 we show the precipitable water vapour (PWV) measured by LHATPRO. We notice that there is no exact synchronicity between the IFS measurements in the telluric H₂O absorption band and the LHATPRO PWV. Some shape similarity of the PWV variations between UT ≈ 8.1 hr and ≈ 9.1 hr, and the IFS variations between UT ≈ 8.3 hr and ≈ 9.3 hr, which might be explained by telluric water vapour fluctuations first crossing the LHATPRO field of view, and ≈ 12 min later the SPHERE/IFS field of view, is most likely coincidental. A strict correlation between IFS and LHATPRO telluric measurements is not expected, as they monitor different parts of the sky (SPHERE/IFS was tracking 51 Eri, and LHATPRO was staring at zenith), and also cover different fields of view (40 mas for SPHERE/IFS vs. 1.4 deg for LHATPRO).

Fig. B.1 Variations in the atmospheric conditions during the IFS 2017-09-28 observations. Top: Telluric absorption as traced by SPHERE/IFS at ≈ 1.4μm. Middle: Sky infrared temperature as traced by LHATPRO. Bottom: PWV as traced by LHATPRO. The vertical dashed line in all figures marks our cut-off time for the first 140 IFS frames, which form the basis of our spectral analysis, out of a total of 154 IFS data frames.

The data stress the importance of a simultaneous monitoring of the telluric absorption along the line of sight for high precision (better than ±2% for the first 140 frames of the present data set) spectro-photometric observations of exoplanets. This can be accomplished either by employing high spectral resolution, which facilitates the monitoring of individual lines in telluric H₂O absorption bands, or - in the case of low- to medium spectral resolution (R ≲20 000) observations – by simultaneous monitoring of the spectro-photometric signal of the host star (employing, e.g. continuous satellite spots).

Appendix C Spectral correlation matrix

The extracted exoplanet spectrum from our SPHERE/IFS data is affected by spectral covariance, which can alter the values of the fitted atmospheric parameters. In order to see by how much our data is affected by this, we followed the methods in Greco & Brandt (2016) to estimate the average spectral correlation ψ_ij within an annulus of 1.5λ/D at the separation of the planet, masking out the planet in a 2λ/D radius. Where

(C1)

where C is the covariance matrix, and 〈I_i〉 is the average intensity within the annulus at wavelength λ_i. The covariance matrix is then used to compute the log-likelihood ln ℒ according to

(C.2)

where S is the observed spectrum, and F is the model spectrum. The correlation matrix for our IFS YH spectrum is shown in Fig. C.1.

We ran a 'nominal' retrieval for our 2018 data with 4000 live points using the covariance matrix to compute the log-likelihood. We observed that the values of the fitted parameters remain within error bars to the ones from the retrieval for which we did not use the covariance matrix. However, the best-fit model has higher and lower flux in the Y and J-bands, respectively, compared to the 'nominal' best-fit model without using the covariance matrix. For this reason, we decided to include the covariance matrix in all our retrievals.

Fig. C.1 Spectral correlation matrix between each pair ofspectral channels in our SPHERE/IFS data.

Appendix D Enforced clouds retrieval and attempt to reproduce previous results

Figure D.1 shows the best-fit spectrum along with 34 randomly sampled posterior PDFs for the 'enforced clouds' retrieval using our new data + photometric points. The best-fit parameters are quoted in Table 5, and the posterior PDFs are shown in Fig. D.2, the values quoted for each parameter correspond to the median of the posterior distribution, the uncertainties show the 16th and 84th percentile, representing a 1σ uncertainty range. In this way the values are not identical to the ones that produce the overall best fit, which are given at the top of Fig. D.1. Being the median, they also do not necessarily correspond to the most probable value that can be seen from the peak of the marginalised posterior distribution shown in Fig. D.2.

Figure D.3 shows the best-fit spectrum along with 100 randomly sampled posterior PDFs for the 'nominal' retrieval using SAM17's data. Figure D.4 shows the corresponding posterior PDFs. The parameters for the best fit spectra are summarised in Table D.1. Values and errors quoted in Table D.1 are derived from the equally weighted posterior distribution produced by the Multinest algorithm for each parameter, i.e. marginalised over all parameters except the one in question. The last line gives the parameters derived by Samland et al. (2017) for their best-fitting 'PTC-C' model. That latter model implies a cloud fraction of 100%. Figures D.5 and D.6 are analogous to the above but for the 'enforced clouds' case.

As can be seen from Table D.1 we reproduce most of the parameters to within the calculated uncertainties, albeit with the major difference that our atmosphere shows no significant trace of clouds, and our metallicity is sub-stellar whereas SAM17 found a strongly super-stellar metallicity. Be reminded that τ_cloud denotes the optical depth of the cloud deck at the location where the atmosphere becomes optically thick due to gas opacities alone, i.e. τ_gas ≈ 1. This implies, that our nominal solution shows essentially no clouds at all (τ_cloud ≈ 10⁻⁵@τ_gas ≈ 1), whereas the enforced clouds solution (τ_cloud ≈ 10@τ_gas ≈ 1) requires a rather unphysical gas giant with a radius of 1.56 R_Jup, which for compensation needs to be unusually cool.

Fig. D.1 petitRADTRANS best-fit spectrum of 51 Eri b for the 'enforced clouds' retrieval run on our new SPHERE spectro-photometric data (teal circles and purple squares) along with the photometric points included in Samland et al. (2017) (shown as squares). The photometric points describe the average flux in the respective filter, the x-error bar represents the filter widths. 34 randomly drawn samples from the posterior probability distribution are shown in gray, to show the spread of model parameter combinations to fit the data. Residuals in multiples of 1σ uncertainties of the data are shown below.

Fig. D.2 Corner plot of the posterior PDFs of the 'enforced clouds' retrieval run on the new data set.

Fig. D.3 Best-fit spectrum (top) of the nominal retrieval run on the original data set from Samland et al. (2017). This is to be compared to Fig. 11 in Samland et al. (2017).

Fig. D.4 Corner plot of the posterior PDFs of the retrieval run on the original data set from Samland et al. (2017). This is to be compared to Fig. 12 in Samland et al. (2017).

Fig. D.5 Best-fit spectrum of the posterior PDFs of the retrieval run on the original data set from Samland et al. (2017) when restricting the range of the τcloud prior to positive values.

Fig. D.6 Corner plot of the posterior PDFs of the retrieval run on the original data set from Samland et al. (2017) when restricting the range of the τcloud prior to positive values.

References

All Tables

All Figures

	Fig. 1 Median combined S/N detection map from ANDROMEDA based on SPHERE/IFS data. The circle indicates the position of 51 Eri b. The azimuthal wings around the planet signal is the characteristic planet signature that ANDROMEDA fits for ADI data.
In the text

	Fig. 2 Newly obtained 51 Eri b YH spectrum and K1K2 photometry. Overplotted for comparison are the YH spectrum and K1K2 photometry from Samland et al. (2017), as well as the updated GPI J and H spectra and the K1 and K2 spectra from Rajan et al. (2017).
In the text

	Fig. 3 Flux uncertainty over wavelength as a function of wavelength for the YH spectra of 51 Eri b from this work compared to the one presented in Samland et al. (2017).
In the text

	Fig. 4 5σ detection limits. Planet contrast (top) with respect to the star as a function of separation to the star, and planet mass detection limits (bottom) as a function of the separation from the star for the SPHERE/IRDIS_EXT data of 51 Eri b. The mass limit for the IFS data is cut to the lowest mass computed by the model grid (1 MJup) for both ANDROMEDA and TRAP.
In the text

Fig. 5 petitRADTRANS best-fit spectrum of 51 Eri b for the 'nominal' retrieval run (shown in black) on our new SPHERE spectro-photometric data (teal circles and purple squares) along with the photometric points included in SAM17 and (Rajan et al. 2017; shown as squares). The photometric points describe the average flux in the respective filter, and the x-error bar represents the filter widths. The 34 randomly drawn samples from the posterior probability distribution are shown in grey in order to show the spread of model parameter combinations to fit the data. Residuals in multiples of 1 σ uncertainties of the data for the best-fit model are shown below.

In the text

	Fig. 6 Corner plot of the posterior PDFs of the 'nominal' retrieval run on the new data set.
In the text

	Fig. 7 Retrieved pressure-temperature profiles in white, with confidence intervals (1σ, 2σ, and 3σ) for our four different retrievals shown in Table 5. Overplotted are the resulting self-consistent P–T structures in pink.
In the text

Fig. A.1 Comparison of different data reduction pipelines for the SPHERE/IFS data. Top: SPHERE/IFS dispersion solutions of the reduced 3D data cubes resulting from the standard ESO pipeline EsoRex, the Data Reduction, and Handling (DRH) of the SPHERE consortium, and the vlt-sphere Python package tracing the separations of satellite spots. Bottom: Satellite spots showing strong spectral gradients in the count rates at the edges of telluric absorption bands.

In the text

Fig. B.1 Variations in the atmospheric conditions during the IFS 2017-09-28 observations. Top: Telluric absorption as traced by SPHERE/IFS at ≈ 1.4μm. Middle: Sky infrared temperature as traced by LHATPRO. Bottom: PWV as traced by LHATPRO. The vertical dashed line in all figures marks our cut-off time for the first 140 IFS frames, which form the basis of our spectral analysis, out of a total of 154 IFS data frames.

In the text

	Fig. C.1 Spectral correlation matrix between each pair ofspectral channels in our SPHERE/IFS data.
In the text

Fig. D.1 petitRADTRANS best-fit spectrum of 51 Eri b for the 'enforced clouds' retrieval run on our new SPHERE spectro-photometric data (teal circles and purple squares) along with the photometric points included in Samland et al. (2017) (shown as squares). The photometric points describe the average flux in the respective filter, the x-error bar represents the filter widths. 34 randomly drawn samples from the posterior probability distribution are shown in gray, to show the spread of model parameter combinations to fit the data. Residuals in multiples of 1σ uncertainties of the data are shown below.

In the text

	Fig. D.2 Corner plot of the posterior PDFs of the 'enforced clouds' retrieval run on the new data set.
In the text

	Fig. D.3 Best-fit spectrum (top) of the nominal retrieval run on the original data set from Samland et al. (2017). This is to be compared to Fig. 11 in Samland et al. (2017).
In the text

	Fig. D.4 Corner plot of the posterior PDFs of the retrieval run on the original data set from Samland et al. (2017). This is to be compared to Fig. 12 in Samland et al. (2017).
In the text

	Fig. D.5 Best-fit spectrum of the posterior PDFs of the retrieval run on the original data set from Samland et al. (2017) when restricting the range of the τcloud prior to positive values.
In the text

	Fig. D.6 Corner plot of the posterior PDFs of the retrieval run on the original data set from Samland et al. (2017) when restricting the range of the τcloud prior to positive values.
In the text