Direct discovery of the inner exoplanet in the HD206893 system. Evidence for deuterium burning in a planetary-mass companion

Long term precise radial velocity (RV) monitoring of the nearby star HD206893, as well as anomalies in the system proper motion, have suggested the presence of an additional, inner companion in the system. Here we describe the results of a multi-epoch search for the companion responsible for this RV drift and proper motion anomaly using the VLTI/GRAVITY instrument. Utilizing information from ongoing precision RV measurements with the HARPS spectrograph, as well as Gaia host star astrometry, we report a high significance detection of the companion HD206893c over three epochs, with clear evidence for Keplerian orbital motion. Our astrometry with $\sim$50-100 $\mu$arcsec precision afforded by GRAVITY allows us to derive a dynamical mass of 12.7$^{+1.2}_{-1.0}$ M$_{\rm Jup}$ and an orbital separation of 3.53$^{+0.08}_{-0.06}$ au for HD206893c. Our fits to the orbits of both companions in the system utilize both Gaia astrometry and RVs to also provide a precise dynamical estimate of the previously uncertain mass of the B component, and therefore derive an age of $155\pm15$ Myr. We find that theoretical atmospheric/evolutionary models incorporating deuterium burning for HD206893c, parameterized by cloudy atmospheres provide a good simultaneous fit to the luminosity of both HD206893B and c. In addition to utilizing long-term RV information, this effort is an early example of a direct imaging discovery of a bona fide exoplanet that was guided in part with Gaia astrometry. Utilizing Gaia astrometry is expected to be one of the primary techniques going forward to identify and characterize additional directly imaged planets. Lastly, this discovery is another example of the power of optical interferometry to directly detect and characterize extrasolar planets where they form at ice-line orbital separations of 2-4\,au.


Introduction
What distinguishes an extrasolar giant planet (EGP) from a brown dwarf (BD) is still a matter of debate. The current IAU definition of a planet is an object that is below the mass required for thermonuclear fusion of deuterium, which we currently believe occurs at 13 Jupiter masses. Identifying objects near this mass limit will clarify how distinct the boundary between massive planets and brown dwarfs really is (e.g., Mollière & Mor-dasini 2012;Bodenheimer et al. 2013). However, what qualifies an object to be an exoplanet may be more related to its formation mechanism rather than solely its mass (e.g., Chabrier et al. 2007). Hybrid exoplanetary systems containing both a brown dwarf as well as an exoplanet, which presumably formed at the same time from the same protoplanetary disk, are therefore ideal systems to study different possible formation pathways using diagnostics like atmospheric atomic ratios (e.g., Öberg et al. 2011; Article number, page 1 of 13 arXiv:2208.04867v5 [astro-ph.EP] 3 Apr 2023 A&A proofs: manuscript no. aa  Grandjean et al. 2019a) along with samples from the posterior for the model fit to the radial velocities and proper motion anomalies. Right: The posterior distributions for our best fit model of the predicted proper motion anomaly (purple histograms) calculated as the difference between the Gaia EDR3 proper motion and the long-term proper motion as calculated by comparing Gaia with Hipparcos (labelled "HG" in the figure) and displayed individually in terms of Right Ascension and Declination. The vertical blue lines show the actual values, with 1σ uncertainties represented by the dashed lines. Madhusudhan et al. 2014;Mordasini et al. 2016;Mollière et al. 2022) or initial entropy (e.g., Spiegel & Burrows 2012;Marleau & Cumming 2014). Finding and characterizing planetary mass companions near this deuterium burning limit (e.g., Bowler et al. 2013;Bonnefoy et al. 2014;Hinkley et al. 2013), or systems containing both a brown dwarf and an exoplanet, will allow us to clarify how we draw a distinction between these two classes of objects (e.g., Mollière & Mordasini 2012).
HD 206893 is a nearby (40.707±0.067 pc, Gaia Collaboration et al. 2018) F5V star with roughly solar metallicity ranging from [Fe/H] = −0.07 to −0.05 (Holmberg et al. 2007;Gáspár et al. 2016) to +0.01 (Netopil 2017), and a highly uncertain age due to the lack of clear membership to any young kinematic moving groups. Most age estimates of this star have been in the range of 100-300 Myr (Zuckerman & Song 2004;Delorme et al. 2017;Milli et al. 2017;Kammerer et al. 2021), but some works have allowed ages as young as 3 or 50 Myr (Kammerer et al. 2021;Ward-Duong et al. 2021), while another has quoted an age as high as 2.1 Gyr (e.g., David & Hillenbrand 2015). A substellar companion, HD 206893B was first identified in Milli et al. (2017) showing a projected orbital separation of ∼10 au using the VLT SPHERE instrument. The somewhat uncertain system age, as well as an unconstrained astrometric orbit, allowed only loose constraints on the mass of HD 206893B, initially placing it in the range of 24-73 M Jup (Milli et al. 2017), but later refined to be as low as ∼5 M Jup (Kammerer et al. 2021) and up to 30-40 M Jup (Delorme et al. 2017;Kammerer et al. 2021;Ward-Duong et al. 2021). This object possesses extraordinarily red infrared colors, and has been characterized in numerous subsequent works (Delorme et al. 2017;Stolker et al. 2020;Meshkat et al. 2021;Ward-Duong et al. 2021;Kammerer et al. 2021).
Evidence has been emerging for an additional, inner companion in the system. Using precise radial velocity (RV) monitoring with the HARPS spectrograph over 1.6 yr, Grandjean et al. (2019a) revealed a significant RV drift that could not be solely due to HD 206893B, and hypothesized the presence of a ∼15 M Jup companion. Similarly, Kammerer et al. (2021) showed that the anomaly in the Gaia EDR3 proper motion cannot be explained by the B companion alone. They constrained the parameter space spanned by the mass and on-sky position of an additional, hypothesized 8-15 M Jup companion.
In this paper, using optical interferometry over several epochs, we present the discovery of HD 206893c, a ∼12-13 M Jup exoplanet orbiting the star at 3.5 au. Our study marks the first direct detection of HD 206893c, providing highly precise astrometry (∼50-100 µarcsec) of the HD 206893c orbit, as well as medium resolution (R ∼ 500) spectroscopy in the K-band. In addition to obtaining a precise dynamical mass for the HD 206893c exoplanet, our study also provides precise astrometry of the orbit of HD 206893B. This increased precision on the HD 206893B orbit, combined with its very well constrained luminosity, allows us to derive a well constrained system age of 155±15 Myr. This discovery of HD 206893c establishes HD 206893 as a hybrid planetary system hosting both an exoplanet and a brown dwarf, and therefore presents itself as a valuable laboratory for studying possible formation pathways of EGPs and BDs.

Precision radial velocities
As described in Grandjean et al. (2019a), HD 206893 has been monitored for several years since 2016 with the High Accuracy Radial velocity Planet Searcher (HARPS, La Silla, ESO) as part of the Young Nearby Stars survey (YNS; Lagrange et al. 2013). The Spectroscopic data via Analysis of the Fourier Interspectrum RVs (SAFIR) software (Galland et al. 2005) was used to compute the RV signal. Figure 1 (left panel) shows the processed historical RV datasets taken from Grandjean et al. (2019a) that were used in this study to predict the semi-major axis of HD 206893c. The residuals between the RV data points and the posterior curves shown in the figure are likely due to stellar pulsations. As shown in Table A.1, there is an additional σ RV = 40 m s −1 amplitude best-fit stellar jitter term not currently represented in the uncertainties on the individual RV points shown in Figure 1. However, as discussed in Section 3.1, such residuals in the RV data will not have a significant impact on the final derived physical parameters for HD 206893c.

GRAVITY observations & data processing
All the GRAVITY observations presented in this paper were acquired using the dual-field/on-axis mode, with medium spectral resolution (R = 500). We closely followed the observing strategy outlined in Gravity Collaboration et al. (2019), where the science fiber is alternately centered on the star to obtain the phase reference, and on the position of HD 206893c. These observations were reduced using the procedure described in Gravity Collaboration et al. (2020) and Nowak et al. (2020). The coherent flux is extracted for the on-star and on-planet exposures using the GRAVITY pipeline (Lapeyrere et al. 2014). The on-planet data are then phase-referenced to the on-star observations, and the data are fitted with a model which includes both a planet component and a low-order polynomial (for speckle deconvolution). The astrometry is obtained by varying the astrometry of the planet component, and looking for the minimum χ 2 or, equivalently, by maximising the periodogram power. More details are provided in Section 2.3 as well as Appendix B of Nowak et al. (2020).
Following the work from Kammerer et al. (2021), we used the previous GRAVITY astrometry of HD 206893B in combination with the new HARPS observations (Grandjean et al. 2019a,b) and Gaia DR3 astrometry (Gaia Collaboration et al. 2016, 2021 to predict a position for HD 206893c in August 2021. We use the proper motion anomaly calculated as the difference between the Gaia EDR3 proper motion and the longterm proper motion as calculated by comparing Gaia with Hipparcos proper motions (labelled "HG" in Figure 1). The previous GRAVITY astrometry measurements were particularly useful for constraining the orbit of HD 206893B, and subsequently remove its signal from the RV data and Gaia astrometry. The RV data allowed us to constrain the period, and thus the semimajor axis of HD 206893c, while the Gaia astrometry provided constraints on its position angle.
During the two nights between 27 and 28 August, we used the GRAVITY instrument (Gravity Collaboration et al. 2017) to try to detect the exoplanet. GRAVITY only probes a small circular region of ∼50 mas radius around the position of the single mode fiber. We searched for the companion over a region on the order of 10 000 mas 2 by offsetting the science fiber of the instrument. At the same time the fringe tracker observed the star to correct atmospheric turbulence ). This allowed an integration time of several tens of seconds. The log of the observations, as well as the position of the single mode fibers, are summarised in Table 1. During the second night, on the last pointing, a weak signal was detected.
The exoplanet position was confirmed one month later, in September. A clear detection was obtained during runs ID 1104.C-0651, part of the ExoGRAVITY large program , and the confirmation during a dedicated run (ID 105.20T0.001). In October, a final observation was carried out, as part of the ExoGRAVITY large program, with the goal of obtaining a good calibration of the spectrum. Figure 2 provides examples of the periodogram power maps obtained from our ob- servations, both for a non-detection (fiber mispointed during the search), and for a detection of the planet. The K-band contrast relative to the host star for HD 206893c is approximately 8.2×10 −5 , compared to 1.6×10 −4 for HD 206893B. The reliability of this contrast ratio is ensured Article number, page 3 of 13  Table 2. GRAVITY astrometry of HD 206893B and HD 206893c. The astrometry with are data published in Kammerer et al. (2021) but reduced with the latest version of the GRAVITY pipeline. The co-variance matrix has σ 2 ∆RA and σ 2 ∆Dec on the diagonal, and ρ × σ ∆RA × σ ∆Dec on the off-diagonal, where ρ is the correlation coefficient.
since the host star is observed simultaneously using the GRAV-ITY fringe tracker, so the ratio of coherent energy is constantly monitored, and is the result of a simultaneous measurement. The main systematic uncertainty that can arise in this measurement is due to the response of the injection of the signal into the fiber, defined by the "lobe" of the interferometer, which can be computed analytically as described in Appendix A of Wang et al. (2021). Previous calibrations reveal that these systematics are small.
The astrometry obtained at the three epochs is given in Table 2. In addition to HD 206893c, we have recomputed the astrometry of HD 206893B from Kammerer et al. (2021), and added an additional position obtained in August 2021. With respect to the theoretical performance of GRAVITY (10 µas, Lacour et al. 2014), the astrometry accuracy is still a factor 5 to 20 worse. The difference during the night of 27 August can be explained by the mispointing of the science fiber while searching for the exoplanet. The error bars still do not decrease considerably during the other runs because of the limitation due to systematics, particularly the optical aberrations within the fiber injection optics (GRAVITY Collaboration et al. 2021).

Statistical significance of the detection
In Fig. 2 we show an example of the periodogram power maps for three observations of HD 206893c: one in August 2021, obtained when the fiber was mispointed during the search, and two obtained in September and October 2021, when the fiber was properly centered on the planet. Each of the three maps represents the field of view of the GRAVITY science fiber. A clear peak with some sidelobes in the periodogram power (not visible in the August 2021 observation) are characteristic of a successful detection, and give the position of the planet.

Orbital fits
To assess the orbital parameters for HD 206893B and c, we employ the method presented in Lacour et al. (2021) to perform Bayesian parameter estimation with the orbitize! software package (Blunt et al. 2020) for jointly fitting the orbits of multiple substellar companions using relative astrometry from GRAVITY, the stellar RVs described in the previous section, and Hipparcos and Gaia absolute stellar astrometry from the EDR3 edition of the HGCA (Brandt 2021;Kervella et al. 2022). This method accounts for shifts in the system barycenter due to multiple planets in the system, but does not account for planet-planet interactions, which are negligible on these timescales. The covariances of the uncertainties in RA/Dec due to u-v coverage are also automatically handled by orbitize!. To find the best fit orbits for the HD 206893B and c objects, the parallel-tempered affine-invariant sampler (Foreman-Mackey et al. 2013;Vousden et al. 2016) was run using 20 temperatures, 1000 walkers per temperature, and 60,000 steps per walker. The walker chains were visually inspected for convergence, and the last 10,000 steps from each walker at the lowest temperature was used to form the posterior distributions for each parameter.
The results of our MCMC orbital fits are summarized in Table 3 and a depiction of the orbits for both HD 206893B and c are shown in Figure 3. The best estimate of the mass and semimajor axis of HD 206893c is 12.7 +1.2 −1.0 M Jup and 3.53 +0.08 −0.06 au (∼5.7 yr orbital period) for a dynamical fit that incorporates the precise RV monitoring from HARPS. Excluding the RV information from the fit results in values of 11.5 +2.4 −2.2 M Jup and 3.68 +0.12 −0.09 au. These values are largely consistent with, but also significantly more precise than, the previous estimations of these parameters for the hypothesized HD 206983c companion of 8-15 M Jup and 5.6 au (Grandjean et al. 2019a;Kammerer et al. 2021). At the same time, our fits reveal a precise and significantly nonzero eccentricity e = 0.41 +0.03 −0.03 for HD 206893c (e = 0.36 +0.05 −0.06 when excluding the RV measurements). For completeness we also summarize the results of our orbital fit for HD 206893B. We find a mass and semi-major axis of M B = 28.0 +2.2 −2.1 M Jup and a B = 9.6 +0.4 −0.3 au corresponding to a ∼26-yr orbital period (26.2 +3.7 −3.6 M Jup and 9.7 +0.8 −0.4 au when excluding the RV measurements). More details on the orbit fitting, including the orbital alignment parameters (e.g., i, ω, Ω) can be found in Appendix A. As for HD 206893c, this mass determination for HD 206893B is consistent with, but also significantly more precise than, previous estimates (e.g. Grandjean et al. 2019a;Meshkat et al. 2021;Ward-Duong et al. 2021;Kammerer et al. 2021) that quoted a mass of 5-30 M Jup , and clearly places HD 206893B in the lowmass brown dwarf regime.
3.2. GRAVITY spectroscopy of HD 206893c: model grid fits to derive atmospheric parameters Fig. 4 shows our collected spectrophotometry for both HD 206893B and c. We estimate the primary atmospheric parameters of HD 206893c (i.e., effective temperature, surface gravity, metallicity) with the methodology outlined in Kammerer et al. (2021). We fit the GRAVITY spectroscopy shown in Fig. 4 with a grid of DRIFT-PHOENIX models (Helling et al. 2008) using the species toolkit (Stolker et al. 2020). We have chosen these models to fit the spectroscopy of HD 206893c for simplicity, as well as to be consistent with those fits carried out in Kammerer et al. (2021) for the HD 206893B companion. Moreover, among the three different model grids used in Kammerer et al. (2021), the DRIFT-PHOENIX grid is the only one that yields a good fit to the extremely red colors of HD 206893B without requiring an additional source of reddening. Even if the DRIFT-PHOENIX grid does not fully capture the full physical picture of the atmosphere, the models empirically fit well, and thus are effective for measuring the bolometric luminosity. We defer a more in-depth characterization of the atmosphere of HD 206893c, involving a larger class of models, to a future work. The best fit parameters for HD 206893B and c are shown in Table 3, and the best fit DRIFT-PHOENIX spectra together with our GRAVITY spectra of both companions as well as other spectrophotometry of B from the literature are shown in Fig. 4. We highlight the bolometric luminosities log(L/L ) = −4.23±0.01 dex and −4.42 +0.02 −0.01 calculated from the best-fit models for HD 206893B and c, respectively. We also find a best-fit metallicity [Fe/H] = 0.27 +0.02 −0.05 and 0.28 +0.02 −0.04 for B and c, respectively, which is notably different from the nearly solar metallicity of the host star ([Fe/H] = 0.04 ± 0.02; Kammerer et al. 2021). However, for the purposes of this paper, the metallicity of the atmospheric model has a negligible effect on the final calculated bolometric luminosity, L bol . Specifically, when the fit is restricted to solar-metallicity models ([Fe/H] = 0.0), the bolometric luminosity of the best-fit model is negligibly smaller (by only 0.02 dex, or 5 %). A similar difference is seen for HD 206893B, in which the luminosity using solar-metallicity models would be log(L/L ) = −4.20 dex, only 0.03 dex higher than in Table 3.
Our bolometric luminosity for HD 206893c is calculated by integrating a model that is fit only to the GRAVITY K-band spectrum at ∼2 µm since there are no photometric measurements at other wavelengths.To test the reliability of this L bol estimate, we have calculated the bolometric luminosity for HD 206893B from a fit only to its K-band spectrum from GRAVITY, ignoring the other available data. This yields log(L/L ) = −4.26 ± 0.01, only 7 % smaller than the −4.23 ± 0.01 value from the best fit to all of the spectrophotometric data at 1-5 µm. Under the assumption that the spectral shapes of the two objects are similar, this suggests that our bolometric luminosity determination for HD 206893c is robust.
The derived bolometric luminosity depends to some extent on the chosen atmospheric model. Therefore, there may be additional systematic uncertainties due to these differences between models. To estimate these systematics, we fit an ensemble of five additional theoretical atmosphere models (i.e., AMES-Dusty, BT-Settl, Exo-REM, petitCODE, as well as a simple blackbody curve) to the K-band GRAVITY spectrum of HD 206893c. As expected, the quality of each fit varied from model to model, and each fit returned a slightly different value for the bolometric luminosity. The rough overall standard deviation in log(L/L ) is 0.20 dex. Thus for the analysis in this paper, going forward we adopt a value of log(L/L ) = −4.42 ± 0.20 for HD 206893c, which should account for systematic uncertainties between models. Some of the returned radii are clearly lower than expected (but not pathologically so) but these deviations are accounted for in our presented uncertainty. We discuss the physical implications of this expanded uncertainty more in Section 4. For the atmospheric model fitting, letting the mass of c be entirely free leads to a clearly incorrect estimation of the mass, with the best-fit surface gravity log g and R implying M c ∼100 M Jup . Therefore, we use the dynamical mass (12.7 +1.2 −1.0 M Jup ; Table 3) as a prior, and obtain 9.9 +1.5 −1.7 M Jup , roughly 20 % smaller than the prior. Nevertheless, the dynamical mass should be treated as much more reliable. For the other atmospheric parameters (R, T eff , log g), we find a similar level of variation of 10-20 % in the values across the models. Nonetheless, the DRIFT-PHOENIX models are the most consistent with the previous model fitting in (Kammerer et al. 2021), and so we report their values in Table 3.

Constraints on the system age and on the cloudiness of HD 206893c
The dynamical mass determination of HD 206893B, precise to 10%, can be combined with the object's bolometric luminosity to derive an age for the system by using cooling tracks. At constant age, the luminosity of hot-start gas giants scales roughly as L bol ∼ M 2 at intermediate luminosities (Arras & Bildsten 2006;Marleau & Cumming 2014), so that the mass ratio M B /M c ≈ 2 A&A proofs: manuscript no. aa Adjusting the initial entropy (Spiegel & Burrows 2012; Marleau & Cumming 2014) of either object does not bring their current luminosities closer. Indeed, an elevated initial entropy will not increase the luminosity of c at its present age because hot starts already assume an initial cooling timescale much shorter than the current age (Marleau & Cumming 2014). Similarly, no reasonable lower initial entropy for B could sufficiently delay the beginning of deuterium burning so that it would be observed in the rising part of a "deuterium flash" (see Fig. 8 of Marleau & Cumming 2014 or Fig. 13 of Bonnefoy et al. 2014). Given that the mass of HD 206893c is close to the deuterium-burning limit M ≈ 11.5-14.5 M Jup (Spiegel et al. 2011;Mollière & Mordasini 2012), nuclear reactions are a likely candidate to bring the luminosity of HD 206893c closer to that of B, and so are clouds, which may play a role in the L/T transition that low-mass objects experience (Dupuy & Liu 2012;Liu et al. 2016). Coincidentally, both can occur at similar effective temperatures (Saumon & Marley 2008), depending on the mass of the object. Therefore, we turn to the models of Saumon & Marley (2008, hereafter SM08). The SM08 models provide luminosities and magnitudes for objects cooling with a cloudy or a clear (cloudfree) atmosphere at all ages, which they show is very similar to respectively COND03 (Baraffe et al. 2003) or DUSTY00 (Chabrier et al. 2000). The cloudiness of the atmosphere influences the cooling rate and thus the luminosity and radius evolution of low-mass objects. The advantage of the SM08 models is that they also provide a "hybrid sequence" using a transition from cloudy at high effective temperature T eff 1400 K to cloudfree at lower T eff 1200 K. This is a simplified model of the L/T transition but has the potential of delaying the cooling differently for HD 206893B and c. Despite the numerous improvements in opacities since 2008, the recent cloud-free Sonora models (Marley et al. 2021) have very similar isochrones to SM08 or COND03. Given this, and since cloudy Sonora models are not yet available, the SM08 models are a good modeling choice, and continue to be applied to modelling substellar objects (e.g., VHS J1256-1257AB b; Dupuy et al. 2022).
In Figure 5 (left panel), we show isochrones for the three varieties of SM08 models: clear (190 Myr), cloudy (155 Myr), and hybrid (155 Myr). They all pass through the well-constrained luminosity of B at its dynamical mass, allowing for relatively tight constraints on the age of the system. A more careful statistical assessment is planned for future work but isochrones at 155 ± 15 Myr (155 ± 20 Myr) for hybrid (cloudy) models pass within one mass errorbar of HD 206893B. Figure 5 shows the important finding that both the cloudy and the hybrid models are able to fit both objects, but that the clear models of SM08 are somewhat inconsistent (to about one sigma) with the luminosity of c. This applies also to the COND03 models, which would require an age of 195 Myr as we verified separately. The 140-170 Myr isochrones that pass through one mass standard deviation of HD 206893B are also easily consistent with the luminosity of HD 206893c (not shown).
Figure 5 also displays cooling curves for the three SM08family models, using the best-fit masses. When using clear atmospheres, the cooling curves of HD 206893c are somewhat underluminous, but are overluminous for HD 206893B. Deuterium burning barely plays a role near the current age of the highermass object, the nuclear fuel being instead spent at an earlier age ∼ 10 Myr (e.g., Burrows et al. 2001). Cloudy atmospheres, whether assumed to be present at all times or only down to T eff = 1400 K (hybrid sequence), change the relative cooling so as to "pinch" the cooling curves together.
Interestingly, according to the hybrid model of SM08, HD 206893c might be in the process of losing its clouds while its luminosity would be the same as an object that is still cloudy. Also for HD 206893B, the transition from a cloudy to a clear atmosphere would change significantly the luminosity compared to a pure cloudy cooling sequence only at later ages t 500 Myr.  Saumon & Marley (2008) matching the dynamical mass and bolometric luminosity of HD 206893B, with the respective age indicated in the legend. The deuterium-burning mass limit, which depends on bulk planetary properties, is shown as a fuzzy pink region (Spiegel et al. 2011). Right: Corresponding cooling curves for the best-fit masses. The system age of 155±15 Myr comes from our isochrone analysis (see text). The thick part of the hybrid cooling curves highlights 1400 K > T eff > 1200 K, with cloudy (clear) atmospheres at earlier (later) times.

HD 206893c
HD 206893B  Table 3. Orbital and physical parameters for HD 206893c, and reevaluated parameters of HD 206893B. M, a and e are obtained from the dynamical fits. Values not listed in parentheses are those which use both the astrometry and RVs during the dynamical orbit fit (and are the most trustworthy), while those in parentheses exclude the RV measurements. R, T eff , log g and log(L/L ) are obtained from fitting the spectrophotometry. The fits use a grid of DRIFT-PHOENIX models assuming no additional dust, and using mass priors of 26 ± 2 and 12 ± 1 M Jup for B and c, respectively. The contrast refers to that measured in the K-band. The luminosity uncertainties for HD 206893c are only the statistical ones within the fits to the DRIFT-PHOENIX models. Further analyses should use instead the more appropriate value of 0.2 dex that reflects the systematic modelling uncertainties (Section 3.2 and 4).

Discussion
The highly precise astrometry delivered by GRAVITY (∼50-100 µarcsec) let us derive dynamical masses for the objects in the HD 206893 system that are precise at the level of 5-10%. Continued astrometric monitoring of both objects could, of course, lead to tighter constraints on the derived dynamical masses, and thus yield a more precise age determination. In combination with SM08 models, the masses and luminosities reveal that both cloudy models and hybrid cloudy models featuring a cloudy-toclear transition describe well the cooling of both objects, and that a clear atmosphere for both objects is not consistent with the observations. Thus already the bolometric luminosities, independently of spectroscopy, disfavour the scenario in which both companions have cooled with a clear atmosphere up to now. However, while a 140-Myr SM08 clear isochrone is overluminous at the mass of of HD 206893B, it does match the luminosity of HD 206893c. Thus both objects could have different degrees of cloudiness.
It is unlikely that unknown systematic errors in our dynamical analysis are leading to calculated values of the mass of HD 206893B and c that are erroneous. Our calculation of the dynamical mass uses a combination of stellar astrometry that may have systematic uncertainties that are not well understood, as well as RV measurements that are impacted by stellar activity (e.g., Grandjean et al. 2019a). However, Table 3 shows that excluding the RV data from our fits gives very similar calculated values for the mass of HD 206893c (12.7 +1.2 −1.0 M Jup versus 11.5 +2.4 −2.2 M Jup with the RVs excluded). It should be noted that in this analysis we did not model the stellar activity in our orbit fit, which may contribute to a slight bias. However, it is unlikely that such unknown systematic errors within our analysis have significantly biased our calculated masses, especially given that our measured dynamical mass of HD 206893c is very much consistent with the predictions presented in Grandjean et al. (2019a) and Kammerer et al. (2021). At most, these systematics could have led to underestimated astrometric uncertainties, but this would not have changed the overall conclusions of this study.
As discussed previously, we compared the spectra from the different atmosphere model families, including SM08 (spectra courtesy of D. Saumon, 2022), to the GRAVITY data. In several cases the shape is reproduced only approximately, and the shape of the SM08 models is not red enough. However, there is a large amount of correlated noise within the GRAVITY spectrum, precluding a simple statistical analysis. To address this, taking a conservative approach to estimate the systematic mod-elling errorbar on the luminosity of HD 206893c, we have therefore adopted the scatter in the L bol value from an ensemble of different models. We allow the K-band shape not to match exactly as a way for accounting for detailed physics currently missing in the atmospheric models, while assuming that they capture the overall luminosity evolution. The spectral energy distribution over a small wavelength range, for instance the K band, is much more model-dependent than the overall bolometric flux.
Finally, in Appendix C we describe our study of the dynamical stability of the HD 206893 system. Our analysis demonstrates that while the architecture of this system may resemble that of the β Pictoris system, the HD 206893 system is much closer to instability.

Conclusions
We have presented the discovery of an inner ∼12 − 13 M Jup exoplanet at 3.5 au in the HD 206893 system, making this one of the first directly imaged "hybrid" planetary systems containing both a brown dwarf (Milli et al. 2017) and a bona fide exoplanet. This is also an example of a discovery of an additional companion following a previous discovery of an outer one (e.g., Lagrange et al. 2019;Nowak et al. 2020). The highly precise astrometry delivered by GRAVITY (∼50-100 µarcsec) has enabled us to derive precise dynamical masses of 12.7 +1.2 −1.0 M Jup for HD 206893c and 28.0 +2.2 −2.1 M Jup for HD 206893B. The dynamical mass of B, combined with its luminosity, allows us to immediately and unambiguously establish a robust age of ∼150 Myr for this system. Further, having both an EGP and a BD in the same system, with a common age, HD 206893 will serve as an extraordinarily valuable system for the purpose of furthering our understanding of planet formation, and highlighting distinctions between the formation pathways of EGPs and BDs.
We have used the theoretical atmospheric models from SM08, which self-consistently treat deuterium burning in substellar objects to show the bolometric luminosities of both objects can be well-modelled by cloudy/hybrid models that allow us to fine tune the age value to 155±15 Myr. Besides residing in a hybrid exoplanetary system, HD 206893c is an object straddling the deuterium-burning limit, and thus is an ideal laboratory to establish the precise mass where this process can occur. In addition to having a dramatic impact on its evolution, the duration and extent of deuterium burning in a substellar object may give valuable clues to its initial conditions and internal structure (e.g., the presence/lack of a rocky core, initial entropy, and initial deuterium abundance, Spiegel & Burrows 2012;Mollière & Mordasini 2012;Bodenheimer et al. 2013;Mordasini 2013). Ours is the first such study to begin addressing these questions, and additional studies like this may help to craft a more robust discrimination between EGPs and BDs.
Finally, in addition to being only the second directly imaged exoplanet whose presence was hinted at using the RV method, our precise dynamical mass for HD 206893c means that this effort potentially represents the first discovery via direct imaging of a bona fide exoplanet that partially relies on precise measurements of the host star astrometry from the Gaia mission. This complements numerous other recent discoveries of substellar companions (e.g., Bonavita et al. 2022;Kuzuhara et al. 2022;Franson et al. 2022a,b;Currie et al. in prep.) whose discovery also used Gaia astrometry. Using precise Gaia astrometry of the host star to point the way to orbiting planets suitable for direct imaging is expected to be one of the primary strategies for direct exoplanet detection and characterization going forward, thereby ending the current era of "blind" direct imaging EGP searches that have had notoriously low detection rates (e.g., Bowler & Nielsen 2018). Secondly, recent exoplanet direct imaging surveys have had limited sensitivity (e.g. Nielsen et al. 2019;Vigan et al. 2021) to the peak of the orbital distribution of giant exoplanets that coincides with the location of the water ice line at 2-3 au for solar-type stars (e.g., Fernandes et al. 2019;Frelikh et al. 2019;Fulton et al. 2021). However, upcoming interferometric observations using JWST Hinkley et al. 2022) are likely to have the combination of sufficient sensitivity and resolution to reach these orbital zones, providing complementary characterization at 3-5 µm. Until then, along with the discovery of β Pic c (Lagrange et al. 2019;Nowak et al. 2020) this discovery of an exoplanet at 3.5 au is more evidence that optical interferometry now enables direct characterization of these planets at the ice-line orbital separations of 2-3 au where they form.

Quantity
Prior Posterior

Appendix A: Detailed orbit fitting results
We use the default orbital basis in Blunt et al. (2020). The only difference is that τ, the epoch of periastron after a given reference epoch in units of orbital period, uses a reference epoch of MJD 59000. To fit the radial velocity data, we also include two free parameters to fit for a bulk offset in the radial velocity (γ RV ) and a jitter term to inflate the error bars on the radial velocity data to better account for stellar activity (σ RV ). We list the median and 68-percentile credible intervals centered about the median as the error bars in Table A.1. We also computed orbital periods (P B and P c ) for both companions as well as the mutual inclination of their orbits and listed those in Table A.1.  Left is the trajectory in (e B , e c ) phase space, starting from the observed values (orange crosses, which are random solutions to the orbital fit using all available data). We exclude the initial conditions with e B > 0.15 and a B < 9.5 au. The colours represent the corresponding eccentricity evolution period (secular period). The right panel represents the aforementioned period with respect to the mass of b.