A spectroscopic multiplicity survey of Galactic Wolf-Rayet stars. III. The northern late-type nitrogen-rich sample

Massive stars are powerful cosmic engines. In the phases immediately preceding core-collapse, massive stars in the Galaxy with $M_i \gtrsim 20$ $M_{\odot}$ may appear as classical Wolf-Rayet (WR) stars. As the final contribution of a homogeneous RV survey, this work constrains the multiplicity properties of northern Galactic late-type nitrogen-rich Wolf-Rayet (WNL) stars. We compare their intrinsic binary fraction and orbital period distribution to the carbon-rich (WC) and early-type nitrogen-rich (WNE) populations from previous works. We obtained high-resolution spectra of the complete magnitude-limited sample of 11 Galactic WNL stars with the Mercator telescope on the island of La Palma. We used cross-correlation to measure relative RVs and flagged binary candidates based on the peak-to-peak RV dispersion. By using Monte Carlo sampling and a Bayesian framework, we computed the three-dimensional likelihood and one-dimensional posteriors for the upper period cut-off, power-law index, and intrinsic binary fraction. Adopting a threshold $C$ of 50 km s$^{-1}$, our Bayesian analysis produced an intrinsic fraction of $0.42\substack{+0.15 \\ -0.17}$ for the parent WNL population alongside distributions for the power-law index and the orbital periods. The observed period distribution of Galactic WN and WC binaries from the literature is in agreement with what is found. The period distribution of Galactic WN binaries peaks at $P{\sim}1$-$10$d and that of the WC population at $P{\sim}5000\,$d. This shift cannot be reconciled by orbital evolution due to mass loss or mass transfer. At long periods, the evolutionary sequence O($\xrightarrow{}$LBV)$\xrightarrow{}$WN$\xrightarrow{}$WC seems feasible. The high frequency of short-period WN binaries compared to WC binaries suggests that they either tend to merge or that the WN components in these binaries rarely evolve into WC stars in the Galaxy.


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. 2019Shenar et al. , 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 corecollapse 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 Article number, page 1 of 15 arXiv:2212.06927v1 [astro-ph.SR] 13 Dec 2022 A&A proofs: manuscript no. mainfile 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 1 (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. 2009Shara et al. , 2012Faherty 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 1). We determined the observed binary fraction ( f WC obs ) to be 0.58 ± 0.14. Using Monte Carlo (MC) simulations with the measured RV uncertainties, we derived the intrinsic binary fraction ( f WC int ) of the Galactic WC population to be at least 0.72 at the 10% significance limit with orbital periods up to 10 5 d. The 16 early-type WN (WNE) stars were analysed in Dsilva et al. (2022, henceforth Paper 2), where the observed binary fraction ( f WNE obs ) 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 WNE int ) to be 0.56 +0.20 −0.15 with orbital periods explored up to 10 5 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 2 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 WC int = 0.96 +0.04 −0.22 and π WC = 1.90 +1.26 −1.25 ), 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 (postmain-sequence, hydrogen-rich) or very massive main-sequence (hydrogen-burning) WR stars (e.g. ?). 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 mainsequence and classical WNL stars. When the known binaries are excluded, the luminosities of the stars in our sample are lower than or around 10 6 L (Hamann et al. 2019). We aim to increase 1 http://pacrowther.staff.shef.ac.uk/WRcat/  N  108  11  1140  125  120  7  385  30  123  8  387  60  124  10  577  50  134  15  2901  170  136  39  2907  280  148  20  1161  110  153  49  2858  200  155  85  3075  140  156  18  797  70  158  13  1031  70 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.

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 narrowband magnitudes (Smith 1968;Massey 1984) and applied the criterion ≤ 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 1. In short, the data were first reduced using the standard HER-MES pipeline (Raskin et al. 2011) and then corrected for the instrumental response (Paper 1, 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 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 pseudocontinuum 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).

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

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

Wind variability and its effect on radial velocity measurements
Comparing the WNL stars in our sample to the WCs (Paper 1) and WNEs (Paper 2), we observe significantly more variability in the WNL line profiles. This is expected because strong lineprofile 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 1 and Paper 2, 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 −1 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;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.  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 Table 2. Radial velocity measurements for our sample of WNL stars together with an overview of their known multiplicity properties. ∆ RV and σ RV are calculated in this work and are used to identify tentative spectroscopic binaries. The spectral types are taken from the GCWR unless indicated otherwise, with 'd.e.l.' implying the dilution of emission lines, and 'a' indicating the presence of absorption lines. The binary status of this work is reported based on the spectroscopic observations.  Demers et al. (2002), (c) Marchenko et al. (1995) fixed period (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

WR#
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 2). 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.

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.
In principle, the chosen threshold should be high enough to avoid false positives due to intrinsic variability. The highcadence 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 −1 . This is similar to what was found for the N v lines in WR 138 (Paper 1). 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 −1 . 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 WNL obs ) of 0.36 ± 0.15.

Multiplicity properties: WNL population
We performed a Bayesian analysis with MC simulations with a method developed in Paper 2 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 −1 (seven objects), 50 ≤ ∆RV < 250 km s −1 (two objects), 250 ≤ ∆RV < 650 km s −1 (two objects), ∆RV ≥ 650 km s −1 (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 WNL min ), maximum period (log P WNL max ), and the powerlaw index (π WNL ). The power-law index π determines the shape of the distribution of binaries with orbital periods P as follows: (2) Along with the period distribution, we included the intrinsic binary fraction, f WNL int , 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 threedimensional 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 2 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 WNL min at 0.15 [d] (approximately 1 d). We explored values of log P WNL max 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 WNL int 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 2 grid point, resulting in almost 4×10 8 populations. We calculated the three-dimensional likelihood of reproducing the observed distribution of objects in ∆RV bins. Assuming flat priors on log P WNL max , π WNL and f WNL int , we computed their marginalised posterior likelihoods (Fig. 3). We computed the 68% credible intervals for each posterior (the highestdensity interval, HDI68 The value of log P WNL max 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.
The derived multiplicity properties agree with those of the Galactic WNE population (Paper 2). An interesting difference is that of the upper-period cutoff, which is better constrained for the WNE population (log P WNE max = 4.60 0.40 −0.77 ) 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 WNL max , which is reconcilable with log P WNE max 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 2) and O-type populations (Sana et al. 2012  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.

Multiplicity properties: Combined WN population
After merging the data sets of the 16 WNE stars in Paper 2 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 WN obs as a function of the threshold C. The characteristic kink in the curve can be seen around 50 km s −1 , 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 WN obs = 0.41 ± 0.09. 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 −1 (16 objects), 50 ≤ ∆RV < 250 km s −1 (6 objects), 250 ≤ ∆RV < 650 (6 objects), ∆RV ≥ 650 km s −1 (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 WN int on the WN population increases due to the contribution from the WNE sample.
The derived multiplicity properties for the WN sample were used to construct a visualisation of the period distribution. Similar to Paper 2, we sampled 10 7 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 bestfit parameters is shown with a solid red line. We also plot the density cloud of distributions for the Galactic WC population in blue (Paper 2). The period distribution for main-sequence O stars is shown with a dashed black line (Sana et al. 2012).

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.

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 2, 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 -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.

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.

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 K. Dsilva , T. Shenar, H. Sana and P. Marchant: A spectroscopic multiplicity survey of Galactic Wolf-Rayet stars 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.

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 shortperiod 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 −1 . 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 WN int 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 3 -10 5 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.

Conclusions
With a homogeneous, magnitude-limited spectroscopic monitoring of 11 northern Galactic WNL stars, we measured relative RVs using cross-correlation and found f WNL obs = 0.36 ± 0.15. Using a Bayesian framework around MC simulations, we derived a fraction corrected for observational biases ( f WNL int = 0.62 +0.16 −0.17 ). We observed that the computed intrinsic multiplicity properties for the Galactic WNL population are congruous with those derived for the WNE population in Paper 2.
We combined the WNE and WNL samples and simulated ∼4 × 10 8 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 shortperiod 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 1), 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.