XMM–Newton observations of the massive colliding wind binary and non-thermal radio emitter Cyg OB2 #8A [O6If + O5.5III(f)] Free

ABSTRACT

We report on the results of four XMM–Newton observations separated by about ten days from each other of Cyg OB2 #8A [O6If + O5.5III(f)]. This massive colliding wind binary is a very bright X-ray emitter — one of the first X-ray emitting O-stars discovered by the Einstein satellite — as well as a confirmed non-thermal radio emitter whose binarity was discovered quite recently. The X-ray spectrum between 0.5 and 10.0 keV is essentially thermal, and is best fitted with a three-component model with temperatures of about 3, 9 and 20 MK. The X-ray luminosity corrected for the interstellar absorption is rather large, i.e. about 10³⁴ erg s⁻¹. Compared to the ‘canonical’L_X/L_bol ratio of O-type stars, Cyg OB2 #8A was a factor of 19–28 overluminous in X-rays during our observations. The EPIC spectra did not reveal any evidence for the presence of a non-thermal contribution in X-rays. This is not unexpected considering that the simultaneous detections of non-thermal radiation in the radio and soft X-ray (below 10.0 keV) domains is unlikely. Our data reveal a significant decrease in the X-ray flux from apastron to periastron with an amplitude of about 20 per cent. Combining our XMM–Newton results with those from previous ROSAT-PSPC and ASCA-SIS observations, we obtain a light curve suggesting a phase-locked X-ray variability. The maximum emission level occurs around phase 0.75, and the minimum is probably seen shortly after the periastron passage. Using hydrodynamic simulations of the wind–wind collision, we find a high X-ray emission level close to phase 0.75, and a minimum at periastron as well. The high X-ray luminosity, the strong phase-locked variability and the spectral shape of the X-ray emission of Cyg OB2 #8A revealed by our investigation point undoubtedly to X-ray emission dominated by colliding winds.

1 INTRODUCTION

The Cyg OB2 (VI Cygni) association has several particularities that stimulated the interest of astronomers. It has a diameter of about 2°, corresponding to about 60 pc at a distance of 1.7 kpc (Knödlseder 2000). It harbours a huge number of early-type stars: about 100 O-type and probably more than 2000 B-type stars (Knödlseder 2000; Comerón et al. 2002). Considering its mass, density and size, Knödlseder (2000) proposed it may be the first object in the Galaxy to be reclassified as a young globular cluster. However, a complete census of the massive star content of Cyg OB2 is not easy to achieve because of the heavy extinction in this direction (Comerón et al. 2002). So far, a spectral classification has only been proposed for its brightest and bluest members (Massey & Thompson 1991).

Another particularity of Cyg OB2 is that it contains some of the brightest OB stars of our Galaxy (see e.g. Herrero, Puls & Najarro 2002), among which we find some of the brightest X-ray emitting early-type stars. Historically, the Einstein X-ray observatory discovered the first X-ray sources whose optical counterparts were known to be massive stars in Cyg OB2, i.e. Cyg OB2 #5, #8A, #9 and #12 (Harnden et al. 1979). The same field was further investigated with various X-ray observatories: ROSAT (Waldron et al. 1998), ASCA (Kitamoto & Mukai 1996; De Becker 2001) and more recently Chandra (Waldron et al. 2004). This paper is the first of a series presenting the XMM–Newton view of Cyg OB2. It will focus on its brightest X-ray emitter, i.e. Cyg OB2 #8A (BD +40° 4227).

Cyg OB2 #8A was recently discovered to be a binary system consisting of an O6If primary and an O5.5III(f) secondary (De Becker, Rauw & Manfroid 2004c; De Becker & Rauw 2006). The system is eccentric (e= 0.24 ± 0.04) with a period of 21.908 ± 0.040 d and a major axis of 142 ± 6 R_⊙. The fact that Cyg OB2 #8A is a binary system could reconcile the high bolometric luminosity reported by Herrero et al. (2002) with its spectral classification, believed at that time to be a single O5.5I(f) star. The analysis of a time-series of the He ii λ 4686 line revealed a phase-locked profile variability likely attributed to a wind–wind interaction (De Becker & Rauw 2006).

In the framework of the campaign devoted to the multiwavelength study of non-thermal radio emitters (see De Becker 2005), Cyg OB2 #8A is a particularly interesting target. The non-thermal radio emission, supposed to be synchrotron radiation (White 1985), requires (i) the presence of a magnetic field and (ii) the existence of a population of relativistic electrons. Although in the past few years the first direct measurements of surface magnetic fields have been performed for a few early-type stars, e.g. β Cep (Donati et al. 2001), θ¹ Ori C (Donati et al. 2002) and ζ Cas (Neiner et al. 2003), the estimation of the strength of the magnetic field of massive stars remains a difficult task. Indeed, the Zeeman effect in spectral lines is concealed by the large rotational broadening of massive stars. Therefore, our knowledge of magnetic fields in early-type stars is at most fragmentary. The relativistic electrons are supposed to be accelerated through the first-order Fermi mechanism described for instance by Bell (1978), and applied to the case of massive stars by Pollock (1987), Eichler & Usov (1993) and Chen & White (1994). This process requires the presence of hydrodynamic shocks. We mention that an alternative scenario was proposed by Jardine, Allen & Pollock (1996), but we will assume here that the first-order Fermi mechanism in the presence of hydrodynamic shocks (the so-called diffusive shock acceleration–DSA–mechanism) is the dominant process. For a discussion of the physical processes involved in the general scenario of the non-thermal emission from massive stars, we refer, for example, to De Becker, Rauw & Swings (2005a) and references therein. The issue to be addressed here is that of the nature of these shocks: are they intrinsic to the stellar winds (see e.g. Feldmeier, Puls & Pauldrach 1997), or are they due to the wind–wind collision in a binary system (see e.g. Stevens, Blondin & Pollock 1992).

In the case of Wolf–Rayet (WR) stars, the non-thermal radio emitters are mostly binary systems: out of 17, 14 are at least suspected binaries (for a general discussion, see De Becker 2005). In some cases the non-thermal emitting region has been resolved and is clearly associated with the colliding-wind region (see e.g. the case of WR 140 Dougherty et al. 2005). But for the O-stars, the situation is less clear even though the fraction of binaries (confirmed or suspected) among non-thermal radio emitters has recently evolved to a value closer to that of WR stars. We know indeed that among the 16 non-thermal O-type radio emitters, 11 are confirmed binaries and two more are suspected binaries (see De Becker 2005). The recent discovery of the binarity of Cyg OB2 #8A lends further support to the second scenario where the population of relativistic electrons, and consequently the non-thermal radio emission, is produced in the interaction zone between the winds of two stars (in this binary system).

In addition, one can wonder whether non-thermal radiation can be produced in the high-energy domain as a counterpart to this non-thermal emission in the radio waveband. The detection of radio synchrotron radiation is an evidence for the presence of relativistic electrons. The existence of such high-energy electrons suggests that other non-thermal emission mechanisms such as inverse Compton (IC) scattering may be at work. The relativistic electrons accelerated at the collision zone may transfer energy to the ultraviolet (UV) photons emitted by the photosphere of the stars, therefore upscattering them to the high-energy domain. As a result, these radio non-thermal emitting binary systems could be non-thermal emitters in the X-ray and soft γ-ray domains as well. In this context, several targets have been investigated in the X-ray domain with XMM–Newton: 9 Sgr (Rauw et al. 2002), HD 168112 (De Becker et al. 2004b) and HD 167971 (De Becker et al. 2005b). Up to now, no unambiguous detection of non-thermal X-ray emission has been revealed by the X-ray observations of non-thermal radio emitters.

Beside the putative non-thermal emission, the X-ray spectrum of massive binaries like Cyg OB2 #8A is expected to be dominated by thermal emission produced by the plasma heated by hydrodynamic shocks due to intrinsic instabilities or to the wind–wind collision. As a colliding wind binary, Cyg OB2 #8A might be compared to other massive binaries where the colliding winds contribute significantly to the thermal X-ray emission (see for instance WR 140, Pollock et al. 2005; and WR 25, Pollock & Corcoran 2005). In this context, the possibility to detect a non-thermal emission component as discussed above will depend strongly on the properties of the thermal emission contributions. De Becker et al. (2005b) discussed the unlikelihood of the simultaneous detection of non-thermal radio and soft X-ray emission and proposed that short period binaries (a few days) were more likely to present a non-thermal X-ray emission below 10.0 keV than wide binaries. A more detailed discussion of this issue can be found in De Becker (2005). To investigate the X-ray emission of the massive members of Cyg OB2, we obtained four pointings with the XMM–Newton X-ray observatory. This paper is devoted to the massive binary Cyg OB2 #8A. The study of the other bright X-ray emitting massive stars, along with that of other fainter sources of the field, is postponed to a forthcoming paper.

The present paper is organized as follows. Section 2 describes the observations and the data reduction procedure. The spectral analysis of EPIC and RGS data of Cyg OB2 #8A is discussed in Section 3, whilst Section 4 is devoted to a discussion of the X-ray luminosity and to its variability. The discussion of archive X-ray data is provided in Section 5. Section 6 is devoted to a general discussion. Finally, Section 7 summarizes the main results of this analysis and presents the conclusions.

2 OBSERVATIONS

We have obtained four observations of Cyg OB2 with the XMM–Newton satellite, with a separation of about ten days between each pointing (see Table 1). The aim-point was set to the position of Cyg OB2 #8A in order to obtain high-resolution RGS spectra of this system. Because of the brightness of the massive stars located in the field of view we used the EPIC medium filter to reject optical light.

2.1 EPIC data

2.1.1 Data reduction

All three EPIC instruments were operated in the full frame mode (Turner et al. 2001; Strüder et al. 2001). We used the version 6.0.0 of the XMM Science Analysis System (sas) for the data reduction. The raw EPIC data of the four pointings were processed through the emproc and epproc tasks. The event lists were screened in the standard way: we considered only events with pattern 0–12 and pattern 0–4, respectively, for EPIC-MOS and EPIC-pn.

We selected the source X-ray events from inside a 60 arcsec radius circular region centred on Cyg OB2 #8A, excluding its intersection with a circular 15 arcsec radius region centred on Cyg OB2 #8C [right ascension (RA) = 20:33:17.9 and declination (Dec.) =+41:18:29.5, Equinox 2000.0]. The background region was defined as an annulus centred on the source and covering the same area as the circular source region, excluding its intersection with a 15 arcsec circular region centred on a point source (RA = 20:33:13.9 and Dec. =+41:20:21.4, Equinox 2000.0). For EPIC-MOS2 data, we excluded the intersection of these two regions (source and background) with a rectangular box to reject a bad column that crosses the central CCD, at slightly more than 30 arcsec away from the centre of the source region. We did the same in the case of EPIC-pn data to avoid a CCD gap located at about 40 arcsec from Cyg OB2 #8A. In each case, the boxes were adjusted after a careful inspection of the relevant exposure maps. Fig. 1 shows the source and background regions used for the three EPIC instruments in the case of Observation 1. The regions for the other observations differ only by the rotation angle. We generated the response matrix files (RMFs) with the rmfgen task for EPIC-MOS data. For EPIC-pn data, because of a problem with rmfgen,¹ we used the canned response matrix for on-axis sources provided by the SOC. The ancilliary response files (ARFs) were generated with the arfgen task. We finally rebinned our spectra to get at least nine and 16 counts per energy bin, respectively, for EPIC-MOS and EPIC-pn. All our spectra were then analysed using the xspec software (see Section 3.1).

2.1.2 High-level background episodes

We extracted a high-energy light curve (pulse invariant–PI–channel numbers >1000, i.e. photon energies above ∼10 keV) from the complete event lists to investigate the behaviour of the background level during the four pointings. High background time intervals are known to occur because of solar soft proton flares (Lumb 2002). We detected high background level episodes mostly in the fourth pointing. Even though such a high background level is not expected to affect significantly the spectral analysis of sources as bright as Cyg OB2 #8A (see e.g. De Becker et al. 2004a,b), we decided to filter our data sets to reject the most affected time intervals. After inspection of the light curves from the four data sets, we selected the time intervals below a threshold of 20 cts s⁻¹ for EPIC-MOS and 75 cts s⁻¹ for EPIC-pn. As a consequence, the effective exposure times are reduced (see Table 2) compared to the values provided in Table 1. However, this allows us to obtain the cleanest possible spectra thereby increasing the reliability of our analysis. As can be seen from the spectra shown for instance by De Becker et al. (2004b), the background correction produces spectra with large error bars on the normalized flux for spectral bins strongly affected by a high background level. In the case of bright sources like Cyg OB2 #8A, the data analysis does not suffer critically from the rejection of a fraction of the exposure time.

2.1.3 Pile-up?

Considering the X-ray brightness of Cyg OB2 #8A, one can wonder whether the EPIC data are affected by pile-up. According to the XMM–Newton User's Handbook, the count rate threshold above which pile-up may occur for point sources in full frame mode are about 0.7 and 8.0 cts s⁻¹, respectively, for EPIC-MOS and EPIC-pn in full frame mode. As will be shown later (see Table 7), the critical value is reached for some EPIC-MOS data sets, and we have to check whether our data are affected.

First, we generated pattern histograms and searched for the presence of patterns 26–29 events expected to be due to pile-up. We did not find such patterns for any of our data sets. Next, we used the epatplot task to draw curves of the singlet and doublet events as a function of PI. We obtained a first series of curves on the basis of event lists filtered using the standard screening criteria and the spatial filter described here above for the source region. A second series of curves was then built on the basis of event lists obtained with a slightly modified spatial filter, where the core of the point spread function (PSF) was excluded. Since pile-up is expected to occur mainly in the core of the PSF, these latter event lists should essentially be unaffected. As the curves built using epatplot are supposed to be pile-up sensitive, we may expect some differences between the two sets of curves if our data are indeed affected. However, no significant differences were found. Consequently, we consider that our data are unaffected by pile-up.

2.2 RGS data

2.2.1 Data reduction

The two RGS instruments were operated in spectroscopy mode during the four observations (den Herder et al. 2001). The raw data were processed with the sas version 6.0.0 through the rgsproc task. The first and second-order spectra of the source were extracted using the rgsspectrum task. We selected the background events from a region spatially offset from the source region. The response matrices were constructed through the rgsrmfgen task for RGS1 and RGS2 data of the four pointings.

2.2.2 High background level episodes

We followed the same procedure as described in Section 2.1.2 to select good time intervals (GTIs) unaffected by soft proton flares. However, as the mean level of the light curves was different according to the data set and also to the instrument, we refrained from adopting the same count rate threshold for all the data. After rejection of the time intervals contaminated by the high background, we obtained the effective exposure times quoted in Table 2. As for the EPIC data, the pointing whose exposure time is the most severely reduced is the fourth one.

3 ANALYSIS OF CYG OB2 #8A DATA

3.1 Spectral analysis

As briefly discussed in Section 1, several physical mechanisms are expected to be responsible for the X-ray emission of massive stars. On the one hand, the heating of the plasma of the stellar winds by hydrodynamic shocks is responsible for a thermal emission. These shocks may occur in stellar winds of individual stars (see e.g. Feldmeier et al. 1997) or in the wind–wind collision zone of binary systems (see e.g. Stevens et al. 1992). These two types of hydrodynamic shocks are able to produce plasma with characteristic temperatures of the order of a few 10⁶ K and of a few 10⁷ K, respectively. To a first approximation, such a thermal emission can be modelled by optically thin thermal plasma models (mekal model: Mewe, Gronenschild & van den Oord 1985; Kaastra 1992). On the other hand, non-thermal emission processes like IC scattering are expected to produce a power-law component in the X-ray spectrum. In this section, we will use composite models made of mekal and power-law models. We note that solar abundances (Anders & Grevesse 1989) are assumed for the plasma throughout this paper.

3.2 Interstellar material and wind absorption

Absorption models are required to account for the fact that both local circumstellar (wind) and interstellar material (ISM) are likely to absorb a significant fraction of the X-rays. The ISM absorption column was fixed to a value of N_H= 0.94 × 10²² cm⁻² obtained from the dust-to-gas ratio given by Bohlin, Savage & Drake (1978), using the colour excess [E(B−V) = 1.6] provided by Torres-Dodgen, Tapia & Carroll (1991).

To account for the fact that the wind material is ionized, an ionized wind absorption model was used for the local absorption component. We adopted the same opacity table as in the case of the multiple system HD 167971 (De Becker et al. 2005b), obtained with the wind absorption model described by Nazé et al. (2004). Briefly, this model considers the absorption by the 10 most abundant elements (H, He, C, N, O, Ne, Mg, Si, S and Fe) fixed at their solar abundances. The stellar fluxes were taken from the Kurucz library of spectra. For details on the model, we refer to Nazé et al. (2004). In De Becker et al. (2005b), we showed that there were no significant differences between opacities derived from various sets of parameters covering at least spectral types from O5 to O8. From the optical data already presented in De Becker et al. (2004c), we can estimate some crucial stellar and wind parameters of the stars in Cyg OB2 #8A. Following the spectral types of the two components, i.e. O6If and O5.5III(f), we adopted typical stellar radii and effective temperatures from Martins, Schaerer & Hillier (2005), allowing us to estimate the bolometric luminosity of the two stars. The mass loss rates and terminal wind velocities were then obtained from the mass loss recipes of Vink, de Koter & Lamers (2000, 2001). These parameters are quoted in Table 3. Provided that the stellar parameters of Cyg OB2 #8A lie within the parameter space discussed by De Becker et al. (2005b), we estimated that the wind absorption model that we used for HD 167971 suits the local absorption of Cyg OB2 #8A as well.

3.3 EPIC spectra

In order to fit the EPIC spectra, we tried different models including mekal and power-law components. The quality of the fits was estimated using the χ² minimization technique and the best-fitting parameter values are quoted in Table 4. We checked the consistency of our results with both the Cash statistic (Cash 1979) and the χ² statistic using a Churazov weighting (Churazov et al. 1996). We did not find any significant differences in the results obtained with the three methods. This was not unexpected as we are dealing with good quality spectra containing rather large numbers of counts per energy bin.

The best fits between 0.5 and 10.0 keV were obtained using a three-temperature thermal model. In the case of EPIC-MOS2, our fits pointed to normalization parameter values that deviated significantly from those of EPIC-MOS1 and EPIC-pn. This might be due to the bad column crossing the source region in the case EPIC-MOS2 (see Fig. 1), resulting in serious problems in obtaining an accurate ARF. For this reason, we discarded the EPIC-MOS2 data from our discussion and we will concentrate on EPIC-MOS1 and EPIC-pn data. As quoted in Table 4, the reduced χ² lies between 1.01 and 1.56 depending on the instrument and on the data set. The average plasma temperature of the three thermal emission components are, respectively, about 3 × 10⁶, 9 × 10⁶ and 20 × 10⁶ K, with average relative emission measures of about 1.0, 0.17 and 0.12. We note the good agreement achieved for the four observations, with a slightly lower temperature for the hard component of the second observation. The upper panel of Fig. 2 shows the EPIC-pn spectrum of Cyg OB2 #8A between 0.5 and 10.0 keV fitted by the three-temperature thermal model. Clearly, the most spectacular feature of this spectrum is the Fe K blend at about 6.7 keV, which is observed for all instruments and in all data sets. The presence of such a strong Fe K blend suggests that the hard X-ray emission component in the Cyg OB2 #8A spectrum has a thermal origin. Indeed, the results obtained with models including a power law were rather poor. Even though in some cases the replacement of the hardest thermal component discussed previously by a power law led to a slightly improved reduced χ², these models were rejected because they failed to fit the iron line at about 6.7 keV. The bottom panel of Fig. 2 shows the result of the fit of the EPIC-pn spectrum of Observation 1 with such a model. We see that the iron blend is poorly fitted, and that the fit of the softer part of the spectrum is less satisfactory than in the case of the upper panel of the same figure. In this case, the two thermal components yield kT of about 0.26 and 1.20 keV, whilst the power law has a photon index of about 3. We note also that we tried to use more sophisticated models with wind absorption columns affected to each emission component. However, this did not improve the quality of the fits and in most cases we obtained similar values (within the 1σ error bars) for every local absorption component. For these reasons, we used only one local absorption column as described in Table 4.

3.4 RGS spectra

As a first step, we combined the first- and second-order spectra in order to inspect the main spectral features and to identify the spectral lines. Above about 17 Å, the spectrum is very absorbed and we concentrated our analysis on the spectral domain below this wavelength. We identified the prominent lines through a comparison with the aped (Smith & Brickhouse 2000) and spex (Kaastra, Mewe & Raassen 2004) line lists (see Table 5). In order to perform a more detailed analysis of the RGS data, we obtained fluxed (RGS1 plus RGS2) spectra for our four observations (see Fig. 3). Individual lines show intensity variations from spectrum to spectrum of a few times 10 per cent and apparent velocity variations of a few hundred km s⁻¹, both of which are at the limit of detectability with the data available, whose exposures were typically only 20 ks. The reality of such apparent changes could be assessed with longer exposures of 60–80 ks. We insist on the fact that the strong absorption undergone by the X-ray spectrum of Cyg OB2 #8A in the RGS bandpass is a critical issue. Our attempts to derive reliable specific information such as line fluxes or velocity shifts did not provide any satisfactory result.

The four RGS spectra were similar in form with no obvious changes in the long-wavelength absorption cut-off near 17 Å. Models with only interstellar absorption at the expected value of 0.94 × 10²² cm⁻² are able to account for the cut-off, although we are not able to exclude a further circumstellar component for some plasma emission models. EPIC and RGS yield an effective combination with the RGS resolving lines down to Fe xvii (15.015 Å) and EPIC up to Fe xxv (18.500 Å), emphasizing the broad range of ionization conditions that exist in the X-ray emitting plasma of Cyg OB2 #8A. A similar situation was found in the X-ray spectrum of the colliding-wind system WR140 (Pollock et al. 2005). We have constructed general models involving a bremsstrahlung continuum absorbed by the expected fixed amount of interstellar material underlying line emission unconstrained by any physical plasma models from H- and He-like ions of Ne, Mg, Si, S, Ar and Ca as well as ions of Fe from Fe xvii to Fe xxv. The best-fitting continuum temperature was 1.9 ± 0.1 keV.² In the absence of high-resolution data of good statistical weight, the same line velocity profile was used for all the lines, allowing them all to be red or blue shifted from the laboratory wavelength by the same velocity and broadened by the same velocity width. Such models are able to provide a good fit simultaneously to EPIC-MOS and RGS data (we did not consider EPIC-pn data here because of the slightly poorer spectral resolution of this latter instrument).

In addition, we used the same kind of composite models as used in Section 3.3. We subtracted the background of individual spectra and then applied the response matrix for a global fitting between 5 and 35 Å. We note that the difference in the RGS1 and RGS2 count rates reported in Table 7 stems for the dead CCD of RGS1 that falls in the 10–14 Å wavelength domain. We obtained the best fit with a two-component mekal model. The χ² are slightly better if we use three thermal emission components but the error bars on the resulting fit parameters increase substantially. As we are dealing here with spectra containing sometimes small numbers of counts per energy bin, we used the Cash statistic (Cash 1979) to compare the results obtained with the χ² statistic. We did not find any significant differences between the results obtained with the two approaches. We obtained typical temperatures of about 2 × 10⁶ and 8–12 × 10⁶ K. These temperatures are close to the values obtained for the two softer thermal components of the 3-T model fitted to the EPIC data (see Section 3.3). The fact that a third emission component is not needed for data from the RGS instruments is explained by their different bandpass, i.e. 0.4 < E < 2.5 keV, whilst the ∼20 × 10⁶ K thermal component is mostly required for higher energies. We checked the consistency of the results obtained between EPIC and RGS data through a simultaneous fit of the spectra from the four instruments, i.e. EPIC-MOS1, EPIC-pn, RGS1 and RGS2. We used the same three-temperature model as in Section 3.3 and we obtained parameter values (see Table 6) very close to those presented in Table 4 for the simultaneous fit of EPIC-MOS1 and EPIC-pn spectra, with similar or slightly larger reduced χ².

4 X-RAY LUMINOSITY OF CYG OB2 #8A

4.1 Variability analysis from XMM–Newton data

The count rates obtained with the five instruments on board XMM–Newton for the four observations are quoted in Table 7. We observe significant variability of Cyg OB2 #8A on a time-scale of about 10 d, i.e. the typical separation between two pointings in our series. The largest variation is found between Observations 1 and 2, with an amplitude of about 20 per cent. We emphasize that the variations we observe for all instruments are correlated. To illustrate the variability observed between the different observations (see Fig. 4), we compared the count rates obtained in several energy bands in Observation 1 with those obtained in Observation 2 (left-hand panels), Observation 3 (middle panels) and Observation 4 (right-hand panels), respectively, for EPIC-MOS1 (upper part) and EPIC-pn (lower part). In each of these plots labelled (a)–(f), the upper section displays the count rates and the lower section shows the relative variability of the observed count rate. Although it is not shown here, we note that the same comparison was performed for the X-ray fluxes, estimated on the basis of the 3-T model with parameters given in Table 4. A plot of the relative variability of the observed X-ray flux of Cyg OB2 #8A for EPIC-pn has been presented by De Becker & Rauw (2005). We clearly see that there is a decrease in the X-ray count rate (flux) between the first and the second observation in the whole EPIC bandpass. The first and third observations appear to be very similar. In the case of the fourth observation, we see that the count rate decreases only in the hard energy band (above about 2.0 keV). We note that all the variability trends discussed here are consistent in both EPIC-MOS1 and EPIC-pn data, either if we consider count rates or observed fluxes.

We finally searched for short-term variability, i.e. within a single exposure. We have binned the event lists into 100, 200, 500 and 1000 s time intervals in four energy bands, respectively, 0.5–10.0, 0.5–1.0, 1.0–2.5 and 2.5–10.0 keV. We calculated the count rates in each time bin, along with their s.d., after subtraction of a background scaled according to the respective surface areas of the source and background regions (the same as used for the spectra extraction, see Section 2.1.1). GTIs were considered to compute the count rates using effective time bin lengths. A first inspection of the light curves does not reveal any significant variability correlated between the EPIC instruments on time-scales shorter than single exposures. This lack of significant variation is confirmed by variability tests applied to every light curve (χ² and pov-test as described by Sana et al. 2004).

4.2 Overall luminosity

On the basis of the best-fitting parameters presented in Table 4 for the three-temperature model, we have evaluated the fluxes between 0.5 and 10.0 keV for the four exposures. The observed, i.e. absorbed, fluxes are provided in the second-to-last column of Table 4. Considering a distance to Cyg OB2 #8A of 1.8 kpc (Bieging, Abbott & Churchwell 1989), we computed its unabsorbed X-ray luminosity, i.e. corrected for the ISM absorption, in the case of the simultaneous fit of EPIC-MOS1 and EPIC-pn data. The results are collected in Table 8. In this table, we also provide the L_X/L_bol ratio. On the basis of the bolometric luminosities given in Table 3, we also computed the expected intrinsic X-ray luminosity using the empirical relation proposed by Sana et al. (2006). Although this latter relation relies on a rather small sample of O-type stars compared to that of Berghöfer et al. (1997), we preferred to use this one because it was established in the same energy domain as for the present analysis, i.e. between 0.5 and 10.0 keV. We therefore obtain X-ray luminosities of 3.06 × 10³² and 2.20 × 10³² erg s⁻¹, respectively, for the primary and the secondary. The sum of these two quantities, i.e. L_X of 5.26 × 10³² erg s⁻¹, allowed us to calculate X-ray luminosity excesses ranging between about 19 and 28 (see column 4 of Table 8).

5 ARCHIVE X-RAY DATA

5.1 ROSAT-PSPC data

Cyg OB2 has been observed twice with the ROSAT-PSPC instrument. A first observation was performed on 1991 April 21 (sequence number rp200109n00, ∼3.5 ks), and the second one between 1993 April 29 and 1993 May 5 (sequence number rp900314n00, ∼19 ks). The latter consisted mainly of four exposures spread over about 5 d. The analysis of Waldron et al. (1998) revealed a significant variation of the soft X-ray flux (below 2 keV) between the 1991 and the 1993 observations, with the highest emission level observed in 1993. We retrieved the screened data from the archive and we used the xselect software to analyse the data of Cyg OB2 #8A. We extracted a light curve of this observation and we split it by applying time filters to obtain four separated data sets with effective exposure times of about 3–4 ks. We selected the source events within a 1 arcmin circular region. The background was selected in an annular region around the source region of the same area, excluding its intersection with a 30 arcsec circular region centred on a point source located to the north relative to Cyg OB2 #8A (RA = 20:33:13.9 and Dec. =+41:20:21.4, Equinox 2000.0). We used the xspec software to analyse the spectra and we obtained reasonable fits with a single temperature mekal model, with a kT of about 0.5–0.7 keV. We determined the count rates for each subexposure in the 0.4–2.5 keV energy band and we collected them in Table 9, along with the time of each exposure.

5.2 ASCA-SIS data

The Cyg OB2 association was observed with ASCA (Tanaka, Inoue & Holt 1994) during the performance verification phase on 1993 April 29 (sequence number 20003000, ∼30 ks). A first analysis of these data was reported by Kitamoto & Mukai (1996). These authors already pointed out the need to use two thermal emission components, with characteristic temperatures of the order of 0.6 and 1.5 keV, respectively, in order to model reasonably the data of Cyg OB2 #8A. We retrieved the raw data and processed them using the xselect software. Even though both gas-imaging spectrometers (GIS) and solid-state imaging spectrometers (SIS) were operated during the observation, we only used SIS data because of their better spatial resolution. We extracted the source events within a 2.21 and 2.95 arcmin radius circular region, respectively, for SIS0 and SIS1. In both cases, we selected the background events from a rectangular box located a few arcmin to the east of Cyg OB2 #8A. As the source region crosses over two CCDs, we constructed two RMFs, i.e. one for each CCD, and we obtained the effective RMF through a weighted sum of the two individual response matrices.

The spectral analysis was performed with the xspec software, and the best-fitting results were obtained with the three-temperature thermal model described in Section 3.3. We note that we obtained a lower value of the reduced χ² by replacing the third thermal component by a power law, but we estimate that this apparently better result is only due to the rather poor quality of the data in the hard part of the spectrum, unlikely to reveal the Fe K line clearly present in our XMM–Newton EPIC spectra. As the quality of the SIS0 data appeared to be significantly poorer than that of SIS1, we considered only the latter in our spectral analysis. The best-fitting parameters obtained with the 3-T thermal model between 0.5 and 10.0 keV are given in Table 10. The SIS1 spectrum and the corresponding model are presented in Fig. 5. From this model, we obtained an absorption corrected L_X of 1.86 × 10³⁴ erg s⁻¹ between 0.5 and 10.0 keV, leading to an X-ray luminosity excess of about 35. We finally note that the observed count rate in the same energy band is 0.331 ± 0.004 cts s⁻¹ for SIS1.

6 DISCUSSION

6.1 Orbital modulation of the X-ray flux

6.1.1 Observational material

As Cyg OB2 #8A is a binary system, one could wonder whether the existing X-ray observations reveal a modulation of the X-ray flux. This issue was first addressed by De Becker et al. (2005a) where the results from several X-ray observations (ROSAT and ASCA) were combined to obtain a phase-folded light curve, on the basis of the ephemeris published by De Becker et al. (2004c). The light curve suggested a phase-locked modulation of the X-ray flux, probably due to the combined effect of the variation of the absorption along the line of sight and of the X-ray emission itself as a function of orbital phase. Because of inconsistencies between ROSAT-HRI and ROSAT-PSPC count rates and because of a poor sampling of the orbital cycle, this preliminary light curve did not allow us to draw a firm conclusion. However, the light curve presented by De Becker et al. (2005a) suggests clearly that the ROSAT-HRI count rates show a phase modulation similar to that of the PSPC data.

Using our four XMM–Newton observations, along with the results from archive ROSAT-PSPC and ASCA-SIS1 data,³ we constructed a new light curve. To compare the count rates from the different instruments in a consistent way, we used the 3-T model with the parameters obtained for the simultaneous fit of EPIC-MOS1 and EPIC-pn data for Observation 1, and we convolved it with the respective response matrices of ROSAT-PSPC and ASCA-SIS1 to obtain faked spectra. We obtained count rates of 0.250 ± 0.004 and 0.270 ± 0.001 cts s⁻¹, respectively, for both instruments. On the basis of these values, and of the count rates obtained in Section 5, we compared the X-ray emission level from all observations after normalization with respect to the XMM–Newton Observation 1. The normalized X-ray count rates obtained this way are plotted as a function of the orbital phase in Fig. 6. We note that this light curve does not suggest any large error on the orbital parameters derived by De Becker et al. (2004c), considering the short period and the large time interval separating some of the observations discussed here. Over the time range between 1991 and 2004, an error of 0.040 d on the period would indeed lead to an error on the orbital phase of the order of 0.4. This suggests that the error of 0.040 d given for the 21.908 d period might be a somewhat conservative value, reinforcing our confidence in the orbital parameters proposed by De Becker et al. (2004c).

The light curve of Fig. 6 presents a maximum at an orbital phase close to 0.75. We note the relative consistency of the various observations that contribute to the rather steep decrease in the X-ray emission between phases 0.75 and 1.0. We note however a discrepancy between the EPIC and PSPC points close to phase 0.9. Unfortunately, the lack of observations between phases 0.0 and 0.4 prevents us from constraining the position of the minimum, likely located shortly after phase zero. A somewhat more detailed view of the variability of the observed X-ray flux can be obtained on the basis of Fig. 7. The X-ray fluxes in the soft (0.5–2.0 keV) and hard (2.0–10.0 keV) bands are individually plotted as a function of the orbital phase in the two upper panels. In order to investigate the variability in the hard part of the X-ray spectrum, we used only XMM–Newton and ASCA data. The separation between the two stars in Cyg OB2 #8A and the position angle (PA) ⁴ are also provided. The decrease of the X-ray flux between apastron and periastron is particularly obvious in the hard energy band. The maximum close to phase 0.75 suggested by Fig. 6 comes mainly from the soft band, as shown by the upper panel of Fig. 7.

In addition to the decrease in the X-ray flux, the comparison between the XMM–Newton Observations 1 and 2 reveals a significant decrease in the characteristic plasma temperature of the hottest thermal component (see Table 4). As these two observations are, respectively, close to phases 0.5 (largest separation) and 0.0 (smallest separation), the decrease in the post-shock plasma temperature could be related to the decrease in the pre-shock velocity at the position of the collision zone. For this rather close binary system, the winds will indeed not yet have reached their terminal velocity before they collide. Using our estimate of the orbital (De Becker et al. 2004c) and wind (Table 3) parameters and adopting a β= 1 wind velocity law, we calculate the distance between the stagnation point (i.e. the location of the wind–wind interaction on the binary axis) and the centre of each star. Following an estimate of the inclination angle of the system,⁵ i.e. about 32°± 5°, we infer absolute distances of about 95 and 81 R_⊙, respectively, for the primary and the secondary at apastron (respectively 58 and 50 R_⊙ at periastron). Using the corresponding pre-shock velocities, we expect the post-shock temperature to vary (on-axis) between about 35 × 10⁶ and 25 × 10⁶ K from phases 0.5 to 0.0⁶ (see Table 11). According to the values quoted in Table 4, the observed temperature for the hottest thermal component is about 21 × 10⁶ and 18 × 10⁶ K, respectively, for Observations 1 and 2. The fact that the predicted and observed temperatures are different (about a factor of 1.4–1.7) may suggest that the effective pre-shock velocity is lower than predicted by about a factor of 1.2–1.3. However, as the EPIC bandpass does not extend further than 10 keV, we should bear in mind that our spectral fits may underestimate the characteristic temperature of the post-shock plasma. Moreover, some fraction of the emission is expected to come from off-axis where the shocks are oblique, leading the averaged observed emission to have a lower characteristic temperature than anticipated at the stagnation point. Finally, some radiative inhibition (Stevens & Pollock 1994) may be at work close to the line of centres of the binary system, therefore reducing the speed of the wind flows before they collide. As a result, considering the numbers given above, we consider that the agreement between the temperatures estimated from the fit of our X-ray spectra and that derived from the geometry and hydrodynamic of the colliding winds in Cyg OB2 #8A is very satisfactory.

Some explanation for the variability described in Section 4.1 (see Fig. 4) and in Figs 6 and 7 can be given on the basis of the different orbital phases of our XMM–Newton observations. In the case of a binary system like Cyg OB2 #8A, we can expect some variability in the X-ray domain mainly for two reasons described as follows.⁷

• The variation of the line-of-sight absorption along the orbital cycle, likely to affect mostly the softer part of the spectrum (i.e. below about 2.0 keV).

• The variation of the separation between the two stars, as the orbit is eccentric, likely to affect the physical conditions in the colliding wind zone. If some variability is observed in the harder part of the spectrum, it should most probably come from this process.

Considering the strong differences between Observations 1 and 2 in the whole EPIC bandpass (parts a and d of Fig. 4), both factors might play a significant role. Considering only XMM–Newton results, this decrease is more spectacular in the hard part of the spectrum (Figs 4 and 7). This should be due to the variation of the separation between the two stars. The fact that Observation 2 occurs when the primary is ‘in front’ of the secondary (PA close to 0°), i.e. at an orbital phase where the absorption should be larger, may explain the decrease in the soft part of the spectrum. The parameters quoted in Table 4 point indeed at a higher local absorption column in the case of Observation 2 as compared to the first one. The same trend is also suggested by the fits of RGS spectra (see Table 6). The lack of significant variability between Observations 1 and 3 (parts b and e of Fig. 4) is compatible with the fact that they fall nearly at the same phase close to apastron. Finally, the decrease in the X-ray emission observed above 2.0 keV in Observation 4 (parts c and f of Figs 4 and 7) might be explained by the decrease in the separation, lowering the X-ray emission from the collision zone. The somewhat higher local absorption in the case of Observation 4 suggested by the wind absorption parameter given in Table 4 may be an artefact of the fit, where the apparent higher absorption is compensated for by the larger normalization parameter of the same observation.

6.1.2 Hydrodynamic simulations

Considering the light curves presented above (see Figs 6 and 7), it is obvious that the X-ray emission from Cyg OB2 #8A presents orbital modulations which are worth to be further investigated. For a better understanding of the physical processes responsible for this orbital modulation, we compared our XMM–Newton observations with the predictions of detailed hydrodynamic simulations using the same approach as Sana et al. (2004). In this method, the VH-1 numerical code based on the piece-wise parabolic method (Colella & Woodward 1984), is used to solve the partial differential equations of hydrodynamics followed by a remap on to a fixed grid after each time-step. A radiative energy loss term is included in the hydrodynamic equations in order to treat the cooling self-consistently. For details see Stevens et al. (1992), Pittard & Stevens (1997) and Sana et al. (2004).

Briefly, two spherically symmetric ionized winds of constant velocity are assumed, resulting in an axisymmetric geometry around the lines of centres where the hydrodynamic problem is reduced to a two-dimensional flow. The orbital motion is neglected. As the acceleration of the wind is not taken into account, the pre-shock velocity is estimated on the basis of a classical velocity law (β= 1.0), for a terminal velocity estimated to be 2.6 times the escape velocity (Vink et al. 2000). We adopted a square grid size of 300 × 300 cells, corresponding to physical distances of 1.5 × 10¹³×1.5 × 10¹³ cm and we let the flow evolve a time long enough so that the system relaxes from the initial conditions. At each step, grids of density, pressure, radial and axial velocities are obtained, and therefore grids of temperature can be calculated. The X-ray emission from the system is evaluated by summing up the emissivity of each cell of the grid at each time-step. The column of absorbing material was computed following a three-dimensional geometry taking into account the inclination of the system and the orientation of the line of sight with respect to the system (see Sana et al. 2004, for details).

We chose three different system configurations to evaluate the evolution of the X-ray flux between 0.5 and 10.0 keV along the orbital cycle: (i) apastron (φ= 0.5), (ii) intermediate (φ= 0.75) and (iii) periastron (φ= 0.0). The pre-shock velocities, the post-shock temperatures and the cooling parameter (χ) estimated in each case are given in Table 11 for both stars of the system. In the three cases, the collision zone relaxes from initial conditions after about 2500 steps, i.e. about 4 × 10⁵ s, and turns out to be highly unstable. This instability comes most probably from the cooling as the shock of the primary wind, and that of the secondary at periastron, are radiative (Stevens et al. 1992), as suggested by the values quoted for the cooling parameter (χ) in Table 11. We mention also that it may also partly result from shear instabilities as the winds in Cyg OB2 #8A have different velocities.

The mean and s.d. of the X-ray luminosity are computed over a large number of time-steps after relaxation from the initial conditions. The averaged predicted X-ray luminosities, absorbed by both the wind and the interstellar material (using N_H= 0.94 × 10²² cm⁻², see Section 3.2), are given in the second column of Table 12. The comparison of the predicted and observed X-ray luminosities reveals that the theoretical modelling leads to X-ray luminosities in reasonable agreement (within about a factor of 2) with the observed luminosities. The simulations predict also a minimum of the X-ray luminosity at periastron that is consistent with what we observe (see Fig. 6). The rather high emission level observed close to phase 0.75 is also predicted. We mention that the observed large amplitude of the variation between phases 0.75 and 0.0 is also predicted, even though this amplitude is larger for the predicted values. In addition, the relative emission levels obtained at the three orbital phases selected for our simulations are in reasonable agreement with the fluxes plotted in Fig. 7 in the hard energy band, i.e. the spectral domain where the colliding winds are mainly expected to produce X-rays. In particular, the rather similar emission levels obtained at phases 0.5 and 0.75 is in agreement with the observed hard X-ray fluxes. In addition, we mention that the results of the simulations are in agreement with the large X-ray overluminosity revealed by the observations. The rather good agreement of the general behaviour between the predicted and observed variabilities confirms that the observed orbital modulations of the X-ray flux are produced by the line-of-sight absorption and by the variation of the separation between the two stars, as these two effects are dominant in our simulations. This result provides strong support to the scenario where the modulations of the X-ray flux come from the colliding winds.

6.2 Non-thermal emission

Cyg OB2 #8A is known as a bright non-thermal radio emitter (Bieging et al. 1989). The fact that it is a confirmed binary system suggests that the non-thermal emission is produced in the wind interaction zone. This has recently been confirmed by Blomme (2005) who presented a phase-folded radio light curve built on the basis of new data showing a strong phase-locked variability of the radio flux density. This is in agreement with the results of the study of Van Loo, Runacres & Blomme (2006) who showed that the production of the observed non-thermal emission from isolated stars is unlikely. In this scenario, electrons are accelerated through the first-order Fermi mechanism up to relativistic velocities (Eichler & Usov 1993). This population of relativistic electrons is expected to interact with the local magnetic field, likely originating from the two stars, to produce synchrotron radiation in the radio domain (see e.g. Dougherty et al. 2003). New developments in the modelling of non-thermal radio emission from massive binaries are described by Pittard et al. (2006). These recent models take into account several physical effects like free–free absorption, the Razin effect and IC cooling.

The latter process is particularly interesting in the sense that the cooling of relativistic electrons by UV photospheric photons through IC scattering is likely to produce a non-thermal high-energy emission component. However, we did not find any evidence for a power-law emission component attributable to a non-thermal emission in our XMM–Newton-EPIC spectra of Cyg OB2 #8A. This was not unexpected as any putative non-thermal emission component would probably be overwhelmed by the strong thermal emission arising from the wind–wind interaction zone. The unlikelihood of the detection of a power-law emission below 10.0 keV in the case of non-thermal radio emitting massive binaries was indeed pointed out by De Becker et al. (2005b). However, such a non-thermal emission could possibly be detected in very hard X-rays and soft γ-rays, where the spectra are not expected to be contaminated by thermal emission. For instance, De Becker et al. (2005a) estimated that the possible contribution of Cyg OB2 #8A to the γ-ray emission from the yet unidentified EGRET source 3EG J2033+4118 could amount up to about 10 per cent. However, these estimations need to be refined using more sophisticated models (see e.g. Pittard & Dougherty 2006). In addition, the observation of the Cyg OB2 region with the IBIS(ISGRI) imager onboard the INTEGRAL observatory did not allow to detect any high-energy emission putatively associated to the massive stars in Cyg OB2 (De Becker 2005).

7 SUMMARY AND CONCLUSIONS

In this paper, we presented the results of four XMM–Newton observations of the massive binary Cyg OB2 #8A separated by about ten days from each other. The best fits of the X-ray spectra were obtained using a three-temperature thermal model, with characteristic plasma temperatures for the three components of about 3 × 10⁶, 9 × 10⁶ and 20 × 10⁶ K. The hottest component is most probably attributed to the thermal X-ray emission from the collision zone between the winds of the two stars. The nature of the X-ray emission appears to be purely thermal, and we failed to fit power-law models to the hard part of the spectrum. This is in agreement with the idea proposed by De Becker et al. (2005b) that non-thermal radio emitters are not likely to display a non-thermal emission component in their X-ray spectrum below 10.0 keV. The X-ray luminosity is very high (about 10³⁴ erg s⁻¹, leading to an X-ray luminosity excess of 19–28 during our observations). This rather high X-ray emission level is in disagreement with the statement by Waldron et al. (2004) that the X-ray emission in Cyg OB2 #8A is probably not different from that expected from isolated O-stars.

The analysis of high spectral resolution RGS data reveals an absorbed spectrum with prominent lines mainly from Si xiv, Si xiii, Mg xii, Mg xi, Ne x, Ne ix and Fe xvii. The comparison of the spectra obtained at different orbital phases suggests a possible variability of some line profiles. A better phase coverage of the orbital cycle with much longer exposure times is needed to address this issue in detail.

We folded the count rates from our XMM–Newton observations, along with those from archive ROSAT-PSPC and ASCA-SIS observations, with the ephemeris given by De Becker et al. (2004c). We observe a variability of the count rate and of the X-ray flux of about 20 per cent between apastron and periastron. The light curve as a function of the orbital phase points to a maximum at about phase 0.75, and suggests a minimum shortly after the periastron passage. The shape of the light curve might be explained by the combined effect of absorption and varying separation along the orbital cycle, this latter factor being responsible for the lower pre-shock velocity reached by the winds at periastron with respect to apastron. Using hydrodynamic simulations, we also find a lower emission level close to phase 0.0. The amplitude of the observed variation is less than predicted. This contrast may be reduced by including shock modification (Pittard & Dougherty 2006), by reducing the value of β for the wind velocity law, or even by slightly changing the inclination angle. As a main conclusion, we state that the strong phase-locked variability — along with the spectral shape — of the X-ray emission of Cyg OB2 #8A revealed by our investigation points undoubtedly to X-ray emission dominated by colliding winds.

In the future, Cyg OB2 #8A appears to be an ideal target for the Wide-band X-ray Imager (WXI) and the Soft Gamma-ray Detector (SGD) onboard the Next-generation X-ray Telescope satellite NeXT (Takahashi et al. 2004), whose sensitivity is expected to be significantly better than that of INTEGRAL. On the other hand, the results from the radio monitoring of Cyg OB2 #8A (Blomme 2005) might be used in parallel with state-of-the-art models to evaluate the non-thermal emission level in the high-energy domain.

ACKNOWLEDGMENTS

Our thanks go to Alain Detal (Liège) for his help in installing the sas. The Liège team acknowledges support from the Fonds National de la Recherche Scientifique (Belgium) and through the PRODEX XMM and INTEGRAL Projects. This research is also supported in part by contract P5/36 ‘Pôle d'Attraction Interuniversitaire’ (Belspo). This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France and NASA's ADS Abstract Service.

REFERENCES