Constraining the properties of HD 206893 B. A combination of radial velocity, direct imaging, and astrometry data

High contrast imaging enables the determination of orbital parameters for substellar companions (planets, brown dwarfs) from the observed relative astrometry and the estimation of model and age-dependent masses from their observed magnitudes or spectra. Combining astrometric positions with radial velocity gives direct constraints on the orbit and on the dynamical masses of companions. A brown dwarf was discovered with the VLT/SPHERE instrument in 2017, which orbits at $\sim$ 11 au around HD 206893. Its mass was estimated between 12 and 50 $M_{Jup}$ from evolutionary models and its photometry. However, given the significant uncertainty on the age of the system and the peculiar spectrophotometric properties of the companion, this mass is not well constrained. We aim at constraining the orbit and dynamical mass of HD 206893 B. We combined radial velocity data obtained with HARPS spectra and astrometric data obtained with the high contrast imaging VLT/SPHERE and VLT/NaCo instruments, with a time baseline less than three years. We then combined those data with astrometry data obtained by Hipparcos and Gaia with a time baseline of 24 years. We used a MCMC approach to estimate the orbital parameters and dynamical mass of the brown dwarf from those data. We infer a period between 21 and 33{\deg} and an inclination in the range 20-41{\deg} from pole-on from HD 206893 B relative astrometry. The RV data show a significant RV drift over 1.6 yrs. We show that HD 206893 B cannot be the source of this observed RV drift as it would lead to a dynamical mass inconsistent with its photometry and spectra and with Hipparcos and Gaia data. An additional inner (semimajor axis in the range 1.4-2.6 au) and massive ($\sim$ 15 $M_{Jup}$) companion is needed to explain the RV drift, which is compatible with the available astrometric data of the star, as well as with the VLT/SPHERE and VLT/NaCo nondetection.


Introduction
Owing to their wide separations from their host stars (>20 au), dynamical masses of directly imaged planet and brown dwarf (BD) companions are difficult to measure. Thus, most of the current mass estimates for these companions rely on evolutionary models. These models still require qualibration, especially at young ages. Such calibrations can be performed using systems for which the companion mass is independently measured.
Indirect measurements of the host star can also provide important constraints on companion properties. In particular, measuring the radial velocity (RV) of a host star enables measurement or at least a constraint on the dynamical mass, M c sin i, of a companion. Calissendorff & Janson (2018), Snellen & Brown (2018), and Dupuy et al. (2019) demonstrated that HIPPARCOS and Gaia measurements of the proper motion of a host star can also be used to better constrain the dynamical mass of its companion.
HD 206893 is a nearby F5V star, located at 40.81 +0.11 −0.11 pc (Gaia Collaboration 2018). Zuckerman & Song (2004)  mated an age between 0.2 and 2 Gyr, which was later refined to 250 +450 −200 Myr by Delorme et al. (2017) using Desidera et al. (2015) method. Additional characteristics of the stars are presented in Table 1. Milli et al. (2017) discovered a BD companion around HD 206893 via direct imaging using the High-contrast Exoplanet REsearch (SPHERE) instrument (Beuzit et al. 2008) and the NaCo instruments at the Very Large Telescope (VLT). The companion has a projected separation of 270 mas (11 au). The star also appears to host a barely resolved disk inclined by 40 ± 10 • from pole-on. From their data taken in 2015-2016, Milli et al. (2017 estimated a period of ∼37 yr, in the case of low eccentricity, and the mass of the companion from models in the range 24 to 73 M J , for a corresponding age in the range 0.2 to 2 Gyr. Using these data and VLT/SPHERE epochs taken in 2016, Delorme et al. (2017) refined the orbital parameters, finding a period P of ∼27 yr and mass M c between 12 and 50 M J for a corresponding age in the range 50−700 Myr. As the age of the system is poorly known, the estimated mass of HD 206893 B, derived from brightness-mass relations, is still poorly constrained.
We report on the High Accuracy Radial velocity Planet Searcher (HARPS) observations of HD 206893, which reveal

Parameter
Value Ref  (Mayor et al. 2003), as part of our Young Nearby Stars survey (YNS; Lagrange et al. 2013). A total of 66 spectra were obtained with a sun altitude below −10 • . We use our Spectroscopic data via analysis of the Fourier Interspectrum RVs (SAFIR) software (Galland et al. 2005) to compute the RVs, which are shown in Fig. 1a.

Direct imaging astrometry
In this paper, we used all epochs obtained for HD 206893 with VLT/SPHERE in Infra-Red Dual Imaging and Spectrograph (IRDIS) mode and with VLT/NaCo. Milli et al. (2017) obtained the first epochs in 2015 and 2016. Then, as part of the SpHere INfrared survey for Exoplanets (SHINE) survey (Chauvin et al. 2017), HD 206893 was observed with SPHERE in September 2016, whose results are reported in Delorme et al. (2017). Finally, HD 206893 was observed two more times with SPHERE in 2017, as part of the follow-up of the SPHERE High-Angular Resolution Debris Disk Survey (SHARDDS; Milli et al. 2017), and in 2018, as part of the SHINE survey. The data were reduced and calibrated as described in Maire et al. (2016) and Delorme et al. (2017). Table 2    HIPPARCOS and Gaia probe the position of the photocenter of the system, which is not the same as the barycenter of the system. Once corrected from the difference in position of the star induced by the companion, the proper motion deduced from the difference of position between HIPPARCOS and Gaia measurements corresponds to the proper motion of the barycenter of the system. The difference between the proper motion as measured by HIPPARCOS and Gaia and the barycenter proper motion provides information on the variation of the tangential velocity over the duration of both datasets (two years for HIPPARCOS, more than two years for Gaia).
These measurements thus constrain both the orientation of the orbit and the dynamical mass of the companion because the impact of the companion on the proper motion of the star should be within the uncertainties of the measurements. Snellen & Brown (2018) and Dupuy et al. (2019) have provided an example of such constraints gained by combining HIPPARCOS and Gaia measurements with direct imaging and RV constraints for the planet β Pic b.
Both the HIPPARCOS and Gaia five-parameter solutions are well behaved for HD 206893 (see Appendix D.1). We present the proper motion measured by HIPPARCOS and Gaia and the proper motion deduced from the difference of position in Table 3.

Radial velocity analysis
The RVs show a strong dispersion due to high frequency variations (hereafter HF). In addition, a linear trend is clearly present (see Fig. 1).
The HF variations seem to have periods slightly less than one day (see Fig. 1a), which is consistent with the photometric variations reported in Delorme et al. (2017). Given the star's spectral type of F5V, the variations could a priori either be due to pulsations or spots. The bisector velocity span (BVS; Queloz et al. 2001) versus RV diagram shows a vertical spread on the 2018 observations and our last epoch shows sinusoidal variations over one hour (see Fig. 2). Such a short timescale indicates that the HF variations are due to pulsations rather than magnetic activity. Additionally, vertical (BVS,RV) spread also indicates pulsations rather than spots because spots usually produce a correlation between RV and BVS (see examples in Lagrange et al. 2009Lagrange et al. , 2013. A linear regression gives a slope of 92 m s −1 yr −1 for the trend. However, because of the very poor data sampling in 2017, coupled with high amplitude stellar variations, the slope is not   well constrained. Binning the RV would also lead to a poor estimation of the slope. To better constrain the slope of the trend, we used a genetic algorithm that fits the RV variations with a sum of 3 sine functions with periods between 0.5 an 2 days to fit the pulsations plus a linear slope to fit the long-term RV variations and an offset to account for the fact that the RV are relative to a reference spectrum (Galland et al. 2005). Fitting three sine function components provides a good overall fit without over-fitting. We used 500 different seeds to estimate the uncertainties associated with the fit. We obtained a median slope of 87 +16 −14 m s −1 yr −1 and a median offset of 1 +0.17 −0.15 km s −1 . Given the resulting (offset,slope) distribution, we can associate to each date a distribution of RV with a linear model. We present this RV distribution in Figure 3. Our RV observations fit broadly into three distinct time epochs. For each of these epochs, we estimated the corresponding RV of the drift and its 1−2σ errorbar from the median of the RV distribution and its 1−2σ limits at the mean date of the epoch. The resulting RVs are presented in red in Fig. 3. We use these RVs in the rest of the study.

Direct imaging astrometry
To constrain the orbital parameters of the BD we fit the astrometry of the companion using a Bayesian approach. We used the model described in Bonnefoy et al. (2014)  We describe the parameters and our priors in Appendix E.1.
We present the results of the MCMC run in Fig. A.1, and the 1−2σ intervals in Table 4. At the end of the MCMC run, the median reduced χ 2 of the chain is 1.1, with a standard deviation of 0.4, ensuring that the data are well fit by the model. We find a period consistent within 1σ with 27 yr period estimated by Delorme et al. (2017); in our results P is in the range 21−33 yr. The orbit is most likely eccentric, between 0 and 0.5, and has a maximum of probability at 0.37. We find an inclination in the range 20−41 • from pole-on, which is consistent with the inclination of the disk, i.e., i Disk = 40 ± 10 • , and the inclination of the equatorial plane of the star, i.e., i * = 30 ± 5 • . The orbit of the BD is then likely coplanar with the disk and in the equatorial plane of its host star. We note a strong degeneracy between period and inclination. We also find a degeneracy between inclination and eccentricity and between the period and eccentricity. We do not discriminate between extended and inclined circular orbits and shorter, eccentric and close to pole-on orbits.

Direct imaging astrometry and RV analysis
To constrain the dynamical mass of the BD, we assumed that the observed companion was responsible for the RV drift then fit the combined RV, presented in Sect. 3, and direct imaging data with a MCMC. Our approach is similar to that used by Bonnefoy et al. (2014) on β Pic b. We used the six orbital parameters, the mass of the BD, and an offset in the RV in our orbital fit. This model was implemented into the likelihood function of our MCMC. We describe the parameters and our priors in Appendix E.2. We present the results of the MCMC run in Appendix B. At the end of the MCMC run, the median reduced χ 2 of the chain is 1.1, with a standard deviation of 0.4, ensuring that the data are well fit by the model. We find similar orbital parameters compared to the fit on the astrometric position of the companion alone. The dynamical mass is now constrained, essentially between 60 and 720 M J (2σ) with a maximum of probability at 140 M J . This mass is incompatible at 2σ with the 12−50 M J mass inferred from its observed photometry and spectra by Delorme et al. (2017).

Direct imaging astrometry, RV, and star proper motion analysis
The proper motion induced by HD 206893 B during HIPPARCOS and Gaia observations should be consistent with the difference between the proper motion measured by HIPPARCOS and Gaia and the barycenter proper motion. The HIPPARCOS and Gaia measurements provide mass constraints complementary to those from the RV observations, as the observed proper motion traces the tangential velocity of the star induced by the BD. HIPPARCOS and Gaia astrometry also yields information on the orientation of the orbit and is thus complementary to direct imaging astrometry as well. None of the solutions at the end of the MCMC on the companion relative astrometry and RV fit HIPPARCOS or Gaia measurements in right ascension and declination, as the amplitude of the tangential velocity is too high (see Fig. B.3).
We add the model to compute the proper motion induced by the companion and the reference proper motion in the likeli-hood function of the MCMC (see Appendices D.3 and D.2). We describe the parameters and our priors in Appendix E.3. Then, we run the MCMC with the direct imaging astrometric data, the three RV epochs presented in Sect. 3 and the proper motion measured by HIPPARCOS and Gaia. The results are presented in Appendix C. At the end of the MCMC run, the median reduced χ 2 of the chain is 3.4; the standard deviation is 0.5, which is significantly worse than the results on relative astrometry and RV. The median χ 2 on the RV data is six times higher than for the solutions of the MCMC on the direct imaging astrometry and RV. We find a mass of 10 +5 −4 M J , which is consistent at 1σ with the mass inferred from the photometry and spectra of HD 206893 B by Delorme et al. (2017); this result favors the hypothesis of young ages. This indicates that HIPPARCOS and Gaia are probing the dynamical influence of the BDs but not the RVs.
However, the median χ 2 on the companion astrometry is slightly worse than for the solutions of the MCMC on the direct imaging alone, but not significantly. Moreover, HIPPARCOS and Gaia proper motion are well fitted except for the Gaia proper motion in RA that is not fitted because the HD 206893 B tangential velocity contribution in RA was null during Gaia observations. We also find longer periods, P between 35 and 54 years with a maximum at 39 years, which are only consistent at 2σ with the periods found in Sect. 4. Adding HIPPARCOS and Gaia constraints favors lower eccentricities, with a maximum at 0.14, and orbits misaligned with the star equator plane, about i = 45 +3 −3 from pole-on. In order to fit HIPPARCOS and Gaia data, the MCMC selected results that fit the companion astrometry less well and therefore were discarded by the MCMC done in Sects. 4 and 5, which shows that HIPPARCOS and Gaia data do not probe the dynamical influence of the BD only.

Discussion
Assuming that the RV drift in the RV of HD 206893 is due to HD 206893 B leads to a mass that is incompatible at 2σ with the 12−50 M J mass inferred from its observed photometry and spectra by Delorme et al. (2017) and with HIPPARCOS and Gaia proper motion measurements.
An additional inner body that would contribute significantly to the observed RV drift and the tangential velocity in RA is needed. Its orbital period should be larger than our time baseline of 1.6 yr (1.4 au). We used the following formula from Lazzoni et al. (2018) for the width of the chaotic zone around a planet to estimate the third body's semimajor axis a crit beyond which the system is not stable: where a B is the first companion semimajor axis, e B is its eccentricity, and µ is the mass ratio of the first companion to the star. Using the characteristics of HD 206893 AB estimated by Delorme et al. (2017), M * = 1.32 M , M c = 22 M J , e = 0.31 and a B between 8 and 10 au, we obtain a critical semimajor axis between 2.1 and 2.6 au, which corresponds to periods between 3 and 4 yr and to an angular separation on the sky between 50 and 60 mas. The critical semimajor axis can be underestimated because the presence of the inner companion induces a wobble of the host star, which produces a bias on the relative positions of the BD toward the host star and thus produces a bias on the eccentricity of its orbit (Pearce et al. 2014). To constrain the properties of the inner body, we used a MCMC to fit the RV alone with a sinusoid (e = 0) and for the inclination of the star: i = 30 • from pole-on, leaving free the phase, mass of the second companion, and period between 1.6 and 4 yr. The fit converges to a mass of 15 M J . Such a mass would lead to a contrast lower than 10 −4 in the near infrared. Given the expected contrast and the upper limit on the separation, it is not possible to detect this body with VLT SPHERE. This third body can be compatible with the HIPPARCOS and Gaia constraint as its period is close to their time baseline; the impact of this third body on the proper motion of the star is then averaged. HIPPARCOS and Gaia data imply a dynamical mass of 10 +5 −4 M J for HD 206893 B. This low mass is compatible with an age for the system lower than 50 Myr. Its probability of membership to the Argus association is higher than previously estimated. It would also explain the low luminosity of the BD compared to the models and its redness reported by Delorme et al. (2017). However, we note that this mass can be under or overestimated as the presence of the second companion could impact on the star tangential velocity. More RV and direct imaging data are needed to further constrain both companions of HD 206893 system.
Acknowledgements. We acknowledge support from the French CNRS and from the Agence Nationale de la Recherche (ANR grant GIPSE ANR-14-CE33-0018).
where y are the observations, σ their corresponding variances, and f (X) the model for parameters X.
To facilitate the convergence of the MCMC we chose the following parameters: -log P: the decimal logarithm of the period in years; f T P : a factor for which the time at the periastron is computed, i.e., T P = f T P × P. In this study the T P is in years and corresponds to the closest periastron passage from the first direct imaging epoch; e: the eccentricity; i: the inclination in radian, 0 • and 90 • correspond to pole-on and edge-on orbit, respectively; -Ω: the longitude of the ascending node in radian; -ω: the argument of periapsis in radian. We used the same model as Bonnefoy et al. (2014), which works on a combination of the parameters presented above.
We used uniform priors in the following intervals: Elsewhere, the prior is null. For the inclination, we took a prior that is the product of a uniform prior in cos(i) and a uniform prior in inclination in the range 0 < i < π.

E.2. MCMC on the direct imaging astrometry and RV
We added to the model a Keplerian in RV to model both the astrometry of the companion and the RV. The likelihood function becomes where y are the observations, σ their corresponding variances, and f (X) the model for parameters X.
To facilitate the convergence of the MCMC we chose additional parameters as -M c : the mass of the companion in M J ; f V 0 : a factor for which the offset in RV are computed, V 0 = f V 0 × K(M c , P, e), with K the semi-amplitude of the RV.
We used the same model as Bonnefoy et al. (2014), which works on a combination of the parameters presented above. We used uniform priors for the additional parameters in the following intervals: Elsewhere, the prior is null.

E.3. MCMC on the direct imaging astrometry and RV and proper motion
We added to the model the tangential velocity model described in Appendix D.3 and the barycenter proper motion model described in Appendix D.2 to fit the difference between the proper motion measured by HIPPARCOS and Gaia and the barycenter proper motion. The likelihood becomes where y are the observations, σ their corresponding variances, and f (X) the model for parameters X. We used the parameters and prior described in Appendix E.2.

E.4. Description of the MCMC runs
We used emcee with the prior and likelihood function as described in Appendices E.1 and E.2. In order to facilitate the convergence, we ran two MCMC for each analysis described in Sects. 1, 5 and 6. We ran the first MCMC with 100 walkers that were initialized randomly with a uniform probability distribution for each parameter in the intervals of their prior for 100 000 iterations. Then a second MCMC was performed with 100 walkers initialized with a uniform probability distribution for each parameter in the intervals where the majority of the walkers of the first run were converging. This approach permits us to check the uniqueness of the solution.
To ensure convergence, as our problem show a multimodal posterior, we verified that the number of iterations was greater than ten times the autocorrelation time of the chain for each parameters (Foreman-Mackey et al. 2013).