Interferometric Detections of sdO Companions Orbiting Three Classical Be Stars

Classical Be stars are possible products of close binary evolution, in which the mass donor becomes a hot, stripped O- or B-type subdwarf (sdO/sdB), and the mass gainer spins up and grows a disk to become a Be star. While several Be+sdO binaries have been identified, dynamical masses and other fundamental parameters are available only for a single Be+sdO system, limiting the confrontation with binary evolution models. In this work, we present direct interferometric detections of the sdO companions of three Be stars—28 Cyg, V2119 Cyg, and 60 Cyg—all of which were previously found in UV spectra. For two of the three Be+sdO systems, we present first orbits and preliminary dynamical masses of the components, revealing that one of them could be the first identified progenitor of a Be/X-ray binary with a neutron star companion. These results provide new sets of fundamental parameters that are crucially needed to establish the evolutionary status and origin of Be stars.


INTRODUCTION
Classical B-emission line stars (hereafter referred to as Be stars) are rapidly rotating main sequence B-type stars surrounded by purely gaseous, ionized, self-ejected disks, and they represent about 15-20% of all B-type stars in the local environment (although this number can increase dramatically at low metallicities, see Rivinius et al. 2013, for a review).A rotation rate close to the critical value is a defining property of Be stars (e.g., Zorec et al. 2016), and it is one of the necessary conditions for the formation of circumstellar decretion disks, from which the eponymous line emission originates.The origin of the rapid rotation among Be stars therefore holds a key to understanding the Be star class as a whole.
Binarity is a natural outcome of the star formation process.Components of close binary systems with periods shorter than a few years undergo periods of mass exchange that has a profound effect on the evolution (Podsiadlowski et al. 1992;de Mink et al. 2014) and the ultimate fate of these systems (e.g., Laplace et al. 2020).Close binarity is observed to be common among massive stars (see Duchêne & Kraus 2013, for a review).The fraction of close binaries among B-type stars in the Galaxy and the Large Magellanic Cloud appears to be at least 50% based on interferometric (Rizzuto et al. 2013) and spectroscopic (Abt et al. 1990;Sana et al. 2012;Kobulnicky et al. 2014;Dunstall et al. 2015;Villaseñor et al. 2021;Hastings et al. 2021) observations corrected for observational and selection biases.
Based on simulations of binary-star populations, de Mink et al. (2013) argued that rapidly rotating massive stars were spun up as a consequence of mass transfer or mergers, while de Mink et al. (2014) came to the conclusion that ∼ 30% of all main-sequence (MS) B-type stars are the products of binary interaction.Binary population synthesis studies specific to Be stars were also performed (e.g., Pols et al. 1991;van Bever & Vanbeveren 1997;Shao & Li 2014;Hastings et al. 2021;Shao & Li 2021), suggesting that if not all, at least a significant fraction is probably made up of binary interaction products.Most Be stars formed by mass transfer in a binary should then have hot, subluminous companions of the subdwarf OB-type (sdOB) type (similar to isolated helium stars) that can further evolve into white dwarfs (WD) or neutron stars (NS) to form Be/X-ray binaries (BeXRB).Although black hole (BH) are also a possibility, the number of Be+BH systems is expected to be small (e.g.Grudzinska et al. 2015).
On the other hand, even single stars could be created by close binary interaction, like the above mentioned case of fast-rotating merger products.The first identified case of a Be star being a merger product might be HD 93521, which is a rapidly rotating runaway late-O star with Hα emission from a faint disk, whose location far from the galactic disk is incompatible with the evolutionary age if only single star evolution is considered.Thus, there is strong evidence that HD 93521 is the result of a merger of two lower-mass stars, which had formed a runaway binary system following a supernova-induced disruption of a triple system or a dynamical ejection (Gies et al., submitted).Supernova explosions of massive stripped companions in close binaries can also lead to systemic disruption, causing the Be stars to become single runaway stars (Berger & Gies 2001).On the basis of the observed fraction of runaway stars in a large sample of Be stars (∼ 13%), Boubert & Evans (2018) conclude that all Be stars could be products of binary mass transfer.
Observationally, confirmed close binary systems comprise a non-negligible fraction of Be stars.Conspicuously, there is a lack of confirmed close MS companions to Be stars (Bodensteiner et al. 2020).As for stripped, evolved companions, the number of confirmed cases remains rather small.Other than the 200 or so known BeXRBs (which are however drawn from a large volume in the Local Group), there are only 15 spectroscopically confirmed cases of Be+sdO systems (no Be+sdB systems) and only a few more candidates (see below).Until recently, there were no confirmed cases of BeXRBs with WD (as opposed to NS or unconstrained) companions, despite the fact that detectable levels of low-luminosity X-ray emission had been expected (Meurs et al. 1992).The elusive X-rays compatible with Be+WD systems, in which the WD accretes the outer Be disk material, were recently detected in the Magellanic Clouds for for a total of six viable candidates (Kennea et al. 2021, and references therein).A related case to Be+WD binaries is the Bn star (rapid rotator without a disk) Regulus, which is orbited by a faint pre-WD star (Gies et al. 2020).The one reported case of a Be star being orbited by a quiescent BH (Casares et al. 2014) relied on possibly overestimated measurements of radial velocities (RVs), and high-resolution spectra resemble the known cases of Be+sdO systems (Rivinius, unpublished).
The low luminosity of sdO stars (compared to the Be primary) makes them difficult to detect spectroscopically in the visible (or at longer wavelengths), but their high T eff makes their contribution more easily detectable in the far ultraviolet (FUV, see below).In the visible, a few sdO companions have been found by detecting antiphased motions of the He II λ4686 line (in absorption or emission) that is associated with the hot sdO star (e.g., Poeckert 1981).By heating up parts of the Be disk, the sdO can also give rise to single-peaked emission line components (e.g., in He I λ6678) that follow the orbital motion but with larger velocity amplitudes than those of the companion (e.g., Rivinius & Štefl 2000).Other spectroscopic features that may indicate the presence of a hot companion include for instance short-lived phases of narrow shell-line absorption, and short-term variations of the (typically double-peaked) strong emission lines like Hα (Maintz et al. 2005).
Importantly, direct detection of sdO companions was shown to be possible with near-IR long-baseline interferometry (so far only for the case of ϕ Per, Mourard et al. 2015, see below), a technique capable of achieving angular resolution equivalent to a (diffraction limited) aperture of hundreds of meters.One should note that recently, an interferometric multiplicity survey among bright (m V ≤ 5.0) Be stars (classical and non-classical) was performed using data from the NPOI optical interferometer, but in that case the dynamical range was insufficient for the detection of companions contributing less than ∼ 4% in the visible (maximum ∆m of 3.5 at λ = 700 nm, Hutter et al. 2021).The study nevertheless provided the first direct detection of the helium/sdO star in the interacting binary υ Sgr (although the NPOI observations would be compatible also with an MS companion with a ≈ B5 spectral type).The helium star contributes ∼ 4% of the total visible flux, despite being the more massive component in this peculiar system.
Another possibility to detect companions of any kind is provided by their tidal influence on the Be star disk.With a close companion in orbit, the outer parts of the disk will be truncated, which translates into a lack of free-free emission at radio wavelengths (Waters et al. 1991;Klement et al. 2017).Analysis of the full sample of Be stars with published radio data led to the conclusion that all of them show the same characteristic turndown in the spectral energy distribution (SED, Klement et al. 2019).The caveats of this method are that (1) the SED turndown could be caused by a hitherto unidentified physical effect that is separate from binarity, and that (2) the nature of the possible companion remains unconstrained.On the other hand, one can estimate the size of the orbit of the possible companion based on the measured size of the Be disk (Klement et al. 2015).
The first identified Be+sdO binary ϕ Per (P = 127 d) was initially found from RV shifts of an emission component in the He II λ4686 line (Poeckert 1981).It was also the first one to be confirmed by cross-correlation techniques in FUV spectra, first using data from the International Ultraviolet Explorer short wavelength prime camera (IUE/SWP, Thaller et al. 1995), and later the Hubble Space Telescope (HST, Gies et al. 1998).The sdO companion of ϕ Per was also the first one -and until this study the only one -to be directly detected by long-baseline interferometry (Mourard et al. 2015), and a combined spectroscopic and astrometric solution revealed that the sdO star has a mass of 1.2 ± 0.2 M .The sdO was found to contribute ∼ 16% of the total flux in the FUV (Gies et al. 1998), but only ∼ 1.5% of the total flux in the H-band (including the contribution from the Be star disk), and ∼ 2.1% relative to the Be star photosphere alone (Mourard et al. 2015).
Several other Be+sdO systems were first detected (or suspected) from variable spectroscopic features before being confirmed in FUV spectroscopy.These are HR 2142 (Peters 1983), 59 Cyg (Rivinius & Štefl 2000), FY CMa (Rivinius et al. 2004;Štefl et al. 2005), and HD 55606 (Chojnowski et al. 2018).The systems (also) found in archival IUE/SWP spectra include FY CMa (Peters et al. 2008), 59 Cyg (Peters et al. 2013), HR 2142 (Peters et al. 2016), and 60 Cyg (Wang et al. 2017).Wang et al. (2018) listed 12 additional IUE candidates, nine of which were subsequently confirmed in new FUV spectra from the HST Imaging Spectrograph (HST/STIS), while the three Be+sdO candidates that were not confirmed are cases where the sdO companions could be obscured by the Be disks during certain orbital phases (Wang et al. 2021).A few additional Be+sdO candidates are suspected from the behavior of He lines in optical spectra but remain unconfirmed due to lack of FUV data: o Pup (Vanzi et al. 2012;Rivinius et al. 2012;Koubský et al. 2012), HD 161306 (Koubský et al. 2014), 7 Vul (Harmanec et al. 2020), andHD 35165 (Rivinius, unpublished).It is worth noting that no sdO companion around a late-type Be star has been found (Rivinius 2019), although the detection in the FUV spectra should be easier due to a lower FUV contribution from the Be star.This could be due to the fact that the search was optimized for hot sdO stars (T eff ∼ 45 kK), and failed to identify possible cooler sdB companions.
While there is still no unambiguously confirmed example of a Be star with a close MS companion (Baade 1992;Gies 2000;Bodensteiner et al. 2020), it should be mentioned that a significant number of Be stars have companions of unconfirmed nature.The SB1 systems 4 Her, 88 Her, Cap, ζ Tau, γ Cas, and Pleione are a few examples of bright Be stars with companions that could be faint, late-type MS stars, thus challenging the binary scenario as the origin of all Be stars.It is worth noting that γ Cas itself is a subject of intense study regarding its hard, variable X-ray emission with a luminosity of less that 1/30th of the values typical in BeXRBs but several times higher than in early B-type stars.These properties are shared by at least several more early-type Be stars, forming a subclass of γ Cas stars (Smith et al. 2016;Nazé & Motch 2018).The origin of the X-ray emission and the connection to a possible prevalence of close binarity among γ Cas stars remains unclear, as scenarios involving a WD, NS, sdO companion, as well as no companion at all, have been suggested (see, e.g., Langer et al. 2020, and references therein).As for Be stars with likely MS companions, the two with the shortest orbital periods appear to be δ Cen (period of ∼5.2 yr, companion spectral type between B4 V and A0 III, Meilland et al. 2008Meilland et al. , 2012) ) and the wider inner pair in the 2+2 quadruple system o And (period of ∼5.7 yr, companion spectral type estimated to be B6 III, Zhuchkov et al. 2010).The orbital periods in both of these are at least an order of magnitude larger that those of the known Be+sdO binaries.Another Be star with a possible MS companion is the highly eccentric δ Sco with an orbital period of ∼ 11 years, which was suggested to be a runaway triple system (Miroshnichenko et al. 2013).
In this paper, we report on the first results from a survey of binarity among Be stars performed at the Center for High Angular Resolution Astronomy (CHARA) Array, which offers sub-milliarcsec angular resolution and is capable of detecting faint companions contributing as little as 0.2% of the total flux in the near-IR (ten Brummelaar et al. 2005;Schaefer et al. 2020).The targets are introduced in Sect. 2 and the observations are described in Sect.3. The new binary detections and preliminary orbits are presented in Sects.4 and 5 before concluding with Sect.6.

TARGETS
Basic information on the observed Be stars is shown in Table 1, where the magnitudes m V and m H and the spectral types were adopted from the Simbad database1 (Wenger et al. 2000).The estimates for angular limb-darkened (LD) diameters of the Be star photospheres were taken from Wang et al. (2021) for 28 Cyg and V2119 Cyg and determined in the same way for 60 Cyg.The distances were taken from Gaia Early Data Release 32 (EDR3).Lastly, Table 1 lists the Be star physical radii R Be computed from the LD diameter and the distance.

28 Cyg
28 Cyg (HR 7708, HD 191610) currently possesses a developed if not particularly dense disk, as evidenced by amateur spectra available from the BeSS database3 (a total of 351 Hα spectra covering the period from 1995 to present, Neiner et al. 2011).At first sight, the double-peaked profile indicates that the circumstellar disk is seen under intermediate inclination.The equivalent width (EW) of the Hα line has mostly stayed between 0 and −10 Å, although there was a period of weakened emission (EW ∼ 0 Å) reflecting only a very weak disk between ∼2007 and ∼2012.Slettebak (1982) reported long-term variability for Hβ and Hγ.While Hβ was showing a symmetric, double-peaked emission in 1971 and 1979, it was in pure absorption in 1964.Interestingly, Hγ showed a 'very weak flanking emission around a central absorption core' (in 1971) and a 'sharp absorption core' (in 1979), suggesting a transient shell spectrum (i.e., disk seen close to edge-on inclination, cf. the case of 60 Cyg below).
28 Cyg has been extensively studied in the context of its short-term spectroscopic and photometric variations caused by non-radial pulsations (NRP), which are believed to be connected to the mass ejection mechanism that supplies material and angular momentum to the circumstellar disk (Rivinius et al. 2003;Baade et al. 2018).Variations attributable to a binary nature have not been reported until the orbiting sdO companion was found and its RVs measured by cross-correlation analysis of 46 archival IUE/SWP spectra (Wang et al. 2018).The sdO signature was detected (and its RVs measured) only in 25 out of the 46 exposures, and the resulting poor phase coverage prevented an estimation of the orbital period or any other orbital parameters.Although the companion was not confirmed by FUV spectroscopy from HST/STIS (three exposures) by Wang et al. (2021), the authors still consider it a likely Be+sdO system and suggest that an orbital phase-dependent obscuration of the companion by the Be star disk might be the reason why a detection was possible only at certain epochs.
V2119 Cyg was placed among Be+sdO binary candidates based on cross-correlation analysis of four IUE/SWP spectra (Wang et al. 2018), which led also to RV measurements of the sdO.The sdO companion has recently been confirmed in three HST/STIS spectra (Wang et al. 2021), making it possible to derive the following sdO parameters: T eff = 43500 K, v sin i < 15 km s −1 , and flux contribution of 4.7 ± 0.7% in the FUV (Wang et al. 2021, with three additional RVs measured for both the sdO and the Be star).Using the observed SED, the authors also derived estimates for the radius and luminosity of the sdO star: R = 0.52 ± 0.07 R , and log L(L ) = 2.94 +0.15  −0.23 .Using the measured RVs (a total of seven measurements for the sdO star), Wang et al. (2021) also presented a preliminary spectroscopic orbit with a period P = 60.286 ± 0.010 d and a velocity semi-amplitude K 2 = 75.5 ± 2.4 km s −1 .

60 Cyg
60 Cyg (HR 8053, HD 200310) shows similar disk properties and long-term variability to both 28 Cyg and V2119 Cyg.The 74 available BeSS spectra covering Hα (1995-present) show a diskless phase in the early 2000s, after which the emission gradually grew to EW of around −10 Å in 2013.After decreasing slightly to -7 Å in the few years after 2013, the EW grew to the present day value of around -12 Å. Slettebak (1982)  Based on spectroscopic analysis, Koubský et al. (2000) suggested that 60 Cyg is a binary with a 146.6 d period.They also derived a preliminary SB1 orbital solution based on RV variations of Hα emission-line wings, with the resulting Be star K 1 -velocity semiamplitude of 10.8 ± 0.1 km s −1 , and a mass function f (m) of 0.0191 M .
60 Cyg was confirmed to be a Be+sdO system from cross-correlation analysis of 23 IUE/SWP spectra (Wang et al. 2017).From the detected sdO signature in the FUV spectra, the authors estimated a mass ratio M sdO /M Be = 0.15 ± 0.02.Using the orbital solution and the Be star parameters derived by Koubský et al. (2000), they also derived the flux contribution of the sdO star to be 3.39±0.15% in the FUV, and its mass and radius to be around M = 1.7 M , and R = 0.48 R , respectively.Like 28 Cyg, 60 Cyg was also studied in the context of short period light and line-profile variations connected to NRP properties (Koubský et al. 2000).

OBSERVATIONS
The interferometric data of the three Be+sdO systems were collected at the CHARA Array (ten Brummelaar et al. 2005;Schaefer et al. 2020) as a part of an ongoing survey of binarity among Be stars that started in late 2020.The results presented here represent the first set of new binary detections and preliminary orbits from this survey.While we were able to derive preliminary orbits for two out of the three targets, we continue to monitor all three targets to arrive at fully constrained orbits that will lead to (more) precise estimates of dynamical masses.The CHARA Array consists of six 1-meter telescopes in a Y-shaped configuration (ten Brummelaar et al. 2005).Two telescopes reside on each arm, and they are labeled according to the direction of the corresponding arm: two East telescopes E1 and E2, and correspondingly S1, S2, and W1, W2, for South and West, respectively.The longest baseline of CHARA -B max ∼330 m -translates to an angular resolution of λ/2B max ∼ 0.5 milliarcsec (mas) in the case of the H-band (λ = 1.6 µm).
The Michigan InfraRed Combiner-eXeter (MIRC-X) installed at the CHARA Array is a highly-sensitive six-telescope interferometric image-plane beam combiner operating in the near-infrared J and H bands (Anugu et al. 2020) with spectral resolution R ∼ 50 or R ∼ 190.Combination of all six telescopes results in a total of 15 visibility-squared (VIS2), 20 closure phases (CP), and 20 triple amplitudes (T3AMP).MIRC-X has demonstrated H-band sensitivity down to a correlated magnitude of 8.2 and VIS2 and CP precision of 1% and 1 • , respectively.It has also demonstrated the ability to detect companions with contrast ratios of up to 500:1 (Anugu et al. 2020;Gallenne et al. 2015).A search in the Simbad database reveals that at the H-band magnitude limit of 8.2, MIRC-X can practically observe a sample of up to 383 Be stars.
The interferometric observing sequence consists of alternating between science targets and calibrators.The calibrators are stars with measured or well-estimated angular diameters, that are used to monitor the transfer function (i.e., the effect of atmospheric and instrumental bias on the interferometric observables) in between science observations, which is then used to calibrate raw science observables.Ideally, the calibrators are as close as possible to the science targets in terms of position on the sky, magnitude, and spectral type, while having a much smaller angular diameter.
In practice, it is usually not possible to meet all of these conditions, so the observer has to compromise between the different factors.
The log of MIRC-X observations for our three science targets and the corresponding calibrators is shown in Table 2.All three targets were observed on at least three separate occasions between 2021 May and 2021 September.We used the H-band (λ = 1.6 µm) and a low spectral resolution of R ∼ 50.The calibrators were selected using the SearchCal software (Bonneau et al. 2006;Bourgés et al. 2014) provided by the Jean-Marie Mariotti Center (JMMC4 ), and their basic characteristics including the uniform disk (UD) diameter are listed in Table 3.We used the official MIRC-X pipeline5 (v.1.3.5) to reduce the data (Anugu et al. 2020).
The interferometric field of view (FoV) is typically limited by the effect of bandwidth smearing, i.e., blurring of fringes with increasing spectral width of the individual wavelength channels.The bandwidth-smearing FoV is given by λ/(B∆λ), where λ is the observed wavelength, B is the baseline length, and ∆λ is the spectral width of an individual wavelength channel (Gallenne et al. 2015).For CHARA in the H-band with λ/∆λ ∼ 50, bandwidth-smearing limits the FoV to maximum offsets from the phase center of ∼ 50 mas.It is possible to correct analytically for the bandwidth smearing to a certain extent, but this was not needed for the present work.

BINARY DETECTIONS
The main interferometric observable VIS2 is a measure of the contrast of interference fringes formed from a combination of two separate beams.VIS2 measurements are sensitive to the geometrical size of the observed light source (along a projected baseline), and a resolved (i.e., larger than the given angular resolution) source will result in a decrease of the fringe contrast (visibility).The CP (phase of the bispectrum), on the other hand, results from a combination of three beams in a closed triangle, and is sensitive to departures from point symmetry in the interferometric FoV.In our case of searching for very faint companions that are offset from the phase center, CP is the more powerful observable, but VIS2 can provide important constraints on the geometrical extents and the flux ratio of the two components.
CP is also more robust to unstable atmospheric conditions compared to VIS2 as it is invariant to atmospheric phase perturbations.An additional and not yet widely used observable is T3AMP (amplitude of the bispectrum), also resulting from a combination of three telescopes in a closed triangle, although this observable can be strongly affected by the atmosphere, similarly to VIS2 (see, e.g., Glindemann 2011; Buscher & Longair 2015, for an introduction to interferometric observables).
A first look at the calibrated VIS2 reveals that the three targets are mostly unresolved (VIS2 > 0.8) even on the longest baselines, corresponding to UD diameters of 0.25 mas (cf.Table 1 with photospheric LD estimates).This prevents us from determining reliable angular extents of the Be stars, the angular orientation and tilt of their disks, or the relative flux contribution of the stellar photospheres and the circumstellar disks.Inspection of the calibrated CPs reveals that the signal remains within around ±5 • , ruling out strong asymmetries or the presence of a low contrast companion.As will be described below, the small CP signals unambiguously point towards the presence of very high contrast companions at separations of a few mas.
We used the dedicated open-source code Companion Analysis and Non-Detection in Interferometric Data6 (CAN-DID, Gallenne et al. 2015) to perform a search for binary companions around the three observed Be stars.As opposed to traditional grid search resulting in a χ 2 map, CANDID uses a more rigorous approach of computing a 2D grid of fits, i.e., a multi-parameter fit is performed for each grid point which represents the initial guess for the position of the companion.With the primary component at the phase center, the parameters are the companion position (∆α, ∆δ), the observed flux ratio of the companion and the primary f O = f (sdO)/f (Be), where the term f (Be) contains both the flux from the stellar photosphere and the surrounding disk (if present), and the UD diameter of the primary.In the fitting process, it is required that multiple starting points reach the same minimum, thus ensuring that the parameter space is fully explored and that the minima are not local.Based on the traveling distance (to convergence) of the fits in the parameter space, CANDID determines a posteriori the optimal grid step, certifying that a global minimum can be found by repeating the fitting process with the correct step size.CANDID also allows one to estimate the detection confidence level as a multiple of the standard deviation σ (assuming that the data uncertainties are uncorrelated) with a maximum level of σ = 8.0.Bootstrap sampling is used to derive the final best-fit parameter values and their uncertainties and to derive the companion position in the form of an error ellipse.CANDID also enables robust determination of a 3σ detection limit on the flux ratio of undetected companions from a given dataset, which is done on the basis of 'injecting' fake companions at different positions and with different flux ratios into the data.
We searched for companions within the bandwidth-smearing FoV, i.e., up to a maximum separation of ∼ 50 mas between the components.In a first step, we performed binary fits for the companion position and the flux ratio using CP only, and with the UD diameter of the primary Be star fixed at the estimated photospheric angular diameter listed in Table 1.This solution resulted in unambiguous and highly significant (> 8.0σ) detections of a high contrast companion at separations < 10 mas in all but one data set, namely 28 Cyg on MJD = 59398.350(with 3.8σ).However, despite having larger uncertainty, the less significant detection agrees well with the measurement from the night that followed (MJD = 59399.337),so that the orbital point should still be valid.
In subsequent steps, we included the VIS2 and T3AMP observables in the binary fits, while introducing the UD diameter of the Be star as a free parameter in the fit.In a few cases (see below), apparent small-scale miscalibration issues in VIS2 and T3AMP (probably caused by variable atmospheric conditions) affected the binary detections.First, the 60 Cyg observation on MJD = 59439.395resulted in unusable VIS2 and T3AMP, but the CP signal remained unaffected.Since the binary detection resulting from CP only is consistent with the other detections for this object, we still consider this to be a reliable orbital point.Second, although including VIS2 confirms the above mentioned weaker detection (despite slightly lowering its significance) for 28 Cyg (MJD = 59398.350),including both VIS2 and in this case a rather noisy T3AMP results in a non-detection.Third, for another two measurements (MJD = 59364.455and MJD = 59479.215),the detection confidence decreased slightly (to 7.2σ and 5.0σ) after the inclusion of VIS2 (and remained similar after including also T3AMP), but the locations of the χ 2 minima remained unambiguous and the flux ratios are consistent with the other measurements.
Overall, all but the two above mentioned measurements (MJD = 59439.396and MJD = 59398.360)resulted in consistent companion locations to within 3σ when using CP only, CP+VIS2, and CP+VIS2+T3AMP.As the final values, we adopt the solution from a simultaneous fit to CP+VIS2, as it results in the most consistent flux ratio across the different epochs (standard deviation of < 0.1%).The one exception is the above mentioned date MJD = 59439.396,where we use the result from fitting only the CP.The numbers are listed in Table 4, where ρ is the angular separation, PA is the position angle measured from North to East, and ∆RA and ∆DEC are the companion offsets.In the remaining columns, the uncertainties are given as 1σ error ellipses with semi-major and semi-minor axes σ-a and σ-b and the PA of the semi-major axis σ-PA (measured from North to East).The uncertainty includes instrumental 0.5% calibration error of the absolute wavelength scale (Anugu et al. 2020).An example plot with the final residuals is shown in Fig. A1.The weighted averages of the observed flux ratios f O = f (sdO)/f (Be) (H-band) are listed in Table 5.Here, the term f (Be) contains a possible contribution from the circumstellar disk of the Be star, and the f O values thus possibly underestimate the flux ratio of the component photospheres only.
Next, we turn to confronting the measured H-band flux ratios with previous results from FUV spectroscopy.The available photospheric flux ratios f O (FUV) were measured through the amplitude of the cross-correlation function (CCF) of the combined and model subdwarf spectra (Wang et al. 2017(Wang et al. , 2021)).Since the FUV studies also led to preliminary estimates of the effective temperatures of the sdO companions, we can use the FUV flux ratio and model spectra of the components to determine the photospheric flux ratio as a function of wavelength and, in particular, find the calculated photospheric flux ratio f C in the H-band.The model spectra used were from the TLUSTY OSTAR2002 grid (Lanz & Hubeny 2003), and from the TLUSTY BSTAR2006 grid (Lanz & Hubeny 2007), for the sdO and Be components, respectively.
Compared to V2119 Cyg and 60 Cyg, the case of 28 Cyg is more problematic, as the spectral signature of the sdO component was absent in recent HST spectra (Wang et al. 2021), even though it was previously detected in 25 of 46 archival spectra from IUE (Wang et al. 2018).This means that no estimate of the effective temperature, and only an upper limit on the flux ratio, could be derived (Wang et al. 2021).Therefore, we assumed the temperature of the sdO component to be T eff = 45 kK (as found in other cases) and we estimated the flux ratio f O in the FUV by comparing the average cross-correlation function peak height for 28 Cyg with those for other cases with a measured f O (see the column with heading S/N in Table 2 of Wang et al. 2018).
The measured and expected flux ratios are summarized in Table 5, which gives the flux ratio f O (FUV) estimated from the FUV spectral analysis, the calculated flux ratio f C in the H-band, and the weighted average of the observed flux ratios f O in the H-band.Comparison of the model spectra with near-IR photometric measurements (2MASS, Cutri et al. 2003) revealed no significant contribution of the circumstellar disks to the total H-band flux at the time of the observations (between the calendar years 1998 and 2000).Inspection of the available BeSS spectra shows that at that time, 28 Cyg had Hα emission similar to the present day, V2119 Cyg had a stronger Hα emission, and 60 Cyg, on the other hand, was found in an almost completely emission-less (and therefore diskless) state.Thus, direct comparison of f C and f O in the H-band should be possible at least for 28 Cyg and V2119 Cyg, while for 60 Cyg, the present-day disk contribution might lower the value of f O , thus underestimating the photospheric flux ratio.The agreement indeed turns out to be excellent for 28 Cyg and V2119 Cyg, but in the case of 60 Cyg, f C (based on IUE spectra only) is ∼ 1.6% compared to f O of ∼ 1.2%, likely revealing the above mentioned contribution of the disk to the H-band flux.c No VIS2 or T3AMP available, the fit was done using CP only.These numbers are similar to what was found for ϕ Per (Sect.1), whose developed disk also contributes significantly to the H-band flux.Thus, our H-band results offer an important verification of the FUV detection methods and confirm the sdO nature of the companions.Lastly, we derived upper limits on the flux contribution of additional undetected companions with the CANDID 'injection' method.In this case, CANDID is used to remove analytically the previously detected companion from the observables before performing the injection, thus obtaining an unbiased detection level map for any potential additional companion that remains undetected.The average values of the limiting companion flux ratios in the H-band f lim were calculated within a separation ρ of 25 and 50 mas from the primary component using the best-quality datasets for each target.Table 6 shows the limiting flux ratio estimates for our three targets.The changing companion positions on different dates are compatible with the expectation for binaries with orbital periods of one to a few months, and we clearly detect orbital motion in datasets separated by at least two nights (Table 4).However, given the still very few orbital points, we can only make preliminary estimates of the orbits.The orbits described below were obtained with the Orbit Fitting Library7 (Schaefer et al. 2006(Schaefer et al. , 2016)).
The full set of parameters for a binary orbit includes the orbital period P , epoch of periastron T (or T RVmax for circular orbits), eccentricity e, semimajor axis a in angular units and a in physical units, orbital inclination i, longitude of the ascending node Ω, longitude of periastron ω, velocity semiamplitudes K 1 and K 2 , and the systemic velocity γ.One can only get the apparent a from astrometric orbits, and the projected a 1 sin i and a 2 sin i (where a = a 1 + a 2 ) only from spectroscopic orbits.Upon combination of astrometric orbits (to derive i) and double-lined spectroscopic orbits (to derive a sin i), one can calculate the total mass using the 3rd Kepler law, and the individual masses using the ratios q = M 1 /M 2 = K 2 /K 1 = a 2 /a 1 .When only a single-lined spectroscopic orbit (and therefore only either a 1 sin i or a 2 sin i) is available, we need to know the distance to obtain a from a , after which we can calculate the total mass, as well as the individual masses from the knowledge of the ratio a 1 /a 2 (see, e.g., Batten 1973 for an introduction to binary orbits).

28 Cyg
For 28 Cyg, no estimates for the orbit or orbital period are available in the literature.Although we have three interferometric detections, two of those were taken on subsequent nights and the results are compatible within the 3σ error bars.The resulting coverage of the astrometric positions as well as that of the RVs measured by Wang et al. (2018) is thus insufficient for period estimation.Assuming that the measurement of angular separation of 6.44 mas on MJD = 59399.337(Table 4) is close to the orbital semimajor axis a , and that the total mass of the system is 10 M (based on the total mass of the ϕ Per system with similar properties, Mourard et al. 2015), the orbital period P resulting from the third Kepler law and Gaia EDR3 distance is 246 d.In Fig. 1, we plot a possible circular astrometric orbit and RV curve with P fixed at this value.The resulting inclination for this orbit is ∼ 118 • , indicating a retrograde (counter-clockwise) orbital motion.Fixing P at lower values results in an increasingly edge-on inclination.Given that spectroscopic shell signatures were detected in the past (Sect.2.1), and assuming that the orbital and disk planes are aligned, a more inclined orbit is a possibility (cf. the case of 60 Cyg described below).The flux ratios are consistent across the three detections, not showing any strong variability that would be analogous to the FUV results, where the companion was detected only at certain epochs.

V2119 Cyg
For V2119 Cyg, Wang et al. (2021) derived a preliminary circular SB1 orbit from the RVs of the sdO companion measured in FUV spectra.We used the available orbital parameters as the starting point in a simultaneous fit to the four astrometric positions and the six available RVs of the sdO companion.The result reveals a circular orbit seen at intermediate inclination with an excellent fit to both the astrometry and the RVs (Fig. 2).The resulting parameters are summarized in Table 7, where they are compared to those of Wang et al. (2021).
With the knowledge of the distance to the system, we can compute the physical size of the orbit as well as the masses of the components.The resulting mass of the sdO companion is 1.62 ± 0.28 M which is slightly above the Chandrasekhar limit (∼ 1.4 M ).Since isolated helium stars with masses above 1.6 M will ultimately explode as core-collapse supernovae (Woosley 2019), V2119 Cyg could be the first identified progenitor of a BeXRB with a NS companion (Podsiadlowski et al. 2004).

60 Cyg
For 60 Cyg, a preliminary single-lined spectroscopic orbit was given by Koubský et al. (2000).We started with a fit to the astrometric positions while keeping the spectroscopic parameters fixed according to the available solution to obtain first guesses of the other orbital parameters.Since we were not able to find a satisfactory circular orbit, we allowed the eccentricity to vary, and arrived at a slightly eccentric orbit (e ∼ 0.2) seen at high inclination (i > 80 • ).Subsequently, we added the 47 RVs of the Be primary (measured from Hα wings by Koubský et al. 2000), for which we arbitrarily adopted errors of ±5 km/s.The scatter of the RV measured is large, showing the difficulty of measuring precise RVs of Be stars.The resulting orbit and RV curve are plotted in Fig. 3.As with V2119 Cyg, we can calculate the physical size of the orbit and preliminary dynamical masses of both components (with knowledge of the distance).The final parameters are listed in Table 8, where they are compared to results from previous studies, namely P , T RVmax , K Be , γ and M Be derived by Koubský et al. (2000), and q and M sdO derived by Wang et al. (2017), revealing a good agreement except for the component masses, which are lower according to our solution.We note that our preliminary Be star dynamical mass of 7.3 ± 1.1 also appears to be low for a B1-type star according to available spectroscopic calibrations based on dynamical masses (e.g.Harmanec 1988), although this might be simply due to a misclassification of the spectral type.
Assuming that the disk orientation is aligned with the binary orbit, the high orbital inclination of 83 • seems to be at odds with the non-shell spectral appearance.Rather, the available spectra in the BeSS database suggest an inclination around 70-75 • (cf.Fig. 1 in Hanuschik et al. 1996).However, signatures of shell absorption have been observed in the past (Sect.2.3), which indicates that either (1) the disk is seen close to edge-on, but shell absorption is for some a Fixed for a circular orbit.1).Left: Preliminary (circular) orbit of the sdO companion resulting from a joint astrometric and RV fit, with the Be star at the center.Right: RV curve of the secondary with RV measurements from Wang et al. (2018Wang et al. ( , 2021)).
reason not generally apparent, or that (2) the orientation of the disk oscillates in time similarly to the cases of Pleione, γ Cas, and 59 Cyg, although much less dramatically (e.g.Tanaka et al. 2007).Considering the first option, the shell absorption could be diluted by emission from an increased amount of disk material at zero RV.This would be similar to the case of the shell star ζ Tau, where a density wave propagates through its disk (manifesting itself via periodic asymmetries of the Hα emission), and as the low density region passes in front of the star (at zero RV), the shell absorption is diluted by emission from the high density region behind the star (also at zero RV, Carciofi et al. 2009).In Be stars without strong density waves in the disk, like 60 Cyg (Koubský et al. 2000), a similar effect would occur if the vertical height of the disk was enhanced, so that there is enough disk material next to the star (or rather below or above) that is not projected onto the stellar disk (cf.Fig. 5 in Hanuschik 1995).
Considering the non-zero eccentricity of the resulting orbit, there is at least one other known case of a slightly eccentric Be+sdO system, for which it was suggested that it might be caused by the dynamical influence of a third wide companion (Peters et al. 2013).Thus, we speculate that also 60 Cyg could have a wide companion, although it is not apparent at the separations and flux ratios probed by CHARA/MIRC-X data (cf.Table 6).

CONCLUSIONS
For the first time, sdO companions of the classical Be stars 28 Cyg, 60 Cyg, and V2119 Cyg were directly detected with optical interferometry, confirming the Be+sdO nature in the case of 28 Cyg through an excellent agreement with the expected near-IR flux ratio (Table 5), and bringing the total of confirmed Be+sdO binaries to 16 (22 when including strong candidates).This represents only the second set of interferometric detections of (classical) Be+sdO binaries, the first one being the single case of ϕ Per (Mourard et al. 2015).These results demonstrate the usefulness of optical interferometry for binary detection and the derivation of astrometric orbits at angular resolution inaccessible to any single-telescope observing technique.
The positions of the detected companions are consistent with orbital periods ranging from one to a few months (Table 4).Preliminary orbital solutions were derived from the new astrometry and available RV measurements for V2119 Cyg and 60 Cyg (Tables 7 and 8), while for 28 Cyg we were only able to suggest a possible orbit due to insufficient number of interferometric measurements.Preliminary dynamical masses for V2119 Cyg and 60 Cyg indicate that in the case of V2119 Cyg, the sdO companion might be massive enough (1.62 ± 0.28) to explode in a core-collapse supernova, and therefore that V2119 Cyg itself could be the first identified progenitor of a BeXRB with a NS component.For 60 Cyg, the resulting orbit is slightly eccentric, making it possibly the second Be+sdO system whose orbit is not circular (after 59 Cyg).Ongoing mapping of the orbits and monitoring in high-resolution spectroscopy are expected to lead to refinement of the orbits and the dynamical masses in the near future, providing crucial and robust test cases for binary evolution models.1).Left: Preliminary (slightly eccentric) orbit of the sdO companion resulting from a joint astrometric and RV fit, with the Be star at the center.Right: RV curve of the primary with RV measurements from Koubský et al. (2000).
The measured flux ratios in the H-band were compared with those derived from previous FUV results and model spectra, revealing a good agreement (Table 5).The interferometric H-band results are thus an important verification of the FUV methods used to detect the companions, and to derive their flux ratios as well as fundamental parameters such as T eff .

Figure 1 .
Figure 1.Relative astrometric orbit and RV curve for the sdO companion of 28 Cyg.Left: Possible (circular) orbit (dashed) of the sdO companion (with fixed P = 246 d) with the Be star at the center.The ellipses are the measured astrometric positions with 3σ uncertainty and the corresponding points on the computed orbit are marked with x signs.The line of nodes is shown in grey.Right: Possible RV curve of the secondary (solid line), corrected for the systemic velocity.RV measurements (black circles with grey error bars) are from Wang et al. (2018).The lower panel shows the (O − C) residuals.

Figure 2 .
Figure2.Relative astrometric orbit and RV curve for the sdO companion of V2119 Cyg (symbols are the same as in Fig.1).Left: Preliminary (circular) orbit of the sdO companion resulting from a joint astrometric and RV fit, with the Be star at the center.Right: RV curve of the secondary with RV measurements fromWang et al. (2018Wang et al. ( , 2021)).

Figure 3 .
Figure3.Relative astrometric orbit and RV curve for the sdO companion of 60 Cyg (symbols are the same as in Fig.1).Left: Preliminary (slightly eccentric) orbit of the sdO companion resulting from a joint astrometric and RV fit, with the Be star at the center.Right: RV curve of the primary with RV measurements fromKoubský et al. (2000).

Table 1 .
Sample Be stars Calculated from the angular diameter and the distance.
reported occasional 'sharp dark cores' in Balmer lines from 1954 to 1976, decreasing Hα emission between 1964 and 1976, and a double-peaked emission in Hβ in 1979-80.As with 28 Cyg, this indicates a transient shell-like spectrum and the possibility of a close to edge-on inclination if the circumstellar disk.

Table 2 .
Log of MIRC-X observations a Total number of individual measurements.

Table 3 .
Calibrators Target UD diam.(H-band) a mH Spectral type

Table 5 .
Measured (O)and calculated (C) sdO to Be flux ratios

Table 6 .
Detection limits on additional components

Table 8 .
Parameters for 60 Cyg a The literature value is TRVmax b Fixed for a circular orbit.