© ESO 2021

1 Introduction

While the number of known exoplanets has grown to over 4000 in the last decade¹, the sample of directly imaged substellar companions remains small (e.g., Bowler 2016). These objects are prime targets for a direct study of their atmospheric properties and composition through imaging and low-resolution spectroscopy (e.g., Biller & Bonnefoy 2018), though. Together with evolutionary tracks and atmospheric models, this enables their effective temperature, radius, surface gravity, age, and mass to be inferred and conclusions to be drawn on their formation history and subsequent evolution (e.g., Bowler 2016; Biller & Bonnefoy 2018). Moreover, astrometric measurements from direct imaging enable deriving the orbital parameters of substellar companions and studying the dynamical interactions with their environment (e.g., Kley & Nelson 2012).

Recently, GRAVITY Collaboration (2019, 2020) used long-baseline interferometry to perform medium-resolution (R ~500) K-band spectroscopy of exoplanets. They were able to demonstrate astrometric measurements with a precision ~ 10 times betterthan what had previously been possible and constrain the atmospheric C/O ratio of the gas giant β Pic b. This, together with a chemical abundance framework of its protoplanetary disk based on Öberg et al. (2011), enabled them to infer a formation by warm-start core accretion (Pollack et al. 1996) between the water and carbon-dioxide icelines. Furthermore, Nowak et al. (2020) were able to directly detect β Pic c, another gas giant in the same system previously discovered with the radial velocity technique (Lagrange et al. 2019). They showed that the precise mass estimate from the radial velocity data together with the K-band spectrum also favors a formation by warm-start core accretion (e.g., Mordasini 2013) for β Pic c. In this regard, near-infrared interferometry significantly advances the field of planet formation and evolution by enabling direct observationsof exoplanets discovered with the radial velocity technique for the first time ever.

Another system for which both direct observations and radial velocity data are available is HD 206893. This debris disk system, at a distance of 40.8 pc (Gaia Collaboration 2018), hosts a directly detected substellar companion at a separation of ~ 10 au, HD 206893 B (Milli et al. 2017). Astronomers are puzzled by the nature of this companion due to its unusually red near-infrared color. Delorme et al. (2017) found that an additional K-band extinction of ~ 0.5 mag is required to match the spectrum of HD 206893 B with those of other dusty, low-gravity, or young brown dwarfs. Furthermore, Delorme et al. (2017) and Ward-Duong et al. (2021) showed that its extremely red color together with its very shallow 1.4 μm water absorption feature are challenging to fit with currently available atmospheric models without an additional extra-photospheric source of extinction.

While Delorme et al. (2017) argue for HD 206893 B being an extremely dusty 15–30 M_Jup L-dwarf, consistent with their age estimate of 50–700 Myr for the host star, Ward-Duong et al. (2021) note that its H- and K-band spectra suggest a lower surface gravity and younger object. Together with the high infrared excess of the disk and a possible Argus moving group membership (membership probability ~ 61%; Ward-Duong et al. 2021), there is a consistent scenario for HD 206893 B being a young (< 50 Myr) gas giant planet. This scenario is also supported by the dynamical mass estimate of ${10}_{-4}^{+5} M_{\text{Jup}}$ from Grandjean et al. (2019) based on radial velocity data and the Gaia proper motion anomaly of the system. Their analysis also points to a second, closer-in companion at a separation of ~ 1.4–2.6 au (HD 206893 C). The picture is further complicated by a gap in the debris disk of the system, which was recently discovered using the Atacama Large Millimeter/Submillimeter Array and could be carved by a third, Jupiter-mass companion at ~ 74 au (Marino et al. 2020). A schematic of the HD 206893 system is shown in Fig. 1 for illustrative purposes.

Here, we present Very Large Telescope Interferometer (VLTI)/GRAVITY K-band spectroscopy of HD 206893 B. From these GRAVITY data, we extract the astrometry of HD 206893 B with a precision of ~ 100 μas and a medium-resolution (R ~ 500) K-band contrast spectrum, which we convert to a spectrum of HD 206893 B with a model spectrum of its host star (cf. Sect. 2). From the astrometry, we improve the constraints on the orbital parameters of HD 206893 B (cf. Sect. 3.1) and on the mass and position of HD 206893 C, also using the Gaia proper motion anomaly of the system. Furthermore, we perform atmospheric model fitting, with and without additional extra-photospheric extinction by high-altitude dust clouds made of enstatite grains (cf. Sect. 4.1). A similar study investigating the dynamical mass and the circumplanetary accretion flux of the PDS 70 b and c protoplanets using GRAVITY data has recently been published by Wang et al. (2021). Moreover, we perform spectral retrievals for HD 206893 B (cf. Sect. 4.2) and check for consistency between our best fit atmospheric parameters and evolutionary tracks (cf. Sect. 4.3). Finally, we discuss our findings in the context of previous works on this system (cf. Sect. 5).

Fig. 1 Schematic of the HD 206893 system with the two inner companions “B” and “C,” the debris disk with its ~ 27 au wide gap, and the planetary-mass companion candidate at ~74 au that could be responsible for clearing this gap. The quoted mass estimates for HD 206893 B and C are explained in Sect. 5.4 of this paper. Distances and sizes are not to scale.

2 Observations and data reduction

As part of the ExoGRAVITY large program (Lacour et al. 2020), we obtained two epochs of GRAVITY (GRAVITY Collaboration 2017) medium-resolution (R ~ 500) K-band spectra of HD 206893 A and B using the Unit Telescopes (UTs) at the VLTI. The observing log is presented in Table 1. The atmospheric conditions varied between average (atmospheric coherence time τ₀ = 5 ms) and below average (τ₀ ≈ 1.5 ms). Both epochs were obtained as bad weather backup targets for the large active galactic nucleus program (PI E. Sturm) based on a time exchange agreement.

From the raw GRAVITY data, we extracted the coherent flux following the standard recipe of the official ESO data reduction pipeline (Lapeyrere et al. 2014). During this first step, the pipeline computes the coherent flux observed on the host star and the companion. However, the coherent flux observed on the faint companion is still contaminated by the halo of the bright host star. This contamination was removed using a Python package developed by our team². The individual steps of this package are outlined in Appendix A of GRAVITY Collaboration (2020) and its output is the decontaminated ratio of the coherent flux between the companion and the host star.

The astrometry for each epoch of data was obtained from the phase of the ratio of the coherent fluxes and is presented in Table 2. The uncertainties were estimated from the scatter of the astrometric values obtained independently for each individual exposure. The typical precision is ~ 100 μas, which is much worse than the theoretical limit of 16.5 μas determined by Lacour et al. (2014). This can be attributed to the low and high frequency phase errors present in our data and introduced by instrumental aberrations (GRAVITY Collaboration 2021). Furthermore, due to the asymmetry of the uv-plane, we used the correlation coefficient ρ to properly describe the confidence intervals and the correlations between the right ascension and the declination offset.

Finally, a spectrum of the companion for each epoch of data was obtained from the amplitude of the ratio of the coherent fluxes. The host star is essentially unresolved in our observations such that we did not need to correct this ratio for the visibility of the host star. The ratio of the coherent fluxes was multiplied by a model spectrum of the host star, which was obtained by interpolating the BT-NextGen stellar model grid (Allard et al. 2012) for the stellar parameters presented in Table 3. This yields a spectrum of HD 206893 B for each epoch of data. The model spectrum of the host star was scaled to match the stellar photometry presented in Table 3, yielding a flux calibrated radius of 1.36 R_⊙. The 2MASS and WISE photometry are sufficient to constrain the stellar spectrum over the GRAVITY wavelength range (cf. Fig. 2). Here, the stellar spectrum is essentially smooth and the weak spectral lines are negligible for our analysis.

Fig. 2 BT-NextGen model spectrum of HD 206893 A (black line), scaled to fit the shown 2MASS and WISE photometry (gray data points). Synthetic photometry based on the model spectrum is also shown with black data points. The top panel shows the transmission curves corresponding to each photometric data point and the bottom panel shows the residuals between the synthetic and the observed photometry.

3 Orbit analysis

3.1 Astrometric analysis

From the interferometric observations with GRAVITY we obtain two new astrometric data points with an unprecedented precision of ~ 100 μas for HD 206893 B (cf. Table 2). Together with astrometric data from the literature (cf. Table A.1), we estimate the orbital parameters of HD 206893 B with orbitize!³ (Blunt et al. 2020), which infers the posterior distribution of the orbital parameters through Markov Chain Monte Carlo (MCMC) sampling with ptemcee⁴ (Foreman-Mackey et al. 2013; Vousden et al. 2016). We initialize the sampler with 20 temperatures, 500 walkers, and 50 000 steps per walker. By visual inspection of the walker chains, we assess convergence and reject the first 40 000 steps before computing the posterior distribution from each walker at the lowest temperature. We note that we use the priors presented in Table 4 for the orbital parameters. These priors are chosen very conservatively in order to not constrain the posterior. The posterior distribution of the orbital parameters and the inferred orbital solutions together with the NACO, SPHERE, GPI, and GRAVITY astrometry are shown in Fig. 3 and the posterior values are quoted in Table 4 (scenario 1). We restrict our orbit analysis to direct imaging astrometry of HD 206893 B because Grandjean et al. (2019) noted that the radial velocities and the Gaia astrometry can only be explained when including an inner companion (HD 206893 C) in the fit. However, including this inner companion could introduce biases in our fit since its orbital parameters and dynamical mass remain poorly constrained.

In general, the GRAVITY astrometry is consistent with those from NACO, SPHERE, and GPI and our orbital solutions are consistent with those from Grandjean et al. (2019) and Ward-Duong et al. (2021). Both of them obtained a double-peaked semimajor axis distribution and an anticorrelation between semimajor axis a and eccentricity e. However, adding the GRAVITY astrometry disfavors small eccentricities e ~ 0 and removes the second peak in the semimajor axis distribution. Instead, a smaller semimajor axis of ${9.28}_{-0.93}^{+1.77} \text{au}$ and a higher eccentricity of ${0.29}_{-0.11}^{+0.06}$ are preferred. Nevertheless, the inclination i is similar to that obtained by Grandjean et al. (2019) and Ward-Duong et al. (2021) with a maximum likelihood value of ~ 145 deg. The mutual inclination i_m between the debris disk of the system reported by Milli et al. (2017) and Marino et al. (2020) and HD 206893 B is $$i_{\text{m}}=\arccos \left(\cos (i)\cos (i_{\text{d}})+\sin (i)\sin (i_{\text{d}})\cos (\Omega -{\Omega }_{\text{d}})\right),$$(1)

where i_d = 140 ± 3 deg is the inclination of the debris disk and Ω_d = 61 ± 4 deg is the longitude of the ascending node of the debris disk (Marino et al. 2020). We compute the mutual inclination and its uncertainty by drawing samples from the posterior distribution of the orbit fit and a normal distribution for i_d and Ω_d. We find $i_{\text{m}}={20.8}_{-11.2}^{+13.6} \text{deg}$, which means that the debris disk and the companion are roughly aligned. We note that we use priors between 0 and 180 deg for the longitude of the ascending node Ω to enforce the debris disk and the companion orbiting in the same direction. However, the direction of rotation of the debris disk is unconstrained and there is a 180 deg ambiguity in the mutual inclination i_m between the debris disk and the companion (Heintz 1978). Therefore, they might as well orbit in opposite directions.

Interestingly, Marino et al. (2020) found that if they enforce coplanarity between the debris disk and the companion, this would lead to the other one of the degenerate solutions being preferred, namely a larger semimajor axis of ~ 11.4 au and a smaller eccentricity of ~ 0.14. To verify their findings, we ran another fit with Gaussian priors of 140 ± 3 deg for the inclination i and 61 ± 4 deg for the longitude of the ascending node Ω. The results are shown in Fig. A.1 and Table 4 (scenario 2) and confirm the findings of Marino et al. (2020), even when adding the GRAVITY astrometry. Therefore, depending on whether coplanarity with the debris disk is assumed or not, degenerate orbital solutions are obtained for HD 206893 B.

3.2 Gaia proper motion anomaly

An independent constraint on the mass of HD 206893 B can be obtained from the proper motion anomaly of its host star (HD 206893 A) measured by Gaia (Kervella et al. 2020). While the proper motion measured by Gaia (Gaia Collaboration 2018) traces the photocenter of the system, the long-term proper motion derived from the difference between the Gaia and HIPPARCOS (van Leeuwen 2007) positions traces the barycenter of the system, once corrected for the influence of the companion at the time of the HIPPARCOS and Gaia measurements. The difference between the proper motion measured by Gaia and the long-term proper motion is the proper motion anomaly that must be caused by one or multiple faint companions.

Grandjean et al. (2019) have already found that the proper motion anomaly measured between the Gaia DR2 and the HIPPARCOS data cannot be explained by HD 206893 B alone and suggested the presence of another closer-in companion of ~15 M_Jup (HD 206893 C). Here, we consider the proper motion anomaly measured between the Gaia EDR3 and the HIPPARCOS data (PMa_RA = −0.102 ± 0.037 and PMa_DEC = −0.612 ± 0.028) to obtain independent constraints on the masses of HD 206893 B and C and the on-sky position of HD 206893 C. Since the orbit of HD 206893 B is known from direct observations (cf. Sect. 3.1) we can compute the proper motion anomaly it causes on its host star. Consistently with Grandjean et al. (2019), we find that the proper motion anomaly measured by Gaia EDR3 cannot be explained by HD 206893 B alone because a proper motion component tangential to the one caused by HD 206893 B is necessary to fit the data. Hence, we assume another companion (HD 206893 C) on a similarly oriented orbit (same inclination and longitude of the ascending node as HD 206893 B), but with zero eccentricity and smaller semimajor axis. We compute the proper motion anomaly it causes on its host star as a function of its on-sky position and its mass, marginalizing over 100 randomly drawn samples of the posterior distribution of the orbit of HD 206893 B and the uncertainties in the proper motion anomaly measured by Gaia EDR3.

Figure 4 shows the predicted masses for HD 206893 B and C. Regions where either of the companions would have negative mass, where the orbital period of HD 206893 C would be smaller than half of the Gaia EDR3 measurement timespan, or where the apparent separation of HD 206893 C would be more than 150 mas are ignored. Our orbital fits for HD 206893 B suggest a minimum orbital separation of ~150 mas and we assume that the orbits of the two companions cannot cross for dynamical stability reasons. For most of the on-sky positions, the mass of HD 206893 C should be between ~ 8–15 M_Jup, a result that is largely consistent with the prediction of Grandjean et al. (2019). Furthermore, depending on the mass of HD 206893 B, the on-sky position of HD 206893 C can be strongly constrained.

Fig. 3 Posterior distribution of the orbital parameters (top) and orbital solutions together with the NACO, SPHERE, GPI, and GRAVITY astrometry (bottom) of HD 206893 B. In the top panel, the values state the 68% confidence intervals around the median. In the bottom panel, the black star highlights the position of HD 206893 A and all error bars show the 1–σ confidence intervals.

Fig. 4 Masses of HD 206893 B (left) and C (right) as a function of the on-sky position of HD 206893 C (at the reference epoch of Gaia EDR3, 2016.0) derived from the proper motion anomaly measured between the Gaia EDR3 proper motion and the HIPPARCOS–Gaia EDR3 long-term proper motion. Physically implausible regions due to negative masses or dynamical instability have been ignored.

Fig. 5 Combined GRAVITY K-band spectrum of HD 206893 B together with the SPHERE Y –H-band spectrum from Delorme et al. (2017) and the GPI J, H, K1, and K2-band spectra fromWard-Duong et al. (2021). The shaded regions highlight the 1–σ confidence intervals. For reference, absorption bands of water and carbon monoxide are indicated.

4 Spectral analysis

Apart from the two astrometric data points, we also obtain two K-band spectra at a resolution of R ~ 500 for HD 206893 B, one for each epoch of GRAVITY data. These two spectra are consistent with each other and we combine them into a single final spectrum, accounting for the covariances, shown in Fig. 5 together with the SPHERE spectrum from Delorme et al. (2017) and the GPI spectra from Ward-Duong et al. (2021). In Sects. 4.1 and 4.2, we use these spectra for atmospheric model fitting and spectral retrievals of HD 206893 B, respectively. To constrain the fits between 3.5 and 5 μm, we supplement the spectra with photometry of HD 206893 B from the literature (cf. Table B.1).

4.1 Atmospheric model fitting

By combining the GRAVITY spectrum with SPHERE and GPI spectra and photometry available in the literature, we reach a broad spectral coverage from ~ 1–5 μm. This spectral region contains absorption bands of water, carbon-monoxide, and methane and is broad enough to estimate the effective temperature, the radius, and the surface gravity of an object. We estimate these parameters for HD 206893 B by fitting its spectra and photometry with atmospheric model grids using species⁵ (Stolker et al. 2020). There is a variety of atmospheric model grids for giant planets and brown dwarfs, all of them being slightly different in terms of underlying physics and complexity. However, all of them assume radiative-convective equilibrium to calculate the temperature structure of the atmosphere self-consistently. Here, we use three different grids: the BT-Settl-CIFIST grid (Allard et al. 2012), the DRIFT-PHOENIX grid (Helling et al. 2008b) (which includes metallicity as an additional free parameter), and the Exo-REM grid (Baudino et al. 2015; Charnay et al. 2018) (which includes both metallicity and C/O ratio as additional free parameters). All three grids include photospheric absorption by dust clouds, but with different approaches to calculate the cloud densities, grain size distributions and compositions. We bin the grids to the spectral resolution of the respective instrument, and use the spectra and filter curves to calculate the synthetic photometry and filter-weighted average flux, respectively.For all grid parameters, we use uniform priors whose boundaries are presented in Table 5. Our atmospheric model fits account for the covariances in the GRAVITY spectrum according to Sect. 2 of Greco & Brandt (2016). While fits to the GRAVITY and the SPHERE spectra alone show good photometric agreement between the two, there seems to be a significant offset between the GPI H-band and the SPHERE spectrum, which may indicate a systematic error in the absolute flux calibration. Given that the SPHERE spectrum agrees well with the GPI J-band spectrum, we decided to fit a separate scaling parameter to each of the GPI spectra while keeping the GRAVITY and the SPHERE spectra fixed in order to preserve the extremely red color of HD 206893 B. Then, we infer the posterior distribution of the model parameters with nested sampling using PyMultiNest⁶ (Buchner et al. 2014; Feroz et al. 2009, 2019).

Several hypotheses have been put forward to explain the extreme redness of HD 206893 B, most notably extinction by local dust, either extra-photospheric or in the form of a circumplanetary disk by Milli et al. (2017), Delorme et al. (2017), and Ward-Duong et al. (2021). Other possibilities like reddening by interstellar dust or extinction by the debris disk could be mostly ruled out. Ward-Duong et al. (2021) could not find any significant interstellar extinction toward the host star based on its photometry, and we can confirm this finding by visual inspection of stellar Ca-lines in high-resolution spectra of HD 206893 A (A.-M. Lagrange, priv. comm.). It seems unlikely that there is an interstellar dust cloud that is obscuring HD 206893 B but not its host star, which is separated by only ~ 250 mas. Moreover, the debris disk of the system would need to be unrealistically optically thick (τ ~ 1.7) to explain the observed reddening, even if viewed edge-on (Ward-Duong et al. 2021). Therefore, extinction by local dust is the most plausible explanation for the extremely red color of HD 206893 B, and we include additional extinction by high-altitude dust clouds made of crystalline enstatite (MgSiO₃) grains in our atmospheric model fits. Since other dust species such as forsterite, corundum, and iron predict similar extinction curves for grain sizes between 0.1–1 μm (e.g., Ward-Duong et al. 2021), we only consider enstatite grains for simplicity here. These grains are described by a log-normal size distribution $$\frac{\text{d}n}{\text{d}a}=\frac{N}{a\sqrt{2\pi }\ln {\sigma }_{\text{a}}}\exp \left(-\frac{{\ln }^2(a/a_{\text{mean}})}{2{\ln }^2{\sigma }_{\text{a}}}\right),$$(2)

where n is the number concentration of grains smaller than a, N is the total number concentration of grains, a is the grain size, a_mean is the geometric mean grain size, and σ_a is the grain size geometric standard deviation (which is dimensionless, Ackerman & Marley 2001). A log-uniform prior between 0.1–10 μm is used for a_mean and a uniform prior between 1.1–10 is used for σ_a. Smaller grains would grow by condensation within timescales of less than a second and are therefore not considered (Charnay et al. 2018). Then, we compute the extinction cross-section σ_ext,λ(a_mean, σ_a) using PyMieScatt⁷ (Sumlin et al. 2018) and scale the flux of the default spectra F_default to that of the reddened spectra $$F_{\text{red}}=F_{\text{default}}\times {10}^{\left(-\frac{A_V}{2.5}\frac{{\sigma }_{\text{ext},\lambda }}{{\sigma }_{\text{ext},\text{V}}}\right)},$$(3)

where A_V is the extinction in the Bessel V -band, another free parameter with a uniform prior between 0–5 mag, and σ_ext,V is the extinction cross-section averaged over the Bessel V -band. In total, our enstatite dust model has three free parameters (a_mean, σ_a, and A_V), which are inferred along with the parameters of the atmospheric model grids.

Table 6 summarizes the atmospheric parameters obtained for HD 206893 B based on the three different atmospheric model grids without (“plain”) and with (“dusty”) additional extinction. The inferred effective temperatures are very similar to those obtained by Delorme et al. (2017) and Ward-Duong et al. (2021), but the surface gravities confirm the trend observed by Ward-Duong et al. (2021), namely that the H and K-band spectra prefer lower surface gravities than those obtained by Delorme et al. (2017) for the SPHERE spectrum at shorter wavelengths. Overall, the parameters inferred from the plain models are spread over a wider range of parameter space than those inferred from the dusty models. Moreover, all dusty models fit the observed data better than the plain models since they have smaller ${\chi }_{\text{red}}^2$. This is not completely surprising given that the dusty models have three more free parameters for describing the additional extinction than the plain models. Most noticeably, for both the plain and the dusty models the DRIFT-PHOENIX (DP) grid predicts a significantly higher surface gravity and mass for HD 206893 B than the BT-Settl-CIFIST (BT) and Exo-REM (ER) grids. However, while the DP grid yields the best fit (i.e., the smallest ${\chi }_{\text{red}}^2$) for the plain models, it yields the worst fit (i.e., the highest ${\chi }_{\text{red}}^2$) for the dusty models.

Another striking difference between the DP grid and the BT and ER grids are the extinction parameters that they predict for the dusty models. While the BT and ER grids consistently prefer small grains with a geometric mean size of ~ 0.33–0.34 μm and a geometric standard deviation of ~ 1.30–1.34, the DP grid prefers large grains with a geometric mean size of ~ 2.29 μm and a geometric standard deviation of ~1.17. This is a difference in geometric mean grain size of almost an order of magnitude. We note that the DP grid uses a different cloud model than the BT and ER grids. With DP, gas is mixed to high altitudes where dust then forms and grows as it rains down. With BT and ER, the cloud model from Ackerman & Marley (2001) is used which assumes that the cloud particles are mixed from the cloud base upward. These fundamentally different approaches could cause the observed difference in predicted dust grain size. Helling et al. (2008a) have further found that the dust to gas ratio predicted by the DRIFT model is larger than the one predicted by the Settl model at small pressures, where the mean grain size is below 1 μm. However, with DP the difference in ${\chi }_{\text{red}}^2$ with and without extinction is small, that is, DP with large grains (dusty) does not fit the data much better than DP without large grains (plain).

Figure 6 shows the best fit model spectra for the dusty models in dark red together with the NACO, SPHERE, GPI, and GRAVITY spectra and photometry of HD 206893 B. It is noteworthy that there are significant differences between the BT and ER grids and the DP grid regarding the depth of the 1.4 μm water absorption feature and the morphology of the H- and K-band peaks. The triangular shaped H-band peak observed by GPI and fit well by the BT and ER grids (but not the DP grids) is typical for a young low-gravity object (e.g., Kirkpatrick et al. 2012; Allers & Liu 2013). Moreover, they deviate significantly at longer wavelengths (> 2.5 μm). There, the available NACO photometry is not precise enough to set meaningful constraints on the model parameters.

The exact same three dusty models are shown in light red, but before applying the additional extinction. Here, the difference in predicted grain size between the DP grid and the BT and ER grids becomes very clear. While for the BT and ER grids, the difference between unextinct (light red) and extinct (dark red) model spectrum decreases with increasing wavelength and approaches zero over the L and M-band, the extinction reaches its maximum near the L-band (where the wavelength is on the order of the geometric mean grain size) for the DP grid. Finally, we note that the ER grid predicts a significantly higher V -band extinction than the BT grid, despite the very similar grain size parameters. This is the case because the ER grid converges toward a significantly larger radius if compared to the BT grid, resulting in a higher bolometric luminosity and therefore requiring higher extinction.

It is also noteworthy that in the absence of additional dust grains (i.e., the plain models), the metallicity is found to hit the upper bound of the DP and ER grids (i.e., significantly super-solar metallicity). It has been observed before that a high metallicity facilitates the formation of dust grains and might be an explanation for the unusually red L dwarf population (Looper et al. 2008; Stephens et al. 2009; Gizis et al. 2012; Marocco et al. 2014). A slight anticorrelation between metallicity and dust extinction (A_V) for the dusty ER grid (cf. Fig. C.6) supports the finding that higher metallicity causes a redder color. However, for the dusty DP grid, where the metallicity also hits the upper bound, we cannot identify such a correlation.

4.2 Spectral retrieval

The atmospheric model fits in Sect. 4.1 have shown that HD 206893 B is difficult to explain with currently available grids and often drives parameters such as the surface gravity or the metallicity to the grid boundaries in order to mimic an extremely red color. Spectral retrievals are better suited for exploring a wide range of parameters and multi-species gas opacities (e.g., Ward-Duong et al. 2021). We thus perform spectral retrievals with petitRADTRANS (Mollière et al. 2019, 2020) and ATMO (Tremblin et al. 2015, 2016) on the same spectra and photometry of HD 206893 B that we also used for the atmospheric model fits in Sect. 4.1.

Fig. 6 Atmospheric models fitted to the observed spectra and photometry of HD 206893 B shown in the background. Dark red lines show the best fit dusty models that include additional extinction by high-altitude dust clouds made of enstatite grains and light red lines show the exact same models before applying the extinction (i.e., without dust). The GPI spectra are not to be taken to face value since they are rescaled during each of the fits. BT = BT-Settl-CIFIST, DP = DRIFT-PHOENIX, and ER = Exo-REM.

4.2.1 petitRADTRANS

For the spectral retrieval with petitRADTRANS (pRT), we follow mostly the implementation including scattering clouds as described in Mollière et al. (2020) for the case of HR 8799 e. We briefly summarize that the P–T structure is parameterized into three different regions. Specifically, we use free temperature nodes at high altitudes, the Eddington approximation for the photospheric region, and a moist adiabat for the deep, radiative-convective part of the atmosphere. Gas abundances are assumedto be in chemical equilibrium, but with an additional parameter for a quenching pressure above which the CO, CH₄, and H₂O abundances are fixed. We include CO, H₂O, CH₄, NH₃, CO₂, H₂S, Na, K, PH₃, VO, TiO, FeH as molecular and atomic line species plus Rayleigh scattering and collision induced absorption (CIA) of H₂ and He. Given the similar temperature to HR 8799 e, we only consider cloud condensates composed of MgSiO₃ and Fe.

Since HD 206893 B is thought to have an unusually cloudy atmosphere, we use the cloud optical depth in the photospheric region as free parameter instead of the mass fractions themselves. The cloud mass fraction above the cloud base, X_cloud, is parameterized as (Mollière et al. 2020) $$X_{\text{cloud}}(P)=X_{\text{cloud},0}{\left(\frac{P}{P_{\text{base}}}\right)}^{f_{\text{sed}}},$$(4)

where X_cloud,0 is the cloud mass fraction at the cloud base, P_base is the pressure at the cloud base, and f_sed is the settling parameter (assumed to be the same for MgSiO₃ and Fe). The cloud parameter, τ_cloud, is then defined as the optical depth at the τ = 1 pressure of the gas-only atmosphere (i.e., neglecting the clouds). To ensure a quasi-physical solution for the cloud properties, we reject samples for which $$X_{\text{cloud},0}>2X_{\text{eq}}(1+f_{\text{sed}}),$$(5)

where X_eq is the equilibrium mass fraction at the cloud base as calculated from the elemental abundances and [Fe/H]. This expression stems from requiring that the surface density of the cloud does not exceed the surface density of condensing species in the atmosphericcolumn above the cloud base, where we allowed for an additional factor of two to slightly relax this condition.

Similar to the atmospheric model fits in Sect. 4.1, we use a nested sampling algorithm (Feroz et al. 2009; Buchner et al. 2014) to sample the posterior distributions of the 20 free parameters with 4000 live points. The model spectra are smoothed to the respective instrument resolution before they are resampled to the data wavelengths, and a separate scaling factor is fitted to each of the GPI spectra to account for systematic calibration offsets. The results of the retrieval are presented in Table 7 and the retrieved spectrum and P–T profile are shown in Figs. 7 and 8, respectively. We note that, for clarity, the scaling factors have not been applied when plotting the spectra, to be able to show the best fit models of pRT and ATMO (described below), the latter of which did not retrieve scaling factors.

Fig. 7 Retrieved spectra of HD 206893 B with the observed spectra and photometry shown in the background. The GPI spectra are not to be taken to face value since they are rescaled during the retrieval with petitRADTRANS (pRT).

Fig. 8 Pressure-temperature profiles of the atmosphere of HD 206893 B retrieved with petitRADTRANS (pRT) and ATMO. The dashed lines show the 1–σ confidence intervals. With pRT we obtain a bimodal posterior for the C/O ratio and the quenching pressure, and the P–T profile for each of the two solutions is shown separately (thick blue line C/O > 0.7, thin blue line C/O < 0.7). All models probe approximately 0.01 to 0.1 bar pressure levels. Dotted lines indicate the condensation curves of various species calculated at 10× solar abundance (Visscher et al. 2010).

4.2.2 ATMO

Another spectral retrieval was also performed with ATMO, a 1D–2D radiative transfer model for planetary atmospheres (Tremblin et al. 2015, 2016). More comprehensive descriptions of the model can be found in Amundsen et al. (2014), Drummond et al. (2016), Goyal et al. (2018), and Phillips et al. (2020). The retrieval aspects of ATMO have been used to fit transit and secondary eclipse data before (e.g., Evans et al. 2017; Wakeford et al. 2017), and are applied here on the HD 206893 B data.

We fit the data assuming chemical equilibrium, and include rainout (Goyal et al. 2020) to account for the depletion of gas phase species due to condensation. The total opacity of the gas mixture is computed using the correlated-k approximation using the random overlap method with resorting and rebinning (Amundsen et al. 2014). The k-coefficients and chemical equilibrium are calculated “on the fly” for each atmospheric layer, spectral band, and iteration such that the derived opacities are physically self-consistent with the P–T profile and chemical composition for each model evaluation in the retrieval. We include spectrally active species H₂ -H₂ and H₂-He CIA opacities, as well as H₂O, CO₂, CO, CH₄, NH₃, K, TiO, VO, FeH, CrH, HCN, C₂H₂, H₂S, and H- (see Goyal et al. 2018 for further details). We fit for the elemental abundances of C and O separately, with the rest of the elements described by a single metallicity parameter, [Fe∕H].

Scattering and absorption effects from clouds on the spectrum were parameterized as follows. We include scattering from small particles using an enhanced Rayleigh-like scattering (Lecavelier Des Etangs et al. 2008) opacity, parameterized as $$\sigma {(\lambda )}_{\text{haze}}={\delta }_{\text{haze}}{\sigma }_{\text{0}}{(\lambda /{\lambda }_0)}^{{\alpha }_{\text{haze}}},$$(6)

where σ(λ) is the total scattering cross-section of the material, δ_haze is an empirical enhancement factor, σ₀ is the scattering cross section of molecular hydrogen at 0.35 μm, and α_haze is a factor determining the wavelength dependence. Condensate cloud absorption is fit separately, and is assumed to have a gray wavelength dependence calculated as $$\kappa {(\lambda )}_{\text{cloud}}={\delta }_{\text{cloud}}{\kappa }_{{\text{H}}_{\text{2}}},$$(7)

where $\kappa {(\lambda )}_{\text{cloud}}$ is the “cloud” absorption opacity, δ_cloud is an empirical factor governing the strength of the gray scattering, and ${\kappa }_{{\text{H}}_{\text{2}}}$ is the scattering opacity due to H₂ at 0.35 μm. $\sigma {(\lambda )}_{\text{haze}}$ and $\kappa {(\lambda )}_{\text{cloud}}$ are added to the total gaseous scattering and absorption, respectively, with a further parameter specifying the pressure level at the top of the gray cloud. To parameterize the P–T profile, we use an analytic radiative equilibrium model by Guillot (2010). Again, we use nested sampling to sample the posterior distribution (Feroz et al. 2009), fitting for a total of 12 free parameters. A flux scaling for calibration offsets was not applied, however. The results of the retrieval are presented in Table 7 and the retrieved spectrum and P–T profile are shown in Figs. 7 and 8, respectively.

Fig. 9 Parameters inferred for HD 206893 B from atmospheric model fits and spectral retrievals compared to AMES-Cond evolutionary tracks. Light red points show the best fit parameters for the plain models and dark red points show the best fit parameters for the dusty models including additional extinction by high-altitude dust clouds made of enstatite grains. The evolutionary tracks are shown for objects with exactly those masses printed below the colorbar. Curves of constant age are shown in dashed gray. BT = BT-Settl-CIFIST, DP = DRIFT-PHOENIX, ER = Exo-REM, and pRT = petitRADTRANS.

4.3 Evolutionary tracks

Similar to Delorme et al. (2017), we compare our best fit atmospheric parameters to evolutionary tracks for giant planets and brown dwarfs to ensure that they correspond to physically plausible substellar companions. Therefore, we use the AMES-Cond evolutionarytracks (Baraffe et al. 2003). Figure 9 shows these tracks for objects of different masses from 1–100 M_Jup in the radius versus effective temperature and surface gravity versus effective temperature planes. The best fit atmosphericparameters are overplotted in light red (plain models), dark red (dusty models), and orange red (free retrievals).

The most significant outliers are the plain ER grid and the ATMO retrieval, which predict unexpectedly large radii of ${2.32}_{-0.02}^{+0.02} R_{\text{Jup}}$ and ${2.23}_{-0.09}^{+0.09} R_{\text{Jup}}$ for a relatively cool (${1049}_{-4}^{+2} \text{K}$ and ${1113}_{-52}^{+51} \text{K}$) and low-mass (${6.58}_{-0.10}^{+0.11} M_{\text{Jup}}$ and ${1.07}_{-0.24}^{+0.82} M_{\text{Jup}}$) object, respectively. With the additional extinction, the radius of the dusty ER grid decreases and the effective temperature increases, leading to an object that is roughly consistent between the atmospheric models and the evolutionary tracks with an extremely young (< 10 Myr) planetary-mass (< 5 M_Jup) companion, given the uncertainties on the surface gravity and the radius from the atmospheric model fits. In the same parameter range, the pRT retrieval can be found, but with a slightly smaller radius and an even lower surface gravity. We note, however, that the dusty ER grid converges toward the lower boundary of the surface gravity (~ 3.5). Other inconsistent parameters are predicted by both the plain BT and dusty DP grids. While the low surface gravity predicted by the plain BT grid is consistent with an extremely young object (< 3 Myr), its predicted small radius is consistent with a rather old object (>300 Myr). The same is observed the other way around for the dusty DP grid. Here, the predicted high surface gravity is consistent with a rather old object and the predicted large radius is consistent with an extremely young object. However, the dusty BT and plain DP grids predict objects that are roughly consistent between the atmospheric parameters and the evolutionary tracks. The dusty BT grid suggests a moderately young (~ 3–300 Myr) object somewhere between ~ 5–30 M_Jup. The plain DP grid suggests a rather old (~ 100–1000 Myr) object somewhere between ~ 15–75 M_Jup. We note that such objects are also in agreement with the age and mass predicted for HD 206893 B by Delorme et al. (2017).

Overall, we find that the BT and ER grids require additional extinction in order to predict physically plausible objects. This is different for the DP grid, which becomes unphysical when additional extinction is included. This could also be related to the much larger grain sizes predicted by the DP grid if compared to the BT and ER grids, which are not expected for high-altitude dust clouds (Hiranaka et al. 2016). Moreover, we find that the three different atmospheric model grids predict objects that populate a wide range of parameter space in age and mass. The (dusty) ER grid predicts an extremely young (< 10 Myr) object of < 5 M_Jup, the (dusty) BT grid predicts a moderately young (~ 3–300 Myr) object of ~ 5–30 M_Jup, and the (plain) DP grid predicts a rather old (100–1000 Myr) object of ~ 15–75 M_Jup.

4.4 Color-magnitude diagram

In a color-magnitude diagram, it can easily be seen that HD 206893 B is the reddest known substellar object (Milli et al. 2017; Delorme et al. 2017; Ward-Duong et al. 2021). Figure 10 shows J–K and H–K color-magnitude diagrams of HD 206893 B and other known planetary-mass companions. For reference, M, L, and T-dwarfs from the SpeX Prism Spectral Libraries⁸ are shown in the background. The apparent magnitudes of the planetary-mass companions were taken from Delorme et al. (2017) (HD 206893 B), Currie et al. (2013) (β Pic b), Skemer et al. (2016) (GJ 504 b), Rajan et al. (2017) (51 Eri b), Zurlo et al. (2016) (HR 8799 b, c, d, and e), Patience et al. (2012) (2M1207 b), and Bohn et al. (2020) (TYC 8998 b). Apparent magnitudes were converted to absolute magnitudes using distances (i.e., parallaxes) from SIMBAD⁹.

The light red lines show the reddening vectors of our best fit enstatite dust model for different V -band extinctions. Here, we assume a_mean = 0.33 μm and σ_a = 1.30, consistent with the best fit dusty BT and ER models. We did not plot the reddening vector for the best fit dusty DP model because it corresponds to a physically implausible object (cf. Sect. 4.3). The shown reddening vectors extend from A_V = 0–3. For the best fit dusty BT and ER models, the predicted V -band extinction of A_V ~ 2.0 and A_V ~ 2.9, respectively,brings HD 206893 B back to the red end of the substellar main sequence, close to where the planetary-mass object β Pic b is located. This implies that HD 206893 B could indeed be a very dusty companion around a young moving group member, such as β Pic b, and marks another similarity between the HD 206893 and the β Pic system besides the very similar system architecture. Finally, compared to the interstellar reddening law applied by Ward-Duong et al. (2021), our enstatite dust model predicts a similar reddening slope in the M_H versus H–K color-magnitude diagram while requiring smaller V -band extinction values of ~ 2–3 instead of ~10 in order to bring HD 206893 B back to the red end of the substellar main sequence (cf. their Fig. 5). The existence of the CT Cha companion (Schmidt et al. 2008) mentioned by Ward-Duong et al. (2021) and suffering from an extreme V -band extinction of ~ 5 magnitudes therefore puts the values of A_V ~ 2–3 obtained for HD 206893 B in a realistic regime.

Fig. 10 J–K (left) and H–K (right) color-magnitude diagrams showing HD 206893 B in red together with the reddening vector of our best fit extra-photospheric enstatite dust model for different V -band extinctions in light red. Other known planetary-mass objects are shown in orange and M, L, and T-dwarfs from the SpeX Prism Spectral Libraries are shown in green. For the dust model, we assume amean = 0.33 μm, σa = 1.30, and AV = 0–3. The red points are in steps of 1 mag.

5 Discussion

5.1 Astrometry

In Sect. 3.1, we infer the orbital parameters of HD 206893 B for two different scenarios. Scenario 1 is constrained by the data only and scenario 2 has an additional constraint on coplanarity with the debris disk of the system observed by Milli et al. (2017) and Marino et al. (2020). By comparison with earlier works from Grandjean et al. (2019) and Ward-Duong et al. (2021), we find that the GRAVITY data resolves the degeneracy between a lower eccentricity, larger semimajor axis and a higher eccentricity, smaller semimajor axis orbit by preferring the latter of these in scenario 1. This is interesting because this orbital solution for HD 206893 B is only roughly aligned with respect to the debris disk of the system (i_m < 34.4 deg). Marino et al. (2020) mention that a misalignment between the orbit of HD 206893 B and the debris disk of the system should be unlikely given its age of at least 50 Myr (Delorme et al. 2017). They argue that HD 206893 B should align with the debris disk of the system due to secular interactions on timescales of only ~10 Myr. In scenario 2, where we enforce alignment between the orbit of HD 206893 B and the debris disk of the system, we clearly find that the lower eccentricity, larger semimajor axis orbit is preferred. This is in agreement with Marino et al. (2020), who obtained the same result without the additional GRAVITY data, and suggests that the GRAVITY data alone prefers a slight misalignment between the orbit of HD 206893 B and the debris disk of the system.

If this misalignment is confirmed by future GRAVITY observations, an explanation for it needs to be found. One possibility would be a significantly younger age (<50 Myr) for the system. We note that such a young age would be in agreement with the age constraint set by comparing our best fit dusty BT and ER models and the spectral retrievals with evolutionary tracks (cf. Sect. 4.3). In addition, an age of ~ 40 Myr would be expected according to Torres et al. (2008) if the system was part of the Argus moving group, for which the probability is ~ 61% (Ward-Duong et al. 2021). We note, however, that the analysis of lithium and barium abundances points to an older age for the host star (Delorme et al. 2017). The discrepancy between the predicted Argus moving group membership from Banyan Sigma (~ 61%, Ward-Duong et al. 2021) and Banyan II (~13.5%, Delorme et al. 2017)can be attributed to the revised kinematics of the Argus association put forth by Zuckerman (2019), though, and another recent and distinct Bayesian methodology by Lee & Song (2019) also provides a membership probability of ~ 63%. Another possibility would be dynamical interactions between HD 206893 B and the other putative companions predicted toexist in this system (e.g., Wu & Lithwick 2011). There is evidence for a second, massive (~ 15 M_Jup) and close (1.4–2.6 au) companion interior to the orbit of HD 206893 B from radial velocity and Gaia data (Grandjean et al. 2019) and a third 0.4–1.7 M_Jup companion further out responsible for carving the gap at ~ 74 au in the debris disk observed with the Atacama Large Millimeter/Submillimeter Array (Marino et al. 2020). However, a profound analysis of potential planet-planet interactions is beyond the scope of this work. Finally, we note that an eccentricity of ~ 0.3 for HD 206893 B preferred by the GRAVITY data is in better agreement with the eccentricity distribution of the brown dwarf population than that of wide-separation (5–100 au) giant planets, the latter of which show a preference for low eccentricities (≈ 0.05–0.25, Bowler et al. 2020). Still, sorting an individual companion into one or the other class of objects based on its eccentricity remains highly speculative, and as highlighted before, there might be other mechanisms responsible for an excitement of HD 206893 B’s eccentricity, such as dynamical interactions with other companions in the system.

5.2 Grid retrievals

In Sect. 4, we analyze the near-infrared spectrum of HD 206893 B. We obtain a dense spectral coverage between 1–2.5 μm by combining data from SPHERE, GPI, and GRAVITY that we extend up to ~ 5 μm with VLT/NACO photometry from Stolker et al. (2020). While Delorme et al. (2017) and Ward-Duong et al. (2021) found that currently available atmospheric models for giant planets and brown dwarfs fail to predict the extremely red color and the very shallow 1.4 μm water absorption feature of HD 206893 B at the same time, they also compared the spectrum of HD 206893 B to those of other dusty, low-gravity or young M and L-dwarfs and found that none of these objects could reproduce all spectral features of HD 206893 B. Hence, both of these authors tried to reconcile the spectrum of HD 206893 B with those of other dusty, low-gravity or young M and L-dwarfs by including additional extinction by high-altitude dust clouds. Therefore, they explored a variety of dust species (forsterite, enstatite, corundum, and iron), grain sizes (0.05–1 μm), and extinction values. They found that reddening by forsterite or enstatite grains with sizes between 0.27–0.50 μm yields the best match to a very low-gravity L3-dwarf. Such a reddening by small (< 1 μm) dust grains in the cool upper atmosphere has been suggested before by Marocco et al. (2014) and Hiranaka et al. (2016) to match the spectra of unusually red L-dwarfs with those of spectroscopic standards. They concluded that scattering by dust clouds with grain sizes between 1–100 μm included in current atmospheric models is not sufficient to describe the peculiarly red L-dwarfs, which seem to require additional scattering by smaller grains in the cool upper atmosphere.

Here, we focused on a slightly different approach by adding extinction by high-altitude dust clouds directly to the atmospheric model grids. Since Ward-Duong et al. (2021) mention that all considered dust species predict similar extinction curves for grain sizes between 0.1–1 μm, we decided to only consider enstatite grains for simplicity here. Our dusty BT, DP, and ER grids fit the extremely red color as well as the very shallow 1.4 μm water absorption feature of HD 206893 B well (cf. Fig. 6). However, both the plain and dusty DP grids fail to predict the pointy H-band peak observed with GPI¹⁰. Such a triangular H-band peak is indicative of a low-gravity object for spectral types between M6–L7, although an older and dusty field brown dwarf with a similar spectral characteristic has been observed, too (Allers & Liu 2013). The H-band morphology of HD 206893 B is therefore at least suggestive of a dusty atmosphere, but likely also a low surface gravity as predictedby the BT and ER grids. Therefore, we obtain the best fits with our dusty BT and ER grids with χ² = 0.751 and 0.757, respectively.In agreement with the comparison between HD 206893 B and other dusty, low-gravity or young M and L-dwarfs by Delorme et al. (2017) and Ward-Duong et al. (2021), we find a grain size distribution with a geometric mean size of ~ 0.3 μm for these grids. Again, we note that the best fit dusty DP model deviates significantly from the best fit dusty BT and ER models at longer wavelengths (>2.5 μm) and better data, for example from the James Webb Space Telescope, is required to obtain tighter constraints.

5.3 Free retrievals

The spectral retrievals, which explore a larger range of atmospheric parameters, yield spectra that are largely consistent with the best fit dusty BT and ER models regarding the spectral morphology between 1–2.5 μm and the predicted flux at longer wavelengths (cf. Fig. 7). The retrieved T_eff from the free retrievals (cf. Table 7) is much below the values from the grid retrievals (cf. Table 6). This is the case because for the grid retrievals T_eff is computed for the interpolated spectrum without the additional extinction applied while for the free retrievals T_eff is obtained by integrating over the final spectrum. Both the free retrievals and the grid retrievals with BT and ER point to a low-gravity atmosphere, which, combined with the retrieved radii, interestingly points to a planetary-mass object. The metallicity and C/O ratio are super-solar in all cases and they are particularly enhanced according to the free retrievals. The high metallicity could also be suggestive of a low-mass object since an anticorrelation between mass and heavy element enrichment has been inferred in the atmospheres of the solar system ice and gas giants and exoplanets (Thorngren et al. 2016), with planet formation models predicting [Fe/H] > 1 only for planets less massive than ~ 1 M_Jup (see Fig. 3 of Mordasini et al. 2016).

The retrieved pressure-temperature profiles from pRT and ATMO (cf. Fig. 8) show different gradients but their photospheres are located at similar pressures (~ 10–100 mbar). The comparison with the condensation profiles in Fig. 8 shows that MgS, MgSiO₃, Mg₂ SiO₄, and Fe clouds are expected to be present in the atmosphere of HD 206893 B. The three latter species would condense out around 0.1 bar when considering the P–T profiles from ATMO although cloud opacities were parametrized during the retrieval without making an assumption about the cloud composition and only limited in their extent by the cloud top pressure for the gray cloud contribution. With pRT, we use the condensation profiles to infer the base of the cloud deck of the considered cloud species (MgSiO₃ and Fe). Figure 8 shows that the cloud base is expected to be deeper in the atmosphere (around 20 bar) when considering the retrieved P–T profiles from pRT.

It is noteworthy that the C/O ratio is constrained with high precision (although a bimodal distribution is obtained with pRT), both with the grid and the free retrievals. This is surprising since the error bars on the C/O ratio are comparable with those obtained for β Pic b by GRAVITY Collaboration (2020) at much better signal-to-noise in the K-band. For the grid retrievals with ER, the best-fit C/O ratios coincide with a discrete point in the model grid. This may indicate that the (linear) interpolation between grid points does not provide sufficiently accurate spectra such that the uncertainties on the C/O ratio from ER are expected to be underestimated. As mentioned above, the posterior distribution of the C/O ratio obtained with pRT is bimodal, which appears related to a correlation with the retrieved quenching pressure (cf. Fig. C.7). With P_quench ~ 1 mbar, the retrieved C/O ratio is comparable to the value from the plain ER fit (i.e., ~ 0.65). For the second solution with P_quench ~ 10 bar, the C/O ratio has a value of ~0.9 and is therefore consistent with ATMO. The retrievals with pRT and ATMO use chemical equilibrium abundances to determine the absorber abundances in the atmosphere. Both models use three free parameters for this, but with different choices on the parametrization: pRT fits for the atmospheric metallicity ([Fe/H]) and then varies the C/O ratio by changing the oxygen abundance. In addition, it retrieves a quenching pressure, above which the abundances of CH₄, H₂O, and CO were held constant. ATMO retrieves the atmospheric metallicity [Fe/H], and allows the C and O abundances to vary freely. The deep quenching pressure that was retrieved with pRT is located below the photosphere. Hence, the combination of P_quench, [Fe/H], andC/O is used as a knob for changing the relative CO/CH₄/H₂O abundances in the photospheric region. In a similar way, the retrieval with ATMO allows for changing the C and O abundances directly, which appears to provide a similar freedom with the abundances retrieval, and, in the retrievals presented here, leads to a similar C/Oratio.

The comparison between the free retrievals presented here shows that a careful vetting of different input model assumptions is important. While retrieving a quenching pressure in the atmosphere appears physically justified (as applied by pRT), as does an independent variation in the oxygen and carbon content in the atmosphere (as applied by ATMO). The fact that both retrievals find consistent C/O constraints could be pure coincidence, and in principle combining both approaches may be best, but it is important to assess whether a retrieval model with a correspondingly increased number of free parameters is justified, for example by considering Bayes factors. Another important difference between the two retrieval models is the parametrization of both the atmospheric P–T structure and the clouds, which are closely linked. The P–T parametrization of pRT can in principle allow for more isothermal P–T structures, while the parametrization of ATMO, which fixes the equilibrium temperature of the planet, will always lead to a non-negligible temperature gradient in the photosphere. As has been shown by, for example, Tremblin et al. (2015, 2016), isothermal atmospheres may mimic clouds, and lead to cloud-free atmospheres even when (synthetic) cloudy spectra are fed into retrievals (Mollière et al. 2020). This may explain why the P–T profile retrieved by pRT is more isothermal than that retrieved by ATMO, but we note that both retrievals constrain the atmospheres to be very cloudy. Additionally, the comparative ease with which the ATMO parametrization can lead to cloudy spectra may be important because in pRT only specific combinations of the physically motivated cloud parameters (settling parameter f_sed, atmospheric mixing strength K_zz, atmospheric P–T profile in comparison to saturation vapor pressure curve of condensates) lead to cloudy solutions.

Apart from atmospheric clouds, another possible source of reddening could be an inclined circumplanetary disk. Evidence for circumplanetary accretion disks around young (< 10 Myr) giant planets has been seen in both planet formation simulations (e.g., Lubow et al. 1999; D’Angelo et al. 2002) and observations (e.g., Bowler et al. 2011; Christiaens et al. 2019). However, recent detections of long-lived accretion disks (“Peter-Pan disks”) around planetary-mass, brown dwarf, and M dwarf objects with ages of ~ 20–55 Myr (Eriksson et al. 2020; Boucher et al. 2016; Murphy et al. 2018; Lee et al. 2020; Silverberg et al. 2020) suggest that a disk around HD 206893 B could be plausible if its age is at the lower end of the range of ~ 50–700 Myr predicted by Delorme et al. (2017). The interferometric visibilities measured by GRAVITY agree with the SPHERE K1- and K2-band photometry within 11.5% and 7.7%, respectively, and we find no dependence of the visibility as a function of the baseline length. Assuming emission from a uniform disk, we thus find an upper limit on the disk diameter of ~89 R_Jup and ~77 R_Jup, respectively, with a longest baseline of 130 m (cf. also Wang et al. 2021, who did a similar exercise for the PDS 70 protoplanets).

5.4 Mass

From the comparison between the best fit atmospheric parameters and the AMES-Cond evolutionary tracks in Sect. 4.3, we find that the dusty ER grid would be consistent with an extremely young object and the plain DP grid would be consistent with a rather old object. In between the two, the dusty BT grid suggests a moderately young (~ 3–300 Myr) object somewhere between ~ 5–30 M_Jup. While the uncertainties on HD 206893 B’s age and mass are large, it is still noteworthy that this solution is also in agreement with both the age of ~ 40–270 Myr estimated for the Argus moving group (Torres et al. 2008; Bell et al. 2015) and the mass of ${10}_{-4}^{+5} M_{\text{Jup}}$ predicted from radial velocity and Gaia data of the system (Grandjean et al. 2019). Compared to Delorme et al. (2017), who report a best fit age and mass of 100–300 Myr and 15–30 M_Jup, respectively, we find that a younger, planetary-mass object would fit the data, too, although this hypothesis hinges on a moderately low Argus moving group membership probability and the dynamical mass estimate from radial velocity and Gaia data is likely biased due to the presence of the inner companion HD 206893 C. A low surface gravity is supported by the spectral retrievals with petitRADTRANS and ATMO, whose best fit atmospheric parameters are roughly consistent with those predicted by the dusty ER grid, but Allers & Liu (2013) mention that older and dusty field brown dwarfs can have very similar H- and K-band spectral features than young low-gravity objects.

By analyzing the Gaia proper motion anomaly of the system (cf. Sect. 3.2), we obtain an independent constraint on the mass of HD 206893 B and find that a second, closer-in companion (HD 206893 C) is required to explain the observed data in accordance with Grandjean et al. (2019). The mass range predicted for HD 206893 B is consistent with the 5–30 M_Jup obtained from the comparison between the best fit atmospheric parameters and evolutionary tracks, except for the case where the two companions are located roughly on the opposite side of their orbits. In this case, their masses remain poorly constrained. The mass of HD 206893 C should be between ~ 8–15 M_Jup and depending on how well the mass of HD 206893 B is known, the on-sky position of HD 206893 C could be narrowed down to a parameter space where searching for it with the GRAVITY instrument could be feasible.

6 Conclusions

We present new VLTI/GRAVITY K-band spectroscopy at a resolution of R ~ 500 of the reddest known substellar companion HD 206893 B. From these observations we obtain two new astrometric data points with a precision of ~ 100 μas as well as a low signal-to-noise K-band spectrum of HD 206893 B. We use the astrometry to update the orbital parameters of HD 206893 B and the spectrum toinfer its atmospheric parameters with atmospheric model fits and spectral retrievals. Given the previously observed difficulties with fitting both the extremely red color as well as the very shallow 1.4 μm water absorption feature of HD 206893 B (Delorme et al. 2017; Ward-Duong et al. 2021), we include additional extinction by high-altitude dust clouds made of enstatite grains in our atmospheric model fits.

From the orbit fits, we find that the GRAVITY data resolves the previously observed degeneracy between a lower eccentricity, larger semimajor axis and a higher eccentricity, smaller semimajor axis orbit by preferring the latter of these. The orbital solution for HD 206893 B preferred by the GRAVITY data suggests a mutual inclination of ${20.8}_{-11.2}^{+13.6} \text{deg}$ between the orbit of HD 206893 B and the debris disk of the system. We argue that a misalignment between them could suggest a significantly younger age for the system or could be caused by dynamical planet–planet interactions with other putative companions in the system. However, such a misalignment needs to be confirmed by future GRAVITY observations and a profound analysis of dynamical planet–planet interactions and their impact on the alignment and eccentricity of HD 206893 B’s orbit is left for future work.

From the grid retrievals, we find that the BT-Settl-CIFIST (BT) and Exo-REM (ER) grids including additional extinction (dusty models) can fit all near-infrared spectral features of HD 206893 B, namely the extremely red color, the very shallow 1.4 μm water absorption feature, and the pointy H-band peak observed with GPI. However, both DRIFT-PHOENIX (DP) grids with (dusty models) and without (plain models) additional extinction fail to reproduce the pointy H-band peak. By comparison to evolutionary tracks, we argue that only the best fit dusty BT and ER models as well as the plain DP model correspond to physically plausible objects. The best fit parameters of these three models are spread over a wide range of ages and masses, though. If the best fit dusty BT and ER models are favored over the best fit plain DP model, due to their ability to fit the pointy H-band peak observed with GPI, we predict an age of ~ 3–300 Myr and a mass of ~ 5–30 M_Jup for HD 206893 B. This mass estimate from atmospheric models and evolutionary tracks is consistent with the mass estimate from the Gaia proper motion anomaly of the system.

From the free retrievals with petitRADTRANS (pRT) and ATMO, we obtain parameters that are roughly consistent with those predicted by the dusty BT and ER models. Most notably are a very low surface gravity, high metallicity, and high C/O ratio. While the atmospheric chemistry and formation of clouds are handled differently between pRT and ATMO, we find consistent C/O ratios of ~ 0.8–0.9. However, their high precision is surprising given the much worse signal-to-noise in the K-band if compared to β Pic b (GRAVITY Collaboration 2020), for which a similar precision was achieved. This might hint at the C/O constraints for HD 206893 B being driven by systematic effects.

The low surface gravity together with the mass estimate of ${10}_{-4}^{+5} M_{\text{Jup}}$ from radial velocity and Gaia data (Grandjean et al. 2019) and a potential Argus moving group membership of the system (membership probability ~ 61%, Ward-Duong et al. 2021) suggests that a planetary nature could be possible for HD 206893 B. While tension exists across the various age indicators, with the extensive host star analysis performed by Delorme et al. (2017) showing both young (rotational period), old (lithium and barium abundance) and ambiguous (chromospheric activity) indicators, the most recent available kinematic analyses from Banyan Sigma and Lee & Song (2019) point to a younger age and potential membership with Argus (Zuckerman 2019). In addition, disentangling youth and low gravity from a dusty atmosphere can be difficult in the L-dwarf regime (Allers & Liu 2013), so that further observations such as more precise L- and M-band photometry from the James Webb Space Telescope and a broader spectral coverage or higher spectral resolution are required to make a robust statement on the nature of HD 206893 B. In agreement with Grandjean et al. (2019), we also find that the Gaia proper motion anomaly of the system suggests a second, closer-in companion (HD 206893 C) whose mass could be in the planetary-mass regime.

Finally, it has been shown that the extreme atmospheric conditions on HD 206893 B responsible for its exceptionally red color cannotbe reproduced by currently available atmospheric models for giant planets and brown dwarfs without further adaptions. The case of HD 206893 B can therefore serve as a benchmark for the further development of such atmosphericmodels which could ultimately lead to a more complete understanding of the objects at the boundary between exoplanets and brown dwarfs. Moreover, future radial velocity or high-contrast imaging observations might confirm the additional companions predicted to exist in this system and improve the mass estimate for HD 206893 B.

Acknowledgements

The authors would like to thank Rob De Rosa for his insights on orbit fitting and the comparison with the GPI results. This research has made use of the Jean-Marie Mariotti Center Aspro service (http://www.jmmc.fr/aspro). P.M. and T.H. acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. T.S. acknowledges support from the Netherlands Organisation for Scientific Research (NWO) through grant VI.Veni.202.230. This work was performed using the ALICE compute resources provided by Leiden University. A.V., G.O., and M.H. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 757561). The manuscript was also substantially improved following helpful comments from an anonymous referee.

Appendix A Orbit fitting and debris disk

Fig. A.1 Same as Fig. 3, but using Gaussian priors of 140 ± 3 deg for the inclination i and 61 ± 4 deg for the longitude of the ascending node Ω in order to enforce coplanarity between the orbit of HD 206893 B and the debris disk of the system.

Appendix B Photometry

Appendix C Atmospheric model fitting posteriors

Fig. C.1 Posterior distribution of the atmospheric parameters of the plain BT-Settl-CIFIST model fitted to the SPHERE, GPI, and GRAVITY spectra and the NACO and SPHERE photometry. The values state the 68% confidence intervals around the median.

Fig. C.2 Same as Fig. C.1, but for the plain DRIFT-PHOENIX model.

Fig. C.3 Same as Fig. C.1, but for the plain Exo-REM model.

Fig. C.4 Same as Fig. C.1, but for the dusty BT-Settl-CIFIST model. We note that the luminosity is calculated from the effective temperature and the radius and is no longer consistent with the luminosity that would be calculated from the reddened spectrum (i.e., the calculated luminosity is larger).

Fig. C.5 Same as Fig. C.4, but for the dusty DRIFT-PHOENIX model.

Fig. C.6 Same as Fig. C.4, but for the dusty Exo-REM model.

Fig. C.7 Posterior distribution of the spectral retrieval with petitRADTRANS.

Fig. C.8 Posterior distribution of the spectral retrieval with ATMO.

References

All Tables

All Figures

	Fig. 1 Schematic of the HD 206893 system with the two inner companions “B” and “C,” the debris disk with its ~ 27 au wide gap, and the planetary-mass companion candidate at ~74 au that could be responsible for clearing this gap. The quoted mass estimates for HD 206893 B and C are explained in Sect. 5.4 of this paper. Distances and sizes are not to scale.
In the text

	Fig. 2 BT-NextGen model spectrum of HD 206893 A (black line), scaled to fit the shown 2MASS and WISE photometry (gray data points). Synthetic photometry based on the model spectrum is also shown with black data points. The top panel shows the transmission curves corresponding to each photometric data point and the bottom panel shows the residuals between the synthetic and the observed photometry.
In the text

	Fig. 3 Posterior distribution of the orbital parameters (top) and orbital solutions together with the NACO, SPHERE, GPI, and GRAVITY astrometry (bottom) of HD 206893 B. In the top panel, the values state the 68% confidence intervals around the median. In the bottom panel, the black star highlights the position of HD 206893 A and all error bars show the 1–σ confidence intervals.
In the text

	Fig. 4 Masses of HD 206893 B (left) and C (right) as a function of the on-sky position of HD 206893 C (at the reference epoch of Gaia EDR3, 2016.0) derived from the proper motion anomaly measured between the Gaia EDR3 proper motion and the HIPPARCOS–Gaia EDR3 long-term proper motion. Physically implausible regions due to negative masses or dynamical instability have been ignored.
In the text

	Fig. 5 Combined GRAVITY K-band spectrum of HD 206893 B together with the SPHERE Y –H-band spectrum from Delorme et al. (2017) and the GPI J, H, K1, and K2-band spectra fromWard-Duong et al. (2021). The shaded regions highlight the 1–σ confidence intervals. For reference, absorption bands of water and carbon monoxide are indicated.
In the text

Fig. 6 Atmospheric models fitted to the observed spectra and photometry of HD 206893 B shown in the background. Dark red lines show the best fit dusty models that include additional extinction by high-altitude dust clouds made of enstatite grains and light red lines show the exact same models before applying the extinction (i.e., without dust). The GPI spectra are not to be taken to face value since they are rescaled during each of the fits. BT = BT-Settl-CIFIST, DP = DRIFT-PHOENIX, and ER = Exo-REM.

In the text

	Fig. 7 Retrieved spectra of HD 206893 B with the observed spectra and photometry shown in the background. The GPI spectra are not to be taken to face value since they are rescaled during the retrieval with petitRADTRANS (pRT).
In the text

Fig. 8 Pressure-temperature profiles of the atmosphere of HD 206893 B retrieved with petitRADTRANS (pRT) and ATMO. The dashed lines show the 1–σ confidence intervals. With pRT we obtain a bimodal posterior for the C/O ratio and the quenching pressure, and the P–T profile for each of the two solutions is shown separately (thick blue line C/O > 0.7, thin blue line C/O < 0.7). All models probe approximately 0.01 to 0.1 bar pressure levels. Dotted lines indicate the condensation curves of various species calculated at 10× solar abundance (Visscher et al. 2010).

In the text

Fig. 9 Parameters inferred for HD 206893 B from atmospheric model fits and spectral retrievals compared to AMES-Cond evolutionary tracks. Light red points show the best fit parameters for the plain models and dark red points show the best fit parameters for the dusty models including additional extinction by high-altitude dust clouds made of enstatite grains. The evolutionary tracks are shown for objects with exactly those masses printed below the colorbar. Curves of constant age are shown in dashed gray. BT = BT-Settl-CIFIST, DP = DRIFT-PHOENIX, ER = Exo-REM, and pRT = petitRADTRANS.

In the text

Fig. 10 J–K (left) and H–K (right) color-magnitude diagrams showing HD 206893 B in red together with the reddening vector of our best fit extra-photospheric enstatite dust model for different V -band extinctions in light red. Other known planetary-mass objects are shown in orange and M, L, and T-dwarfs from the SpeX Prism Spectral Libraries are shown in green. For the dust model, we assume amean = 0.33 μm, σa = 1.30, and AV = 0–3. The red points are in steps of 1 mag.

In the text

	Fig. A.1 Same as Fig. 3, but using Gaussian priors of 140 ± 3 deg for the inclination i and 61 ± 4 deg for the longitude of the ascending node Ω in order to enforce coplanarity between the orbit of HD 206893 B and the debris disk of the system.
In the text

	Fig. C.1 Posterior distribution of the atmospheric parameters of the plain BT-Settl-CIFIST model fitted to the SPHERE, GPI, and GRAVITY spectra and the NACO and SPHERE photometry. The values state the 68% confidence intervals around the median.
In the text

	Fig. C.2 Same as Fig. C.1, but for the plain DRIFT-PHOENIX model.
In the text

	Fig. C.3 Same as Fig. C.1, but for the plain Exo-REM model.
In the text

	Fig. C.4 Same as Fig. C.1, but for the dusty BT-Settl-CIFIST model. We note that the luminosity is calculated from the effective temperature and the radius and is no longer consistent with the luminosity that would be calculated from the reddened spectrum (i.e., the calculated luminosity is larger).
In the text

	Fig. C.5 Same as Fig. C.4, but for the dusty DRIFT-PHOENIX model.
In the text

	Fig. C.6 Same as Fig. C.4, but for the dusty Exo-REM model.
In the text

	Fig. C.7 Posterior distribution of the spectral retrieval with petitRADTRANS.
In the text

	Fig. C.8 Posterior distribution of the spectral retrieval with ATMO.
In the text