<figure><img src="image_1.png" /></figure>

HAL
open science

GRAVITY K -band spectroscopy of HD 206893 B

J. Kammerer, S. Lacour, T. Stolker, P. Mollière, D. Sing, E. Nasedkin, P.

Kervella, J. Wang, K. Ward-Duong, M. Nowak, et al.

► To cite this version:

J. Kammerer, S. Lacour, T. Stolker, P. Mollière, D. Sing, et al. GRAVITY K -band spectroscopy of HD 206893 B. Astronomy and Astrophysics - A&A, 2021, 652, pp.A57. 10.1051/0004-6361/202140749 . insu-03365088

HAL Id: insu-03365088

https://insu.hal.science/insu-03365088

Submitted on 17 Jun 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# GRAVITY K-band spectroscopy of HD 206893 B*

## Brown dwarf or exoplanet

J. Kammerer<fnref n="1" />,²,²¹, S. Lacour³,<fnref n="1" />, T. Stolker⁴, P. Mollière⁵, D. K. Sing²⁶,²⁷, E. Nasedkin⁵, P. Kervella³, J. J. Wang⁶, K. Ward-Duong²¹, M. Nowak⁷, R. Abuter<fnref n="1" />, A. Amorim⁸,⁹, R. Asensio-Torres⁵, M. Bauböck¹⁰, M. Benisty¹¹, J.-P. Berger¹¹, H. Beust¹¹, S. Blunt⁶, A. Boccaletti³, A. Bohn⁴, M.-L. Bolzer¹⁰, M. Bonnefoy¹¹, H. Bonnet<fnref n="1" />, W. Brandner⁵, F. Cantalloube⁵, P. Caselli¹⁰, B. Charnay³, G. Chauvin¹¹, E. Choquet¹², V. Christiaens¹³, Y. Clénet³, V. Coudé du Foresto³, A. Cridland⁴, R. Dembet<fnref n="1" />, J. J. Dexter¹⁴, P. T. de Zeeuw⁴,¹⁰, A. Drescher¹⁰, G. Duvert¹¹, A. Eckart¹⁵,¹⁶, F. Eisenhauer¹⁰, F. Gao¹⁷, P. Garcia⁹,¹⁸, R. Garcia Lopez¹⁹,⁵, E. Gendron³, R. Genzel¹⁰, S. Gillessen¹⁰, J. Girard²¹, X. Haubois²², G. Heißel³, T. Henning⁵, S. Hinkley²³, S. Hippler⁵, M. Horrobin¹⁵, M. Houllé¹², Z. Hubert¹¹,³, L. Jocou¹¹, M. Keppler⁵, L. Kreidberg⁵, A.-M. Lagrange¹¹,³, V. Lapeyrère³, J.-B. Le Bouquin¹¹, P. Léna³, D. Lutz¹⁰, A.-L. Maire²⁴,⁵, A. Mérand<fnref n="1" />, J. D. Monnier²⁰, D. Mouillet¹¹, A. Müller⁵, T. Ott¹⁰, G. P. P. L. Otten¹²,²⁸, C. Paladini²², T. Paumard³, K. Perraut¹¹, G. Perrin³, O. Pfuhl<fnref n="1" />, L. Pueyo²¹, J. Rameau¹¹, L. Rodet²⁵, G. Rousset³, Z. Rustamkulov²⁷, J. Shangguan¹⁰, T. Shimizu¹⁰, J. Stadler¹⁰, O. Straub¹⁰, C. Straubmeier¹⁵, E. Sturm¹⁰, L. J. Tacconi¹⁰, E. F. van Dishoeck⁴,¹⁰, A. Vigan¹², F. Vincent³, S. D. von Fellenberg¹⁰, F. Widmann¹⁰, E. Wieprecht¹⁰, E. Wiezorrek¹⁰, J. Woillez<fnref n="1" />, S. Yazici¹⁰, and the GRAVITY Collaboration

*(Affiliations can be found after the references)*

Received 8 March 2021 / Accepted 2 June 2021

## ABSTRACT

**Context.** Near-infrared interferometry has become a powerful tool for studying the orbital and atmospheric parameters of substellar companions.

**Aims.** We aim to reveal the nature of the reddest known substellar companion HD 206893 B by studying its near-infrared colors and spectral morphology and by investigating its orbital motion.

**Methods.** We fit atmospheric models for giant planets and brown dwarfs and perform spectral retrievals with petitRADTRANS and ATMO on the observed GRAVITY, SPHERE, and GPI spectra of HD 206893 B. To recover its unusual spectral features, first and foremost its extremely red near-infrared color, we include additional extinction by high-altitude dust clouds made of enstatite grains in the atmospheric model fits. However, forsterite, corundum, and iron grains predict similar extinction curves for the grain sizes considered here. We also infer the orbital parameters of HD 206893 B by combining the ~100 µas precision astrometry from GRAVITY with data from the literature and constrain the mass and position of HD 206893 C based on the Gaia proper motion anomaly of the system.

**Results.** The extremely red color and the very shallow 1.4 µm water absorption feature of HD 206893 B can be fit well with the adapted atmospheric models and spectral retrievals. By comparison with AMES-Cond evolutionary tracks, we find that only some atmospheric models predict physically plausible objects. Altogether, our analysis suggests an age of ~3–300 Myr and a mass of ~5–30 $M_{Jup}$ for HD 206893 B, which is consistent with previous estimates but extends the parameter space to younger and lower-mass objects. The GRAVITY astrometry points to an eccentric orbit ($e = 0.29^{+0.06}_{-0.11}$) with a mutual inclination of <34.4 deg with respect to the debris disk of the system.

**Conclusions.** While HD 206893 B could in principle be a planetary-mass companion, this possibility hinges on the unknown influence of the inner companion on the mass estimate of $10^{-5} M_{Jup}$ from radial velocity and Gaia as well as a relatively small but significant Argus moving group membership probability of ~61%. However, we find that if the mass of HD 206893 B is <30 $M_{Jup}$, then the inner companion HD 206893 C should have a mass between ~8–15 $M_{Jup}$. Finally, further spectroscopic or photometric observations at higher signal-to-noise and longer wavelengths are required to learn more about the composition and dust cloud properties of HD 206893 B.

**Key words.** planets and satellites: atmospheres – planets and satellites: detection – planets and satellites: gaseous planets – techniques: interferometric

## 1. Introduction

While the number of known exoplanets has grown to over 4000 in the last decade<fnref n="1" />, 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

<footnote>* Based on observations made with ESO telescopes at Paranal Observatory under program IDs 1103.B-0626 and 1104.C-0651.</footnote>

<footnote n="1">¹ http://exoplanet.eu/</footnote>
<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>

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 \sim 500$) K-band spectroscopy of exoplanets. They were able to demonstrate astrometric measurements with a precision ~10 times better than what had previously been possible and constrain the atmospheric C/O ratio of the gas giant $\beta$ 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 $\beta$ 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 $\beta$ Pic c. In this regard, near-infrared interferometry significantly advances the field of planet formation and evolution by enabling direct observations of 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^{+5}_{-4} M_{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 \sim 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).

## 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 \sim 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 $\tau_0 = 5$ ms) and below average ($\tau_0 \approx 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
<table><caption>**Table 1.** Observing log.</caption><thead><tr><th rowspan="2">Date</th><th colspan="2">UT time</th><th colspan="2">NEXP/NDIT/DIT</th><th rowspan="2">Airmass</th><th rowspan="2">τ<sub>0</sub></th><th rowspan="2">Seeing</th></tr><tr><th>Start</th><th>End</th><th>HD 206893 B</th><th>HD 206893 A</th></tr></thead><tbody><tr><td>2019-07-17</td><td>08:52:56</td><td>09:56:06</td><td>5/12/60 s</td><td>6/64/1 s</td><td>1.17–1.54</td><td>1.5–2.2 ms</td><td>1.11–1.71"</td></tr><tr><td>2019-08-13</td><td>03:21:16</td><td>04:21:09</td><td>5/12/60 s</td><td>7/64/1 s</td><td>1.03–1.12</td><td>3.2–4.9 ms</td><td>0.80–1.00"</td></tr></tbody></table>

**Notes.** NEXP, NDIT, and DIT denote the number of exposures, the number of detector integrations per exposure, and the detector integration time, respectively, and $τ_0$ denotes the atmospheric coherence time.

<table><caption>**Table 2.** Relative astrometry of HD 206893 B.</caption><thead><tr><th>MJD (days)</th><th>ΔRA (mas)</th><th>ΔDec (mas)</th><th>σ<sub>ΔRA</sub> (mas)</th><th>σ<sub>ΔDec</sub> (mas)</th><th>ρ</th></tr></thead><tbody><tr><td>58 681.396</td><td>130.73</td><td>198.10</td><td>0.04</td><td>0.06</td><td>-0.58</td></tr><tr><td>58 708.165</td><td>127.03</td><td>199.27</td><td>0.09</td><td>0.13</td><td>-0.88</td></tr></tbody></table>

**Notes.** The covariance matrix can be obtained using $σ^2_{ΔRA}$ and $σ^2_{ΔDec}$ on the diagonal and $ρσ_{ΔRA}σ_{ΔDec}$ off-diagonal, where $ρ$ is the correlation coefficient.

package developed by our team<fnref n="2" />. 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.

### 3. Orbit analysis

#### 3.1. Astrometric analysis

From the interferometric observations with GRAVITY we obtain two new astrometric data points with an unprecedented precision

<table><caption>**Table 3.** Stellar parameters and 2MASS and WISE photometry of HD 206893 A from the literature.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>T<sub>eff</sub></td><td>K</td><td>6500 ± 100</td><td>D17</td></tr><tr><td>log g</td><td>–</td><td>4.45 ± 0.15</td><td>D17</td></tr><tr><td>[Fe/H]</td><td>dex</td><td>0.04 ± 0.02</td><td>D17</td></tr><tr><td>R</td><td>R<sub>☉</sub></td><td>1.26 ± 0.02</td><td>D17</td></tr><tr><td>π</td><td>mas</td><td>24.51 ± 0.06</td><td>G18</td></tr><tr><td>J<sub>2MASS</sub></td><td>mag</td><td>5.869 ± 0.023</td><td>S06</td></tr><tr><td>H<sub>2MASS</sub></td><td>mag</td><td>5.687 ± 0.034</td><td>S06</td></tr><tr><td>K<sub>S2MASS</sub></td><td>mag</td><td>5.593 ± 0.021</td><td>S06</td></tr><tr><td>W1<sub>WISE</sub></td><td>mag</td><td>5.573 ± 0.176</td><td>W10</td></tr><tr><td>W2<sub>WISE</sub></td><td>mag</td><td>5.452 ± 0.052</td><td>W10</td></tr><tr><td>W3<sub>WISE</sub></td><td>mag</td><td>5.629 ± 0.015</td><td>W10</td></tr><tr><td>W4<sub>WISE</sub></td><td>mag</td><td>5.481 ± 0.043</td><td>W10</td></tr></tbody></table>

**Notes.** D17 = Delorme et al. (2017), G18 = Gaia Collaboration (2018), S06 = Skrutskie et al. (2006), W10 = Wright et al. (2010).

<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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!<fnref n="3" /> (Blunt et al. 2020), which infers the posterior distribution of the

<footnote n="2">² Python package available on GitHub upon request.</footnote>

<footnote n="3">³ https://github.com/sblunt/orbitize</footnote>
<table><caption>**Table 4.** Orbital parameters inferred for HD 206893 B.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th><th>This work <sup>(a)</sup></th><th>G19 <sup>(a)</sup></th><th>W21 <sup>(a)</sup></th><th>This work <sup>(b)</sup></th><th>M20 <sup>(b)</sup></th></tr></thead><tbody><tr><td>a (au)</td><td>Semimajor axis</td><td>LogUniform(1, 100)</td><td>9.28<sup>+1.77</sup><sub>-0.93</sub></td><td>9.74<sup>+1.46</sup><sub>-1.41</sub></td><td>10.4<sup>+1.8</sup><sub>-1.7</sub></td><td>11.37<sup>+1.09</sup><sub>-0.75</sub></td><td>11.35<sup>+1.13</sup><sub>-0.77</sub></td></tr><tr><td>e</td><td>Eccentricity</td><td>Uniform(0, 1)</td><td>0.29<sup>+0.06</sup><sub>-0.11</sub></td><td>0.25<sup>+0.13</sup><sub>-0.17</sub></td><td>0.23<sup>+0.13</sup><sub>-0.16</sub></td><td>0.13<sup>+0.05</sup><sub>-0.03</sub></td><td>0.14<sup>+0.05</sup><sub>-0.04</sub></td></tr><tr><td>i (deg)</td><td>Inclination</td><td>Sine(0, 180) <sup>(a)</sup><br/>N(140, 3) <sup>(b)</sup></td><td>154<sup>+12</sup><sub>-9</sub><br/>—</td><td>146<sup>+13</sup><sub>-7</sub><br/>—</td><td>146<sup>+14</sup><sub>-7</sub><br/>—</td><td>—</td><td>—</td></tr><tr><td>ω (deg)</td><td>Argument of per.</td><td>Uniform(0, 360)</td><td>123<sup>+89</sup><sub>-44</sub></td><td>181<sup>+223</sup><sub>-135</sub></td><td>173<sup>+79</sup><sub>-109</sub></td><td>71<sup>+36</sup><sub>-25</sub></td><td>74<sup>+36</sup><sub>-29</sub></td></tr><tr><td>Ω (deg)</td><td>Lon. of asc. node</td><td>Uniform(0, 180) <sup>(a)</sup><br/>N(61, 4) <sup>(b)</sup></td><td>55<sup>+49</sup><sub>-28</sub><br/>—</td><td>267<sup>+44</sup><sub>-29</sub><br/>—</td><td>77<sup>+34</sup><sub>-33</sub><br/>—</td><td>—</td><td>—</td></tr><tr><td>τ</td><td>Time of per. since τ<sub>ref</sub><sup>(c)</sup></td><td>Uniform(0, 1)</td><td>0.58<sup>+0.10</sup><sub>-0.20</sub></td><td>—</td><td>—</td><td>60<sup>+3</sup><sub>-4</sub></td><td>62<sup>+4</sup><sub>-4</sub></td></tr><tr><td>π (mas)</td><td>Parallax</td><td>N(24.51, 0.06) <sup>(d)</sup></td><td>24.50<sup>+0.06</sup><sub>-0.06</sub></td><td>24.51<sup>+0.06</sup><sub>-0.06</sub></td><td>24.51<sup>+0.06</sup><sub>-0.06</sub></td><td>24.51<sup>+0.06</sup><sub>-0.06</sub></td><td>24.51<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><td>M<sub>tot</sub> (M<sub>☉</sub>)</td><td>Total mass</td><td>N(1.32, 0.02) <sup>(e)</sup></td><td>1.32<sup>+0.02</sup><sub>-0.02</sub></td><td>1.32<sup>+0.02</sup><sub>-0.02</sub></td><td>1.32<sup>+0.02</sup><sub>-0.02</sub></td><td>1.32<sup>+0.02</sup><sub>-0.02</sub></td><td>1.32<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>P (yr)</td><td>Orbital period</td><td>—</td><td>24.59<sup>+7.38</sup><sub>-3.60</sub></td><td>26.45<sup>+6.18</sup><sub>-5.54</sub></td><td>29.1<sup>+8.1</sup><sub>-6.7</sub></td><td>33.36<sup>+4.96</sup><sub>-3.30</sub></td><td>33.28<sup>+5.09</sup><sub>-3.33</sub></td></tr></tbody></table>

**Notes.** The posterior states the 68% confidence interval around the median. $N(μ, σ)$ denotes a normal distribution with mean $μ$ and standard deviation $σ$. For comparison, the orbital parameters inferred by other works are shown as well. Scenario 1 is constrained by the data only and scenario 2 has an additional constraint on coplanarity with the debris disk. <sup>(a)</sup>Scenario 1, <sup>(b)</sup>Scenario 2, <sup>(c)</sup>Measured in fractions of orbital period with $τ_{ref} = 55800$ MJD, <sup>(d)</sup>Gaia Collaboration (2018), <sup>(e)</sup>Delorme et al. (2017), G19 = Grandjean et al. (2019), W21 = Ward-Duong et al. (2021), M20 = Marino et al. (2020).

orbital parameters through Markov Chain Monte Carlo (MCMC) sampling with ptemcee<fnref n="4" /> (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^{+1.77}_{-0.93}$ au and a higher eccentricity of $0.29^{+0.06}_{-0.11}$ 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_m = \arccos(\cos(i) \cos(i_d) + \sin(i) \sin(i_d) \cos(\Omega - \Omega_d)), \quad (1)$$

where $i_d = 140 \pm 3$ deg is the inclination of the debris disk and $\Omega_d = 61 \pm 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 $\Omega_d$. We find $i_m = 20.8^{+13.6}_{-11.2}$ 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 $\Omega$ 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 $\Omega$. 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

<footnote n="4">⁴ https://github.com/willvousden/ptemcee</footnote>
<figure><img src="image_4.png" /><figcaption>**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-\sigma$ confidence intervals.</figcaption></figure>

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 ($PM_{RA} = -0.102 \pm 0.037$ and $PM_{DEC} = -0.612 \pm 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
<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**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 from Ward-Duong et al. (2021). The shaded regions highlight the 1–σ confidence intervals. For reference, absorption bands of water and carbon monoxide are indicated.</figcaption></figure>

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.

## 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
<table><caption>**Table 5.** Prior boundaries for the atmospheric model grids used in this work.</caption><thead><tr><th>Model</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>log <i>g</i></th><th>[Fe/H]<br>(dex)</th><th>C/O</th><th><i>R</i><br>(<i>R</i><sub>Jup</sub>)</th></tr></thead><tbody><tr><td>BT-Settl-CIFIST</td><td>1000–2000</td><td>2.5–5.5</td><td>–</td><td>–</td><td>0.8–6.0</td></tr><tr><td>DRIFT-PHOENIX</td><td>1000–2000</td><td>3.0–5.5</td><td>-0.6–0.3</td><td>–</td><td>0.8–6.0</td></tr><tr><td>Exo-REM</td><td>1000–2000</td><td>3.5–4.5</td><td>-0.5–0.5</td><td>0.3–0.75</td><td>0.8–6.0</td></tr></tbody></table>

**Notes.** The boundaries for the effective temperature and the radius were chosen based on previous works on HD 206893 B, while those for the other parameters exploit the maximum allowed range.

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<fnref n="5" /> (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<fnref n="6" /> (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 circum-planetary 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 ($\tau \sim 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{dn}{da} = \frac{N}{a\sqrt{2\pi}\ln\sigma_a} \exp\left(-\frac{\ln^2(a/a_{\text{mean}})}{2\ln^2\sigma_a}\right), \quad (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*<sub>mean</sub> is the geometric mean grain size, and *σ*<sub>a</sub> 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*<sub>mean</sub> and a uniform prior between 1.1–10 is used for *σ*<sub>a</sub>. 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 *σ*<sub>ext,l</sub>(*a*<sub>mean</sub>, *σ*<sub>a</sub>) using PyMieScatt<fnref n="7" /> (Sumlin et al. 2018) and scale the flux of the default spectra *F*<sub>default</sub> to that of the reddened spectra

$$ F_{\text{red}} = F_{\text{default}} \times 10^{\left(-\frac{a_V}{2.5} \frac{\sigma_{\text{ext,l}}}{\sigma_{\text{ext,V}}}\right)}, \quad (3) $$

where $A_V$ is the extinction in the Bessel V-band, another free parameter with a uniform prior between 0–5 mag, and $\sigma_{\text{ext,V}}$ is the extinction cross-section averaged over the Bessel V-band. In total, our enstatite dust model has three free parameters ($a_{\text{mean}}$, $\sigma_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^2_{\text{red}}$. This

<footnote n="5">⁵ https://github.com/tomasstolker/species</footnote>

<footnote n="6">⁶ https://github.com/JohannesBuchner/PyMultiNest</footnote>

<footnote n="7">⁷ https://github.com/bsumlin/PyMieScatt</footnote>
<table><caption>**Table 6. Atmospheric parameters inferred for HD 206893 B using grid retrievals.**</caption><thead><tr><th>Model</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>log <i>g</i></th><th>[Fe/H]<br>(dex)</th><th>C/O</th><th><i>R</i><br>(<i>R</i><sub>Jup</sub>)</th><th><i>M</i><br>(<i>M</i><sub>Jup</sub>)</th><th><i>a</i><sub>mean</sub><br>(μm)</th><th>σ<sub>a</sub></th><th><i>A</i><sub>V</sub><br>(mag)</th><th><i>χ</i><sup>2</sup><sub>red</sub></th></tr></thead><tbody><tr><td colspan="11"><strong>Plain</strong></td></tr><tr><td>BT</td><td>1600<sup>+1</sup><sub>-1</sub></td><td>3.50<sup>+0.00</sup><sub>-0.01</sub></td><td>—</td><td>—</td><td>0.98<sup>+0.01</sup><sub>-0.01</sub></td><td>1.17<sup>+0.02</sup><sub>-0.02</sub></td><td>—</td><td>—</td><td>—</td><td>0.957</td></tr><tr><td>DP</td><td>1431<sup>+9</sup><sub>-9</sub></td><td>5.14<sup>+0.20</sup><sub>-0.13</sub></td><td>0.27<sup>+0.02</sup><sub>-0.04</sub> (↑)</td><td>—</td><td>1.20<sup>+0.03</sup><sub>-0.03</sub></td><td>77<sup>+45</sup><sub>-21</sub></td><td>—</td><td>—</td><td>—</td><td>0.841</td></tr><tr><td>ER</td><td>1049<sup>+2</sup><sub>-4</sub></td><td>3.50<sup>+0.00</sup><sub>-0.00</sub> (↓)</td><td>0.49<sup>+0.01</sup><sub>-0.01</sub> (↑)</td><td>0.65<sup>+0.00</sup><sub>-0.00</sub></td><td>2.32<sup>+0.02</sup><sub>-0.02</sub></td><td>6.58<sup>+0.11</sup><sub>-0.10</sub></td><td>—</td><td>—</td><td>—</td><td>1.024</td></tr><tr><td colspan="11"><strong>Dusty</strong></td></tr><tr><td>BT</td><td>1589<sup>+13</sup><sub>-22</sub></td><td>3.83<sup>+0.38</sup><sub>-0.14</sub></td><td>—</td><td>—</td><td>1.17<sup>+0.07</sup><sub>-0.05</sub></td><td>3.55<sup>+5.67</sup><sub>-1.08</sub></td><td>0.33<sup>+0.05</sup><sub>-0.06</sub></td><td>1.30<sup>+0.15</sup><sub>-0.12</sub></td><td>1.99<sup>+0.44</sup><sub>-0.50</sub></td><td>0.751</td></tr><tr><td>DP</td><td>1444<sup>+12</sup><sub>-11</sub></td><td>5.02<sup>+0.14</sup><sub>-0.14</sub></td><td>0.26<sup>+0.03</sup><sub>-0.05</sub> (↑)</td><td>—</td><td>1.75<sup>+0.25</sup><sub>-0.16</sub></td><td>123<sup>+56</sup><sub>-35</sub></td><td>2.29<sup>+0.11</sup><sub>-0.10</sub></td><td>1.17<sup>+0.07</sup><sub>-0.05</sub></td><td>0.82<sup>+0.23</sup><sub>-0.18</sub></td><td>0.774</td></tr><tr><td>ER</td><td>1347<sup>+6</sup><sub>-7</sub></td><td>3.55<sup>+0.06</sup><sub>-0.04</sub> (↓)</td><td>0.06<sup>+0.09</sup><sub>-0.07</sub></td><td>0.75<sup>+0.00</sup><sub>-0.01</sub> (↑)</td><td>2.03<sup>+0.08</sup><sub>-0.08</sub></td><td>5.71<sup>+1.01</sup><sub>-0.65</sub></td><td>0.34<sup>+0.04</sup><sub>-0.04</sub></td><td>1.34<sup>+0.11</sup><sub>-0.08</sub></td><td>2.87<sup>+0.36</sup><sub>-0.30</sub></td><td>0.757</td></tr></tbody></table>

**Notes.** The plain models are unmodified whereas the dusty models include additional extinction by high-altitude dust clouds made of enstatite grains. The values state the 68% confidence intervals around the median. Arrows (↑ or ↓) indicate if a parameter converges toward the upper or lower grid boundary.

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^2_{red}$) for the plain models, it yields the worst fit (i.e., the highest $\chi^2_{red}$) 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^2_{red}$ 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
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 7.** Atmospheric parameters inferred for HD 206893 B using free retrievals.</caption><thead><tr><th>Model</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>log <i>g</i></th><th>[Fe/H]<br>(dex)</th><th>[C/C<sub>☉</sub>]<br>(dex)</th><th>[O/O<sub>☉</sub>]<br>(dex)</th><th>C/O</th><th><i>R</i><br>(<i>R</i><sub>Jup</sub>)</th><th><i>M</i><br>(<i>M</i><sub>Jup</sub>)</th></tr></thead><tbody><tr><td><i>petitRADTRANS</i></td><td>1216<sup>+13</sup><sub>-17</sub></td><td>2.87<sup>+0.63</sup><sub>-0.47</sub></td><td>1.62<sup>+0.35</sup><sub>-0.41</sub></td><td>—</td><td>—</td><td>0.82<sup>+0.04</sup><sub>-0.19</sub></td><td>1.70<sup>+0.05</sup><sub>-0.05</sub></td><td>0.83<sup>+2.71</sup><sub>-0.55</sub></td></tr><tr><td><i>ATMO</i></td><td>1113<sup>+51</sup><sub>-52</sub></td><td>2.72<sup>+0.24</sup><sub>-0.10</sub></td><td>0.74<sup>+0.16</sup><sub>-0.21</sub></td><td>1.45<sup>+0.10</sup><sub>-0.14</sub></td><td>1.24<sup>+0.09</sup><sub>-0.13</sub></td><td>0.90<sup>+0.03</sup><sub>-0.04</sub></td><td>2.23<sup>+0.09</sup><sub>-0.09</sub></td><td>1.07<sup>+0.82</sup><sub>-0.24</sub></td></tr></tbody></table>

of HD 206893 B that we also used for the atmospheric model fits in Sect. 4.1.

### 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 assumed to be in chemical equilibrium, but with an additional parameter for a quenching pressure above which the CO, CH$_4$, and H$_2$O abundances are fixed. We include CO, H$_2$O, CH$_4$, NH$_3$, CO$_2$, H$_2$S, Na, K, PH$_3$, VO, TiO, FeH as molecular and atomic line species plus Rayleigh scattering and collision induced absorption (CIA) of H$_2$ and He. Given the similar temperature to HR 8799 e, we only consider cloud condensates composed of MgSiO$_3$ 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_{\text{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}}}, \qquad (4)$$

where $X_{\text{cloud},0}$ is the cloud mass fraction at the cloud base, $P_{\text{base}}$ is the pressure at the cloud base, and $f_{\text{sed}}$ is the settling parameter (assumed to be the same for MgSiO$_3$ and Fe). The cloud parameter, $\tau_{\text{cloud}}$, is then defined as the optical depth at the $\tau = 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}}), \qquad (5)$$

where $X_{\text{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 atmospheric column 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
<figure><img src="image_8.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>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).</figcaption></figure>

the best fit models of pRT and ATMO (described below), the latter of which did not retrieve scaling factors.

### 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_0(\lambda/\lambda_0)^{\alpha_{\text{haze}}}, \qquad (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}_2}, \qquad (7) $$

where κ(λ)cloud is the “cloud” absorption opacity, δcloud is an empirical factor governing the strength of the gray scattering, and κH₂ is the scattering opacity due to H₂ at 0.35 µm. σ(λ)haze and κ(λ)cloud are added to the total gaseous scattering and absorption, respectively, with a further parameter specifying
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

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.

### 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 evolutionary tracks (Baraffe et al. 2003). Figure 9 shows these tracks for objects of different masses from 1–100 $M_{\text{Jup}}$ in the radius versus effective temperature and surface gravity versus effective temperature planes. The best fit atmospheric parameters 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}$ K and $1113_{-52}^{+51}$ 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_{\text{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_{\text{Jup}}$. The plain DP grid suggests a rather old (~100–1000 Myr) object somewhere between ~15–75 $M_{\text{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_{\text{Jup}}$, the (dusty) BT grid predicts a moderately young (~3–300 Myr) object of ~5–30 $M_{\text{Jup}}$, and the (plain) DP grid predicts a rather old (100–1000 Myr) object of ~15–75 $M_{\text{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<fnref n="8" /> 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).

<footnote n="8">⁸ http://pono.ucsd.edu/~adam/browndwarfs/spexprism/library.html</footnote>
<figure><img src="image_11.png" /><figcaption>**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 *a*<sub>mean</sub> = 0.33 µm, *σ*<sub>a</sub> = 1.30, and *A*<sub>V</sub> = 0–3. The red points are in steps of 1 mag.</figcaption></figure>

Apparent magnitudes were converted to absolute magnitudes using distances (i.e., parallaxes) from SIMBAD<fnref n="9" />.

The light red lines show the reddening vectors of our best fit enstatite dust model for different V-band extinctions. Here, we assume *a*<sub>mean</sub> = 0.33 µm and *σ*<sub>a</sub> = 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*<sub>V</sub> = 0–3. For the best fit dusty BT and ER models, the predicted V-band extinction of *A*<sub>V</sub> ~ 2.0 and *A*<sub>V</sub> ~ 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*<sub>H</sub> 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*<sub>V</sub> ~ 2–3 obtained for HD 206893 B in a realistic regime.

# 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

<footnote n="9">⁹ http://simbad.u-strasbg.fr/simbad/</footnote>
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 to exist 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<fnref n="10" />. 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 predicted by the BT and ER grids. Therefore, we obtain the best fits with our dusty BT and ER grids with $\chi^2 = 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

<footnote n="10">¹⁰ We note that since we allow for different scaling parameters for each of the GPI spectra during each of the fits, we cannot use them to obtain information about the absolute flux, but they are still useful to compare the spectral morphology between data and model.</footnote>
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_{\text{quench}} \sim 1$ mbar, the retrieved C/O ratio is comparable to the value from the plain ER fit (i.e., $\sim 0.65$). For the second solution with $P_{\text{quench}} \sim 10$ bar, the C/O ratio has a value of $\sim 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$_4$, H$_2$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_{\text{quench}}$ [Fe/H], and C/O is used as a knob for changing the relative CO/CH$_4$/H$_2$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/O ratio.

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_{\text{sed}}$, atmospheric mixing strength $K_{\text{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 circumbinary disk. Evidence for circumbinary 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_{\text{Jup}}$ and ~77 $R_{\text{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_{\text{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^{+5}_{-4} 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_{\text{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_{\text{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_{\text{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 \sim 500$ of the reddest known substellar companion
HD 206893 B. From these observations we obtain two new astro-
metric 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 to infer its atmospheric parameters with atmo-
spheric 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 eccentric-
ity, 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 sug-
gests a mutual inclination of $20.8_{-11.2}^{+13.6}$ 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_{\text{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} 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 (chromo-
spheric 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 pre-
cise L- and M-band photometry from the James Webb Space
Telescope and a broader spectral coverage or higher spectral res-
olution 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 sug-
gests 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 con-
ditions on HD 206893 B responsible for its exceptionally red
color cannot be 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 atmospheric
models which could ultimately lead to a more complete under-
standing 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.

References

Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872

Allard, F., Homeier, D., & Freytag, B. 2012, Phil. Trans. Roy. Soc. London Ser. A, 370, 2765

Allers, K. N., & Liu, M. C. 2013, ApJ, 772, 79

Amundsen, D. S., Baraffe, I., Tremblin, P., et al. 2014, A&A, 564, A59

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Baudino, J. L., Bézard, B., Boccaletti, A., et al. 2015, A&A, 582, A83

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Biller, B. A., & Bonnefoy, M. 2018, Exoplanet Atmosphere Measurements from Direct Imaging, eds. H. J. Deeg, & J. A. Belmonte, 101

Blunt, S., Wang, J. J., Angelo, I., et al. 2020, AJ, 159, 89

Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2020, MNRAS, 492, 431

Boucher, A., Lafrenière, D., Gagné, J., et al. 2016, ApJ, 832, 50

Bowler, B. P. 2016, PASP, 128, 102001

Bowler, B. P., Liu, M. C., Kraus, A. L., Mann, A. W., & Ireland, M. J. 2011, ApJ, 743, 148

Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, AJ, 159, 63

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172

Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJ, 877, L33

Currie, T., Burrows, A., Madhusudhan, N., et al. 2013, ApJ, 776, 15

D'Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647

Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79

Drummond, B., Tremblin, P., Baraffe, I., et al. 2016, A&A, 594, A69

Eriksson, S. C., Asensio Torres, R., Janson, M., et al. 2020, A&A, 638, L6

Evans, T. M., Sing, D. K., Kataria, T., et al. 2017, Nature, 548, 58

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, Open J. Astrophys., 2, 10

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration 2018, VizieR Online Data Catalog: I/345

Gizis, J. E., Faherty, J. K., Liu, M. C., et al. 2012, AJ, 144, 94

Goyal, J. M., Mayne, N., Sing, D. K., et al. 2018, MNRAS, 474, 5158

Goyal, J. M., Mayne, N., Drummond, B., et al. 2020, MNRAS, 498, 4680

Grandjean, A., Lagrange, A. M., Beust, H., et al. 2019, A&A, 627, L9

GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94

GRAVITY Collaboration (Lacour, S., et al.) 2019, A&A, 623, L11

GRAVITY Collaboration (Nowak, M., et al.) 2020, A&A, 633, A110

GRAVITY Collaboration (Abuter, R., et al.) 2021, A&A, 647, A59

Greco, J. P., & Brandt, T. D. 2016, ApJ, 833, 134

Guillot, T. 2010, A&A, 520, A27

Heintz, W. D. 1978, Double stars, 15

Helling, C., Ackerman, A., Allard, F., et al. 2008a, MNRAS, 391, 1854

Helling, C., Dehn, M., Woitke, P., & Hauschildt, P. H. 2008b, ApJ, 675, L105

Hiranaka, K., Cruz, K. L., Douglas, S. T., Marley, M. S., & Baldassare, V. F. 2016, ApJ, 830, 96

Kervella, P., Arenou, F., & Schneider, J. 2020, A&A, 635, L14

Kirkpatrick, J. D., Gelino, C. R., Cushing, M. C., et al. 2012, ApJ, 753, 156

Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211

Lacour, S., Eisenhauer, F., Gillessen, S., et al. 2014, A&A, 567, A75

Lacour, S., Wang, J. J., Nowak, M., et al. 2020, in SPIE Conf. Ser., 11446, 1144600

Lagrange, A. M., Meunier, N., Rubini, P., et al. 2019, Nat. Astron., 3, 1135

Lapeyrere, V., Kervella, P., Lacour, S., et al. 2014, in SPIE Conf. Ser., 9146, Optical and Infrared Interferometry IV, eds. J. K. Rajagopal, M. J. Creech-Eakman, & F. Malbet, 91462D

Lecavelier Des Etangs, A., Pont, F., Vidal-Madjar, A., & Sing, D. 2008, A&A, 481, L83

Lee, J., & Song, I. 2019, MNRAS, 486, 3434

Lee, J., Song, I., & Murphy, S. 2020, MNRAS, 494, 62

Looper, D. L., Kirkpatrick, J. D., Cutri, R. M., et al. 2008, ApJ, 686, 528

Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001

Marino, S., Zurlo, A., Faramaz, V., et al. 2020, MNRAS, 498, 1319

Marocco, F., Day-Jones, A. C., Lucas, P. W., et al. 2014, MNRAS, 439, 372

Milli, J., Hibon, P., Christiaens, V., et al. 2017, A&A, 597, L2

Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67

Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131

Mordasini, C. 2013, A&A, 558, A113

Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41

Murphy, S. J., Mamajek, E. E., & Bell, C. P. M. 2018, MNRAS, 476, 3290

Nowak, M., Lacour, S., Lagrange, A. M., et al. 2020, A&A, 642, A2

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16

Patience, J., King, R. R., De Rosa, R. J., et al. 2012, A&A, 540, A85

Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62

Rajan, A., Rameau, J., De Rosa, R. J., et al. 2017, AJ, 154, 10

Schmidt, T. O. B., Neuhäuser, R., Seifahrt, A., et al. 2008, A&A, 491, 311

Silverberg, S. M., Wisniewski, J. P., Kuchner, M. J., et al. 2020, ApJ, 890, 106

Skemer, A. J., Morley, C. V., Zimmerman, N. T., et al. 2016, ApJ, 817, 166

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Stephens, D. C., Leggett, S. K., Cushing, M. C., et al. 2009, ApJ, 702, 154

Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182

Sumlin, B. J., Heinson, W. R., & Chakrabarty, R. K. 2018, J. Quant. Spectr. Rad. Transf., 205, 127

Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64

Torres, C. A. O., Quast, G. R., Melo, C. H. F., & Sterzik, M. F. 2008, Young Nearby Loose Associations, ed. B. Reipurth, 5, 757

Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17

Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19
van Leeuwen, F. 2007, A&A, 474, 653

Visscher, C., Lodders, K., & Fegley, Bruce, J. 2010, ApJ, 716, 1060

Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919

Wakeford, H. R., Sing, D. K., Kataria, T., et al. 2017, Science, 356, 628

Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148

Ward-Duong, K., Patience, J., Follette, K., et al. 2021, AJ, 161, 5

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Wu, Y., & Lithwick, Y. 2011, ApJ, 735, 109

Zuckerman, B. 2019, ApJ, 870, 27

Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57

¹ European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany
e-mail: j.kammerer.bw@googlemail.com

² Research School of Astronomy & Astrophysics, Australian National University, Canberra, ACT 2611, Australia

³ LESIA, Observatoire de Paris, PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Janssen, 92195 Meudon, France

⁴ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands

⁵ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

⁶ Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

⁷ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁸ Universidade de Lisboa - Faculdade de Ciências, Campo Grande, 1749-016 Lisboa, Portugal

⁹ CENTRA - Centro de Astrofísica e Gravitação, IST, Universidade de Lisboa, 1049-001 Lisboa, Portugal

¹⁰ Max Planck Institute for extraterrestrial Physics, Giessenbacherstraße 1, 85748 Garching, Germany

¹¹ Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

¹² Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

¹³ School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Melbourne, Australia

¹⁴ JILA and Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA

¹⁵ Institute of Physics, University of Cologne, Zülpicher Straße 77, 50937 Cologne, Germany

¹⁶ Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, 53121 Bonn, Germany

¹⁷ Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany

¹⁸ Universidade do Porto, Faculdade de Engenharia, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal

¹⁹ School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

²⁰ Astronomy Department, University of Michigan, Ann Arbor, MI 48109, USA

²¹ Space Telescope Science Institute, Baltimore, MD 21218, USA

²² European Southern Observatory, Casilla 19001, Santiago 19, Chile

²³ University of Exeter, Physics Building, Stocker Road, Exeter EX4 4QL, UK

²⁴ STAR Institute/Université de Liège, Belgium

²⁵ Center for Astrophysics and Planetary Science, Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

²⁶ Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD, USA

²⁷ Department of Earth & Planetary Sciences, Johns Hopkins University, Baltimore, MD, USA

²⁸ Academia Sinica, Institute of Astronomy and Astrophysics, 11F Astronomy-Mathematics Building, NTU/AS campus, No. 1, Sect. 4, Roosevelt Rd., Taipei 10617, Taiwan
**Appendix A: Orbit fitting and debris disk**

<table><caption>**Table A.1.** Relative astrometry of HD 206893 B from the literature.</caption><thead><tr><th>MJD<br>(days)</th><th>Sep.<br>(mas)</th><th>PA<br>(deg)</th><th>Inst.</th><th>Band</th><th>Source</th></tr></thead><tbody><tr><td>57 300</td><td>270.0 ± 2.6</td><td>69.95 ± 0.55</td><td>SPHERE</td><td><i>H</i></td><td>M17</td></tr><tr><td>57 608</td><td>269.53 ± 12.15</td><td>62.76 ± 2.16</td><td>NACO</td><td><i>L'</i></td><td>S20</td></tr><tr><td>57 647</td><td>265 ± 2</td><td>62.25 ± 0.11</td><td>SPHERE</td><td><i>K1/K2</i></td><td>D17</td></tr><tr><td>57 653</td><td>267.6 ± 2.9</td><td>62.72 ± 0.62</td><td>GPI</td><td><i>H</i></td><td>W21</td></tr><tr><td>57 682</td><td>265.0 ± 2.7</td><td>61.33 ± 0.64</td><td>GPI</td><td><i>K1</i></td><td>W21</td></tr><tr><td>57 948</td><td>260.3 ± 2.0</td><td>54.2 ± 0.4</td><td>SPHERE</td><td><i>H</i></td><td>G19</td></tr><tr><td>58 066</td><td>256.9 ± 1.1</td><td>51.01 ± 0.35</td><td>GPI</td><td><i>K2</i></td><td>W21</td></tr><tr><td>58 277</td><td>246.51 ± 21.34</td><td>42.80 ± 2.24</td><td>NACO</td><td>NB4.05</td><td>S20</td></tr><tr><td>58 289</td><td>249.11 ± 1.60</td><td>45.50 ± 0.37</td><td>SPHERE</td><td>H2/H3</td><td>G19</td></tr><tr><td>58 385</td><td>251.7 ± 5.4</td><td>42.6 ± 1.6</td><td>GPI</td><td><i>J</i></td><td>W21</td></tr><tr><td>58 415</td><td>239.12 ± 17.55</td><td>42.53 ± 2.17</td><td>NACO</td><td><i>M'</i></td><td>S20</td></tr></tbody></table>

**Notes.** “Sep.” and “PA” denote the angular separation and the position angle, respectively, and “Inst.” denotes the instrument with which the data were acquired. M17 = Milli et al. (2017), D17 = Delorme et al. (2017), G19 = Grandjean et al. (2019), S20 = Stolker et al. (2020), W21 = Ward-Duong et al. (2021).
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>
## Appendix B: Photometry

<table><caption>Table B.1. Photometry of HD 206893 B from the literature.</caption><thead><tr><th>Inst.</th><th>Band</th><th>λ<br>(μm)</th><th>Value<br>(mag)</th><th>Source</th></tr></thead><tbody><tr><td>SPHERE</td><td><i>J</i></td><td>1.245</td><td>18.33 ± 0.17</td><td>D17</td></tr><tr><td>SPHERE</td><td><i>H</i></td><td>1.626</td><td>16.79 ± 0.06</td><td>D17</td></tr><tr><td>SPHERE</td><td><i>K1</i></td><td>2.104</td><td>15.20 ± 0.10</td><td>D17</td></tr><tr><td>SPHERE</td><td><i>K2</i></td><td>2.255</td><td>14.88 ± 0.09</td><td>D17</td></tr><tr><td>NACO</td><td><i>L'</i></td><td>3.805</td><td>13.80 ± 0.31</td><td>S20</td></tr><tr><td>NACO</td><td>NB405</td><td>4.056</td><td>13.17 ± 0.55</td><td>S20</td></tr><tr><td>NACO</td><td><i>M'</i></td><td>4.781</td><td>12.78 ± 0.51</td><td>S20</td></tr></tbody></table>

Notes. “Inst.” denotes the instrument with which the data were acquired. D17 = Delorme et al. (2017), S20 = Stolker et al. (2020).

## Appendix C: Atmospheric model fitting posteriors

<figure><img src="image_13.png" /><figcaption>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.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. C.2. Same as Fig. C.1, but for the plain DRIFT-PHOENIX model.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. C.3. Same as Fig. C.1, but for the plain Exo-REM model.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**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).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. C.5. Same as Fig. C.4, but for the dusty DRIFT-PHOENIX model.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. C.6. Same as Fig. C.4, but for the dusty Exo-REM model.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. C.7. Posterior distribution of the spectral retrieval with petitRADTRANS.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. C.8. Posterior distribution of the spectral retrieval with ATMO.</figcaption></figure>