© The Authors 2023

1. Introduction

Wolf-Rayet (WR) stars are a spectroscopic class of objects with strong, broad emission lines of helium and Hβ. They are thought to be the descendants of main-sequence O stars with initial masses greater than about 20–30 M_⊙, as well as immediate progenitors of black holes. Wolf-Rayet binary systems are therefore natural progenitors of black-hole binaries, which could coalesce to produce observable gravitational waves. Based on the chemical composition of their spectra, WR stars are divided into nitrogen- (WN), carbon- (WC), and oxygen-rich (WO) WR stars. These elements are thought to be the exposed products of hydrogen, helium, and carbon burning, respectively. Conti (1976) hypothesised that main-sequence O stars could strip themselves through their stellar winds and appear as WR stars before ending their lives as compact objects (called the Conti scenario). However, most massive stars live in binary systems that will interact over the course of their lifetime (e.g., Sana et al. 2012), and in this way, envelope stripping via a companion can also form WR stars. The dominant formation channel, and hence subsequent evolution, of WR stars is still an open question (Paczyński 1967; Neugent & Massey 2014; Shenar et al. 2019, 2020). Constraining the intrinsic multiplicity properties of the WN and WC populations allows us to understand the consequences of binary interaction on the final stages of stellar evolution. Moreover, orbital analysis of WR binaries directly allows us to measure masses of black-hole progenitors prior to core-collapse in a model-independent fashion.

The observed binary fraction of Galactic WR stars has been reported to be ∼0.4 in the literature, with a total of 227 objects (van der Hucht 2001; henceforth VDH01). This well-studied sample consisted of 127 WN stars, 87 WC stars, 10 WN/WC stars, and 3 WO stars. The reported binary fraction of ∼0.4 includes photometric, spectroscopic, and visual binaries, alongside candidates flagged as binaries due to the dilution of the emission lines in their spectra. As VDH01 is a compilation of studies with different data sets collected with different instruments and analysed using various techniques, it is non-trivial to correct for the observational biases. As a result, determining the intrinsic multiplicity properties of the Galactic WR population using these observational statistics is problematic.

The Galactic Catalogue of WR stars¹ (henceforth GCWR; originally Appendix 1 in Rosslowe & Crowther 2015) has a total of 667 Galactic WR stars (v1.25). Most of the recent additions to this catalogue were discovered through near- and mid-infrared surveys (Mauerhan et al. 2009; Shara et al. 2009, 2012; Faherty et al. 2014; Rosslowe & Crowther 2018). Unfortunately, not all of these WR stars have been investigated for possible companions. In order to obtain a reasonable observational-bias correction using spectroscopy, for example, it is of paramount importance to accurately quantify the uncertainties on the radial velocity (RV) measurements. With this in mind, we undertook a magnitude-limited (V ≤ 12 mag) spectroscopic survey of northern Galactic WR stars with the Mercator telescope on La Palma in 2017. Using the High-Efficiency and high-Resolution Mercator Echelle Spectrograph (HERMES; Raskin et al. 2011), we obtained at least six epochs of data for 39 northern Galactic WR stars.

The 12 WC stars in our sample were analysed in Dsilva et al. (2020, henceforth Paper I). We determined the observed binary fraction ($f_{\mathrm{obs}}^{\mathrm{WC}}$) to be 0.58 ± 0.14. Using Monte Carlo (MC) simulations with the measured RV uncertainties, we derived the intrinsic binary fraction ($f_{\mathrm{int}}^{\mathrm{WC}}$) of the Galactic WC population to be at least 0.72 at the 10% significance limit with orbital periods up to 10⁵ d. The 16 early-type WN (WNE) stars were analysed in Dsilva et al. (2022, henceforth Paper II), where the observed binary fraction ($f_{\mathrm{obs}}^{\mathrm{WNE}}$) was found to be 0.44 ± 0.12. After adopting a Bayesian framework, our MC simulations revealed the intrinsic binary fraction of WNE stars ($f_{\mathrm{int}}^{\mathrm{WNE}}$) to be $0.{56}_{-0.15}^{+0.20}$ with orbital periods explored up to 10⁵ d. The orbital period distribution of the parent WNE population revealed that the majority of binaries have periods shorter than 10 d, similar to what is found for main-sequence O stars (Sana et al. 2012).

Using the same Bayesian framework, Paper II also re-derived the multiplicity properties for the WC population. We found that the underlying orbital period distribution was significantly different from that derived for the WNE population ($f_{\mathrm{int}}^{\mathrm{WC}}=0.{96}_{-0.22}^{+0.04}$ and π^WC $=1.{90}_{-1.25}^{+1.26}$), with a majority of WC binaries having orbital periods greater than 1000 d. By comparing the distributions for the orbital period of the WNE and WC populations, we found a lack of short-period WC binaries (P <  10 d), questioning the evolutionary link between WNE and WC binaries. This deficit is backed by the literature, where the number of short-period WC binaries is much smaller than that of the corresponding WNE binaries, even though they should be easiest to detect given the large Doppler motion exhibited by their spectral lines.

As the third paper in this series, the focus of this work is to analyse the 11 WNL stars in our sample and to investigate the multiplicity properties of the underlying late-type nitrogen-rich (WNL) WR population. WNL stars can either be classical (post-main-sequence, hydrogen-rich) or very massive main-sequence (hydrogen-burning) WR stars (e.g., de Koter et al. 1997). It is challenging to differentiate between the two, although the main-sequence WNL stars are generally more luminous (and more massive) than classical WNL stars (Hamann et al. 2019). While attempting to link the WNL population with the WNE and WC populations through binary evolution, it is important to differentiate between main-sequence and classical WNL stars. When the known binaries are excluded, the luminosities of the stars in our sample are lower than or around 10⁶ L_⊙ (Hamann et al. 2019). We aim to increase the analysed sample size and explore the connection between the multiplicity properties of main-sequence O stars, WN, and WC stars. The sample and data are presented in Sect. 2. The method of measuring RVs is explained in Sect. 3. Section 4 presents our results. The discussion is presented in Sect. 5, and our conclusions are presented in Sect. 6.

2. Sample and data reduction

From the GCWR, we selected WN stars that have spectral types later than WN5. Entries listed as ‘WN5-6’ were considered to be WN5 for simplicity. For objects to be visible with the Mercator telescope, we applied a magnitude cut of V ≤ 12 and δ ≥ −30°. For stars with missing V-band magnitudes, we used their narrow-band v magnitudes (Smith 1968; Massey 1984) and applied the criterion v  ≤  13. This resulted in a sample of 11 WNL stars. We discuss the selection biases in Sect. 5.2.

Over the course of the RV monitoring campaign, we obtained at least six epochs with the 1.2 m Mercator telescope on La Palma using the HERMES spectrograph. HERMES covers the optical regime with a wavelength range from 3800 Å to 9000 Å with a resolving power of R = λ/Δλ ∼ 85 000. We also used archival HERMES data in our analysis when available, resulting in a time baseline of two to eight years. The number of spectra and time coverage for the 11 WNL stars in our sample is shown in Table 1.

The data reduction and normalisation is described in Paper I. In short, the data were first reduced using the standard HERMES pipeline (Raskin et al. 2011) and then corrected for the instrumental response (Paper I, Royer et al., in prep.). In addition to the bias and flat-field correction, the HERMES pipeline also corrects for barycentric motion. Following this, the spectra were corrected for telluric contamination using Molecfit (Smette et al. 2015; Kausch et al. 2015). After correcting for these effects, the spectra are only affected by interstellar (and maybe circumstellar) reddening. In order to normalise them in a homogeneous fashion and to minimise the human systematics, we used a normalised WNL model spectrum from the Potsdam Wolf-Rayet code (Gräfener et al. 2002; Todt et al. 2015) to identify pseudo-continuum regions around 8100 Å in the red and 5100 Å in the blue. We then anchored a continuum WNL model to the red and applied a reddening to fit the spectrum in the blue. After exploring various reddening laws, we concluded that they did not affect the homogeneity of the normalisation process. We therefore decided to use the reddening law from Fitzpatrick (2004).

3. Radial velocity measurements

As the spectra of WR stars are dominated by strong and broad emission lines, measuring RVs is a challenging task. We decided to use cross-correlation to measure RVs for the stars in our sample. An implicit assumption of cross-correlation is that the template is an accurate representation of the data. Therefore, the choice of template is important. Because stellar atmospheric models usually do not accurately reproduce the complex line profiles of WR spectra and the line profiles cannot be reproduced by Gaussian or Lorentzian functions, we decided to construct and use a spectrum with a high signal-to-noise ratio (S/N) as a template. We briefly describe the method and assumptions below.

3.1. Cross-correlation for measuring radial velocities

The cross-correlation function (CCF) is a convolution of the template with the data, and it is described using a log-likelihood function (Zucker 2003, see Paper I for more details). The RV of an epoch is measured by maximising this function, and the uncertainty can be derived from maximum log-likelihood theory. The measured RVs are then used to create a co-added spectrum, which is in the rest frame of the initial template. The co-added spectrum is then used to re-measure the RVs. This iterative process is continued until the measured RVs do not change within their measurement errors. One caveat of this method is that the measured RVs are relative to an epoch, and not absolute. This is not a problem when probing RV variations to identify binary candidates (as in this work), but an absolute shift must be applied when combining them with RVs from the literature.

3.2. Wind variability and its effect on radial velocity measurements

Comparing the WNL stars in our sample to the WCs (Paper I) and WNEs (Paper II), we observe significantly more variability in the WNL line profiles. This is expected because strong line-profile variability has been observed in WNL stars before (St-Louis et al. 2009; Chené & St-Louis 2011; Michaux et al. 2014; Chené et al. 2020). This directly increases the measurement uncertainties on the RVs, as the assumption of the template reproducing the data partially breaks down. In addition to correcting for observational biases, it is also important to understand the effect of this variability on the RV measurements.

In order to reduce the effect of wind variability on RV measurements, we selected spectral lines that were least affected by spectral variability for each star. Since the line-forming regions in WR stars are spatially extended, spectral lines of different elements probe different physical regions in the outflowing wind. Therefore, the behaviour of these lines is sensitive to local wind physics and can exhibit disparate RVs. To address this, we used lines of He II, N III, N IV, and in some cases, N V, to measure RVs. Additionally, for each star, we used the combination of all the above-mentioned masks and called it ‘full spec’. The final masks for each object are indicated in Appendix B. In Paper I and Paper II, we obtained high-cadence data of the long-period binaries WR 137 and WR 138 away from their periastron passages. This allowed us to quantify the amplitude of this variability on the RV measurements in WC and WNE stars, which were found to be ∼5 and ∼15 km s⁻¹ respectively.

In this work, we focus on WR 136, which has a spectral type of WN6b(h). Although it is classified as an SB1? system in the GCWR with a period of 4.5 d (Koenigsberger et al. 1980; Aslanov & Cherepashchuk 1981), there is no derived orbit for this system. Further studies with photometry (Moffat & Shara 1986) and polarimetry (Robert & Moffat 1989) failed to confirm this period. We chose WR 136 for this purpose due to a lack of a long-period binary in our sample. Assuming it is a single star, we obtained high-cadence data with HERMES over two periods. We first procured 18 spectra over 16 days in August 2019. Second, we obtained between 1 and 6 epochs within a night on seven occasions between August and September of 2019. Finally, we obtained between 6 and 15 epochs per night over ten nights in October 2019. This resulted in a total of 87 high-cadence spectra.

We used lines of He II, N III, N IV, and N V and the Balmer lines to measure RVs for WR 136. The RV dispersion is shown in Fig. 1. The observed RV dispersion comprises a statistical part (i.e. the measurement errors, σ_p) and contribution from wind variability (σ_w). Assuming they are uncorrelated, the weighted standard deviation of the RVs for each mask can be represented by

Fig. 1. High-cadence RV measurements for WR 136 over a variety of masks. In brackets at the bottom are values of (σw, σRV) (km s−1). As there are 87 data points for each mask, the legend has been omitted for clarity.

$$\begin{array}{c}\hfill {\sigma }_{\text{RV}}=\sqrt{{\sigma }_p^2+{\sigma }_w^2}.\end{array}$$(1)

By measuring σ_RV and σ_p, we obtain an upper limit on the wind variability. For each mask, values of (σ_RV, σ_w) are indicated in Fig. 1. We find that lines of He II are less affected by wind variability than N V λ4945. One possible reason is that the N V λ4945 line is a very weak line in the spectrum of WR 136, limiting its RV accuracy. In any case, the trend of He II lines showing greater stability than lines of nitrogen is unusual (see e.g., Paper II). A plethora of reasons might explain this, such as different wind structures or simply different wind physics across the various line-forming regions. This good example demonstrates that each WR star is unique. The RV amplitudes for the masks are similar to what was found in Koenigsberger et al. (1980), although we did not find the same period as they reported after analysing the data with a Lomb-Scargle periodogram (Lomb 1976; Scargle 1982). We used diagnostic plots such as Fig. 1 for each object to select an appropriate mask for measuring RVs. These are indicated in Appendices A and B.

4. Results

4.1. Observed binary fraction and detection probability

After measuring the RVs, we computed the peak-to-peak amplitude (ΔRV) for each object. We chose a threshold C above which we attributed the observed ΔRV to Doppler motion in a binary system. Figure 2 is a non-parametric threshold plot that shows the observed binary fraction as a function of the chosen threshold. By moving right to left, we attribute the RV amplitude of these stars to Doppler motion, and classify them as candidate spectroscopic binaries.

Fig. 2. Non-parametric threshold plot showing the observed binary fraction as a function of the threshold C. The adopted threshold of 50 km s−1 is shown by the dot-dashed vertical line. The entries in bold are confirmed spectroscopic binaries in the literature.

In principle, the chosen threshold should be high enough to avoid false positives due to intrinsic variability. The high-cadence study on WR 136 demonstrated that with the appropriate mask, the measured peak-to-peak RV amplitude due to intrinsic variability is ∼15 km s⁻¹. This is similar to what was found for the N V lines in WR 138 (Paper I). Therefore, we adopted a similar strategy and chose a threshold that is at least three times higher than the peak-to-peak intrinsic variability, at C = 50 km s⁻¹. This resulted in the classification of 4 out of 11 stars as RV variable objects for which binarity is likely the cause. This corresponds to an observed binary fraction ($f_{\mathrm{obs}}^{\mathrm{WNL}}$) of 0.36 ± 0.15.

4.2. Multiplicity properties: WNL population

We performed a Bayesian analysis with MC simulations with a method developed in Paper II to formally account for the observational biases. The ideal solution is to calculate the likelihood of the observed RV time series for each object based on the sampling and possible binary properties of the system. However, this is too computationally intensive, and instead, we divided the objects into multiple ΔRV bins and tried to reproduce them with our simulations. The ΔRV bins were ΔRV < 50 km s⁻¹ (seven objects), 50  ≤   ΔRV  <   250 km s⁻¹ (two objects), 250  ≤   ΔRV  <   650 km s⁻¹ (two objects), ΔRV ≥ 650 km s⁻¹ (no object).

The parameters included in the model for the population with binaries are the period distribution, intrinsic binary fraction, inclination, eccentricity, mass ratio, time of periastron passage, and the orientation of the binary system in three-dimensional space. The period distribution is characterised by the minimum period ($\log P_{\min }^{\mathrm{WNL}}$), maximum period ($\log P_{\max }^{\mathrm{WNL}}$), and the power-law index (π^WNL). The power-law index π determines the shape of the distribution of binaries with orbital periods P as follows:

$$\begin{array}{c}\hfill p(logP)\sim {(logP)}^{\pi }.\end{array}$$(2)

Along with the period distribution, we included the intrinsic binary fraction, $f_{\mathrm{int}}^{\mathrm{WNL}}$, as a model parameter. The eccentricities were drawn from a flat distribution between 0.0 and 0.9, with random inclinations, times of periastron passage, and three-dimensional orientations of the orbital planes. We also assumed a flat mass-ratio distribution from 0.1–2, with the mass of the WNL star to be 20 M_⊙. We refrained from using the observed distributions of eccentricity and mass-ratio since they show observational biases, which we aim to correct for. Moreover, we explored various mass-ratio distributions in Paper II and noted that they did not affect the posteriors significantly. Given the presence of very short-period binaries in the WNL sample, we fixed the value of $\log P_{\min }^{\mathrm{WNL}}$ at 0.15 [d] (approximately 1 d). We explored values of $\log P_{\max }^{\mathrm{WNL}}$ from 1.5 to 5.0 [d] in steps of 0.1 dex, π^WNL from −2.0 to 1.0 in steps of 0.1, and $f_{\mathrm{int}}^{\mathrm{WNL}}$ from 0.30 to 1.00 in steps of 0.02. We also used prior information on the known orbital periods for those classified as binaries (Table 2), and required at least three objects with 1 <  P <  10 d.

With the framework set up, we simulated 10 000 sets of 11 WNL stars for each P − M₂ grid point, resulting in almost 4 × 10⁸ populations. We calculated the three-dimensional likelihood of reproducing the observed distribution of objects in ΔRV bins. Assuming flat priors on $\log P_{\max }^{\mathrm{WNL}}$, π^WNL and $f_{\mathrm{int}}^{\mathrm{WNL}}$, we computed their marginalised posterior likelihoods (Fig. 3). We computed the 68% credible intervals for each posterior (the highest-density interval, HDI68). For the three model parameters, our posteriors yield values of $f_{\mathrm{int}}^{\mathrm{WNL}}$ $=0.{42}_{-0.12}^{+0.15}$, π^WNL $=-0.{70}_{-1.02}^{+0.73}$ and $\log P_{\max }^{\mathrm{WNL}}=4.{90}_{-3.40}^{+0.09}$ [d] for the parent WNL population. The value of $\log P_{\max }^{\mathrm{WNL}}$ is set by the upper limit of the grid due to the lack of long-period constraints introduced by the orbital bins because no WNL binaries are confirmed at P >  10 d.

Fig. 3. Marginalised posteriors of $\log P_{\max }^{\mathrm{WNL}}$, $f_{\mathrm{int}}^{\mathrm{WNL}}$, and πWNL for the Galactic WNL population. The one-dimensional posteriors are calculated assuming flat priors. For each posterior, the solid red lines show HDI68, and the dashed red line shows the mode.

The derived multiplicity properties agree with those of the Galactic WNE population (Paper II). An interesting difference is that of the upper-period cutoff, which is better constrained for the WNE population ($\log P_{\max }^{\mathrm{WNE}}=4.{60}_{-0.77}^{0.40}$) given the prior information on the presence of long-period binaries. We lack these constraints for the WNL population, where the only orbital periods that are constrained in the literature are shorter than 10 days (Table 2). This results in a large HDI68 for $\log P_{\max }^{\mathrm{WNL}}$, which is reconcilable with $\log P_{\max }^{\mathrm{WNE}}$ within the credible intervals.

The negative power-law index for the WNL sample indicates a preference for binary systems with shorter orbital periods. This is similar to what was found for the WNE (Paper II) and O-type populations (Sana et al. 2012). The lower cutoff for the orbital period distribution ($\log P_{\min }^{\mathrm{WNL}}$) was fixed at 0.15 [d] in order to account for binaries with periods as short as 1.6 d in our sample. The similarities with respect to the WNE population indicate that the overall multiplicity properties of the Galactic WN population are consistent. We therefore combined the WNE and WNL populations and performed a similar Bayesian analysis.

4.3. Multiplicity properties: Combined WN population

After merging the data sets of the 16 WNE stars in Paper II and the 11 WNL stars in this paper, we ended up with 27 Galactic WN stars. Based on their ΔRV values, we created Fig. 4, similar to Fig. 2, depicting $f_{\mathrm{obs}}^{\mathrm{WN}}$ as a function of the threshold C. The characteristic kink in the curve can be seen around 50 km s⁻¹, differentiating RV variation due to Doppler motion and intrinsic variability, thereby reaffirming our choice of threshold. We therefore constrained the observed binary fraction for the WN sample to be $f_{\mathrm{obs}}^{\mathrm{WN}}$ = 0.41 ± 0.09.

Fig. 4. Same as Fig. 2, but for the 27 WN stars (16 WNE + 11 WNL) in our sample. The threshold C = 50 km s−1 is shown by the vertical dot-dashed line. Entries in bold represent confirmed spectroscopic binaries in the literature.

We performed a Bayesian analysis on the combined WN sample with the same underlying assumptions for the distributions stated in Sect. 4.2. We divided the sample into various ΔRV bins as follows: ΔRV < 50 km s⁻¹ (16 objects), 50 ≤ ΔRV < 250 km s⁻¹ (6 objects), 250 ≤ ΔRV < 650 (6 objects), ΔRV ≥ 650 km s⁻¹ (no object). Similarly, we also defined minimum orbital period bins based on what is known for the sample: P  >   1 d (10 objects), P  >   10 d (3 objects), P  >   100 d (2 objects) and P  >   1000 d (1 object). The one- and two-dimensional marginalised posteriors (computed by assuming flat priors) are shown in Fig. 5. As compared to Fig. 3, the power-law index (π^WN) becomes more negative, and $f_{\mathrm{int}}^{\mathrm{WN}}$ on the WN population increases due to the contribution from the WNE sample.

Fig. 5. Same as Fig. 3, but for the combined (WNE + WNL) Galactic WN population.

The derived multiplicity properties for the WN sample were used to construct a visualisation of the period distribution. Similar to Paper II, we sampled 10⁷ sets of parameters from the posteriors in Fig. 5 and used them to create orbital period distributions. The density cloud representing these distributions is shown in red in Fig. 6. The period distribution constructed using the best-fit parameters is shown with a solid red line. We also plot the density cloud of distributions for the Galactic WC population in blue (Paper II). The period distribution for main-sequence O stars is shown with a dashed black line (Sana et al. 2012).

Fig. 6. Visualisation of the 107 period distributions created by sampling the posteriors from Fig. 5. The light and dark red regions depict 95% and 68% of all distributions. The solid red line is the distribution created with the best-fit values from the posteriors. The blue regions and lines are the same for the WC population from Paper II. The dashed black line is the period distribution for O stars from Sana et al. (2012).

5. Discussion

5.1. Comparison with the literature

In order to verify if what is observed agrees with our results, we investigated the literature to consider WR stars beyond those included our sample. Since the objects from VDH01 have been thoroughly analysed for binarity, we focused on this subset of WR stars. The catalogue contains 87 WC stars and 127 WN stars, of which 78 are WNL stars. Of these 78 objects, 21 are classified as spectroscopic binaries and 14 currently have orbital solutions. Of these, only one has a period longer than 100 days.

The total number of Galactic WN binaries with known orbital periods is 28. Combining this with the 18 known orbital periods for Galactic WC binaries, we constructed the cumulative observed orbital period distribution along with the observed distributions from our sample (Fig. 7). While a few WC binaries have periods shorter than 10 d, this fraction is still much smaller than that of short-period WN binaries. This shows that the analysed samples are representative of the Galactic WR population.

Fig. 7. Cumulative distribution of the reported orbital periods from the literature for the Galactic WN (solid red) and WC (dash-dotted blue) populations. The observed distributions based on spectroscopic periods from our sample are also plotted, with the WN (dashed red) and WC (dotted blue) binaries.

5.2. Observational biases due to the limiting magnitude

Another potential observational bias introduced by the conditions of the RV campaign is that of the limiting magnitude. Because WR stars are generally faint in the optical, the selection criteria of V ≤ 12 inherently favour binaries. Vanbeveren & Conti (1980) demonstrated this clearly for the Large Magellanic Cloud, although this effect is not dominant in the Galaxy due to the large spread in the distances.

Similar to what was done in Paper II, we assumed that the companion contributes to around 50% of the flux in a WR + OB system (e.g., Shenar et al. 2019). Therefore, WR binaries would be twice as bright as single WR stars on average. To account for this bias, we need to account for single WR stars between 12.0 and 12.7 mag. When the V-band magnitude was missing, we searched for v-band magnitudes between 13.0 and 13.7.

We found four single WNL entries in that range: WR 58 (WN4b/CE), WR 82 (WN7(h)), WR 130 (WN8(h)), and WR 131 (WN7h+absorption). Although WR 131 could be considered a candidate binary, no evidence of a companion has been found in the literature, and hence we conservatively considered it to be a single star. Including these stars in our sample would change the intrinsic binary statistics from 5 binaries out of 11 WNL stars to 5 binaries out of 15 WNL stars. This would result in a binary fraction of 0.33, which is within our credible interval.

5.3. Consequences on binary evolution

Depending on their mass, main-sequence O stars are thought to evolve into either red supergiants or luminous blue variables (LBVs) before becoming WR stars, after which they are expected to collapse to form compact objects (Conti 1976; Meynet & Maeder 2003; Crowther 2007; Langer 2012). In this section, we use our derived multiplicity properties for the Galactic WN and WC populations to infer possible connections between the evolutionary stages before collapse. As the multiplicity properties of the WNL and WNE populations are compatible, we discuss the WN population as a whole.

5.3.1. Linking the O and WN populations

As indicated in Fig. 6, the orbital period distributions of O stars and WN stars is in excellent agreement. This is not surprising because WN stars can be formed by the envelope stripping of the primary by the secondary in an OB binary. Using the analytical expressions from Soberman et al. (1997), we can calculate the change in the orbital period due to mass transfer. Depending on the mass ratio of the OB binary, mass transfer does not increase the orbital period by more than a factor of three, regardless of the conditions for mass transfer. For example, conservative mass transfer in a 40 + 30 M_⊙ binary where the primary transfers 20 M_⊙ would have an increase in the period (P_f/P_i) of a factor of 1.7. Therefore, the increase in the orbital period would not shift the distribution drastically, which is in accordance with what we find.

A different process that can significantly harden a binary is common-envelope evolution, triggered by mass transfer in systems with extreme mass ratios or after the donor has developed a deep convective envelope (see e.g., the review by Ivanova et al. 2013). However, envelope ejection is energetically disfavoured unless interaction occurs after the star becomes a red supergiant (e.g., Klencki et al. 2021), but this phase is potentially not reached by the most massive stars as it would require evolution beyond the Humphreys-Davidson limit (Humphreys & Davidson 1979; Davies et al. 2018; Gilkis et al. 2021).

Depending on their initial mass, some massive stars are also thought to go through the LBV phase during their transition from main-sequence O stars to the WN phase (see e.g., Lamers & Nugis 2002; Crowther 2007; Langer 2012). Multiplicity studies of the LBV population indicate a binary fraction of 0.70 ± 0.09, with the majority of them at periods of a year or longer (Mahy et al. 2022). The multiplicity properties of the LBVs seem inconsistent with the large number of short-period WN binaries that are found both in our sample and in the GCWR (Fig. 7). However, this could just be a consequence of shorter-period systems being too compact to fit an LBV star, leading an O star to fill its Roche lobe before reaching that evolutionary phase. This would imply that the progeny of the LBV binary population are the WN binaries with P >  1 yr.

5.3.2. Evolution from the WN to WC populations

The largest discrepancy between the orbital period distributions for both the observed (Fig. 7) and intrinsic (Fig. 6) populations is the glaring lack of WC binaries at short periods. Because WN stars are expected to evolve into WC stars by further stripping their envelope, the change in their orbital parameters is governed by wind or binary mass loss. This results in an increase of the orbital period by a factor of ∼1.5–2. Therefore, this implies that only WN binaries with P >  100 d would preferentially evolve into WC binaries. In the short-period regime of P <  10 d, only a small fraction of WN binaries might become WC binaries.

One possible factor preventing WNE stars in short-period binaries from evolving into WC binaries is that they are unable to strip their envelopes further before they reach core collapse, for example, due to reduced mass-loss rates (possible implications from Neijssel et al. 2021). Why this would only affect WN binaries with short periods and not the entire population is unclear. Another possibility is that short-period WN binary systems undergo interaction before they can reach the WC phase. The cause of this interaction might be the expansion of the WN star due to inflation (Gräfener et al. 2012; Sanyal et al. 2015; Grassitelli et al. 2018; Ro 2019). While it is unclear whether inflation occurs in nature and how an inflated donor star would respond to mass transfer, an unstable mass-transfer phase resulting in a merger could prevent short-period WN stars from evolving into WC binaries. Conversely, if these mass transfer episodes are stable, they would result in further stripping, leading to more short-period WC stars and worsening the discrepancy in the period distributions.

Additionally, in systems undergoing wind-wind collisions, slowed outflows may potentially interact with the binary and be ejected with additional angular momentum, leading to a contraction of the orbit and possibly a merger. For example, MacLeod & Loeb (2020) showed how the interaction of colliding winds with the stars in equal-mass main-sequence binaries could lead to a shrinkage of the orbit. Depending on the ratio of the orbital velocity to the wind velocity, this effect could also occur in WN binaries. Whether this scenario provides a valid explanation remains to be shown.

A final possible scenario could be that the O star companion fills its Roche lobe before the WN evolves into a WC. The likelihood of this occurring during the short-lived WNE phase is too small for this to be the predominant explanation for the lack of short-period WC systems, however.

Due to the duration and sampling of the RV campaign, it is also possible that there is a population of long-period binaries (P >  1000 d) that we do not detect. This is because the intrinsic variability of WN stars is higher than that of the WC, making it challenging to detect low RV amplitudes. To elaborate, the progenitor binary population would be composed of WNL plus mid- to late-type O star pairs. The observed RV variations for these systems across our observing campaign would be about 15 to 20 km s⁻¹. Therefore, the progenitor population of long-period WC binaries might simply be undetected. However, this would imply that for the Galactic WN population, the value of $f_{\mathrm{int}}^{\mathrm{WN}}$ is close to 1.00, as most (all?) of the apparently single WN stars would then reside in long-period binary systems.

A large fraction of the main-sequence O population are in hierarchical triple systems (e.g., Sana et al. 2014; Moe & Di Stefano 2017). Mass loss due to stellar winds or binary stripping, among other things, will change the ratio of the inner and outer separation. If the stripping of the primary by the secondary resulted in a WN + OB binary, then the nature of mass transfer and evolution of the inner orbit determines the stability of the triple. When this ratio changes beyond a critical point, the system becomes dynamically unstable. This could result in the disruption of the triple, or in a merger of the inner binary (Toonen et al. 2020).

In the case of mass loss through stellar winds, the dynamically unstable phase can last from thousands to millions of years (Toonen et al. 2022). Although the most common predicted outcome is the preservation of the triple hierarchy, it is also common for the system to be disrupted, with one of the stars ejected. In case of an ejection, the resulting orbit has a period of the order of 10³–10⁵ d (Toonen et al. 2022). As an example in a hierarchical triple, if the inner WN + OB binary was formed through a combination of stripping via stellar winds and mass transfer, it could be dynamically unstable. An ejection could result in the formation of a WN + OB binary in a wide orbit, where the WN star would proceed to evolve and become a wide WC + OB binary. Detailed simulations involving the formation of the WR star, mass transfer, and stellar winds are required to investigate the frequency and feasibility of this channel.

6. Conclusions

With a homogeneous, magnitude-limited spectroscopic monitoring of 11 northern Galactic WNL stars, we measured relative RVs using cross-correlation and found $f_{\mathrm{obs}}^{\mathrm{WNL}}$ =0.36  ±  0.15. Using a Bayesian framework around MC simulations, we derived a fraction corrected for observational biases ($f_{\mathrm{int}}^{\mathrm{WNL}}$ $=0.{62}_{-0.17}^{+0.16}$). We observed that the computed intrinsic multiplicity properties for the Galactic WNL population are congruous with those derived for the WNE population in Paper II.

We combined the WNE and WNL samples and simulated ∼4  ×  10⁸ populations of 27 WN stars. We found the intrinsic multiplicity properties for the Galactic WN population to be similar to those derived for main-sequence O stars in Sana et al. (2012). We constructed the orbital period distribution for Galactic WNs and noted the peak at P <  5 d. The abundance of short-period WN binaries can be explained through case A or case B mass transfer in O binaries. Short-period WN binaries can also form via common envelope evolution through late case B or case C mass transfer, regardless of the mass ratio. However, it is observationally non-trivial to distinguish between the multitude of channels.

From the intrinsic multiplicity properties of the Galactic WC population (Paper I), we see a clear lack of short-period WC binaries, with the orbital period distribution peaking at P  ∼  5000 d. The lack of short-period WC binaries is also supported by the literature. While a few are known, their number is smaller than for their WN counterparts, even though the number of WNs and WCs in the Galaxy is similar. The discrepancy in distributions is not only found in our simulations (Fig. 6), but can also be seen in the observed period distribution of Galactic WR binaries (Fig. 7). Orbital evolution via mass loss from the WN to WC phase cannot explain the shift towards longer periods.

The multiplicity properties of Galactic LBVs (Mahy et al. 2022) appear to agree with that of the long-period O, WN, and WC binary populations. This might imply that at long periods, an evolutionary channel of O (→LBV) → WN → WC might exist. However, the evolution of the short-period WN binaries is still an enigma. Invoking exotic channels such as common-envelope evolution and dynamical interactions in triples could explain their evolution, but these channels are highly uncertain. Further observational and theoretical studies are required to understand their importance in the grand scheme of binary evolution.

As the immediate progenitors of black holes, the multiplicity properties of WR stars directly affect those for black-hole binaries (e.g., through black-hole kicks in a core-collapse scenario), and hence the predictions for gravitational-wave progenitors. Expanding the sample size of magnitude-limited studies is critical for understanding the binary evolution of massive stars.

Acknowledgments

This work was published with funds from the European Research Council under European Union’s Horizon 2020 research programme (grant agreement No. 772225). T.S. also acknowledges funding received from the European Union’s Horizon 2020 under the Marie Skłodowska-Curie grant agreement No. 101024605. P.M. acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N.

References

Appendix A: Comments on individual objects

WR 108: According to the GCWR, the spectra of WR 108 have diluted emission lines and show absorption. While it is quite common for WNL stars to show absorption, diluted emission lines are often attributed to a companion star. We measured a ΔRV of ∼16 km s⁻¹ over 1140 days, and classified WR 108 as a single star.

WR 120: The GCWR classifies WR 120 as a single WN7 star. The few absorption lines in the spectra exhibit moderate line-profile variability. However, the S/N of the data in the extreme blue regions makes it challenging to identify whether this is simply due to intrinsic variability or due to Doppler motion. With seven epochs spanning 385 days, we find a ΔRV of 40.5 km s⁻¹ using N III, and classified it as a single star.

WR 123: WR 123 is notorious for its strong photometric and spectral line-profile variability (Marchenko & Moffat 1998). It was found to have a periodicity of ∼10 hours using photometric and spectroscopic monitoring (Lefèvre et al. 2005; Dorfi et al. 2006; Chené et al. 2011). Similar to WR 120, we find the absorption lines in the spectra to vary quite significantly. By monitoring N Vλλ4604, 4620 for eight epochs over 387 days, we measured a ΔRV of 57.1 km s⁻¹ and classified it as a binary system. The sampling and number of spectra are too sparse to find any significant periodicity, and so further monitoring is required to verify its binary status.

WR 124: The GCWR classifies WR 124 as a SB1? with no dilution of emission lines. It is surrounded by a spectacular nebula M1-67 and is the fastest runaway in the Galaxy (Moffat et al. 1982). Given its 2.4 d and 2.7 d periodicity (Moffat et al. 1982; Moffat & Shara 1986), it was considered to be a binary with a neutron star. Toalá et al. (2018) were unable to rule out the presence of an embedded neutron star when studying the system with X-rays. If this were the case, the RV amplitude of the WR star would be below our threshold. In any case, we measured ΔRV to be 19.3 km s⁻¹ with ten epochs over 577 days, and classified WR 124 as a single star.

WR 134: Similar to WR 123 and WR 124, WR 134 is classified as a SB1? with no dilution of emission lines. Its binary status has been widely debated in the literature, and the observed periodicity has been attributed to a low-mass companion (Marchenko & Moffat 1998; Rustamov & Cherepashchuk 2012). However, Aldoretta et al. (2016) demonstrated that large-scale structures in the wind, particularly in He IIλ5412, exist with a period of 2.255 ± 0.008 d. We obtained 15 epochs over 2901 days and measured ΔRV to be 29.6 km s⁻¹ using N V λ4945 and classified it as a single star.

WR 136: WR 136 is surrounded by a ring nebula (NGC 6888), and is classified as a SB1? with no dilution of emission lines in the GCWR. A period of 4.554 d was found using line-profile variations (Koenigsberger et al. 1980) and RV measurements (Aslanov & Cherepashchuk 1981). However, this period was not verified later with polarisation (Robert & Moffat 1989) and photometry (Moffat & Shara 1986). A neutron star companion was suggested as a companion, which would again result in small RV variations of the WR star. With 39 epochs over a time baseline of 2907 days, we measured a ΔRV of 15.4 km s⁻¹ using He II lines and classified it as a single star.

WR 148: WR 148 is one of the few WN8 stars in the Galaxy that is confirmed to be in a binary system. The GCWR classifies it as a SB1 system, as it has a short period of 4.3 d and was thought to have a compact companion (Drissen et al. 1986a; Bracher 1979; Marchenko et al. 1996). However, Munoz et al. (2017) found a period of 4.317336 ± 0.000026 and were able to disentangle the spectra to find a O4-6 companion. Hamann et al. (2019) found a luminosity (logL (L_⊙)) of 6.2, indicating that it might be a main-sequence star. However, they also reported that 15% of the flux is contributed by its companion, so that further analysis taking this into account is required. We obtained 20 epochs over 1161 days and measured RVs using N IVλλλ7103, 7109, 7123. We found WR 148 to have a ΔRV of 171.8 km s⁻¹ and classified it as a binary system.

WR 153: The GCWR classifies WR 153 to be a quadruple system (WN6o/CE+O3-6 + B0:I+B1:V-III) where the WR binary has a period of 6.6887 days and a circular orbit (Massey 1981; Demers et al. 2002). With 49 epochs over 2858 days, we measured a ΔRV of 576.6 km s⁻¹ using N Vλλ4604, 4620 and classified it as a binary system.

WR 155: The GCWR classifies WR 155 as a SB2 (WN6o + O9II-Ib) system. This binary has the shortest known period (1.6412436 d; Moffat & Shara 1986; Marchenko et al. 1995). Drissen et al. (1986b) studied the system with polarimetry and found masses of 42 and 30 M_⊙ for the WN and O stars, respectively. We obtained 85 epochs over 3075 days, measured ΔRV to be 609.5 km s⁻¹ using N V λ4945 and λ4058 and classified it as a binary system.

WR 156: WR 156 is another WN8 star with large photometric and spectral line-profile variability, with diluted emission lines. Photometric studies did not find any conclusive periodicity (Moffat & Shara 1986; Marchenko et al. 1998). It is reported to have a luminosity of 6.01 (Hamann et al. 2019), which questions its status as a massive WNL versus a classical WNL star. We observed WR 156 over 797 days and obtained 18 epochs. Using N IVλλλ7103, 7109, 7123 and λ5737 lines, we measured a ΔRV of 16.1 km s⁻¹ and classified it as a single star.

WR 158: The GCWR classifies WR 158 as a WN7h star with diluted emission lines. Andrillat & Vreux (1992) discovered O Iλ8446 in its spectrum, hypothesising that it is produced in the slowly expanding cool shocked region located outside the ionisation front of the bubble created by the WR wind, or in a binary with a Be companion. Hamann et al. (2019) reported a luminosity (logL (L_⊙)) of 6.06, which might indicate its status as a main-sequence WNL star. We obtained 13 spectra over 1031 days, measured a ΔRV of 15.0 km s⁻¹ using N IVλλλ7103, 7109, 7123 and classified it as a single star.

Appendix B: Relative radial velocity measurements

Relative RVs for the 11 WNLs in our sample. The reference epoch has a RV of 0.0 km s⁻¹. We refrain from providing absolute RV measurements as this is highly method dependent, especially for WR stars. The barycentric Julian date (BJD) is given as the middle of the exposure. The average S/N is given in Table 1. Along with the measured RVs, we indicate the measurement error, that is, the statistical uncertainty σ_p (Eq. 9 in Paper I).

All Tables

All Figures

	Fig. 1. High-cadence RV measurements for WR 136 over a variety of masks. In brackets at the bottom are values of (σw, σRV) (km s−1). As there are 87 data points for each mask, the legend has been omitted for clarity.
In the text

	Fig. 2. Non-parametric threshold plot showing the observed binary fraction as a function of the threshold C. The adopted threshold of 50 km s−1 is shown by the dot-dashed vertical line. The entries in bold are confirmed spectroscopic binaries in the literature.
In the text

	Fig. 3. Marginalised posteriors of $\log P_{\max }^{\mathrm{WNL}}$, $f_{\mathrm{int}}^{\mathrm{WNL}}$, and πWNL for the Galactic WNL population. The one-dimensional posteriors are calculated assuming flat priors. For each posterior, the solid red lines show HDI68, and the dashed red line shows the mode.
In the text

	Fig. 4. Same as Fig. 2, but for the 27 WN stars (16 WNE + 11 WNL) in our sample. The threshold C = 50 km s−1 is shown by the vertical dot-dashed line. Entries in bold represent confirmed spectroscopic binaries in the literature.
In the text

	Fig. 5. Same as Fig. 3, but for the combined (WNE + WNL) Galactic WN population.
In the text

Fig. 6. Visualisation of the 107 period distributions created by sampling the posteriors from Fig. 5. The light and dark red regions depict 95% and 68% of all distributions. The solid red line is the distribution created with the best-fit values from the posteriors. The blue regions and lines are the same for the WC population from Paper II. The dashed black line is the period distribution for O stars from Sana et al. (2012).

In the text

	Fig. 7. Cumulative distribution of the reported orbital periods from the literature for the Galactic WN (solid red) and WC (dash-dotted blue) populations. The observed distributions based on spectroscopic periods from our sample are also plotted, with the WN (dashed red) and WC (dotted blue) binaries.
In the text