© ESO 2020

1. Introduction

The radiation fields of massive stars drive powerful stellar winds associated with huge mass-loss rates (typically 10⁻⁷ to 10⁻⁵ M_⊙ yr⁻¹) and highly supersonic velocities (typically 2000 km s⁻¹). In binary systems consisting of two massive stars, the collision of their stellar winds produces an interaction zone contained between two oppositely faced hydrodynamical shocks separated by a contact discontinuity (e.g. Stevens et al. 1992). The presence of this interaction zone leads to a number of observational signatures that can range from radio waves into the high-energy domain (see e.g. Rauw 2013). Thermal X-ray emission is created by the plasma from the interaction zone which may be heated up to 10⁷ K through the kinetic energy of the winds (Stevens et al. 1992). In some systems, this results in a prominent X-ray emission in the 0.5−10 keV band that varies with orbital phase, as a result of the changing optical depth along the line of sight and/or of the changing orbital separation in eccentric systems (for a review, see Rauw & Nazé 2016, and references therein).

Aside from the heating of post-shock plasma, hydrodynamic shocks in colliding wind binaries (CWBs) can also produce a population of relativistic particles via diffusive shock acceleration through the first order Fermi mechanism. A subset of the CWBs indeed display synchrotron radio emission (Benaglia 2010, and references therein), which is produced via the interaction of relativistic electrons with a magnetic field (Eichler & Usov 1993; Pittard et al. 2006; Pittard & Dougherty 2006). This non-thermal radio emission is often variable as a result of changing line-of-sight optical depth and, in eccentric systems, also of changing shock strength (Blomme et al. 2010, 2013).

The presence of a population of relativistic electrons, along with the enormous supply of stellar photospheric ultraviolet (UV) photons, especially in short-period massive binaries, should result in a strong inverse Compton (IC) scattering emission in hard X-rays and soft γ-rays (Pollock 1987; Chen & White 1991; Pittard & Dougherty 2006; Reimer et al. 2006; del Palacio et al. 2016). X-ray observations of O-star binaries that display non-thermal radio emission however revealed no clear evidence of non-thermal X-ray emission in the 0.5−10 keV energy domain (Rauw et al. 2002; De Becker et al. 2006). At these energies, the putative non-thermal X-ray emission is overwhelmed by the strong thermal emission from the wind interaction which involves plasma at temperatures up to kT = 2 keV. However, at energies above 10 keV, where the thermal emission rapidly declines, the conditions for detections of a non-thermal emission are more favourable. Indeed, a detection of the very massive CWB η Carinae above 20 keV has been reported with INTEGRAL and Suzaku, and most recently with NuSTAR (Leyder et al. 2008, 2010; Sekiguchi et al. 2009; Hamaguchi et al. 2018)¹. The CWB was also detected up to giga-electronvolt energies with the AGILE and Fermi telescopes (Tavani et al. 2001; Abdo et al. 2010). However, detection of the hard X-ray emission between 10 keV and 20 keV with NuSTAR is somewhat uncertain (Hamaguchi et al. 2018; Panagiotou & Walter 2018).

In the present work, we consider the CWB Cyg OB2 #8A (≡BD+40° 4227A, Schulte 8A), an eccentric (e ∼ 0.2) O6 I + O5.5 III binary with an orbital period of 21.9 days (De Becker et al. 2004). It was one of the first massive stars found to be an X-ray emitter (Harnden et al. 1979). Since that pioneering detection, the system has been shown to display phase-locked variations in its X-ray emission (Cazorla et al. 2014, and references therein). Cyg OB2 #8A is also known for its strong non-thermal radio emission (Bieging et al. 1989; Blomme et al. 2010). Owing to the compactness of its orbit, we expect IC scattering to be particularly efficient in this system. De Becker et al. (2007) analysed the INTEGRAL/IBIS data of the Cygnus region obtained during the two first Announcements of Opportunity, but failed to detect any emission directly associated with Cyg OB2. These authors inferred 3σ upper limits on the count rates at the position of the unidentified EGRET source 3EG 2033+4118; these count rates, assuming a photon index Γ = 1.5, convert into upper limits on the flux of 6.1 × 10⁻¹², 4.2 × 10⁻¹², and 4.0 × 10⁻¹¹ erg s⁻¹ cm⁻² in the 20−60, 60−100, and 100−1000 keV energy bands, respectively. In view of the spatial resolution of IBIS, these upper limits apply to the combined emission from all putative emitters within a radius of 6′ around the position of 3EG 2033+4118, which includes Cyg OB2 #8A as well as Cyg OB2 #9, another non-thermal radio emitting O-star binary (Rauw 2004; Blomme et al. 2013, and references therein). In this work, we analyse the overall set of available INTEGRAL observations performed up to now and refine the upper limit on the hard X-ray emission from Cyg OB2 #8A. The NuSTAR satellite is able to observe in hard X-rays with a much better angular resolution than INTEGRAL. We also analyse the first NuSTAR observations of Cyg OB2 #8A. We first present the observations (Sect. 2) and revise the orbital solution of the system (Sect. 3). We then present the X-ray spectral analysis (Sect. 4) and discuss the results (Sect. 5). The summary and conclusions of this study are provided in Sect. 6.

2. Observations

2.1. NuSTAR

Cyg OB2 #8A was observed with the NuSTAR satellite on 25 and 26 August 2018. NuSTAR (Harrison et al. 2013) features two focussing telescopes that operate in the energy domain from 3 keV to 79 keV. The focussing optics provide a superior sensitivity in the hard X-ray domain compared to previous missions that relied on coded masks. We analysed these observations with the NuSTAR Data Analysis Software package (nustardas v1.8.0) built in HEASOFT software (v6.24). Stray light contamination from the off-axis source Cygnus X-3 is high in the field of view of both observations (see Fig. 1). Moreover, part of the observations were taken close to the South Atlantic Anomaly (SAA). Great caution must thus be taken to correctly take the high X-ray background into account. The images reveal two sources: Cyg OB2 #8a and Cyg OB2 #9. We focus on the former target because the latter source is weaker and is more heavily affected by the stray light from Cyg X-3.

Fig. 1. NuSTAR RGB image (red = 1.6−5.6 keV, green = 5.6−13.6 keV, blue = 13.6−21.6 keV) of the Cyg OB2 region taken on 2018 August 26 (FPMA data only).

The events from Cyg OB2 #8A were extracted over a 30″ radius circular area centred on the optical position of Cyg OB2 #8A (Gaia Collaboration 2018). As recommended by the NuSTAR team², the background region is defined as an annulus of 50 and 80″ inner and outer radius centred on Cyg OB2 #8A. Based on the configuration of the observations, we reduced the level-1 data using nupipeline (v0.4.6) with the keywords tentacle=yes and saamode=strict. We then extracted the events and spectra from focal plane modules (FPMs) A and B for the source and background regions. The corresponding response files (ARF and RMF) were created with numkarf and numkrmf. The source spectra were grouped with a minimum signal-to-noise ratio of 4 in the ISIS software (Houck 2013).

2.2. XMM-Newton

Twelve XMM-Newton observations were used in this study (Table A.1). Seven of these observations were previously analysed by Cazorla et al. (2014).

For all observations but 0677980601, the European Photon Imaging Cameras (EPIC; Turner et al. 2001; Strüder et al. 2001) were operated in full-frame mode. For observation 0677980601, the metal oxide semi-conductor (MOS) cameras were operated in large window mode with the field of view centred on Cyg OB2 #9. Cyg OB2 #8A was thus only observed by the pn camera for this observation. Given the optical brightness of our target (V = 8.98), the medium optical filter was used to reject optical and UV photons. The EPIC data were processed with the emchain and epchain tasks from the Science Analysis Software (SAS) package (version 17.0; Current Calibration files as of 2018 June 22) to extract the event lists for the MOS and pn cameras. We then selected the good time intervals (GTI) defined when the total count rate in the 0.2−10 keV energy range over the full detector is lower than 0.009 for pn and 0.004 count s⁻¹ arcmin⁻² for MOS. We selected the X-ray events by keeping only the single and double events (PATTERN≤4) for the pn camera and the single, double, triple, and quadruple events (PATTERN≤12) for the MOS cameras. Finally, we rejected the dead columns and bad pixels using the bit masks FLAG==0 for pn and #XMMEA_SM for MOS. The source extraction region was the same as for the NuSTAR analysis. This region allows us to extract 90% and 85% of the flux at 1.5 keV on-axis for the pn and MOS cameras, respectively. The background events were extracted from an annulus of 50 and 100″ inner and outer radius centred on Cyg OB2 #8A. The X-ray point sources detected in the background region using the SAS task edetect_chain were filtered out. We extracted the spectra for the source and background regions and built the corresponding ARF and RMF. The source spectra were grouped with a minimum signal-to-noise ratio of 7 using the task specgroup.

2.3. Swift

Cyg OB2 #8A was observed 73 times with Swift (Gehrels et al. 2004) in photon counting mode from 2006 December 15 to 2018 August 27 (Table A.1). Six of these observations were previously analysed by Cazorla et al. (2014).

We used the HEASOFT task XRTPIPELINE (v0.13.4) and the calibration files released on 2018 July 10 to reject the hot and bad pixels and select the grades between 0 and 12. We used XSelect (v2.4e) to extract events from Cyg OB2 #8A. The definition of the source region is the same as for the XMM-Newton and NuSTAR data. For the spectral analysis, the ARFs corresponding to the source position on the detector were computed with xrtmkarf (v0.6.3). The spectra were grouped with a minimum of 20 counts per bin using grppha. The spectra having less than three bins were rejected.

2.4. INTEGRAL

The INTEGRAL instrument observed the Cygnus region thousands of times with the IBIS imager (Ubertini et al. 2003). We selected all public observations performed with ISGRI up to now (i.e. from rev. 11 to 1977), where the Cygnus region is located in the fully coded field of view (i.e. with an off-axis angle lower than 480′), leading to 1736 observations corresponding to 3.1 Ms of exposure time³. We analysed these observations using the off-line scientific analysis (OSA) software (v11.0) provided by the Integral Science Data Center (Courvoisier et al. 2003). We built a mosaic image of the overall INTEGRAL observations in the three following energy bands: 20−60, 60−100, and 100−200 keV. Following De Becker et al. (2007), we then extracted the spectra by considering a Gaussian with a half width at half maximum of 6′ at the position of the EGRET source 3EG J2033+4118 corresponding to the Cygnus OB2 region for each of the mosaic image with the mosaic_spec task from the OSA software. The response files and background spectra were created using the spe_pick command to average the files proportionally to the exposures from the different observing periods.

2.5. Optical spectroscopy

Optical spectroscopy of Cyg OB2 #8A was obtained with the Aurélie spectrograph (Gillet et al. 1994) at the 1.52 m telescope of the Observatoire de Haute Provence (OHP; France). The data were taken during two observing campaigns of six nights each in June 2016 and September 2017. Aurélie was equipped with a 2048 × 1024 CCD with a pixel size of 13.5 μm squared. We used a 600 l mm⁻¹ grating providing a reciprocal dispersion of 16 Å mm⁻¹. The resolving power, measured on the Thorium-Argon calibration exposures, was 7000 over the wavelength range from 4440 to 4890 Å. Typical integration times were 1.5 h. The data were reduced using the MIDAS software (version 17FEBpl 1.2). In this way, we obtained eight new spectra (see Table 1).

3. Revised orbital solution

We measured the radial velocities (RVs) of the He Iλ 4471 absorption line on our new optical spectra. Because of the severe blending of the primary and secondary components at most orbital phases, the MIDAS two-Gaussian fit routine deblend/line only provided a reliable fit for two observations (HJD 2 458 002.538 and HJD 2 458 004.516). We therefore adopted the same approach as De Becker et al. (2004) to measure the RVs on the other spectra: we fitted each spectrum with a combination of two Gaussians in which the widths and relative intensities were fixed to the values obtained via deblend/line for the observation taken on HJD 2 458 002.538. In this procedure, only the positions of the lines and the overall line intensity were varied. The resulting RVs are listed in Table 1. The typical errors on these new RVs are around 10 km s⁻¹ for the primary star, but can exceed 20 km s⁻¹ for the secondary at those phases where the lines are most heavily blended.

We then combined the new RV points with those of De Becker et al. (2004) to compute a revised orbital solution with the Liège Orbital Solution Package (LOSP) code (Sana et al. 2006). The results are given in Table 2 and illustrated in Fig. 2. Compared to the solution of De Becker et al. (2004), we find a lower eccentricity (0.18 vs. 0.24) and a slightly larger mass ratio (1.26 vs. 1.16).

Fig. 2. Revised orbital solution of Cyg OB2 #8A. Top panel: primary and secondary RVs are indicated by the filled and open symbols, respectively. The blue dots correspond to the RVs from De Becker et al. (2004), whereas the red squares indicate the new RVs. The solid and dashed lines illustrate the best-fit RV curve (see Table 2). Bottom panels: relative orbital separation and position angle as a function of orbital phase. A value of 0° for the position angle corresponds to the conjunction with the primary star in front.

Comparing the minimum masses that we inferred from our revised orbital solution with the typical masses of stars of same spectral type quoted by Martins et al. (2005), we estimated an orbital inclination of (33.8 ± 1.0)° for the primary and (32.9 ± 1.0)° for the secondary. We thus adopted i = (33.4 ± 0.7)° as a reasonable estimate of the inclination. This then leads to a semi-major axis of the orbit of a = 142.8 R_⊙ = 9.9 × 10¹² cm. Bailer-Jones et al. (2018) computed a distance for Cyg OB2 #8A of 1.5 ± 0.1 kpc. With an optical brightness ratio (primary/secondary) close to 2 and the bolometric corrections corresponding to the spectral types derived by De Becker et al. (2004), we obtained luminosities of (7.2 ± 0.1) × 10⁵L_⊙ for the primary star and (4.3 ± 0.1) × 10⁵L_⊙ for the secondary star. Using the spectral type – effective temperature calibration of Martins et al. (2005), these bolometric luminosities imply stellar radii of (20.6 ± 0.2) R_⊙ for the primary and (14.0 ± 0.1) R_⊙ for the secondary. Both stars therefore remain well inside their Roche lobes at all orbital phases.

4. X-ray spectral analysis

The goal of this study is to provide constraints on the properties of a putative hard non-thermal X-ray emission from Cyg OB2 #8A. To do so, we analysed the NuSTAR and INTEGRAL spectra to test the presence of a power-law component at high energies. This requires a good characterisation of the X-ray spectrum at lower energies as well because some of the thermal emission that dominates below 10 keV could extend into the higher energy domain. All our spectral fits were done under the X-ray spectral fitting package (XSPEC; v12.10.0e).

The X-ray spectrum of Cyg OB2 #8A is heavily absorbed at energies below 0.8 keV by the large column density of interstellar material and by the material of the stellar winds. For the interstellar medium, we followed Cazorla et al. (2014) in adopting a total neutral hydrogen column density of N_H = 0.91 × 10²² cm⁻² corresponding to the Bohlin et al. (1978) relation with E(B−V) = 1.56. This value is higher than the N_H I value of (0.63 ± 0.07) × 10²² cm⁻² inferred by Herrero et al. (2001) from the interstellar Ly α line measured on HST/STIS spectra of Cyg OB2 #8A. However, the Ly α line only traces H I, whereas the total neutral hydrogen column density must also account for H₂. Assuming that the average N_H2/N_H I ratio of Bohlin et al. (1978) holds for the sightline towards Cyg OB2 #8A, the results of Herrero et al. (2001) translate into N_H = (0.84 ± 0.10) × 10²² cm⁻² in reasonable agreement with the value we adopted.

4.1. Improvement of the spectral model of Cazorla et al. (2014)

Cazorla et al. (2014) fitted seven of the XMM-Newton spectra with an absorbed optically thin thermal plasma model consisting of three components produced by the astrophysical plasma emission code (APEC, Smith et al. 2001). The parameters of each APEC component are the temperature of the plasma (kT) and the normalisation of the emission (n). The three APEC models are absorbed by both the interstellar medium reproduced by wabs and the stellar winds modelled by phabs.

The wabs model uses the abundances determined by Anders & Grevesse (1989). However, these abundances were updated by Wilms et al. (2000) leading to lower metal abundances. We thus improved the Cazorla et al. (2014) model with the TBnew model, the abundances of Wilms et al. (2000), and the updated cross sections of Verner et al. (1996). These abundances were also used for phabs and the APEC emission. We also took the effects of the dust scattering along the line of sight into account by adding the dustscat component (Predehl & Schmitt 1995) whose hydrogen column density was set to the hydrogen column density of the TBnew component divided by 1.5 (Nowak et al. 2012).

We fitted the spectra using the Markov chain Monte Carlo (MCMC) method. We used the XSPEC_emcee⁴ programme designed by Jeremy Sanders, which allows the MCMC analysis of X-ray spectra in XSPEC using emcee⁵, an extensible, pure-Python implementation of the Goodman & Weare (2010) MCMC ensemble sampler. This method uses a number of “walkers”, which evolve independently from each other in the parameter space reducing the autocorrelation time. We followed the same methodology as adopted by Cazorla et al. (2014). We first fitted the XMM-Newton spectra for each observation and checked the constancy of the temperatures of the three APEC components with phase. Since they agree within the errors, we then set the temperatures to their mean values (0.34, 0.83 and 2.13 keV) and fitted the XMM-Newton spectra again to determine the normalisation of the APEC models. The normalisation of the softest component (which is probably mainly due to X-ray emission from the individual stars) is almost constant with phase and has very large error bars. We thus fixed it to its mean value: 2.11 × 10⁻² cm⁻⁵. We finally fitted XMM-Newton spectra adopting as free parameters the phabs hydrogen column density and the normalisations n₂ and n₃.

We also grouped the Swift spectra by phase interval; we used an interval of 0.1 beginning at ϕ = 0.05. We simultaneously fitted the Swift spectra from each phase group with the same model as for XMM-Newton. The variation in the spectral parameters with phase is shown in Fig. 3. The best-fitting parameters of all the models we tested are summarised in Table B.1.

Fig. 3. Evolution with phase of the absorption and normalisation factors from the fit with three APEC components whose temperatures and n1 component were fixed. The red asterisks, blue circles, and green squares indicate the spectral parameters from the Swift, XMM-Newton, and NuSTAR fits, respectively.

4.2. Search for a non-thermal hard X-ray component

The putative non-thermal photons emitted by the wind shock region mainly contribute to the flux observed at high energies, i.e. above 10 keV. To test the presence of such a component, we thus investigated the NuSTAR spectra. The two NuSTAR observations were taken close to periastron passage (phase 0). We thus interpolated the values of the wind hydrogen column density and the n₂ and n₃ components to the phases of the NuSTAR observations and applied the so-constructed models to the NuSTAR spectra. For the ϕ = 0.028 spectra, the interpolated model predicts a hard X-ray flux exceeding what is observed with NuSTAR. We thus let the n₃ component free, leading to a best-fit value of n₃ = (3.3 ± 0.3) × 10⁻³ cm⁻⁵ (Table B.1 and top panel of Fig. 4). This value is small compared to what is observed at phases 0.019 and 0.042 with XMM-Newton. On the other hand, for the ϕ = 0.085 NuSTAR spectra, the observed flux above 10 keV is underestimated by the interpolated model. We thus added a power-law component to the model and fitted NuSTAR spectra letting only the power-law parameters free. This leads to a negative photon index and a flux consistent with zero. The F-statistic shows that the addition of the power-law component does not significantly improve the fit (p = 0.25). We thus let the n₃ component of the thermal model free leading to a best fit of n₃ = (3.7 ± 0.4) × 10⁻³ cm⁻⁵, i.e. close to the value observed at phase 0.056 with XMM-Newton (bottom panel of Fig. 4).

Fig. 4. Best-fitting dustscat × TBnew × phabs × (apec+apec+apec) models of the ϕ = 0.028 (top panel) and ϕ = 0.085 (bottom panel) NuSTAR spectra.

4.3. Upper limits on the flux of a non-thermal component

As becomes obvious from the various trials above, there is currently no unambiguous detection of a hard, non-thermal X-ray emission in Cyg OB2 #8A. We can nevertheless use the INTEGRAL and NuSTAR spectra to establish upper limits on the flux of such a component.

Using the spectra extracted from the mosaic images of the 3.1 Ms of INTEGRAL observations, we constrained the 3σ upper limit on the count rate of the 3EG 2033+4118 source in the three energy bands to 3.5 × 10⁻² count s⁻¹ in 20−60 keV, 6.3 × 10⁻³ count s⁻¹ in 60−100 keV, and 4.4 × 10⁻³ count s⁻¹ in 100−200 keV. We then determined the flux in the 20–30 keV energy band reproducing the total count rate of 4.5 × 10⁻² count s⁻¹ over the 20−200 keV energy band. Because the hard X-rays in this band should mostly stem from IC scattering, the total count rate was converted to flux by modelling the emission with a power-law characterised by different values of spectral index Γ. The resulting 3σ upper limits on the non-thermal flux of 3EG 2033+4118 in the 20−30 keV energy band are reported in Table 3 for each value of Γ. As mentioned in Sect. 1, these fluxes contain the emission from all putative emitters within a radius of 6′ around the position of 3EG 2033+4118. They are thus very large upper limits for the non-thermal emission from Cyg OB2 #8A.

We then deduced directly the 3σ upper limit on the non-thermal flux from Cyg OB2 #8A observed during the two NuSTAR observations in the same energy band (20−30 keV). To do so, we applied an absorbed 3T plasma plus power-law model on the four NuSTAR spectra. We constructed the thermal part of the model with the best-fitting parameters determined in Sect. 4.2 for each spectrum (see Table B.1). The non-thermal part of the model is constructed with a power-law whose spectral index is fixed to the different values of Γ reported in Table 3. For each value of Γ, we determined a 3σ upper limit on the flux of the power-law in the 20−30 keV energy band. These upper limits (also reported in Table 3) are more stringent than those determined with INTEGRAL observations. In particular, for the value Γ = 1.5 prevailing for the IC emission, the upper limit on the non-thermal emission is twice better when using NuSTAR observations.

5. Discussion

5.1. Light curve in 2−10 keV

From the best-fitting models presented in Sect. 4.1, we computed the Cyg OB2 #8A light curve showing the evolution of its observed (absorbed) flux over the entire orbit (Fig. 5). This new X-ray light curve is more densely sampled than in the previous studies of De Becker et al. (2006), Blomme et al. (2010), Yoshida et al. (2011), and Cazorla et al. (2014). This allows us to get further insight into the properties of the wind interaction. Our light curve clearly reveals a maximum emission near ϕ ∼ 0.8, whilst the minimum occurs shortly after periastron (De Becker et al. 2006; Cazorla et al. 2014). Our spectral fits indicate that the emission measure of the hardest plasma component follows the same trend as the global flux.

Fig. 5. Cyg OB2 #8A light curve in the 2−10 keV band. The black, red, and blue lines indicate the observed flux and error bars for the Swift, XMM-Newton, and NuSTAR spectra.

Among the CWBs that have been monitored in detail in X-rays, the system that comes closest to Cyg OB2 #8A is probably HD 166 734 (Nazé et al. 2017). This system consists of two O-type supergiants (O7.5 I + O9 I) on a highly eccentric orbit (e = 0.62) with a period of 34.5 days (Mahy et al. 2017). The light curve of this system displays a minimum around the periastron which is attributed to a disruption of the shock (Nazé et al. 2017). Indeed, the residual X-ray emission of HD 166 734 at this phase range is consistent with purely intrinsic emission from the two supergiants. In Cyg OB2 #8A, the situation is different: even at minimum, the system remains overluminous in X-rays by about a factor 20 (log(L_X/L_bol) =  −5.4), indicating that the shock does not get disrupted at periastron.

Finally, we note that the shape of the observed light curve in the 2−10 keV domain is very close to the model cwb4 computed by Pittard (2010). The authors modelled the wind shock region in 3D taking the Coriolis deflection and the radiative cooling into account. They considered two O6V stars with mass-loss rates of 2 × 10⁻⁷ M_⊙ yr⁻¹ and terminal velocities of 2500 km s⁻¹ orbiting each other with a period of 6.1 days and an eccentricity of 0.36. This eccentricity is higher than that measured for Cyg OB2 #8A explaining the larger amplitude of the variations in the cwb4 light curve. The decay of the light curve of cwb4 near the periastron is explained by a change in the cooling regime of the wind shock. Indeed, close to periastron, the shock, which is adiabatic over most parts of the orbit, becomes radiative leading to a sudden cooling of the shocked gas. A similar process may thus explain the shape of the light curve of Cyg OB2 #8a.

Another result that can be noticed is the complexity of the dependence of the observed (absorbed) X-ray flux on orbital separation. Indeed, the flux versus orbital separation plot describes a strong hysteresis (see Fig. 6 and Cazorla et al. 2014), as theoretically predicted by Pittard & Parkin (2010): at the same orbital separation, the emission of Cyg OB2 #8A is different depending on whether the primary is approaching or receding from the secondary.

Fig. 6. Observed X-ray flux of Cyg OB2 #8A in the 2−10 keV energy bands as a function of orbital separation computed with the revised orbital solution (see Table 2). Blue filled circles correspond to XMM-Newton spectra, whilst red asterisks stand for Swift observations binned in orbital phase (0.1 phase bins).

5.2. Evolution of the wind column density

The variations in the X-ray light curves of CWBs can be due to either phase-dependent, line-of-sight wind column densities as the stars move around each other and/or to the variations in the orbital separation with phase. In view of the relatively low orbital inclination inferred in Sect. 3, we expect the strongest influence to come from the changing orbital separation. Nevertheless, we observe some variations in the best-fit wind column density, in which the wind absorption is stronger at phases near periastron than at apastron (see top panel of Fig. 3). This situation most likely reflects the fact that the X-ray emitting part of the wind interaction zone is somewhat compact (i.e. dominated by the regions near the shock apex), is more deeply buried inside the winds at periastron, and is seen through the denser winds at that time. Independent pieces of evidence that the X-ray emission zone is somewhat compact comes from a Suzaku observation analysed by Yoshida et al. (2011). These authors found the flux above 3 keV to vary on timescales of 20 ks. In the 3T best-fit model of Yoshida et al. (2011), this modulation reflects temperature changes (at the 20% level) of the hottest plasma component. Interpreting these variations as hydrodynamical instabilities of the wind interaction zone, Yoshida et al. (2011) estimated a volume of about 2.4 × 10³⁶ cm³ for the hottest plasma, corresponding to a radius (for a spherical plasma volume) of 8.3 × 10¹¹ cm = 11.9 R_⊙ (i.e. 0.1 times the orbital separation at periastron). The actual shape of the wind contact surface is closer to a plane or a cone, and the hottest plasma is likely confined to a layer close to this surface. Therefore, the size of the X-ray emitting region is most probably larger than the stellar radii, but somewhat smaller than the orbital separation, consistent with theoretical expectations.

Whilst no clear trend can be detected in the absorptions fitted to the noisier Swift data, the results obtained for the XMM-Newton spectra indicate a much clearer trend (Fig. 3). The XMM-Newton column density is largest at phases shortly after periastron when the primary is in front. If the shock cone is turned towards the secondary, the denser primary wind would be in front of the X-ray emitting plasma near the shock apex at phases between 0.95 and 0.39 (see bottom panel of Fig. 2), which is consistent with the observations.

We adapted the model of Williams et al. (1990) to evaluate the line-of-sight column density as a function of orbital phase and attempt to constrain the orbital inclination. However, whereas the trend with phase appears well reproduced, the inclination $i=({20}_{-4}^{+6})°$ fitted with this model is significantly lower than our estimate based on the typical masses of the stars, and may be explained by the limitations of this model (see Appendix C).

5.3. Predictions on the inverse Compton emission

Before we can address the issue of the presence or absence of an IC hard X-ray emission, it is important to review the evidence for the presence of non-thermal electrons in the winds of Cyg OB2 #8A. The radio emission of Cyg OB2 #8A is thought to arise from a synchrotron emission process because the emission level is well above that expected from free-free processes in the wind of the stars given their wind density. Moreover, the analysis of the centimetre and millimetre emission shows phase-locked variations clearly revealing that a significant part of this emission arises from the wind interaction region (Blomme et al. 2010, 2017). The thermal free-free emission of the hot high-density gas of the colliding wind region may also give rise to phase-locked variability in excess compared to that produced by the winds of the individual stars. However, the predicted spectral indices in the centimetre wavelength domain for such a purely thermal emission are significantly positive (Pittard 2010) and would thus be difficult to reconcile with the observed spectral index which was found to be close to 0 (Blomme et al. 2010). The radio emission of Cyg OB2 #8A is thus most probably emitted by the non-thermal radiation from relativistic electrons.

In a relatively close CWB such as Cyg OB2 #8A, a significant population of relativistic electrons is required to produce an observable synchrotron radio emission despite the large optical depths of the winds in the radio domain (Blomme et al. 2010). In the innermost parts of the binary these electrons should suffer a substantial loss in energy via IC scattering of the intense stellar UV radiation field. The corresponding emission should take the form of a hard X-ray emission with a power-law spectral energy distribution (Pittard & Dougherty 2006).

The energy of IC scattered photons is given by

$$\begin{array}{c}\hfill h{\nu }_{\mathrm{IC}}\simeq \frac{4}{3}{\gamma }^{2}h{\nu }_{\ast },\end{array}$$(1)

where γ stands for the Lorentz factor of the relativistic electrons. Blomme et al. (2017) showed that electron energies of about 450 MeV, corresponding to γ ≃ 880, are required to explain the 3 mm emission via the synchrotron mechanism. Adopting the spectral type–effective temperature calibration of Martins et al. (2005), we find that the spectra of the two stars in Cyg OB2 #8A peak near 16 eV. Therefore, adopting a conservative value for the Lorentz factor of the relativistic electrons of the order of 100, we expect the IC emission to extend out to energies of at least 100 keV.

To predict the flux of the IC emission from Cyg OB2 #8A we used the theoretical model of the wind-wind collision region and the resulting non-thermal emission created by Pittard et al. (2019). The main characteristics of this model are summarised here. The non-thermal particle distribution is calculated assuming that the shocks are coincident with the contact discontinuity. The immediate post-shock distribution of the primary particles (i.e. directly accelerated at the shocks) is computed at each point along the shocks using the semi-analytic model of Blasi et al. (2005). The diffusive shock acceleration leads to a spectral index of the non-thermal particle distribution that can be energy dependent. The non-thermal particle distribution is then evolved downstream of the shocks by solving the kinetic equation. Energy losses due to IC, synchrotron and relativistic bremsstrahlung emission, and ionic and adiabatic cooling are calculated assuming a black-body distribution of incoming stellar photons and electron to proton number density ratio of 0.01.

We created several models differentiating by the values of the mass-loss rates Ṁ (varying the wind velocities has less impact on the value of the predicted IC emission) and pre-shock magnetic field strength ζ_B which is very uncertain. The last parameter is defined as the ratio between the pre-shock magnetic and the kinetic energy densities. Table 4 lists these values for the four models presented in Fig. 7.

Fig. 7. Predicted anistotropic IC emission from Cyg OB2 #8A. Top panel: model 1 along the orbital phase. Bottom panel: spectral energy distribution of model 1 (purple), 2 (green), 3 (blue), and 4 (orange) at orbital phase 0.028.

The top panel of Fig. 7 shows the predicted anisotropic IC emission along the orbital phase for model 1. The predicted IC emission is maximum prior to periastron and decreases by about 25% near phase 0.2. The variation is similar for the three other models. It is caused by the changing separation of the stars and the orientation of the stars and the wind-wind collision region to the observer.

The bottom panel of Fig. 7 shows the spectral energy distribution of the predicted anisotropic IC emission for the four different models listed in Table 4. Decreasing ζ_B reduces the maximum momentum of the particles, but increases the normalisation of the non-thermal particles at lower energies. Increasing the mass-loss rate of the primary star by a higher factor than those of the secondary star leads to a colliding wind region which is closer to the secondary. This leads to an increased IC cooling and thus a higher total normalisation of the emission.

The largest flux of the IC emission at ϕ = 0.028 between 20 keV and 30 keV is reached by model 4 with a flux of 1.81 × 10⁻¹⁴ erg s⁻¹ cm⁻². It is still about six times lower than the upper limit on the flux obtained with INTEGRAL (F_INTEGRAL) and about 2.5 times lower than the upper limit on the flux obtained with NuSTAR (F_NuSTAR) for Γ = 1.5.

Compared to other works in the literature, the predicted emission at 10 keV from model 1 is about 3 orders of magnitude lower than that for WR 140 (Pittard & Dougherty 2006; Reimer et al. 2006). This large difference is due to the powerful WR wind in the WR 140 system.

Finally, we estimated the kinetic luminosity available in the shock as $L=0.5\dot{M}{\mathit{v}}_{\infty }^2$. For model 1, the kinetic luminosity of the primary wind is L₁ = 1.27 × 10³⁶ erg s⁻¹ while that of the secondary wind is L₂ = 1.27 × 10³⁵ erg s⁻¹. Considering the wind momentum ratio of this model (η = 0.1), Pittard & Dawson (2018) determined the half-opening angles θ of the shock region (as measured from the secondary star) of the primary and secondary shocks to be 73° and 30°, respectively. The fraction of wind that is shocked is f₁ = 0.5(1 − cosθ) = 0.35 for the primary and f₂ = 0.5(1 + cosθ) = 0.93 for the secondary. The total kinetic luminosity is thus L ∼ f₁ L₁ + f₂ L₂ = 5.6 × 10³⁴ erg s⁻¹. The fraction of the wind luminosity that is thermalised is somewhat lower than f₁ and f₂ since the majority of the pre-shock wind is not normal to the shock. Considering the above computed fractions thus leads to an upper limit of the IC to kinetic luminosity ratio. The total kinetic luminosity has to be compared to the full IC emission from model 1 at phase 0.028: 1.2 × 10³¹ erg s⁻¹. The IC emission is thus about three orders of magnitude lower than the available power.

The radiated luminosity from thermal particles is L ∼ (f₁ L₁/χ₁)+(f₂ L₂/χ₂), where $\chi ={\mathit{v}}_8^4{d}_{12}/{\dot{M}}_7$ the cooling parameter for each wind with v₈ the terminal velocity in 10⁸ cm s⁻¹, the distance separating the stars in 10¹² cm, and Ṁ₇ the mass-loss rate in 10⁻⁷ M_⊙ yr⁻¹. This leads to χ₁ = 13 for the primary and χ₂ = 130 for the secondary of model 1. The primary wind thus dominates the total luminosity, which is roughly L ∼ 3.4 × 10³⁴ erg s⁻¹, leading to a flux of f = 1.25 × 10⁻¹⁰ erg s⁻¹ cm⁻². The unabsorbed flux (i.e. corrected from the ISM hydrogen column density) computed from the XMM-Newton spectrum at phase 0.019 in the 0.5−10 keV energy range is about 5 × 10⁻¹¹ erg s⁻¹ cm⁻². The theoretical and observational fluxes are thus in good agreement because the factors f₁ and f₂ overestimate the wind luminosity that is thermalised.

6. Summary and conclusions

We used X-ray data taken with Swift and XMM-Newton to improve the phase coverage of the 2−10 keV light curve of Cyg OB2 #8A. The improved fitting of the X-ray spectra with a 3T plasma model allowed us to achieve a better description of the characteristics of the thermal X-ray emission along the orbit. Our results indicate a light curve with an almost constant level between phases 0.2 and 0.6 followed by an increase until a maximum around phase 0.8 and then a sharp decay with a minimum around phase 0.1. This light curve shape is similar to that of other eccentric systems with comparable orbital periods. However, unlike the case of HD 166 734, the X-ray emission of Cyg OB2 #8A preserves a luminosity well above the level of the intrinsic emission of O-type stars, even at minimum. This indicates that the colliding wind region in Cyg OB2 #8A is not disrupted during the periastron passage.

We also analysed 3.1 Ms of INTEGRAL observations of the EGRET source 3EG 2033+4118 corresponding to the Cygnus OB2 region. This allowed us to refine the 3σ upper limit on the hard X-ray count rate (between 20 and 200 keV) produced by the sources located within less than 6′ from 3EG 2033+4118. We then converted this upper limit on the count rate to an upper limit on the flux of the non-thermal emission by modelling a power-law with different spectral indices Γ. For Γ = 1.5 (expected for the IC emission), the 3σ upper limit on the 20−30 keV emission from the Cygnus OB2 region is 1.1 × 10⁻¹³ erg s⁻¹ cm⁻².

Two NuSTAR observations of Cyg OB2 #8A taken shortly after periastron (phases 0.028 and 0.085) were then investigated to directly search for evidence of an IC hard X-ray emission from the binary thanks to its better angular resolution. Combining these spectra with XMM-Newton data at neighbouring orbital phases, we found no compelling evidence of strong excess above the thermal emission out to energies of 30 keV. Meanwhile, we used the NuSTAR spectra to infer more stringent upper limits on the flux of the putative non-thermal component. For Γ = 1.5, the 3σ upper limit on the 20−30 keV emission from Cyg OB2 #8A is 4.6 × 10⁻¹⁴ erg s⁻¹ cm⁻².

We used the new theoretical model of Pittard et al. (2019) to predict the anisotropic IC emission from the wind shock region. The predicted non-thermal fluxes are lower than our observational upper limits, indicating that another quantum leap in sensitivity is required before the non-thermal X-ray emission of O-type star binaries can possibly be detected.

Acknowledgments

The OHP observing campaigns were supported by the Direction Générale de l’Enseignement Non-Obligatoire et de la Recherche Scientifique of the Fédération Wallonie-Bruxelles. The Liège team thanks the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme (contracts XMaS and HERMES). This work was supported by the Fonds National de la Recherche Scientifique – FNRS, notably under Grant n° T.0192.19. This research has made use of the software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. This research has also made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration and of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). We thank the PIs who obtained the Swift, INTEGRAL and XMM-Newton observations used in this work.

References

Appendix A: log of the X-ray observations of Cyg OB2 #8A

Appendix B: Results of the spectral fitting of the X-ray spectra

Appendix C: Modelisation of the wind column density

We adapted the model of Williams et al. (1990) to evaluate the line-of-sight column density of the stellar winds as a function of orbital phase. For an orbital inclination i ≤ 90°, we have

$$\begin{array}{c}\hfill N_{\mathrm{wind}}(\varphi )=\frac{N_0}{\frac{r}{a}|cos(v+\omega )|sini\theta \sqrt{A}}(arctan\frac{u_2}{\theta }-arctan\frac{u_1}{\theta }),\end{array}$$(C.1)

where N₀ is a scaling parameter depending on the wind density, v is the true anomaly, $u_1=\frac{B}{2\sqrt{A}}$, $u_2=\sqrt{A}tan\psi +\frac{B}{2\sqrt{A}}$, with

$$\begin{array}{cc}& A=1+{cot}^{2}i{tan}^2(v+\omega )\hfill \end{array}$$(C.2)

$$\begin{array}{cc}& B=2|tan(v+\omega )|{cot}^{2}i\hfill \end{array}$$(C.3)

$$\begin{array}{cc}& C=1+{cot}^{2}i\hfill \end{array}$$(C.4)

$$\begin{array}{cc}& {\theta }^2=C-\frac{B^2}{4A}\hfill \end{array}$$(C.5)

and ψ is given by

$$\begin{array}{ccc}\hfill \psi & =\frac{\pi }{2}+(v+\omega )\hfill & \hfill \mathrm{if}(v+\omega )\in [0,\frac{\pi }{2}]\end{array}$$(C.6)

$$\begin{array}{ccc}\hfill \psi & =\frac{3\pi }{2}-(v+\omega )\hfill & \hfill \mathrm{if}(v+\omega )\in [\frac{\pi }{2},\frac{3\pi }{2}]\end{array}$$(C.7)

$$\begin{array}{ccc}\hfill \psi & =(v+\omega )-\frac{3\pi }{2}\hfill & \hfill \mathrm{if}(v+\omega )\in [\frac{3\pi }{2},2\pi ]·\end{array}$$(C.8)

For Cyg OB2 #8A all parameters but i and N₀ are known from the orbital solution. The value N₀ can be determined requiring that the mean value of the computed N_wind equals the mean value of the observed values in the relevant phase interval. We then computed the χ² as a function of orbital inclination by comparing the model to the column density values determined from the XMM-Newton spectra. The lowest χ² is reached for $i=({20}_{-4}^{+6})°$ and the best fit is shown in Fig. C.1. Whereas the trend with phase appears well reproduced, this inclination is significantly lower than our estimate based on the typical masses of the stars and would imply unrealistically high masses of 190 and 150 M_⊙ for the stars.

Fig. C.1. Variations in the column density in 1022 cm−2 obtained in the fits to the XMM-Newton spectra (triangular symbols with error bars) fitted with the model of Williams et al. (1990) for i = 20° (dashed line).

The most likely explanation for this discrepancy is that the model underestimates i because it assumes that the emission arises from a single point at the shock apex. In reality, however, the X-ray emitting region is extended and the wind column density inferred from the observations represents the mean of the columns along the various sightlines towards this extended region. Another limitation of the model is that this does not account for Coriolis deflections of the shock cone, which might also affect the phase-dependence of N_wind.

All Tables

All Figures

	Fig. 1. NuSTAR RGB image (red = 1.6−5.6 keV, green = 5.6−13.6 keV, blue = 13.6−21.6 keV) of the Cyg OB2 region taken on 2018 August 26 (FPMA data only).
In the text

Fig. 2. Revised orbital solution of Cyg OB2 #8A. Top panel: primary and secondary RVs are indicated by the filled and open symbols, respectively. The blue dots correspond to the RVs from De Becker et al. (2004), whereas the red squares indicate the new RVs. The solid and dashed lines illustrate the best-fit RV curve (see Table 2). Bottom panels: relative orbital separation and position angle as a function of orbital phase. A value of 0° for the position angle corresponds to the conjunction with the primary star in front.

In the text

	Fig. 3. Evolution with phase of the absorption and normalisation factors from the fit with three APEC components whose temperatures and n1 component were fixed. The red asterisks, blue circles, and green squares indicate the spectral parameters from the Swift, XMM-Newton, and NuSTAR fits, respectively.
In the text

	Fig. 4. Best-fitting dustscat × TBnew × phabs × (apec+apec+apec) models of the ϕ = 0.028 (top panel) and ϕ = 0.085 (bottom panel) NuSTAR spectra.
In the text

	Fig. 5. Cyg OB2 #8A light curve in the 2−10 keV band. The black, red, and blue lines indicate the observed flux and error bars for the Swift, XMM-Newton, and NuSTAR spectra.
In the text

	Fig. 6. Observed X-ray flux of Cyg OB2 #8A in the 2−10 keV energy bands as a function of orbital separation computed with the revised orbital solution (see Table 2). Blue filled circles correspond to XMM-Newton spectra, whilst red asterisks stand for Swift observations binned in orbital phase (0.1 phase bins).
In the text

	Fig. 7. Predicted anistotropic IC emission from Cyg OB2 #8A. Top panel: model 1 along the orbital phase. Bottom panel: spectral energy distribution of model 1 (purple), 2 (green), 3 (blue), and 4 (orange) at orbital phase 0.028.
In the text

	Fig. C.1. Variations in the column density in 1022 cm−2 obtained in the fits to the XMM-Newton spectra (triangular symbols with error bars) fitted with the model of Williams et al. (1990) for i = 20° (dashed line).
In the text