© The Authors 2024

1 Introduction

Wolf–Rayet (WR) stars are evolved massive stars with typical current masses of 10–25 M_⊙. The minimum initial mass at solar metallicity for a single star to become a WR star is ~25 M_⊙ (Crowther 2007, and references therein). WR stars usually fuse helium or heavier elements in their cores, having lost part or all of their outer hydrogen layer. However, the most massive among them, which can reach up to 80 M_⊙ or even higher, can still be fusing hydrogen and exhibit substantial hydrogen content in their atmospheres (for a review, see: Crowther 2007). WR stars are considered to mainly descend from O type stars, with the WR phase accounting for ~5% of their entire lifetime. They exhibit high mass-loss rates and dense winds that can reach terminal velocities of the order of ~1–2 × 10³ km s⁻¹ (Crowther 2007). Due to their high-velocity winds, they are characterised by broad emission line spectra, which, depending on their spectral type, show lines of N (WN spectral type), C (WC spectral type), or a combination of C and O (rare WO spectral type) in addition to strong lines of He. Studying the WR phase in the lifetime of massive stars is crucial for better understanding the final stages of stellar evolution. Moreover, through their high mass-loss rates and their ultimate death as core-collapse supernova, they dominate the feedback to the local interstellar medium (ISM), can trigger or pause the formation of new stars, and overall contribute to the chemical evolution of galaxies.

In WR stars, due to the high density and higher metal abundances of their winds, it is anticipated that X-rays would be either entirely or highly absorbed. This absorption is expected unless the X-rays are generated at locations farther out from the star, where the wind density is comparatively low. Indeed, theoretical simulations make evident the higher opacity of the WR winds – especially for the WC type – compared to O-type stars, with the winds of the WR stars remaining optically thick up to a few thousand stellar radii (Figure 2 in Oskinova 2015). Consequently, the presence of X-rays in WR winds necessitates mechanisms capable of generating hot gas at significantly larger radii compared to the lower-density winds observed in OB stars.

X-rays from WR stars have been observed both for presumed single stars and colliding wind binary (CWB) systems, the latter being usually X-ray brighter. For single WR stars, the X-ray emission could be produced from a distribution of shocks embedded in their stellar winds, although alternative scenarios have been proposed (Oskinova 2015; Rauw 2022). The WN type is generally brighter than the WC type, whose wind opacity to soft X-rays is greater. In CWB systems (usually WR+OB), a part of the kinetic energy of the winds injected into the wind–wind interaction region is converted into heat and the shock-heated plasma emits additional X-rays (Stevens et al. 1992).

The X-ray spectra of WR stars are characteristic of optically thin thermal plasmas with emission lines of abundant elements (Mg, Si, S, Ar, Ca, Fe, etc.) and a temperature that ranges up to a few million degrees for single stars to tens of millions of degrees for CWBs. In particular, the Fe XXV line (at ~6.7 keV) is observed for plasma temperatures of > 10M degrees (>1 keV) that cannot be generally reached in the winds of single nonmagnetic massive stars and could be treated as a diagnostic of CWBs (Rauw 2022).

In our Galaxy, 669 WR stars¹ have been identified so far (July 2023), with the vast majority located within the Galactic plane with a confirmed binary fraction of ~30–40% (e.g. van der Hucht 2001; Crowther 2007; Rosslowe & Crowther 2015; Dsilva et al. 2020, 2022, 2023). The young massive star cluster Westerlund1 (Wd1; l = 339.55°, b = −00.40°) is an ideal laboratory with which to perform a comprehensive study of a large coeval population of WR stars. The Wd1 WR population is the richest known in a single star cluster, amounting to an impressive number of 24 WR stars thought to have been created in a single burst of star formation 3–5 Myr ago (Clark & Negueruela 2002; Negueruela & Clark 2005; Groh et al. 2006; Crowther et al. 2006a; Gennaro et al. 2017; Andersen et al. 2017)². For comparison, the Arches and Quintuplet clusters have 13 and 19 known WR stars, respectively (e.g. Clark et al. 2018a,b, 2019a). Moreover, Wd1 is considered to be the most massive young stellar cluster known within the Milky Way (e.g. Clark et al. 2005; Brandner et al. 2008; Gennaro et al. 2011; Lim et al. 2013; Andersen et al. 2017), and it is the closest super star cluster to the Sun (~4 kpc Davies & Beasor 2019; Navarete et al. 2022; Negueruela et al. 2022). The Wd1 WR population comprises 16 WN type and 8 WC type stars, with binarity being firmly confirmed through optical and infrared studies for 4 out of the 24 WR stars (e.g. Crowther et al. 2006a; Bonanos 2007; Clark et al. 2011, 2020; Koumpia & Bonanos 2012; Ritchie et al. 2022). However, the true binary fraction is likely higher based on indirect signs typical for CWBs, such as infrared excess, hard X-ray colours, and flat radio spectra (Crowther et al. 2006a; Skinner et al. 2006; Clark et al. 2008, 2019b; Dougherty et al. 2010).

Twelve WR stars in Wd1 have already been detected in two pre-EWOCS Chandra observations (~60 ks; Skinner et al. 2006; Clark et al. 2008). These WR stars exhibited hard X-ray spectra and X-ray colours overall, suggesting a binary nature. While the early shallow X-ray observations provide very important insights into the WR population, only deeper X-ray observations can provide crucial information on the nature of all the known WR stars in Wd1, through detailed spectral and timing analysis. The EWOCS observations (Extended Westerlund 1 and 2 Open Clusters survey), along with the pre-EWOCS Chandra archival observations, reach a total exposure time of ~1.1 Ms, and provide the unique opportunity to obtain a very deep X-ray look into the Wd1 WR population.

This is the second paper in the series of EWOCS papers, the aim of which is to study the WR population in detail. The paper is organised as follows: In Sec. 2 we provide a very short review of the data, in Sec. 3 we present our analysis and results, including the WR source detection, variability (short- and long-term), X-ray colours, and spectral analysis. In Sec. 4 we discuss our main findings regarding the spectral analysis and the implications for binarity. In Sec. 5 we summarise our results. We assume a distance to Wd1 of 4.23 kpc (Negueruela et al. 2022).

2 Observations

This work is based on 44 Chandra observations of Wd1, 36 of which have been taken as part of the EWOCS Chandra Large Program (Proposal number: 21200267, P.I. Guarcello; ACIS-I; 967.80 ks), while 8 are pre-EWOCS archival observations (ACIS-S; 151.93 ks). All EWOCS observations are pointed towards Wd1, while two of the pre-EWOCS observations are pointed at Wd1 and six at the magnetar CXO J164710.20-455217. The EWOCS observations were performed between June 2020 and August 2021, while the pre-EWOCs observations occurred between June 2005 and February 2018. A detailed description of all the observations, exposure times, technical characteristics, and data reduction is presented in Guarcello et al. (2024). The list of Chandra observations used in this paper are contained in the Chandra Data Collection (CDC) 153³.

3 Analysis and results

3.1 X-ray detection and photometry of Wolf–Rayet stars

Due to source confusion, and the bright and irregular background in Wd1, a rather sophisticated strategy using four different methods was adopted for source detection in five energy bands. The source validation and photometry was performed using the ACIS-Extract software (Broos et al. 2010). A detailed description of the source detection and photometry of the 5963 validated sources in Wd1 can be found in Guarcello et al. (2024).

We detected 20 out of the 24 WR stars known in Wd1. Among them, 8 stars (J, X, P, I, Q, V, H, and S) were detected in X-rays for the first time. The exact locations of the X-ray detected WR stars within the star cluster are shown in the Chandra colour image of Wd1 (Fig. 1), while a close-up view of each X-ray detected WR star can be seen in Fig. 2. The extreme crowding characterising the core of the cluster is evident in both Fig. 1 and Fig. 2. Moreover, the close-up colour images of each WR star in Fig. 2 reveal that 12 out of the 20 WR stars appear to be predominately hard (F, N, W, J, R, X, G, P, Q, D, and V) while 1 WR star (H) appears to be predominately soft.

The vast majority of the X-ray sources in Wd1 are very faint (Guarcello et al. 2024), and since WR stars are expected to be significantly brighter than X-ray-emitting low-mass stars, even in the case of contamination, the low-mass sources would contribute only a few photons to the signal we attribute to the WR stars. Thus, no impact on our study is expected, except for the very faint WR stars. We therefore inspected NIRCam JWST images in order to evaluate if the X-rays are generated indeed from the WR stars or from another nearby source. The WR stars are very bright and saturated and therefore it was impossible to detect very close sources. However, another much fainter infrared source is detected within the extraction region of R and J, while within the extraction region of H we find another very bright IR source.

In Fig. 2 we also show the extraction regions of the WR stars and nearby sources with magenta colours. The sizes of the extraction regions were determined by the ACIS-Extract software, which takes into account the local point spread function (PSF) and reduces the region size where necessary in order to avoid overlap with nearby sources (Guarcello et al. 2024). For all WR stars, except for stars V and A, the extraction regions encircle 86–91% of the PSF. Inside the yellow circle that illustrates the position of star V (Fig. 2), ACIS-Extract finds two sources of similar brightness. The source to the east corresponds to the optical position of star V. Due to overlap, the extraction region for star V is reduced to encircle 60% of the integrated PSF. Similarly, for star A inside the yellow circle that illustrates its position (Fig. 2), ACIS-Extract finds two more sources, one to the east and one to the west of star A, which contain 130.7 and 233.7 net broad-band counts, respectively. Both sources are much fainter than star A, but since the source in the east is very close to the WR star, ACIS-Extract reduces the size of the star A extraction region to 58% of the integrated PSF (hereafter referred to as the small extraction region). Since star A is the brightest among the WR stars in the sample, and the nearby stars are comparatively very faint, we decided to also use a larger extraction region corresponding to ≥90% of the integrated PSF (hereafter referred to as the large extraction region) and contains 10295 net counts. The contribution of the neighbouring sources to the large extraction region is 3% for the broad-band net counts (most of it originating from energies of higher than 2 keV). In the following analysis (long- and short-term periodicity, colours, and spectra), we used both small and large extraction regions of star A. However, as no significant differences were observed, we present only the results from the large extraction region.

Basic properties of the detected WR stars are listed in Table 1, where sources are sorted according to their X-ray brightness. In the broad band (0.5–8.0 keV), four WR stars have more than 1000 net counts and are all of spectral type WN. Eight WR stars have between 100 and 1000 net counts (two WC and six WN), and another eight WR stars have between 20 and 100 net counts (two WC and six WN). The undetected sources are presented in Table 2 and are all of WC spectral type.

Fig. 1 Smoothed colour Chandra image of the central 4 × 4 arcmin of Westerlund 1. Soft (0.5–2.0 keV), medium (2.0–4.0 keV), and hard (4–8.0 keV) X-ray photons are shown with red, green, and blue colours, respectively. Yellow circles show 19 out of the 20 detected WR stars. The WR star N is not depicted, as it is located ~4.5 arcmin south of the centre of the star cluster.

Fig. 2 X-ray-detected WR population in Wd1. We present the X-ray-brightest WR star from the top left panel and moving rightwards the X-ray brightness decreases. Each panel has a size of 10″ × 10″ and the colour coding is the same as in Fig. 1. The yellow regions are shown just for illustration purposes centred at the X-ray coordinates of the WR stars, while the magenta regions correspond to the extraction regions of detected sources used for photometry and spectral analysis, except for star A where a larger extraction region was used (Sec. 3.1).

3.2 Short-term variability

In order to explore any possible short-term variability, which could indicate physical changes for the WR stars within a timescale of few hours, we extracted the probability-weighted light curves of the single observations using the CIAO tool glvary. The tool utilizes the Gregory-Loredo algorithm variability test (Gregory & Loredo 1992), through the glvary implementation within the CIAO tool, and is optimal for detecting variability in a source since it corrects for instrumental variations that are not accounted for by the dmextract tool⁴. The glvary tool splits the events of a single observation into multiple time bins, looks for significant deviations and then assigns a variability index (varindex) to each observation based on the probability for a source to be variable. It also outputs the average standard deviation, and the average count-rate of the light curve. The variability index ranges from 0 to 10, and sources with varindex≥6 are categorised as definitely variable. We used a minimum time bin set at 1000 s (mintime=1000 s), and explored the variability for the brightest WR stars (>200 net counts) in the broad (0.5–8 keV), soft (0.5–2 keV), and hard (2–8 keV) energy bands.

We detected short-term variability for four WR stars, namely star A (OBSID 22988), star N (OBSID 23287), star U (OBSIDs 22983 and 22986) and star L (OBSID 22977). This is quite a rare event (2–5% occurrence) taking into account that there are overall 41 observations for star A, 39 for stars U and L, and 38 for star N. The variability detected for stars A, U, and N is exclusively observed within the broad and/or hard X-ray bands, with no corresponding variation in the soft band. This suggests that the observed variability is more likely linked to temperature variations of the thermal plasma responsible for the X-ray emission. Temperature variations assuming a constant local absorption would also result in variations of the soft flux, albeit at much lower levels not detected by the glvary tool. On the contrary, star L shows short-term variability in the soft band alone, pointing to possible changes in the local absorption of the binary system within that observation. The short-term light curves of all observations where variability has been detected are presented in Fig. A.1 in the Appendix.

3.3 Long-term variability and periodicity search

Given their long baseline spanning more than one year, the EWOCS observations provide a unique opportunity to search for long-term variability of the WR stars. The identification of long-term periodic signals through X-ray data reveals modulations in the flux of the wind of the WR star. This modulation could be intrinsic to the wind, for example a corotating interaction region (Chené & St-Louis 2011; Massa et al. 2014) as suggested in the case of WR6 (Ignace et al. 2013), or could be due to orbital modulation revealing the presence of a companion.

In the search for variability and periodic signals we did not include the pre-EWOCS observations since they are separated by several years during which the sensitivity of the ACIS detector to soft X-rays has declined owing to contamination build up⁵, and we cannot be sure that any periodicity was stable for such a long time. For the construction of the long-term light curve we used the probability-weighted light curve with the optimal time binning for each observation provided by the glvary tool (see Sec. 3.2), as well as the raw light curve, with a single point per observation, provided by the dmextract tool. We initially performed a chi-square test for the probability-weighted light curves of all sources except for the five faintest, to access the significance of the long-term variability. We find that all sources show long-term variability at the 99% confidence level except for star U which we do not include in the following analysis. We then constructed the Lomb-Scargle (LS) periodogram which is designed to detect periodic signals in unevenly spaced observations using the astropy class LombScargle⁶ (Zechmeister & Kürster 2009; VanderPlas et al. 2012; VanderPlas & Ivezić 2015). For the construction of the LS periodogram we used the probability-weighted light curves provided by the glvary tool. The advantage of using the individual light curves instead of the single count-rate value for each observation, is that we can account for variability within the same observation. In addition, by utilising the glvary light curves we can avoid instrumental variations affecting the light curve (i.e. bad pixels, chip gaps, etc.). One aspect that requires extra attention when using the glvary light curves to construct periodograms, is the fact that the tool assigns multiple time bins even when the light curve is considered constant. This affects only the calculation of the false alarm probability (FAP) levels which result in lower power values. In order to address this issue, for observations with identical count-rate values in consecutive time bins, we retained only one value and the average time stamp.

In the four panels of Fig. 3, we display for the brightest WR star A, from top to bottom, the raw light curve produced by the dmextract tool (one point per observation), the probability-weighted light curve produced by the glvary tool (multiple points per observation), the LS periodogram, and the overlaid window function (WF) associated with the time series produced by the LombScargle astropy class. Inspecting the light curves in the first two panels, we notice that for a few observations, the raw light curve exhibits very low count-rates compared to that of the probability-weighted light curve. This happens as the source passes through the detector’s chip gaps, and the corresponding decrease of the effective area is not accounted for by the dmextract tool in the production of the raw light curve. In the third panel along with the LS periodogram, we indicate the highest significance period and its 1σ error produced by the fit, as well as, with horizontal dashed, solid, and dotted lines, the 0.01%, 0.1%, and 1% FAP levels respectively, which are calculated through the LombScargle class of astropy using the bootstrap method with 100 000 iterations. In the bottom panel, we overplot on the LS periodogram the WF (cyan colour) which is centred and normalised at the frequency of the main periodicity. All the peaks of the WF (apart from the highest) correspond to aliases produced by the uneven sampling of the underlying “true” light curve. Therefore, the WF helps us understand which of the peaks present in the LS periodogram could correspond to real periodicities. In Fig. 4, we show the phase-folded probability-weighted light curve for star A based on the highest significant period measured. For all the phase-folded light curves presented in this study, phase 0 is defined as the minimum flux in the broad band.

We performed the same analysis for all X-ray detected WR stars except for the five faintest ones (D, V, E, H, and S) since their signal-to-noise ratios are too low and no firm conclusions could be drawn. In order to be conservative, we consider significant periodicity signals those that exceed at least 0.1% of the FAP, and whose periodograms do not exhibit many peaks with similar significance. These criteria, apart from star A, are fulfilled by another eight sources (B, L, F, N, J, X, P and Q), for which we present light curves and LS periodograms in Fig. B.1 in the Appendix. The phase-folded probability-weighted light curves of the same eight WR stars are presented in the Appendix in Fig. B.2. In Figs. C.1 and C.2 in the Appendix, we retain, for future reference, the light curves, periodograms, and phase-folded probability-weighted light curves for the sources that did not meet our criteria for significant periodicity detection (O, W, R, G, and I).

For the nine sources with measured periods we also searched for possible secondary periodicity. To achieve this, we fitted a sinusoidal function to the light curve of the source. The period of the sinusoidal was set at the main period of the LS periodogram. We then computed the residual light curve which resulted after the subtraction of the best-fit sinusoidal from the initial light curve. As a final step, we constructed the LS periodogram of the residual light curve and searched for a secondary periodicity. We found secondary periodicities for three (A, B, and F) out of the nine sources with measured periods. The residual LS periodograms are presented in Fig. D.1 in the Appendix. For stars A and F, even though the significance of the measured period is not high (FAP<1%) we present them since they coincide with optical periodicities.

In the following, we report and discuss the main results for the WR stars with measured main and secondary periodicities. The light curve of star A shows significant long-term broad-band variability (first two panels of Fig. 3), with a clear periodicity at 81.75 ± 4.63 days (power=0.58; third panel of Fig. 3). We constructed the light curves and periodograms for the soft and hard bands and we find that they are similar and always agree within the errors. The residual time series after the removal of the main signal shows a shorter period of 7.70 ± 0.04 days (Fig. D.1; FAP<1%), which is compatible with the optical modulation of 7.63 days reported by Bonanos (2007). These authors were not able to detect the ~81 day period since their observations roughly covered 20 days.

For star B, which is also a confirmed binary, we find a period of 3.514±0.009 days (power=0.47), which coincides with the period calculated from optical observations (3.52 days; Koumpia & Bonanos 2012). Even though we observe other highly significant peaks in the periodogram that may not be adequately described by the WF, we consider this result robust due to its agreement with the optical data. The search of the residual time series resulted in a highly significant secondary period at ~2.2 days.

For star L, which is a confirmed binary, we find a main periodicity at 46.41 ± 1.38 days (top right panel in Fig. B.1). If we zoom in on the highest significant peak in Fig. D.2 of the periodogram, we find that another significant peak (FAP > 0.01%) at approximately 53 days corresponds to the period of 53.95 days found by Ritchie et al. (2022) in the optical band. The peaks at approximately 46 days and 53 days are aliases of each other (Fig. D.2 in the Appendix).

Star F is a confirmed binary with a period of 5.05 days and a suspected triple system (Clark et al. 2011). We find a clear peak at 2.263 ± 0.004 days while the time series of the residuals reveals a second period at 5.15 ± 0.02 days (Fig. D.1), which surpasses the FAP=1%. We report it since it is compatible with the period found in the optical band by Clark et al. (2011).

For star N, which is long suspected to be a binary due to its infrared excess and spectroscopy (Crowther et al. 2006a; Ritchie et al. 2022), we recover a clear signal at FAP=0.01% in the periodogram, with a relatively long-term periodicity of 58.61 ± 2.79 days. The inspection of the residual time series does not reveal any significant peak.

Stars J and X exhibit their highest peaks at periods of ~5.15 and ~9.01 days, respectively. No secondary period is detected for these two stars. However, as the main periods barely pass the 0.1% FAP threshold and a series of other peaks are located at relatively high power we characterise these periods as uncertain. Star P shows a peak slightly more significant than the 0.1% FAP threshold at ~5.19 days, and a clear periodicity in its phase-folded light curve (Fig. B.2). Therefore, we mark also this period measurement as uncertain. Finally, for star Q we measure a main period at 7.4 days. However, another significant peak that is not described by the WF is present at ~130 days. The search of the residual time series for a secondary period did not result in any significant peaks. Due to the presence of these two significant peaks, we mark the main period of star Q as uncertain.

We present the highest significance period measurements of the nine WR stars in Column 10 of Table 1 and we mark as uncertain (with the letter u) the period of three stars according to the previous discussion.

Fig. 3 EWOCS light curves and periodogram of star A. Top panel: broad-band raw light curve produced by the dmextract tool. Second panel: broad-band probability-weighted light curve produced by the glvary tool. Third panel: resulting LS periodogram and the FAP levels of 0.01, 0.1, and 1% are marked with the dashed, solid, and dotted horizontal lines. A clear peak at about 81 days is detected. Bottom panel: Same as third panel, with the normalised window function overlaid on the periodogram.

Fig. 4 Phase-folded probability-weighted light curve of star A in the broad band.

3.4 X-ray colours

We computed the X-ray colours of the WR stars, both from the single and the combined observations, using only the EWOCS data. Our aim is to find possible variability of the spectral properties of the WR stars, even for those stars which are too faint for a detailed spectral analysis.

We define the X-ray colours as C₁ = log₁₀(S/M), C₂ = log₁₀(M/H) and C₃ = log₁₀(S/H) where S, M, and H are the net counts in the soft (0.5–2.0 keV), medium (2.0–4.0 keV), and hard (4.0–8.0 keV) bands, respectively. We calculate the X-ray colours and their corresponding uncertainties for the 20 detected WR stars using the BEHR tool (Bayesian Estimation of Hardness Ratios; Park et al. 2006). The tool provides reliable estimates and confidence limits even when one or both bands have very low counts, since it evaluates the posterior probability distribution of the X-ray colours taking into account the number of source and background counts. Moreover, in order to estimate the spectral parameters of the WR stars (especially the faint ones) from their X-ray colours, we created grids by simulating absorbed thermal spectra.

3.4.1 Combined observations

We present the C₂ over C₃ diagram for the combined observations, along with the simulated absorbed thermal grid, in Fig. 5. Confirmed binary systems from optical studies are marked with open symbols, while WN and WC stars are marked with filled stars and circles, respectively. Errors are calculated with the BEHR tool at the 90% confidence level.

In this colour–colour diagram the WR stars are concentrated mainly in two groups. The first group, hereafter referred to as “medium colour group”, comprises seven WR stars (A, B, U, L, O, E, and S) located within the simulated grid of estimated spectral parameters at N_H=1–2.5 × 10²²atoms cm⁻² for the absorption and kT=1.8–5.0 keV for the temperature of the emitting plasma. These thermal plasma temperatures are typical for CWBs, and indeed three out of four confirmed optical binary systems belong to this group.

The second group, hereafter referred to as the “hard colour group”, is located at the bottom-left of the colour-colour diagram and contains ten WN (J, R, W, G, X, P, I, Q, V, and D) and two WC (F and N) stars. This group contains the faintest non-confirmed binary WR stars of WN spectral type and the two brightest WC spectral type WR stars, all appearing hard in the colour Chandra image (see Fig. 2). From their colours alone we cannot confirm whether their spectra are intrinsically hard, or whether they appear hard due to high absorption. The WC type stars N and F, are known to be dust-forming stars (Crowther et al. 2006a), and therefore absorption could play a significant role.

Only one WR star, the WC type star H, appears isolated at the softest side of the colour-colour diagram. Star H is among the two faintest X-ray detected WR stars and as a consequence its colours exhibit quite large errors. In fact, H is the only WR star appearing soft in the colour Chandra image (see Fig. 2). Taking the large errors into account, the source could fall on the medium colour group or remain at the soft side of the diagram.

Fig. 5 X-ray colour–colour diagram of the WR stars. WN stars are marked with star symbols and WC stars with circles, while open symbols correspond to optically/infrared-confirmed binary systems. Error bars are reported at the 90% confidence level. The grid shown with dotted grey lines represents the values of NH and kT obtained from simulated absorbed thermal spectra. The grey arrows and the grey numbering indicate the direction in which the parameters of the two spectral components increase and their corresponding values, respectively.

3.4.2 X-ray colour variations among the individual observations

In order to explore possible spectral variation experienced by the WR stars over different timescales, we calculated their X-ray colours in each single observation. We have used observations that have at least 10 net counts in the broad band in order to impose statistical confidence. In Fig. 6 we show the colour-colour diagrams for the two brightest sources (A and B), while in Fig. E.1 we show the same diagram for other four bright WR stars. The remaining WR stars are too faint in the individual observations to extract meaningful results. We always present with small symbols the individual observations while the larger symbols refer to the combined one.

For star A, the colour–colour measurements overlaid on the simulated grid are symmetrically distributed around the 2.1 keV value obtained from the combined observation. They exhibit a rather constant absorption of N_H ~2.0×10²² atoms cm⁻², and a temperature varying between 1.5 keV and 3.0 keV. Individual observations of star B show variations in both N_h (1.0–3.5 × 10²²atoms cm⁻²) and kT (1.3–5.0 keV) that cannot be attributed to colour uncertainties. Again, the symbol representing the combined observation falls amidst the individual observations. The same is true for the rest of the WR stars, although the error bars on the colours of the individual observations are quite large.

3.5 Spectral analysis

To gain further insight into the X-ray properties of the WR population, the spectra of all available observations (EWOCS and pre-EWOCS) were fitted together. Since the spectral properties of the WR stars can vary over time, particularly if they are in binaries, any fitting results obtained correspond to average measurements. The details of the photon extraction and background region selection method are reported in Guarcello et al. (2024).

To fit the spectra of the WR stars we used absorbed optically thin plasma models in thermal equilibrium, as have been typically used to describe the shock heated plasma observed in WR stars. We introduced a model of the form: tbabs×phabs× ∑ (vapec). The model component tbabs is used to account for the ISM absorption fixed to the value of 2.05 × 10²² atoms cm⁻², which is the weighted average absorption towards Westerlund 1 using the N_H tool of HEASARC (HI4PI Collaboration 2016). When we convert this value to optical extinction using the relation N_H=1.9×10²¹ A_V (Zombeck 2006), we find an optical extinction of A_V ~10.8, which is the exact same value as the average extinction towards Wd1 calculated by Negueruela et al. (2010). This value is also in agreement with the value A_V =11.40±2.40 reported by Damineli et al. (2016). The second absorption parameter, phabs, was let free to vary in order to account for local variations of the absorption of mainly the stellar winds and/or of the ISM. The model vapec corresponds to collisionally ionised diffuse gas with the possibility to vary the abundances. The plasma temperature is reported in keV, and its normalisation is given in units of $\frac{{10}^{-14}}{4\pi D^2}\int n_e{n}_{H}dV$, where n_e and n_H are the electron and hydrogen densities integrated over the volume V of the emitting region, and D is the distance to the source in cm (Smith et al. 2001). In all cases we used the abundance table provided by Wilms et al. (2000) as it is the most updated table to be used with the tbabs model for ISM absorption.

There has been evidence for non-equilibrium plasma in wide WR+O systems, while in shorter period WR binaries electron and ion temperatures were found to be equal (see review by Rauw 2022). For this reason, we have also tested the vpshock model which allows one to account for non-thermal equilibrium (Borkowski et al. 2001). The model resulted in good fits and in many cases fewer spectral components were necessary in comparison to the vapec model. However, the ionisation parameter was in most cases unconstrained, and the best-fit plasma temperature was either very high and/or only constrained by a lower limit particularly for the low-luminosity sources. In addition, the plasma temperature was almost in all cases degenerate with the abundances parameter. Despite these problems, for two cases (stars L and O) the non-equilibrium model worked slightly better by fitting the strengths of the emission lines with one thermal component instead of two. Regarding particular cases (stars A and F) where three thermal components were required in order to fit the equilibrium case (described later on in this section), the non-equilibrium case required two thermal plasma components with the second thermal component constrained only by a lower limit for star A, while for star F three plasma components were necessary with the ionisation parameter unconstrained and the softer thermal component frozen. Therefore, the need for up to three thermal components cannot necessarily be attributed to the use of a specific model but rather to the spectral complexity of these bright sources. Overall, a non-equilibrium model could work well in some cases but the degeneracies and/or low counts of most of the spectra do not allow to distinguish between the two scenarios (equilibrium and non-equilibrium). Therefore, for clarity and uniformity, we only present the thermal equilibrium emission model vapec.

Abundances in WR stars are expected to deviate strongly from solar, and therefore we also performed the fitting using typical WN and WC abundances reported in van der Hucht et al. (1986)⁷ for the local absorption component and the thermal components (model: tbabs× vphabs× ∑ vapec). In the case of a binary system, the two sets of abundances (solar and WR) are expected to represent the two extreme situations: one in which the entire X-ray emission is dominated by the wind of the WR star and the second by a possible OB stellar companion. The reality is much more complicated, and it is likely a mixture of both, and for that reason we present both cases in Tables 3 and 4.

In XSPEC abundances are reported with respect to hydrogen and solar values, and in addition the hydrogen abundance cannot be set for the local absorption. Therefore, the hydrogen abundance cannot be set to zero. For this reason, and following similar procedure reported in Nazé et al. (2021), we set the helium abundance in XSPEC to an arbitrary value of 100 and then scale all the other element abundances relative to helium. This results in XSPEC WN abundances of 5.2, 809.7, 5.8, 73.7, 85.4, 113.5, 40.6, and 46.5 for elements C, N, O, Ne, Mg, Si, S, and Fe respectively, and WC abundances of 16 336, 0.0, 3871.6, 208.3, 1058, 359.5, 120.8, and 138.7 for elements C, N, O, Ne, Mg, Si, S, and Fe respectively. For Al, Ar, and Ca, for which we generally have only weak or no spectral information, we set them at values of 40 and 100 for the WN and WC abundances, respectively. In one or two cases of the brightest stars, lines from these elements are more prominent; we address these in Sec. 4. Choosing another arbitrary value for the helium abundance affected only the absorbing column and the normalisation of the thermal components without affecting the other parameters, the fluxes, or the fit statistic.

For the six brightest stars, we grouped the spectra using a minimum number of 15 total counts per bin, and we perform χ² spectral fitting using the XSPEC v.12.13.0c software (Arnaud 1996). For the faintest WR stars (net counts<370), we used the wstat statistics of the CIAO sherpa tool (Freeman et al. 2001). Using the wstat statistics, we grouped the spectra in bins of 5 total counts and calculated the null hypothesis probability in order to estimate the quality of the fit⁸. In Tables 3 and 4, we present the best-fit spectral parameters and their corresponding 1σ errors, along with the fluxes and luminosities of the best-fit models for solar and WR abundances respectively. We find that the best-fit values of the parameters as well as the fluxes agree between the two abundances. Fluxes were calculated using the flux command of XSPEC, and the sample_flux of Sherpa, and the flux errors are reported at the 68% confidence level. Throughout this work, whenever we refer to corrected fluxes or luminosities, we are always referring to the corrected values solely for the ISM absorption. We are not presenting the corrected values for the entire column density (i.e., for the ISM and local absorption) primarily because, in WR stars the wind itself is the source of the X-ray emission (through shocks embedded in the wind). Therefore, only the energy that escapes the wind and is released into the ISM corresponds to the intrinsic luminosity of the star (for more details: Sect. 2.2 and references therein of Rauw 2022). Moreover, due to the moderate spectral resolution of the ACIS detector, we cannot constrain well the abundances. Therefore, these were generally kept frozen, except for specific cases and for elements where a significant improvement in the fit was noted when allowed to vary. These cases are reported for each WR star in the corresponding Tables. The spectra along with their bestfit models for the brightest WR star, A, are presented in Fig. 7, while for the remaining sources they are presented in Fig. F.1.

Fig. 6 Colour–colour diagrams for the two brightest WR stars in Wd1. Black large symbols correspond to the combined exposure while small symbols to the individual observations. Error bars are shown at the 90% confidence level.

3.5.1 Description of spectral fitting results

In Table 5, we present a summary of the emission lines visible in the X-ray spectra of 17 WR stars in order of X-ray brightness. The number of lines we observe is correlated with the luminosity of the WR star and all 17 sources exhibit the iron line at 6.7 keV. In more detail, regarding the individual spectra, for star A, three thermal plasma components were necessary in order to obtain a good fit (${\chi }_{\nu }^2$ =1.05). We present the spectrum along with its best-fit model in the left panel of Fig. 7. The spectrum shows characteristic lines of a >10 MK thermal plasma, with the lines of Si, S, and Fe being quite prominent. The iron line shows significant residuals at lower energies (~6.4 keV), which we investigate further in Sect. 3.5.2. Moreover, we observe a soft flux excess with respect to the best-fit model below 1 keV. We had difficulty fitting this excess even after introducing an additional very soft (<0.5 keV) thermal component. Soft excess of this kind has been seen, for example, during the low state of the wind colliding system gamma-2 Vel = WR 11 (WC+O) where a radiative recombination continuum could be present (Schild et al. 2004). This system has a similar orbital period (78.5d) to that of star A, and shows a strong iron line. However, we noticed that most of the soft excess disappears if we adopt a smaller spectral grouping (5 counts per bin) in an attempt to highlight more details and see the emission lines in the soft part of the spectrum (right panel in Fig. 7). The observed bump at ~0.7 keV roughly corresponds to the O VIII emission line, while a small remaining excess at ~0.9 keV could correspond to the Fe-L line. We conclude that part of the observed soft excess may be due to emission lines that cannot be fitted well when the spectra are binned to provide higher counts (i.e. 15 counts) per spectral bin.

WR star B exhibits a very similar spectrum to that of star A (panel 1 in Fig. F.1). Its spectrum was fitted well with two thermal components. However, some residuals remain close to the iron line, which may indicate the existence of a lower ionisation Fe line, which is discussed in more detail in Sect. 3.5.2.

Similar to star A is the combined spectrum of star U (panel 2 in Fig. F.1). The spectrum is fitted well by two thermal plasma components. The WR star U has not been confirmed as a binary by optical/infrared studies, but it exhibits a hard spectrum (kT~2 keV), and a hint of the Fe line above 6 keV. When we group the spectrum using 5 total counts per bin (panel 3 in Fig. F.1) we clearly see the line at ~6.7 keV.

Star L is a known long period binary, and apart from the emission lines seen in other WN bright stars, shows a S XVI at ~3.5 keV and a hint of the Ca XX line at ~5.4keV (see panel 4 in Fig. F.1). The best-fit model consists of two thermal components. However, residuals remain around the Ar XVII (3.1 keV), S XVI (3.5 keV) and Ca XX (5.4 keV) lines.

Star F is the brightest WC star and it is a confirmed binary. Its X-ray spectrum is very different from those of the bright WN stars. It is harder with more prominent lines of Ar and Ca, compared to those of Si and S (see panel 5 in Fig. F.1). This could be the result of high local absorption, or an actual hardening of the spectrum due to intrinsic properties of the system, since the S line is much weaker than the Si line although less affected by absorption. The best-fit model is obtained using a three-temperature plasma model, where the elements of Ar and Ca for solar composition are free to vary, while for WR composition in addition to Ar and Ca also the abundances of Fe and Si are free to vary. However, residuals remain around the emission line of Ar (~3.1 keV; panel 5 in Fig. F.1).

The spectrum of star O (panel 6 in Fig. F.1 is fitted well by two thermal plasma components, with the strong lines requiring the abundances of the elements S, and Ar to be variable. We can only obtain a lower limit on the temperature of the hotter component. The star does not show an Fe line; however, when we use a spectral grouping of 5 counts per bin, we see a hint of a line at 6.7 keV (panel 7 in Fig. F.1).

WR star N, the second brightest source in the hard colour group, is of WC spectral type. It is highly suspected to be a binary (Crowther et al. 2006a; Ritchie et al. 2022) and its X-ray spectrum shows a strong emission line of Fe (panel 8 in Fig. F.1). Stars W and J show similar spectra to that of N, (panel 9 and 10 in Fig. F.1). All three sources require two thermal components in order to fit well the soft part of the spectrum, but stars N and W also require a high local absorption.

In the spectrum of star R we see for the second time a line at ~5.3 keV (as in L). The best-fit model consists of two thermal components plus a Gaussian line at 5.3 keV to fit the Ca XX line which is very strong, while the S XVI line, which is weaker, is not fitted well by the model (panel 11 in Fig. F.1). The spectrum of star X is hard, and is fitted by a single thermal model (panel 12 in Fig. F.1).

WR star G shows a very absorbed spectrum with a strong iron line (panel 13 in Fig. F.1). All attempts to fit this spectrum with our default model have failed. Only a single thermal component, seen through a partial fraction absorber pcfabs can adequately fit its spectrum. However, the partial absorber components were kept fixed at the best-fit values of an initial attempt to fit the spectrum. Due to poor statistics, letting these parameters free in order to calculate the errors resulted in none of the spectral parameters to be constrained. Since even in this case, the null hypothesis probability is very low, these spectral fitting results, except for the best-fit absorbed flux, should be treated with caution.

For star P we obtain the best fit when using a single thermal model with a high temperature (kT~4–5 keV) and high local absorption (panel 14 in Fig. F.1). The spectrum of star I is well fitted adopting two thermal plasma model (panel 15 in Fig. F.1), with a low and high temperature component (kT~0.1 and ~5.0 keV) and high local absorption. For star Q, a single thermal model is required to provide good fitting results. However, the temperature is not constrained, giving only a lower limit of ~3 keV (panel 16 in Fig. F.1).

The spectrum of the WR star D is fitted by a single thermal model that yielded only a lower temperature limit of ~1 keV (panel 17 in Fig. F.1). The spectrum of star V appears highly absorbed, flat, and with only the Fe line visible (panel 18). One thermal component leaves significant residuals at soft energies, and therefore, we introduced a second thermal soft component. The soft component after a initial fit converged to the value of 0.1 keV. However, when we calculated the uncertainties, and due to low statistics, the overall fit was unstable, and for that reason the temperature of the soft component was frozen at the value of 0.1 keV.

We fitted the spectra of the faintest sources E, H, and, S, with a single thermal component. Except for the spectrum of star H, which resulted in a very low null hypothesis probability, for the other stars we obtain reasonable statistics. The best-fit model of star E (panel 19 in Fig. F.1), indicates low local absorption and a high temperature (kT~3.0 keV). The spectrum of star H is soft and moderately absorbed (kT~0.28 keV), with no emission above ~3 keV (panel 20 in Fig. F.1). The spectrum of star S is flatter (i.e. harder) than the other two faint sources. It requires a moderate local absorption, and a relatively high temperature (kT~ 1.3–1.7 keV; panel 21 in Fig. F.1).

Fig. 7 Spectral modelling of star A. Left: combined spectrum of star A along with the best-fit model consisting of three vapec components (for solar composition). The fit residuals are shown in terms of sigma with error bars of size 1σ. Right: same spectrum as in the left panel but at low energies (<1.5 keV), using a grouping of 5 counts per bin.

3.5.2 Fluorescent Fe Kα line

In the spectrum of stars A, B, and J, together with the Fe line at ~6.7 keV, we observe a line at ~6.4 keV. The wide separation is not compatible with a redshifted ~6.7 keV line, but is consistent with a fluorescent Fe line. In Fig. 8, we present the spectra of these stars between 5 and 7.5 keV, using spectral groupings where the lower energy Fe line is clearly visible. The continuum is fitted with a bremsstrahlung model (model in XSPEC: bremss; magenta dotted line), and the Fe lines with two Gaussian line models (model in XSPEC: gauss; red dashed and blue dash dotted lines). The spectral fitting results for the two lines are presented in Table 6.

The fitting results confirm that the Gaussian line at lower energies is consistent with a fluorescent ~6.4 keV line in all cases, with higher signal to noise for A. The line energies indicate that the line could arise from any ionisation state of Fe below Fe XVII (for an example on the expected centroid of FeKα line on the ionisation level see: Yamaguchi et al. 2014). The presence of a fluorescent line indicates that cool and dense gas coexists with the hot gas responsible for the ~6.7 keV. The line might arise from photoionised cool material located within the wind of the system, possibly close to the wind–wind collision zone, when it is illuminated by the energetic photons produced in the shock-heated plasma. This is the first clear evidence of fluorescent Fe Kα emission for WR binary systems. Fluorescent lines have already been observed for the putative single WR star WR6 (Oskinova et al. 2012), while both thermal and fluorescent iron lines have been observed for colliding wind binaries such as η Carina with similar equivalent width measurements, and also for other massive stars such as γ Cas stars (e.g. Corcoran et al. 1998, 2001; Panagiotou & Walter 2018; Rauw 2022).

Fig. 8 Candidate fluorescent Fe line in the spectra of stars A, B, and J. We show the spectra between 5 and 7.5 keV, fitting the continuum with a bremsstrahlung model (magenta dotted line) and two Gaussian line models to represent the fluorescent (red dashed line) and hot gas emission (blue dashed dotted line).

4 Discussion

4.1 Overall X-ray properties

We observe that overall the spectral fitting results do not show any significant differences, regardless of whether we use solar or WR abundances (Tables 3 and 4). Moreover, all sources (except for H) exhibit hard spectra. Our fitting results, similarly to a number of other works on WR stars (e.g. Skinner et al. 2010; Zhekov 2012; Nazé et al. 2021), show for the majority of the sources that two temperature models can provide acceptable fits requiring a cool and a hot plasma component (kT< 0.8 keV and kT>1.5 keV). However, for the brightest WN and WC sources an additional thermal component at kT<0.2 keV is necessary, while for the fainter sources just one thermal component, usually hard (kT>1.0 keV) is adequate.

The stars in the hard colour group (F, N, J, R, X, G, P, I, Q, D, and V) show systematically higher local absorption than the stars in the medium colour group (Tables 3, and 4). However, we should stress that for fainter sources obtaining physically relevant results was challenging due to the quality of the spectra. For example, degeneracies were often present among the temperature of the hotter thermal plasma component, the local absorption, and the abundance parameters leading in some cases to poorly constrained plasma temperatures (low null hypothesis probability as shown in Tables 3 and 4). For that reason these results should be treated with caution.

All optically known WN stars, but only half of the WC stars, are detected in X-rays. This is expected, due to the different chemical composition and higher opacity of WC stellar winds (e.g. Oskinova 2015). When we fit the spectra using the generic WR abundances reported in van der Hucht et al. (1986), we find that in some cases a higher S abundance of a factor of 2–3 is required (Table 4). This is in agreement with previous analyses, indicating that the S abundance might need to be higher (e.g. Skinner et al. 2010; Zhekov 2012). In addition, the spectra of stars O and F presented very strong emission lines that could only be fitted with very high abundances of specific elements for both solar and WR abundances (notes in Tables 3 and 4).

Non-thermal emission for colliding wind binaries has been so far detected only above 20 keV (e.g. η Car and Apep Leyder et al. 2010; Hamaguchi et al. 2018; del Palacio et al. 2023). We explored for possible signs of a non-thermal component in the spectrum of star A in the Chandra energy range using the EWOCS data, both because it is the brightest WR star in Wd1 (~10000 net counts) and radio observations have shown signs of non-thermal emission (Dougherty et al. 2010). However, we found no signs of non-thermal emission. This indicates that any non-thermal X-ray emission, if present, must be significantly weaker in this energy range than the thermal emission observed, and difficult to detect with the current instrumentation (De Becker 2007).

The 1 Ms deep Chandra EWOCS observations reveal for the first time that the vast majority of the WR stars in Wd1 exhibit an iron line at ~6.7 keV. In order to present the properties of the iron line for the entire WR star population, we fit all the spectra between 5.0 and 8.0 keV with a simple bremsstrahlung spectrum and a Gaussian line model (bremms+gauss). When another line is present near the 6.7 keV line, as discussed in Sect. 3.5.2, we include an additional gauss component. We calculated the flux of the 6.7 keV line by including the cflux component in the model. The line width values were always set to zero. We present the results for the flux and the energy of the line in Fig. 9. The brightest lines are observed in stars A, B, F, and N, which are the two brightest WN and WC stars respectively. The line energies for the four brightest stars are consistent with Fe XXV emission (~6.67 keV). For almost all the remaining WR stars, the line energies agree within the errors with being an Fe XXV line, except for X, whose line centroid is found at a lower (~6.4 keV) energy.

Overall, we measure integrated absorbed fluxes (0.5–8.0 keV) that span from ~3×10⁻¹⁶ erg cm⁻² s⁻¹ (for star H) to 1.7×10⁻¹³ erg cm⁻² s⁻¹ (for star A). The corresponding observed luminosities are 6×10²⁹ erg s⁻¹ and 3.6×10³² erg s⁻¹. For all WR stars, the mean values of the observed and corrected luminosity distribution are log L_X=31.19, and log L_X=31.32 respectively (with L_X in units of erg s⁻¹). The WC stars, as expected due to the higher opacity of their dense metal rich winds, are generally fainter than WN stars, with 50% detection compared to 100%, and mean absorbed luminosities of log L_X=31.25 (WN) and log L_X=30.92 (WC) and mean corrected luminosities of log L_X=31.36 (WN) and log L_X=31.16 (WC). For different WN spectral types we see no particular differences between the observed mean luminosity and the corrected mean luminosity. However, the ISM-corrected luminosity distribution of the WNE WR stars is shifted to higher values with median corrected luminosity log L_x(WNL)=31.11 and log L_x(WNE)=31.35. This difference is small but nevertheless agrees with previous results where WNE are found to be more luminous than WNL (e.g. Oskinova 2005, 2016).

Additionally, the EWOCS observations provide a unique opportunity to investigate potential variations in colour and spectra parameters along the orbits of known binaries. This analysis, combining X-ray and optical/infrared spectroscopy, will be presented in future papers of the EWOCS series. Here we briefly mention some results on star A. The colours of the individual observations show that although N_H seems rather constant, kT values vary (see Fig. 6). However, no specific correlation was found between orbital phase and colours/spectral parameters. We also performed Lomb-Scargle analysis to investigate whether the colours and median-photon energy vary over time, but still we did not obtain conclusive results. Moreover, since we detected a ~6.4 keV fluorescent line in the spectrum of star A (see Sect. 3.5.2), we investigate how the Fe lines change at four different orbital phases using a period of 81.75 days. We fitted the combined spectra corresponding to the four different phases between 5 and 8.5 keV using a bremss+gaussian+gaussian model in XSPEC and a spectral grouping of five counts per bin. We have always detected two Fe lines, and found that they are centred at 6.40–6.42 keV and at 6.61–6.74 keV for the low flux observations, while for the highest flux observations the lines are centred at slightly higher energies 6.50–6.57 keV and at 6.75–6.80 keV.

Fig. 9 Fe XXV line properties for WR population. Fe XXV line flux over the central line energy.

4.2 X-ray luminosity versus bolometric luminosity

For O type stars (both single and binaries) a linear relationship between the X-ray luminosity (L_X) and the bolometric luminosity (L_bol) has been observationally established, and it is of the form L_X ~ 10⁻⁷ L_bol (e.g. Long & White 1980; Pallavicini et al. 1981; Chlebowski et al. 1989; Sana et al. 2006; Nazé et al. 2011; Rauw et al. 2015). For WR stars a number of studies analysing data mostly for WN stars have not found a clear connection between X-ray and bolometric luminosities, with the relation exhibiting a large scatter of more than an order of magnitude around the L_X − L_bol of O type stars (e.g. Wessolowski 1996; Oskinova 2005; Skinner et al. 2010, 2012; Nazé et al. 2021; Crowther et al. 2022). In some cases, studies have found that there is a connection between the kinetic properties of the wind and X-ray luminosity for single and binary stars (e.g. Skinner et al. 2010, 2012; Zhekov 2012), but in some other cases no clear correlation between these properties is evident (e.g. Wessolowski 1996; Nazé et al. 2021). The situation is also not very clear in terms of binarity, since binary WR stars are generally found to exhibit values of L_X/L_bol ≥ 10⁻⁷ (e.g. Oskinova 2005; Crowther et al. 2022), while a recent study of 18 confirmed binary systems revealed that not all binary systems are necessarily X-ray bright, and binaries exhibiting the same L_bol have X-ray luminosities that can differ by an order of magnitude or more (Nazé et al. 2021).

In this section, we examine the L_X − L_bol for the WR stars in Wd1 at an initial level and in light of binarity content. More thorough analysis on the connection of X-ray and wind properties will follow in a subsequent publication including all massive stars in Wd1. We use the bolometric luminosities provided in Rosslowe (2016) scaled to the adopted distance for Wd1. For the X-ray luminosities we use the ISM corrected L_X in the energy band 0.5–8.0 keV provided in Table 3 for the detected stars, while for the undetected stars we have calculated upper lim-its⁹. In the top panel of Fig. 10 we show the L_X versus L_bol plot for all the WR stars in Wd1. The L_X ~10⁻⁷ L_bol is shown with a solid line while the dashed lines represent the L_X ~10⁻⁸ L_bol and L_X~10⁻⁶ L_bol relations. The nomenclature is the same as in Fig. 5, while the upper limits are shown with black triangles. We observe that the entire detected population with the exception of the very faint stars E, S, and H, lies within the L_X ~ 10⁻⁸ L_bol and L_X ~ 10⁻⁶ L_bol lines, in agreement with what has been observed in previous works (e.g. Wessolowski 1996; Skinner et al. 2010; Crowther et al. 2022). Three of the brightest WR stars (A, B, and F) and firmly confirmed binaries show values higher than 10⁻⁷, with star A being the only one clearly surpassing L_X ~ 10⁻⁷ L_bol. We also present in the bottom panel of Fig. 10 the same relation for the 4–8 keV X-ray luminosity range. This provides us with the L_X/L_bol relation free of any absorption effects (foreground or local). In this energy range there is a clearer separation between the WR stars with A, B, and F all being a bit below the 10⁻⁷ line, and the rest of the WR stars (apart from the faintest E, S, and H) nicely following the 10⁻⁸ line.

Overall, apart from stars A, B, and F, all other WR stars show suppressed values of the L_X/L_bol relation compared to what has been found for O type stars. According to Owocki et al. (2013) lower L_X/L_bol is expected for optically thick winds and for increasing mass-loss rate. Owocki et al. (2013) also discuss how orbital separation affects their results, suggesting that close, short (day to week) period binaries should exhibit lower L_X/L_bol values. Our results are consistent with these findings since for star A, which is the only WR star clearly surpassing L_X~10⁻⁷ L_bol, we measure a long period (~81 days). This finding aligns with results presented for a number of Galactic long period WR binaries (i.e. WR25: P~208 days; WR138: P~1521 days; WR133: P~112 days; Oskinova 2005), as well as in the case of Melnick 34, an X-ray bright WR system (WN5h + WN5h, P~155 days) in the 30 Dor region of the Large Magellanic Cloud (Tehrani et al. 2019).

Fig. 10 X-ray versus bolometric luminosity for the WR stars in Wd1. The top panel depicts the “classical” correlation using the 0.5–8.0 keV band in the X-rays, while the bottom panel shows the same correlation using the 4–8.0 keV X-ray band. Upper limits on X-ray-undetected WR stars are shown with black triangles. X-ray luminosities are obtained using solar abundances.

4.3 Comparison to previous X-ray studies

Using the first two Chandra observations of Wd1 (~60 ks), Skinner et al. (2006) detected 12 of the 22 known WR stars at the time. One of them, the WC star K, remains undetected through the EWOCS data. However, the detection of stars K and E was considered marginal with 9±3 and 6±3 counts in the 0.5–7.0 keV band, respectively. Moreover, star K was only detected in the first Chandra observation (18.8 ks). Skinner et al. (2006) also do not report one of the brightest sources, star U, which was then not identified as WR star. For the stars present also in our catalogue, the unabsorbed luminosities they list in their Table 1 are higher by up to an order of magnitude compared to the ones we report. This can be explained by the fact that for 10 out of the 12 detected sources Skinner et al. (2006) use PIMMs conversion assuming a spectrum softer (kT=1.0 keV) than the actual spectrum (see Tables 3 to 4). Moreover, they report the spectral properties of star A detecting only the strong lines of S and Si. Their results agree within the errors with the spectral properties we obtain. They conclude that probably all their detections are binaries, particularly the WC stars.

Clark et al. (2008) used the same Chandra observations as Skinner et al. (2006), and examined the global properties of the X-ray sources along with photometric and spectroscopic information from optical and radio datasets. They consider the marginal detection of star K by Skinner et al. (2006) to be in error, and they present the results of the then newly discovered WR star U. They present the spectral fit results for the five brightest sources (A, B, U, L, and F) where in none of those the iron line at 6.7 keV was observed. The observed fluxes they present agree well with the ones we obtain.

Our deeper EWOCS observation produces spectra with much more detail including many more lines, the Fe line at 6.4 keV and 6.7 keV for the first time (see Sec. 3.5), and a better characterisation of the softer part of the spectrum (<1 keV). For this reason in most cases more than one thermal plasma model is required to fit adequately the spectra. Moreover, we constrain better and at lower values the temperatures for stars L and F using thermal emission components in thermal equilibrium, while for star B we always obtain higher temperatures for the thermal plasma.

In addition Clark et al. (2008) present the intensity colour diagram for all the detected sources, and find that all WR stars exhibit hard colours especially in comparison to other massive stars. Only E has a soft colour (kT~0.6 keV) and it is the faintest source in their sample. Our results agree for all the sources except for source E. We find instead that star E has hard colours similar to that of the brightest sources of the ‘medium colour group’, which is also confirmed by its spectrum (kT~3.0 keV).

4.4 The nature of stars N and T

Stars N and T are the most distant WR stars from the centre of the cluster, at a projected distance of about 4.5 arcmin. Both are WC9 dusty stars, with star N also being an SB2 candidate, and are therefore, considered to be binaries with O-type companions (Crowther et al. 2006a; Ritchie et al. 2022). X-ray emission is detected only from star N, and the EWOCS data confirm its binary nature (hard spectrum, strong Fe 6.7 keV emission line, P = 58 days). Crowther et al. (2006a) discussed the location of the two stars outside of the cluster and suggested that they were ejected due to either dynamical interactions within the cluster or following a supernova kick. Moreover they commented that neither mechanism can explain the observed binarity since dynamical interactions favour single stars while the kick would require that the companion has already undergone a supernova explosion.

We plot the proper motions of the two stars relative to the cluster members using Gaia EDR3 astrometry (Fig. 11), where positions and proper motions have been plotted relative to the cluster members in Wd1. The red colour points show the cluster members from Negueruela et al. (2022), while blue indicates the WR stars T and N along with their proper motion vectors over a 0.2 Myr baseline. Neither proper motion vector points back towards the cluster, suggesting that neither WR star has been recently ejected. In the reference frame of the cluster, the proper motions of these stars are (PM_α, PM_δ) = (0.5 ± 0.1, −0.04 ± 0.08) and (−0.28 ± 0.07, −0.08 ± 0.06) mas yr⁻¹ for N and T respectively. Relative to the inherent velocity dispersion of the cluster (with the effect of uncertainties removed) of σ_α = 0.252 and σ_δ = 0.280 mas yr⁻¹, these proper motions are not very large and are both within 2σ of the central velocity of the cluster.

For these two binary systems to lie outside of the cluster core and to have proper motions that do not point back towards the cluster makes it very unlikely that they were formed at the center of the cluster and then ejected. It seems more likely that these stars formed outside the cluster core in a more dispersed population surrounding Wd1. Many star clusters are known to be surrounded by dispersed populations (Wright 2020). Detecting this dispersed population in the optical and infrared, particularly the more-numerous low-mass component, is challenging due to the typical large contamination of memberships in the outskirts of the clusters. However, in the X-rays with the EWOCS data, a dispersed population around the cluster core is observed (Guarcello et al. 2024). In fact, out of the 24 WR stars only 4 (N, T, X, and W) are located farther than 2 arcmin from the cluster core, while the number of X-ray validated sources more distant than 2 arcmin from the core is 2559 out of the 5963. If we restrict the sample to an energy between 2 and 3.5 keV (where most of the members are) there are 1678 sources more distant than 2 arcminutes out of the 3253. While 83% of the WR stars are located within 2 arcmin, this percentage decreases to approximately 49%-57% for the general population of the X-ray sources. This reveals that a significant population of dispersed, low to intermediate-mass young stars exists in the surroundings of Wd1.

Fig. 11 Spatial distribution of Wd1 cluster members (red) and WR stars N and T (blue) with proper motion vectors shown for the latter (black lines) over a 0.2 Myr baseline. The proper motions shown are those relative to that of the cluster, which was calculated from the membership of Negueruela et al. (2022).

4.5 Binarity

A large percentage of WR stars are members of binary systems with massive stellar companions (typically OB stars). The presence of a companion can be directly revealed through optical photometry, when eclipses or periodic signals are present in the light curve of the star, or from large periodic variations of radial velocity. In spectroscopic data, the companion’s existence can also be detected through specific absorption features (e.g. Bonanos 2007; Clark et al. 2011; Ritchie et al. 2022). While these methods offer direct and unambiguous evidence of a companion’s existence, they also pose various challenges in detecting complete samples. These challenges encompass factors such as the orientation of the orbit, the duration of the data collection period, and, particularly for WR stars, the hindering of the companion’s spectral information by their broad line spectra.

When optical spectroscopy and photometry are insufficient to reveal the existence of a companion, indirect methods for detecting the companion come into play. For example in X-rays, since isolated massive stars have typically soft X-ray spectra coming from self-shocks developed in the stellar wind, it is generally assumed that hard colours or spectra (≥2 keV) are indication of high temperatures developing in the wind–wind collision zone of binary systems. Moreover, the presence of the Fe line at 6.7 keV in the X-ray spectrum is considered a diagnostic for the identification of CWBs, since it is formed only at temperatures higher than about 10 million degrees. Such high temperatures cannot be developed in the winds of single WR stars unless they are highly magnetised, which is quite rare since very few cases have been discovered (i.e. de la Chevrotière et al. 2014; Hubrig et al. 2016, 2020). At radio frequencies, binary systems can exhibit significant deviations from the expected thermal emission from single winds due to the presence of additional non-thermal emission, leading to flat or even negative spectral indexes (Williams et al. 1997; Dougherty & Williams 2000; Dougherty et al. 2005; De Becker & Raucq 2013). Excess near to mid-IR emission is also considered as an indication of binarity in persistent (i.e. WR104; Tuthill et al. 1999) and episodic (i.e. WR140; Han et al. 2022) dusty WC systems. The formation of dust is favoured in WR+O binary systems since they provide the necessary ingredients (elements C and H) and the required high density in the wind–wind collision zone. In addition the wind–wind collision zone may provide protection over ionising photons that tend to destroy the dust (e.g. Williams et al. 2005; Crowther et al. 2006b; Crowther 2007). However, the absence of dust emission in WR stars does not exclude binarity, since an episodic dust production could have prevented us from observing it.

Based on the above criteria, previous studies of the WR population in Wd1 have estimated a high binary fraction (62–70%) that could even reach unity (Crowther et al. 2006a; Clark et al. 2008; Dougherty et al. 2010). Along the same lines, an updated estimate on the binary WR fraction in Wd1, based on the new Chandra results from the EWOCS observations, is presented in this work. We are currently analysing the dust production from Wd1 JWST data which will be presented in a dedicated study as part of the EWOCS series.

Below we summarise all signs of binarity for the WR stars in Wd1 with a stress on the new indications provided by the EWOCS data. Star A is unambiguously a binary system, through optical photometric, infrared spectroscopic, radio and previous X-ray analysis (Bonanos 2007; Crowther et al. 2006a; Dougherty et al. 2010; Skinner et al. 2006; Clark et al. 2008). With the present work we determine a main period of 81.75 days and a secondary of 7.7 days, with a hard X-ray spectrum and the identification, for the first time, for this source of a strong Fe line at 6.7 keV. The identification of two periodic signals could indicate either the presence of a third body or the existence of corotating interaction regions within the stellar wind (e.g. Chené & St-Louis 2011; Massa et al. 2014). Given the specific characteristics of this system, and assuming a Keplerian orbit, the period of 7.63 days would imply that, the WR wind directly hits the surface of the companion star and therefore, no wind–wind interaction could develop, which increases the likelihood of the longer period being the orbital period. However, a colliding wind scenario could still develop in the case of the shorter period, if mutual radiative effects slow down the WR wind and prevent it from hitting the stars surface. Therefore,, more detailed modelling is necessary to understand the nature of this system.

WR star B is a confirmed binary through photometric and previous X-ray data (Bonanos 2007; Koumpia & Bonanos 2012; Skinner et al. 2006; Clark et al. 2008). We detect the iron line at 6.7 keV in its spectrum, and two periodicities, the main at 3.51 days (compatible with the optical), and the secondary at ~2.2 days. Star L is a confirmed binary through X-ray and radial velocity measurements (Clark et al. 2008; Ritchie et al. 2022). The EWOCS data reveal its hard X-ray spectrum with an Fe line at 6.7 keV, with a long period of ~46 days, although a significant peak at ~53 days is also detected (Sec. 3.3). Star F, the X-ray brightest WC star, is a confirmed binary through optical, infrared, radio and X-ray data (Clark et al. 2011, 2008; Crowther et al. 2006a; Dougherty et al. 2010), while long-term IR variability suggests the presence of a third body (Clark et al. 2011). We find a very hard spectrum with a strong Fe line at 6.7 keV and two measured periods of ~2 and ~5 days, which either strengthen the third body scenario or support a corotating interaction region origin in either star of the binary.

WR star N shows signs of binarity through its hard X-ray colours, mid to NIR excess, and infrared spectroscopy (Skinner et al. 2006; Clark et al. 2008; Crowther et al. 2006a; Ritchie et al. 2022). With the EWOCS data we confirm this finding, since the X-ray spectrum of star N is hard and shows a particularly strong Fe line at ~6.7 keV. Moreover, we find that star N exhibits a long period of ~58 days.

WR stars U, O, W, R, G, and D were considered most probably binaries due to their hard X-ray colours (Clark et al. 2008). We confirm these identifications through our colour and spectral analysis which reveal their hard spectra typical of CWBs. Moreover, all of them show the iron line at 6.7 keV.

We find that star E has colours and spectra (see Sects. 3.4 and 3.5) compatible with a hard X-ray source (kT~3.0 keV). No iron line is visible but this could be due to the low count-rate of the source (~56 net counts) in the EWOCS observations. Its nature was debated since the source appeared soft (kT~0.6 keV) in Clark et al. (2008, 2019b). However, we deem that result highly uncertain because of the very low number of net counts (~8) based only on two Chandra observations (~60 ks). Due to its hard spectrum our data indicate that E is most likely a binary.

The only indication of binarity for WR star V comes from radio observations (Dougherty et al. 2010). The source is detected for the first time in the X-rays through the EWOCS data that reveal the very hard X-ray nature of the source (kT> 1.5 keV; hard colour group), with a strong iron line at 6.7 keV, making it a strong CWB candidate. The WR stars J, X, P, I, and Q are detected for the first time in X-rays through EWOCS. We consider all of them to be binaries due to their very hard spectra and iron lines (kT>3.0 keV; hard colour group).

Star S has an intriguing nature. It has been suggested that it is the remnant of a short period binary system disrupted in a supernova explosion, and now appears to be single (Clark et al. 2014). With EWOCS it is detected for the first time in X-rays, and it is the faintest WR star (23 net counts) in the sample. Its colour and spectrum reveal its hard X-ray nature although the number of counts do not allow for further characterisation.

Star H is considered a possible CWB through its infrared excess (Crowther et al. 2006a). However, in the X-rays it presents very soft (kT~0.5 keV) colours and spectra, so its nature is quite puzzling. Within the Chandra extraction region of H, two bright infrared sources are observed through NIRCam JWST (Sec. 3.1). Therefore,, if the X-ray detection truly corresponds to emission from the WR star and not its nearby source, H could be a single star that is an unusual dust emitter or it could be a binary system with a very particular projected geometry that at the epoch of observation allows us to see only the soft wind of one of the two stars.

Three out of four X-ray non-detected WR stars, C, M, and T, are all dusty WC stars (Crowther et al. 2006a). In addition, large radial velocities have been detected for star M (Ritchie et al. 2022) further confirming its binary nature. Therefore, we consider all of them to be binaries. They are probably not detected in X-rays because their very opaque and dense winds do not allow for any of the X-ray emission to escape (at least down to our detection limit). Star K is also not detected in X-rays, and among the undetected WR stars is the only non-dusty WC star. Based on the available information, we cannot reach a conclusion regarding its single or binary nature.

To conclude, in our sample, given the hard X-ray spectra (≥ 1.5 keV) observed for 19 out of the 20 detected Wr stars, in addition to the iron line at 6.7 keV for 17 of them, we infer a binary fraction of 80% only on the basis of the X-ray data. Combining our results with infrared observations (Crowther et al. 2006a), we conclude that all these WR stars, except probably H, and K, are binaries, resulting in a binary fraction in the WR star population of ~92% (100% for WN and 75% for WC). This fraction could indeed approach unity as also predicted by previous works (Crowther et al. 2006a; Clark et al. 2008, 2019b).

A WR binary fraction potentially reaching unity in Wd1, is higher than current estimates for our Galaxy (about 30–40%; van der Hucht 2001; Dsilva et al. 2020, 2022, 2023). However, when corrected for biases, the estimated binary fraction is significantly higher, nearly reaching unity for WC stars (0.96+0.04%) and being somewhat lower for WN stars (0.52+0.14%); the higher frequency of short-period WN stars suggests different formation mechanisms (Dsilva et al. 2022, 2023). The discrepancy between the overall Galactic binary fraction and that measured for Wd1 may be attributed to the unique nature of Wd1, recognised as the most massive young star cluster in the Milky Way (Clark et al. 2005; Brandner et al. 2008; Gennaro et al. 2011; Lim et al. 2013; Andersen et al. 2017). On the other hand, the high binary WR fraction in Wd1 would require that either all O-type stars, considered to be the progenitors of WR stars, are in binary systems or that this fraction is reached only for the O stars that happen to evolve into WR stars.

In any case, a very high fraction of binaries aligns with the widely accepted notion that a significant proportion of massive stars exist in binary or multiple systems, with evidence suggesting that the binary fraction increases with the primary star’s mass. O-type stars, for instance, exhibit a binary fraction exceeding 70% (see the review by Marchant & Bodensteiner 2023). Confirming a high binary fraction for massive stars carries significant implications for interpreting their stellar evolution and the nature of young star clusters. For instance, different dynamical evolution scenarios may have different impacts on stellar multiplicity, and these should be considered in any modelling study. Furthermore, since many physical properties of young massive star clusters are currently calculated mainly assuming single-star evolution, factors such as age, mass, and consequently, the shape of the high-end of the initial mass function could be miscalculated. This highlights the significance of incorporating binarity as a parameter not only in studies calculating these quantities but also in extragalactic studies utilising population synthesis models.

5 Summary

In this paper, we present the X-ray analysis of the WR population in the Galactic young and very massive star cluster Westerlund1 using 36 new Chandra observations from the EWOCS project (Extended Westerlund 1 and 2 Open Clusters Survey), in addition to 8 archival shallow observations, reaching a total depth of ~1.1 Ms. The WR population in Westerlund 1 is the richest in our Galaxy, comprising 24 WR stars, 8 carbon type, and 16 nitrogen type. Our main findings, made possible by the exceptional depth and long time baseline of the EWOCS dataset, are summarised below:

• We detect 20 out of the 24 WR stars in X-rays, 8 of which are detected for the first time in X-rays, reaching down to an observed luminosity of ~7×10²⁹ erg s⁻¹ and a maximum observed luminosity of 3.6×10³² erg s⁻¹ for an adopted distance to Wd1 of 4.23 kpc;

• The 4 undetected WR stars are all of carbon spectral type and are probably not detected due to their very dense, opaque stellar winds;

• We find evidence for short-term X-ray variability in 4 WR stars, and find clues in one case (star L) that suggest that the variability is connected to variations in local absorption;

• Most WR stars show clear signs of long-term variability, with defined periodicities measured for nine stars. For three systems, we measure two periods, while for five we measure a periodical signal for the first time through the X-ray data. We retrieve periods consistent with the optical period for three stars;

• The majority of the WR stars have hard X-ray colours. Only one star (the carbon type star H) appears to be soft, although contamination from a nearby infrared source is possible;

• Our spectral analysis reveals details seen for the first time in the spectra of WR stars in Westerlund 1, such as emission lines of Mg, S, Si, Ca, and Fe, and it confirms the hard-X-ray nature (kT>2 keV) for all but one WR star);

• The iron emission line at ~6.7 keV (Fe XXV) is detected for the first time in the spectra of 17 out of the 20 X-ray-detected WR stars in Wd1. The presence of the line is typically associated with the very hot plasma located in the wind–wind collision zone of colliding wind binaries;

• We observe the fluorescent Fe line at ~6.4 keV in the spectra of three WR stars. This is the first clear detection of this line in WR binary systems. The presence of the fluorescent line points to the existence of cold gas that could be photoionised by the hot continuum produced in the wind–wind collision zone.

Overall, the hard X-ray spectra of the WR population and the presence of the iron line in the majority of them, along with other binarity diagnostics from optical, infrared, and radio observations, allow us to estimate a binarity fraction of 92%, which could also reach unity with more detailed and continuous observations in other wavelengths. This result has important implications for advancing our understanding of the evolution of massive stars and the formation of massive star clusters.

Acknowledgements

K. A. acknowledges support of INAF-OAPa Director Grant n. 72/2022. K. A. thanks Vinay L. Kashyap, Anthony Moffat, and Sean Gunderson for the helpful recommendations/discussions. M. G. G., C. A., E. F., L. P., and S. S. acknowledge the INAF grant 1.05.12.05.03. J. F.A.C. is a researcher of CONICET and acknowledges their support. Support for this work was also provided by the National Aeronautics and Space Administration through Chandra Proposal 21200267 issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. I. N. is partially supported by the Spanish Government Ministerio de Ciencia e Innovación (MCIN) and Agencia Estatal de Investigación (MCIN/AEI/10.130 39/501 100 011 033/FEDER, UE) under grant PID2021-122397NB-C22, and also by MCIN with funding from the European Union NextGenerationEU and Generalitat Valenciana in the call Programa de Planes Complementarios de I+D+i (PRTR 2022), project HIAMAS, reference ASFAE/2022/017. M. G. G. also acknowledge partial support from the Grant INAF 2022 YODA. The scientific results reported in this article are based on observations made by the Chandra X-ray Observatory. We also made use of NASA’s Astrophysics Data System Bibliographic Services.

Appendix A Short-term light curves of the WR stars that show variability

Fig. A.1 Light curves of the single observations where short-term variability is detected. The raw light curves produced with the dmextract tool are shown with the black line, with a time binning of 5ks. The probability-weighted light curves and their 3σ errors obtained with the glvary tool are shown with the shaded area.

Appendix B Long-term light curves, periodograms, and folded light curves of the Wolf–Rayet stars with defined period measurements

Fig. B.1 Long-term light curves and periodograms of WR stars with detected periods. The FAP levels of 0.01, 0.1, and 1 percent are marked with the dashed, solid, and dotted horizontal lines.

Fig. B.2 Phase-folded probability-weighted light curves of WR stars in Wd1 with detected periods.

Appendix C Long-term light curves, periodograms, and folded light curves of the Wolf–Rayet stars with non-detected periods

Fig. C.1 Long-term light curves and periodograms of WR stars with non-detected periods. The FAP levels of 0.01, 0.1, and 1 percent are marked with the dashed, solid, and dotted horizontal lines. In some cases, in order to show the significance of the highest peak an additional FAP level of 10 percent is shown with a dash-dotted horizontal line.

Fig. C.2 Phase-folded probability-weighted light curves of WR stars in Wd1 with non-detected periods.

Appendix D Residuals of long-term light curves

Fig. D.1 LS periodograms of the residual time series of stars A, B, and F, after the removal of the main period. The FAP levels of 0.01, 0.1, and 1 percent are marked with the dashed, solid, and dotted horizontal lines.

Fig. D.2 Star L normalised window function (cyan line) overlaid on the periodogram (black line) to reveal the aliases using as main period the ~46 (left panel) and the ~53 days (right panel) respectively. The FAP levels of 0.01, 0.1 and 1 percent are marked with the dashed, solid, and dotted horizontal lines.

Appendix E X-ray colour variability of four of the brightest WR stars

Fig. E.1 Colour–colour diagrams for 4 of the brightest Wolf–Rayet stars. Large symbols correspond to the combined exposure while small symbols to the individual observations. WN type stars are marked with star symbols, WC type stars are marked with circles, while optically confirmed binaries are marked with open symbols.

Appendix F X-ray spectra and fitting of the WR stars in Wd1

Fig. F.1 Spectra of the WR stars in Wd1 along with their best-fit models. Spectra are presented from the X-ray brightest to the X-ray faintest source.

References

All Tables

All Figures

	Fig. 1 Smoothed colour Chandra image of the central 4 × 4 arcmin of Westerlund 1. Soft (0.5–2.0 keV), medium (2.0–4.0 keV), and hard (4–8.0 keV) X-ray photons are shown with red, green, and blue colours, respectively. Yellow circles show 19 out of the 20 detected WR stars. The WR star N is not depicted, as it is located ~4.5 arcmin south of the centre of the star cluster.
In the text

Fig. 2 X-ray-detected WR population in Wd1. We present the X-ray-brightest WR star from the top left panel and moving rightwards the X-ray brightness decreases. Each panel has a size of 10″ × 10″ and the colour coding is the same as in Fig. 1. The yellow regions are shown just for illustration purposes centred at the X-ray coordinates of the WR stars, while the magenta regions correspond to the extraction regions of detected sources used for photometry and spectral analysis, except for star A where a larger extraction region was used (Sec. 3.1).

In the text

Fig. 3 EWOCS light curves and periodogram of star A. Top panel: broad-band raw light curve produced by the dmextract tool. Second panel: broad-band probability-weighted light curve produced by the glvary tool. Third panel: resulting LS periodogram and the FAP levels of 0.01, 0.1, and 1% are marked with the dashed, solid, and dotted horizontal lines. A clear peak at about 81 days is detected. Bottom panel: Same as third panel, with the normalised window function overlaid on the periodogram.

In the text

	Fig. 4 Phase-folded probability-weighted light curve of star A in the broad band.
In the text

Fig. 5 X-ray colour–colour diagram of the WR stars. WN stars are marked with star symbols and WC stars with circles, while open symbols correspond to optically/infrared-confirmed binary systems. Error bars are reported at the 90% confidence level. The grid shown with dotted grey lines represents the values of NH and kT obtained from simulated absorbed thermal spectra. The grey arrows and the grey numbering indicate the direction in which the parameters of the two spectral components increase and their corresponding values, respectively.

In the text

	Fig. 6 Colour–colour diagrams for the two brightest WR stars in Wd1. Black large symbols correspond to the combined exposure while small symbols to the individual observations. Error bars are shown at the 90% confidence level.
In the text

	Fig. 7 Spectral modelling of star A. Left: combined spectrum of star A along with the best-fit model consisting of three vapec components (for solar composition). The fit residuals are shown in terms of sigma with error bars of size 1σ. Right: same spectrum as in the left panel but at low energies (<1.5 keV), using a grouping of 5 counts per bin.
In the text

	Fig. 8 Candidate fluorescent Fe line in the spectra of stars A, B, and J. We show the spectra between 5 and 7.5 keV, fitting the continuum with a bremsstrahlung model (magenta dotted line) and two Gaussian line models to represent the fluorescent (red dashed line) and hot gas emission (blue dashed dotted line).
In the text

	Fig. 9 Fe XXV line properties for WR population. Fe XXV line flux over the central line energy.
In the text

	Fig. 10 X-ray versus bolometric luminosity for the WR stars in Wd1. The top panel depicts the “classical” correlation using the 0.5–8.0 keV band in the X-rays, while the bottom panel shows the same correlation using the 4–8.0 keV X-ray band. Upper limits on X-ray-undetected WR stars are shown with black triangles. X-ray luminosities are obtained using solar abundances.
In the text

	Fig. 11 Spatial distribution of Wd1 cluster members (red) and WR stars N and T (blue) with proper motion vectors shown for the latter (black lines) over a 0.2 Myr baseline. The proper motions shown are those relative to that of the cluster, which was calculated from the membership of Negueruela et al. (2022).
In the text

	Fig. A.1 Light curves of the single observations where short-term variability is detected. The raw light curves produced with the dmextract tool are shown with the black line, with a time binning of 5ks. The probability-weighted light curves and their 3σ errors obtained with the glvary tool are shown with the shaded area.
In the text

	Fig. B.1 Long-term light curves and periodograms of WR stars with detected periods. The FAP levels of 0.01, 0.1, and 1 percent are marked with the dashed, solid, and dotted horizontal lines.
In the text

	Fig. B.2 Phase-folded probability-weighted light curves of WR stars in Wd1 with detected periods.
In the text

	Fig. C.1 Long-term light curves and periodograms of WR stars with non-detected periods. The FAP levels of 0.01, 0.1, and 1 percent are marked with the dashed, solid, and dotted horizontal lines. In some cases, in order to show the significance of the highest peak an additional FAP level of 10 percent is shown with a dash-dotted horizontal line.
In the text

	Fig. C.2 Phase-folded probability-weighted light curves of WR stars in Wd1 with non-detected periods.
In the text

	Fig. D.1 LS periodograms of the residual time series of stars A, B, and F, after the removal of the main period. The FAP levels of 0.01, 0.1, and 1 percent are marked with the dashed, solid, and dotted horizontal lines.
In the text

	Fig. D.2 Star L normalised window function (cyan line) overlaid on the periodogram (black line) to reveal the aliases using as main period the ~46 (left panel) and the ~53 days (right panel) respectively. The FAP levels of 0.01, 0.1 and 1 percent are marked with the dashed, solid, and dotted horizontal lines.
In the text

	Fig. E.1 Colour–colour diagrams for 4 of the brightest Wolf–Rayet stars. Large symbols correspond to the combined exposure while small symbols to the individual observations. WN type stars are marked with star symbols, WC type stars are marked with circles, while optically confirmed binaries are marked with open symbols.
In the text

	Fig. F.1 Spectra of the WR stars in Wd1 along with their best-fit models. Spectra are presented from the X-ray brightest to the X-ray faintest source.
In the text