# Investigating 39 Galactic Wolf-Rayet stars with VLTI/GRAVITY

## Uncovering a long-period binary desert*

K. Deshmukh<sup>1,2,★★</sup>①, H. Sana<sup>1,2</sup>①, A. Mérand<sup>3</sup>, E. Bordier<sup>4</sup>①, N. Langer<sup>5,6</sup>, J. Bodensteiner<sup>3</sup>, K. Dsilva<sup>7</sup>①, A. J. Frost<sup>8</sup>, E. Gosset<sup>9</sup>, J.-B. Le Bouquin<sup>10</sup>, R. R. Lefever<sup>11</sup>①, L. Mahy<sup>12</sup>①, L. R. Patrick<sup>13</sup>①, M. Reggiani<sup>1</sup>, A. A. C. Sander<sup>11</sup>①, T. Shenar<sup>14</sup>①, F. Tramper<sup>13</sup>①, J. I. Villaseñor<sup>15</sup>①, and I. Waisberg<sup>16,17</sup>

*(Affiliations can be found after the references)*

Received 23 September 2024 / Accepted 4 November 2024

## ABSTRACT

**Context.** Wolf-Rayet stars (WRs) represent one of the final evolutionary stages of massive stars and are thought to be the immediate progenitors of stellar-mass black holes. Their multiplicity characteristics form an important anchor point in single and binary population models for predicting gravitational-wave progenitors. Recent spectroscopic campaigns have suggested incompatible multiplicity fractions and period distributions for N- and C-rich Galactic WRs (WNs and WCs) at both short and long orbital periods, in contradiction with evolutionary model predictions.

**Aims.** In this work, we employed long-baseline infrared interferometry to investigate the multiplicity of WRs at long periods and explored the nature of their companions. We present a magnitude-limited ($K < 9$; $V < 14$) survey of 39 Galactic WRs, including 11 WN, 15 WC, and 13 H-rich WN (WNh) stars.

**Methods.** We used the *K*-band instrument GRAVITY at the Very Large Telescope Interferometer (VLTI) in Chile. The sensitivity of GRAVITY at spatial scales of ~1 to 200 milliarcseconds and flux contrast of 1% allowed an exploration of periods in the range $10^2-10^5$ d and companions down to ~5 $M_\odot$. We carried out a companion search for all our targets, with the aim of either finding wide companions or calculating detection limits. We also explored the rich GRAVITY dataset beyond a multiplicity search to look for other interesting properties of the WR sample.

**Results.** We detected wide companions with VLTI/GRAVITY for only four stars in our sample: WR 48, WR 89, WR 93, and WR 115. Combining our results with spectroscopic studies, we arrived at observed multiplicity fractions of $f_{\text{obs}}^{\text{WN}} = 0.55 \pm 0.15$, $f_{\text{obs}}^{\text{WC}} = 0.40 \pm 0.13$, and $f_{\text{obs}}^{\text{WNh}} = 0.23 \pm 0.12$. The multiplicity fractions and period distributions of WNs and WCs are consistent in our sample. For single WRs, we placed upper limits on the mass of potential companions down to ~5 $M_\odot$ for WNs and WCs, and ~7 $M_\odot$ for WNh stars. In addition, we also found other features in the GRAVITY dataset, such as (i) a diffuse extended component contributing significantly to the *K*-band flux in over half the WR sample; (ii) five known spectroscopic binaries resolved in differential phase data, which constitutes an alternative detection method for close binaries; and (iii) spatially resolved winds in four stars: WR 16, WR 31a, WR 78, and WR 110.

**Conclusions.** Our survey reveals a lack of intermediate- (a few hundred days) and long- (a few years to decades) period WR systems. The 200d peak in the period distributions of WR+OB and BH+OB binaries predicted by Case B mass-transfer binary evolution models is not seen in our data. The rich companionship of their O-type progenitors in this separation range suggests that the WR progenitor stars expand and interact with their companions, most likely through unstable mass transfer, resulting in either a short-period system or a merger.

**Key words.** techniques: interferometric – binaries: general – stars: massive – stars: Wolf-Rayet

## 1. Introduction

Massive stars with initial masses of greater than $8 M_\odot$ 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_\odot$, 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

<footnote>* Based on observations collected at the European Southern Observatory under ESO program IDs 111.24JN (PI: H. Sana) and 109.23CN (PI: I. Waisberg).</footnote>

<footnote>** Corresponding author; kunalprashant.deshmukh@kuleuven.be</footnote>
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<sup>+0.04</sup><sub>-0.22</sub>) as compared to WNs (0.52<sup>+0.14</sup><sub>-0.12</sub>), 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 binary 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

<footnote n="1">¹ http://pacrowther staff shef ac uk /WRcat /</footnote>
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 \sim 500$) or high ($R \sim 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 $\Delta 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<fnref n="2" /> (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 molecule, 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

<figure><img src="image_1.png" /><figcaption>**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.</figcaption></figure>

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.

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 Bry 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

<footnote n="2">² https://github.com/amerand/PMOIRED</footnote>
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^2_{red}$) 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^2_{red}$ 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^2_{red}$ 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\sigma$ 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.

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.

<figure><img src="image_2.png" /></figure>

**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.

## 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 $\Delta K \sim 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 $\Delta$East or $\Delta$E and $\Delta$North or $\Delta$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 $\Delta K < 5$, finding a $\leq 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
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

colliding wind binary based on its properties, but could not con-
firm 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_{\text{red}}^2$ at $\Delta\text{E} = 0.02 \pm 0.08$ mas and $\Delta\text{N} = -5.18 \pm 0.10$ mas
relative to the central primary. The flux of the companion is
$0.24 \pm 0.02$ times that of the WR star. Based on a Gaia DR3
parallax of $0.27 \pm 0.02$ mas, the projected separation of the com-
panion is about $19 \pm 1$ au.

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_{\text{red}}^2$ at $\Delta\text{E} = -124.80 \pm 0.04$ mas and $\Delta\text{N} = 154.93 \pm 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 \pm 0.01$ times that of the WR star. The Gaia parallax for WR 115 is $0.58 \pm 0.46$ mas, and cannot be used to constrain the physical projected separation due to its large error-bar.

**WR 48.** Also known as $\theta$ Mus A, WR 48 was initially clas-
sified as a spectroscopic binary with a 19-day orbital period

<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

(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 compan-
ion at $\Delta\text{E} = 1.59 \pm 0.03$ mas and $\Delta\text{N} = -12.56 \pm 0.06$ mas rela-
tive to the central binary. It has a flux $1.27 \pm 0.10$ times that of the
central binary. Based on the *Gaia* parallax of $0.46 \pm 0.07$ mas,
the projected separation is $28 \pm 4$ au. One more epoch of speckle
interferometry data (Hartkopf et al. 2012) and two more epochs
of archival GRAVITY data are available for $\theta$ 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 $\Delta\text{E} = 5.60 \pm 0.03$ mas and $\Delta\text{N} = 13.19 \pm 0.03$ mas. The flux of the companion is $0.20 \pm 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 \pm 0.03$ mas for WR 93, the projected separation is about $29 \pm 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
<figure><img src="image_5.png" /></figure>

**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.

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 ~($\lambda$/B) × $e_{DPHI}$/360° radian is necessary for a detection in DPHI, where $\lambda$ 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} \sim 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 (Bry + 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.

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 (Bry + 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 µm.
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Line formation strengths $\xi$ (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 $\xi$ values are normalized for better line-to-line comparison. The radii where most line formation occurs $R_{\text{max}}$ are highlighted by the dashed black lines.</figcaption></figure>

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_{\odot}$.

Using a PoWR atmosphere model of WR 78, the line formation strengths per radius $\xi$ (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 $\xi$ 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.

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
<table><caption>**Table 1.** Combined multiplicity properties from spectroscopy (literature) and interferometry (this work) for our sample of 39 WRs.</caption><thead><tr><th>WR class</th><th>Total stars</th><th>Single Stars</th><th>Unresolved spectroscopic companions</th><th>Resolved spectroscopic companions</th><th>Resolved wide companions</th><th>Total WRs with companions</th><th>Observed multiplicity fraction</th></tr></thead><tbody><tr><td><strong>WNh</strong></td><td>13</td><td>10</td><td>2</td><td>0</td><td>1</td><td>3</td><td><strong>0.23 ± 0.12</strong></td></tr><tr><td><strong>WN</strong></td><td>11</td><td>5</td><td>3</td><td>2</td><td>1</td><td>6</td><td><strong>0.55 ± 0.15</strong></td></tr><tr><td>WNL</td><td>6</td><td>3</td><td>0</td><td>2</td><td>1</td><td>3</td><td>0.50 ± 0.20</td></tr><tr><td>WNE</td><td>5</td><td>2</td><td>3</td><td>0</td><td>0</td><td>3</td><td>0.60 ± 0.22</td></tr><tr><td><strong>WC</strong></td><td>15</td><td>9</td><td>1</td><td>4</td><td>2<sup>(*)</sup></td><td>6</td><td><strong>0.40 ± 0.13</strong></td></tr><tr><td><strong>cWR</strong></td><td>26</td><td>14</td><td>4</td><td>6</td><td>3</td><td>12</td><td><strong>0.46 ± 0.10</strong></td></tr><tr><td><strong>All WR</strong></td><td>39</td><td>24</td><td>6</td><td>6</td><td>4</td><td>15</td><td><strong>0.38 ± 0.08</strong></td></tr></tbody></table>

**Notes.** Rows are divided into WNh, WN, and WC, with WN further divided into WNL and WNE. The last column lists the observed multiplicity fraction along with errors calculated using binomial statistics (\*(*) as WR 48 is a triple candidate, it is counted twice).

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.

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):

$$E_{B-V} = 1.21 \times E_{b-v} \quad (1)$$

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

$$A_V = R_V \times E_{B-V}, \quad (2)$$

$$A_K = 0.12 \times A_V. \quad (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

$$M_K = m_K - A_K - 5 \times \log_{10} \left( \frac{d}{10} \right), \quad (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.
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

### 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_\odot$. 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_\odot$.

## 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_\odot$. 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
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

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 $\Delta H \lesssim 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
<figure><img src="image_10.png" /></figure>

**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_\odot$ (red) and 40 $M_\odot$ (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_\odot$. Furthermore, the black dotted line represents an interaction boundary for the relatively extreme case of $M_1 = 40 M_\odot$, $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.

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_\odot$ in the Galaxy, although it is highly uncertain with estimates ranging from 18 to 33 $M_\odot$ (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_\odot$ (Shenar et al. 2020). To keep a conservative approach, we set a common threshold at 25 $M_\odot$ 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_\odot$ 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_\odot$. 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_\odot$ primary, mass ratio of 1 and eccentricity of 0.3. We equated the maximum radial extent of a 40 $M_\odot$ 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.

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^3 R_\odot$).

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_\odot$) 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 compan-
ions should be preserved, likely still on the main sequence. Addi-
tionally, 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 chan-
nel. Based on our results, the mass threshold to produce cWR
stars via the single-star channel might be even higher than 40 $M_{\odot}$
(see e.g. Fig. 10 in Langer 2012). Alternatively, these wide, less
gravitationally bound companions might be lost due to dynam-
ical 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 distri-
butions 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 dif-
ference of these predicted period distributions between Milky
Way and LMC (see their Fig. B.1).

In literature, the lack of observed wide WR binaries pre-
dicted 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 system-
atical 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 com-
panions to the cWR stars cannot be excluded based on our sur-
vey. 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 bina-
ries 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 signif-
icant 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 gen-
eral, and of long-period ones in particular (see again Fig. 10).

How to produce so many single cWR stars is not imme-
diately 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. Ini-
tial 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 calcula-
tions 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 proper-
ties. We detected wide companions for only four WRs in our
sample, namely WR 48, WR 89, WR 93, and WR 115. In addi-
tion, 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_\odot$, and that for WNh stars to be around 7 $M_\odot$.

The evolutionary implications of our results span the differ-
ent spectral classes of WRs, and their connection to the predeces-
sors 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²⁻¹⁰⁵ 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 bina-
ries are expected to evolve into OB+BH binaries via super-
nova or collapse of the cWR star. Binary evolution models for
OB+BH stars in the LMC predict a peak in their period dis-
tribution at around 200 days. Assuming minimal mass-loss and
no momentum kick during BH formation, this period distribu-
tion 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 popula-
tion 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).

Acknowledgements. The research leading to these results has received fund-
ing 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 Forschungs-gemeinschaft*, 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-1377790B-C41 funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe”. F.T. gratefully acknowledges support by grant PID2022-1377790B-C41, funded by the Spanish Ministry of Science, Innovation and Universities/State Agency of Research MICIU/AEI/10.13039/501100011033.

References

Abbott, D. C., & Conti, P. S. 1987, ARA&A, 25, 113

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 031040

Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Phys. Rev. X, 11, 021053

Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Phys. Rev. X, 13, 041039

Absil, O., Le Bouquin, J. B., Berger, J. P., et al. 2011, A&A, 535, A68

Aguilera-Dena, D. R., Langer, N., Antoniadis, J., et al. 2022, A&A, 661, A60

Aguilera-Dena, D. R., Müller, B., Antoniadis, J., et al. 2023, A&A, 671, A134

Antokhin, I. I., Bertrand, J. F., Lamontagne, R., & Moffat, A. F. J. 1995, IAU Symp., 163, 62

Bartzakos, P., Moffat, A. F. J., & Niemela, V. S. 2001, MNRAS, 324, 18

Chené, A. N., & St-Louis, N. 2011, ApJ, 736, 140

Chené, A. N., Moffat, A. F. J., Cameron, C., et al. 2011, ApJ, 735, 34

Clark, J. S., Lohr, M. E., Patrick, L. R., et al. 2018, A&A, 618, A2

Clark, J. S., Patrick, L. R., Najarro, F., Evans, C. J., & Lohr, M. 2021, A&A, 649, A43

Cohen, M., Barlow, M. J., & Kuhi, L. V. 1975, A&A, 40, 291

Conti, P. 1976, in Proc. 20th Colloq. Int. Astrophys. (Univ. Liege), 193

Crowther, P. A. 2007, ARA&A, 45, 177

Crowther, P. A., & Bohannan, B. 1997, A&A, 317, 532

Crowther, P. A., Smith, L. J., & Willis, A. J. 1995, in Wolf-Rayet Stars: Binaries; Colliding Winds; Evolution, eds. K. A. van der Hucht, & P. M. Williams, 163, 152

Crowther, P. A., Hadfield, L. J., Clark, J. S., Negueruela, I., & Vacca, W. D. 2006, MNRAS, 372, 1407

Davis, A. B., Moffat, A. F. J., & Niemela, V. S. 1981, ApJ, 244, 528

Dray, L. M., & Tout, C. A. 2003, MNRAS, 341, 299

Drissen, L., Robert, C., & Moffat, A. F. J. 1992, ApJ, 386, 288

Dsilva, K., Shenar, T., Sana, H., & Marchant, P. 2020, A&A, 641, A26

Dsilva, K., Shenar, T., Sana, H., & Marchant, P. 2022, A&A, 664, A93

Dsilva, K., Shenar, T., Sana, H., & Marchant, P. 2023, A&A, 674, A88

Eldridge, J. J., & Stanway, E. R. 2016, MNRAS, 462, 3302

Eldridge, J. J., Izzard, R. G., & Tout, C. A. 2008, MNRAS, 384, 1109

Fahed, R., & Moffat, A. F. J. 2012, MNRAS, 424, 1601

Faherty, J. K., Shara, M. M., Zurek, D., Kanarek, G., & Moffat, A. F. J. 2014, AJ, 147, 115

Figer, D. F., McLean, I. S., & Najarro, F. 1997, ApJ, 486, 420

Fitzpatrick, E. L. 1999, PASP, 111, 63

Frost, A. J., Bodensteiner, J., Rivinius, T., et al. 2022, A&A, 659, L3

Frost, A. J., Sana, H., Mahy, L., et al. 2024, Science, 384, 214

Gallenne, A., Mérand, A., Kervella, P., et al. 2015, A&A, 579, A68

Gamen, R. C., & Niemela, V. S. 2003, IAU Symp., 212, 184

Gamow, G. 1943, ApJ, 98, 500

Gaposchkin, S. 1949, Peremennye Zvezdy, 7, 36

Georgy, C., Ekström, S., Hirschi, R., et al. 2015, in Wolf-Rayet Stars, eds. W. R. Hamann, A. Sander, & H. Todt, 229

Gosset, E., Remy, M., Manfroid, J., et al. 1991, Inf. Bull. Var. Stars, 3571, 1

Gosset, E., Nazé, Y., Sana, H., Rauw, G., & Vreux, J. M. 2009, A&A, 508, 805

Haiman, Z., & Loeb, A. 1997, ApJ, 483, 21

Hamann, W. R., Gräfener, G., Liermann, A., et al. 2019, A&A, 625, A57

Hartkopf, W. I., Mason, B. D., Gies, D. R., et al. 1999, AJ, 118, 509

Hartkopf, W. I., Tokovinin, A., & Mason, B. D. 2012, AJ, 143, 42

Hendrix, T., Keppens, R., van Marle, A. J., et al. 2016, MNRAS, 460, 3975

Higgins, E. R., Sander, A. A. C., Vink, J. S., & Hirschi, R. 2021, MNRAS, 505, 4874

Hill, G. M., Moffat, A. F. J., & St-Louis, N. 2002, MNRAS, 335, 1069

Hill, G. M., Moffat, A. F. J., & St-Louis, N. 2018, MNRAS, 474, 2987

Hillier, D. J. 1987, ApJS, 63, 965

Hillier, D. J. 2015, in Wolf-Rayet Stars, eds. W. R. Hamann, A. Sander, & H. Todt, 65

Hopkins, P. F., Kereš, D., Oňorbe, J., et al. 2014, MNRAS, 445, 581

Janssens, S., Shenar, T., Sana, H., et al. 2022, A&A, 658, A129

Justham, S., Podsialdowski, P., & Vink, J. S. 2014, ApJ, 796, 121

Koenigsberger, G., Georgiev, L., Hillier, D. J., et al. 2010, AJ, 139, 2600

Lachaume, R., Rabus, M., Jordán, A., et al. 2019, MNRAS, 484, 2656

Lamberts, A., Millour, F., Liermann, A., et al. 2017, MNRAS, 468, 2655
Lamers, H. J. G. L. M., Maeder, A., Schmutz, W., & Cassinelli, J. P. 1991, ApJ, 368, 538

Langer, N. 2012, ARA&A, 50, 107

Langer, N., Hamann, W. R., Lennon, M., et al. 1994, A&A, 290, 819

Langer, N., Schürmann, C., Stoll, K., et al. 2020, A&A, 638, A39

Lapeyrere, V., Kervella, P., Lacour, S., et al. 2014, SPIE Conf. Ser., 9146, 91462D

Lau, R. M., Eldridge, J. J., Hankins, M. J., et al. 2020a, ApJ, 898, 74

Lau, R. M., Hankins, M. J., Han, Y., et al. 2020b, ApJ, 900, 190

Lau, R. M., Hankins, M. J., Han, Y., et al. 2022, Nat. Astron., 6, 1308

Lefever, R. R., Sander, A. A. C., Shenar, T., et al. 2023, MNRAS, 521, 1374

Li, Z., Zhu, C., Lü, G., et al. 2024, ApJ, 969, 160

Lima, E. F., Bica, E., Bonatto, C., & Saito, R. K. 2014, A&A, 568, A16

Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13

Luehrs, S. 1997, PASP, 109, 504

Mahy, L., Lanthermann, C., Hutsemékers, D., et al. 2022, A&A, 657, A4

Marchant, P., & Bodensteiner, J. 2024, ARA&A, 62, 21

Martins, F., & Plez, B. 2006, A&A, 457, 637

Menon, A., Ercolino, A., Urbaneja, M. A., et al. 2024, ApJ, 963, L42

Mérand, A. 2022, SPIE Conf. Ser., 12183, 121831N

Millour, F. 2014, EAS Publ. Ser., 69, 17

Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15

Moffat, A. F. J., & Seggewiss, W. 1977, A&A, 54, 607

Monnier, J. D., Tuthill, P. G., Danchi, W. C., Murphy, N., & Harries, T. J. 2007, ApJ, 655, 1033

Monnier, J. D., Zhao, M., Pedretti, E., et al. 2011, ApJ, 742, L1

Montes, G., Pérez-Torres, M. A., Alberdi, A., & González, R. F. 2009, ApJ, 705, 899

Montes, G., González, R. F., Cantó, J., Pérez-Torres, M. A., & Alberdi, A. 2011, A&A, 531, A52

Nazé, Y., Rauw, G., Sana, H., & Corcoran, M. F. 2013, A&A, 555, A83

Nazé, Y., Rauw, G., & Gosset, E. 2021a, MNRAS, 502, 5038

Nazé, Y., Gosset, E., & Marechal, Q. 2021b, MNRAS, 501, 4214

Nazé, Y., Rauw, G., Johnson, R., Gosset, E., & Hoffman, J. L. 2023, MNRAS, 526, 2167

Neugent, K. F., & Massey, P. 2014, ApJ, 789, 10

Niemela, V. S., Mandrini, C. H., & Mendez, R. H. 1985, Rev. Mex. Astron. Astrofis., 11, 143

Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457

Paczyński, B. 1967, Acta Astron., 17, 355

Parkin, E. R., & Gosset, E. 2011, A&A, 530, A119

Pauli, D., Langer, N., Aguilera-Dena, D. R., Wang, C., & Marchant, P. 2022, A&A, 667, A58

Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9

Rajagopal, J. 2010, Rev. Mex. Astron. Astrofis. Conf. Ser., 38, 54

Rauw, G., Gosset, E., Manfroid, J., Vreux, J. M., & Claeskens, J. F. 1996, A&A, 306, 783

Richardson, N. D., Shenar, T., Roy-Loubier, O., et al. 2016, MNRAS, 461, 4115

Richardson, N. D., Lee, L., Schaefer, G., et al. 2021, ApJ, 908, L3

Richardson, N. D., Daly, A. R., Williams, P. M., et al. 2024, ApJ, 969, 140

Rosslowe, C. K., & Crowther, P. A. 2015, MNRAS, 447, 2322

Rosslowe, C. K., & Crowther, P. A. 2018, MNRAS, 473, 2853

Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444

Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15

Sander, A. A. C., Hamann, W. R., Todt, H., et al. 2019, A&A, 621, A92

Sander, A. A. C., Vink, J. S., & Hamann, W. R. 2020, MNRAS, 491, 4406

Sander, A. A. C., Lefever, R. R., Poniatowski, L. G., et al. 2023, A&A, 670, A83

Schootemeijer, A., & Langer, N. 2018, A&A, 611, A75

Schootemeijer, A., Shenar, T., Langer, N., et al. 2024, A&A, 689, A157

Shara, M. M., Howell, S. B., Furlan, E., et al. 2022, MNRAS, 509, 2897

Shenar, T., Sablowski, D. P., Hainich, R., et al. 2019, A&A, 627, A151

Shenar, T., Gilkis, A., Vink, J. S., Sana, H., & Sander, A. A. C. 2020, A&A, 634, A79

Shenar, T., Wade, G. A., Marchant, P., et al. 2023, Science, 381, 761

Shylaja, B. S. 1990, Ap&SS, 164, 63

Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593

Skinner, S. L., Zhekov, S. A., Güdel, M., Schmutz, W., & Sokal, K. R. 2010, AJ, 139, 825

Skinner, S. L., Zhekov, S. A., Güdel, M., Schmutz, W., & Sokal, K. R. 2012, AJ, 143, 116

Skinner, S. L., Schmutz, W., Güdel, M., & Zhekov, S. A. 2021, Res. Notes Am. Astron. Soc., 5, 125

Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77

Smith, L. F. 1968, MNRAS, 138, 109

Smith, N. 2014, ARA&A, 52, 487

Smith, L. J., Crowther, P. A., & Prinja, R. K. 1994, A&A, 281, 833

Smith, L. F., Shara, M. M., & Moffat, A. F. J. 1996, MNRAS, 281, 163

Smith, N., E Andrews, J., Moe, M., et al. 2020, MNRAS, 492, 5897

Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014, ApJS, 211, 10

Soulain, A., Lamberts, A., Millour, F., Tuthill, P., & Lau, R. M. 2023, MNRAS, 518, 3211

St-Louis, N., Chené, A. N., Schnurr, O., & Nicol, M. H. 2009, ApJ, 698, 1951

Turner, D. G. 1982, IAU Symp., 99, 57

van der Hucht, K. A. 2001, New Astron. Rev., 45, 135

van der Hucht, K. A., Hidayat, B., Admiranto, A. G., Supelli, K. R., & Doom, C. 1988, A&A, 199, 217

Vanbeveren, D., & Conti, P. S. 1980, A&A, 88, 230

Vanbeveren, D., De Donder, E., Van Bever, J., Van Rensbergen, W., & De Loore, C. 1998, New Astron., 3, 443

Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015

1 Institute of Astronomy, KU Leuven, Celestijnlaan 200D, 3001 Leuven, Belgium

2 Leuven Gravity Institute, KU Leuven, Celestijnenlaan 200D, box 2415, 3001 Leuven, Belgium

3 European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany

4 I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany

5 Argelander Institut für Astronomie, Auf dem Hügel 71, DE-53121 Bonn, Germany

6 Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, DE-53121 Bonn, Germany

7 Institut d'Astronomie et d'Astrophysique, Université Libre de Bruxelles CP 226, Boulevard du Triomphe, 1050 Brussels, Belgium

8 European Southern Observatory, Santiago, Chile

9 Groupe d'Astrophysique des Hautes Energies, STAR, Université de Liège, Quartier Agora (B5c, Institut d'Astrophysique et de Géophysique), Allée du 6 Août 19c, B-4000 Sart Tilman, Liège, Belgium

10 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

11 Zentrum für Astronomie der Universität Heidelberg, Astronomisches Rechen-Institut, Mönchhofstr. 12-14, 69120 Heidelberg, Germany

12 Royal Observatory of Belgium, B-1180 Brussels, Belgium

13 Departamento de Astrofísica, Centro de Astrobiología, (CSIC-INTA), Ctra. Torrejón a Ajalvir, km 4, 28850 Torrejón de Ardoz, Madrid, Spain

14 The School of Physics and Astronomy, Tel Aviv University, Tel Aviv 6997801, Israel

15 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

16 Independent Researcher, Rehovot, Israel

17 Department of Particle Physics and Astrophysics, Weizmann Institute of Science, Rehovot 76100, Israel
# Appendix A: Summary of observations

<table><caption>Table A.1. GRAVITY observations of WRs.</caption><thead><tr><th>WR #</th><th>Calibration Star ID</th><th>MJD of Observation</th><th>Spectral Resolution</th></tr></thead><tbody><tr><td colspan="4">P111.24JN</td></tr><tr><td>8</td><td>HD 64545</td><td>60058.062</td><td>Medium</td></tr><tr><td>9</td><td>HD 61103</td><td>60295.214</td><td>Medium</td></tr><tr><td>12</td><td>HD 75027</td><td>60088.029</td><td>Medium</td></tr><tr><td>14</td><td>HD 75027</td><td>60058.013</td><td>Medium</td></tr><tr><td>15</td><td>HD 80347</td><td>60264.312</td><td>Medium</td></tr><tr><td>16</td><td>HD 85779</td><td>60064.099</td><td>Medium</td></tr><tr><td>18</td><td>HD 88687</td><td>60066.125</td><td>Medium</td></tr><tr><td>21</td><td>HD 88687</td><td>60326.230</td><td>Medium</td></tr><tr><td>22</td><td>HD 90274</td><td>60089.074</td><td>High</td></tr><tr><td>23</td><td>HD 93808</td><td>60060.064</td><td>Medium</td></tr><tr><td>24</td><td>HD 91652</td><td>60351.197</td><td>High</td></tr><tr><td>31</td><td>HD 96010</td><td>60332.294</td><td>Medium</td></tr><tr><td>31a</td><td>HD 95096</td><td>60351.397</td><td>Medium</td></tr><tr><td>42</td><td>HD 98677</td><td>60089.029</td><td>Medium</td></tr><tr><td>47</td><td>HD 107837</td><td>60089.140</td><td>Medium</td></tr><tr><td>48</td><td>HD 116703</td><td>60065.182</td><td>High</td></tr><tr><td>52</td><td>HD 116140</td><td>60088.192</td><td>Medium</td></tr><tr><td>55</td><td>HD 116471</td><td>60332.327</td><td>Medium</td></tr><tr><td>57</td><td>HD 116007</td><td>60332.375</td><td>Medium</td></tr><tr><td>66</td><td>HD 135378</td><td>60146.056</td><td>Medium</td></tr><tr><td>75</td><td>HD 147155</td><td>60180.008</td><td>Medium</td></tr><tr><td>79</td><td>HD 151790</td><td>60109.202</td><td>High</td></tr><tr><td>79b</td><td>HD 152370</td><td>60109.227</td><td>Medium</td></tr><tr><td>81</td><td>HD 153951</td><td>60147.186</td><td>Medium</td></tr><tr><td>85</td><td>HD 156641</td><td>60176.170</td><td>Medium</td></tr><tr><td>87</td><td>HD 155808</td><td>60176.116</td><td>Medium</td></tr><tr><td>89</td><td>HD 155808</td><td>60175.149</td><td>Medium</td></tr><tr><td>92</td><td>HD 156291</td><td>60146.097</td><td>Medium</td></tr><tr><td>93</td><td>HD 158021</td><td>60116.186</td><td>High</td></tr><tr><td>97</td><td>HD 159488</td><td>60146.181</td><td>Medium</td></tr><tr><td>98</td><td>HD 159488</td><td>60179.001</td><td>Medium</td></tr><tr><td>108</td><td>HD 166031</td><td>60182.147</td><td>Medium</td></tr><tr><td>110</td><td>HD 166240</td><td>60148.158</td><td>Medium</td></tr><tr><td>111</td><td>HD 166031</td><td>60177.158</td><td>Medium</td></tr><tr><td>113</td><td>HD 169783</td><td>60181.999</td><td>High</td></tr><tr><td>114</td><td>HD 167811</td><td>60179.153</td><td>Medium</td></tr><tr><td>115</td><td>HD 169783</td><td>60180.176</td><td>Medium</td></tr><tr><td colspan="4">P109.23CN</td></tr><tr><td>78</td><td>HD 96568</td><td>59715.111</td><td>High</td></tr><tr><td>79a</td><td>HD 160263</td><td>59794.043</td><td>High</td></tr></tbody></table>

Notes: The columns list the WR number, corresponding calibration star, time of observation and GRAVITY spectral resolution used (Medium R~500 or High R~4000). Also mentioned are the ESO Programs 111.24JN and 109.23CN, under which all WRs were observed.

# 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

$$V = \frac{F_1 V_1 + F_2 V_2}{F_1 + F_2} \qquad (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 \sim 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.

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
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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

$$|V| = \sqrt{\frac{1 + f^2 + 2f \cos(\vec{\rho} \cdot \vec{x}_0)}{1 + f^2}}, \quad (C.1)$$

where $f$ is the flux ratio of the secondary to the primary, $\vec{\rho}$ is the spatial frequency ($\rho = B/\lambda$) and $\vec{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
<figure><img src="image_12.png" /><figcaption>**Fig. C.1.** The fine grid search for a companion for WR 89. We found the global minimum for $\chi_{\text{red}}^2$ at $\Delta E = 0.02 \pm 0.08$ mas and $\Delta N = -5.18 \pm 0.10$ mas with the primary star set at the origin. The color scale denotes $\chi_{\text{red}}^2$ across the grid.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. C.2.** The fine grid search for a companion for WR 115. We found the global minimum for $\chi_{\text{red}}^2$ at $\Delta E = -124.80 \pm 0.04$ mas and $\Delta N = 154.93 \pm 0.02$ mas with the primary star set at the origin. The color scale denotes $\chi_{\text{red}}^2$ across the grid.</figcaption></figure>

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 $\Delta E = -100$ to 100 mas and $\Delta N = -100$ to 100 mas, with steps of 20 mas each. According to Eq. C.1, the wider a binary (higher the $x_0$), 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 $\Delta E = -15$ to 15 mas and $\Delta N = -15$ to 15 mas, with steps of 3 mas each. As an extension to the discussion about $\chi_{\text{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_{\text{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 $\Delta E = 0.02 \pm 0.08$ mas and $\Delta N = -5.18 \pm 0.10$ mas relative to the central primary, as shown in Fig. C.1.

### 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_0$) as well as its orientation with respect to the baseline ($\vec{p} \cdot \vec{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 $\Delta E = -200$ to $-100$ mas and $\Delta N = 100$ to 200 mas; and $\Delta E = 100$ to 200 mas and $\Delta N = -200$ to $-100$ mas, with steps of 20 mas each. We found an approximate detection around $\Delta E = -120$ mas and $\Delta N = 150$ mas. Consequently, we performed a fine grid search closer to the approximate solution, and detected the companion at $\Delta E = -124.80 \pm 0.04$ mas and $\Delta N = 154.93 \pm 0.02$ (see Fig. C.2).

### 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 $\Delta E = -1.59 \pm 0.03$ mas and $\Delta N = 12.56 \pm 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.

### 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 $\Delta E = 5.60 \pm 0.03$ mas and $\Delta N = 13.19 \pm 0.03$ mas relative to the central primary, as shown in Fig. C.4.

More details about the companions such as flux ratios, physical separations and potential spectral types are discussed in the main text of the paper.
<figure><img src="image_14.png" /><figcaption>**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_{\text{red}}^2$ at $\Delta E = 1.59\pm0.03$ mas and $\Delta N = -12.56\pm0.06$ mas with the wide companion set at the origin. The color scale denotes the $\chi_{\text{red}}^2$ across the grid.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. C.4.** The fine grid search for a companion for WR 93. We found the global minimum for $\chi_{\text{red}}^2$ at $\Delta E = 5.60\pm0.03$ mas and $\Delta N = 13.19\pm0.03$ mas with the primary star set at the origin. The color scale denotes $\chi_{\text{red}}^2$ across the grid.</figcaption></figure>

## 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 D.1 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 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 con-
sistent 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.

**WR 57** 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).

**WR 79** 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:1afpe 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.
<table><caption>**Table D.1.** Parametric modeling of all 39 WRs in our GRAVITY sample using the PMOIRED tool.</caption><thead><tr><th>WR #</th><th>Spectral Type</th><th>Point Source<br>fwhm (mas)</th><th>Fully Resolved<br>Component (flux % w.r.t WR)</th><th>Resolved<br>Line-emitting Region (mas)</th><th>Wide Binary<br>Separation (mas)</th><th>Data<br>Quality</th></tr></thead><tbody><tr><td>24</td><td>WN6ha-w</td><td>0.12 ± 0.06</td><td>2 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>85</td><td>WN6h</td><td>0.46 ± 0.19</td><td>18 ± 2</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>22</td><td>WN7h + O9III-V</td><td>0.12 ± 0.02</td><td>0</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>78</td><td>WN7h</td><td>0.002 (fixed)</td><td>0</td><td>0.51 ± 0.02</td><td>-</td><td>Good</td></tr><tr><td>87</td><td>WN7h+abs</td><td>0.00</td><td>10 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>12</td><td>WN8h + OB</td><td>0.74 ± 0.05</td><td>18 ± 1</td><td>-</td><td>-</td><td>Poor</td></tr><tr><td>16</td><td>WN8h</td><td>0.002 (fixed)</td><td>26 ± 9</td><td>0.98 ± 0.11</td><td>-</td><td>Poor</td></tr><tr><td>66</td><td>WN8h</td><td>0.78 ± 0.08</td><td>21 ± 2</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>89</td><td>WN8h+abs</td><td>0.65 ± 0.32</td><td>15 ± 7</td><td>-</td><td>5.18 ± 0.10</td><td>Good</td></tr><tr><td>79a</td><td>WN9ha/O8:Ifape</td><td>0.00</td><td>0</td><td>-</td><td>-</td><td>Poor</td></tr><tr><td>79b</td><td>WN9ha/O6:Ifape</td><td>0.00</td><td>16 ± 1</td><td>-</td><td>-</td><td>Poor</td></tr><tr><td>108</td><td>WN9ha</td><td>0.49 ± 0.10</td><td>6 ± 1</td><td>-</td><td>-</td><td>Poor</td></tr><tr><td>31a</td><td>WN11h</td><td>0.002 (fixed)</td><td>23 ± 3</td><td>0.91 ± 0.07</td><td>-</td><td>Average</td></tr><tr><td>47</td><td>WN6 + O5V</td><td>0.53 ± 0.05</td><td>12 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>75</td><td>WN6</td><td>0.51 ± 0.05</td><td>13 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>115</td><td>WN6-w</td><td>1.02 ± 0.16</td><td>31 ± 4</td><td>-</td><td>198.94 ± 0.03</td><td>Good</td></tr><tr><td>8</td><td>WN6</td><td>0.81 ± 0.07</td><td>24 ± 2</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>55</td><td>WN7</td><td>0.63 ± 0.02</td><td>12 ± 1</td><td>-</td><td>-</td><td>Good</td></tr><tr><td>98</td><td>WN8 + O8-9</td><td>0.30 ± 0.06</td><td>7 ± 1</td><td>-</td><td>-</td><td>Good</td></tr><tr><td>18</td><td>WN4</td><td>0.70 ± 0.02</td><td>16 ± 1</td><td>-</td><td>-</td><td>Good</td></tr><tr><td>31</td><td>WN4 + O8V</td><td>0.64 ± 0.15</td><td>32 ± 4</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>21</td><td>WN5 + O4-6</td><td>0.70 ± 0.06</td><td>20 ± 2</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>97</td><td>WN5b + O7</td><td>0.75 ± 0.04</td><td>18 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>110</td><td>WN5-6b</td><td>0.002 (fixed)</td><td>0</td><td>Not Constrained</td><td>-</td><td>Good</td></tr><tr><td>52</td><td>WC4</td><td>0.31 ± 0.15</td><td>19 ± 2</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>9</td><td>WC5 + O7</td><td>0.51 ± 0.07</td><td>35 ± 8</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>111</td><td>WC5</td><td>0.12 ± 0.10</td><td>9 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>114</td><td>WC5</td><td>0.72 ± 0.06</td><td>14 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>15</td><td>WC6</td><td>0.70 ± 0.06</td><td>13 ± 1</td><td>-</td><td>-</td><td>Good</td></tr><tr><td>23</td><td>WC6</td><td>0.45 ± 0.01</td><td>5 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>48</td><td>WC6 + O6-7V</td><td>0.57 ± 0.31</td><td>11 ± 7</td><td>-</td><td>12.66 ± 0.06</td><td>Poor</td></tr><tr><td>14</td><td>WC7</td><td>0.38 ± 0.04</td><td>8 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>42</td><td>WC7 + O7V</td><td>0.30 ± 0.08</td><td>12 ± 1</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>79</td><td>WC7 + O5-8</td><td>0.00</td><td>0</td><td>-</td><td>-</td><td>Average</td></tr><tr><td>93</td><td>WC7 + O7-9</td><td>0.62 ± 0.09</td><td>3 ± 1</td><td>-</td><td>14.33 ± 0.03</td><td>Good</td></tr><tr><td>57</td><td>WC8</td><td>0.40 ± 0.38</td><td>123 ± 13</td><td>-</td><td>-</td><td>Good</td></tr><tr><td>113</td><td>WC8d + O8-9IV</td><td>1.32 ± 0.03</td><td>335 ± 7*</td><td>-</td><td>-</td><td>Good</td></tr><tr><td>81</td><td>WC9</td><td>0.62 ± 0.14</td><td>21 ± 2</td><td>-</td><td>-</td><td>Good</td></tr><tr><td>92</td><td>WC9</td><td>0.20 ± 0.17</td><td>36 ± 3</td><td>-</td><td>-</td><td>Average</td></tr></tbody></table>

**Notes:** Shown are the WR number (column 1), spectral type (column 2), fwhm of the point source used to model the WR (column 3) and contribution of the fully resolved component (column 4). For WRs with resolved line-emitting regions, column 5 provides the fwhm of the resolved component, while column 6 provides projected separations for resolved wide binaries. Lastly, column 7 gives a purely qualitative evaluation of GRAVITY data for each WR determined on case-by-case inspection. (* In case of WR 113, a Gaussian disk was used instead of the fully resolved component)