© The Authors 2024

1. Introduction

Massive stars with initial masses of greater than 8 M_⊙ live eventful lives that often conclude with powerful explosions and leave behind a compact object such as a neutron star or stellar-mass black hole (Woosley et al. 2002). They enrich their surroundings with gas, dust, nucleosynthesis products, and ionizing radiation, strongly influencing the local and galactic environment throughout their lives and after their death (Haiman & Loeb 1997; Nomoto et al. 2013; Hopkins et al. 2014). Most massive stars reside in binary or multiple star systems and undergo interaction with a companion in the form of mass transfer or merger (Sana et al. 2012, 2014; Moe & Di Stefano 2017). The multiplicity of massive stars therefore plays an important role in their evolution (Langer 2012; Marchant & Bodensteiner 2024). Understanding the multiplicity properties of massive stars at different evolutionary stages can provide strong constraints on the nature of past interactions and predictions of their future evolution.

Wolf-Rayet stars (WRs) represent one of the final stages of stars with initial masses of greater than ∼20 M_⊙, directly preceding black hole formation (Conti 1976; Higgins et al. 2021; Pauli et al. 2022; Aguilera-Dena et al. 2023). Observationally, WRs are characterized by strong, broad emission lines in their spectra, originating from strong optically thick winds. A large fraction of WRs are hydrogen-free, and these are referred to as classical WRs or cWRs. Wolf-Rayet stars are also classified as nitrogen-rich (WN), carbon-rich (WC) or, very rarely, oxygen-rich (WO) based on the dominating spectral lines, and in particular the emission line strengths and their ratios (Smith 1968; Crowther 2007). WN stars are often subdivided into early- and late-type WNs, namely WNE (earlier than WN6) and WNL (WN6 and later), respectively. Another subtype of WRs are the hydrogen-rich WN stars. The suffix ‘h’ was added to the spectral classification to denote the presence of hydrogen in WR spectra (Smith et al. 1996). WNh stars tend to be similar to late-type WN stars, and show a continuity with very massive O stars in terms of spectral properties. Unlike cWRs, a large fraction of WNh stars are likely very massive main sequence stars that have strong winds that give them a WR-like spectroscopic appearance (Crowther 2007). In this work, we refer to and treat WNh stars separately from cWRs.

Classical WRs are evolved, post-main sequence objects (typically core-He burning) that have been stripped of their hydrogen envelope. It has long been thought that this stripping happens either through stellar winds (Conti 1976; Abbott & Conti 1987; Smith 2014) or via interaction with a close companion (Paczyński 1967; Vanbeveren et al. 1998). In the single-star channel, also referred to as the modified Conti scenario, the formation and evolution of WRs is thought to follow the path O → red supergiant (RSG)/ luminous blue variable (LBV) → WN → WC (→ WO in rare cases). Despite large uncertainties, single-star evolution models generally cannot reproduce the luminosity distribution of observed WR populations, suggesting that the binary channel plays a significant role in producing WRs (Schootemeijer & Langer 2018; Shenar et al. 2019; Hamann et al. 2019). The contribution of each of the two scenarios to the Galactic WR population remains to be determined (Neugent & Massey 2014; Shenar et al. 2019, 2020; Schootemeijer et al. 2024). Understanding the multiplicity of massive stars in the Wolf-Rayet phase can help anchor the role of binary evolution in the formation of WRs.

Over the years, many studies have addressed the multiplicity of WRs in the Galaxy (Vanbeveren & Conti 1980; van der Hucht 2001). The latter studied a sample of 227 WRs and arrived at an observed binary fraction of 0.4 based on photometry, spectroscopy, visual binaries, and dilution of emission lines due to a companion. The heterogeneity of the instruments and methods used was a limiting factor for appropriate bias correction and subsequent derivation of the intrinsic binary or companion fraction. Presently, the Galactic WR Catalog v1.30¹ (originally Appendix 1 in Rosslowe & Crowther 2015) contains 679 stars. Taking a more systematic approach than previous studies, Dsilva et al. (2020, 2022, 2023) considered a homogeneous, magnitude-limited sample of 39 Galactic WRs (27 WNs and 12 WCs, where WNs also include H-rich WNs) to establish their intrinsic binary fraction. In particular, these authors used spectroscopic data to determine radial velocity variation in WRs and employed well-understood bias correction to derive their intrinsic binary fraction.

WN binaries were found to have a period distribution that dominates at short periods (P < 100 d), similar to that of the OB stars (Dsilva et al. 2022, 2023). This is not compatible with the intermediate RSG/LBV phase, wherein objects can have radii large enough to initiate interactions (Mahy et al. 2022), while common-envelope transition toward short-period WNs can be excluded from theoretical considerations. Dsilva et al. (2020) found an observed binary fraction of 0.58 ± 0.14 for WCs, while Dsilva et al. (2023) found 0.41 ± 0.09 for WNs (1σ errors from binomial statistics). Based on a Bayesian inference analysis of the additional radial velocity (RV) noise, Dsilva et al. (2023) predicted that the period distributions for WN and WC stars should differ significantly, with WNs skewed toward shorter periods and WCs skewed toward longer periods. These authors also predicted a higher intrinsic binary fraction for WCs ($0.{96}_{-0.22}^{+0.04}$) as compared to WNs ($0.{52}_{-0.12}^{+0.14}$), although the difference between the two could originate from small number statistics. These results are based on modeling RV variations seen in the WR spectral lines, which are heavily affected by wind variability. This confines the ability of spectroscopy to only detect binaries with short periods, high mass ratios, and high inclinations. The binary detection probability drops at large orbital periods (P > 1 year) and the Bayesian framework developed for the bias correction at long periods is subject to assumptions as to the shape of the period and mass-ratio distributions (Dsilva et al. 2022, 2023). With optical long-baseline interferometry, it is possible to detect binaries with longer periods that are usually beyond the reach of spectroscopy. Interferometry has previously been used to study hot dust in WRs (Rajagopal 2010), the wind collisions in WR+O star binaries (Lamberts et al. 2017), and to determine the orbits of such binaries (Monnier et al. 2011; Richardson et al. 2021, 2024). Interferometry can also serve as a powerful and less time-consuming complementary method for searching for new binary companions at long periods (e.g., Rajagopal 2010; Richardson et al. 2016), better constraining the multiplicity and resolving the potential period distribution discrepancy between WNs and WCs.

The evolutionary connection between WNs and WCs is thought to follow the path WNL → WNE → WC. This evolution is a direct consequence of the strong winds of WRs, which continuously strip the outer layers at high mass-loss rates, revealing the products of hydrogen burning (WNs) and helium burning (WCs) (Gamow 1943; Lamers et al. 1991). The orbital period of WR binaries changes during this phase due to a combination of mass loss and angular momentum loss via winds, the former increasing the period by a factor of ∼2 at most (Dsilva et al. 2023). Such an increase is not sufficient to explain the predicted period distribution discrepancy of WNs and WCs, indicating that WC binaries should preferably evolve from WN binaries with longer orbital periods, and that the short-period WN binaries should avoid evolving into WC binaries altogether. While the latter is a problem that is best tackled with spectroscopy, the binarity at longer periods can be better studied with interferometry.

To investigate the multiplicity of WRs at long periods, we turned to the Very Large Telescope Interferometer (VLTI) operated by the European Southern Observatory (ESO). Equipped with several infrared instruments, the VLTI is capable of resolving objects with angular separations in the range of ∼1 to 200 milliarcseconds (mas). For a typical ∼2 kpc distance for Galactic WRs, this corresponds to orbital periods in the range of 10² − 10⁵ days, extending the period parameter space well beyond that of spectroscopic measurements (< 10³ days). A census of WR multiplicity at long periods is an important step toward understanding not only the origin of such systems, but also their evolution into potential compact object binaries.

In this paper, we present our results on the multiplicity of 39 Galactic WRs based on interferometric observations with the VLTI, and discuss their implications on the formation and evolution of WRs. The paper is organized as follows: Sect. 2 describes the observations and data reduction of our WR sample, and Sect. 3 details our data analysis approach. In Sect. 4 we describe the results of our analysis, including detected binaries as well as other features found in the dataset. Section 5 focuses on the multiplicity of our sample as well as the limitations of our observations. Finally, we present the possible interpretations of our results and our conclusions in Sect. 6.

2. Sample, observations, and data reduction

We chose VLTI/GRAVITY and its Auxiliary Telescopes (ATs) to investigate the multiplicity of WRs, owing to the unique spectro-interferometric capabilities of the instrument in the K-band (1.95–2.45 μm). We selected all WRs with K <  9 visible during the observing run, and excluded targets that were previously proposed/observed with GRAVITY. An exception to this was WR 48, which has been observed with GRAVITY before and is a candidate triple system. This added up to a total of 49 WRs – 23 WCs and 26 WNs+WNhs.

A total of 39 WRs out of the requested 49 were observed. Observations for the 10 remaining targets were hampered owing to technical issues or poor observing conditions. Additionally, of the WRs excluded from our proposal, archival GRAVITY data for 4 WRs were also included in our sample for this work (1 WC and 3 WN+WNh stars). All data were obtained in the medium (R ∼ 500) or high (R ∼ 4000) resolution mode of GRAVITY, depending on the K-band magnitudes of the WRs (high resolution for brighter targets).

We performed the data reduction and visibility calibration using the standard GRAVITY Pipeline v1.6.0 (Lapeyrere et al. 2014). The calibrator stars selected for program 111.24JN were FIII stars within a 2° radius and ΔK <  2, ideally brighter than the WR. For the targets from the archival program 109.23CN, the calibrators chosen were early-type A dwarf stars. A total of 4 targets in our sample were not reduced successfully because of bad data quality. To summarize, a total of 14 targets were excluded from our analysis due to issues with observations/data reduction, comprising of 9 WC stars and 5 WN+WNh stars. Nevertheless, our sample still covers all spectral types of WRs uniformly, as listed in Table D.1.

Our final reduced sample was hence comprised of 39 WRs. This also included 4 stars analyzed by Dsilva et al. (2020, 2022, 2023) as part of their study. A full list of targets along with their corresponding calibration stars, times of observation and spectral resolution used are summarized in Table A.1. The data quality varied significantly across the sample. The last column of Table D.1 lists a qualitative data evaluation for every target as good, average or poor based on visual inspection and signal-to-noise ratio (S/N) considerations. Although not quantified, this classification is meant to convey the overall quality of observations.

3. Data analysis

The GRAVITY instrument provides a K-band spectrum of the target as well as interferometric observables such as visibility amplitude (|V|), closure phases (T3PHI) and differential phases (DPHI). We used all the available data types to qualitatively assess every system. A brief description of this qualitative approach is discussed in Appendix B. This was followed by geometric modeling to look for any resolved companions or other features. PMOIRED² (Mérand 2022) was our primary analysis tool for modeling the spectro-interferometric data. PMOIRED enables parametric modeling where one can add individual components such as circular and elliptical uniform disks, Gaussian point sources, and many more complex structures, and fit such models to the data given some prior parameter values or grids. We also used the telluric correction feature in the module (based on a grid of atmospheric spectra generated from molecfit, Smette et al. 2015) to obtain the continuum normalized flux (NFLUX) for every target.

WRs typically have strong optically thick winds that can extend up to several stellar radii. In the K-band, the Potsdam Wolf-Rayet Models (PoWR) predict the extent to be few tens of solar radii. An example for WR 78 is given in Fig. 1, illustrating that in the K-band the continuum opacity already starts to extend to slightly larger radii than at optical wavelengths. Lines can be optically thick to much larger radii, but as evident from Fig. 1 this effect is not prominent in the K-band. Considering that WRs in our sample are typically ∼2 kpc away, this extent corresponds to an angular size much less than 1 mas, and thus mostly unresolvable for the VLTI. Consequently, our initial expectation was for all WRs to be approximately point sources. Nevertheless, we used a Gaussian component to model the WR and left the full width at half maximum (fwhm) as a free parameter in case the WR winds were more extended than expected. To search for potential resolved companions to WRs, we assumed main sequence stars which can be treated as point sources.

Fig. 1. Radius in solar radii where the optical depth is equal to unity based on a PoWR atmosphere model for WR 78. The blue curve only takes the continuum opacity into account, while the light-red curve also accounts for the line opacities.

Some WRs are also known dust producers, which is apparent through the additional flux contribution of dust seen in infrared photometry. The spatial scales for such dusty regions around WRs can reach several thousands of astronomical units (au), making them resolvable with the VLTI. In our sample of 39 WRs, only one star is a known dust producer (WR 113; WC8d+O8-9IV). Nevertheless, such emission in the context of GRAVITY data can be factored in as a fully resolved flux spanning the whole field of view. Such an emission component was included agnostically for all targets and ignored in case it did not contribute significantly.

The K-band spectra for WRs are characterized by emission lines from certain elements corresponding to their spectral type, which are thoroughly discussed in literature (Figer et al. 1997; Crowther et al. 2006; Rosslowe & Crowther 2018; Clark et al. 2018). For WN stars the Brγ line (often blended with He I and He II) is almost always present in emission, along with a component of the He I + N III 2.11 μm and the He II 2.189 μm features. WC stars, on the other hand, typically show the C IV 2.07–2.08 μm and C III 2.114 μm (blended with He I 2.113 μm) features. Some WCs also show a C IV emission line around 2.43 μm (Clark et al. 2018). Additionally, WRs show a further variety of spectral lines varying from star to star depending on the composition and velocities of their winds, sometimes even showing different spectral features for the same spectral type (Clark et al. 2021). Nevertheless, we mainly focus on the multiplicity aspect of WRs in our sample in this work, and only comment on the spectral lines wherever necessary.

As described in Appendix B and shown in Fig. B.1, the presence of spectral lines may or may not manifest itself in |V| data depending on what components the system is made up of. Moreover, fitting |V| and T3PHI data alone is typically sufficient to model the presence (or absence) of a resolved companion. We therefore used |V|-T3PHI data for companion search assuming flat spectra, and only invoked NFLUX and DPHI for targets showing spectral features in any of |V|, T3PHI or DPHI data.

We converged on the final model for every WR using reduced chi-squared (${\chi }_{\mathrm{red}}^2$) minimization in PMOIRED. We examined the data visually first to look for any obvious signatures of a companion. This was followed by finding best-fit models without and with a companion to look for any significant improvement in the ${\chi }_{\mathrm{red}}^2$ value of the latter. When including NFLUX and DPHI data, we only modeled the most prominent spectral features that were sufficient to explain the data. A side-effect of this was much higher values for ${\chi }_{\mathrm{red}}^2$ due to unmodeled features, particularly for targets observed in the high resolution mode. This primarily concerns the binary detections and is further described in Appendix C. Furthermore, the residuals from our model fits were not always flat at zero, often showing some baseline-dependent signals. This could arise due to the complex spatial properties of the target beyond our modeling capability, or due to imperfect visibility calibration. We only modeled companion(s) and fully resolved flux components; exploring more detailed structure in GRAVITY snapshot data is beyond the scope of this work.

The errors obtained by PMOIRED do not account for systematics and are overly optimistic. Consequently, once we converged on the final model, we performed bootstrapping on the model to get realistic error-bars on our best-fit parameters (Lachaume et al. 2019). As an additional note, all error-bars noted in this paper, obtained from model fitting or otherwise, are 1σ errors.

4. Results from interferometry

In this section we discuss the results of our interferometric data analysis along with more details about modeling for the different categories of WRs found in our GRAVITY sample.

4.1. Unresolved stars

Single stars are expected to be mostly unresolved point sources in GRAVITY data due to their small angular sizes. Similarly, stars in binaries with angular separations smaller than the VLTI resolution limit (∼1 mas) can also resemble single stars in interferometric data. We combined these two categories into one, termed “unresolved stars” since they were indistinguishable and hence treated identically in our analysis. Here, unresolved refers to the lack of a resolved companion in |V|-T3PHI. In terms of interferometric observables, these stars were expected to show mostly flat, close to unity visibility amplitudes and zero closure and differential phases. A total of 29 WRs in our sample fall under this broad category. A good example of this is WR 22 (Fig. 2), which is a WN7h + O9III-V spectroscopic binary (period ∼80 days) showing point-source-like properties. The best-fit model was a point source component alone.

Fig. 2. |V|-T3PHI data for the WN7h + O9III-V spectroscopic binary system WR 22. The different colors in T3PHI represent four combinations of 3 ATs each, while those in |V| represent the six baselines of 2 ATs each. The top panels show the best fit model (gray) along with the data, while the bottom panels show the residuals. The flat-at-zero T3PHI and flat-at-unity |V| data are typical of a single unresolved point source.

We expected all unresolved stars to show point-source-like features just as WR 22. Yet several of them show a significant downward shift in |V| at all baselines, indicative of a second, fully resolved component. WR 18, a single WN4 star, is a good representative for such cases (Fig. 3). It was modeled using a point source and a fully resolved component uniformly spanning the whole field of view. Such a component could not be constrained in size and its exact spatial extent remains uncertain.

Fig. 3. |V|-T3PHI data for the WN4 star WR 18 (similar to Fig. 2). The flat-at-zero T3PHI and flat but vertically offset |V| data are best explained by a combination of an unresolved point source and a fully resolved component.

4.2. Resolved wide binaries in visibility amplitude and closure phase

In our GRAVITY data, the presence of a wide companion with an angular separation of about 1–200 mas was likely to produce a strong signal in |V| and T3PHI, down to a contrast of ΔK ∼ 5. In particular, a non-zero T3PHI implies asymmetry in the target and motivates the search for a companion. For the central WR, we used a model similar to the unresolved stars described in Sect. 4.1. To search for a companion, we added a second point source with the ΔEast or ΔE and ΔNorth or ΔN coordinates (with respect to the primary) searched over a grid and flux ratio set as a free parameter (Gallenne et al. 2015). The search was done in two steps – an initial larger and coarser grid to get an approximate solution, followed by a smaller finer grid to constrain the companion well. Details of the grid search are summarized in Appendix C.

Only four targets were found to have potential wide companions. It is important to note that our observations represent a single epoch for each target, and follow-up observations would be necessary to confirm whether a companion is indeed gravitationally bound to the central WR. That being said, the probability of chance alignment of a random star in a 200 mas radius around the WR is very low. For a quantitative estimate, we followed the method described for Eq. (8) from Sana et al. (2014) to estimate the probability of spurious association for companions within 200 mas and brighter than ΔK <  5, finding a ≲0.1% probability of chance alignment. Below we describe the four systems in more detail:

WR 89: This object is classified as a WN8h+abs single star (Smith et al. 1996), where abs stands for absorption lines of unknown origin, though they could in principle also be intrinsic to the WR. It was identified as a thermal radio source by Montes et al. (2009) and a strong X-ray source by Nazé et al. (2013). The latter supported the possibility of WR 89 being a colliding wind binary based on its properties, but could not confirm it using archival optical spectra. This hinted at a potential long-period binary (years to decades) and necessitated long-term spectroscopic monitoring to confirm binarity.

In GRAVITY data, WR 89 shows clear deviation from point-source-like behavior (see Fig. 4). We performed a grid search for a potential binary companion and found a global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = 0.02  ±  0.08 mas and ΔN = −5.18  ±  0.10 mas relative to the central primary. The flux of the companion is 0.24  ±  0.02 times that of the WR star. Based on a Gaia DR3 parallax of 0.27  ±  0.02 mas, the projected separation of the companion is about 19  ±  1 au.

Fig. 4. |V|-T3PHI data for the WN8h+abs star WR 89 (similar to Fig. 2). The T3PHI and |V| data deviate drastically from symmetric/point-source-like behavior, and are best modeled by a wide binary system with two distinct point sources.

WR 115. This source is classified as a WN6-w single star (Hamann et al. 2019), where the -w stands for weak lines. It also shows spectral line variability potentially caused by corotating interaction regions (St-Louis et al. 2009). A speckle-imaging companion search by Shara et al. (2022) identified a candidate early B dwarf companion with a projected separation of 0.20 arcseconds (200 mas) and position angle of 320 degrees East of North. Whether or not this supposed companion is bound to the central WR could only be confirmed with follow-up astrometric observations.

In GRAVITY data, we see clear signs of a companion in |V|, T3PHI as well as DPHI. We performed a grid search for a potential binary companion and found a global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = −124.80  ±  0.04 mas and ΔN = 154.93  ±  0.02 mas relative to the central primary. In terms of projected separation and position angle, this translates to ∼199 mas and 321 degrees East of North respectively. This is largely consistent with Shara et al. (2022), although follow-up observations in the future would be necessary to confirm the orbit. The flux of the companion is 0.05  ±  0.01 times that of the WR star. The Gaia parallax for WR 115 is 0.58  ±  0.46 mas, and cannot be used to constrain the physical projected separation due to its large error-bar.

WR 48. Also known as θ Mus A, WR 48 was initially classified as a spectroscopic binary with a 19-day orbital period (Moffat & Seggewiss 1977; Hill et al. 2002). This was based on a composite spectrum (WR+O at least) that showed features of a WC6 star moving at the orbital period in short-period binary, and features of an O star that were virtually stationary, possibly unrelated to the binary. Hartkopf et al. (1999) performed a speckle interferometric survey of WRs, in which they detected an O supergiant wide companion to WR 48, although with some uncertainty over their relative positions.

In GRAVITY data, we detected this potential outer companion at ΔE = 1.59  ±  0.03 mas and ΔN = −12.56  ±  0.06 mas relative to the central binary. It has a flux 1.27  ±  0.10 times that of the central binary. Based on the Gaia parallax of 0.46  ±  0.07 mas, the projected separation is 28  ±  4 au. One more epoch of speckle interferometry data (Hartkopf et al. 2012) and two more epochs of archival GRAVITY data are available for θ Mus A. A more comprehensive analysis of the object using all available data will be covered in Gosset et al. (in prep.).

WR 93. This source was classified as a WC7 + O7-9 spectroscopic binary by van der Hucht et al. (1988) based on its spectra, although no orbital period has been determined to date. It belongs to one of the embedded clusters in the H II complex NGC 9357 (Lima et al. 2014).

We detected a strong binary signal in GRAVITY data for WR 93. We performed the routine grid search for a companion, finding the best fit for ΔE = 5.60 ± 0.03 mas and ΔN = 13.19 ± 0.03 mas. The flux of the companion is 0.20 ± 0.01 times that of the primary, and is likely the O7-9 companion known from spectroscopy. Based on the Gaia parallax of 0.49  ±  0.03 mas for WR 93, the projected separation is about 29  ±  2 au.

More details about the companion grid search along with plots for these WRs are presented in Appendix C.

4.3. Resolved spectroscopic binaries in differential phase

A total of 12 stars in our sample have been previously classified as confirmed spectroscopic binaries with OB companions and measured periods ranging from 6–80 days. For a typical parallax of 2 kpc, this translates to separations < 0.5 mas which is too small to be resolved with the VLTI when using |V| and T3PHI. However, since WRs have emission lines in the K-band and OB stars do not, we could potentially get a non-zero signal in the differential phase data (DPHI) at the emission line wavelengths. This is because DPHI is sensitive to spatial asymmetries in the wavelength regime – in this case, owing to the presence of emission lines in WRs and lack thereof in OB stars. This has been shown in the fitting of other binary systems. For example, in Frost et al. (2022) the properties of a decretion disk around a Be star were determined from fits to DPHI, despite its size being below the VLTI resolution limit and there being no drop in the |V| to fit.

This method of companion detection was applied to the whole sample regardless of their current classification. For the purpose of modeling in this case, the NFLUX and DPHI data were also included and models were updated to include spectral lines. A minimum binary separation of ∼(λ/B)×e_DPHI/360° radian is necessary for a detection in DPHI, where λ is the wavelength of observation, B is the baseline and e_DPHI is the error level in DPHI data. For our GRAVITY observations assuming e_DPHI ∼ 1°, the minimum separation comes out to be ≳0.01 mas. To model such a companion, a second point source with a flat spectrum was added with a grid search region spanning 2  ×  2 mas centered on the primary star. Only five stars out of the sample of 39 were found to have a DPHI signal corresponding to any spectral line(s). WR 98, a WN8/WC7 + O8-9 spectroscopic binary, is an example of such a binary detection as shown in Fig. 5. Three broad emission lines at 2.11 μm (He I + N III blend), 2.16 μm (Brγ + He I + He II blend) and 2.19 μm (He II), typical for WN stars as discussed in Sect. 3, can be seen in NFLUX as well as DPHI. The features in DPHI indicate an asymmetry in all spectral lines, hinting towards a close binary as discussed above.

Fig. 5. GRAVITY data for WR 98, a WN8/WC7 + O8-9 spectroscopic binary, zoomed in at 2.10–2.22 μm along with the best-fit model shown in red. The top left panel shows the u-v coverage for the observation. The remaining panels, from left to right, show NFLUX, T3PHI, DPHI and |V| respectively, color-coded according to baselines. The data are zoomed in to a wavelength range of around 2.10–2.22 μm, with three prominent emission lines present at 2.11, 2.16 and 2.19 μm. The DPHI data (along with some of |V| data) show slight but clear features corresponding to the spectral lines. The best-fit model for such a system is a very close binary consisting of the line-emitting WR and a flat-spectrum companion, along with a fully resolved component.

The remaining four stars were WR 9, WR 42, WR 47 and WR 79, which are discussed in more detail in Appendix D. All five binaries detected using this method are known spectroscopic binaries. The remaining seven binaries were not detected likely due to their smaller projected separation at the time of observation, higher flux contrast, and/or poor S/N of GRAVITY data. Nevertheless, we demonstrate the potential of this method for short-period binary detection in WRs, in addition to the typical |V|-T3PHI detections for wide binaries.

4.4. Resolved line-emitting regions

WRs are characterized by strong optically thick winds that dominate their spectra, resulting in broad emission lines. The different emission lines are produced by a multitude of processes, including resonance scattering, collisional excitation, (dielectronic) recombination, and continuum fluorescence, with complex interdependencies (see, e.g., Hillier 2015), for a more detailed discussion and examples). The dominating process for the origin of a specific line further depends on the detailed, non-local thermodynamic equilibrium (nonLTE) population numbers and can change for different conditions and thus different WR subtypes. In the dense winds of WR stars, ionization changes in the outer wind are common, giving rise to the appearance of emission lines of many different ionization stages in the spectrum. For the purpose of this study, the recombination of helium from the fully ionized He III to He II is particularly important as it is responsible for the prominent appearance of He I lines in the K-band.

The spectro-interferometric capabilities of VLTI/GRAVITY offer a unique insight into extended line-emitting regions. For a target with multiple components contributing to the flux, spectroscopy alone can be limited in distinguishing the flux coming from different components. Combining with interferometry, however, it is possible to do so. Taking a WR star for example, if the total flux is a sum of the continuum originating from the (unresolved) photosphere and emission lines originating from the extended (resolved) winds, the |V| data will reveal their presence as two different components. Such spectral disentangling can provide a physical measurement of the extent of line-emitting regions, that has only been estimated from radiation hydrodynamical models before.

Figure 6 shows the GRAVITY data for the WN7h star WR 78. The spectrum shows WN-like broad emission lines at 2.11 μm (He I + N III blend), 2.16 μm (Brγ + He I + He II blend) and 2.19 μm (He II) μm. The T3PHI and DPHI data are flat at zero, implying no point asymmetry in the object. The overall level of |V| is essentially unity, characteristic of an unresolved point source. However, two spectral features can be seen corresponding to the lines at 2.11 μm and 2.16 μm, while such a feature is absent for 2.19 μm. The likely explanation for such |V| data is the presence of two components in the system – (i) an unresolved point source responsible for the continuum and the line 2.19 μm (uniform disk with fixed radius = 0.001 mas); (ii) a slightly resolved component only emitting the lines 2.11 and 2.16 μm (gaussian with fwhm = free parameter). On fitting such a model to the data, we obtain an fwhm of 0.51 ± 0.02 mas. WR 78 has a parallax of 0.56 ± 0.03 mas in Gaia DR3, based on which, the fwhm or “size” of the resolved component is 0.89 ± 0.05 au or 191 ± 11 R_⊙.

Fig. 6. GRAVITY data for WR 78, a single WN7h star, along with the best-fit model shown in red. The top left panel shows the u-v coverage for the observation. The remaining panels, from left to right, show NFLUX, T3PHI, DPHI and |V| respectively, color-coded according to baselines. The data are zoomed in to a wavelength range of around 2.09–2.21 μm, with three strong emission lines present around 2.11, 2.16 and 2.19 μm. The |V| data show clear features corresponding to the first two spectral lines, decreasing in prominence from longer to shorter baselines. The best-fit model for such a system is an unresolved point source emitting the 2.19 μm spectral line and a slightly resolved Gaussian (fwhm ≈ 0.51 ± 0.02 mas) emitting the 2.11 and 2.16 μm spectral lines.

Using a PoWR atmosphere model of WR 78, the line formation strengths per radius ξ (defined in Hillier 1987) are shown in Fig. 7 for the observed lines in Fig. 6. The values are normalized, as the absolute values of ξ carry little physical meaning, but the figure illustrates where the line formation is highest in the stellar wind. It is evident that the peak of the formation for the He I lines occurs just a bit outside of the region where He recombines from He III to He II as discussed above. Notably, the main formation further happens well within 1 au for all the considered lines, including those we resolve in our data. The reason for this is likely twofold: firstly, imperfections in the model for WR 78 might slightly underestimate the necessary stellar radius by incorrectly estimating the optically thick wind onset region and the velocity stratification (see also Lefever et al. 2023). This is not unlikely as our model is not a dynamically consistent one (Sander et al. 2020, 2023) and could affect the radial scale of the line-formation calculations. Secondly, even if the peak of line formation is unresolved, the emission strength can still be strong enough farther out in the wind to be resolved with GRAVITY, which seems to be the case for He I lines.

Fig. 7. Line formation strengths ξ (Hillier 1987) at 2.11 (top), 2.16 (middle), and 2.19 μm per radius for the PoWR atmosphere model used in Fig. 1. The ξ values are normalized for better line-to-line comparison. The radii where most line formation occurs Rmax are highlighted by the dashed black lines.

In our sample of 39 stars, a total of four WRs (all of them late WN/WNh stars) were found to have resolved line emitting regions. It is important to note here that our observations were taken in the snapshot mode of GRAVITY and higher S/N observations with more uv-coverage can provide an even more detailed look into the line-forming regions. Such observations could be undertaken and discussed in a future study.

The wide variety of interferometric data for WRs in our sample could further motivate separate studies for each category described above. A target-by-target summary for all stars in our sample is provided in Appendix D. WR 113 being the only dust producing WR in our sample shows a distinct signature, also discussed in detail in Appendix D. In the following sections, we focus on the multiplicity analysis and subsequent interpretation for our WR sample.

5. Multiplicity fraction

5.1. Observed multiplicity fraction

The multiplicity of WRs as noted in literature is largely based on spectroscopic and photometric observations. These methods are less likely to detect companions to WRs that are lower in mass, inclination and/or in wider orbits due to their smaller effect on the Doppler motion of the WR. With our GRAVITY observations, we have taken a complementary approach of detecting the wider companions with interferometry. By doing so, we opened up a new parameter space covering two orders of magnitude in projected separations (three orders in orbital period) down to K-band flux contrast as low as ∼1%.

Table 1 summarizes the multiplicity status for our sample classified according to their spectral classes – WNh (H-rich WN), WNL (late-type WN), WNE (early-type WN) and WC. Here, WNL refers to WRs with spectral types WN6 or later, while WNE refers to earlier spectral types.

For every spectral class, we have reported the number of stars that were – (a) classified as single; (b) confirmed double-lined spectroscopic binaries (SB2) but the companion was unresolved in GRAVITY data; (c) confirmed SB2s with the secondary resolved in GRAVITY DPHI data; (d) found to have a wide companion resolved in GRAVITY |V|-T3PHI data. For the special case of WR 48, which is likely a triple system, it is both a confirmed short-period binary and has a wide companion, and is therefore included twice. The observed multiplicity fraction, which is calculated as the total number of WRs that have companions divided by the sample size, is listed in the last column, along with error-bars derived using binomial statistics.

The observed multiplicity fraction for WRs altogether was found to be 0.38 ± 0.08. Based on spectral classes, the observed fraction for WNL (0.50 ± 0.20), WNE (0.60 ± 0.22) and WC (0.40 ± 0.13) were found to be consistent with each other within the error-bars. WNh stars, however, turned out to have a relatively lower multiplicity fraction (0.23 ± 0.12) than their classical counterparts. In the next subsection we discuss the important physical limitations of our interferometric survey before proceeding to Sect. 6, where we interpret the multiplicity fractions from Table 1 in an evolutionary context.

5.2. Determining physical parameters for detected companions and detection limits

Our magnitude-limited GRAVITY observations of 39 WRs form a homogeneous sample in terms of the parameter space explored, making it possible to uniformly assess the overall properties and limitations of our detection methods. A large fraction of WRs in our sample was unresolved in GRAVITY data, including some known close binaries as well as stars classified as single. For all such targets, we calculated a 3σ detection limit for a wide companion using PMOIRED. This was done by injecting a point source randomly within 200 mas around the central primary and estimating the flux required for a 3σ detection (as described in Absil et al. 2011). A graphical example of this for WR 18 is shown in Fig. 8. The detection limit was expressed as magnitude contrast with respect to the central WR. It tends to vary across different positions with respect to the central star, resulting in an approximate Gaussian distribution as seen in the second plot of Fig. 8. As a conservative estimate for the limiting magnitude contrast, we picked the 99 percentile value at the lower end of the distribution. For example, in case WR 18, the 99 percentile value is 5.057 magnitude.

Fig. 8. Companion detection limit plots for WR 18. The first plot shows a randomly created dense grid about 200 mas in radius around the central primary star that was used to inject a second point source and calculate the flux ratio needed for a 3σ detection. The scale for flux ratio is given in magnitude with reference to the primary. The second plot shows a histogram of limiting magnitudes over all the grid points. Lastly, the third plot shows detection limit as a function of radial distance, along with a rolling average in red.

This contrast in K-band magnitudes or corresponding flux ratios could be interpreted physically, based on certain assumptions and provided we have a good understanding of the fundamental stellar parameters of the stars. This could be done for both detection limits as well as detected wide companions to some extent, as described in the following text.

Hamann et al. (2019) and Sander et al. (2019) have reported fundamental stellar parameters for a large number of Galactic WN (including WNh) and WC stars respectively using updated Gaia DR2 distances. This includes 24 WRs from our sample, most of which are classified as single stars. For WR binaries in which the companion contributes significantly to the spectrum, making it composite, the determination of stellar parameters was not possible. We used all information available from the two studies for our WR sample, with the primary goal of obtaining corresponding absolute magnitudes for detected companions or detection limits. A brief description is given below.

Hamann et al. (2019) have listed the following stellar parameters for WN and WNh stars – the narrow-band color excess (E_b − v), relative visibility (R_V) associated with the applied reddening law, and absolute v-band magnitude (M_v). Firstly, we converted the narrow-band color excess to broadband color excess using the equation from Turner (1982):

$$\begin{array}{c}\hfill E_{B-V}=1.21\times E_{b-v}.\end{array}$$(1)

We used this to calculate the V-band extinction and subsequently the K-band extinction following the equations from Fitzpatrick (1999):

$$\begin{array}{cc}\hfill A_V& =R_V\times E_{B-V},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill A_K& =0.12\times A_V.\hfill \end{array}$$(3)

Finally, we employed the apparent K-band magnitudes listed on the Galactic Wolf Rayet Catalog v1.30 and Gaia DR3 distances to compute absolute K-band magnitudes using the equation

$$\begin{array}{c}\hfill M_K=m_K-A_K-5\times {\text{log}}_{10}\left(\frac{d}{10}\right),\end{array}$$(4)

where d is the distance to the star expressed in parsec. We performed a similar exercise based on Sander et al. (2019) for the WC stars in our sample. In this case, however, the R_V parameter was not reported. The uncertainty in R_V only propagates into A_K with a factor of 0.12. Consequently, we assumed the standard value of R_V = 3.1 for all WRs as a reasonable approximation, and determined their absolute K-band magnitudes.

5.2.1. Detected companions

As described in Sect. 4.2, we detected wide companions in four WRs in our sample – WR 48, 89, 93 and 115. WR 89 and 115 have been classified and analyzed by Hamann et al. (2019) as single stars, owing to the lack of any companion contribution to their spectra. Using the K-band flux ratio we obtained from interferometry along with the absolute K-band magnitudes of the WRs, we calculated the absolute K-band magnitudes of their wide companions. Subsequently, assuming the companion is a main sequence OB star, we used synthetic photometry from Martins & Plez (2006), Pecaut & Mamajek (2013) to estimate the spectral types of the companions. Following is a short summary:

• WR 89: The absolute K magnitude for this WR was found to be –6.607 ± 0.179, while that for its companion was −5.057 ± 0.179. Based on the latter, we estimate the companion to be approximately an O5III giant.

• WR 115: The absolute K magnitude for this WR was found to be –4.595 ± 1.736, while that for the wide companion was –1.342 ± 1.736. The large error-bars in this case are primarily due to the uncertain Gaia parallax of WR 115 (0.58 ± 0.46). Consequently, it is not possible to constrain the spectral type of the companion. Nevertheless, based on the reported parallax, we estimated the companion to be a B2V dwarf.

• WR 48, WR 93: Both of these WRs have a composite spectrum and were excluded from the analysis of Sander et al. (2019). As a result, the above mentioned method could not be applied to these stars to determine the nature of their companions.

It is important to note that our estimates are based on the assumption that the companion is on the main sequence. Since we are basing the spectral type solely on one data point (absolute K magnitude), we also explored the possibility of a late-type pre-main sequence companion that is equally bright. Based on the pre-main sequence tracks computed by Siess et al. (2000), F/G/K-type stars can have absolute K magnitudes similar to those calculated above for our WR companions, provided they are very young (age < 1 Myr). Such stars are unlikely to reside in binaries with WRs, which are likely a few Myr old. We therefore exclude the possibility of pre-main sequence companions, although future observations in a different wavelength band could potentially help confirm the nature of the companions.

5.2.2. Detection limits

For WRs with no resolved companions, we computed the upper limits for companion detection in GRAVITY data. The limits were determined in terms of magnitude contrast with respect to the unresolved primary. We converted this contrast to a limiting flux ratio in the K-band. With the knowledge of absolute K magnitudes of the WRs and limiting flux ratios, we derived the corresponding absolute K magnitudes for the companion. Essentially, any companion brighter than this limiting magnitude would have yielded at least a 3σ detection. Therefore, we use this limiting absolute magnitude to infer upper limits on the mass of the companion.

Taking a conservative approach, we assumed the companion to be a main sequence dwarf. We again employed synthetic photometry provided by Pecaut & Mamajek (2013) to estimate the spectral type of the limiting companion. Figure 9 illustrates the absolute K magnitudes of the WRs, limiting magnitude contrasts and the absolute K magnitudes of the companions along with their spectral types. There is a scatter in the limiting magnitudes owing to the scatter in WR magnitudes as well as limiting contrasts, the latter being dependent on data quality. Nevertheless, the sample as a whole has a median detection limit around spectral type B4-5V, or about 5 M_⊙. The WNh stars have, on average, higher limits than cWRs due to their higher luminosity, and a median detection limit around spectral type B2V or mass of about 7 M_⊙.

Fig. 9. Absolute K-band magnitudes for WRs determined using stellar parameters (top points), limiting magnitude contrast from GRAVITY data (arrows), and the absolute K-band magnitudes for limiting companions along with corresponding spectral types (bottom). WNh stars (blue) have, on average, higher luminosities and consequently earlier limiting companion spectral types. WN (orange) and WC (red) stars show lower luminosities and later limiting companion spectral types.

6. Evolutionary implications

The multiplicity fraction of WRs, along with the period distribution of binaries is an informative indicator of their evolutionary history. In this work, we extended the discovery parameter space far beyond spectroscopic studies and were able to observationally constrain the multiplicity of WRs up to orbital periods of ∼10⁵ days and down to companion masses of ∼5 M_⊙. Many WRs in literature have been classified as candidate binaries based on the presence of absorption features or the dilution of emission lines. We only classified a WR as a binary if it either has a known orbital period or it is detected in our interferometric survey. Based on this definition, Table 1 summarizes the observed multiplicity fraction of WRs across different spectral classes. Our results therefore represent lower limits on WR multiplicity in the Galaxy based on conservative binary detection criteria. Here, we discuss the physical interpretation of our results, and corresponding evolutionary implications.

6.1. The curious case of WNh stars

The multiplicity fraction of the 13 WNh stars in our sample is particularly low (0.23 ± 0.12) compared to cWR stars (0.46 ± 0.10). As mentioned in the introduction, a significant fraction of WNh stars are thought to be very massive main sequence stars displaying the WR phenomenon due to their high luminosities (log(L/L_⊙)≳6) and subsequent strong winds. In such a scenario, the formation of WNh stars does not require binary interaction, allowing for a lower multiplicity fraction compared to cWRs. For example, WR 12 and WR 22 are known spectroscopic binaries with high luminosities (Hamann et al. 2019). The newly detected wide binary WR 89 shares similar properties, consisting of a high luminosity WNh star and an O star companion unlikely to have interacted before.

However, for the WNh stars that are truly single, the need of a high luminosity is even more pronounced. For the stars classified as single in our sample that were also reported in Hamann et al. (2019), we compared the WNh stars to their hydrogen-poor counterparts (WNs). While the WNh stars are on average more luminous, there is a significant overlap between the two populations. In particular, not all WNh stars classified as single are significantly more luminous than WN stars. Langer et al. (1994) discussed the possibility of WNh stars, H-rich or H-poor, being pre- or post-LBV evolutionary phase of very massive stars respectively. There are also observed connections with late-type WNh stars being seen after LBV variability stopped, e.g. Romano’s star in M33 (Smith et al. 2020) or HD5980 in the Small Magellanic Cloud (Koenigsberger et al. 2010). Alternative formation channels could also be contributing to the WNh population, e.g., Li et al. (2024) discuss the formation of WNh stars via binary mergers, which is indeed a form of binary interaction but results in single WNh stars. While the signatures of massive binary mergers have been reported before (Shenar et al. 2023; Frost et al. 2024), analogs for WNh stars are not known, potentially due to their violent winds.

WNh stars are most likely progenitors of cWR through the wind stripping channel. As such, it is surprising that their multiplicity is lower than that of the WNs. The higher the luminosity of the WR star, the higher is the companion detection limit. As discussed in Sect. 5.2.2, companions later than spectral type B2 V are unlikely to be detected for WNh stars in our survey. Such companions can also avoid spectroscopic detection due to their extreme mass ratios and lower luminosity. A combination of these factors might possibly explain the lower multiplicity fraction of WNh stars but would require quantitative modeling for confirmation.

6.2. Comparing WN and WC multiplicity

The spectroscopic survey of Dsilva et al. (2020, 2022, 2023) focused on a sample of 39 Galactic WRs to investigate multiplicity properties and evolutionary link between WN and WC stars. Since, spectroscopy is most efficient at short periods, the short-period WR binaries were most easily detected. They found a lack of short-period (≲100 days) WC binaries compared to their WN counterparts. In our sample however, the number of short-period WNs and WCs is comparable (see the literature review in Appendix D). Consequently, the inconsistency between WNs and WCs at short periods might be of statistical origin.

As discussed in the introduction, Dsilva et al. (2023) also used Bayesian inference to model the excess RV noise and to predict binary properties of WRs well beyond the spectroscopic regime. They inferred the multiplicity fraction of WCs to be much higher than WNs, especially at longer periods (≳100 days), which formed the rationale of our interferometric detection campaign. However, as seen in Table 1, we find no such overabundance of wide binaries for WCs compared to WNs. It is therefore reasonable to conclude that the evolutionary link of WNs evolving into WCs does not face any tension based on our results.

Furthermore, the lowest luminosity WR stars spend little to no time in the WC phase (Crowther 2007; Aguilera-Dena et al. 2022), suggesting a slightly higher occurrence of WN binaries. The values in Table 1 might reflect this expectation, although the overlapping error-bars do not establish it conclusively.

Finally, 14 stars were excluded from our analysis due to poor observations/data reduction, 9 of which are WCs (including 6 WCd stars) while only 3 are WNs. The exclusion of these stars from our sample potentially introduces a bias in the multiplicity statistics, particularly due to the link between dust-production and binarity (Lau et al. 2020a). While the different spectral types are represented in our sample, relative sample sizes of different subcategories of WRs might not be representative of the total WR population.

6.3. Comparing cWR binaries with O star binaries

With GRAVITY, we explored the angular separation range from ∼1 to 200 mas, corresponding projected separations from ∼1 to 300 au. The binaries in this range were detected based on just one epoch, and do not have a known orbit yet. Sana et al. (2014) presented a similar high angular resolution multiplicity survey of O-type stars using the VLTI/PIONIER and NACO/SAM instruments. They detected binaries with angular separations of ∼1 to 250 mas and magnitude contrasts ΔH ≲ 4 − 5. Tramper et al. (in prep.) further analyzed the O star sample to derive the projected separations (au) and mass ratios of the detected binaries. They also report the spectral types and corresponding masses for the O stars as well as their companions.

Massive O-type binaries are the likely progenitors of cWR binaries and a comparison between the two populations can be insightful to investigate their evolutionary connection. Only the massive O stars are thought to evolve into cWRs, with the single-star channel suggesting a mass threshold of ∼25 M_⊙ in the Galaxy, although it is highly uncertain with estimates ranging from 18 to 33 M_⊙ (Eldridge et al. 2008; Eldridge & Stanway 2016; Georgy et al. 2015; Limongi & Chieffi 2018; Dray & Tout 2003). On the other hand, the binary channel allows for lower mass progenitors, down to ∼20 M_⊙ (Shenar et al. 2020). To keep a conservative approach, we set a common threshold at 25 M_⊙ for any O star to be considered a potential cWR progenitor regardless of the channel.

To fix the parameter space for comparison of cWR and O star binaries, we converted the angular separation (mas) to physical projected separation (au) for WR binaries using their Gaia DR3 parallax measurements. Since our study probes the ∼1 to 300 au projected separation range, we only selected the O star binaries detected in this range. Furthermore, all companions to O stars detected in this separation range have derived spectral types earlier than B3V (Tramper et al., in prep.), quite similar to the detection limits for the cWRs in our sample as shown in Fig. 9. The possibility of pre-main sequence companions can be ruled out on similar grounds as in Sect. 5.2.1. We note here that WNh stars are excluded from the comparison since they do not resemble H-depleted WN or WC stars in an evolutionary context.

Figure 10 shows the cumulative companion fraction (total number of companions/sample size) plotted against the projected separation for our cWR sample and the massive O star sample. Also included is an additional line for a higher mass cut of 40 M_⊙ for the O star sample. While the bottom horizontal-axis shows the projected separation – a direct observable from interferometry; the top horizontal-axis provides an approximate period scale, assuming a total system mass of 40 M_⊙. The figure also shows an interaction boundary, which was calculated as the Roche radius of the primary for an extreme case of a 40 M_⊙ primary, mass ratio of 1 and eccentricity of 0.3. We equated the maximum radial extent of a 40 M_⊙ star (Marchant & Bodensteiner 2024) to the Roche radius and found the corresponding binary separation. This was also done for eccentricities of 0 and 0.5 to show an approximate spread in the extent of this boundary. Finally, a few cWR binaries not in our sample are indicated with arrows in the plot at approximated separations for illustration. These include WR 98a (Hendrix et al. 2016), WR 104 (Soulain et al. 2023), WR 112 (Lau et al. 2020b), WR 137, WR 138 (Richardson et al. 2016) and WR 140 (Lau et al. 2022), all of which have known orbital periods.

Fig. 10. Cumulative companion fraction and corresponding binomial error regions plotted against the projected binary separation for cWR stars in our sample (blue) and O star primaries more massive than 25 M⊙ (red) and 40 M⊙ (orange) from the SMASH+ survey (Tramper et al., in prep.). The bottom horizontal-axis range is bound by the limits of the GRAVITY detection regime. The top horizontal-axis shows the orbital period for corresponding binary separations assuming a total system mass of 40 M⊙. Furthermore, the black dotted line represents an interaction boundary for the relatively extreme case of M1 = 40 M⊙, q = 1 and e = 0.3. The gray shaded region spans the eccentricity range from e = 0 to e = 0.5. Also shown are a few WR binaries from literature approximately placed on the plot for illustration.

As the bottom horizontal-axis of the plot is based on the GRAVITY regime, all unresolved binaries (known from spectroscopic studies) lie below the projected separation of 1 au. Up to this separation, the companion fraction for cWRs is ∼38% while that of O stars is ∼50%. As such, this does not necessarily cause any evolutionary tension between the two populations, considering that a significant fraction of close O star binaries might end up merging – thus “destroying” the binary. However, as we uncover the new parameter space of separations ranging from 1 to 300 au, the multiplicity fractions of the two populations seem to diverge even more. Up to 300 au, the companion fraction is ∼90% for O stars and only ∼50% for WRs. This divergence can be discussed in two parts separated by the binary interaction limit defined by the maximum extent of an evolved massive star, thus about 25–50 au (≈5 − 10 × 10³ R_⊙).

Up to the interaction boundary, the multiplicity fraction of O stars steadily increases while there are virtually no detected cWR binaries in this range. In other words, most O star binaries in this regime seem to not have WR binary successors. Such dearth of WRs in wide binaries has also been reported in the Small (Schootemeijer et al. 2024) and Large Magellanic Cloud (LMC) (Bartzakos et al. 2001; Shenar et al. 2019), although only from spectroscopy. For our Galactic sample, a possible explanation for this can be that massive O stars might tend to expand considerably in their later evolutionary stages, initiating binary interaction. In case of stable mass transfer, the growing cumulative O star multiplicity fraction with period should qualitatively be preserved to the cWR stage since the orbital period is modified gradually. However, a drastic change such as that seen in Fig. 10 could reflect unstable mass transfer, with the interaction resulting in either a short period WR+OB binary or in a merger that may or may not evolve into the WR stage, depending on the mass of the merger product and the efficiency of the wind-stripping mechanism. In either case, the result of such interactions could be a lack of long-period cWR binaries. As mentioned earlier and indicated in Fig. 10, some cWR binaries are found in literature with orbital separations in the range of 2–20 au. However, some of them have non-zero eccentricities, raising doubts about their post-interaction status.

Beyond the interaction boundary, we again see a steady increase in the multiplicity fraction of O stars (for mass thresholds of 25 as well as 40 M_⊙) while there is virtually no increase for cWRs. If such O stars indeed do not interact and form cWRs by themselves via the single-star channel, their companions should be preserved, likely still on the main sequence. Additionally, their orbital separations should only change by a factor of two even in the extreme case of the system losing half of its total mass via primary wind mass-loss. The lack of cWR binaries in this regime is therefore concerning for the single-star channel. Based on our results, the mass threshold to produce cWR stars via the single-star channel might be even higher than 40 M_⊙ (see e.g. Fig. 10 in Langer 2012). Alternatively, these wide, less gravitationally bound companions might be lost due to dynamical interactions with their environment. Detailed simulations, accounting for the local birth density and dynamical history of each system would be needed to investigate the likelihood of such a scenario and whether it can explain the lack of very wide cWR binaries compared to their O star counterparts.

6.4. Implications for black hole binaries

The growing catalog of merging black hole (BH) binaries (Abbott et al. 2019, 2021, 2023) has garnered significant interest in their possible formation channels. Massive binary evolution is considered a strong potential contributor to such BH binaries. Langer et al. (2020) investigated the properties of OB star-black hole (OB+BH) systems at LMC metallicity using detailed binary evolution models as these represent a crucial intermediate stage in the formation of BH+BH binaries.

The evolutionary stage immediately preceding an OB+BH system is thought to be an OB+cWR binary, where the cWR undergoes a supernova/collapse to form the BH (see e.g. Fig. 1 in Langer et al. 2020). For example, Langer et al. (2020) assume the BH mass to be equal to the He-core mass of the cWR, and no momentum kick during its formation. As a result, the OB+BH systems are virtually OB+cWR systems at the time of collapse. Under these assumptions, our cWR binary sample offers a glimpse into the initial orbital properties of the OB+BH phase and allows for a direct comparison between observations and theoretical computations.

According to Langer et al. (2020), the period distribution of OB+BH binary models at LMC metallicity is expected to peak around 200 days (Fig. 6 in Langer et al. 2020), and is dominated by systems that have undergone Case B mass transfer. Pauli et al. (2022) conducted a similar population study for WRs in the LMC, finding peaks around 10 and 30 days in the period distributions of WN and WC binaries respectively. However, they also found a significant population of WR binary models at periods up to few 100 days (Figs. 6 and 7 in Pauli et al. 2022). While these simulations have been computed with LMC metallicity, Janssens et al. (2022) showed that one expects only a small difference of these predicted period distributions between Milky Way and LMC (see their Fig. B.1).

In literature, the lack of observed wide WR binaries predicted by such simulations has been so far attributed to the limitations of spectroscopic campaigns at longer periods. With GRAVITY being sensitive to periods in the range of months to decades, the present observational campaign provides a systematical exploration of the period range of OB+BH progenitors predicted from simulations. However, our survey failed to reveal any OB+cWR binaries in the proximity of the 200-day period peak predicted by evolutionary models. Similarly, no OB+cWR binaries were found in a projected separation range of ∼1 to 20 au, corresponding to a period range of few tens of days to a few thousand days. Of course, neutron star or black hole companions to the cWR stars cannot be excluded based on our survey. While a few cWR binaries have been found in this range by other studies (see Fig. 10), there is no sign of a period peak around 200 days.

The OB+BH binaries predicted in this period range were all products of Case B mass transfer. The lack of OB+WR binaries that we report might thus warrant a revision of this mass transfer phase. In particular, the absence of OB+cWR binaries suggests that most of the long-period interacting O star binaries undergo unstable mass transfer, potentially leading to mergers (Justham et al. 2014; Menon et al. 2024).

6.5. Single cWRs and others caveats

The high percentage (∼54%) of single cWRs is another significant result of our study. Even after expanding the discovery parameter space by two orders of magnitude in separation and flux contrast thanks to optical long-baseline interferometry, we could not find any new binaries aside from WR 89 (which was strongly suspected to be a binary from X-ray observations). It thus seems that there is a significant lack of WR binaries in general, and of long-period ones in particular (see again Fig. 10).

How to produce so many single cWR stars is not immediately evident. As there are few single O stars, many of the single cWRs could have a binary evolution history. However, binary mergers are not necessarily helpful in all cases, since they increase the H-rich envelope mass in the initial primary stars. Initial secondary stars could evolve into cWRs as single stars after their companions have been ejected from the orbit following a supernova. Such systems, however, will likely evolve through a previous O+cWR phase. Tailored population synthesis calculations are required to investigate whether a consistent solution is possible.

One possibility is that our sample is affected by a significant selection bias. For example, the VLTI-ATs required optically bright sources for the STRAP injections (V <  14) which introduces a bias against significantly reddened objects and indeed some known WC dust emitters (WR 98a, WR 112) were not included in this sample because of this limitation. While it would be good to complement the current survey with a dedicated VLTI+UTs survey that can observe reddened targets, it remains to be seen whether this possible bias alone is sufficient to explain the lack of long-period cWR binaries. Indeed, the numbers of dusty WR that we missed seems insufficient to reconcile observations and theory. In addition, their immediate WN progenitors, should not suffer from such biases and do not reveal further binaries in that period range.

7. Conclusion

We conducted an interferometric survey of 39 Galactic WRs with VLTI/GRAVITY to investigate their multiplicity properties. We detected wide companions for only four WRs in our sample, namely WR 48, WR 89, WR 93, and WR 115. In addition, we found four WRs with spatially resolved line-emitting regions and five WRs with resolved companions in differential phase data. We also detected a fully resolved component that contributes significantly to the K-band flux for most WRs in our sample.

Combining with spectroscopic studies, we determined the observed multiplicity fraction for our sample to be 0.38 ± 0.08; with 0.23 ± 0.12 for WNh, 0.55 ± 0.15 for WNs, and 0.40 ± 0.13 for WCs. Here, the classification of binaries was based on a known orbital period and/or a wide companion detected in our study. Companion detection limits were also calculated for the GRAVITY data and were converted from flux contrasts to absolute magnitudes assuming main sequence dwarf companions. Subsequently, we found the median detection limit for cWRs to be around 5 M_⊙, and that for WNh stars to be around 7 M_⊙.

The evolutionary implications of our results span the different spectral classes of WRs, and their connection to the predecessors and descendants of WR binaries (OB binaries and OB+BH binaries, respectively). The low multiplicity fraction of WNh stars could be potentially attributed to their higher luminosities, enabling them to form via the single-star channel, or making it more difficult to detect companions as compared to cWRs. WNh stars could also be products of binary mergers that inevitably now appear as single stars.

Within the Northern cWR population, Dsilva et al. (2023) reported a potential tension between the period distributions of WNs and WCs based on a modeling of the RV excess noise. If confirmed, cWR stars cannot be explained by the conventional WN → WC evolutionary link between them. However, we find the WN and WC multiplicity properties to be consistent within the error bars, and suggest that the tension reported in Dsilva et al. (2023) is either due to limited sample sizes or an additional RV variability signal not well modeled in Dsilva et al. (2023).

We also compared our cWR population with the massive O star population from the SMASH+ survey of Sana et al. (2014), and find a stark contrast in their multiplicity properties at long periods (10² − 10⁵ days). There is a significant lack of cWR wide binaries as compared to massive O stars, indicating that binary interactions might be “destroying” wide O star binaries instead of evolving them into wide cWR binaries. On the other hand, at very long periods beyond the regime of binary interaction, many of the wide O star binaries are expected to evolve into wide cWR binaries based on single star physics. The lack of such cWR binaries might raise doubts as to the effectiveness of the Conti scenario, although no robust conclusions can be drawn at present.

With regard to their eventual fate, many OB+cWR binaries are expected to evolve into OB+BH binaries via supernova or collapse of the cWR star. Binary evolution models for OB+BH stars in the LMC predict a peak in their period distribution at around 200 days. Assuming minimal mass-loss and no momentum kick during BH formation, this period distribution can then also represent the parent period distribution, that is, the distribution for the OB+cWR stars. With GRAVITY, we find no OB+cWR binaries in or around this period range, or even up to periods of a few years. We believe this implies that long-period interacting O star binaries might not always evolve into cWR+OB stars, and instead might merge through unstable mass transfer, possibly impeding the common-envelope channel toward double-BH merger at solar metallicity. Finally, we also discuss the high observed fraction of single cWRs stemming from observational biases or their true single nature. Future work on higher sensitivity observations and refined theoretical population models can lead the way to understand these objects better.

Data availability

All reduced GRAVITY data may be made available upon request to the authors. Supplementary material for Appendix D can be found on Zenodo (https://doi.org/10.5281/zenodo.14042165).

Acknowledgments

The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 and Horizon Europe research and innovation programme (grant agreement numbers 772225: MULTIPLES; and 101054731: Stellar-BHs-SDSS-V). A.A.C.S., and R.R.L. are supported by the German Deutsche Forschungsgemeinschaft, DFG in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander). L.R.P. acknowledges support by grants PID2019-105552RB-C41 and PID2022-137779OB-C41 funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe”. F.T. gratefully acknowledges support by grant PID2022-137779OB-C41, funded by the Spanish Ministry of Science, Innovation and Universities/State Agency of Research MICIU/AEI/10.13039/501100011033.

References

Appendix A: Summary of observations

Appendix B: Qualitative analysis of GRAVITY data

Three interferometric observables of interest captured by GRAVITY are |V|, T3PHI and DPHI, each capable of storing information regarding the size and structure of the target to varying extents. In addition, GRAVITY also provides a K-band spectrum, which is normalized and referred to as NFLUX. For a qualitative view, just inspecting the general features in these observables can motivate our modeling approach.

The visibility amplitude (|V|) is indicative of the extent to which the target is resolved. At its extremes, we have |V| = 1 at all baselines for an unresolved source, whereas |V| = 0 at all baselines for a fully resolved source. Partially resolved sources typically show intermediate behavior between the two extremes. In case the target is comprised of two components, the |V| data is sensitive to their flux ratio. In a simplified case of two components, the combined complex visibility follows the equation

$$\begin{array}{c}\hfill V=\frac{F_1{V}_1+F_2{V}_2}{F_1+F_2}.\end{array}$$(B.1)

Figure B.1 shows four simple examples of two-component systems to get a qualitative understanding of how |V| manifests for such targets if they are composed of a combination of unresolved or partially resolved (i.e., V∼1) and a fully resolved (V=0) component. Because the latter has zero visibility, we can ignore the phase term and the linearity of the complex visibility applies to the amplitude directly.

Fig. B.1. Four simple examples of |V| plots for two-component targets where one component is unresolved or slightly resolved and the other component is fully resolved. The first column describes the two components along with their flux properties, while the second column shows how the components add up together resulting in the combined |V|. The boxes in each plot show the spectrum. In interferometric data, we only have access to the final (green) plots and this table describes a qualitative approach to guessing the components present in the target based on |V|. Although the four examples serve as a good representative set, actual data and corresponding models can be much more complex, warranting the addition of more parameters.

It is important to note here that the examples shown are for simple cases, while it is possible for |V| to be much more complex owing to different spatial and/or spectral components in the system. When it comes to searching for appropriate models, a more comprehensive inspection is necessary. For example, in case of a binary system consisting of two distinct point sources, |V| depends on their flux ratio, separation as well as orientation. Such a system can only be modeled by keeping the unknown parameters as free variables and searching a grid in the parameter space to find the best-fit solution.

The closure phase (T3PHI) data is sensitive to departure from point asymmetry within the target being observed. A point-symmetric source typically displays zero closure phase. A non-zero closure phase will be observed for an object with a resolved non-point like symmetry, such as a non-equal binary. T3PHI is a strong indicator of a target harboring a binary system, and is fit along with |V| data using the grid search described above.

The differential phase (DPHI) data can help detect asymmetry at a certain wavelength relative to a different wavelength at scales much smaller than those accessible with |V|-T3PHI. A non-zero differential phase (DPHI) in a line indicates a wavelength-dependent shift in the photocenter. For example, in a very close binary system in which one star has a certain spectral line which the other does not, DPHI data might detect a signal at the spectral line wavelength, even though |V|-T3PHI are consistent with a single unresolved point source. Although DPHI is a powerful tool to detect such asymmetry, detailed modeling requires one to know the properties of the system beforehand, especially when it is unresolved in |V|-T3PHI data. Consequently, DPHI is mostly used for qualitative evaluation of targets, and only invoked in modeling when a system is resolved in |V|-T3PHI as well.

When modeling interferometric data, it is ideal to have a consistent fit across all observables - |V|, T3PHI, DPHI and NFLUX. In case there is no spectral feature in NFLUX and hence in |V|, and DPHI = 0, it is likely that |V| and T3PHI alone can be sufficient for modeling - the first two examples in Fig. B.1 fall in this category. If we do see a spectral feature in NFLUX as well as |V|, it is important to include NFLUX and DPHI for modeling the system - the third example in Fig. B.1 falls in this category. Lastly, if we see a spectral feature in NFLUX but not in |V|, and DPHI = 0, it is likely that the spectral feature is present in multiple components and is somehow suppressed in |V| due to the flux ratio following Eq. B.1. Such cases can be approached in two ways - 1. assuming flat spectra for all components without taking NFLUX into account, and simply modeling |V|-T3PHI to obtain general properties of the system; 2. performing detailed modeling including spectral lines and finding a solution to replicate NFLUX as well as |V| consistently, along with T3PHI and DPHI. The first approach is often sufficient for a companion search since |V|-T3PHI data alone can clearly indicate the presence (or lack) of a resolved binary or multiple system. It does not require introducing unknown/uncertain parameters to force a better fit to the data. The second approach is more informative, but also more challenging due to a higher number of parameters and potential degeneracies introduced with them. Consequently, the second approach is favored only for specific cases where we have physically motivated parameters to add to the model.

Appendix C: Grid search for resolved companions

The prominent signatures of a binary in interferometric data include a non-zero closure phase indicating asymmetry, and an approximately sinusoidally modulating visibility amplitude. The latter can be expressed analytically as

$$\begin{array}{c}\hfill |V|=\sqrt{\frac{1+f^2+2fcos\left(\frac{\overrightarrow{\rho }.\overrightarrow{x_0}}{\lambda }\right)}{1+f^2}},\end{array}$$(C.1)

where f is the flux ratio of the secondary to the primary, $\overrightarrow{\rho }$ is the spatial frequency (ρ = B/λ) and $\overrightarrow{x_0}$ is the vector pointing from the primary to the secondary (Millour 2014). While it is not possible to find an accurate binary model for interferometric data using a purely analytical approach, it can still approximately inform us on the orientation and separation of the binary. Consequently, a qualitative look at GRAVITY data can motivate the approach we take to find the best binary model.

Only four WRs in our sample (WR 89, WR 115, WR 48 and WR 93) displayed potential signatures of a resolved binary. Nevertheless, we performed a grid search for companions on all 39 WRs in our GRAVITY sample. The flux ratio of the secondary to the primary was left as a free parameter. We started with a coarse initial grid spanning the range ΔE = –100 to 100 mas and ΔN = –100 to 100 mas, with steps of 20 mas each. According to Eq. C.1, the wider a binary (higher the x₀), the more sinusoidally modulated its |V| data. In case the |V| data for a WR were mostly flat and the coarse grid search did not result in a good model, we performed a second, finer grid search much closer to the primary. This second search spanned the range ΔE = –15 to 15 mas and ΔN = –15 to 15 mas, with steps of 3 mas each. As an extension to the discussion about ${\chi }_{\mathrm{red}}^2$ in Sect. 3, we note here that including NFLUX and DPHI data for modeling while not taking every single spectral feature into account can lead to large values of ${\chi }_{\mathrm{red}}^2$. While we ensured that our solutions were consistent even when using |V|-T3PHI data alone, we report the results including NFLUX and DPHI for a more complete picture of the systems.

Based on the two searches described above, 35 WRs in our sample were found to have no wide companion detections that were significant. Here, we discuss the remaining four WRs.

C.1. WR 89

For WR 89, we found the best-fit model in the fine grid search. The companion was detected at ΔE = 0.02 ± 0.08 mas and ΔN = −5.18 ± 0.10 mas relative to the central primary, as shown in Fig. C.1.

Fig. C.1. The fine grid search for a companion for WR 89. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = 0.02 ± 0.08 mas and ΔN = −5.18 ± 0.10 mas with the primary star set at the origin. The color scale denotes ${\chi }_{\mathrm{red}}^2$ across the grid.

C.2. WR 115

The GRAVITY data for WR 115 displayed a high amount of sinusoidal modulation in |V| (see SM Fig. 51). The coarse as well as fine grid searches failed to obtain an accurate binary model for WR 115. We turned to Eq. C.1 for a qualitative estimate of a suitable model. The sinusoidal modulation in |V| is sensitive to the binary separation (x₀) as well as its orientation with respect to the baseline ($\overrightarrow{\rho }.\overrightarrow{x_0}$), with the highest modulation seen for baselines aligned with the binary. From SM Fig. 51, it is evident that the highest modulation is present in the J2A0 and G1A0 baselines.

Based on these inputs, we performed another coarse grid search in the ranges ΔE = –200 to –100 mas and ΔN = 100 to 200 mas; and ΔE = 100 to 200 mas and ΔN = –200 to –100 mas, with steps of 20 mas each. We found an approximate detection around ΔE = –120 mas and ΔN = 150 mas. Consequently, we performed a fine grid search closer to the approximate solution, and detected the companion at ΔE = −124.80 ± 0.04 mas and ΔN = 154.93 ± 0.02 (see Fig. C.2).

Fig. C.2. The fine grid search for a companion for WR 115. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = −124.80 ± 0.04 mas and ΔN = 154.93 ± 0.02 mas with the primary star set at the origin. The color scale denotes ${\chi }_{\mathrm{red}}^2$ across the grid.

C.3. WR 48

For WR 48, as mentioned in the main text of the paper, a known outer O supergiant companion is reported in Hartkopf et al. (1999), which is brighter than the central object. Although our GRAVITY data for WR 48 is incomplete, we see clear indications of this companion (see SM Figs. 14 and 42). In this case, we fixed the O supergiant at the origin and performed a grid search for the WR. We found the best-fit model in a fine grid search, where the WR was detected at ΔE = −1.59 ± 0.03 mas and ΔN = 12.56 ± 0.06 mas relative to the central O supergiant (see Fig. C.3). The main text of the paper reports the position of the O supergiant companion relative to the WR.

Fig. C.3. The fine grid search for a companion for WR 48. Here, the wide companion is brighter than the WR, so we fixed the companion at the origin and searched the WR position across the grid. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = 1.59 ± 0.03 mas and ΔN = −12.56 ± 0.06 mas with the wide companion set at the origin. The color scale denotes the ${\chi }_{\mathrm{red}}^2$ across the grid.

C.4. WR 93

The GRAVITY data for WR 93 showed strong indications of a wide companion (see SM Figs. 28 and 48), and were analyzed similarly to WR 89. We found the best-fit model in the fine grid search. The companion was detected at ΔE = 5.60 ± 0.03 mas and ΔN = 13.19 ± 0.03 mas relative to the central primary, as shown in Fig. C.4.

Fig. C.4. The fine grid search for a companion for WR 93. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = 5.60 ± 0.03 mas and ΔN = 13.19 ± 0.03 mas with the primary star set at the origin. The color scale denotes ${\chi }_{\mathrm{red}}^2$ across the grid.

More details about the companions such as flux ratios, physical separations and potential spectral types are discussed in the main text of the paper.

Appendix D: Comments on individual objects

Here we provide some details for each star in our sample, including the known properties and new results from interferometry. Following the descriptions for all WRs is Table 2 summarizing the best-fit geometric models for all stars. All corresponding figures can be found as part of the supplementary material uploaded on Zenodo (https://doi.org/10.5281/zenodo.14042165), and are referred to in this section as SM Fig. #.

WR 8 is classified as a single WN6/WC4 star by Crowther et al. (1995). It displays corotating-interaction-region-type (CIR-type) variability in its spectral lines (Chené & St-Louis 2011). GRAVITY data for WR 8 did not show signs of binarity in any of the interferometric observables. The data were fit reasonably well with a simple point source and a fully resolved component, as shown in SM Fig. 1.

WR 9 is classified as a WC5 + O7 spectroscopic binary with a 14.3-day orbital period (van der Hucht et al. 1988). The GRAVITY data were best modeled using a point source and a fully resolved component (see SM Fig. 2). Although the T3PHI residuals seemed to have a signal, we did not a find a binary or even a triple system model that could explain it. However, features corresponding to spectral lines could be seen in the |V| and DPHI data, likely indicating a close binary as discussed in Sect. 4.3. SM Fig. 37 shows the complete spectro-interferometric data for WR 9, along with the best-fit model consisting of two close point sources and a fully resolved component. The position of the companion is essentially unconstrained.

WR 12 is classified as a WN8h + OB spectroscopic binary with a 23.9-day orbital period, also known to harbor colliding winds (Fahed & Moffat 2012). The |V| and T3PHI data were modeled sufficiently well with a point source and a fully resolved component as shown in SM Fig. 3.

WR 14 is a WC7 initially classified as a potential SB1 spectroscopic binary (Shylaja 1990), which was later classified as single (Drissen et al. 1992). The GRAVITY data were best modeled using a central point source and fully resolved component, the latter not being significant (SM Fig. 4).

WR 15 is a WC6 classified as a single star (Sander et al. 2019). The GRAVITY data were best modeled using a point source and a fully resolved component, although leading to significant residuals (SM Fig. 5). A companion search for a binary companion or even an additional tertiary companion was not able to fit the data. We therefore speculate the origin of the residuals to be either some structure within the fully resolved component or poor visibility calibration.

WR 16 is a WN8h classified as a single star (Hamann et al. 2019), also faintly detected in X-rays (Skinner et al. 2012). The GRAVITY data for WR 16 largely resembled a point source with a slight contribution from a fully resolved flux (see SM Fig. 6). Although the data were noisy, on close inspection, we found features in |V| data corresponding to the spectral lines in NFLUX data.

As discussed in Sect. 4.4, the case of WR 16 can be similarly explained by a two component model for the central star - an unresolved component for the continuum and 2.06 μm spectral line; and a partially resolved component for the 2.11 and 2.16 μm spectral lines. The partially resolved component was modeled as a Gaussian with an fwhm of 0.98 ± 0.11 mas. SM Fig. 38 shows the best-fit model for the complete spectro-interferometric GRAVITY data for WR 16.

WR 18 is a WN4 classified as a single star (Hamann et al. 2019), also known to emit X-rays (Skinner et al. 2010). The interferometric data were best modeled using a point source primary and a fully resolved diffuse component, as shown in Fig. 3 in the main text.

WR 21 is a WN5 classified as a spectroscopic binary with an O4-6 companion and a 8.25-day orbital period (Fahed & Moffat 2012), showing the colliding wind phenomenon (Nazé et al. 2023). GRAVITY data (particularly |V|) for WR 21 show the presence of a significant fully resolved component. The |V|-T3PHI were modeled adequately with a point source and a fully resolved component (see SM Fig. 7). The data show slight deviation from point source behavior, although not significant enough to consider binarity. Higher SNR data might resolve this ambiguity.

WR 22 is a WN7h classified as an eclipsing spectroscopic binary with an O9III-V companion and an orbital period of about 80 days (Gosset et al. 1991). It has been studied extensively as an X-ray source as well (Gosset et al. 2009; Parkin & Gosset 2011). As discussed in Sect. 4.1 and shown in Fig. 2, WR 22 was consistent with a simple point source model, with a fully resolved component essentially absent.

WR 23 is a WC6 classified as a single star (Sander et al. 2019). The GRAVITY data were best modeled by a point source and a fully resolved emission component (see SM Fig. 8).

WR 24 is a WN6ha-w classified as a single star (Hamann et al. 2019), also showing X-ray emission (Skinner et al. 2010). The GRAVITY data were consistent with a point source as shown in SM Fig. 9, with negligible contribution from a fully resolved component.

WR 31 is a WN4 classified as a spectroscopic binary with an O8V companion and a 4.8-day orbital period (Niemela et al. 1985), also showing the colliding wind phenomenon (Fahed & Moffat 2012; Nazé et al. 2023). The GRAVITY data for WR 31 were best modeled by a point source and a fully resolved flux contributing significantly (see SM Fig. 10).

WR 31a is a WN11h/candidate LBV classified as a single star (Smith et al. 1994). The GRAVITY data for WR 31a were largely consistent with a point source and fully resolved component model as seen in SM Fig. 11. On close inspection however, similar to WR 16, it shows spectral features in |V|, indicating resolved line-emitting regions (see Sec. 4.4). In particular, the emission lines around 2.06 and 2.16 μm showed significant features in |V| data.

We modeled WR 31a using two components: an unresolved component for the continuum; and a partially resolved component for the 2.06 and 2.16 μm spectral lines. The partially resolved component was modeled as a Gaussian with an fwhm of 0.91 ± 0.07 mas. SM Fig. 39 shows the best-fit model for WR 31a.

WR 42 is a WC7 classified as a spectroscopic binary with an O7V companion and a 7.89-day orbital period (Davis et al. 1981). The GRAVITY data for WR 42 resembled a point source with a slight contribution from a fully resolved component (see SM Fig. 12). The |V| and DPHI data, however, did show some spectral features corresponding to the emission lines in the NFLUX data (see SM Fig. 40).

Similar to WR 9, this behavior can be explained by the presence of a very close companion that is unresolved in a |V|-T3PHI binary search, but slightly resolved in DPHI. This is consistent with the already confirmed binary nature of WR 42, but cannot be reliably used to determine the position of the companion.

WR 47 is a WN6 classified as a spectroscopic binary with an O5V companion and a 6.2-day orbital period, also displaying the colliding wind phenomenon (Fahed & Moffat 2012). The GRAVITY data for WR 47 were modeled sufficiently well with a point source and a fully resolved component (see SM Fig. 13).

With a more detailed approach, we found that the DPHI data show a spectral feature. The data are considerably noisy, and we tried an approximate model similar to WR 42, which could explain the DPHI feature. SM Fig. 41 shows the spectro-interferometric data for WR 47, along with an approximate model over-plotted on the data.

WR 48 is discussed in detail in the main text of the paper. SM Fig. 14 and SM Fig. 42 show the |V|-T3PHI and complete GRAVITY data respectively for WR 48, along with the best-fit model. Despite the data being noisy (and, for some baselines, absent), we found a good fit with a model consisting of two point sources and a negligible fully resolved component. A more detailed analysis of WR 48 will be covered in Gosset et al. (in prep.).

WR 52 is a WC4 classified as a single star (Sander et al. 2019). The GRAVITY data were consistent with a point source as shown in SM Fig. 15, with some contribution from a fully resolved component. The |V| data do show some features, although no simple model was found to perfectly describe them. Nevertheless, a wide binary was effectively ruled out after a thorough grid search.

WR 55 is a WN7 classified as a single star (Hamann et al. 2019). The interferometric data were best modeled using a point source primary and a fully resolved diffuse component, as shown in SM Fig. 16.

WR57 is a WC8 classified as a single star (Sander et al. 2019). GRAVITY data for WR 57, |V| data in particular, reveal peculiar features. SM Fig. 17 shows the |V|-T3PHI data modeled using a point source and a fully resolved component, the latter contributing to a much larger extent than other WRs. However, the closure and differential phases were both flat at zero. SM Fig. 43 shows the complete interferometric data with the model including spectral lines of the WR.

WR 113 (later in this appendix) was the only other WR in our sample that approximately resembles the |V| data for WR 57. The former is known to produce dust (Cohen et al. 1975), while the latter is not. More observations to improve the SNR, or even perform image reconstruction, were deemed necessary to understand the true origin of the significant fully resolved component in WR 57.

WR 66 is a WN8h classified as a single star (Hamann et al. 2019). It is known to show a 3.5-4 hr periodic photometric variability (Antokhin et al. 1995; Rauw et al. 1996; Nazé et al. 2021a) and faint X-ray emission (Skinner et al. 2021). The source of variability is not well understood, with the presence of a compact companion being one possibility among few others. The GRAVITY data for WR 66 were best modeled by a point source and a fully resolved component (see SM Fig. 18).

WR 75 is a WN6 classified as a single star (Hamann et al. 2019). The GRAVITY data were consistent with a point source as shown in SM Fig. 19, with some contribution from a fully resolved component.

WR 78 is a WN7h classified as a single star (Sander et al. 2019), also known to emit X-rays (Skinner et al. 2012). The GRAVITY data for WR 78 show features of spatially resolved line-emitting region for certain spectral lines, as described in detail in the main text of the paper. In context of companion search using |V|-T3PHI data, WR 78 was best modeled as a point source with no contribution from a fully resolved component (see SM Fig. 20).

WR79 is a WC7 classified as a spectroscopic binary with an O5-8 companion and a 8.9-day orbital period, also showing the colliding wind phenomenon (Luehrs 1997). The |V|-T3PHI data were largely consistent with a point source and a slightly contributing fully resolved component, as shown in SM Fig. 21. On closer inspection of the complete spectro-interferometric data, we found small spectral features in DPHI as well as |V| corresponding to the WR spectral lines in NFLUX. This can be explained by a close binary as described in Sect. 4.3. SM Fig. 44 shows the best-fit model for WR 79 consisting of a very close binary and a slight fully resolved component. Similar to other WRs of this kind, the position of the companion is not constrained.

WR 79a is a WN9ha/O8:Iafpe classified as a single star (Crowther & Bohannan 1997), also known to emit X-rays (Skinner et al. 2010). The visibilities are more than unity due to poor visibility calibration. The best-fit model for GRAVITY data of WR 79a was an unresolved point source by itself (SM Fig. 22). A binary search with or without the bad |V| data did not fit the data well, and we classified it as a single star.

WR 79b is a WN9ha/O6:Iafpe classified as a single star (Sota et al. 2014). Although the GRAVITY data are noisy, they were best explained with an unresolved point source and a fully resolved component (SM Fig. 23). The visibility residuals seem to contain some signal, but similar to WR 15, it can be attributed to the resolved component or poor visibility calibration.

WR 81 is a WC9 classified as a single star (Sander et al. 2019). The GRAVITY data for WR 81, particularly |V|-T3PHI, were well explained by a point source and a fully resolved component (see SM Fig. 24). On a closer look, the |V| data have a prominent spectral feature corresponding to the 2.06 μm line of the WR. Such a feature was only seen in WR 81 and WR 92 (see below). We modeled it by including the same spectral features in the fully resolved component (see SM Fig. 45), something that explains the data but is currently non-trivial to reconcile with any physical interpretation. Higher S/N data in the future can help constrain the properties of the fully resolved component in greater detail.

WR 85 is a WN6h classified as a single star (Hamann et al. 2019). The interferometric data were best fit by a point source and a fully resolved component (see SM Fig. 25). The |V| data show some residual signal which can be attributed to some structure in the fully resolved component.

WR 87 is a WN7h+abs classified as a single star, where abs stands for an absorption component (Smith et al. 1996). It is also a known X-ray source, although the origin of X-rays is uncertain (Nazé et al. 2013). The GRAVITY data for WR 87 were best fit by a point source and a slight fully resolved component, as shown in SM Fig. 26.

WR 89 is discussed in detail in the main text of the paper. SM Fig. 46 shows the complete spectro-interferometric data for WR 89 along with the best-fit model.

WR 92 is a WC9 classified as a single star (Sander et al. 2019). The |V|-T3PHI data for WR 92 were best modeled using a point source and a fully resolved component (see SM Fig. 27). The |V| data, however, have a prominent spectral feature corresponding to the 2.06 μm line of the WR similar to WR 81, and were also modeled similarly (see SM Fig. 47).

WR 93 is discussed in detail in the main text of the paper. SM Fig. 28 and SM Fig. 48 show the |V|-T3PHI and complete GRAVITY data respectively for WR 93, along with the best-fit model.

WR 97 is a WN5b classified as a spectroscopic binary with an O7 companion and a 12.6-day orbital period (Smith et al. 1996), also known to exhibit the colliding wind phenomenon (Nazé et al. 2021b). The GRAVITY data were best modeled by a central point source and a fully resolved component, as shown in SM Fig. 29.

WR 98 is a WN8/WC7 classified as a spectroscopic binary with an O8-9 companion and a 47.8-day period (Gamen & Niemela 2003; Faherty et al. 2014). It is also known to emit in the radio band, potentially due to colliding winds (Montes et al. 2011). The |V|-T3PHI data for WR 98 could be explained well by a point source model, with only a slight contribution from a fully resolved component (see SM Fig. 30). Including DPHI and NFLUX does, however, reveal spectral features indicating the presence of a close binary, as described in Sect. 4.3. Figure 5 shows the complete spectro-interferometric data along with the best-fit model.

WR 108 is a WN9ha classified as a single star Smith et al. (1996). The GRAVITY data for WR 108 were best fit by a point source and a fully resolved component, as shown in SM Fig. 31.

WR 110 is a WN5-6b classified as a single star, potentially hosting corotating interaction regions in its wind (Chené et al. 2011). The GRAVITY data for WR 110 were consistent with a point source model and a negligible fully resolved component (see SM Fig. 32). The |V| data, however, show features of a spatially resolved line-emitting region for multiple spectral lines, similar to WR 78 described in the main text of the paper. In this case however, the spectral features could not be constrained as simply as other similar cases. SM Fig. 49 shows the complete spectro-interferometric data for WR 110, along with an approximate model similar to WR 78. More observations with higher SNR were deemed necessary to better constrain its properties.

WR 111 is a WC5 classified as a single star (Sander et al. 2019). The GRAVITY data were best fit by a point source and a slight fully resolved component, as shown in SM Fig. 33.

WR 113 is a WC8d classified as a spectroscopic binary with an O8-9IV companion and a 29.7-day orbital period (Gaposchkin 1949), also known to produce dust (Hill et al. 2018). Shara et al. (2022) resolved a wide companion candidate at a separation of 1.16 arcseconds, although whether or not it is gravitationally bound the central binary is not known.

Monnier et al. (2007) studied the near-infrared sizes of dusty WRs including WR 113, finding its K-band “size” to be 27 ± 3 mas. In our GRAVITY data, we see strong evidence of this dusty disk reflecting in the strong |V| offset at all baselines (SM Fig. 34). The DPHI data shows spectral features, although in this case they could potentially arise due to the presence of a dusty disk. Constraining the size of the dusty disk was not possible based on the current data. Nevertheless, our best-fit model WR 113 was a central point source and a slightly off-center gaussian disk with an fwhm of around 15 mas (SM Fig. 50). To constrain the latter at a higher precision, more GRAVITY data would be necessary.

WR 114 is a WC5 classified as a single star (Sander et al. 2019). The GRAVITY data for WR 114 were best modeled with a point source and a fully resolved component (see SM Fig. 35).

WR 115 is discussed in detail in the main text of the paper. SM Fig. 36 and SM Fig. 51 show the |V|-T3PHI and complete GRAVITY data respectively for WR 115, along with the best-fit model.

All Tables

All Figures

	Fig. 1. Radius in solar radii where the optical depth is equal to unity based on a PoWR atmosphere model for WR 78. The blue curve only takes the continuum opacity into account, while the light-red curve also accounts for the line opacities.
In the text

Fig. 2. |V|-T3PHI data for the WN7h + O9III-V spectroscopic binary system WR 22. The different colors in T3PHI represent four combinations of 3 ATs each, while those in |V| represent the six baselines of 2 ATs each. The top panels show the best fit model (gray) along with the data, while the bottom panels show the residuals. The flat-at-zero T3PHI and flat-at-unity |V| data are typical of a single unresolved point source.

In the text

	Fig. 3. |V|-T3PHI data for the WN4 star WR 18 (similar to Fig. 2). The flat-at-zero T3PHI and flat but vertically offset |V| data are best explained by a combination of an unresolved point source and a fully resolved component.
In the text

	Fig. 4. |V|-T3PHI data for the WN8h+abs star WR 89 (similar to Fig. 2). The T3PHI and |V| data deviate drastically from symmetric/point-source-like behavior, and are best modeled by a wide binary system with two distinct point sources.
In the text

Fig. 5. GRAVITY data for WR 98, a WN8/WC7 + O8-9 spectroscopic binary, zoomed in at 2.10–2.22 μm along with the best-fit model shown in red. The top left panel shows the u-v coverage for the observation. The remaining panels, from left to right, show NFLUX, T3PHI, DPHI and |V| respectively, color-coded according to baselines. The data are zoomed in to a wavelength range of around 2.10–2.22 μm, with three prominent emission lines present at 2.11, 2.16 and 2.19 μm. The DPHI data (along with some of |V| data) show slight but clear features corresponding to the spectral lines. The best-fit model for such a system is a very close binary consisting of the line-emitting WR and a flat-spectrum companion, along with a fully resolved component.

In the text

Fig. 6. GRAVITY data for WR 78, a single WN7h star, along with the best-fit model shown in red. The top left panel shows the u-v coverage for the observation. The remaining panels, from left to right, show NFLUX, T3PHI, DPHI and |V| respectively, color-coded according to baselines. The data are zoomed in to a wavelength range of around 2.09–2.21 μm, with three strong emission lines present around 2.11, 2.16 and 2.19 μm. The |V| data show clear features corresponding to the first two spectral lines, decreasing in prominence from longer to shorter baselines. The best-fit model for such a system is an unresolved point source emitting the 2.19 μm spectral line and a slightly resolved Gaussian (fwhm ≈ 0.51 ± 0.02 mas) emitting the 2.11 and 2.16 μm spectral lines.

In the text

	Fig. 7. Line formation strengths ξ (Hillier 1987) at 2.11 (top), 2.16 (middle), and 2.19 μm per radius for the PoWR atmosphere model used in Fig. 1. The ξ values are normalized for better line-to-line comparison. The radii where most line formation occurs Rmax are highlighted by the dashed black lines.
In the text

Fig. 8. Companion detection limit plots for WR 18. The first plot shows a randomly created dense grid about 200 mas in radius around the central primary star that was used to inject a second point source and calculate the flux ratio needed for a 3σ detection. The scale for flux ratio is given in magnitude with reference to the primary. The second plot shows a histogram of limiting magnitudes over all the grid points. Lastly, the third plot shows detection limit as a function of radial distance, along with a rolling average in red.

In the text

Fig. 9. Absolute K-band magnitudes for WRs determined using stellar parameters (top points), limiting magnitude contrast from GRAVITY data (arrows), and the absolute K-band magnitudes for limiting companions along with corresponding spectral types (bottom). WNh stars (blue) have, on average, higher luminosities and consequently earlier limiting companion spectral types. WN (orange) and WC (red) stars show lower luminosities and later limiting companion spectral types.

In the text

Fig. 10. Cumulative companion fraction and corresponding binomial error regions plotted against the projected binary separation for cWR stars in our sample (blue) and O star primaries more massive than 25 M⊙ (red) and 40 M⊙ (orange) from the SMASH+ survey (Tramper et al., in prep.). The bottom horizontal-axis range is bound by the limits of the GRAVITY detection regime. The top horizontal-axis shows the orbital period for corresponding binary separations assuming a total system mass of 40 M⊙. Furthermore, the black dotted line represents an interaction boundary for the relatively extreme case of M1 = 40 M⊙, q = 1 and e = 0.3. The gray shaded region spans the eccentricity range from e = 0 to e = 0.5. Also shown are a few WR binaries from literature approximately placed on the plot for illustration.

In the text

Fig. B.1. Four simple examples of |V| plots for two-component targets where one component is unresolved or slightly resolved and the other component is fully resolved. The first column describes the two components along with their flux properties, while the second column shows how the components add up together resulting in the combined |V|. The boxes in each plot show the spectrum. In interferometric data, we only have access to the final (green) plots and this table describes a qualitative approach to guessing the components present in the target based on |V|. Although the four examples serve as a good representative set, actual data and corresponding models can be much more complex, warranting the addition of more parameters.

In the text

	Fig. C.1. The fine grid search for a companion for WR 89. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = 0.02 ± 0.08 mas and ΔN = −5.18 ± 0.10 mas with the primary star set at the origin. The color scale denotes ${\chi }_{\mathrm{red}}^2$ across the grid.
In the text

	Fig. C.2. The fine grid search for a companion for WR 115. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = −124.80 ± 0.04 mas and ΔN = 154.93 ± 0.02 mas with the primary star set at the origin. The color scale denotes ${\chi }_{\mathrm{red}}^2$ across the grid.
In the text

Fig. C.3. The fine grid search for a companion for WR 48. Here, the wide companion is brighter than the WR, so we fixed the companion at the origin and searched the WR position across the grid. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = 1.59 ± 0.03 mas and ΔN = −12.56 ± 0.06 mas with the wide companion set at the origin. The color scale denotes the ${\chi }_{\mathrm{red}}^2$ across the grid.

In the text

	Fig. C.4. The fine grid search for a companion for WR 93. We found the global minimum for ${\chi }_{\mathrm{red}}^2$ at ΔE = 5.60 ± 0.03 mas and ΔN = 13.19 ± 0.03 mas with the primary star set at the origin. The color scale denotes ${\chi }_{\mathrm{red}}^2$ across the grid.
In the text