© The Authors 2025

1. Introduction

HR 6819 (HD 167128, QV Tel; G = 5.3) is among the first discovered binary systems composed of a classical Be star and a bloated pre-subdwarf (pre-sd) stripped star (Bodensteiner et al. 2020a; El-Badry & Quataert 2021), and the first whose nature was unambiguously confirmed by optical interferometry (Frost et al. 2022). This binary system is a product of recent mass transfer between two B-type components, which originally orbited each other on a close orbit before the more massive component evolved to fill its Roche lobe. Most of the mass of this component, currently observed as the pre-sd star, was lost to its companion, which was rejuvenated and spun up to near-critical rotation in the process. After the mass transfer ceased, the spun-up component formed a viscous Keplerian disk, thus becoming a classical Be star (Rivinius et al. 2013). The pre-sd star is currently found in a short-lived evolutionary stage that lasts only a few hundred thousand years (El-Badry & Quataert 2021), which precedes its settling into a much fainter helium-burning subdwarf OB-type (sdOB) star and the spectrum becoming dominated by the Be star. HR 6819 is of great importance to stellar astrophysics in general and to the evolution of binary stars in particular as it represents a newly identified, brief snapshot in the evolution of interacting binaries (Marchant & Bodensteiner 2024). The immediate post-mass-transfer stage of HR 6819 also offers a window into the physics of the mass-transfer process, the surface chemical composition of the newly formed Be star, and the surface chemical composition of the exposed core of the former mass donor.

HR 6819 belongs among the so-called black hole (BH) impostors, which are multiple stellar systems that were initially erroneously claimed to contain dormant (non-accreting) stellar mass BHs (Bodensteiner et al. 2022). HR 6819 specifically was reported to be a hierarchical triple system (Rivinius et al. 2020), only to be shown later to consist solely of a classical Be star and a bloated pre-sd star of similar brightness, orbiting each other on a ∼40-day orbit (Bodensteiner et al. 2020a; Gies & Wang 2020; El-Badry & Quataert 2021). Several other systems similar to HR 6819 have been interferometrically confirmed more recently (Rivinius et al. 2024), including another BH impostor and the first suggested Be + pre-sd binary, LB-1 (Liu et al. 2019; Abdul-Masih et al. 2020; Shenar et al. 2020).

The key to the unambiguous confirmation of the nature of HR 6819 and similar systems is the employment of high-angular-resolution observations. They include the single-telescope techniques of speckle interferometry and adaptive-optics-assisted imaging and, more importantly, multiple-telescope optical/near-IR interferometry (Eisenhauer et al. 2023). In the triple-star BH scenario, two widely separated, similarly bright components would be resolved, while in the binary pre-sd scenario, no such outer component would be seen. Despite the tentative report of a speckle companion with unconstrained brightness (Klement et al. 2021), no bright wide companion was detected with follow-up observations by VLT/MUSE, thus confirming the Be + pre-sd scenario (Frost et al. 2022). Furthermore, given the proximity and brightness of HR 6819, optical interferometry obtained with the Very Large Telescope Interferometer (VLTI) was able to resolve the two binary components themselves, enabling the detection of the orbital motion on the basis of two observations separated by ∼13 days (Frost et al. 2022).

Binary systems such as HR 6819 are cornerstone objects for our understanding of the formation of classical Be stars (Shao & Li 2014, 2021; Hastings et al. 2021) and rapid rotators in general (de Mink et al. 2013). Be stars represent up to 40% of B-type stars depending on the metallicity of the environment, and several observational clues indicate that they may be dominated by descendants of HR 6819-like systems. Most apparent is the lack of close, main-sequence companions to (early-type) Be stars (Bodensteiner et al. 2020b), and the majority of Be stars thus appear to be single. When companions are found, they are either evolved products of stripping in mass-transferring binaries, such as neutron stars in Be X-ray binaries (Reig 2011), or their nature remains uncertain (e.g., Gies et al. 2000). The hardest to detect and characterize are stripped sdOB and (pre-)white dwarf (WD) companions, although they are expected to be particularly common if most Be stars originate from mass transfer in binaries (Pols et al. 1991; Shao & Li 2014). Currently, there are only ∼25 confirmed Be + sdOB binaries (Wang et al. 2018, 2021; Klement et al. 2022), and there are only candidate Be + WD systems, among them the class of peculiar X-ray sources with the prototype γ Cas (Smith et al. 2016; Nazé et al. 2022; Gies et al. 2023).

Other observational clues suggest that many Be stars may possess low-mass companions at small orbital separations. Among these, the prevalence of spectral energy distribution (SED) turndown (i.e., the steepening of the spectral slope toward radio wavelengths observed in all Be stars with sufficient radio data) is indicative of disk truncation possibly caused by low-mass companions (Waters et al. 1991; Klement et al. 2019). Furthermore, the structure of the radio SED for the several Be stars with radio SED measurements at multiple wavelengths, including the known binary γ Cas, suggests the presence of a circumbinary extension of the Be disk (Klement et al. 2017). In several Be stars, signs of circum-companion gas – sometimes in the form of a disk – were observed as well (e.g., Bjorkman et al. 2002; Chojnowski et al. 2018; Klement et al. 2024, de Amorim et al., in prep.), including in HR 6819 itself (Bodensteiner et al. 2020a) as well as in LB-1 (Shenar et al. 2020; Hennicker et al. 2022). Recent smoothed particle hydrodynamic simulations of Be disks in binary systems strongly support the notion laid out above, that the Be disks are not truncated by the binary orbit but instead form complex structures, including a bridge region connecting the inner Be disk to the companion, circum-companion structures that may take the form of an accretion disk, and a large-scale circumbinary spiral arm that may contribute to the radio fluxes (Rubio et al. 2025).

The aim of this work is to advance our knowledge of HR 6819 and its place among classical Be stars by obtaining model-independent dynamical masses of the binary system. They were derived with the help of a time series of high-angular-resolution and high-spectral-resolution interferometric observations from VLTI/GRAVITY, which are described in Sect. 2. The fitting of the relative astrometry of both components, as well as of the Be star and its disk parameters, is described in Sect. 3, while the combined three-dimensional orbital solution using previously published radial velocity (RV) measurements from high-resolution FEROS spectroscopy is outlined in Sect. 4. The results are discussed in Sect. 5 in the context of our current understanding of the origin and the general properties of classical Be stars.

2. VLTI/GRAVITY observations

A total of 12 snapshots (1 h observing time including calibrator data and overheads) of HR 6819 were collected between 2021 September and 2023 September by the VLTI (Eisenhauer et al. 2023) four-telescope beam combiner GRAVITY (GRAVITY Collaboration 2017) operating in the near-IR K band (1.98 to 2.40 μm). The 1.8 m Auxiliary Telescope (AT) Array, with each AT using the dedicated adaptive optics system NAOMI (Woillez et al. 2019), was used in the “Large” configuration (maximum baseline of ∼130 m) in order to achieve the highest angular resolution that was possible at VLTI at the time (λ/2B_max ∼ 1.75 mas in the near-infrared K band).

The GRAVITY instrument was used in the single-field on-axis mode, which splits the light equally between the fringe tracker (Lacour et al. 2019) and the science beam combiner. The spectral resolution of the science spectrograph was set to high (λ/Δλ ∼ 4000), enabling a detailed analysis of spectro-interferometric features. Owing to its brightness, HR 6819 was observed in relaxed atmospheric conditions (seeing < 1.4″ and coherence time in the V band τ₀ > 1.4 ms). The data were reduced and calibrated using the official European Southern Observatory pipeline version 1.6.6 in the EsoReflex environment (Freudling et al. 2013). This resulted in the following observables covering the entire K band for each snapshot: spectrum, four closure phases (T3PHI, one per telescope triangle), six absolute visibilities (|V|) and six differential phases (DPHI, one of each per telescope pair).

Each observation of HR 6819 consisted of an object–sky–object sequence, resulting in two sets of data points, which may differ from each other due to the changes in baseline projection caused by the Earth rotation (particularly in |V|; see Fig. 1). Every observation of HR 6819 was followed by an observation of the calibrator HD 161420 (K = 5.4, uniform disk diameter = 0.321 ± 0.008 mas, spectral type F0 IV) in order to enable the absolute calibration of |V| and T3PHI. The calibrator diameter is small considering the angular resolution of the data and the calibrator data thus represent the instrumental response of a point-like source, which is ideal for calibration purposes. The first two calibrated snapshot observations from 2021 were taken in the split polarimetric mode with a Wollaston prism in the beam path, while the rest of the data were obtained in the combined mode. The GRAVITY data (Table A.1) are of very good quality, although for the snapshot taken on 2023 May 01, we had to discard the |V| part of the dataset due to a miscalibration issue (calibrated |V| > 1), probably caused by rapidly changing atmospheric conditions. The full set of the reduced and calibrated data will be made available in the Optical Interferometry Database¹ (Haubois et al. 2014).

Fig. 1. PMOIRED Keplerian geometrical model fit to the calibrated GRAVITY snapshot taken on 2023 July 17, representing the average quality of the model fit (χ2 = 3.1) as well as the typical quality of the data. The upper-left panel shows the (u, v) coverage, the lower-left panel the normalized Brγ line profile (NFLUX), and the remaining panels the interferometric observables closure phase (T3PHI), differential phase (DPHI), and the absolute visibility (|V|). The data (black with gray error bars) are overplotted with the model fit (red). The blue dots in the NFLUX and DPHI panels represent the continuum region used for normalization. The baselines and telescope triangles are given in the corner of each panel, and their colors correspond to the (u, v) coverage panel. For the layout of the VLTI telescope locations, see Fig. 10 in the VLTI manual (https://www.eso.org/sci/facilities/paranal/telescopes/vlti/documents/VLT-MAN-ESO-15000-4552_v115.pdf).

To visualize and model the reduced GRAVITY data, we used the Python code Parametric Modeling of Optical InteRferomEtric Data² (PMOIRED; Mérand et al. 2022). PMOIRED enables the fitting of multicomponent geometrical models composed of, for example, uniform disks (UDs) and Gaussians, to interferometric data, including the possibility to fit normalized line profiles and associated spectro-interferometric signatures when the spectral resolution is sufficient. PMOIRED uses least-squares minimization to obtain best-fit parameters, and contains features such as telluric correction of the K-band spectrum for GRAVITY data, grid search for the detection of binary companions, and a bootstrapping resampling method to obtain realistic uncertainties of the fitted parameters. The grid search method is the same as implemented in the CANDID³ code dedicated specifically to interferometric companion detection (Gallenne et al. 2015).

The fitting of spectro-interferometric features enables the derivation of RVs, whose accuracy depends on the wavelength calibration of the data. The accuracy of GRAVITY in the high spectral resolution mode is known to be at least ∼0.02% (Gallenne et al. 2023), with the pipeline delivering the data at vacuum wavelengths. To check the wavelength calibration precision for each snapshot, we used PMOIRED to fit and subtract telluric features from the K-band spectrum. This was done by interpolating a one-dimensional grid of Molecfit (Smette et al. 2015) models parameterized by the precipitable water vapor while correcting for the shape of the continuum. This procedure resulted in a telluric-corrected K-band spectrum (tellurics do not affect the interferometric observables), the fitted precipitable water vapor, and finally the calibrated spectral dispersion model, which was then applied to the data to facilitate the highest possible accuracy of the fitted RVs. The accuracy of the wavelength calibration resulting from the spectral dispersion model is plotted in Fig. A.1 for each snapshot, revealing that it is better than expected and no worse than 0.004% (∼12 km s⁻¹) for all epochs except for 2021 September 05, where it is 0.007% (∼21 km s⁻¹). To obtain realistic uncertainties of the fitted RVs from GRAVITY data (Sect. 4.3), we thus conservatively chose to quadratically add the uncertainty of 12 km s⁻¹ to all epochs except for 2021 September 05, for which we added the uncertainty of 21 km s⁻¹.

3. Geometrical model fitting

3.1. Relative astrometry and flux ratios of the two binary components

The relative astrometry and flux ratios of the two components were derived using geometrical model fitting with the PMOIRED code to the spectral region around the hydrogen Brγ line, which is the most prominent spectral line in the K band. We selected a 16 nm-wide spectral region centered at the Brγ vacuum wavelength of λ_Brγ,vac = 2.1661178 μm, and normalized the line profile and DPHI using the region edges. Brγ shows a double-peaked emission line profile compatible with what is expected from a Keplerian disk surrounding the Be star and seen at intermediate inclination (see Fig. 1). The spectro-interferometric features in Brγ (T3PHI, DPHI, and |V|) contain a combination of the signal from the Keplerian Be disk, and a much more pronounced signal from the photocenter displacement caused by the bright companion. The companion itself does not appear to show any features in the Brγ region. We fitted the normalized Brγ line profile, and the interferometric observables T3PHI, DPHI, and |V|, simultaneously.

The first two GRAVITY snapshots taken in 2021 were already analyzed by Frost et al. (2022), first using a simple model approximating the Brγ line profile with a Gaussian, and subsequently with a geometrical model mimicking a more realistic contribution from the Be star disk with receding and approaching spatially offset components, which result in a double-peaked line profile similar to the observed one. Our results for the relative astrometry for these two epochs, using a more complex geometrical model (see below), are consistent with these preliminary results.

To obtain a realistic representation of the Be star disk in our geometrical model, we used the prescription for a Keplerian disk that was recently implemented in PMOIRED (Delaa et al. 2011; Meilland et al. 2012; Klement et al. 2024). The full model used to obtain the final fits of the Brγ spectral region consists of the following components: UDs for the two stellar components (with diameters UD_Be and UD_pre-sd) with a relative astrometric offset (ΔRA, ΔDec of the pre-sd star relative to the Be star), and with a relative flux ratio (parameterized as the contribution of the pre-sd star to the total continuum flux f_pre-sd), and a Keplerian disk whose position coincides with that of the Be star. The parameters of the Keplerian disk, contributing solely in the spectral line, are the following: the angular full width at half maximum of the Brγ line emission region (FWHM_disk, Brγ), disk inclination (i_disk, assuming a geometrically thin disk so that the flattening ratio is r = cos i_disk), disk position angle (PA_disk, measured from north to east), the inner diameter of the disk, the Keplerian rotational velocity v_orb at the inner disk diameter (re-parameterized through v_orb sin i_disk), the power law exponent of the radial rotation law (β), the central wavelength of the disk Brγ signature (re-parameterized as RV of the Be star RV_Be), and the equivalent width of the Brγ line profile (EW_Brγ). All parameters were kept free in the fitting process with the following exceptions: both UD_Be and UD_pre-sd were fixed at 0.15 mas (see below), the inner diameter of the Keplerian disk was kept equal to UD_Be, and β was fixed at 0.5 to represent Keplerian rotation. The remaining free parameters were therefore v_orb sin i_disk, RV_Be, ΔRA, ΔDec, f_pre-sd, FWHM_disk, Brγ, i_disk, PA_disk, and EW_Brγ.

As for the choice of the UD values of both components, we ran an initial fit with the value of 0.12 mas, corresponding to the absolute radii given by Bodensteiner et al. (2020a) and El-Badry & Quataert (2021) scaled to their adopted distances. As the distance resulting from the dynamical parallax of our orbital fits (Sect. 4.2) is significantly lower than the adopted distances in the previous works, we repeated the fitting process with the value of 0.15 mas for the UDs, which resulted from rescaling the radii to the smaller distance. We note that the exact values of the UDs do not influence the astrometric results, as in the considered range they represent angularly unresolved components at the angular resolution of the data. However, the value of UD_Be does influence the results for the Keplerian disk, as it corresponds to the disk’s inner diameter. Namely, the higher value of UD_Be results in slightly lower fitted values of the Keplerian velocity at the base of the disk, which is in our model given by v_orb sin i_disk. As we could not fit the value of UD_Be, the exact choice of the best-estimated value of 0.15 mas thus introduces a possible systematic bias into the fitted parameters of the Keplerian disk.

The model including the Keplerian disk described above successfully converged for all epochs with χ_red.² ranging from 1.8 to 5.3, and with the resulting astrometric positions showing a clear orbital motion according to the expectations. A representative example of one of the datasets and the quality of our model fit is shown in Fig. 1, along with synthetic images and component fluxes for the resulting geometrical model shown in Fig. 2. The final values of the fitted parameters and their uncertainties were determined with a bootstrapping algorithm implemented in PMOIRED, the results of which are shown in Fig. A.2 for the same representative epoch. The relative positions and flux ratios for each epoch are listed in Table A.1. We note that the errors do not include contributions from possible systematics such as small-scale miscalibration of the data due to rapidly changing atmospheric conditions during observations.

Fig. 2. PMOIRED synthetic image resulting from the Keplerian model fit to the 2023 July 17 snapshot (left) and the fluxes of the model components (right). The model components are the Be star primary (green) and Keplerian disk (orange) in the image center, and the pre-sd secondary (purple) offset in the image by ∼1.3 mas toward the northeast, with the sum of their fluxes normalized to 1.0 in the continuum (black).

The weighted average of the flux ratio is f_pre-sd = 0.439 ± 0.013, that is to say, the pre-sd star is the slightly fainter component, contributing ∼44% of the total K band continuum flux on average (the total flux includes a possible small contribution from the Be star disk). Some variability in the flux ratio is apparent with the minimum and maximum values from individual snapshots being 0.376 ± 0.018 and 0.540 ± 0.028, respectively, although most measurements cluster around the average value. The resulting parameters for the Be star and its disk are given in Tables 1 and A.2 and are discussed in Sect. 3.2.

3.2. The Be star and its Keplerian disk

The measurement of the Keplerian disk parameters is complicated by the fact that the spectro-interferometric signals in Brγ (T3PHI, DPHI, and |V|) are dominated by the photocenter displacement caused by the bright companion rather than by the typical signature of a resolved Keplerian disk (see Fig. A.3). It also needs to be considered that the changing (u, v) coverage of the data across the epochs means that they do not constrain the parameters in exactly the same way. Nevertheless, the overall consistency of the results across the epochs (listed in Table A.2) as well as the generally very good quality of the fit to the observed data (Fig. 1) shows that we are partly resolving the Brγ line-emitting region of the Keplerian disk, even though its size is clearly close to the limits of the angular resolution of the data. The resulting Keplerian disk parameters are given in Table 1, which lists the weighted average of the values across all epochs as well as their full range, including errors.

The resulting value of v_orb sin i_disk is dependent on the fixed value of UD_Be = 0.15 mas (corresponding to Be star radius R_Be ∼ 4.7 R_⊙ at the distance of 296 pc; see Sect. 4.2 and El-Badry & Quataert 2021). Relying on this value of UD_Be, we can estimate the critical rotation fraction of the Be star from comparison with a previously derived v_rot sin i, where v_rot is the surface rotational velocity of the Be star, while assuming that i = i_disk (i.e., that the rotational planes of the Be star and its disk are aligned). The spectroscopic analyses of both Bodensteiner et al. (2020a) and El-Badry & Quataert (2021) arrived at a consistent value of v_rot sin i ∼ 180 km s⁻¹. Adopting the larger uncertainty of ±20 km s⁻¹ from El-Badry & Quataert (2021), the resulting rotational parameter W = v_rot/v_orb is 0.744 ± 0.091 (corresponding to v_rot/v_crit = 0.81 ± 0.08 where v_crit is the critical rotation velocity). This value is typical for classical Be stars, which generally rotate close to but below their critical velocity (Rivinius et al. 2013).

The remaining parameters of the Keplerian disk are realistic, albeit not particularly precise, given their spread across the individual epochs. Still, they provide a consistency check against the orbital solution, as for a disk in a post mass-transfer binary, i_disk should be equal to the orbital inclination i, while PA_disk should be parallel to the line of nodes given by Ω. While the average PA_disk is indeed equal to Ω, the average i_disk is ∼2σ below i (Sect. 4, Table 1). Inspection of the individual epochs in which the resulting i_disk is significantly below i (these include the 2023 July 17 epoch shown in Fig. 1) reveals that the corresponding best-fit models appear to underestimate the central reversal of the Brγ line profile in each case. On the other hand, for the epochs when i_disk is consistent with i, the Brγ central reversal is deeper and more consistent with the observations. Thus, we conclude that the disagreement between the average i_disk and the orbital i is not significant and the overall results are consistent with the Be disk being aligned with the orbit. We note that fixing the values of i_disk at different values in the full range (Table 1) during the fitting process has negligible influence on the resulting astrometric positions and flux ratios.

4. Orbital solution for HR 6819

4.1. Astrometric orbital solution

To fit the astrometric positions in order to obtain the purely astrometric as well as the combined orbital solutions (Sect. 4.2), we used the IDL code orbfit-lib⁴ and the implemented Newton-Raphson method, which performs a formal orbital fit by calculating a first-order Taylor expansion for the equations of orbital motion to minimize the residuals (Schaefer et al. 2016). Unlike recent studies concerning HR 6819 (Bodensteiner et al. 2020a; Gies & Wang 2020; El-Badry & Quataert 2021), we considered the Be star as the primary component, as it is both (slightly) brighter in the K band on average (Sect. 3.1) and currently much more massive than the pre-sd star.

The 12 sets of relative astrometry determined above were used to derive astrometric orbital solutions, assuming both circular and eccentric orbits. In both solutions (as well as in the combined solution derived in Sect. 4.2; see Fig. 3), there are three astrometric points that do not fit the orbit within 3σ. One of these is the lower quality dataset from 2023 May 01, which may be affected by a miscalibration issue that had forced us to completely discard the |V| part of that dataset (Sect. 2). The other two are the datasets obtained in 2021 in split polarimetric mode, but it is unclear if any possible bias introduced by the polarimetric mode could be responsible for the lower quality of the orbital fit. To test how these three astrometric points influence the orbital solutions, we tried excluding them from the orbital fits, but the resulting parameters remained well within their respective 1σ errors in each individual solution.

Fig. 3. Astrometric orbit for the eccentric combined solution using Be star RVs (dashed line) with the Be star at the center (star symbol), the line of nodes (gray line), the ascending node (large plus sign), and the periastron position (small plus sign coinciding with the astrometric position on HJD = 60116.6). The relative astrometric positions of the pre-sd star (error ellipses corresponding to 3σ uncertainties) coinciding with the calculated positions on the orbit (x symbols) are shown alongside the average HJD of each observation. For the secondary axes in units of AU, a distance of 296 pc is assumed.

The eccentric orbital solution results in considerably smaller residuals, with χ_red.² decreasing from 10.6 to 5.5. Although the obtained eccentricity is small at e = 0.0306 ± 0.0071, it appears to be significant relative to its uncertainty at a 3σ level, as well as according to the statistical F-test comparing the circular and eccentric solutions. Nevertheless, with the exception of e and the associated epoch and longitude of the inferior conjunction or periastron, respectively, for the circular and eccentric solutions (T₀ and ω_pre-sd), the resulting orbital parameters are compatible within 1σ between the two solutions. The results are also all consistent with the previous works that analyzed FEROS spectra (Rivinius et al. 2020; Bodensteiner et al. 2020a; El-Badry & Quataert 2021), while the angular semimajor orbital axis (a″), and the position angle of the line of nodes (Ω) are both in a rough agreement with the preliminary values determined by Frost et al. (2022). The orbital inclination i ∼ 39° resulting from our solution lies in between the values estimated indirectly by Bodensteiner et al. (2020a) and El-Badry & Quataert (2021) and that of the preliminary astrometric fit of Frost et al. (2022). The resulting orbital elements for the purely astrometric solution are listed in Table 2.

4.2. Combined orbital solution

The secondary pre-sd star presents a multitude of narrow absorption lines in its optical spectrum, which enables accurate and precise measurements of its RV (e.g., Bodensteiner et al. 2020a, their Fig. 7). The pre-sd RVs were recently measured from archival FEROS spectra with different procedures on four separate occasions, leading to consistent results with the velocity semi-amplitude K_pre-sd ∼ 61 km s⁻¹ (Rivinius et al. 2020; Bodensteiner et al. 2020a; Gies & Wang 2020; El-Badry & Quataert 2021). For our combined orbital solution, we specifically utilized the RVs from Bodensteiner et al. (2020a), obtained by fitting several lines simultaneously with Gaussian profiles, and RVs from Gies & Wang (2020), obtained from cross-correlation with model spectra. The resulting orbital elements ended up being consistent within the respective 1σ errors when using either of these RV sets, with the only exception being the systemic velocity, which has no bearing on the physical parameters. The final combined solutions presented in Table 3 were obtained using the pre-sd RVs from Gies & Wang (2020) since they contain one additional measurement from an archival IUE spectrum, which helped increase the time baseline.

As for RV measurements of the Be star, the situation is complicated by the rotationally broadened and disk-emission-contaminated spectral lines, as well as the extreme mass ratio resulting in a very small semi-amplitude K_Be. From the recent works, only Gies & Wang (2020) published measurements for the individual FEROS epochs, obtained from fitting the broad wings of the Hα emission line, which originate from the very inner parts of the Be star disk, and therefore trace the radial motion of the Be star itself. Using the RV curve of the pre-sd star, Gies & Wang (2020) further refined the Be star RV star measurements by first subtracting the estimated photospheric profile of the pre-sd star from the composite Hα profiles at each epoch, resulting in the RV semi-amplitude K_Be = 3.9 ± 0.7. This value is consistent within 1σ with those obtained by Bodensteiner et al.(2020a, K_Be = 4.0 ± 0.8) and El-Badry & Quataert (2021, K_Be = 4.5 ± 1.0) using different methods. We used both K_Be and the individual Be star RVs from Gies & Wang (2020) in subsequent combined orbital solutions (see below).

For the combined orbital solution, the uncertainties of the astrometric and RV datasets were scaled in order to balance their contribution to the overall residuals. We considered both circular and eccentric orbits, similar to the purely astrometric solution described above. For the first set of combined solutions, we fitted the astrometry and pre-sd RVs, and used K_Be to calculate the masses and the distance (first and third columns in Table 3). In a slightly different approach, we also fitted the astrometry along with the RVs for both components (second and fourth columns in Table 3). These two approaches returned consistent results, as expected, with the only difference being the distribution of the resulting parameter uncertainties. Specifically, when fitting RVs of both components, the resulting K_pre-sd has a significantly smaller uncertainty than given by Gies & Wang (2020), while for most of the other parameters, the uncertainties are slightly higher. This then results in a better-constrained M_pre-sd, and slightly higher uncertainties for the other computed parameters. As for the purely astrometric solution, allowing for nonzero eccentricity results in overall lower residuals, with the eccentricity remaining small but significant at e ∼ 0.029. Both the masses and the distance are slightly higher in the eccentric solution than in the circular solution, but remain consistent within their respective errors. The astrometric orbit and RV curve corresponding to the eccentric solution using Be star RVs (fourth column in Table 3) are plotted in Figs. 3 and 4, respectively.

Fig. 4. RV curves for the combined eccentric solution using Be star RVs. The Be star and pre-sd star RV curves (solid and dashed curves, respectively) are overplotted with the RVs measured by Gies & Wang (2020) for the Be star and the pre-sd star (x symbols and empty circles, respectively), RVs measured from the Keplerian disk signature in Brγ in GRAVITY data (green stars), and RVs from He Iλ20587 and He Iλ21126 DPHI signatures in GRAVITY data (blue triangles and red diamonds, respectively). The errors of the RVs derived from GRAVITY data are dominated by the wavelength calibration uncertainty (Sect. 2). The lower panel shows the residuals in units of σ for the main RV dataset of Gies & Wang (2020).

The final results are listed in Table 3. Given the significance of the small but nonzero eccentricity, we took only the eccentric solutions into account for our final values and obtain d = 296.0 ± 8.0 pc, M_Be = 4.24 ± 0.31 M_⊙, M_pre-sd = 0.270 ± 0.056 M_⊙, and q = M_pre-sd/M_Be = 0.064 ± 0.012 (M_Be/M_pre-sd = 15.6 ± 3.0). If we considered more conservatively the full error ranges from both the circular and eccentric solutions, the resulting parameters would be d = 291 ± 13 pc, M_Be = 4.12 ± 0.43 M_⊙, M_pre-sd = 0.265 ± 0.061 M_⊙, and the same q as given above for the eccentric solutions only.

The distance resulting from the dynamical parallax of 296.0 ± 8.0 pc is significantly different from the Gaia DR3 distance of 368 ± 16 pc (calculated from the measured parallax while using photogeometric priors; Bailer-Jones et al. 2021). Without the knowledge of the Be star RVs, using the Gaia distance to compute K_Be would result in 18 km s⁻¹ or more for the four orbital solutions in Table 3, a value that is clearly overestimated by more than a factor of 4. This confirms that Gaia DR3 parallaxes can be strongly biased by binary motions such as those of HR 6819. HR 6819 was not included among non-single stars in the recent release of the Gaia DR3 catalog (Gaia Collaboration 2023b,a). However, this is not surprising given the stringent acceptance criteria for astrometric binary star processing in the DR3 release (Halbwachs et al. 2023), and both the relative error of the measured parallax and the orbital period of HR 6819 being rather low.

4.3. Radial velocities from GRAVITY data

Radial velocity measurements for the Be star disk, which should trace the Be star itself, are a direct by-product of the fitting of the Keplerian disk signature to the Brγ region described in Sect. 3. Unfortunately, the weakness of the Be disk signal (see Sect. 3.2) limits the precision of these RVs, despite the high wavelength calibration accuracy of the data. Nevertheless, the Brγ RVs do follow the Be star orbit determined from the combined solution, but with an rms error of ∼9 km s⁻¹, which is higher than the actual K_Be of ∼4 km s⁻¹.

We searched the entire K-band spectral range of the GRAVITY data for any additional spectral features that would enable us to obtain further information about the binary system. We found weak and narrow DPHI signatures in two helium lines He Iλ20587 (Fig. A.4) and He Iλ21126. While the He Iλ20587 line profiles show a hint of possible emission at < 5% of the continuum level, those of He Iλ21126 appear to show an absorption profile at a similar level. The DPHI signature is clearer for He Iλ21126, while at certain epochs it is barely discernible in He Iλ20587. Fitting of Gaussian profiles to the two lines shows that they originate from (the vicinity of) the pre-sd star rather than from the disk of the Be star. The weakness of the DPHI signal limits the accuracy and precision of the extracted RVs, but they still clearly trace the RV curve of the pre-sd star (Fig. 4). The rms residual is ∼19 km s⁻¹ for the weaker He Iλ20587, and ∼10 km s⁻¹ for He Iλ21126.

If the He Iλ20587 is indeed in emission (Fig. A.4), we might be detecting small amounts of circumstellar gas surrounding the pre-sd star (cf. Hanson et al. 1996), while the He Iλ21126 line probably originates in the photosphere of the pre-sd star. Given that the Be star is already surrounded by a well-developed self-ejected disk, parts of which might be reaching all the way to the secondary pre-sd star (Gies & Wang 2020), any circumstellar material around the secondary probably originates from accreted material from the Be disk. A possible presence of circum-companion gas was already suggested by Bodensteiner et al. (2020a) based on a small emission peak in Hα apparently following the orbit of the stripped star. The presence of an extended “bridge” region connecting the Be disk to the companion, as well as the possible presence of circum-companion structure formed by accreted material, is in line with recent predictions based on state-of-the-art hydrodynamical modeling (Rubio et al. 2025).

5. Discussion and conclusions

Owing to the high spatial and spectral resolution and the overall good quality of the VLTI/GRAVITY data, we were able to obtain the astrometric orbit of the binary HR 6819, composed of a newly formed classical Be star and a bloated, recently stripped pre-sd companion. Previously published RV measurements enabled us to obtain a combined three-dimensional orbital solution, leading to the first precise model- and geometrical parallax-independent dynamical masses for this rare type of astrophysical object, which represents a very short stage in the evolution of interacting intermediate-mass binaries. Our orbital analysis indicates that the distance derived from the Gaia DR3 parallax is biased due to the binary motions; the distance of 296.0 ± 8.0 pc derived here is significantly lower than the Gaia one (368 ± 16 pc; Bailer-Jones et al. 2021) and agrees well with the estimate of 310 ± 60 pc based on SED fitting (Rivinius et al. 2020). The orbit appears to be slightly eccentric, with e = 0.0289 ± 0.0058, and thus for the final values of the distance and the dynamical masses listed in this section and in the abstract we take only the eccentric solutions into account (Table 3).

The dynamical mass of the Be star (M_Be = 4.24 ± 0.31 M_⊙) is lower than what was estimated in recent works based on spectroscopy (Bodensteiner et al. 2020a; El-Badry & Quataert 2021). The Be star may be hotter than expected; while Bodensteiner et al. (2020a) assumed the Be star to have a spectral type of B2-3 V from the fitted T_eff, the dynamical mass is more consistent with B5 V to B7 V (following the mean parameters of main-sequence dwarf stars based on the literature and catalog survey⁵; Pecaut & Mamajek 2013, version 2022.04.16). The Be star is rotating near-critically and is surrounded by a self-ejected Keplerian disk that is aligned with the binary orbit. These findings are consistent with the expectation that mass and angular momentum transfer from the originally more massive component is responsible for the rapid rotation of the Be star and the subsequent formation of its disk. The dynamical mass of the pre-sd star is at the lower end of the estimated error ranges from Bodensteiner et al. (2020a) and El-Badry & Quataert (2021).

Of the known Be + (pre-)sdOB binaries, HR 6819 is only the second system with similarly well-constrained dynamical masses, after φ Per (Mourard et al. 2015), and only the ninth system with known dynamical masses when counting the seven parallax-dependent masses from Klement et al. (2022, 2024). In addition, dynamical masses are available for OGLE-LMC-T2CEP-211, which is a post-mass-transfer eclipsing system composed of a likely classical Be star and a Type II Cepheid (Pilecki et al. 2018); this binary system might be in a very short progenitor stage of HR 6819-like systems when the stripped star is passing through an instability strip while evolving toward lower temperatures shortly after the mass transfer.

HR 6819’s pre-sd star is expected to settle into a faint sdOB star within a few hundred thousand years (Bodensteiner et al. 2020a; El-Badry & Quataert 2021), at which point its spectrum will become dominated by the Be star. Compared to the known Be + sdOB systems with well-constrained parameters, HR 6819 has the most extreme mass ratio, the lowest stripped star mass, and the smallest orbital separation (Mourard et al. 2015; Klement et al. 2024). It is most similar to the first confirmed Be + subdwarf B-type (sdB) system, κ Dra, which has a slightly lower total mass but a ∼2 times higher mass ratio and a ∼30% larger orbital separation (Klement et al. 2022).

Given its similarity to κ Dra, HR 6819’s pre-sd star may become a relatively cool sdB star with a T_eff comparable to that of HR 6819’s Be star, which would make it extremely hard to detect with current observing techniques, even compared to the already challenging sample of known Be + sdOB binaries. Using the estimated pre-sd radius of ∼4.5 R_⊙ (Bodensteiner et al. 2020a; El-Badry & Quataert 2021), a future sdB radius of ≲0.7 R_⊙ (sdB radius for κ Dra; Klement et al. 2022), and a constant T_eff for simplicity, we can expect the pre-sd star’s total flux contribution to decrease to ≲1%. At the same angular separation, the higher contrast would put the stripped companion detection beyond the capabilities of VLTI/GRAVITY, and its spectral lines would not be visible in optical or far-UV spectra. Furthermore, the extreme mass ratio would make the measurement of the RV shifts of the Be star extremely elusive without prior knowledge of the orbit of the pre-sd companion, and the system would thus not be detectable as a spectroscopic binary. Only the low contrast resulting from the bloated evolutionary phase of the stripped star enabled both the angular resolution of its orbit with GRAVITY and the detection of the spectral lines of the future sdB star in the optical spectra. The angular resolution limit is less challenging for the CHARA array in the northern hemisphere, which has a longer maximum baseline than the VLTI (HR 6819 itself is not observable by CHARA).

It has been suggested by population synthesis studies that many, and possibly the majority of, Be stars are binary interaction products with faint evolved companions such as sdOBs or low-mass WDs (Pols et al. 1991; Shao & Li 2014, 2021). Thus, if the configuration of HR 6819 (once the stripped star contracts) is common among Be stars, the important observational implication is that most present-day Be stars will appear as single. The relatively small number (≈25) of confirmed descendants of HR 6819-like systems (Be + sdOB binaries; Wang et al. 2018; Klement et al. 2024) might thus represent only the most easily detectable tip of the iceberg of this population. The incidence of spectroscopic Be binaries with low-mass companions of unknown nature (e.g., Gies et al. 2000) as well as the prevalence of SED turndown among bright Be stars (see Sect. 1) are other observational clues suggesting that (mostly) invisible companions abound among Be stars.

As a by-product of the GRAVITY data analysis, we were able to measure the size of the Brγ line-emitting region of the Be disk, which appears to be smaller than the orbit (a″ ∼ 1.3 mas) at FWHM_Brγ,disk ∼ 0.5 mas. On the other hand, the disk kinematics derived from the Hα line suggest that the Hα line-emitting region probably extends all the way to the pre-sd companion (Gies & Wang 2020). This is not completely unexpected, as the Hα line emitting-region is generally the largest among all the emission lines. The angular size measurement of the Brγ region should prove useful for future modeling of Be disk hydrodynamics in close binary systems (Rubio et al. 2025), with a similar measurement already available for HR 2142 (Klement et al. 2024).

While Bodensteiner et al. (2020a) report a small emission peak in Hα probably originating from circum-companion gas accreted from the outer parts of the Be star disk, we were not able to detect a similar signature in the Brγ line. We did, however, detect weak signatures of the helium lines He Iλ20587 and He Iλ21126 that do originate from the pre-sd star, with He Iλ20587 possibly showing a small level of emission (Sect. 4.3, Fig. A.4). This emission may originate from circum-companion gas accreted from the Be star disk. Thanks to the highly accurate wavelength calibration of our GRAVITY data (< 0.01%; see Fig. A.1), we could measure RVs from the two helium lines, which were found to agree well with the pre-sd RV curve based on FEROS spectra. This confirms that high-spectral-resolution GRAVITY data can be used to measure precise RVs, which may be extremely difficult to obtain by other means (see also the case of HR 2142; Klement et al. 2024).

Acknowledgments

TS acknowledges support by the Israel Science Foundation (ISF) under grant number 0603225041. This work is based on observations collected at the European Southern Observatory under ESO programmes 107.22R6.002, 109.23FX.001, and 111.24PG.001.

Note added in proofs.In the abstract and later in the paper, we compared the distance derived from the dynamical parallax (296.0 ± 8.0 pc) to the Gaia DR3 distance given by Bailer-Jones et al. (2021, 368 ± 16 pc), which in this case is also almost identical to the distance calculated directly from the Gaia DR3 parallax. However, according to Maíz Apellániz (2022), the uncertainty of the Gaia distance is severely underestimated due to HR 6819 being a bright star. Using the calibration of Maíz Apellániz (2022) along with an OB-type prior for the Gaia parallax results in a distance of ${375}_{-58}^{+85}$, which is still much higher than the distance derived from the dynamical parallax, but considering the large uncertainty, it is only 1.4σ from our resulting value of 296 pc.

References

Appendix A: Detailed results from the fitting of GRAVITY data

The best-fit parameters are given for each epoch in Tables A.1 and A.2. The wavelength calibration accuracy for each epoch is plotted in Fig. A.1. Bootstrapping results from the geometrical model fitting for a representative epoch are shown in Fig. A.2. Averaged model for the Be star and its Keplerian disk as it would appear without a companion is compared to the dataset from the same epoch in Fig. A.3. Finally, an example of the fit to the He Iλ20587 line is shown in Fig. A.4.

Fig. A.1. Wavelength calibration accuracy, Δλ/λ, as a function of wavelength for each GRAVITY snapshot, as determined by PMOIRED fitting of telluric features. R is the spectral resolution and < err> is the error in the wavelength calibration averaged across the full K-band range. An accuracy of 0.001% in the K band corresponds to ∼3 km s−1.

Fig. A.2. Probability distributions resulting from PMOIRED bootstrapping giving the final values of the model parameters for the 2023 July 17 snapshot. The quantities and units are the same as those shown in Tables A.1 and A.2 except for RVBe, which is displayed here before barycentric correction and before including the wavelength calibration error (which ends up dominating the uncertainty). Orange represents the original fit to the data, and blue shows the bootstrapping results.

Fig. A.3. Same as Fig. 1, but showing the geometrical model with the Keplerian disk parameters averaged across all epochs (Table 1) and without any contribution from the companion, i.e., as it would be expected from a single Be star.

Fig. A.4. Same as Fig. 1, but showing the spectral region around He Iλ20587. The line profile (NFLUX) shows a possible small emission component, while DPHI shows a clear deviation from zero in the spectral line. There is no discernible signal in T3PHI or |V|. Only DPHI was used for the fitting of the helium line RVs.

All Tables

All Figures

Fig. 1. PMOIRED Keplerian geometrical model fit to the calibrated GRAVITY snapshot taken on 2023 July 17, representing the average quality of the model fit (χ2 = 3.1) as well as the typical quality of the data. The upper-left panel shows the (u, v) coverage, the lower-left panel the normalized Brγ line profile (NFLUX), and the remaining panels the interferometric observables closure phase (T3PHI), differential phase (DPHI), and the absolute visibility (|V|). The data (black with gray error bars) are overplotted with the model fit (red). The blue dots in the NFLUX and DPHI panels represent the continuum region used for normalization. The baselines and telescope triangles are given in the corner of each panel, and their colors correspond to the (u, v) coverage panel. For the layout of the VLTI telescope locations, see Fig. 10 in the VLTI manual (https://www.eso.org/sci/facilities/paranal/telescopes/vlti/documents/VLT-MAN-ESO-15000-4552_v115.pdf).

In the text

Fig. 2. PMOIRED synthetic image resulting from the Keplerian model fit to the 2023 July 17 snapshot (left) and the fluxes of the model components (right). The model components are the Be star primary (green) and Keplerian disk (orange) in the image center, and the pre-sd secondary (purple) offset in the image by ∼1.3 mas toward the northeast, with the sum of their fluxes normalized to 1.0 in the continuum (black).

In the text

Fig. 3. Astrometric orbit for the eccentric combined solution using Be star RVs (dashed line) with the Be star at the center (star symbol), the line of nodes (gray line), the ascending node (large plus sign), and the periastron position (small plus sign coinciding with the astrometric position on HJD = 60116.6). The relative astrometric positions of the pre-sd star (error ellipses corresponding to 3σ uncertainties) coinciding with the calculated positions on the orbit (x symbols) are shown alongside the average HJD of each observation. For the secondary axes in units of AU, a distance of 296 pc is assumed.

In the text

Fig. 4. RV curves for the combined eccentric solution using Be star RVs. The Be star and pre-sd star RV curves (solid and dashed curves, respectively) are overplotted with the RVs measured by Gies & Wang (2020) for the Be star and the pre-sd star (x symbols and empty circles, respectively), RVs measured from the Keplerian disk signature in Brγ in GRAVITY data (green stars), and RVs from He Iλ20587 and He Iλ21126 DPHI signatures in GRAVITY data (blue triangles and red diamonds, respectively). The errors of the RVs derived from GRAVITY data are dominated by the wavelength calibration uncertainty (Sect. 2). The lower panel shows the residuals in units of σ for the main RV dataset of Gies & Wang (2020).

In the text

	Fig. A.1. Wavelength calibration accuracy, Δλ/λ, as a function of wavelength for each GRAVITY snapshot, as determined by PMOIRED fitting of telluric features. R is the spectral resolution and < err> is the error in the wavelength calibration averaged across the full K-band range. An accuracy of 0.001% in the K band corresponds to ∼3 km s−1.
In the text

Fig. A.2. Probability distributions resulting from PMOIRED bootstrapping giving the final values of the model parameters for the 2023 July 17 snapshot. The quantities and units are the same as those shown in Tables A.1 and A.2 except for RVBe, which is displayed here before barycentric correction and before including the wavelength calibration error (which ends up dominating the uncertainty). Orange represents the original fit to the data, and blue shows the bootstrapping results.

In the text

	Fig. A.3. Same as Fig. 1, but showing the geometrical model with the Keplerian disk parameters averaged across all epochs (Table 1) and without any contribution from the companion, i.e., as it would be expected from a single Be star.
In the text

	Fig. A.4. Same as Fig. 1, but showing the spectral region around He Iλ20587. The line profile (NFLUX) shows a possible small emission component, while DPHI shows a clear deviation from zero in the spectral line. There is no discernible signal in T3PHI or |V|. Only DPHI was used for the fitting of the helium line RVs.
In the text