A spectroscopic multiplicity survey of Galactic Wolf-Rayet stars: II. The northern WNE sequence

Most massive stars reside in multiple systems that will interact over the course of their lifetime. Classical Wolf-Rayet (WR) stars represent the final end stages of stellar evolution at the upper-mass end. As part of a homogeneous, magnitude-limited ($V\leq12$) spectroscopic survey of northern Galactic WR stars, this paper aims to establish the observed and intrinsic multiplicity properties of the early-type nitrogen-rich WR population (WNE). We obtained high-resolution spectroscopic time series of the complete magnitude-limited sample of 16 WNE stars observable with the 1.2 m Mercator telescope at La Palma, typically providing a time base of about two to eight years. We measured relative radial velocities (RVs) using cross-correlation and used RV variations to flag binary candidates. Adopting a peak-to-peak RV variability threshold of 50 km/s as a criterion found an observed multiplicity fraction of 0.44$\pm$0.12. Using an updated Monte Carlo method with a Bayesian framework, we calculated the three-dimensional likelihood for the intrinsic binary fraction, the maximum period, and the power-law index for the period distribution for the WNE population. We also re-derived multiplicity parameters for the Galactic WC population. We found an intrinsic multiplicity fraction of $0.56\substack{+0.20 \\ -0.15}$ for the parent WNE population. For the Galactic WC population, we re-derive an intrinsic multiplicity fraction of $0.96\substack{+0.04 \\ -0.22}$. The derived multiplicity parameters for the WNE population are quite similar to those derived for main-sequence O binaries but differ from those of the WC population. The significant shift in the WC period distribution towards longer periods is too large to be explained via expansion of the orbit due to stellar winds, and we discuss possible implications of our results.


Introduction
Stars with initial masses M i 8 M will end their lives through core-collapse supernovae or through direct collapse (Heger et al. 2003), leaving behind compact objects such as neutron stars or black holes. Through their strong stellar winds and energetic end-Based on observations made with the Mercator Telescope, operated on the island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.
Based on observations obtained with the HERMES spectrograph, which is supported by the Research Foundation -Flanders (FWO), Belgium, the Research Council of KU Leuven, Belgium, the Fonds National de la Recherche Scientifique (F.R.S.-FNRS), Belgium, the Royal Observatory of Belgium, the Observatoire de Genève, Switzerland and the Thüringer Landessternwarte Tautenburg, Germany.
of-life explosions, these stars dominate the transport of energy and momentum into the interstellar medium (see Mac Low et al. 2005;Crowther 2019).
At the upper mass end, stars with M i 20 M may enter the classical Wolf-Rayet (WR) phase, where after stripping (most of) their hydrogen-rich envelope, they reveal deep layers that are strongly enriched by products of nucleosynthesis (for a review, see Crowther 2007). The spectra of WR stars are characterised by strong, broad emission lines that result from strong, optically thick, radiation-driven winds. WR stars are generally understood to be the evolved, core He-burning counterparts of massive O stars and are immediate progenitors of black holes. WR stars play an important role in constraining the evolution of massive stars (Vanbeveren et al. 1998;Meynet & Maeder 2003;Shenar et al. 2016;Hamann et al. 2019;Sander et al. 2019 A&A proofs: manuscript no. aanda as well as the properties of supernovae and compact objects (de Mink et al. 2014;Marchant et al. 2016;Hainich et al. 2018;Langer et al. 2020).
WR stars are divided into three main spectral classes: the nitrogen (WN), carbon (WC), and oxygen (WO) sequences, which reflect the nitrogen-, carbon-, or oxygen-rich chemical composition of their atmospheres. Conti (1976) first proposed that over the course of their main-sequence evolution, O stars that are massive enough may lose sufficient mass via stellar winds to reveal, first, the products of H-burning, followed by those of He-burning (see also Rublev 1965). These products are nitrogen and carbon, respectively, and the evolutionary counterparts are thought to be WN and WC stars. Ever since, this picture of evolution has been called the "Conti scenario". In reality, there is a fluid transition between O and WN stars (Of, Of/WN, e.g. Conti et al. 1995;Crowther & Bohannan 1997), as well as WN and WC stars (WN/WC stars, Conti & Massey 1989).
Wolf-Rayet stars can also be produced in close binary systems where the envelope of the star is stripped through Roche-lobe overflow by a companion, or through a combination of stellar winds and mass transfer (Paczyński 1967). The dominant formation channel for WR stars as a function of mass and metallicity is still widely debated (Vanbeveren et al. 1998;Neugent & Massey 2014;Shenar et al. 2019Shenar et al. , 2020. As the immediate progenitors of black holes and neutron stars, the characteristics of WR stars directly determine the observed properties of supernovae (Zapartas et al. 2017;Eldridge et al. 2018) as well as the expected distribution of gravitational wave systems (Eldridge et al. 2019).
Radial velocity (RV) surveys in the last few decades have shown that the majority of main-sequence O stars are in multiple systems with periods of up to a few years, with a significant fraction at periods shorter than a few tens of days (Mason et al. 2009;Sana et al. 2012Sana et al. , 2013Sana et al. , 2014Kobulnicky et al. 2014; Barbá et al. 2014;Maíz Apellániz et al. 2016;Almeida et al. 2017;Maíz Apellániz et al. 2019). As they are likely to interact during their lifetime, comparing the multiplicity properties of O stars to that of WR stars during the WN and WC phases is crucial for calibrating the physics of binary interaction. Moreover, orbital analyses of WR binary systems allow us to measure black hole progenitor masses in a model-independent way.
The Galactic Catalogue of WR stars 1 (henceforth GCWR) is maintained by Paul Crowther (originally, Appendix 1 in Rosslowe & Crowther 2015), with 667 objects currently (v1.25). Until recently, the observed multiplicity fraction of WR stars was reported to be ∼0.4 in the literature (van der Hucht 2001, henceforth VDH01). However, there are multiple caveats to this value. Firstly, this fraction includes photometric, spectroscopic (with and without derived orbits), and visual WR binaries. Secondly, the number of WRs with known orbits remains small with respect to the number of WRs listed as binaries in the GCWR. Thirdly, VDH01 is a compilation of decades of work and involves different techniques, each with their own sensitivity and limitations. This means that correcting for observational biases is non-trivial. As a consequence, the intrinsic multiplicity fraction of WR stars and the distributions of orbital elements and mass ratios are insufficiently constrained. These are important resources that can be used to calibrate the physics in binary stellar evolution codes (e.g. efficiency of mass transfer and angular momentum transport) as well as recipes used in stellar population synthesis (e.g. BPASS, Eldridge et al. 2017).
In order to obtain a reliable correction for observational biases, accurate RV measurements, and their uncertainties, are of  Dsilva et al. (2020, henceforth Paper 1), the survey targets all WR stars brighter than V = 12 that are visible from La Palma. In Paper 1, we focused on the multiplicity properties of the WC population. Based on the 12 WC stars in our sample, we determined an observed binary fraction ( f WC obs ) of 0.58 ± 0.14, where uncertainties due to the sample size have been taken into account. Using Monte Carlo (MC) simulations and the measured RV uncertainties, we deduced the intrinsic multiplicity fraction ( f WC int ) of the Galactic WC population to be at least 0.72 at the 10% significance limit. Moreover, most of these systems had ∆RV values that are indicative of orbital periods longer than a few hundred days, which is consistent with the small fraction (∼10 %) of known short-period WC binaries in VDH01.
As a second paper in this series, we focus on 16 northern Galactic early-type WN (WNE) stars. Because WNE stars are believed to be the progenitors of WC stars, a direct connection between their multiplicity properties is expected. The sample, observations, and data reduction are presented in Sect. 2. Section 3 elaborates on the RV measurements and relevant masks used for cross-correlation. In Sect. 4 we discuss the spectroscopic binary fraction and compare it to the literature. In Sect. 5 we discuss the correction for observational biases and the intrinsic binary fraction. In Sect. 6 we discuss the evolutionary consequences. Section 7 presents our conclusions.

Sample
This work focuses on the WNE stars in our sample, which we define as WN stars with spectral types equal to or earlier than WN5. In cases where objects are classified as 'WN5-WN6' in the GCWR, we consider them as WN5 for simplicity. From the GCWR we selected the abovementioned stars with V ≤ 12 and declination δ ≥ −30 • . For objects without a V -band magnitude, we used narrow-band magnitudes (Smith 1968;Massey 1984), with the condition ≤ 13. This resulted in a magnitude-limited sample of 16 WNE stars whose known properties are summarised in Table 2. Their spectral-type distribution is shown in Fig. 1. Observations were carried out with the HERMES spectrograph mounted on the 1.2 m Mercator telescope at the Roque de los Muchachos observatory in La Palma, Spain. HERMES covers the wavelength range from 3800 Å to 9000 Å with a spectral resolving power of R = λ/∆λ ∼ 85000. We obtained at least six epochs of each of the target stars. When available, we used archival HERMES observations in the RV analysis, which reach a typical time base of two to eight years. The number of spectra and the sampling time of the observations are shown in Table 1.

Data reduction and normalisation
The data reduction procedures for our spectra are comprehensively described in Paper 1. After the standard reduction from the HERMES pipeline (Raskin et al. 2011), a correction for the instrumental response function is applied. Molecfit Kausch et al. 2015) is then used to correct the spectra for telluric contamination. After this step, what remains is the effect of interstellar reddening superimposed on the continuum from the wind. The spectra are then normalised in a homogeneous way that minimises its effect on the RV measurements. We used a model of a normalised spectrum of a WNE star from the Potsdam Wolf-Rayet code (PoWR: Gräfener et al. 2002;Hamann & Gräfener 2003, 2004Todt et al. 2015) to identify pseudo-continuum regions, in the blue around 5300 Å and in the red around 8100 Å. We then anchored the stellar continuum model for the same WNE star to the red point and applied a reddening to fit the slope. We explored various reddening laws and R V values and realised that they did not significantly impact the homogeneity of the normalisation process. Therefore, we used the reddening law from Fitzpatrick (2004) and R V of 3.1, the average for the galaxy.

Measuring radial velocities
Deriving RVs of WR stars is not a trivial process given their strong, broad emission lines that form in their outflowing winds. Of the methods that are used extensively in the community, we used cross-correlation as our tool of choice, with a statistical framework that allows us to derive meaningful uncertainties on our RV measurements (Zucker 2003, see Paper 1 for more details). As mentioned in Paper 1, we focus on the RVs of the WR star even in SB2 systems to maintain homogeneity in the analysis and to focus on binary detection instead of orbital analysis. We briefly describe the formulation and assumptions of the cross-correlation function (CCF).

Cross-correlation and mask selection
The CCF is the convolution of a chosen template to the data to measure the velocity shift. It is represented by a log-likelihood function, where the maximum is the RV of the spectrum. The uncertainty on the RV measurement can then be derived from maximum log-likelihood theory. An implicit assumption of the method is that the template accurately reproduces the data. Because WR models often do not to accurately represent the shape of both strong and weak lines simultaneously, we chose the highest signal-to-noise (S/N) epoch as the template. The measured RVs were then used to construct a weighted co-added spectrum of a higher S/N. This process was iterated until the derived RVs no longer changed significantly within the measurement errors, at which point the final template and measurements were converged. The S/N of the template contributes significantly to the measured uncertainty, and this was minimised by creating a high S/N coadded spectrum. A caveat is that the RVs derived in this fashion are relative RVs. While this is sufficient for binary detection, an absolute shift must be applied when comparing to RVs from the literature.
When compared with the WC stars analysed in Paper 1, the WNE stars in our sample exhibit stronger line-profile variability in their spectra. Therefore, the assumption that the template accurately represents the data often breaks down for variable stars. This directly results in larger uncertainties on the derived RVs. In order to ensure that the templates can represent the data as accurately as possible, for each star we selected encompassing spectral lines that are least affected by line-profile variability. Depending on how variable the stars are, we used lines of N iv, N v, He ii, and when possible, Balmer lines. In cases of strong line-profile variability, we tend to only use N v lines and occasionally N iv lines because they are thought to form in the hotter, inner regions of the stellar wind. We used a combination of all the different lines of various strengths and ions used to measure RVs and call this mask 'full spec' (e.g. a combination of N v, N iv, and He ii lines). At the bottom of the plot are the values of (σ w , σ RV ) (km s −1 ) for each mask. As mentioned in the text, 'full spec' is a combination of all the masks to its left and not the entire spectrum. The N iv line at 5205 Å is very weak, which explains its high scatter.

Studying the long-period binary system WR 138 to quantify wind variability
In Paper 1, we analysed high-cadence data of the long-period binary system WR 137 in order to understand the effect of wind variability and other systematic effects on the RV measurements. As a nitrogen-rich analogue, we studied the long-period binary WR 138 (WN5o + O9V) with a period of 4.16 yr and a low eccentricity of 0.3 (Palate et al. 2013;Annuk 1990;Richardson et al. 2016). Palate et al. (2013) also studied the system in X-rays, concluding that its emission is normal for a colliding-wind WR + OB binary system, although the sampling of their data did not allow them to check for the existence of phase-locked variability of the X-ray emission. With this in mind, we collected high-cadence HERMES data in August of 2019 (with the periastron passage in June 2020). In total, we obtained 18 spectra over 16 consecutive nights and 6 spectra over one night a few days later (Table B.11). The phased RV plot of WR 138 is shown in Fig. 2 with an emphasis on the high-cadence data over different masks. Because the high-cadence data were collected at the orbital phase of 0.8, the effect of windwind collision (e.g. Lührs & Moffat 2002) on the measured RVs could be significant. Using the systemic and stellar parameters (Annuk 1990;Richardson et al. 2016), we computed the cooling parameter χ (Stevens et al. 1992) and found the shocked region to be adiabatic. As the cooling timescale is longer than the escape timescale, the contribution to recombination lines in the optical is minimal. Therefore, the contribution of wind-wind collision to the RV measurements should not change significantly over the timescale of the high-cadence observations. Even if wind-wind collision were to contribute significantly to the RV measurements, the derived scatter in the high-cadence study would be an upper limit on the intrinsic wind variability of the WR star.
For the high-cadence data, the contributions to the uncertainties in the measured RVs are statistical (i.e. from the method, σ p ) and physical (i.e. from the wind physics, σ w ). The weighted stan-  Fig. 4. Non-parametric threshold plot showing the inverse cumulative distribution of peak-to-peak RV variability in our sample, which is equivalent to the observed binary fraction ( f WNE obs ) as a function of the adopted peak-to-peak RV variability threshold C. Names in bold are known binaries with established spectroscopic periods for the WR components. dard deviation of the high-cadence measurements are indicative of the true error, which we denote by σ RV (eq. 12 in Paper 1). Assuming σ p and σ w are uncorrelated, they can simply be added quadratically as For each mask, σ p is calculated from the co-variance matrix of the CCF (eq. 9 in Paper 1). Using Eq. 1 along with the calculated values of σ p and σ RV , we obtain an estimate for σ w . The high-cadence data therefore allow us to isolate and quantify the effect of wind variability on the RV measurements. The lines that are least affected by wind variability are N v weak lines (App. A, σ w ∼ 5 km s −1 ). Lines of He ii and N iv lines are found to be more variable (σ w ∼ 10 km s −1 ). As their formation regions can be spatially different, these measurements can be helpful probes of turbulence and other physical processes. The difference in dispersion for lines of N iv arise from the different line strengths and since the S/N in the blue is much lower than in the red. The mask 'full spec' is a combination of the N iv, He ii, and N v lines. The mask 'N4 7109' covers the triplet N iv λλλ7103, 7109, 7123. Ultimately, we chose to measure RVs for WR 138 using weak N v lines, namely a combination of N v λλ4604, 4620, and N v λ4945.
Although we observed that weak N v lines showed the least intrinsic variability for WR 138, this cannot necessarily be generalised for all WNE stars. However, for all the objects in our sample, we realised that the N v lines indeed had the least RV dispersion, and we therefore used them to measure RVs. For stars with strong line-profile variability (WR 2, WR 6, WR 7, and WR 110), we used the combination of all masks ('full spec') to measure RVs (App. A). The full journal of RV measurements for each object with their respective masks is given in App. B.

Observed binary fraction
For each star in our sample, we measured the RVs and obtained their peak-to-peak variability amplitude (∆RV). We then selected Table 2. Overview of the known multiplicity properties of our sample of WNE stars together with the results of our RV measurements. ∆ 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. The binary status of this work is reported based on the spectroscopic observations.

WR# Spectral Type Binary Status
WN5h VB VB (l) --9.9 3.6 no 110 WN5b WN4o   a threshold to determine the boundary beyond which the RV amplitude would be dominated by orbital motion in a binary system. We plot the binary fraction as a function of this threshold in a non-parametric threshold plot (Fig. 4). By moving the threshold from higher to lower RV values (right to left in Fig. 4), we classify more stars as RV variable, hence putative spectroscopic binaries. A summary of the characteristics of the WNE stars in our sample can be found in Table 2. As we cannot verify the status of objects classified as visual binaries (VB) in the GWCR, we simply carry over their status in our classification.
The next step was to choose an adequate value for the threshold C such that stars with ∆RV≥ C were classified as binaries. It is important that the chosen threshold avoids false positives, in particular, those that are erroneously classified as binaries due to strong intrinsic variability. The high-cadence study of WR 138 is a useful guideline to ensure this. The peak-to-peak RV amplitude measured due to wind variability was observed to be 15 km s −1 for the N v lines (Fig. 3). We therefore chose a conservative threshold larger than three times this value, that is, 50 km s −1 . This results in 7 of the 16 observed binaries in the sample, f WNE obs = 0.44 with a binomial error of ± 0.12. Interestingly, all the known spectroscopic binaries in our sample with an established orbital solution are properly identified with the adopted threshold. For objects with ∆RV < 50 km s −1 and a large enough RV time series (WR 1 and WR 2), we investigated possible periodicities using a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982), but found no coherent periodicity.

Detection probability
In order to determine the intrinsic multiplicity fraction of Galactic WNE stars, f WNE int , it is important to understand the observational biases of the RV campaign. The number of detected binaries is merely a minimum estimate of the true number of WNE binaries in our sample. Physical and orbital properties of the systems, geometrical effects, time coverage, and sampling of the RV study all affect the binary detection probability. As in Paper 1, we modelled the detection probability, estimated the number of undetected binaries, and constrained some of the properties of the parent WNE population.
Using MC simulations, we determined the sensitivity of our survey with the method explained in Patrick et al. (2019). We constructed a grid with different orbital periods (P) and secondary masses (M 2 ) between 1 and 10 5 d and 1 and 40 M , respectively. The grid was evenly spaced with steps of ∆ log(P/d) = 0.1 and ∆M 2 = 1 M . For each grid point of log P and M 2 , we simulated RV time series for 40 000 binary systems assuming a primary mass of 20 M , which is similar to the typical masses around 10-30 M reported for WNs in the Galaxy ). The eccentricities were drawn from a flat distribution between 0.0 and 0.9. We assumed random orientations of the orbital plane in three-dimensional space and random times of periastron passage.
The RV time series at each grid point were simulated according to the observational sampling (number of epochs, time base) for each object. This allowed us to compute the probability that a given orbital configuration in log P and M 2 would give us an RV signal that would pass our detection threshold (i.e. ∆RV ≥ C = 50 km s −1 ). An example of such a probability detection plot is shown in Fig. 5. After we determined the probability detection plots for each object, we computed the average over the grid of log P and M 2 . Assuming a flat mass-ratio between 0.1 and 2.0 with a primary mass of M 1 = 20 M , we averaged over the secondary mass to calculate the average detection probability as a function of orbital period (Fig. 6). In Paper 1, we adopted C = 10 km s −1 , so that the efficiency of detecting WC binaries dropped below ∼80% at orbital periods of 1000 d (Fig. 8 in Paper 1). Here, due to a larger adopted threshold C for our WNE sam-Article number, page 5 of 17 A&A proofs: manuscript no. aanda ple, this efficiency drops below 80% already at periods of ∼100 d (Fig.5). Figure 6 shows the global detection probability of our RV campaign as a function of the orbital period. As a first approach, the observed binaries in our sample can be corrected with their respective detection probability, and thus the number of true binaries in our sample can be estimated. We ignored WR 3 in this exercise as we did not classify it as a RV variable system (Fig. 4), and we ignored WR 6 since we are unable to confirm the orbital period with our data (App. A). After defining a quantity p(P) to be the detection probability for each binary detected with a given orbital period P, we expect our sample to contain 1/p(P) binary systems. The largest correction factor comes from WR 138, which has a detection probability of ∼0.33. We therefore expect three similarly long period binaries (1/0.33) to be present in our sample on average. Given the six detected binaries in our sample, we expect about nine binaries in the sample, corresponding to an f WNE int of 0.56. This straightforward counting correction is only possible because the orbital periods of the binary systems are known, but it can be strongly affected by small sample statistics.

Adopting a Bayesian approach
Because p(P) is a strong function of P, the bias-corrected binary fraction depends on the underlying period distribution. Here, we adopt a more formal approach to determine both within a Bayesian framework. The model parameters of a population with binaries include those that describe the period distribution, binary fraction, inclination, eccentricity, mass ratio, orientation of the orbital plane, time of periastron passage, and so on. Ideally, we would like to calculate the likelihood of observing the RV time series for each object given these model parameters. However, this is computationally too intensive, and thus we adopted a different strategy. We divided the 16 WNE stars into four ∆RV bins: ∆RV ≤ 50 km s −1 (nine objects), 50 ≤ ∆RV ≤ 250 km s −1 (four objects), 250 < ∆RV ≤ 650 km s −1 (three objects) and ∆RV ≥ 650 km s −1 (no objects). The shape of the period distribution is parametrised by its power-law index π as follows: p(log P) ∼ (log P) π . (2) The period distribution further has lower and upper bounds, called log P min and log P max . Along with f WNE int , we explored π, log P min , and log P max as model parameters. Because very short period WNE binaries are present in our sample, we fixed log P min = 0.15 [d]. We explored log P max in the range of 3.0 to 5.0 in steps of 0.1. We explored values of π from −1.0 to 1.0 in steps of 0.1, and f WNE int from 0.3 to 1.0 in steps of 0.02. For each bin in this three-dimensional parameter space, we simulated 40 000 sets of 16 WNE stars, that is, over 6 × 10 8 populations. The assumptions for the eccentricity, inclination, and so on are the same as the MC simulations above (Sect. 5.1). Based on our prior knowledge of the known orbital periods (Table 2), when possible, we enforced that our simulations had at least seven binaries in period bins: P > 1 d (seven systems), P > 10 d (four systems), P > 100 d (two systems), and P > 1000 d (one system).
We then proceeded to compute the three-dimensional likelihood of observing the different numbers of objects in ∆RV bins as a function of log P max , π and f WNE int (Fig. 7). Assuming flat priors for each of the three parameters, we computed the marginalised posterior likelihood for the individual parameters. For each posterior, we report the mode and the 68% highest-density interval (HDI68). Our computations yield HDI68s of f WNE int = 0.56 +0.20 −0.15 , log P max = 4.60 0.40 −0.77 , and π = −0.30 +0.55 −0.53 . The upper bound of log P max is poorly constrained and is set by the maximum of the grid (5.00). However, an increase in log P max would simply result in an increase in f WNE int in order to maintain the fraction of binaries at shorter periods, which is required to reproduce our observations (Fig. 4). A negative value for π indicates a slight preference for an overabundance of shorter-period systems. The obtained values are consistent with the values estimated for mainsequence O stars ( f bin = 0.69 ± 0.09, π = −0.55 ± 0.22: Sana et al. 2012). Similarly, our best-fit results suggest a long-period tail in the period distribution that extends to at least 20 yr and probably more.  One of the caveats of our simulations is the assumption on the mass ratio distribution. We have assumed a flat distribution between 0.1 and 2. However, this need not be the case in nature, as the outcome of binary interaction scenarios is an important variable. Table 2 shows that known companions around WNE stars are main-sequence O stars (i.e. more massive companions). However, this result may suffer from a bias because binaries with low-mass companions will be much harder to detect. This is because the WR companion will not show much RV variation, and the brightness ratio will make it hard to detect any spectral features of the companion. We therefore decided to explore various power-law indices for the mass ratio distribution (κ). With the same setup as above, we assumed κ = +1.0 and -1.0 for the mass-ratio distribution and computed the same three-dimensional likelihood (Figs. C.1 and C.2, respectively). The conclusions from Fig. 7 did not change significantly, with κ = −1.0 yielding more short-period binaries (π = −0.6) alongside a larger f WNE int of 0.66, and κ = +1.0 yielding more massive companions with a slightly lower f WNE int of 0.54. In all cases, the constraint on the minimum extent of the period distribution remained similar with log P max > 3.8.

Revisiting the WC sample
We reanalysed our results from Paper 1 using the improved statistical framework presented in Sect. 5.2. We calculated the four-dimensional likelihood and posteriors for log P min , log P max , f WNE int and π in order to provide a fair comparison. Unlike in the WNE sample, we did not fix P min because we do not have many observational constraints. We explored log P min from 0.15 up to values of 2.95 and log P max from 3.0 to 5.0, both in steps of 0.1. In order to fit the absence of short P systems, a low number of intermediate P and a large number of long P binaries, we explored values of π from −1.0 to +4.0 in steps of 0.1. We also explored values of f WC int from 0.30 to 1.00 in steps of 0.02. As in Sects. 5.1 and 5.2, we assumed a flat distribution for the eccentricity between 0.0 and 0.9, random orientations of the orbital plane in three-dimensional space, and random times of periastron passage.
Article number, page 7 of 17 A&A proofs: manuscript no. aanda The WC sample was separated into four ∆RV bins: 10 ≤ ∆RV ≤ 30 km s −1 (six objects), 30 ≤ ∆RV ≤ 250 km s −1 (no objects), 250 < ∆RV ≤ 300 km s −1 (one object) and ∆RV ≥ 300 km s −1 (no objects). We then determined the likelihood of observing the same binned distribution for given values of the four model parameters. We also enforced that our simulations have binaries in the following period bins: P > 20 d (three systems), P > 2000 d (two systems). Assuming a primary mass of 10 M , we simulated 10 000 populations of 12 WC stars for each bin in this four-dimensional parameter space, hence ∼ 1.1 × 10 10 populations. Fig. 8 shows the four-dimensional likelihood plots along with the one-dimensional posteriors assuming flat priors. The white areas of the plot have a computed likelihood of zero, meaning that none of the 10 000 populations computed for this set of parameters matched the observations. The HDI68 values are f WC int = 0.96 +0.04 −0.22 , log P min = 0.75 +0.26 −0.60 , log P max = 4.00 +0.42 −0.34 , and π = 1.90 +1.26 −1. 25 . The values of f WC int and π are in agreement with the results of Paper 1, where we found the majority of WC binaries with P > 100 d and f WC int > 0.72. Based on the above analyses, the visualisation of the intrinsic period distributions of the WC and WNE populations is presented in Fig. 9. The period distribution for the WNE population using the best-fit values of log P max and π, with P min fixed at 1 d, is shown in red. The solid blue line shows the period distribution for the WC population using the best-fit values of log P min , log P max , and π. In an attempt to express the uncertainties of these parameters on the period distribution, we sampled their posteriors 10 000 times and over-plot the resulting constructed period distributions. The dark and light shaded areas show 68% and 95% of the highest density covered by the distributions, respectively. For the sake of comparison, the distribution for main-sequence O stars from Sana et al. (2012) is also shown (dashed black line).
The derived multiplicity properties might be affected by spectral variability in WR stars on the RV measurements and by the relatively small sample of stars we used. While wind variability in WR stars has a significant effect on the measured RVs, we thoroughly quantified this using our high-cadence studies (Sect. 3.2 and Paper 1) and statistically accounted for it, using a variability threshold that is high enough to avoid false positives. Similarly, the Bayesian framework presented here accounts for the small sample size. Unless our sample is significantly biased and does not (in a statistical sense) reflect the properties of the parent WNE and WC populations, our conclusions that the WNE and WC stars have different period distributions should hold. We further elaborate on these two aspects in Sects. 6.1 and 6.2.

Comparison of the binary properties with the literature
In order to assess whether our results might be biased by the small sample, we explored the literature for statistics on the WNE binary fraction. The GCWR has a total of 387 WN stars, of which 97 are classified as WNE 2 . Most of the WR stars in the past decades have been detected through infrared surveys as they are embedded in dust and/or suffer from substantial reddening or extinction, and their multiplicity properties are poorly studied. In order to discuss f WNE obs in the context of a comparable sample that has been well studied for multiplicity, we considered the stars in VDH01. The catalogue reports a total of 127 WN stars, of which 49 are of the WNE subtype. Of this subset of 49 stars, 13 are WNE binaries with derived spectroscopic orbital solutions. The observed spectroscopic multiplicity fraction for WNE stars in VDH01 is then 0.27. Table 2 shows that most of the WNE binary systems in our sample with known periods have periods of a few days to a hundred days. WR 138 is the only long-period binary system with a period of about four years. In VDH01, all of the other WNE binaries with derived spectroscopic orbital solutions have orbital periods shorter than 15 d. Therefore, the observed period distribution of the WNEs in our sample seems to be confirmed with the threefold larger literature-based sample. In contrast to the WC stars studied in Paper 1 and a similar literature search (VDH01), there seems to be a systematic skew in the observed period distribution, with a higher frequency of short-period WNE binaries than long-period ones, even in logarithmic space. There is thus no indication that our sample is biased, nor that the results of Sect. 5 are contradicted by known WR properties from the existing literature at large.

Observational biases due to the magnitude limit
Considering a population of Galactic WR stars within a certain volume, it is reasonable to assume that the WR stars in binaries will be brighter than the single WR stars (Vanbeveren & Conti 1980). This implies a potential bias due to our selection criterion, which will be larger or smaller depending on the contrast in brightness between the binaries and the single stars. Because most known WR binaries comprise OB-type companions, we made the simplifying assumption that WR binaries are twice as bright as single WR stars on average (e.g. Shenar et al. 2019). This implies that to avoid this bias, we should include single WR stars in the V-band magnitude range 12.0 to 12.7 mag.
In order to investigate this, we explored the GCWR for WNE stars within this magnitude range. For stars with missing V-band magnitudes, we searched for -band magnitudes between 13.0 and 13.7, similar to our approach in Sect. 2. For the WNE sample, we found two entries, both presumably single stars, WR 126 (WC5/WN) and WR 129 (WN4o). When we conservatively adopt these two targets as single WNE stars, f WNE int changes from 0.56 (9 out of 16) to 0.50 (9 out of 18). This difference is smaller by a factor of three than our errors and is hence negligible.
For the WC sample, we found four new entries in the GCWR. Of these, WR 114 (WC5 + OB?: van der Hucht 2001) and WR 132 (WC6 + ?: Bisiacchi et al. 1983) are candidate binaries, WR 125 (WC7ed + O9III: Williams 2019; Arora et al. 2021) is a confirmed binary, and WR 150 (WC5) is a presumably single star. As was the case with the WNEs, accounting for these objects results in f WC int that is well within the derived uncertainties. Therefore, the effect of this bias is not significant given our derived values. This can be understood through multiple reasons. First, there is a void of non-extinct massive star-forming regions at distances between about 3 and 8 kpc. While we may be able to detect WRs in the Carina star-forming region (∼2.7 kpc), we are unable to detect them in Westerlund 1 (∼3.8 kpc) because they suffer from large interstellar extinction (Clark et al. 2005). Second, a notable fraction of WC binaries are dust-producing systems (Williams 1995) and hence could experience increased circumstellar reddening. This would cause them to be fainter in the V band, which could balance out the increase in brightness due to additional companions. K. Dsilva , T. Shenar, H. Sana and P. Marchant: A spectroscopic multiplicity survey of Galactic Wolf-Rayet stars Fig. 8. Four-dimensional likelihood over log P min , log P max , f WNE int , and π for the WC sample. Assuming flat priors, the one-dimensional posteriors are also shown. For each posterior, the solid red lines show the HDI68, and the dashed red line shows the mode.

Implications for binary evolution
According to our current understanding of stellar evolution (Conti 1976;Meynet & Maeder 2003;Crowther 2007), depending on their initial mass, main-sequence O stars evolve into either cool supergiants or luminous blue variables (LBVs) before becoming WR stars by losing their envelope. Because the spectra of WR stars are thought to reflect the products of fusion in the stellar interiors, the Conti scenario (Conti 1976;Crowther 2007) proposes a spectroscopic evolution from O stars to the WNL, WNE, WC, and, if they are massive enough, WO phases before these stars end their lives as compact objects. The focus of this section is to use our new multiplicity constraints on the WC (Paper 1) and WNE stars to investigate the evolutionary connections between these various categories of objects.
Following the Conti scenario, WNE stars are thought to have evolved from main-sequence O stars via the WNL phase by wind stripping. As discussed before, the multiplicity parameters of the populations of WNE and O stars are congruous within errors (Sect. 5.2 and Fig. 9). If the population of binary WNE stars is the product of mass stripping in OB star binaries, similarities between their distributions are to be expected. For example, if a 40 M star becomes a 20 M WNE after transferring its envelope to a companion, the ratio of the pre-and post-interaction period (P f /P i ) can be determined analytically given assumptions on mass-transfer efficiency and its companion (Soberman et al. 1997). For the case of a 20 M accretor and conservative mass transfer, P f /P i = 1, while for fully non-conservative mass transfer (assuming material is ejected with the specific angular momentum of the accretor), we have P f /P i 0.9. As an example of a system with an initial mass ratio closer to unity, if the accretor has Article number, page 9 of 17 A&A proofs: manuscript no. aanda 35 M at the onset of mass transfer, then P f /P i 2.1 for conservative mass transfer, while P f /P i 2.7 for the non-conservative case. Therefore, mass transfer will not significantly modify the period distribution of O stars as they evolve into WNEs, which is compatible with the results of Fig. 9.
If WC stars are direct descendants of WNE stars, the orbital evolution from the WNE to the WC phase is expected to be mainly governed by mass loss due to their stellar winds. This results in the widening of the orbit, so that a shift to longer orbital periods is expected. To quantify this shift, we considered a 20 M WNE star in a binary system with a 30 M main-sequence O star. Even if the WNE star loses 15 M by the time it becomes a WC star, the orbital period of the system would only change by a factor of two. Therefore, only for WNE binaries in the long-period tail (P > 500 d) of the inferred distribution can mass loss lead to orbital periods compatible with the peak of the WC distribution (P > 1000 d). If indeed WC stars originate from WNE stars losing mass, then explaining the distributions shown in Fig. 9 requires that preferentially longer-period WNE binaries evolve into WC stars, while the shorter-period ones avoid that outcome.
A possibility for the longer-period regime is that we simply do not detect WNE binaries with periods greater than a few thousand days in this current RV campaign. This is shown in Fig.  6, where our detection probability drops to 40% at P ∼ 1000 d and 0% at P > 10000 d. The presence of an undetected longperiod population of WNE binaries would indeed provide natural progenitors to the long-period WC binary population. The value of f WNE int would then be much higher than reported here, however. Almost all of the apparently single WNE stars would be in longperiod binaries, resulting in a f WNE int close to 1.00.

Conclusions
We have established the observed and intrinsic multiplicity properties of a complete magnitude-limited sample of 16 northern Galactic WNE stars using data from the HERMES spectrograph since 2017. We measured the RVs using cross-correlation with a statistical framework that allowed us to derive meaningful uncertainties. Adopting a peak-to-peak variability threshold C of 50 km s −1 , which is more than three times the observed short-term variability for WR 138, we derived an observed spectroscopic binary fraction, f WNE obs , of 0.44 ± 0.12. Using MC simulations, we determined constraints on a parametrised model of the distribution of orbital parameters. In particular, we considered a distribution described by a power-law index π and upper and minimum values log P min and log P max for the period distribution, as well as the intrinsic binary fraction f WNE int . Assuming flat priors for these model parameters, we found f WNE int to be 0.56 +0.20 −0.15 . The power-law index for the period distribution, π, was found to be −0.30 +0.55 −0.53 . Both these values are consistent with what was derived for main-sequence O stars by Sana et al. (2012). Furthermore, the period distribution also favours the majority of systems at orbital periods shorter than 100 d. This is in stark contrast to what we found for the Galactic WC population ( f WC int = 0.96 +0.04 −0.22 , π = 1.90 +1.26 −1.25 .), where the majority of binaries had orbital periods of a few 1000 days. The 68% highest-density interval values of π for the WC and WNE populations do not overlap.
Taking the discrepancy in the period distributions of the WNE and WC populations at face value in our analyses and in the literature (VDH01) questions the presence of a systematic evolutionary connection between WNE and WC stars. A population of WNE binaries at P ∼ 1000 d appears to be missing. These would be ideal progenitors of the observed WC binaries given the orbital evolution due to mass loss. Secondly, the evolved counterparts of the observed short-period WNE binaries seem to be missing, despite the high sensitivity of the RV campaign in Paper 1. These results, combined with a similar deficit from literature studies, indicate the absence of a corresponding population of short-period WC binaries.
Considering the populations of main-sequence O stars, WNE and WC stars are linked from an evolutionary standpoint, it is possible for main-sequence O binaries to evolve into WNE binaries through Case A or early Case B mass transfer. However, depending on the initial period, the system would not survive if the initial mass ratio diverged drastically from unity. It is also possible for short-period WNE binaries to form via a common-envelope scenario, which could occur due to either late Case B or Case C mass transfer in a wide main-sequence O binary, regardless of the initial mass ratio. The long-period tail of the WNE period distribution is consistent with the current multiplicity statistics for LBVs (Mahy et al. 2022), possibly indicating an evolutionary sequence from long-period main-sequence O binaries → LBV → WNE → WC.
The orbital evolution of WNE binaries is mainly governed by mass loss, leading to an expansion of the orbital period by a factor of ∼1.5-2. This change is insufficient to explain the observed differences between the WNE and WC period distributions, which peak at P < 10 d and P∼5000 d, respectively. It thus appears that short-period WN binaries have a low chance of becoming WC binaries, for reasons that are not yet understood.
The existing multiplicity properties gathered from the literature on a larger sample and the smaller but statistically bettercharacterised sample presented here seem to indicate that the underlying period distributions for WC and WNE populations are different. This can offer substantial diagnostics on the evolution of these systems. However, owing to the significant evolutionary implications from the above discrepancy, it is critical to increase the sample size of such magnitude-limited studies.
A&A proofs: manuscript no. aanda Appendix A: Comments on specific objects WR 1: According to the GCWR, WR 1 is classified as an SB1?. Marchenko et al. (1998b) found a variety of photometric periods with Hipparcos, but settled on the 'best' value of 11.68 ± 0.14 d. As a follow-up, Morel et al. (1999) studied the variability of the different line profiles and measured centroid velocities of He ii λ4686. They concluded that a non-degenerate companion causing this line-profile variability on short timescales would have to be luminous enough to be seen in the spectrum (either in spectral lines or through line dilution, which is not the case). They were unable to rule out the possibility of a compact object as a companion. Even though the X-ray emission from WR 1 is variable (Wessolowski & Niedzielski 1996), the amount of X-rays is normal for its spectral type.
Extensive studies of the line-profile variability were undertaken by Flores et al. (2007) and Chené & St-Louis (2010), who found periods of 7.684 d and 16.9 +0.6 −0.3 d, respectively. They both concluded that the periodicity was intrinsic to the stellar wind and not due to orbital motion. Chené & St-Louis (2010) proposed corotating interaction regions (CIRs) as the best scenario to explain their photometric and spectroscopic results. This was further supported by a spectropolarimetric study by St-Louis (2013), who found continuum polarisation at the level of 0.5% for WR 1, indicating a large-scale structure in the wind. In this work, we find WR 1 to have a ∆RV of 37.5 km s −1 over ∼1265 d, which is significantly below the threshold of 50 km s −1 , and hence we classify it as a presumably single star.
WR 2: WR 2 is the only WN2 star in the Galaxy. It is classified as a visual binary in the GCWR and van der Hucht (2001). The optical spectrum clearly shows absorption lines of a B star in it. With the help of optical spectroscopy, Chené et al. (2019) did not find any significant RV variations over a period of ∼10 yr. Using direct imaging with assumptions on the stellar parameters, the authors showed that the contribution of the B star to the optical spectrum is smaller than expected by a factor of ten (5% instead of 50% for its distance and absolute magnitude) and is hence likely to be a background object. We measure a ∆RV of 36.9 km s −1 for WR 2 over ∼1085 d and classify it as a single star.
WR 3: According to van der Hucht (2001), WR 3 is classified as a WN3 + O4 system.  found a period of 46.85 ± 0.02 d with a mean semi-amplitude of the RV curve, K, of 33 km s −1 , using both He II 4686 (emission) and Hγ (absorption) lines. However, the RV curves of the absorption and the emission line were found to be roughly in phase (relative phase shift of 0.15 ± 0.03), implying that the absorption lines are intrinsic to the WR star. Given that the line profiles of WR 3 are fairly stable and triangular, Marchenko et al. (2004) hypothesised it to be a single WN3 star with a high hydrogen abundance and were able to model the optical spectrum. The GCWR classifies it as a WN3ha (hydrogen rich with absorption lines). In this work, we find WR 3 to have a ∆ RV of 13.1 km s −1 over ∼ 120 d measured using the N V line at 4945 Å. We therefore reject the period found by  and classify it as a single star.
WR 6: The GCWR classifies WR 6 as a SB1? system. The system demonstrates spectacular photometric (Marchenko et al. 1998b, and references therein) and spectral (Morel et al. 1997;Flores et al. 2007) variability with a period of 3.77 d. Additionally, this periodicity has been shown to vary on longer timescales Robert et al. 1992). The line-profile variability was thought to be caused by wind-wind collisions with a (compact) binary companion, but Morel et al. (1997) argued that this is rather due to a structured wind.
With uninterrupted photometry over 136 d from the BRITE satellite, Schmutz & Koenigsberger (2019) modelled the light curve and showed that the epoch-to-epoch variability was due to rapid apsidal motion in a binary. Koenigsberger & Schmutz (2020) used archival UV and X-ray data collected over several decades to model the emission peaks of various line profiles and also computed RVs. They were able to explain both the RVs and light curves with rapid apsidal motion in an eccentric binary (WR+B, e=0.1). Moreover, the authors also presented arguments for a third companion. In this work, we measure a ∆RV of 170.0 km s −1 over ∼1025 d and classify WR 6 to be a binary system, in agreement with other findings.
WR 7: The GCWR classifies WR 7 as a single star. With Hipparcos photometry, Marchenko et al. (1998b) observed a long-term trend but did not find any periods. We find that WR 7 exhibits strong line-profile variability over short timescales, although we were unable to discern any conclusive periodicity. With observations over ∼640 d, we measure a ∆RV of 46.5 km s −1 and classify it as a single star. However, recent observations from TESS indicate the presence of two short periods, about 0.3 d and 2.6 d. A detailed analysis of these periods is in progress (Toalá et al. in prep).
WR 10: In the GCWR, WR 10 is classified as a single WN5h star. It is a visual binary with an A2V star as a companion, but it lacks RV variation (Niemela et al. 1999), which we confirm. In this work, we find WR 10 to have a ∆RV of 10.2 km s −1 over ∼680 d, thus classifying it as a single star.
WR 110: WR 110 is classified as a single WN5b (broadlined) star in the GCWR. Over a time span of ∼680 d, we do not find significant RV variations with a ∆RV of 31.8 km s −1 and classify it as a single star.
WR 127: According to the GCWR, WR 127 is a binary system with spectral type WN3b + O9.5V. A spectroscopic orbit was derived by Massey (1981) with a period of 9.5550 ± 0.0002 d. Lamontagne et al. (1996) and Marchenko et al. (1998b) found shallow eclipses in the light curve over the same period. de La Chevrotière et al. (2011) improved the orbital solution, modelled the wind-wind collision zone, and revised the classification of the O-star, reclassifying WR 127 as a WN5o+O8.5V system. We find WR 127 to have a ∆RV of 332.0 km s −1 over ∼2650 d and classify it as a binary system. WR 128: According to the GCWR, WR 128 is classified as a WN4(h) star with a binary status of SB2?. It was found to be a photometric variable by Antokhin & Cherepashchuk (1985) with a period of 3.871 d in the V band. The depth of the eclipses led the authors to conclude that the companion was a neutron star with a bright accretion disk.  could not verify this period in their study, and the nature of the companion has also not been confirmed yet. In this study, we find WR 128 to have a ∆RV of 9.7 km s −1 over ∼1120 d and classify it as a single star.
WR 133: WR 133 is classified as a WN5o + O9I system. In addition to being a SB2 system, it also shows photometric (Marchenko et al. 1998b) and polarimetric  variability, although the photometric variability is not coherent with the binary orbit. The spectroscopic orbit was first derived by Underhill & Hill (1994). Richardson et al. (2021) combined interferometric and spectroscopic data to derive the first visual orbit for a WN star and also to improve upon the parameters of the system. In this study, we measure RVs for WR 133 with N v λ4945. We find a ∆RV of 111.6 km s −1 over a timescale of ∼2900 d and classify it as a binary system. WR 138: According to the GCWR, WR 138 is classified as a WN5o + B? system. The companion was classified as a O9 star by Annuk (1990), who found a period of 1538 d. This was further confirmed and improved upon by Palate et al. (2013), who derived a period of 1521 ± 35 d with optical spectroscopy. They also found evidence of wind-wind collision with X-ray observations. Finally, Richardson et al. (2016) were able to resolve the system with interferometry and, with spectroscopy, they classified the companion as an O9V star with a fairly high sin i of ∼350 km s −1 . We measured the RVs with weak N v lines, particularly a combination of N v λλ4604, 4620, and N v λ4945. Over a timescale of ∼2900 d, we find WR 138 to have a ∆RV of 92.1 km s −1 and classify it as a binary system.
WR 139: Also known as V444 Cyg, WR 139 is an eclipsing short-period binary system with a spectral classification of WN5o + O6III-V with a period of 4.212435 d. It is a well-studied object both with photometry Marchenko et al. 1998b, and references within) and spectroscopy (Marchenko et al. 1994). Marchenko et al. (1994) also found evidence of a stable, short-period signal (P∼0.36 d) from the He ii λ4686 and He ii λ5412 radial velocities that they attributed to pulsations. The RVs were measured with N v λ4945. With a ∆RV of 618.3 km s −1 over ∼2950 d, we classify WR 139 as a binary system.
WR 141: According to the GCWR, WR 141 is a SB2 binary system with a spectral type of WN5o + O5V-III and a period of 21.7 d. Orbital solutions were derived by both Marchenko et al. (1998a) and Ivanov et al. (1999), but we list the former in Table 2 as it is an SB2 solution and as both derived orbital parameters are consistent within errors. In this study, we find WR 141 to have a ∆RV of 209.8 km s −1 over ∼2040 d and classify it as a binary system. The RVs were measured with N v λ4945.
WR 151: WR 151 is another short-period eclipsing binary with a period of 2.12691 d that is classified as a WN4o + O5V system in the GCWR. Hutton et al. (2009) combined RV measurements from Lewis et al. (1993) and photometry to fit a simultaneous orbital solution. We measure RVs with N v λ4945, find a ∆RV of 650 km s −1 over ∼950 d and classify it as a binary.
WR 152: The GCWR classifies WR 152 as a WN3(h). We measure RVs using N v λ4945 and find a ∆RV of 8.3 km s −1 over ∼920 d, thus classifying it as a single star.
WR 157: According to the GCWR, WR 157 is a visual binary with a spectral type WN5o (+B1II). As part of a larger sample, Maíz Apellániz et al. (2021) observed the system with spectroscopy and reclassified the companion as a B0.7 II star. In this study, we measure RVs for the WR star with N v λλ4604, 4620 and note that it barely moves. We measure a ∆RV of 11.8 km s −1 over ∼2170 d. However, He i λ4471 shows an SB2 behaviour, implying that the system is a triple (Fig. A.1). He ii λ4541 also appears to show the same SB2 behaviour, although it is less certain due to the blend with the WR emission. We conclude that the system has a tight OB + OB binary and a WR star whose connection to the OB binary is yet to be established. This is similar to BAT99 126, the WR star in the Large Magellanic Cloud that was found to be a quadruple system (Janssens et al. 2021).  simulated 10 000 sets of 12 WC stars for each bin in this fourdimensional parameter space, hence ∼ 7 9 populations.     , and π for κ = −1.0. For each posterior, the solid red lines show HDI68, and the dashed red line shows the mode.