S 308 and other X-ray emitting bubbles around Wolf-Rayet stars

S 308 is an X-ray emitting bubble that surrounds the Wolf-Rayet star WR6. The structure shines in the optical as well and is thus known as the Dolphin Nebula. Due to its large angular extent, it has been covered at only 90% with past XMM-Newton observations. Thanks to the unique dataset provided by the all-sky survey performed in X-rays by SRG/eROSITA, we can show for the first time the image of the bubble in its entire extent in this band, together with its spectral characterization. Moreover, we have tried to apply the same procedure for other wind-blown bubbles detected in the optical/IR and we searched for X-ray extended emission around them. We first analyzed the diffuse emission of S308, providing a detailed spectral analysis. We then considered a sample of 22 optical/IR selected wind-blown bubbles from a previous study based on WISE data, providing an estimate of the X-ray flux for the first time. We obtained the best fit for S308 with a two-temperature non-equilibrium plasma model (kT$_{1}=0.8_{-0.3}^{+0.8}$ keV and kT$_{2}=2_{-1}^{+3}$ keV) showing super-solar N abundance and low absorption. We did not detect any of the 22 optical/IR emitting bubbles in X-rays, but using our best fit model, we estimated the 3$\sigma$ flux upper limits for each bubble. We demonstrate the new possibility offered by SRG/eROSITA to study known wind-blown bubbles and look for other ones. A two-temperature plasma description seems to fit the data quite well for S308. Since all of the 22 bubbles studied still remain undetected by SRG/eROSITA, it is very likely that absorption effects and spatial compactness are responsible for the challenges standing in the way of detecting these bubbles in soft X-rays.


Introduction
Wolf-Rayet (WR) stars are massive stars, responsible for carrying considerable amounts of material in the interstellar medium (ISM) via strong winds.Among the different phases of a massive star's life, the WR phase is the very late stage when the outer layers have already been expelled and what is shining is just a very hot core.There are different types of WR stars, depending on which elements dominate their optical spectra.The main types are WN, WC, and WO, which have their spectra dominated by nitrogen, carbon, and oxygen, respectively.Some WN stars show hydrogen in their spectra, while WO and WC do not show hydrogen.For a modern overview, see for example Sander et al. (2019).As a consequence of this very powerful activity and expulsion of material from their envelope, WR stars play a key role in determining the composition of the ISM.Some authors even consider the possibility that the Solar System originated in the wind-blown bubble of a WR star (Dwarkadas et al. 2017).As a consequence, WR stars also play a key role in the star formation process.
Since WR stars are known to be in the very late stage of their lives, they are considered to be the closest example of stars that are going to explode as core-collapse supernovae (Woosley et al. 2002;Sander & Vink 2020) in a relatively short amount of time.
⋆ e-mail: fcam@mpe.mpg.deThis makes studies of these objects particularly interesting, not just for stellar physics, but also relating to our understanding of supernova explosion mechanisms and supernova remnant studies.At the time of writing (summer 2023), the Wolf-Rayet Star Catalogue 1 (Rosslowe & Crowther 2015) lists 667 WR stars.
WR6 is a prototypical WR star which is visible as a naked eye object in the Canis Major constellation (also known as EZ Canis Maior, HD50896 van der Hucht (2001)).Given it is very luminous also in the X-ray band, it has been covered by deep grating observations with Chandra and XMM-Newton.Making use of this rich dataset, Oskinova et al. (2012) and Huenemoerder et al. (2015) showed the presence of an intriguing plasma emitting structure.While most of the spectrum is well described by collisional equilibrium thermal plasma components, the very high quality of the spectra shows another component in the hard X-rays, whose origin is still unknown.It is well known that WR stars are surrounded by dense material, so it is possible that this additional component originates there.One noticeable example of hard X-ray emitting plasma around massive stars is found in colliding-wind binaries (CWB).However, works such as Pradhan et al. (2021) and Abaroa et al. (2023) demonstrate the description of the dense wind around WR stars is still far from being complete.On top of this unclear theoretical scenario, Koenigsberger & Schmutz (2020) provided evidence for an un-seen companion around WR6.This clearly could alter the interpretation given by the earlier studies of Oskinova et al. (2012) and Huenemoerder et al. (2015), with the additional cold plasma component possibly due to reprocessed material around or by the unseen companion.In the context of having a possible double system in WR6, we notice how del Palacio et al. (2023) highlighted the presence of a hard X-ray emitting component around the so called Apep system, which is composed of two WR stars.
Moreover, the WR6 system is particularly interesting also because it is surrounded by an extended circumstellar bubble called S308.Detected for the first time with ROSAT (Wrigge 1999), this structure was extensively studied in X-rays, employing several XMM-Newton observations (Chu et al. 2003;Toalá et al. 2012).Given the large extent of the remnant, this dataset was barely sufficient to cover 90% of its extent in X-rays.As a consequence, also due to its very low surface brightness, this circumstellar bubble has never been imaged in its entire extent in the X-rays.In addition, the impossibility of finding a background that is fully independent from the source in the XMM-Newton dataset led Toalá et al. (2012) to conclude that their spectral analysis was likely to be biased for this reason.
Therefore, the launch of the X-ray telescope eROSITA in 2019 (Predehl et al. 2021) on board the Spectrum-Roentgen-Gamma Observatory (SRG, Sunyaev et al. (2021)) has provided a unique opportunity, thanks to its large field of view and XMM-Newton-like spectral resolution.These factors provide the ideal conditions for studying extended sources, such as galaxy clusters, supernova remnants, and bubbles.Given the unique dataset provided by the first four all-sky survey scans of eROSITA (with the merged dataset called eRASS:4), in this paper, we release the image and spectral analysis of the bubble in its entire extent, for the first time.Moreover, the eRASS:4 dataset allows us to look for X-ray emission around 22 similar bubbles that are known from their optical/IR emission (Toalá et al. 2015).In Section 2, we describe the data reduction process, while in Section 3, we present the image and the spectral analysis carried out on S308 together with the flux estimate on the others bubbles.Finally, in Section 4, we discuss our results and in Section 5, we present our conclusions.

Data reduction
We started the data analysis with the reduction of the eROSITA dataset that passed over S 308 for four consecutive all-sky scans (April 11-19, 2020, October 15-22, 2020, April 10-21, 2021, and October, 12-27, 2021) summing up to a total observing time of 1860s.To carry out this operation, we employed the eROSITA science analysis software (eSASS) version 2112142 (Brunner et al. 2022).With the command evtool, we merged the event files from tiles 102114 and 105114, using the events from all the four scans (dataset called 'eRASS:4').With the same command, we cleaned the merged event file from solar flaring, selecting all photon patterns available (PATTERN=15) in order to maximize the signal.To produce the image shown in Figure 1, we corrected for the vignetting over the whole field of view in the 0.2-10 keV band, obtaining a final vignetted-corrected exposure time of 586s.
For the spectral analysis, we analyzed the entire nebula, masking the point sources.After this masking operation, we finally extracted the spectrum, background, redistribution matrix file (RMF), and ancillary response file (ARF) with the command srctool.

Spatial Analysis of S308
To analyze the diffuse emission from the X-ray emitting bubble, we started by creating the image shown in Figure 1.The image was smoothed with CIAO X-ray Data Analysis software (Fruscione et al. 2006) tool csmooth, using as the parameters a minimum significance of S/N=2 and maximum of 7, along with a minimum smoothing scale of 5 pixels and maximum smoothing scale of 20 pixels.As can be seen from Figure 1, the inner part of the bubble seems to be considerable less luminous than the outer parts.The straightforward justification for this would be that the intense winds coming from the bright central Wolf-Rayet star have swept out the material expelled from the star itself in a previous phase.We also created an optical/X-ray composite image as shown in Figure 2.

Spectral analysis of S308
To perform the spectral fitting analysis, we employed the Python interface of XSPEC (Arnaud 1996), called PyXSPEC, using the Cash statistics (Cash 1979) version implemented in XSPEC (see Kaastra 2017, for a discussion).We will refer to this quantity later as CSTAT.We fitted the spectra from 0.3 keV to 10 keV.The errors on the parameters are expressed in 1 σ confidence level (68%).We extracted the spectra from the entire nebula and from the background using the regions shown in Figure 1.
After removing the point sources, we started to characterize the background following the modeling described in Okon et al. (2021).We employed three different spectral components: one power law to describe the cosmic X-ray background (CXB) with the slope fixed at 1.4 plus three thermal components (VAPEC, Smith et al. (2001)) representing the galactic halo and the Local Hot Bubble (LHB).The galactic halo is modeled by one component with kT 1 =0.658 keV and a second component with kT 2 =1.22 keV.The third component representing the LHB is modeled with a thermal model of temperature kT=0.105keV.The normalization of all these four models are left free to vary.Before fitting this model to the background spectrum, we modeled the instrumental background using the best-fit model obtained from the filter wheel closed data, as described in Yeung et al. (2023).In this way, we obtained the best-fit model for the background, accounting for both sky and instrumental contributions.
We then fit the best-fit background model and the source model together to the source spectrum.As presented in Camilloni et al. (2023), we took into account the geometric area of the extraction region (expressed in sr) in our spectral modeling ("constant" factor in front of the source model), which was kept fixed.While fitting simultaneously source and background model, we kept all the parameters frozen, except an overall global normalization in the background model.To better estimate the errors in the spectral fitting analysis, we used a code based on the Markov chain Monte Carlo (MCMC) library emcee (Foreman-Mackey et al. 2013).We initialized our walkers using a Gaussian distribution centered on the parameters of a first fit run.We then explored the posterior using logarithmically uniform priors.We left the code run for 40000 steps, selecting only the last 2000 steps to ensure most of the walkers were converged.All the abundances are expressed in solar units.
We started modeling the diffuse emission with a single absorbed (TBabs, Wilms et al. (2000)) VAPEC model, obtaining a CSTAT/DOF (where DOF=degree of freedom) close to keV.This second component displays a clearly supersolar N abundance (N/N ⊙ =9.1 +0.6 −1.1 ).The best-fit parameters are shown in Table 1.
The result is quite different if two non-equilibrium collisional shock models are employed (VPSHOCK, Borkowski et al. (2001)) and the same abundance ratios are adopted.Again, we considered the galactic absorption of these components using the TBabs model.Despite the fact that the CSTAT statistics was improved, whilst obtaining the best-fit value so far, the temperatures seem a bit unrealistic, being well above 1 keV for both components (Table 2).
We then tested again a double non-equilibrium model with the abundances of N, O, Ne, Mg, and Fe free to vary in one component.We have chosen to let vary these specific ions since these are the lines that are more prominent in our spectra.In the second model component, we left only C and N free to vary.Our aim here was testing the assumption made by Chu et al. (2003) and Toalá et al. (2012), using two non-equilibrium models instead.In the first instance, we obtained a fit with the two components inverted, due to the overly high degeneracy among the parameters.To solve this issue, we started from the best fit obtained with the fixed abundance and we launched another fit run, allowing the abundances listed above to be free to vary.
We retrieved abundances similar to those originally adopted by Chu et al. (2003).The final picture is a two-phase gas: one presenting lower temperature and very low ionization timescale, while the other presents a high temperature and a slightly higher ionization timescale.The plasma is far from reaching equilibrium in both components, as per the definition of ionization timescale given in Arnaud (1996); Borkowski et al. (1994Borkowski et al. ( , 2001)).From Figure 3, we observe that the hot and high ionization timescale component is much weaker.In performing the same exercise of letting free the abundances in the double VAPEC model, we obtained CSTAT/DOF=1.02.Nevertheless, with the double VAPEC we recover the temperatures of Toalá et al. (2012), but obtaining a much higher N abundance (N/N ⊙ = 9.3 +0.5 −0.9 ).Given the difference of CSTAT/DOF is around 1% between the two models, we suggest that deeper data with much higher spectral resolution would be needed to finally determine whether the difference between the two models is real or the effect is mainly related to the lack of spectral resolution.We also ran again the fits starting from 0.2 keV instead of 0.3 keV, but the results are substantially the same.We opted to continue using the results from the spectra starting from 0.3 keV, given the calibration of eROSITA is more stable and known above 0.3 keV.
In this context, given the limited spectral resolution of present-day X-ray instruments (especially on extended sources), Kavanagh (2020) discusses how both APEC and nonequilibrium models would not serve as the most suitable models to describe the complex environments around bubbles and superbubbles.Since, with both models, we have a CSTAT/DOF value that is very close to 1 and given the very small difference of statistics between the two models tested, we conclude that the argument of Kavanagh (2020) should be valid also around windblown bubbles, as in the case of S 308.

Galactic search for other WR X-ray shining bubbles
Given the result obtained for WR6, we looked also at the WR nebulae listed in Table 3 of Toalá et al. (2015) and we calculated the flux per each of them obtained with the eRASS:4 dataset.We proceeded as follows.
We started by downloading the WISE image from the IRSA3 archive for each one of the WR nebulae.Next, for each WISE image, we defined a circular region in DS9 matching the boundaries of the nebula.The region size depends on the dimension and the shape of the bubble, which can be irregular and not always clearly visible.In this case, we tried to match the brightest boundary of the bubble visible in IR, keeping the circle centered on the WR star.Then, we masked the point sources inside the extraction region and we extracted a spectrum for each one of these regions.We used a region of the same size to extract the background from a nearby region.In this way, we extracted an independent background spectrum for each bubble.Assuming the double VPSHOCK model with fixed abundance employed on S308 (described in Section 3.2), we measured the X-ray fluxes at the position of the optical/IR bubbles listed in Table 3 of Toalá et al. (2015).The results are shown in Table 4 and the extraction regions in Figures 6 and 7. Energy (keV) -5 -2.5 0 2.5 5 data model error with the abundances left free to vary, as described in the text.For all the plots, the source model is shown as two red lines (two components).The other dash-dotted lines represent the background, demonstrating that the flat powerlaw is the most significant component of the particle background, while the two thermal dashdot components represent the galactic halo.Finally, the powerlaw shows the unresolved Cosmic X-ray Background.The orange line is the combined model.The models displayed are extracted from the median values of the last 2000 steps of the chain run of emceee.The spectra have also been rebinned for a clearer presentation.
For simplicity and given the minimal differences among Table 2 and Table 3, we proceeded with the double VPSHOCK model with fixed abundance to extract the flux of the 22 bubbles.In addition, thanks to explicit background modeling, we took into account the background in the flux estimation.Therefore, the fluxes reported in Table 4 at 3σ confidence level are only from the source: this allows us to put upper limits on the source flux for each bubble.The values at 3σ confidence level are obtained taking 99.8% of the posterior probability distribution of the last 2000 steps of the chains.We computed the significance checking whether S/N>3, where S is the counts in the source region while and N is total count in the background region.Both S and N values are normalized for the region area.In cases where the ratio is higher than 3, we claimed a 3 sigma detection.Otherwise, we estimated the 3σ flux upper limit.We do not detect any of the 22 bubbles above the background level.

Morphological and spectral analysis of S308
Thanks to the all-sky dataset provided by eROSITA, we present the first complete image of the wind-blown bubble S308 in Xrays, together with the first spectral analysis of the object in its entire extent.In addition, we were able to independently select the background region, contrary to what has been possible so far with only XMM-Newton pointed observations.From our analysis, the emission of the nebula can be described by a two component non-equilibrium thermal model (kT 1 = 0.8 +0.8 −0.3 keV and kT 2 = 2 +3 −1 keV), presenting very low column density N_H= 0.14 +0.05 −0.05 10 22 cm −2 and an enrichment of N in the first component.The abundance values are remarkably in accordance with those found by Chu et al. (2003).However, our temperatures do not agree with those of Toalá et al. (2012) when the double VPSHOCK model is employed, but they agree if the double VAPEC shown in Table 1 is used instead.Given the , where n e is the electron density of the plasma (cm −3 ), n H , is the hydrogen density (cm −3 ) and D (cm) is the distance to the source.Normalization is expressed as 10 −14 n e n H dV 4πD 2 , where n e is the electron density of the plasma (cm −3 ), n H , is the hydrogen density (cm −3 ) and D (cm) is the distance of the source Ne/Ne ⊙ 0.5  Normalization is expressed as 10 −14 n e n H dV 4πD 2 where n e is the electron density of the plasma (cm −3 ), n H is the hydrogen density (cm −3 ) and D (cm) is the distance of the source rather unrealistic abundance of N in the second component and the marginal improvement in the statistics given by the double VPSHOCK, we conclude this description is a valid alternative to the double VAPEC first presented by Toalá et al. (2012).We would like to stress that our results are not biased by the background region choice, as opposed to those of Toalá et al. (2012); this has allowed us to leave the column density and the nitrogen abundance free to vary.Looking at the cold, very soft X-ray peaked emission displayed in Figure 1, to explain the high temperature we found using a non-equilibrium plasma description, we notice how the ionization timescales are small, suggesting the gas is far from the equilibrium (Borkowski et al. 1994;Arnaud 1996;Borkowski et al. 2001).One possible explanation for such high temperatures and non-equilibrium in both components could be that the shell is expanding in a very low-density medium.This would explain the very low emissivity of the hot component shown in the spectra of Figure 3.However, we do agree with Kavanagh (2020) in that most of these models are probably not best suited to describe the complex plasma environment around bubbles, suggesting that new and different models Article number, page 6 of 11 F.Camilloni, W. Becker & M. Sasaki: S 308 and other X-ray shining bubbles Table 4. 3σ X-ray fluxes measured with eROSITA for the bubbles listed in Table 3 of Toalá et al. (2015).The bubbles are located between 180 and 360 degrees of galactic longitudes.The first column displays the identifier of the WR star and the second the flux from the source region.
WR identifier Flux (10 should be employed.For instance, the VAPEC results shown in Table 1 have temperature values similar to those of Toalá et al. (2012), but the very high enrichment of N seems a bit unrealistic, as described in Section 3.2.In general, models such as VP-SHOCK and VAPEC are more or less sufficient to describe most of the extended sources with the spectral and spatial resolution available with the actual X-ray instruments, but with the advent of forthcoming observatories such as Arcus (Smith et al. 2022), XRISM (XRISM Science Team 2020), and Athena (Nandra et al. 2013) there will be a need for more refined models.For instance, none of the models considered above take into account photoionization, making the description of this scientific case unlikely to be very accurate given the presence of a very intense illuminating source such as WR6.An intense radiation field provides a viable explanation for finding a plasma more out of equilibrium in the cold inner component than in the hot one, which we would expect to be farther out with respect to our interpretation.As visible in Figure 1, and as previously described in Chu et al. (2003) and Toalá et al. (2012), the structure of the nebula is characterized by a limb of brightened emission, in accordance with our suggestion of hot material expanding in a very low density circumstellar medium.This structure is probably shaped by the intense radiation field of WR6 and powerful winds which push the material further from the star.Looking at the compo- sition of the gas, from our fits, we find a clear signature of N, given a supersolar abundance value in one of the two plasma components, in accordance with Toalá et al. (2012).If the double VAPEC model is adopted, the value is highly supersolar.The presence of a considerable amount of N in the spectrum of S 308 is a clear signature of the enrichment of the outer layers of WR6 that are being progressively expelled during the Wolf-Rayet phase of the star.

Search for new X-ray emitting bubbles
Moreover, we decided to exploit the unique dataset offered by eROSITA to scan the sky in the galactic longitude range of 180 -360 degrees, with an aim to look for other bubbles similar to S308.We recall that so far X-ray emission has been claimed to be observed only from three other wind-blown bubbles beside S308.Except for NGC 6888 (Toalá et al. 2016, and references therein) around WR136, the other two wind-blown bubbles that are known are located within our search area, specifically: NGC 2359 around the Wolf-Rayet star WR 7 (Zhekov 2014) and NGC 3199 around WR 18 (Toalá et al. 2017).In this work, we did not detect any of the 22 bubbles at 3σ confidence level.However, adopting the double VPSHOCK as the best-fit model, we are able to provide the flux upper limits for each bubble.Several details that have been observed for a number of bubbles are displayed in Figures 6 and 7, while the explicit background determination in the flux estimation allows us to reject many false positives.The most intriguing cases are still WR75 and WR95 and from our images, some diffuse-like emission can be spotted around the position of the WR star.Looking at the spectrum of WR75 (Figure 4), another thermal component (source model) is definitely needed below 1 keV, where the contribution of the background is lower.Further investigations with XMM-Newton could be very useful to confirm the presence of an eventual X-ray emitting bubble around this WR star.
One of the possible reasons for the lack of detections for such bubbles in X-rays is that in some cases, they could be quite compact and not spatially resolved by eROSITA.We tried a simple experiment to test this hypothesis, comparing the radial profile of WR85 with that of one nearby known star (HD322835), which is clearly a point source.We selected this particular WR star as candidate since in the WISE image shown by Toalá et al. (2015) its bubble-shaped emission is quite evident and in X-rays, it appears to be very compact (the radial profiles are presented in Figure 5).WR85 is embedded in the diffuse emission of the supernova remnant RX J1713.7-3946,resulting in a higher background level; however, looking at the shape of the radial profile one can spot a small deviation from the profile of a point source.The superior angular resolution of Chandra would be really helpful in ultimately understanding whether the source is extended in X-rays or not.
Another explanation for non detecting most of these bubbles comes from a spectroscopic argument.Looking at the image of S308 (Figure 1) and its spectrum (Figure 3), the emission is highly peaked in very soft X-rays and our spectral fits retrieve very low values for the column density.Therefore, it is very likely that absorption effects should play an important role with increasing distance and this would explain why most of the bubbles remain undetected in X-rays.An interesting case is represented by WR22 and WR85, where the emission of the Wolf-Rayet star is quite clear; however, this is not necessary true for many of the other stars we studied, suggesting that absorption might play a significant role.We recall how we have masked the point source for the significance and flux estimation.As an additional indication, the relatively high galactic latitude of WR6 provides an explanation for the lack of absorption and the consequent bright emission observed in X-rays: the mechanism is clearly described in Meyer et al. (2021).Checking this conclusion with the optical extinction database 4 of Lallement et al. (2019), at 1.5 kpc, we find A V = 0.12 in the direction of WR6.From using the relation of Predehl & Schmitt (1995), this optical extinction value is equivalent to a column density in X-rays of ≈ 10 20 cm −2 , which is fully in accordance with our data.
On top of these physical arguments, we recall that the total exposure provided by eROSITA in scanning mode is quite short comparing to XMM-Newton pointed observation employed to discover the other X-ray detected bubbles.This provides a possible explanation for why we do not have a significant detection for the bubbles surrounding WR7 and WR18.

Conclusion
In this work, we provide the first complete spectral characterization in X-rays of the wind-blown bubble S308 surrounding the WR star WR6, together with the first image in its entire extent.We find that collision equilibrium and non-equilibrium thermal models provide substantially similar fit statistic, suggesting new observations with much higher spectral resolution are needed to assess the real nature of these bubbles.Moreover, we also searched for new wind-blown bubbles around other WR stars and we find none are significant above the background.However, we were able to provide flux upper limits.As discussed in Section 4, the main reasons for these non detections are likely to be absorption and compactness of the nebulae.Nevertheless, the theoretical picture to explain why these bubbles have lower temperature and luminosities than expected remains unclear.Specifically, local absorption, instabilities, and clumps with low X-ray luminosity have been advocated as a viable explanation to explain such observational facts (Toalá & Arthur 2011).However, a recent discussion presented by Dwarkadas (2023) shows how the debate is still open, with momentum-driven solution and asymmetries that could play a more significant role than considered before.New data with a much higher spatial and spectral resolution would be very helpful to finally clear the picture about these intriguing objects.4).In yellow, we show the WISE contours only for the most representative bubbles.A magenta cross indicates the position of the WR star.

Fig. 1 .
Fig. 1.False color image (Red: 0.2-0.7 keV, Green: 0.7-1.1 keV, Blue: 1.1-10 keV) of the X-ray emitting bubble S308.The star WR6 shines in the center of the bubble.The red line indicates the source extraction region, while the blue line indicates the background one.The smoothing algorithm is based on the work of Ebeling et al. (2006) and is applied as described in Section 3.1.The image scaling has been considerably stretched to highlight the faint soft X-ray diffuse emission of the nebula.

Fig. 2 .
Fig. 2. Composite eRASS:4 RGB and optical image of the S308 nebula.The eROSITA image was smoothed using a Gaussian filter with 1 σ standard deviation value to calculate the Gaussian kernel.The colors in the X-ray image correspond to R: 0.2-0.7 keV, G: 0.7-1.2keV and B:1.2-8.0 keV.The optical image was overplotted onto the eRASS:4 image with a 50% transparency.It was taken using an ASA 50cm reflecting telescope and Hα and O[III] filters.Optical image courtesy of Nik Szymanek & Telescope Live.

Fig. 3 .
Fig.3.Spectrum of the diffuse emission surrounding WR6.Top-left: tbabs*(vapec+vapec), Top-right: tbabs*(vpshock+vpshock) with fixed abundances, Bottom: tbabs*(vpshock+vpshock) with the abundances left free to vary, as described in the text.For all the plots, the source model is shown as two red lines (two components).The other dash-dotted lines represent the background, demonstrating that the flat powerlaw is the most significant component of the particle background, while the two thermal dashdot components represent the galactic halo.Finally, the powerlaw shows the unresolved Cosmic X-ray Background.The orange line is the combined model.The models displayed are extracted from the median values of the last 2000 steps of the chain run of emceee.The spectra have also been rebinned for a clearer presentation.

Fig. 4 .
Fig. 4. Spectrum of the potential X-ray emitting bubble surrounding WR75.The spectra have been obtained and the lines are labeled as in Figure3.

Fig. 7 .
Fig. 7. Exposure corrected images in X-rays (Red: 0.2-0.7 keV.Green: 0.7-1.1 keV.Blue: 1.1-10 keV) of the regions surrounding the WR stars listed in Table4.Top left: Corresponding image from WR40 to WR101.We excluded the images of WR68 whose emission is contaminated by surrounding unrelated diffuse emission.In green and red, the source and background extraction regions are shown, respectively, employed for the flux estimation (Table4).In yellow, we show the WISE contours only for the most representative bubbles.A magenta cross indicates the position of the WR star.

Table 1 .
Best-fit values obtained with a double VAPEC model.The normalization factors ("factor") of the instrumental and sky background are also expressed.

Table 2 .
Best-fit values obtained with a double VPSHOCK model.The normalization factors ("factor") of instrumental and sky background are also expressed.

Table 3 .
Best-fit values obtained with a double VPSHOCK model, with N, O, Ne, Mg and Fe free to vary in the first component and only C and N are free to vary in the second.The normalization factors ('factor') of instrumental and sky background are also expressed.