© C. Desgrange et al. 2022

1 Introduction

Direct imaging has proven to be successful at imaging and characterizing the properties of young planetary system architectures, and young (≤100 Myr) gaseous giant planets orbiting at large distances from their host stars (a ≳ 10 au). Large-scale surveys of several hundreds of young nearby stars (<150 pc), such as SHINE (Chauvin et al. 2017; Desidera et al. 2021; Langlois et al. 2021; Vigan et al. 2021) and GPIES (Nielsen et al. 2019), have now extended our vision of the giant planet demographics down to 10 au. The goal is to access the bulk of the giant planet population close to the snowline and to bridge the gap with complementary indirect detection methods such as radial velocity, transit, micro-lensing, and soon astrometry (with the Gaia Data Release 4) that are sensitive to planets closer to their stars (≲ 10 au). Over the last two decades discoveries of emblematic planetary systems such as HR 8799 (Marois et al. 2008, 2010), β Pictoris (Lagrange et al. 2010), HD 95086 (Rameau et al. 2013a), 51 Eri (Macintosh et al. 2015), and PDS 70 (Keppler et al. 2018; Müller et al. 2018) have offered rich opportunities to explore the diversity of young Solar System analogs, containing giant planets and circumstellar disks shaped with cavities and belts. HR 8799 hosts four exoplanets (Marois et al. 2008, 2010), whereas 51 Eridani hosts one exoplanet (Macintosh et al. 2015). A second planet has been recently detected by the radial velocity technique in the β Pictoris system (Lagrange et al. 2019) and confirmed with interferometric observations (Nowak et al. 2020), while a second forming planet has been detected with MUSE in the PDS 70 system (Haffert et al. 2019). The mass of these protoplanets is highly uncertain, with estimates ranging from 1 to 17 M_Jup (Müller et al. 2018; Christiaens et al. 2019; Mesa et al. 2019; Isella et al. 2019; Stolker et al. 2020b; Wang et al. 2021). These young, planetary systems are benchmark laboratories for exploring the formation and evolution of young giant planets with the current large telescopes and instruments. They are prime targets for upcoming telescopes, such as the James Webb Space Telescope (JWST, first light 2022) and the Extremely Large Telescope (ELT, first light 2027). From that perspective, the hunt for additional planets in the young system HD 95086 is very interesting.

Since 2013, and since the discovery of a 4–5 M_Jup exoplanet HD 95086 b in thermal imaging using the NaCo instrument (Rousset et al. 2003; Lenzen et al. 2003) at the Very Large Telescope (VLT), the HD 95086 planetary system has become a reference to investigate the processes of planetary formation and evolution, and to characterize young planetary architectures. The star is an A-type star with an approximate effective temperature of 7750 K and a mass of 1.6M_⊙. Until very recently, the star was identified as a young star located at the border of the Lower Centaurus Crux (LCC) association and thus with an age of 17 ± 2 Myr (Pecaut et al. 2012). However, based on the Gaia Data Release 2, Booth et al. (2021) showed that the star could instead belong to the Carina association, which according to their age estimation would make it a few million years younger, ${13.3}_{-\mathrm{0..6}}^{+1.1}$ Myr. By using a self-consistent Bayesian analysis, Swastik et al. (2021) derived anew estimation of the stellar parameters consistent with the literature but more accurate, and in particular a stellar metallicity of ${0.14}_{-0.04}^{0.05}$. The physical parameters of the star are summarized in Table 1.

HD 95086 hosts a double-belt debris disk architecture, very similar to that of our Solar System, with an outer belt resolved by the Atacama Large Millimeter Array (ALMA) in the continuum at 1.3 mm (Su et al. 2017; see Fig. 1). The inner warm belt is located at 8 ± 2 au (187 ± 26 K), and the large outer, colder belt between 106 ± 6 au and 320 ± 20 au (57 ± 2 K). Their existence was originally identified from the analysis of Herschel observations, in combination with the characterization of the spectral energy distribution (SED) of HD 95086 (Moór et al. 2013). Based on SED modeling from Herschel, Spitzer, WISE, and APEX observations, the existence of a third belt at 2 au (300 K), has been also proposed by Su et al. (2015), together with a disk halo component that could extend up to 800 au (Su et al. 2017), but this innermost belt has not been confirmed to date. Recently, Zapata et al. (2018) added new constraints on the structure of the outer belt at submillimeter and millimeter wavelengths with ALMA observations at 0 9 and 1 3 mm, and derived a dust-to-gas ratio ≥50. Moreover, Zapata et al. (2018) and Booth et al. (2019) did not detect CO (J = 2–1) and (J = 3–2) emissions, excluding the possibility of HD 95086 being an evolved gaseous primordial disk. By using spectro-spatial filter on ALMA observations, Booth et al. (2019) found tentative evidence of CO (J = 2–1) emission with an integrated line flux of 9.5 ± 3.6 mJy kms⁻¹. It corresponds to a CO mass of (1.4–13) × 10⁻⁶ M_⊕, which they determined to be consistent with second-generation production of gas through collisional cascade (Kral et al. 2017). According to Su et al. (2015), the collisions in the HD 95086 disk might also explain their detection of a 69 µm crystalline olivine feature from the outer disk with the Spitzer telescope as the crystallization of olivine requires a high temperature, as is the case for instance in the core of planetary bodies after their disruption. Finally, the outer belt has also been marginally detected in polarized scattered light in the near-infrared (J band) by SPHERE differential polarimetric imaging (DPI) observations (Chauvin et al. 2018) colocated with the thermal emission seen by ALMA. The physical parameters of the debris disk architecture are given in Table 1.

At infrared wavelengths, following the discovery of HD 95086 b with VLT/NaCo (Rameau et al. 2013a) in the L′ band (3.8 µm), the planet was re-imaged using Gemini/GPI in the H band (1.5–1.8 µm) and in the K1 band (1.9–2.2 µm) (Rameau et al. 2016; De Rosa et al. 2016), and using VLT/SPHERE with IRDIS in H2H3 filters (λ_H2 = 1.593 µm, λ_H3 = 1.667 µm) and in K1K2 filters (λ_K1 = 2.103 µm, λ_K2 = 2.255 µm), and with the integral field spectrograph (IFS) in the YJ (0.95-1.35 µm), and YJH (0.97–1.66 µm) settings (Chauvin et al. 2018). The combination of different photometric measurements in the infrared enabled De Rosa et al. (2016) and Chauvin et al. (2018) to confirm the late L spectral type of HD 95086 b, which is consistent with a dusty atmosphere of about 800-1300 K.

The first orbital fitting of HD 95086 b was performed by Rameau et al. (2016) from previous NaCo astrometric data (epochs 2012 to 2013 from Rameau et al. 2013a) and GPI astrometric monitoring between 2013 and 2016 (published partially in Galicher et al. 2014). They found a semimajor axis of ${62}_{-8}^{+21}$ au, an eccentricity less than 0.21, and an inclination of ${153}_{-14}^{+10\circ }$ at the 68% confidence interval by using Monte Carlo methods. Chauvin et al. (2018) updated the orbital solution using a larger orbital coverage, this time combining NaCo and SPHERE astrometric measurements. This recent MCMC analysis gave consistent results and showed that the planet is orbiting with a period of about ${289}_{-177}^{+12}$ years, a semimajor axis of ${52}_{-24}^{+13}$ au, a relatively low eccentricity (${0.2}_{-0.2}^{+0.3}$), and with an inclination of ${141}_{-13}^{+15\circ }$ at the 68% confidence interval, compatible with a coplanar orbit with the debris disk plane. The physical parameters of the exoplanet are summarized in Table 1. In addition, Chauvin et al. (2018) used the High Accuracy Radial velocity Planet Searcher (HARPS) high-resolution optical spectrograph to search for additional exoplanets with the radial velocity (RV) technique, and could exclude the presence of a very massive (>10 M_Jup), coplanar inner giant planets at less than 1 au.

Given the large cavity seen in ALMA images inside the cold outer belt (see Fig. 1), the system HD 95086 very likely hosts at least one or perhaps two additional planets closer to the star than HD 95086 b, which would explain the architecture of the two debris belts. Su et al. (2015), Rameau et al. (2016), and Chauvin et al. (2018) investigated this possibility by considering various locations, eccentricities, and masses for the inner planets, the physical properties of b, and the characteristics of the inner and outer belts together with the detection performance of current planet imagers. A configuration with one or two additional inner planets between 10 and 30 au, dynamically stable with b carving the outer belt, and participating in the replenishment of the inner belt is possible and worth investigating.

In this paper we extend the study of Chauvin et al. (2018) to revisit the orbital and atmospheric properties of HD 95086 b, and the presence of additional inner giant planets, considering a total of ten epochs acquired with the VLT/SPHERE instrument (Beuzit et al. 2019) between February 2015 and May 2019. These datasets include five new unpublished epochs covering January 2018 to May 2019. In Sect. 2 we present the data acquired, together with the archival data used for this analysis. In Sect. 3 we describe the image processing methods used, along with our data selection and the decision to combine the different datasets considering the individual epoch contrast performance and adaptive optics (AO) correction quality. In Sect. 4 we present the updated astrometry for the exoplanet HD 95086 b based on VLT/NaCo and VLT/SPHERE data, covering a total of seven years of monitoring between 2012 and 2019, and determine the best orbital solution. In Sect. 5 we present for the first time the SPHERE-JH (1 2–1 6 µm) spectroscopic observations of HD 95086 b. Using the MCMC special code (Christiaens et al. 2021) applied to the combined spectrum, we re-analyze the physical parameters of the planet and investigate the presence of a circumplanetary disk around it. In Sect. 6 we finally look for the hypothetical exoplanets c and d in the system by updating the HARPS and SPHERE combined observations, and also by applying the K-Stacker algorithm (Le Coroller et al. 2015) to the SPHERE multi-epoch datasets.

Fig. 1 Composite ALMA-continuum (at 1.3 mm) and SPHERE/IRDIS (at 2.1 µm) observations of HD 95086 (see Su et al. 2017; Chauvin et al. 2018). The exoplanet b is detected at the K1 band (red dot). The white dashed ring at 180 au represents the peak location of the outer cold belt located from 106 ± 6 au to 320 ± 20 au. The inner warm belt is not resolved with ALMA. The pink ellipse represents the ALMA synthetic beam.

2 Observations

2.1 VLT/SPHERE data

The HD 95086 system was monitored during the SHINE survey at 13 different epochs between February 2015 and May 2019 (see Table 2) using the VLT/SPHERE high-contrast instrument (Beuzit et al. 2019). The observations were obtained with the modes IRDIFS (3 epochs) and IRDIFS-EXT (11 epochs) that combine simultaneously the IRDIS (Dohlen et al. 2008) and IFS instruments (Claudi et al. 2008). The IRDIFS-EXT mode combines IRDIS in dual-band imaging (DBI; Vigan et al. 2010) mode with the K1K2 filter doublet λ_K1 = 2.103 ± 0.102 µm, λ_K2 = 2.255 ± 0.109 µm, and IFS in the YJH (0.97–1.66 µm) setting. The IRDIFS mode combines IRDIS in DBI with H2H3 filters (λ_H2 = 1.593 ± 0.055 µm, λ_H3 = 1.667 ± 0.056 µm), and IFS in the YJ (0.95-1.35 µm) setting. Each observing sequence was performed with the pupil-tracking mode. This combination enables the use of angular (Marois et al. 2006) and/or spectral differential imaging techniques (Racine et al. 1999; Sparks et al. 2002) to reach higher contrast at subarcsecond separations. The details of the observations are reported in Table 2.

In this work we focused our analysis on data acquired with the IRDIFS-EXT mode as HD 95086 b, and any expected inner planet in the system, are L-type planets, expected to be particularly red and therefore easier to detect at longer wavelengths. A total of ten epochs are considered as the data acquired on May 3, 2017, are not exploitable owing to very poor observational conditions. The observing conditions are summarized in Table 2 and Fig. A.1. The Strehl ratio (SR) and the wind parameters are measured by the SPHERE eXtreme AO (SAXO, Petit et al. 2014) real-time computer named Standard Platform for Adaptive optics Real Time Applications (SPARTA, Fedrigo et al. 2006), while the seeing (є) and the atmospheric coherence time parameters (τ₀) were obtained by the Differential Image Motion Monitor (DIMM) and the Multi-Aperture Scintillation Sensor (MASS, Kornilov et al. 2007) turbulence monitor at the Paranal Observatory.

2.2 Archival data: VLT/NaCo and Gemini-South/GPI

To revisit the orbital and spectral properties of HD 95086 b, we analyzed archival data from the VLT/NaCo imager obtained in 2012 and 2013 (Rameau et al. 2013a,b), together with Gemini-S/GPI obtained between 2013 and 2016 (Rameau et al. 2016; De Rosa et al. 2016). A summary of these observations and the astrometric and spectro-photometric results used for this work are reported in Table 3.

3 Data reduction and analysis

3.1 Pre-processing

All SPHERE observations of HD 95086 were reduced by the SPHERE Data Center (Delorme et al. 2017a), using the SPHERE Data Reduction and Handling pipeline (Pavlov et al. 2008), following the same approach as described by Chauvin et al. (2018). To summarize, the pre-processing corrects the non-coronagraphic point spread function (PSF) and the coronagraphic image cube for bad pixels, dark current, flat non-uniformity; the sky background for both IRDIS and IFS; and the wavelength and cross-talk between spectral channel calibration for IFS. A normalization is applied to calibrate the coronagraphic images in intensity relative to the star (i.e., in terms of contrast). The coronagraphic images are centered by using the four satellite spots to accurately determine the position of the star behind the coronagraphic mask, as also described in Chauvin et al. (2018). To calibrate the astrometry of both IRDIS and IFS on the sky, a star-crowded field (47 Tuc) is regularly observed as part of the long-term analysis of the SPHERE guaranteed time observation (GTO) astrometric calibration described in Maire et al. (2016, Maire et al. 2021a) to measure the detector plate scale, true north, and distortion. The plate scale and true north solutions at each epoch are reported in Table 3.

3.2 Image processing

To detect and characterize potential planetary signals in the images, we used two dedicated pipelines, namely ANDROMEDA (Cantalloube et al. 2015) and SpeCal (Galicher et al. 2018). Both are based on the angular differential imaging (ADI) technique, which removes the starlight residuals in the coronagraphic images. ANDROMEDA is a forward-modeling approach based on a maximum likelihood estimator (Mugnier et al. 2009). It first performs a simple pair-wise subtraction. Then, it searches for the specific signature that would appear in the presence of an unresolved point-source in the residual image, and estimates its probability, jointly for all pairs of subtracted images. The SpeCal pipeline combines a set of different algorithms like classical ADI (cADI, Marois et al. 2006), locally optimized combination of images (LOCI/TLOCI, Lafrenière et al. 2007), and principal component analysis (PCA, Amara & Quanz 2012; Soummer et al. 2012). To exploit the spectral diversity given by the IFS and the IRDIS-DBI modes, in addition to the temporal dimension as done in ADI, SpeCal has been developed to apply LOCI/TLOCI and PCA in angular and spectral differential imaging (ASDI; Mesa et al. 2015). For the analysis of the SHINE survey, the reference algorithms benchmarked with various blind tests are the TLOCI-ADI and the PCA-ASDI algorithms (Langlois et al. 2021). For IRDIS, the contrast performance showed that ANDROMEDA-ADI performs better than TLOCI-ADI (see Figs. 2 and 3), and in addition provides a more robust estimate of the statistical threshold for the candidate detection. The signal-to-noise ratio maps (S/N maps) obtained with TLOCI-ADI for IRDIS and PCA-ASDI for IFS using SpeCal are shown in Fig. B.1. The S/N maps obtained for ANDROMEDA in ADI for IRDIS and IFS (consisting of a combination of the IFS channels after processing individual each channel) are shown in Fig. 2.

3.3 Multi-epoch selection and combination

Given the very red spectrum of HD 95086 b, the planet is clearly detected for each individual epoch (with S /N > 10) in the K1 and K2 bands with SPHERE, and in the K band with GPI. This band (K1, where the background noise is lower than in K2) is therefore used to extract the planet’s relative astrometry at each epoch, and derive its updated orbital properties (see below). The detection becomes more challenging at the H band (S/N ~ 3–7) for both SPHERE and GPI, but remains possible for good-quality observation (January 2018, February 2018, March 2018, April 2019 and May 2019). At the J band the planet is currently detected only by stacking reduced SPHERE-IFS images taken at various epochs to optimize the speckle cancelation, as done by Chauvin et al. (2018). Following a similar multi-epoch strategy, we first selected the best observing epochs obtained by SPHERE between 2015 and 2019, then re-aligned each final IFS datacubes correcting for the planet’s orbital motion, and stacking them to extract the JH spectrum of the planet.

To investigate the contrast performance at each epoch, we considered as the first criterion the contrast curves determined with TLOCI-ADI and ANDROMEDA-ADI for IRDIS, and PCA-ASDI and ANDROMEDA-ADI for IFS (Fig. 3). The observing conditions (atmospheric turbulence conditions and AO telemetry, summarized in Table 2 and in Fig. A.1) were used as a consistency check. From a total of ten epochs, a subsample clearly emerges of six very good epochs. The best epoch is January 2018, followed by February 2018, March 2018, April 2019, May 2019, and May 2017. The worst epochs are January 2016 and February 2015. Despite the very good Strehl ratio (81%) and seeing (0.33 as) for the January 2016 epoch, the atmospheric coherence time is very fast (1.8 ms), producing a wind-driven halo, which degrades the contrast performance by about one order of magnitude at the planet location (Cantalloube et al. 2020).

The 5σ detection limits are 1 × 10⁻⁶ and 2 × 10⁻⁶ at 0.6″ respectively with the IFS and IRDIS for the best epoch (January 2018). IFS performs better in contrast at close separations as PCA-ASDI exploits the spectral diversity of 39 spectral channels to remove the speckles, while TLOCI-ADI does not for the IRDIS dual-band K12. We note that the IFS and IRDIS contrast curves are coherent with each other: the best contrast curves with IFS are the same as with IRDIS (see Fig. 3). In summary, the best epochs for extracting the spectrum of the exo-planet HD 95086 b are the last six: May 2017, January 2018, February 2018, March 2018, April 2019, and May 2019. To extract the astrometry, we preferred the epochs acquired with the satellite spots for the whole sequence of observation (and not only at the beginning and/or the end of the observation) to ensure a correct astrometric calibration for the whole sequence of observation (May 2016, May 2017, January 2018, and April 2019), as well as the two epochs in 2015 (February 2015 and May 2015), as in 2015 no observation with satellite spots for the whole sequence was acquired. We note that the position of the star on the detector can marginally evolve during the sequence of observation due to optomechanical and thermal variations that slightly move the optics. The Differential Tip Tilt Sensor (DTTS, Baudoz et al. 2010) aims to ensure its centering, but a centroid variation can be expected. We adopted the conservative error on the centroid position of 0.2 pixel corresponding to 2.5 mas for the epochs acquired without the continuous satellite spots mode, as done in Chauvin et al. (2018). This value is taken into account in the astrometric calibration error budget, as well as the errors on the true north, plate scale, and pupil offset, following Maire et al. (2021b). In addition, it should be noted that the epoch of March 2018 suffers from a higher centroid variation, which is a rare problem, but even so occurred in a few observations over the whole seven years of exploitation of the SPHERE instrument. This higher centroid variation is about 5 mas both in declination and right ascension based on the expected position of HD 95086 b in March 2018 considering the epochs better calibrated in astrometry, as of January 2018 and April 2019.

The relative astrometry and photometry of the exoplanet HD 95086 b to its host star are gathered in Table 3 for all the epochs obtained with SPHERE-IRDIS in the Kl band, using ANDROMEDA and SpeCal-TLOCI image processing, along with the measurements from archival NaCo and GPI observations. The ANDROMEDA and SpeCal-TLOCI astrometric and photometric measurements from 2015 to 2019 are consistent with each other (see Table 3).

Fig. 2 Signal-to-noise ratio maps for all the SPHERE epochs. At the top, SPHERE-IRDIS in the K1 band and at the bottom, SPHERE-IFS in the YJH bands, both reduced by the pipeline ANDROMEDA. The color bar corresponds to the signal-to-noise ratio, 0–15 for IRDIS and 0–5 for the IFS. Planet b is located in the white square. The region below the inner working angle of the coronagraph is masked. In the top right corner of the IFS data, the one (or two) point-like feature(s) correspond to the remanent of the star on the IFS detector.

Fig. 3 Detection limits at 5σ for IRDIS in the K1 band on the left and IFS in the YJH bands on the right. Data are reduced at the top by the pipeline SpeCal-TLOCI (for IRDIS) or SpeCal-PCAPad (for IFS) and at the bottom by the pipeline ANDROMEDA. The gray area corresponds to the area hidden by the larger coronagraph used for all the epochs (i.e., N_ALC_Ks) which has an inner working angle of 0.107″ and 0.116″ for IFS and IRDIS, respectively. The blue vertical line corresponds to the position of the planet HD 95086 b at the best epoch January 2018 (cyan, 0.625″).

3.4 The 1.2-3.8µm spectrum of HD 95086 b

For the first time, we have extracted the spectrum of the planet HD 95086 b in the J and H bands. This was achieved on the SpeCal-PCA ADI reduced images by measuring the contrast of the planet in an aperture of 1 FWHM of the six-epoch averaged image in which we corrected for the exoplanet’s orbital motion by using its precise astrometry through time before stacking. The exoplanet is still not detectable in the Y band, hardly in the J band, and with a signal-to-noise ratio of about 5 per spectral channel in the H band. By using the BT-NEXTGEN synthetic spectrum of the star HD 95086 A, we converted the contrast values to flux values. We completed the SPHERE-IFS YJH (0.96–1.64 µm, resolution 30) and SPHERE-IRDIS K12 photometric points (2.10 and 2.26 µm) with archival data from GPI, which provides the spectrum of the exoplanet HD 95086 b in the K1 band at low resolution (1.95–2.20 µm, resolution 66), and from NaCo, which provides a photometric measurement in the L′ band (3.80 µm, bandwidth 0.62 µm). The complete spectrum of HD 95086 b is shown in Fig. 4; a zoomed-in image of the J and H bands is shown in Fig. 5. The increasing slope in the spectroscopic K1 and photometric K12 and L′ points acknowledges the redness of exoplanet b (De Rosa et al. 2016) and seems to be verified in our H-band spectrum as well, even though the uncertainties are significant. We empirically estimated the spectral correlation matrix as in Greco & Brandt (2016) and De Rosa et al. (2016) for the measurements obtained with the IFS of SPHERE and GPI (see Appendix C), and used it to compute the covariance matrix used for spectral fitting of HD 95086 b (see Sect. 5).

4 Orbital analysis

We ran the Markov chain Monte Carlo (MCMC) orbital fit (Ford 2005, Ford 2006), as done in Chauvin et al. (2012). We used previous NaCo astrometric data (epochs 2012 to 2013 from Rameau et al. 2013a) and SPHERE measurements (epochs February 2015, May 2015, May 2016, and May 2017 from Chauvin et al. 2018) including two new astrometric points from the SPHERE-IRDIS images obtained in the K1 band (January 2018 and April 2019) with the updated distance value of 86.2 pc from Bailer-Jones et al. (2018). In this run the stellar mass was fixed to its mean value 1.6 M_⊙ and to the stellar distance (86.2 pc). It would have been technically possible to leave the stellar mass free and to let the code redetermine it, but this turned out to be inaccurate. The reason is that only a tiny part of the orbital period is covered by the observations, leading to a degeneracy between the central mass and the inclination.

The priors assumed for this run were logarithmic between 1 and 4000 yr; linear for the eccentricity e between 0 and 1; ∝ sin i for the inclination i between 0 and 180°; and linear between −180° and 180° for the longitude of ascending node Ω, the argument of periastron ω, and for the mean anomaly at the time of the first observation epoch (related to the time of periastron passage T_p). We note that the MCMC run is represented by taking the orbital period P as variable instead of the semimajor axis a. Both approaches are equivalent as P and a are linked via Kepler’s third law.

In addition, when dealing with pure relative astrometric data as we do here, it is well known that there is a ±180° degeneracy between solutions in the longitude of ascending node Ω and the argument of periastron ω. Each solution with (Ω, ω) yields exactly the same projected orbit as the same solution, but with (Ω + 180°, ω + 180°). To overcome this difficulty, as explained in Chauvin et al. (2012), the code actually fits Ω + ω and ω − Ω rather than Ω and ω directly. The former angles are indeed unambiguously determined contrary to ω and Ω. Then, each root solution with fitted values for Ω + ω and ω − Ωis declined as two final separate solutions, one with (Ω, ω) and the other with (Ω + 180°, ω + 180°).

The updated orbital parameters and confident regions from the MCMC orbital fit with the two new astrometric measurements in 2018 and 2019 are shown in Table 4. A sample of the orbit solutions is displayed in Fig. 6, as is the orbit given in Table 4, which corresponds to the maximum a priori probability (MAP). By using both the prior and the reduced chi-squared (${\chi }_r^2$) information, the MAP maximizes the probability given in Eq. (1) used in the MCMC: $$\text{probability}=\frac{\sin \left(i\right)}{P}\times \exp \left(-{\chi }_{\text{r}}^2/2\right).$$(1)

As the MAP gives the peak of probability in the 6D orbital parameter space, the MAP may not correspond to the peak of each 1D distribution (see black vertical lines in diagonal panels in Fig. 7), due to the correlation between the orbital fitted elements. Nevertheless, the 2D distributions in Fig. 7 indicate that the MAP solution (black star) corresponds better to the 2D peaks, which is closer to the 6D reality.

The results are consistent with the previous analysis led by Chauvin et al. (2018). This is expected, first because the fit corresponds to a linear part of the orbit, and second because it only covers a small percent of its whole orbit. In Fig. 7 we can see that a subsample of very eccentric solutions are found (e ≥ 0.4), but they are correlated with lower inclinations (i ≤ 140°) than the inclination of the outer belt (147-153°). Hence, these orbital solutions could still be consistent with the double-belt architecture of the system, even though they represent a small subsample of the MCMC orbital fit solutions.

All in all, these solutions confirm that the exoplanet HD 95086 b, located at a semimajor axis 51–73 au and with a low eccentricity (e ≤ 0.18), is likely sculpting the inner edge of the outer ring, and cannot alone sustain the large cavity observed between 10 and 106 au (Su et al. 2015; Rameau et al. 2016).

By using only the four epochs imaged with SPHERE-IRDIS in the K1 band and acquired with the satellite spots enabling a better astrometric calibration (see Table 2), we also obtained the orbital parameters of HD 95086 b as a by-product of the K-Stacker algorithm (Le Coroller et al. 2015), which is used below in the search for one or two additional inner planets in the system (see Sect. 5). K-Stacker is an optimization algorithm that takes advantage of the Keplerian motion of exoplanets on several epochs to then recenter the images according to their Keplerian motion. Thus, it differs from other classical orbital fitting methods, such as MCMC, as it takes as input images from several epochs instead of derived astrometric positions, and consequently can detect objects otherwise unreachable in each observation considered individually. Hence, the coherence of the orbits found by K-Stacker (within 2σ) illustrates its ability to constrain objects with small orbital motion between epochs (50 mas, i.e., ~ 1 FWHM in the K1 band) and very faint objects as well; about the same orbit has been found in the H band, where HD 95086 b is hardly detectable (i.e., S/N ≤ 5 for each epoch taken individually, see Fig. 2). Within the error bars from orbital parameters similar results to the MCMC solutions are found, although the MCMC fit uses the NaCo and SPHERE data, which covers a longer timescale 2012–2019 (instead of 2016–2019 for K-Stacker). In addition, K-Stacker allows us to constrain the star mass between 1.56 and 1.59 M_⊙ (as described in Le Coroller et al. 2020, Sect. 3.4).

Fig. 4 Spectroscopic and photometric values of the exoplanet HD 95086 b expressed in contrast relative to the host star from 1.2 µm to 3.8 µm with SPHERE-IFS (blue triangles), SPHERE-IRDIS (yellow squares), GPI (pink diamonds), and NaCo (red circle) in the bands JH, K, K12, and L′ at the top. At the bottom, the transmission curves of the filters K1, K2, and L′ are shown1. Results are extracted from an averaged stacking of the six best epochs for the SPHERE data, and for one single epoch for the GPI and NaCo data.

Fig. 5 Spectrum of the exoplanet HD 95086 b expressed in contrast relative to the host star from SPHERE-IFS data in the J band (1.21–1.32 µm) and in the H band (1.46–1.62 µm) at a 68% confidence level. The contrast is given in units of 10−6. The data are reduced by the pipeline SpeCal-PCA ADI where a six-epoch averaged stack was made before extracting the spectrum to boost the S/N.

Fig. 6 Astrometric positions of the planet HD 95086 b between 2012 and 2020 with three different instruments: SPHERE (yellow squares), GPI (pink diamonds), and NaCo (red circles) with 1σ error bars. The astrometric positions from SPHERE were computed with the SpeCal-TLOCI pipeline. The letter “w” in the legend indicates the observations imaged with satellite spots that enable finding the exact position of the star during the whole sequence and later recenter the frames if necessary. The green and orange solid lines correspond to a sample of the orbital solutions found by the MCMC orbital fit and the K-Stacker tools, respectively. The black dashed and dotted lines respectively represent the MCMC and K-Stacker orbits for which the corresponding orbital parameters are given in Table 4. For the MCMC tool, the orbit corresponds to the MAP (see text), whereas for K-Stacker it is the orbit closest to the mean of the orbital solutions. Both of these orbits are by construction true orbital solutions. Notes: The Feb2015 point from SPHERE is hidden by the two 2015 points from GPI; the two 2019 points are very close to each other. The insert (in the top right corner) shows the location of the planet relative to the cold outer belt and warm inner belt from Su et al. (2017).

5 Spectral characterization

The first studies of the infrared colors of HD 95086 b rapidly showed that the spectral properties of the planet fall at the late L to L/T transition and that the planet is underluminous compared to the field dwarfs of similar spectral types (Galicher et al. 2014; De Rosa et al. 2016; Chauvin et al. 2018). The red colors and underluminosity, as shown in Fig. 8, are characteristic of young L/T objects, and are often associated with the inhibited settling of dust in the upper parts of low surface gravity atmospheres. Based on photometric H (1.5–1.8 µm) and spectroscopic K1 (1.9–2.2 µm) observations of HD 95086 b with GPI, De Rosa et al. (2016) confirmed the L-type dusty atmosphere, as evidenced by a featureless low-resolution spectrum and a monotonically increasing pseudo-continuum in the K1 band consistent with a cloudy atmosphere. Considering the 1.2–1.6 µm spectrum extracted in this work, we propose below to reinvestigate the spectral properties of HD 95086 b. We will consider the best atmosphere models fitting current observations from 1.2 to 3.8 µm, as well as the possibility of having a circumplanetary disk around the planet b, as seen for instance in the younger Solar System analog PDS 70, also a member of Sco-Cen for the planet c (Keppler et al. 2018; Isella et al. 2019).

5.1 Comparison to models

To constrain the physical properties of the exoplanet HD 95086 b with atmospheric and circumplanetary disk models, we used the special package (Christiaens et al. 2021), first known as specfit, a module of the open-source python package VIP² (Gomez Gonzalez et al. 2017) and now available as a distinct package³. The special package is compatible with any atmospheric grid if a snippet function which reads the input grid files is provided. It can also fit for blackbody components (either alone or as additional component(s) to the atmosphere), the optical extinction A_v, the optical-to-selective extinction ratio R_v, or the intensity of emission lines that are provided in a dictionary, special utilizes the MCMC sampler emcee (Foreman-Mackey et al. 2013) to retrieve in a Bayesian framework the most likely physical parameters of any stellar or substellar object based on its spectral energy distribution. Models are linearly interpolated between grid points. The log-likelihood expression provided to the sampler is $$\log ℒ\left(D/M\right)=-\frac{1}{2}\left[W{\left(F_{\text{obs}}-F_{\mathrm{mod}}\right)}^T\right]C^{-1}\left[W^T\left(F_{\text{obs}}-F_{\mathrm{mod}}\right)\right],$$(2)

where F_obs and F_mod are the observed and model fluxes; C is the spectral covariance matrix (see Appendix C); and W is a vector of normalized weights that are proportional to the relative width δλ/λ of each spectral channel or photometric filter. The last prevents the fit from putting too much emphasis on the IFS points (higher density of measurements) at the expense of the photometric points, which cover a wider spectral range (e.g., Ballering et al. 2013; Olofsson et al. 2016). Nevertheless, to test their effect on the fits to the HD 95086 b spectrum, we show in Appendix D a comparison of the best-fit models obtained with and without these additional weighting coefficients.

For each fit the MCMC was run with 100 walkers until the number of steps met a criterion based on autocorrelation time for convergence. The adopted criterion is that the number of steps is 50 times larger than the autocorrelation time for all free parameters (see, e.g., documentation of emcee; Foreman-Mackey et al. 2013). This resulted in 5000–20000 steps for atmospheric models and 100000–250000 steps for the atmospheric+circumplanetary disk models. We then used a “burn-in” factor of 0.5.

Fig. 7 Results of the MCMC orbital fitting of HD 95086 based on the NaCo archive data (epochs 2012-2013) and our SPHERE-IRDIS astrometric results obtained in the K1 band with the SpeCal-TLOCI pipeline (epochs 2015-2019). The orbital parameters are: the orbital period (P (yr)), the eccentricity (e), the inclination (i (°)), the longitude of ascending node (ω (°)), the argument of periastron (ω (°)), the time of periastron passage (Tp (yr AD)). The best orbital fit solution is given by the MAP as black stars in the non-diagonal panels, and as the solid black vertical lines in the diagonal panels. The 1σ confidence region defined as the shortest interval comprising 68% of the probability around the MAP solution is shown with the vertical dashed red lines in the diagonal panels. As ω and ω are degenerated; confidence regions of 34% are given for both. If (ω, ω) is a solution, then (ω + 180°, ω + 180°) is a solution as well. The color bar (on the right) indicates the number of solutions corresponding to a given color for each subplot.

Fig. 8 Color-magnitude diagram considering the SPHERE-IRDIS K1 photometry, and the JIFS photometry from 1.2 to 1.32 µm as extracted from the SPHERE-IFS datacubes.

5.1.1 Atmospheric models

We provided the following public grids of atmospheric models as input to special: BT-SETTL (Allard et al. 2012; Allard 2014), DRIFT-PHOENIX (Woitke & Helling 2003, 2004; Helling & Woitke 2006; Helling et al. 2008), and the grids of A and AE forsterite cloud models with 60 µm modal grain size distribution from Madhusudhan et al. (2011, hereafter M11).

These grids have different prescriptions for clouds, opacity, and dust. The BT-SETTL and DRIFT-PHOENIX models both rely on the PHOENIX atmosphere model (Hauschildt 1992), while the M11 models use a modified version of the TLUSTY atmosphere model (Hubeny 1988; Hubeny & Lanz 1995). The BT-SETTL models consider dust formation using cloud microphysics (condensation and sedimentation mixing timescales) to compute dust grain sizes in a self-consistent way. The DRIFT-PHOENIX models rely on the non-equilibrium cloud model DRIFT and account for the formation of dust grains through a kinetic approach (grain formation, growth, settling, advection, and evaporation). Seven different solids are considered for dust and cloud formation: MgSiO₃[s], Mg₂SiO₄[s], MgO[s], SiO₂[s], SiO[s], Al₂O₃[s], and TiO₂[s]. In the M11 grids, different distributions of the grain sizes and vertical extent are considered. The A and AE cloud models correspond to clouds extending to the top of the atmosphere and half the pressure scale height, respectively. In this work we only consider the M11 grids calculated using forsterite (Mg₂SiO₄[s]) clouds owing to the completeness of the grids. Given the very red slope of the spectrum of HD 95086 b and the better fits obtained with the cloud-A grid (see below), we did not consider the grids with thinner and/or lower cloud distributions presented in M11.

All grids consider effective temperature (T_eff) and surface gravity (log(g)) as free parameters, except for the M11 forsterite cloud-A models (log(g) fixed to 4.0). The DRIFT-PHOENIX and M11 forsterite cloud-A grids also include metallicity (log(Z/Z_⊙)) as a free parameter. The photometric radius R_phot was allowed to take values between 0.1 R_Jup and 5 R_Jup, the lower bound allowing for the possibility of a fraction of the emission not reaching the observer (e.g., in the presence of a large amount of dust around the planet). Table 5 summarizes the sampling of the free parameters for each grid. For each model grid, we carried out two types of fits; the optical extinction A_V was either fixed to 0 (the value estimated for the star; Chen et al. 2012) or left as a free parameter to account for the possible presence of dust around the planet. We considered the Cardelli et al. (1989) extinction law and allowed the value of A_V to span from 0 to 20 mag.

The best-fit models retrieved by special for each grid of atmospheric models are presented in the top four panels of Fig. 9. The solid and dashed green lines correspond to the maximum-likelihood samples among all posterior samples, when A_V is included as a free parameter and set to 0 mag, respectively. In order to visualize the uncertainties on the retrieved physical parameters, we also plotted 200 random samples from the posterior distribution inferred by the MCMC in the free A_V case (cyan curves in Fig. 9). Table 6 reports the most likely parameters for each type of model based on their marginalized distribution (some of which are shown in Figs. D.1–D.4).

For each type of model we computed the Akaike information criterion (AIC; Akaike 1974) to evaluate which models reproduced the observed spectrum better. The definition of the AIC takes into account the number of free parameters in order to avoid overfitting. Low values of AIC correspond to a good fit with a relatively small number of free parameters. For each type of model we report ∆AIC = AIC_model − AIC_BB in Table 6, where AIC_BB is the AIC obtained for a fit to a single blackbody component (no atmospheric model).

We show in Figs. D.1–D.4. the corner plots retrieved by special for some of the models with minimum AIC values, to show the degeneracy between some parameters. When log(g) is a free parameter, an additional panel is shown for the mass posterior distribution estimated from the log(g) and R_phot posterior distributions (i.e., not a free parameter of the fit).

Inspection of Figs. 9 and D.2 and Table 6 reveal two types of models with high support. Depending on the assumed atmospheric grid, either the planet is found to have a relatively low effective temperature (800-1200 K), small to medium amount of extinction (A_V ≲ 10 mag), and super-solar metallicity (DRIFT-PHOENIX and M11 forsterite cloud-A models; e.g., Figs. D.2 and D.3) or it has a higher effective temperature (1200-1600 K) and a high level of extinction by surrounding dust (A_V ≳ 10 mag) for a solar metallicity (BT-SETTL and M11 forsterite cloud-AE models; e.g., Fig. D.1). In particular, we find that only the DRIFT-PHOENIX and M11 forsterite cloud-A models can reproduce the red slope of the spectrum when the extinction is set to A_V = 0 mag (dashed lines in Fig. 9), while the BT-SETTL and M11 cloud-AE models are unable to account for the observed spectrum without extinction (∆AIC > 10; Burnham & Anderson 2002). These two solutions can also be seen in the DRIFT-PHOENIX posterior samples (top left panel of Fig. 9) and corner plot (Fig. D.2) obtained when AV is set as a free parameter: two clusters of solutions can be seen corresponding to low T_eff, high log(g), high log(Z/Z_⊙), and low A_V on the one hand, and high T_eff, unconstrained log(Z/Z_⊙), and high A_V on the other hand.

Our results suggest that a large amount of dust is present, either in the upper part of the atmosphere (super-solar metallicity; see details in Sect. 5.2) and/or around the planet (to account for the extinction). Whether circumplanetary dust could emit an additional thermal component detectable in our spectrum is further investigated in the next section.

Fig. 9 Measured spectrum of HD 95086 b (see Fig. 4 for legend information) compared to the best-fit models retrieved by special for each atmospheric model grid considered in this work (top four panels: DRIFT-PHOENIX, BT-SETTL, Mil cloud-A, Mil cloud-AE), and for different circumplanetary disk models (bottom two panels). All solid lines correspond to extinction AV considered as a free parameter, while the dashed lines correspond to a fixed extinction of AV = 0 mag. The circumplanetary disk models consist of either mixed atmosphere (BT-SETTL or Mil c-A) + extra blackbody models (solid lines) or CPD-only models (dotted lines). The latter correspond to either a debris CPD model (bottom left panel) or a viscous CPD model (bottom right panel; Zhu et al. 2015), where the parameter $\log \left(MM/M_J^2\right)$ corresponds to the mass accretion rate. All types of models have a similar level of support, except for the BT-SETTL and Mil c-AE models without extinction, and the viscous CPD model (Table 6).

5.1.2 Circumplanetary disk models

HD 95086 b is located between two debris disk belts. As it is one of the reddest substellar object known (De Rosa et al. 2016), and a large amount of dust appears to be necessary to account for the observed spectrum (see previous section), we also investigated the possible presence of a circumplanetary disk (CPD) signature in our spectrum. We considered two types of CPDs: a circumplanetary primary viscous disk, in which the exoplanet b still accretes material, and a circumplanetary debris disk consisting of heated grains and modeled by a blackbody component.

We used the SED predictions from the grid of accreting CPD models presented in Zhu (2015), which are characterized by two free parameters: the inner truncation radius of the CPD and the mass accretion rate, spanning 1 to 4 R_Jup and 10⁻⁴ to 10⁻⁷ $M_{\text{Jup}}^2$ yr^-1, respectively. For debris CPDs we considered either a single blackbody component (without atmospheric model), as performed recently for PDS 70 b (Wang et al. 2020; Stolker et al. 2020a), or an additional blackbody component besides the emission from the atmosphere. We allowed the values of the blackbody temperature Tbb to range between 100 K and 2000 K and a blackbody radius between 0.1 and 10 R_Jup, with the condition that Tbb is lower than or equal to T_eq, where T_eq is the equilibrium temperature corresponding to a distance of R_bb (i.e., the extreme case of a spherical shell of optically thick hot dust). If the condition is not met for a particular sample, its log-likelihood is set to minus infinity.

For viscous CPD models the MCMC converged at the edge of the parameter space in terms of inner truncation radius (1.0 R_Jup) for a mass accretion rate of ~10^−6.5 $M_{\text{Jup}}^2$ yr⁻¹. However, the slope of the viscous CPD models is too red to reproduce the observed spectrum on its own (dotted line in bottom right panel of Fig. 9), leading to a poor fit (∆AIC ~ 636). This does not prevent the possibility of a combination of atmospheric and viscous CPD emission. Future measurements in the Hα filter are required to constrain the accretion rate and definitely rule out this hypothesis.

The single-blackbody model leads to a satisfactory fit for a temperature of 1032 ± 45 K and a photometric radius of 1.1 ± 0.2 R_Jup (bottom left panel of Fig. 9 and Fig. D.1). It is one of the models that minimizes the ∆AIC. For atmospheric+debris CPD models, we find that the addition of a blackbody component can also reproduce the observed spectrum without significantly increasing the ∆AIC value. The minor improvement for the maximum-likelihood atmospheric+debris CPD models is such that ∆AIC is still lower than 10, hence implying a similar level of support as models with fewer free parameters (see bold values of ∆AIC in Table 6), including super-solar metallicity atmosphere models. The corner plots associated with atmosphere+debris CPD models (e.g., Fig. D.4 for M11 c-A+debris CPD) show that a number of solutions correspond to negligible contribution from the additional blackbody. However, some high-likelihood solutions involve values of T_bb ~ 800-1300 K and R_bb ~ 0.3–0.7 R_Jup that are comparable to T_eff and R_phot, respectively. This is the case in particular for the maximum-likelihood BT-SETTL+CPD and M11c-A+CPD models, shown as a green curve in the bottom two panels of Fig. 9. Considering the black-body temperature to be the equilibrium temperature of the dust, T_bb ≈ 1319 K (resp. 871 K) for T_eff ≈ 1438 K (resp. 933 K) would imply that the heated dust is located near the top of the atmosphere in either case.

5.1.3 Evolutionary models

Atmosphere modeling of exoplanets and brown dwarfs consists in describing the physical and chemical processes at play in substellar atmospheres using radiative-convective equilibrium models, which can include non-equilibrium chemistry processes, the effect of stellar irradiation, cloud formation, dust settling, and/or mixing to simulate spectra. Even so, their results must be compared to predictions of evolutionary models, which give the evolution of the internal structure of exoplanets and brown dwarfs in time, to exclude non-physical solutions. To do so, based on the apparent photometry, age, and distance of HD 95086 b, we used the Bern EXoplanet cooling tracks (BEX, Marleau et al. 2019) with the AMES-COND atmospheres (Baraffe et al. 2003), corresponding to hot or warm start initial conditions, to derive the predicted bulk properties of the planet (luminosity, mass, effective temperature, surface gravity, and radius). Following a similar approach to Delorme et al. (2017b), we compared the regime of solutions between the two atmosphere and evolutionary models for the predicted surface gravity and radius, as shown in Fig. 10. At the age of HD 95086, BEX models predict for effective temperatures between 800 and 1600 K, typical values of 1.2–1.5 R_Jup and log(g) between 3.6 and 4.2, with a predicted mass of 3–12 M_Jup for HD 95086 b.

From the best-fit solutions of the atmosphere models shown in Fig. 10 and reported in Table 6, the only model consistent with the BEX predictions is M11 with forsterite cloud-A models. The spectrum of HD 95086 b is well reproduced by a super-solar metallicity 800–1200 K atmosphere, with only a small to medium amount of additional extinction required. The atmospheric fits lead to values of radii (${132}_{-0.31}^{+0.16}$ R_Jup) comparable to the physical radius predicted by the BEX models. For all other grids of models, the favored values of photometric radius (0.6–1.0 R_Jup) appear smaller than expected, which may either suggest that the M11 forsterite cloud-A models are the most appropriate for the case of HD 95086 b or that a fraction of the atmospheric flux is obscured by circumplanetary dust, as we discuss below.

5.2 Origin of the red spectral slope of b

5.2.1 Super-solar metallicity atmosphere

The good fit to high-metallicity atmospheric models suggests that HD 95086 b could be somewhat similar to the unusually red L dwarfs (e.g., Gizis et al. 2012; Marocco et al. 2014, and references therein). Looper et al. (2008) and Stephens et al. (2009) suggested that a high metallicity was indeed responsible for the unusually red slope of the field L dwarfs in their samples (for which a low gravity appeared unlikely). A high metallicity facilitates the production of dust grains in the atmosphere, hence clouds. To reproduce the spectrum of unusually red L dwarfs, Marocco et al. (2014) tested different extinction laws corresponding to different dust compositions, namely corundum (Al₂O₃), enstatite (MgSiO₃), iron, and ISM-like (R_v = 3.1). They found that dereddening with any of these extinction laws (including ISM) makes the spectra consistent with that of standard field L dwarfs. Their findings also corroborate our two categories of best-fit models. It appears observationally difficult from the near-infrared spectrum alone to constrain where the dust is located, either in the upper atmosphere (super-solar metallicity enhancing cloud formation) or around the planet (high circumplanetary extinction).

5.2.2 Circumplanetary disk

The first evidence for the presence of a viscous circumplanetary disk was presented in Christiaens et al. (2019) and Isella et al. (2019) for the case of protoplanet PDS 70 b. This is consistent with the estimated young age of the system (~5 Myr) and the presence of a large amount of gas in the protoplanetary disk. Since the HD 95086 system is older (~13.3 Myr), and a low amount of CO has been observed (Booth et al. 2019), the viscous disk is not favored. Our best-fit viscous CPD model also appears too red on its own to account for the observed near-infrared spectrum. Nonetheless, considering that Chinchilla et al. (2021) and Eriksson et al. (2020) find evidence for ongoing accretion onto a planetary-mass object with a main estimated age of 25 Myr and 30–40 Myr, respectively, future measurements in the Hα filter are required to constrain the accretion rate of HD 95086 b, and definitively rule out the possibility of a combined atmosphere+viscous CPD model. We note that the comparison with the studied system by Eriksson et al. (2020) could be nuanced as they studied a binary system of two M stars that could keep the gas-rich disk longer (e.g., known as the Peter Pan disks, Silverberg et al. 2020), and whose age is poorly constrained and debated. It could be younger, with a recent averaged age estimated at about 20 Myr (and an age range of 3–65 Myr, Ujjwal et al. 2020). Finally, if H-alpha measurements are necessary to measure the accretion rate, they could also be a sign of chromospheric activity as it is a common feature in late M and early L objects (Chinchilla et al. 2021).

An alternative explanation for the red slope of HD 95086 b is the presence of circumplanetary dust causing high extinction. In this scenario we can expect the dust located the closest to the planet to be heated to high enough temperatures to show a signature corresponding to an IR excess comparable to the atmospheric emission alone. We find in Sect. 5 that the addition of a second blackbody component to model a circumplanetary debris disk around HD 95086 b can also reproduce the observed spectrum with a similar level of support to models with fewer free parameters (i.e., without an extra blackbody component). In particular, the solution from the model of atmosphere M11 with cloud-A and an extra blackbody component and optical extinction A_V is consistent with predictions from evolutionary models. This model favors a super-solar metallicity (${0.5}_{-0.3}^{0.4}$) with a medium level of extinction ($A_{\text{v}}={8.9}_{-0.3}^{+0.4}$ mag). Hence, both a circumplanetary disk and a super-solar metallicity could account for the red spectral slope of b, and both could be present together; if there is a debris CPD, some of this debris might have been accreted onto the planet’s atmosphere, which naturally increases the metallicity.

A small subset of the debris disk CPD solutions from models of atmosphere is expected to have a non-negligible signature at near-IR wavelengths (solutions with large T_bb and R_bb similar to R_phot, as shown in the bottom panels of Fig. 9). In these cases, considering the blackbody temperature to correspond to the equilibrium temperature at the separation of the dust implies that the heated dust would be located near the top of the atmosphere (within ~1 R_Jup distance above the atmosphere). However, since a significant fraction of posterior samples corresponds to no significant excess at near-IR wavelengths, it is also possible that only cold circumplanetary dust is present. We note that Pérez et al. (2019) used ALMA 1.3 mm observations from Su et al. (2017) to search for the presence of a circumplanetary disk around HD 95086 b, and derived an upper limit of 30 µJy at the planet location. All the models retrieved by special are consistent with this non-detection, being over three orders of magnitude fainter than their upper limit.

To test whether longer infrared wavelength observations may allow us to distinguish between high-metallicity cold circumplanetary dust and hot+cold circumplanetary dust, we show in Fig. 11 the predictions at longer wavelengths for the highest-likelihood models of each type. The predicted spectra for six of the most likely models are reported up to 11 µm: two models with high metallicity only and no extinction (in cyan), two models with mid to high extinction (i.e., cold circumplanetary dust), and two models with both extinction and an extra blackbody component (i.e., both hot and cold circumplanetary dust). We see that multiple accurate measurements at longer wavelengths, possibly including a spectrum, would be required to distinguish between the different scenarios since the predicted fluxes are still relatively similar for the different scenarios.

On the other hand, observations at short wavelengths in the visible and the near-infrared may also reveal the presence of a CPD. If there is a similar disk around HD 95086 b, it should be polarized, and also detectable at short wavelengths based on Mie’s theory, for example by ZIMPOL and SPHERE-IRDIS in dual-polarization imaging mode (van Holstein et al. 2021). In particular, we note the excess of flux at the shorter wavelengths (below 1.4 µm), with respect to the best-fit models (within 1–2σ) in Fig. 11, which would be consistent with a debris disk. However, this CPD should be modeled with a more tuned model than the additional blackbody component used in this work to account for the observed spectrum.

Future ground-based observations with ERIS at the VLT, METIS at the ELT, and JWST in space, should soon enable us to unambiguously confirm the origin of the very red spectrum of HD 95086 b. This is indeed highlighted by synthetic observations of CPD from Szulágyi et al. (2019, for the ERIS instrument) and Chen & Szulágyi (2021, for the JWST and ELT telescopes).

Fig. 10 Comparison of the solutions for the different models of atmosphere (left: BT-SETTL, middle left: DRIFT-PHOENIX, middle right: M11 cloud-A, and right: M11 cloud-AE) and models of evolution (top panel: BEX-Hot and bottom panel: BEX-Warm). As for the atmospheric models, the best-fit solutions are shown in pink, with a pentagon if the optical extinction AV = 0, a square if the extinction is a free parameter, and a triangle if there is an additional blackbody fitted; the samples of solutions are in red if the extinction AV = 0 and in blue if the extinction is a free parameter. For the evolutionary models, the solid curves represent the surface gravity as a function of the planetary radius given different effective temperature (blue for the coolest to red for the hottest). The colored circles represent the expected log g and R at a given age (10 or 20 Myr): hence lower and upper values expected for HD 95086 aged 13.3 Myr.

6 Searching for planet c

The architecture of the young planetary system HD 95086 offers an interesting comparison case with two emblematic systems HR 8799 (Götberg et al. 2016; Su et al. 2015) and PDS 70 (Keppler et al. 2018), and more generally with the interpretation that these multiple-belt debris disks are young analogs to our Solar System. The imaged giant planets would be responsible for the dynamical clearing of the debris disks and the formation of observed multiple-belt architecture, as suggested by Kennedy & Wyatt (2014) and Shannon et al. (2016). For HD 95086 the observed planet-belt architecture composed of a warm and relatively narrow inner belt at ~8 au, a broad cavity from typically 10-100 au inside which the massive (4–5 M_Jup, a ~ 53 ± 16 au and e ≤ 0.17, see Table 4) planet HD 95086 b orbits, and finally, a cold outer belt lying between 106 and 320 au, suggests that probably more than one giant planet is orbiting in this system (Su et al. 2017; Rameau et al. 2016; Chauvin et al. 2018).

Applying Eqs. (4) and (5) of Shannon et al. (2016) to the case of HD 95086 (1.6 M_⊙, 14.3 ± 2 Myr, cavity from 10 to 100 au), Chauvin et al. (2018) derived a minimum mass of the planets in the cavity of 0.35 M_Jup and a typical number of required planets of 2.4 (i.e., 2–3 giant planets depending on their respective separation). Comparing these results to the outcome of HARPS and SPHERE combined detection limits (up to 2017) in the context of a planet–disk coplanar configuration, they found that there might still be room for two additional stable planets c and d in the cavity in addition to b with typical masses between 0.35 M_Jup (dynamical clearing constraint) and 6 M_Jup for a semimajor axis between 10 and 30 au or 0.35 M_Jup and 5 M_Jup beyond 30 au.

Considering the most recent SPHERE multi-epoch observations (up to May 2019) and including the deepest high-contrast images obtained so far with SPHERE (in early 2018), we revisit here the search for planet c (and d) beyond the inner warm narrow belt at ~8 au, that is beyond ~100 mas for the HD 95086 system. We first push the current exploration of the close environment to search for c using the K-Stacker algorithm. With no clear detection, we then set new upper limits on the potential masses of inner giant planets in the system.

Fig. 11 Spectrum of HD 95086 b compared to the most likely models retrieved by special with the highest-likelihood atmospheric and CPD models, extended up to 11 µm. Cyan curves show models with high metallicity but no extinction: DRIFT-PHOENIX (solid line) and Mil cloud-A (dashed line). Green curves show models with mid to high extinction: BT-SETTL (solid line) and Mil cloud-AE (dashed line). Black curves show models with both extinction and an extra blackbody component: BT-SETTL (solid line) and Mil cloud-A (dashed line). There is no atmospheric model+viscous CPD reported here as their best model does not represent a good enough fit to the data with respect to the other models.

6.1 K-Stacker exploration

The existence of HD 95086 c has been investigated by Le Coroller et al. (2020) using the K-Stacker algorithm (Nowak et al. 2018). The code was applied to a combination of six IFS observations (median over Y-H wavelengths) obtained at different epochs (from February 2015 to May 2017). Following a similar strategy, the complete set of observations reported in Table 2 is now considered for both IFS and IRDIS. They represent a total of ten epochs between February 2015 and May 2019. Based on the detection limits achieved with SpeCal-PCAPad and SpeCal-TLOCI in the H and K1 bands with IFS and IRDIS, respectively (see Fig. 3), we selected the eight best ones, rejecting February 2015 and January 2016. For IFS, the choice of H band (median over the channels between 1.47 and 1.59 µm) over Y and / bands was motivated by deeper sensitivities down to masses of 2 M_Jup given the very red colors (mid-L to L/T types) expected for the detectable planets according to the SPHERE detection limits. As fixed parameters for K-Stacker, the recent Gaia-DR2 distance and the primary stellar mass derived in this work (see Sect. 4) were used, together with an exploration range of orbital parameters for HD 95086 c compatible with stable dynamical orbits considering the system architecture and the presence of HD 95086 b (see Table 7).

For both K-Stacker runs, IRDIS at the K1 band and ΓFS at the H band, we do not detect any clear point-source signal with a signal-to-noise ratio higher than 5, which could indicate the probable presence of a closer-in planet. Further characterization to improve the detection limits in the 100–300 mas regime using either reference differential imaging with star-hopping (Wahhaj et al. 2021) or molecular mapping techniques (Petrus et al. 2021) will be needed.

6.2 Detection probabilities

Without any clear detection of additional planets orbiting HD 95086 on individual epochs or using K-Stacker and all the best epochs available, we now explore the completeness of previous HARPS and new SPHERE observations using the pyMESS2 code (A.-M. Lagrange, priv. comm.), a Pythonized version of the MESS2 code (Lannier et al. 2017), together with the Exoplanet Detection Map Calculator (Exo-DMC Bonavita 2020)⁴. Both codes are the latest versions of the Multi-purpose Exoplanet Simulation System (MESS; Bonavita et al. 2012), a Monte Carlo tool for the statistical analysis of direct imaging survey results. In a similar fashion to its predecessors, both codes combine the information on the target stars with the instrument detection limits to estimate the probability of detection of a given synthetic planet population, ultimately generating detection probability maps.

For each star in the sample, they generate a grid of masses and physical separations of synthetic companions, then estimate the probability of detection given the provided detection limits at each epoch. The default setup uses a flat distribution in log space for both the mass and semimajor axis, but in a similar fashion to their predecessors, they allow for a high level of flexibility in terms of possible assumptions on the synthetic planet population to be used for the determination of the detection probability. For each point in the mass-semimajor axis grid, pyMESS2 and DMC generate a fixed number of sets of orbital parameters. By default all the orbital parameters are uniformly distributed except for the eccentricity, which is generated using a Gaussian eccentricity distribution with µ = 0 and σ = 0.3 (constraint: e ≥ 0), following the approach by Hogg et al. (2010) (see Bonavita et al. 2013, for details). This allows us to properly take into account the effects of projection when estimating the detection probability using the contrast limits in Fig. 3. They calculate the projected separations corresponding to each orbital set for all the values of the semimajor axis in the grid (see Bonavita et al. 2012, for a detailed description of the method used for the projection). This allows us to estimate the probability that each synthetic companion is truly in the instrument field of view and therefore that it will be detected if the value of the mass is higher than the limit. In this specific case, we chose to restrict the inclination and the longitude of the node of each orbital set to make sure that all companions in the population would lie in the same orbital plane as the disk (i.e., with i = −30 ± 3° and ω = 97 ± 3°) for one run, and consider no a priori information on the system’s orientation for a second run. We also assumed a uniform distribution for the eccentricity, with a maximum value of 0.6.

While MESS was limited in its use to direct imaging, both DMC and pyMESS2 can also be used to draw similar constraints using other kinds of datasets, including radial velocity (RV) data. Given the provided RV time series, both codes use the local power analysis (LPA) approach described by Meunier et al. (2012) to estimate, for each mass and separation in the grid, for what fraction of the generated orbital sets the signal generated by the companion would be compatible with the data. DMC at this stage is limited to an independent determination of the detection probability for each technique that is later combined. The two sets of maps are merged by considering, for each point in the grid, the best value of the probability. With pyMESS2 (and MESS2), the detection probability directly combines the detectability of each generated planet by checking if each planet can be detectable at least in one DI epoch or in the RV epochs. The probability is then derived by counting, for each [mass,a], how many of the generated planets are detected with either technique. In the end, combining the two methods allows a more accurate determination of the detection probability.

Figure 12 (left) illustrates the advantages of combining HARPS radial velocity observations with direct imaging, as shown by Chauvin et al. (2018), but here updating their results with the latest SPHERE measurements. The detection probability gain with the new SPHERE epochs is shown in the middle (coplanar case) and right (no orbital constraints) panels for the specific region of interest between 5 and 35 au, where the presence of inner giant planets is suspected. The contrast gain of about 1.0 mag at typically 100–200 mas with both IFS and IRDIS presented in Fig. 3, with the multi-epoch combination, enables us to nail down the detection probability map by a fraction of Jupiter mass and au. A further relevant gain would likely imply using alternative observing strategies, such as star-hopping (Wahhaj et al. 2021) or molecular mapping (Hoeijmakers et al. 2018; Petrus et al. 2021), before the arrival of the extremely large telescopes (Chauvin 2018).

Fig. 12 Probability of companion detection as a function of the distance to the star from pyMESS2. On the left, overview of the detection limit results based on radial velocity (≤ 5 au) and direct imaging (all SPHERE+NaCo observations, ≥5 au) data for the whole HD 95086 system. In the middle and right, direct imaging results are zoomed between 5 and 35 au. In the left and middle, coplanarity is assumed between the researched exoplanets and the outer belt (i.e., with i = −30 ± 3° and ω = 97 ± 3°), but not on the right. The results without the last five epochs imaged with SPHERE (in 2018 and 2019) are shown in gray to black contours for a given probability of 10%, 50%, 68%, 90%, and if coplanarity is assumed 99%. The color bars show the levels of probability detection.

7 Conclusion

In this work we presented and analyzed five new observations from SPHERE (2018–2019) on the young Solar System analog HD 95086, as well as a re-analysis of the five previous SPHERE observations (2015–2017), and the results from archival data from the GPI and NaCo instruments. We reported an in-depth characterization of the system HD 95086:

1. Regarding the exoplanet HD 95086 b, we extracted for the first time its spectrum in the J and H bands by combining six epochs imaged with SPHERE-IFS to maximize the signal-to-noise ratio of the planet as it is hardly detectable in these bands.

2. We constrained the physical properties of HD 95086 b providing spectroscopic and photometric measurements in J, H, K1, K2, and L′ bands from the SPHERE-IFS, GPI, SPHERE-IRDIS, and NaCo instruments. We obtained two types of solutions: for a surface gravity log(g) between 3.3–5.0, either the exoplanet seems to have a high effective temperature between 1400–1600 K and a significant extinction (≳10), or lower temperatures between 800–1200 K with a small to medium level of extinction (≲10) and a super-solar metallicity. Both of these solutions reveal the presence of significant dust, which can explain the redness of the exoplanet b. The dust could be present either in the upper layers of the atmosphere, explaining the super-solar metallicity atmosphere found for some atmospheric model best solutions, or around the exoplanet HD 95086 b with the presence of a circumplanetary disk, explaining the high extinction found for the other atmospheric models.

3. Additional modeling combining at the same time models of atmospheres and circumplanetary disks confirm the possibility of the presence of a debris circumplanetary disk around the exoplanet HD 95086 b since the solution is as likely as the other solutions found without it based on the Akaike information criterion. Considering the age of the system, the nature of this possible circumplanetary disk suggested by our modeling is more likely to be a debris disk than a viscous disk. Nevertheless, future measurements sufficiently precise in the Hα filter are required to constrain the accretion rate of HD 95086 b and rule out definitively the possibility of a viscous circumplanetary disk combined with specific atmospheric models. Other future measurements, particularly in M and N bands, are necessary to discriminate between our remaining best-fit atmospheric models.

4. We updated the orbital parameters for the exoplanet b, adding two additional monitoring years from SPHERE. Our best orbital solutions are consistent with previous published orbital parameters.

5. As for additional exoplanets in the system, we pushed the detection performance by combining the best epochs using the K-Stacker algorithm, which combines the best observations through different Keplerian motions to correct for the orbit of any additional exoplanet. We also applied several post-processing algorithms, but we did not find any robust candidates. Nonetheless, we put new constraints on the masses and locations of putative additional exoplanets in the system, and we ruled out any other 5 M_Jup inner planet in the system located at a distance greater than 17 au at a 50% confidence level (or 9 M_Jup inner planet at a distance greater than 10 au at a 50% confidence level).

Future observations with the JWST (GTO target), at the VLT/I with GRAVITY, ERIS, SPHERE and its potential upgrade, and with the first light instruments of the ELT, should enable us to understand the global architecture and origin of HD 95086, and its commonality with our own Solar System.

Acknowledgements

We would like to thank Kate Su for sharing the ALMA continuum image used in Fig. 1 together with the SPHERE-IRDIS high-contrast image at K1-band. We acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU. The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). This work is partly based on data products produced at the SPHERE Data Center hosted at OSUG/IPAG, Grenoble. 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 Zurich (Switzerland), NOVA (Netherlands), ONERA (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 project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COBREX; grant agreement n° 885593). V.C. acknowledges funding from the Australian Research Council via DP180104235. A.V. and G.P.P.L.O. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 757561). Finally, we would like to thank the anonymous referee and the editor for their helpful comments.

Appendix A Observational conditions

Figure A.1 shows the observation conditions for all the available SPHERE epochs used.

Fig. A.1 Distribution of the observational conditions (when available) for all the epochs: Strehl ratio, atmospheric coherence length r0, wind, seeing, and atmospheric coherence time τ0. The Strehl ratio, atmospheric coherence length, and wind were all measured by SPARTA, the computer of the adaptive system of SPHERE, while both the seeing and atmospheric coherence time were measured by the DIMM telescope at Paranal. A wind ≤ 3 m · s−1 indicates the presence of low wind effect (Milli et al. 2017), while a coherence time ≤ 3 ms indicates the presence of wind-driven halo (Cantalloube et al. 2020).

Appendix B Reduced images with SpeCal

Figure B.1 shows the images reduced with the pipelines SpeCal-TLOCI and SpeCal-PCAPad for the IRDIS and IFS data, respectively.

Fig. B.1 Signal-to-noise ratio maps for all the SPHERE epochs (top: IRDIS in the K1 band reduced by the pipeline SpeCal-TLOCI, and bottom: IFS in the YJH bands) reduced by the pipeline SpeCal-PCAPad (PCA ADI+SDI). The color bar corresponds to the signal-to-noise ratio. The region below the inner working angle of the coronagraph is masked.

Appendix C Spectral correlation matrix

Figure C.1 shows the spectral correlation matrix (ϕ) empirically estimated from our IFS datasets, as in Greco & Brandt (2016) and De Rosa et al. (2016). The spectral covariance matrix C used for spectral fits is obtained as C_ij = ϕ_ijσ_iσ_j, where σ_i and σ_j correspond to the reported uncertainties in flux measurements for points i and j.

Fig. C.1 Spectral correlation matrix computed from our SPHERE/IFS (first 15 points) and GPI/IFS (points 15–43) data. For IRDIS and NACO, Cij = δij, where δij is the Kroenecker symbol.

Appendix D Corner plots of the posteriors of the atmospheric models

In this appendix we show the posteriors of the atmospheric models retrieved by special for some of the best models at reproducing the observed spectrum: the single-blackbody model and the BT-SETTL model with free extinction (both in Fig. D.1); the DRIFT-PHOENIX model with free extinction (Fig. D.2); and the M11 cloud-A model with free extinction, without and with an extra blackbody component (Figs. D.3 and D.4, respectively). In each case the mass is not a free parameter of the fit, but is calculated from the log(g) and R_p posterior distributions.

We note that in Eq. (2), which represents the log-likelihood expression provided to the sampler, in addition to the spectral covariance matrix, additional weighting coefficients are used. These coefficients can account for the fact that absolute flux calibration is performed instrument per instrument (i.e., it is not specific to each spectroscopic point), whose individual uncertainties and covariances are not properly accounted for in the fit. In the hypothetical case where the uncertainty on absolute flux calibration is the same for the spectrograph and the photometers involved in acquiring the combined spectrum, but dominates over individual statistical uncertainties, the absence of additional weights would typically pin the absolute flux calibration to that of the spectrograph, given the significantly larger number of data points.

This absolute flux calibration uncertainty can be difficult to assess. It was noted when comparing Spitzer/IRS versus Spitzer/MIPS measurements, a case for which the use of additional weights was suggested (e.g., Ballering et al. 2013; Chen et al. 2014). However, it is known to still be an issue when comparing measurements obtained by the latest generation of high-contrast imagers. The amplitude of the systematic bias between the overlapping SPHERE and GPI spectra of HD 206893 B is on the same order as the total estimated (statistical+systematic) uncertainties (Kammerer et al. 2021). This kind of bias has also been noted when comparing the spectrum of PDS 70 b obtained by SINFONI and SPHERE (Christiaens et al. 2019).

Assigning weights such that their sum is 1 for all measurements obtained by a given instrument (i.e., 1 per photometric measurement and 1/N_spec for each spectrograph measurement, where N_spec is the number of channels) would correspond to assigning the same level of confidence for the absolute flux calibration of each instrument. However, this would too severely undermine the significance of the spectroscopic data in the fit. It would indeed make all measurements of a given spectrograph as significant as a single photometric measurement, while spectrographs provide more than just an absolute flux measurement; they also provide information about the shape of the spectrum.

Therefore, as a compromise we opted for weights that are proportional to the spectral bandwidth, as in Ballering et al. (2013) and Olofsson et al. (2016). Some estimated model parameters are dependent on the absolute flux calibration (e.g., photometric radius), on the features, and/or shape of the spectrum (e.g., surface gravity), or on both together (effective temperature and extinction). Therefore, a certain choice of weights may in principle change the estimated best-fit parameters.

To test in practice the effect of the additional weighting coefficients on the fits to the HD 95086 b spectrum, we show in Figs. D.4 and D.5 the best-fit solutions obtained with and without the additional weighting coefficients, respectively. The best-fit parameters show only minor differences, which may have been expected given the good match between most near-IR spectroscopic and photometric points of HD 95086 b (Fig. 4). The largest relative differences correspond to the estimated radius and mass of the planet, although all parameters still remain within the 1σ error bars of the estimates obtained with and without weights.

Fig. D.1 Corner plots retrieved by special for the single-blackbody model at the top right, BT-SETTL with free extinction AV model at the bottom left.

Fig. D.2 Corner plot retrieved by special for the DRIFT-PHOENIX model with free extinction AV.

Fig. D.3 Corner plot retrieved by special for the M11 cloud-A model with free extinction AV.

Fig. D.4 Corner plot retrieved by special for the M11 cloud-A model with free extinction AV and with an extra blackbody component.

Fig. D.5 Corner plot retrieved by special for the M11 cloud-A model with free extinction AV and with an extra blackbody component without additional weighting coefficients.

References

All Tables

All Figures

Fig. 1 Composite ALMA-continuum (at 1.3 mm) and SPHERE/IRDIS (at 2.1 µm) observations of HD 95086 (see Su et al. 2017; Chauvin et al. 2018). The exoplanet b is detected at the K1 band (red dot). The white dashed ring at 180 au represents the peak location of the outer cold belt located from 106 ± 6 au to 320 ± 20 au. The inner warm belt is not resolved with ALMA. The pink ellipse represents the ALMA synthetic beam.

In the text

Fig. 2 Signal-to-noise ratio maps for all the SPHERE epochs. At the top, SPHERE-IRDIS in the K1 band and at the bottom, SPHERE-IFS in the YJH bands, both reduced by the pipeline ANDROMEDA. The color bar corresponds to the signal-to-noise ratio, 0–15 for IRDIS and 0–5 for the IFS. Planet b is located in the white square. The region below the inner working angle of the coronagraph is masked. In the top right corner of the IFS data, the one (or two) point-like feature(s) correspond to the remanent of the star on the IFS detector.

In the text

Fig. 3 Detection limits at 5σ for IRDIS in the K1 band on the left and IFS in the YJH bands on the right. Data are reduced at the top by the pipeline SpeCal-TLOCI (for IRDIS) or SpeCal-PCAPad (for IFS) and at the bottom by the pipeline ANDROMEDA. The gray area corresponds to the area hidden by the larger coronagraph used for all the epochs (i.e., N_ALC_Ks) which has an inner working angle of 0.107″ and 0.116″ for IFS and IRDIS, respectively. The blue vertical line corresponds to the position of the planet HD 95086 b at the best epoch January 2018 (cyan, 0.625″).

In the text

Fig. 4 Spectroscopic and photometric values of the exoplanet HD 95086 b expressed in contrast relative to the host star from 1.2 µm to 3.8 µm with SPHERE-IFS (blue triangles), SPHERE-IRDIS (yellow squares), GPI (pink diamonds), and NaCo (red circle) in the bands JH, K, K12, and L′ at the top. At the bottom, the transmission curves of the filters K1, K2, and L′ are shown1. Results are extracted from an averaged stacking of the six best epochs for the SPHERE data, and for one single epoch for the GPI and NaCo data.

In the text

	Fig. 5 Spectrum of the exoplanet HD 95086 b expressed in contrast relative to the host star from SPHERE-IFS data in the J band (1.21–1.32 µm) and in the H band (1.46–1.62 µm) at a 68% confidence level. The contrast is given in units of 10−6. The data are reduced by the pipeline SpeCal-PCA ADI where a six-epoch averaged stack was made before extracting the spectrum to boost the S/N.
In the text

Fig. 6 Astrometric positions of the planet HD 95086 b between 2012 and 2020 with three different instruments: SPHERE (yellow squares), GPI (pink diamonds), and NaCo (red circles) with 1σ error bars. The astrometric positions from SPHERE were computed with the SpeCal-TLOCI pipeline. The letter “w” in the legend indicates the observations imaged with satellite spots that enable finding the exact position of the star during the whole sequence and later recenter the frames if necessary. The green and orange solid lines correspond to a sample of the orbital solutions found by the MCMC orbital fit and the K-Stacker tools, respectively. The black dashed and dotted lines respectively represent the MCMC and K-Stacker orbits for which the corresponding orbital parameters are given in Table 4. For the MCMC tool, the orbit corresponds to the MAP (see text), whereas for K-Stacker it is the orbit closest to the mean of the orbital solutions. Both of these orbits are by construction true orbital solutions. Notes: The Feb2015 point from SPHERE is hidden by the two 2015 points from GPI; the two 2019 points are very close to each other. The insert (in the top right corner) shows the location of the planet relative to the cold outer belt and warm inner belt from Su et al. (2017).

In the text

Fig. 7 Results of the MCMC orbital fitting of HD 95086 based on the NaCo archive data (epochs 2012-2013) and our SPHERE-IRDIS astrometric results obtained in the K1 band with the SpeCal-TLOCI pipeline (epochs 2015-2019). The orbital parameters are: the orbital period (P (yr)), the eccentricity (e), the inclination (i (°)), the longitude of ascending node (ω (°)), the argument of periastron (ω (°)), the time of periastron passage (Tp (yr AD)). The best orbital fit solution is given by the MAP as black stars in the non-diagonal panels, and as the solid black vertical lines in the diagonal panels. The 1σ confidence region defined as the shortest interval comprising 68% of the probability around the MAP solution is shown with the vertical dashed red lines in the diagonal panels. As ω and ω are degenerated; confidence regions of 34% are given for both. If (ω, ω) is a solution, then (ω + 180°, ω + 180°) is a solution as well. The color bar (on the right) indicates the number of solutions corresponding to a given color for each subplot.

In the text

	Fig. 8 Color-magnitude diagram considering the SPHERE-IRDIS K1 photometry, and the JIFS photometry from 1.2 to 1.32 µm as extracted from the SPHERE-IFS datacubes.
In the text

Fig. 9 Measured spectrum of HD 95086 b (see Fig. 4 for legend information) compared to the best-fit models retrieved by special for each atmospheric model grid considered in this work (top four panels: DRIFT-PHOENIX, BT-SETTL, Mil cloud-A, Mil cloud-AE), and for different circumplanetary disk models (bottom two panels). All solid lines correspond to extinction AV considered as a free parameter, while the dashed lines correspond to a fixed extinction of AV = 0 mag. The circumplanetary disk models consist of either mixed atmosphere (BT-SETTL or Mil c-A) + extra blackbody models (solid lines) or CPD-only models (dotted lines). The latter correspond to either a debris CPD model (bottom left panel) or a viscous CPD model (bottom right panel; Zhu et al. 2015), where the parameter $\log \left(MM/M_J^2\right)$ corresponds to the mass accretion rate. All types of models have a similar level of support, except for the BT-SETTL and Mil c-AE models without extinction, and the viscous CPD model (Table 6).

In the text

Fig. 10 Comparison of the solutions for the different models of atmosphere (left: BT-SETTL, middle left: DRIFT-PHOENIX, middle right: M11 cloud-A, and right: M11 cloud-AE) and models of evolution (top panel: BEX-Hot and bottom panel: BEX-Warm). As for the atmospheric models, the best-fit solutions are shown in pink, with a pentagon if the optical extinction AV = 0, a square if the extinction is a free parameter, and a triangle if there is an additional blackbody fitted; the samples of solutions are in red if the extinction AV = 0 and in blue if the extinction is a free parameter. For the evolutionary models, the solid curves represent the surface gravity as a function of the planetary radius given different effective temperature (blue for the coolest to red for the hottest). The colored circles represent the expected log g and R at a given age (10 or 20 Myr): hence lower and upper values expected for HD 95086 aged 13.3 Myr.

In the text

Fig. 11 Spectrum of HD 95086 b compared to the most likely models retrieved by special with the highest-likelihood atmospheric and CPD models, extended up to 11 µm. Cyan curves show models with high metallicity but no extinction: DRIFT-PHOENIX (solid line) and Mil cloud-A (dashed line). Green curves show models with mid to high extinction: BT-SETTL (solid line) and Mil cloud-AE (dashed line). Black curves show models with both extinction and an extra blackbody component: BT-SETTL (solid line) and Mil cloud-A (dashed line). There is no atmospheric model+viscous CPD reported here as their best model does not represent a good enough fit to the data with respect to the other models.

In the text

Fig. 12 Probability of companion detection as a function of the distance to the star from pyMESS2. On the left, overview of the detection limit results based on radial velocity (≤ 5 au) and direct imaging (all SPHERE+NaCo observations, ≥5 au) data for the whole HD 95086 system. In the middle and right, direct imaging results are zoomed between 5 and 35 au. In the left and middle, coplanarity is assumed between the researched exoplanets and the outer belt (i.e., with i = −30 ± 3° and ω = 97 ± 3°), but not on the right. The results without the last five epochs imaged with SPHERE (in 2018 and 2019) are shown in gray to black contours for a given probability of 10%, 50%, 68%, 90%, and if coplanarity is assumed 99%. The color bars show the levels of probability detection.

In the text

Fig. A.1 Distribution of the observational conditions (when available) for all the epochs: Strehl ratio, atmospheric coherence length r0, wind, seeing, and atmospheric coherence time τ0. The Strehl ratio, atmospheric coherence length, and wind were all measured by SPARTA, the computer of the adaptive system of SPHERE, while both the seeing and atmospheric coherence time were measured by the DIMM telescope at Paranal. A wind ≤ 3 m · s−1 indicates the presence of low wind effect (Milli et al. 2017), while a coherence time ≤ 3 ms indicates the presence of wind-driven halo (Cantalloube et al. 2020).

In the text

	Fig. B.1 Signal-to-noise ratio maps for all the SPHERE epochs (top: IRDIS in the K1 band reduced by the pipeline SpeCal-TLOCI, and bottom: IFS in the YJH bands) reduced by the pipeline SpeCal-PCAPad (PCA ADI+SDI). The color bar corresponds to the signal-to-noise ratio. The region below the inner working angle of the coronagraph is masked.
In the text

	Fig. C.1 Spectral correlation matrix computed from our SPHERE/IFS (first 15 points) and GPI/IFS (points 15–43) data. For IRDIS and NACO, Cij = δij, where δij is the Kroenecker symbol.
In the text

	Fig. D.1 Corner plots retrieved by special for the single-blackbody model at the top right, BT-SETTL with free extinction AV model at the bottom left.
In the text

	Fig. D.2 Corner plot retrieved by special for the DRIFT-PHOENIX model with free extinction AV.
In the text

	Fig. D.3 Corner plot retrieved by special for the M11 cloud-A model with free extinction AV.
In the text

	Fig. D.4 Corner plot retrieved by special for the M11 cloud-A model with free extinction AV and with an extra blackbody component.
In the text

	Fig. D.5 Corner plot retrieved by special for the M11 cloud-A model with free extinction AV and with an extra blackbody component without additional weighting coefficients.
In the text