© ESO, 2012

1. Introduction

Nonthermal radio emission, bright and variable hard X-ray emission, or periodic changes in the profiles of optical emission lines are the typical signatures of colliding winds in massive binary systems. Indeed, only massive stars have strong enough stellar winds for hosting these kinds of energetic phenomena. However, the reverse is not true; i.e., (many) massive binaries exist that show no signature of colliding winds. In fact, nonthermal radio emission is a rare feature (fewer than 40 cases in our Galaxy, see De Becker 2007; Benaglia 2010), and recent X-ray surveys have shown that most binaries are not significantly overluminous at high energies (Nazé 2009). Understanding why it is so will improve our knowledge of both collision mechanisms and stellar wind properties since such emissions heavily depend on them. To this aim, it is important to collect as much information as possible, notably on the rarest cases showing conspicuous collision signatures like Cyg OB2 #9.

As suggested by its name, this star belongs to the Cyg OB2 association, one of the richest OB associations of our Galaxy since it contains about one hundred O-stars (Knödlseder 2000; Negueruela et al. 2008; Wright et al. 2010). As for other members of this association, Cyg OB2 #9 is strongly extinguished, which renders its study quite difficult (E [B − V]  = 2.11, Massey & Thompson 1991). It was nevertheless soon identified as an early-type star (O5If⁺, Morgan et al. 1955; Walborn 1973) and as a nonthermal radio emitter (Abbott et al. 1984). After these first studies, the spectrum of Cyg OB2 #9 was not studied in more detail, so that it was only proven to be a binary in recent years. A hint towards a multiple nature came from the conic shape of the nonthermal emission, a configuration revealed by the Very Large Baseline Array (Dougherty & Pittard 2006) and only possible in wind-wind collisions. Two years later, two important findings were reported. A first, indirect evidence of binarity came from the periodic modulations of the radio emission (van Loo et al. 2008), occurring with a timescale of 2.355 yr. Soon after that, direct evidence of binarity, i.e., clear doubling of stellar lines, was unveiled by a dedicated spectroscopic campaign (Nazé et al. 2008). The continued monitoring of the system then led to computation of a preliminary orbital solution (Nazé et al. 2010): Cyg OB2 #9 appeared as a long-period, very eccentric (e = 0.7) binary composed of two nearly equal-mass objects. In parallel, the analysis of several XMM-Newton observations of this object revealed no strong overluminosity or peculiar hardness in the X-ray spectrum, but did show some unusual long-term variability (by about 10%, Nazé et al. 2010). Clearly, additional data are needed to pinpoint the properties of the system: improving the orbital solution, enlarging the orbital coverage of the X-ray data (the available ones being all taken far from periastron), and extending the knowledge of the radio emission. Only then can a full modeling be undertaken, to better understand massive stars and their winds.

Since the discovery of the binarity of Cyg OB2 #9, two periastron passages have occurred. The first one took place in the beginning of 2009, when the system was in conjunction with the Sun and therefore unobservable. The second one was expected to occur in June − July 2011, i.e., at the best time for observing the Cygnus constellation. We thus organized a large observing campaign of Cyg OB2 #9 to observe all aspects of this event: radio monitoring of the nonthermal emission using the Expanded Very Large Array, interferometric monitoring to disentangle the astrometric orbit of the two stars using CHARA, spectroscopic monitoring in the optical to improve the orbital solution, and X-ray monitoring to detect and understand the high-energy emission from the wind-wind collision. The results of this multiwavelength campaign will be published as a series of papers, which will be concluded by a modeling of the system. With all these observations in hand, Cyg OB2 #9 will become one of the best known nonthermal radio emitters, with detailed constraints on the wind characteristics and on the collision properties. Cyg OB2 #9 will thus represent a true stepping stone in our understanding of wind-wind collisions and their associated emissions.

This first paper provides the results obtained through optical spectroscopy and X-ray monitoring. Section 2 presents the observations, Sect. 3 analyzes the spectroscopic data with the aim to improve the orbital solution, Sect. 4 discusses the variations observed in X-rays, and Sect. 5 summarizes our findings and concludes this paper.

2. Observations

2.1. The optical domain

2.1.1. Sophie

Our monitoring of Cyg OB2 #9 at the Haute-Provence Observatory (France) began in 2007. It mainly used the high-resolution echelle spectrograph Sophie installed on the 1.93 m telescope. It yields an homogeneous dataset covering a broad wavelength domain with the highest quality available.

The Sophie spectrograph provides coverage of the 3900 − 6900 Å domain with a resolving power R = λ/Δλ of 35 000. In total, 75 echelle spectra were obtained with this instrument. The exposure time was typically 40 min, and the derived signal-to-noise ratios vary from 30 to 80. Whenever the exposures were close (i.e., they were taken within a month and no large radial-velocity [RV] variations were detected), these data were combined, providing in the end 31 independent spectra spanning a period of four years. As for all Sophie spectra, these data were reduced by the automatic Sophie pipeline¹: it locates the orders, performs an optimal extraction with cosmic-ray rejection, calibrates the spectra, and finally reconnects the orders after barycentric correction and correction for the blaze function (derived from a Tungsten lamp exposure – for massive stars, it provides only a rough, first-order normalization). These data were smoothed further using a sliding window of 15 px size (with 1 px  =  0.01 Å), and corrected by hand for any remaining cosmic ray hits. The journal of observations is provided in Table 1.

Owing to the low signal-to-noise and the high extinction of the star, it is difficult to perform a clean and coherent normalization. To ensure homogeneity, we first performed the best normalization possible, using low-order polynomials, on the data from June 2010: these data are not only among the best quality data, but they were also taken at a time when lines are most blended, making the definition of continuum windows easier. We then used this normalized set as a reference for the other data. It enables us to remove stellar lines (through a division of the spectra by the reference), allowing better modeling with low-order polynomials of the continuum oscillations due to the poor blaze correction and poor order merging.

In a final step, we tried to correct for some annoying nonstellar features. To this aim, we have fitted with Gaussians the diffuse interstellar bands (DIBs) in the 5412 Å and 5800 Å regions (see Fig. 1) on the June 2010 data (i.e., when lines are most blended). In addition, a DIB near 5812 Å was unveiled when the stellar lines were deblended – this particular DIB apparently has a similar rest wavelength to the C iv line. It was thus fitted on the deblended spectra and added to the set of DIBs already defined. No good fitting could be found for the DIB at 5797 Å, so that we kept it in the spectra. The fitted interstellar features were then removed from the data, leading to cleaner stellar lines. We have also used the table from Hinkle et al. (2000) to correct the Hα, He i λ 5876 Å, and C iii λ 5696 Å regions for pollution by the telluric lines.

Fig. 1 The spectrum from June 2010 (black thin line) and the DIB correction (red thick line).

2.1.2. WIRO

Our team also monitored Cyg OB2 #9 at the Wyoming Infrared Observatory (WIRO) 2.3 m telescope many times between 2011 May 31 and 2011 July 25 (see Table 1), which provides a more intense coverage of the periastron event, as well as an independent check of the stellar behavior at that time. All new datasets were reduced using standard IRAF reduction routines as outlined in Kiminki et al. (2007). Sixteen observations were obtained using the WIRO-longslit spectrograph with the 1800 line mm^-1 grating in first order over the spectral range 5250–6750 Å. Exposure times varied from 720 s to 5400 s depending on weather conditions; the resulting signal-to-noise was about 200 near 5876 Å. Copper-Argon lamp exposures were taken after each observation to wavelength-calibrate the spectra to an RMS of 0.03 Å. The typical resolution was 1.5 Å FWHM across the chip. Spectra were Doppler corrected to the heliocentric frame and then shifted further by small amounts (typically  < 6 km s^-1) so that Na i interstellar absorption lines in each exposure fell at the mean Na i velocity averaged over the ensemble of 16 measurements. The shift also ensures that interstellar features, e.g. the Na i λ 5889 line, display similar velocities in WIRO and Sophie spectra, enabling us to directly compare them. Finally, continuum normalization was performed using low-degree polynomials.

2.2. X-rays

2.2.1. Swift

The Swift facility observed Cyg OB2 #9 during 5 ks at five carefully selected dates in 2011 (PI Nazé, see Table 1), to follow the star during ten months around the date of periastron. The data were then reduced following the recommendations of the UK Swift center². First, the pipeline (FTOOLS task xrtpipeline) was applied. Events were then extracted in a circular region of 47′′ radius for the source and in a nearby circle of 119′′ radius for the background. The chosen region is as devoid of X-ray sources as possible (see Rauw 2011). The adequate response matrix file (RMF, swxpc0to12s6_20010101v013.rmf) is provided by the Swift team, but the ancillary response file (ARF) was specifically calculated for Cyg OB2 #9 by the task xrtmkarf with the inclusion of an exposure map so that bad columns can be taken into account. Count rates measured in these observations vary from 0.04 to 0.1 cts s^-1, so there is no pile-up.

Fig. 2 Spectrum observed in June 2010 (after DIB and telluric correction). The labels indicate the stellar lines (other features are interstellar). The He ii λ 6560 Å line next to Hα is strongly blended with a DIB.

2.2.2. XMM-Newton

XMM-Newton had observed the Cygnus OB2 region six times in previous years, mostly for studying the X-ray emission of Cyg OB2 #8A (see Table 1). Since Cyg OB2 #9 was in the field-of-view, we analyzed its X-ray emission and found some unexplained variations (Nazé et al. 2010). A full monitoring of Cyg OB2 #9 with XMM-Newton in 2011 was not possible, but an observation was scheduled near the predicted time of the periastron passage (25 ks, ObsID = 0677980601, PI Nazé, see also Table 1).

This new dataset was reduced in the standard way using SAS v. 10.0.0. After the pipeline processing (SAS tasks emproc, epproc), we filtered the data to only keep the best events, as recommended by the SAS team³. A flare in the X-ray background affects the middle of the observation. It was eliminated, so that the effective exposure times are reduced to 24 ks for MOS and 15 ks for pn.

Source events were extracted in a circular region of 35′′ radius around the target, and the background region was defined as a nearby circle of the same size. An annulus could not be used, as numerous X-ray sources can be found in the region (see Rauw 2011). The count rate measured for this new observation is 0.6 cts s^-1 for MOS and 1.8 cts s^-1 for pn, which is five times higher than in previous observations. Despite this brightening, the source is not affected by pile-up, as checked by a run of the task epatplot.

3. Cyg OB2 #9, as seen in the optical

3.1. The spectra

The high extinction of Cyg OB2 #9 has two direct consequences on its optical spectrum. First, the blue part of the spectrum (<4400 Å) is very noisy, hence not really usable. Second, most lines recorded in the spectra actually are interstellar features! Only a few stellar lines can be distinguished and used for further study (see Fig. 2): Balmer lines (Hα, Hβ), He i λ 4471, 5876 Å, He ii λ 4542, 4686, 5412, 6683 Å, N iii λ 4634, 4641 Å (in emission), C iii λ 4650, 5696, 6723, 6729 Å (in emission), C iv λ 5801, 5812 Å. Some of these stellar lines are, however, contaminated by interstellar features, such as C iv λ 5812 Å.

When one compares the spectra taken at different phases (see Fig. 3), it is obvious that the C iii λ 5696 Å emission line follows the motion of one component of the system. The He ii λ 4686 Å emission line instead displays the opposite motion, that of the other component, though with a shift of about 90 km s^-1. This shift implies that the line most probably does not arise from the photosphere, but somewhere close to it, in the wind. The exact position may be influenced by the presence of a wind-wind collision in the system (see below). These two lines are the sole ones “associated” to a single component of the system. Other lines show two components in the 0.9–1.1 phase interval (Fig. 3). Another problem is obvious when one looks at the spectra where both components are clearly separated: the lines of both objects have similar strengths in He ii, N iii, and C iv, which does not help the derivation of the RVs and their attribution to one of the components.

Fig. 3 Comparison of the spectral appearance of Cyg OB2 #9 in June 2010 (φ = 0.571, solid line) and in the beginning of July 2011 (φ = 0.005, thick dotted line). The spectra have been corrected for some DIBs (see text) and the telluric lines (the latter correction being imperfect around 5696 Å).

Fig. 4 RVs determined by simple fits (top) and by finding the best shift of fixed Gaussians (bottom left). Small cyan symbols correspond to WIRO data, other symbols are the same in all panels, with one colored symbol for each line. The RV measured for a narrow interstellar feature is superimposed on the He ii λ 4686 Å subpanel using asterisks, to show the quality of the wavelength calibration. The bottom right panel shows the RVs determined by disentangling. The best-fit orbital solution (Table 3) found from disentangling is superimposed on each panel.

Fig. 5 Disentangling results. Each panel shows, from top to bottom, the spectra of the primary star, the secondary star, and the interstellar features (with continuum shifted to +1) in four different spectral domains. The normalization is imperfect, as is often the case for disentangling results when the orbital RV shifts are modest.

Fig. 6 Top: variation with phase of the RVs (solid line = primary, dashed line = secondary), separation, and position angle (defined to be zero when the primary is in front and along the line-of-sight). Bottom: the orbit of the primary star around the center of mass (see Table 3 for physical parameters). The position of periastron is shown with a big dot, of phases 0.1 to 0.9 with open squares, of phases 0.05 and 0.95 with crosses, and of phase 0.99 with a star. The line of nodes corresponds to the x-axis and the line-of-sight to the y-axis, with an Earthly observer towards the bottom. At periastron, a (arbitrary) sketch of the collision zone is shown as a pink dashed line, with the position of the secondary star shown as an open blue circle.

3.2. RV measurements and orbital solution

As already explained in Nazé et al. (2010), measuring the RVs in Cyg OB2 #9 is difficult because the broad stellar lines are faint, sometimes even contaminated by interstellar lines, and also blended for most of the (highly eccentric) orbit. In a first step, we measured the RVs by fitting one (or two, if the line displayed two components or appeared asymmetrical, which happens within Δφ = 0.1 of periastron) Gaussian(s) to the stellar lines. Results of these simple fits are shown at the top of Fig. 4. It is obvious from that figure that deblending the two components is not easy, since they can only be measured at maximum separation. The shift of He ii λ 4686 Å with respect to the true stellar motion is also obvious, as is the attribution to the primary star of C iii λ 5696 Å⁴ and the main component of He i lines. WIRO data agree well with the Sophie observations, when taking their larger uncertainty (about 20–30 km s^-1) due to a lower resolution into account.

To improve over these simple fits, we fixed the component’s shapes to the Gaussians fitted in the most separated cases, and we cross-correlated them with the blended lines until the best χ² was found, as we did in Nazé et al. (2010). The new RV curves are better, as can be seen in Fig. 4. However, a trend must be noted: the χ² mapping has a tendency to separate the components more at phases of maximum blending. This may be due to slight changes in line shapes (due to noise and imperfect normalization), but allowing the line shapes to vary slightly produces erratic results of lower quality.

Finally, we used these preliminary RVs as input for a disentangling program based on the method of González & Levato (2006). In this program, the spectra of each component are iteratively determined, first on the basis of the approximate input RVs, then by using improved RVs found from cross-correlating the components’ spectra with the observed spectra (for details see Mahy et al. 2010). In our case, the disentangling tries to fit the spectra of three components: the fixed interstellar features, and the primary and secondary stars (whose RVs vary from one observation to the next). The continuum was arbitrarily fixed at 0, 0.5, and 0.5 for these three components (i.e., the true light ratios are not taken into account in this method). The cross-correlation was performed using masks, containing the He i λ 4471 Å, He ii λ 4542 Å, He ii λ 5412 Å, and He i λ 5876 Å lines for the primary and the He ii λ 4542 Å, He ii λ 5412 Å, and He i λ 5876 Å lines for the secondary. The spectra where both components are clearly visible (i.e., maximum separation) received a weight 50 times greater than in the blended cases. This process yields the RV curves seen in the last panel of Fig. 4 (for values, see Table 2), and the disentangled spectra are shown in Fig. 5.

The RVs derived from the disentangling can be used to find the orbital solution. To this aim, we used the Liège Orbital Solution Package (LOSP)⁵, which is based on the SB1 method of Wolfe et al. (1967) generalized to SB2 cases as described, e.g., in Rauw et al. (2000a). The period was fixed to 858.4 d, as found from an extensive radio dataset covering three decades (Blomme et al., in prep.)⁶, and the derived parameters are shown in Table 3 and Fig. 6. It agrees well with the preliminary solutions from Nazé et al. (2010).

3.3. Spectral features

While photospheric He ii absorption lines display similar strength in both components (Fig. 3), the He i lines appear much stronger in the primary component (i.e., the one that appears blueshifted at periastron). In fact, He i λ 4471 Å is so weak in the other component that it is difficult to distinguish it. Comparing the relative strength of He ii λ 4542 Å and He i λ 4471 Å (Conti 1973), we found that they correspond to an O5 − 5.5 type and O3 − 4 type for the primary and secondary, respectively (Walborn et al. 2010, reported the blended spectrum as O4.5Ifc). That one component is hotter than the other one may explain why C iii λ 5696 is best seen in the cooler component while C iv lines have similar strengths in both stars. Since the He ii lines display similar intensities, it also suggests that the two stars should have different luminosity classes, the coldest object being also bigger. This conclusion is supported by the fact that identical luminosity classes would also be incompatible with the coldest star being the most massive component (see previous section).

To investigate this issue further, we turned to other spectral features. Both components display N iii and the neighboring C iii lines in emission. The latter lines are not variable and are much stronger in the primary, which agrees well with the finding that these features reach their maximum strength in O5 stars (Walborn et al. 2010). In our disentangled spectra, the He ii λ 4686 Å line appears in emission in both components, but the width of the feature in the primary spectrum suggests that this line may be contaminated. Indeed, the disentangling program is not fit to study such features since it considers the line shapes to be constant throughout the orbit and to arise solely in the stellar atmospheres. The shift observed for the He ii λ 4686 Å line (see above) clearly indicates that it is not the case: even if the stars might truly display some He ii λ 4686 Å in emission, most of it should be formed in the wind, close to the secondary star and with a potential influence of the wind-wind collision. With N iii + C iii lines (and possibly He ii λ 4686 Å) in strong emission, the primary spectrum thus appears quite compatible with a supergiant class of type fc.

Determining the luminosity class of the secondary star is more difficult. We checked the spectra for signs of N iv λ 4058 Å (Walborn et al. 2002): this line appears as a very weak emission, much weaker than N iii λ 4634, 4641 Å, but the spectra are so noisy in the blue range that it is difficult to perform disentangling, hence to say whether the line belongs to one component or to both. In either case, this suggests that the hotter component (the secondary) is O3.5V or O4III-V and not a supergiant, whereas the presence in its spectrum of N iii and He ii λ 4686 Å lines in emission suggests the contrary!

Fig. 7 Evolution of the Hα line in the Sophie (top) and WIRO (bottom) datasets. Phases are quoted near each spectrum and, in the top panel, the distance between spectra is proportional to the phase difference. Note the two absorptions that disappear before periastron, the line profile then being almost triangular, and the increased absorption and emission components at periastron.

We used the typical stellar parameters of giant O3.5 and O5.5 stars (Martins et al. 2005, temperature and gravity are similar for giants and supergiants in this reference) as input for the atmosphere modeling code CMFGEN. This enabled us to determine the equivalent widths of the He i λ 4471 Å and He ii λ 4542, 5412 Å lines on the synthetic spectra, which we compared to the values found by fitting Gaussians to the observed spectra (at maximum velocity separation) or to the disentangled spectra. This led to nearly equal luminosity ratios, while the orbital solution further suggests equal masses. Therefore, our favored spectral types are O5–5.5I for the primary and O3–4III for the secondary.

Our monitoring of Cyg OB2 #9 allows us to search for spectral signatures of a wind-wind collision. To this aim, one generally turns to Balmer and He ii λ 4686 Å emissions. As already mentioned, the He ii λ 4686 Å line has a peculiar behavior: its shape is relatively constant, but its RVs are shifted relative to the stars (Fig. 4). The recorded velocities instead suggest the line arises close to the secondary star (but not close to its photosphere, because of the velocity shift), and a contamination by the collision zone is possible but difficult to assess.

In contrast, the broad and intense Hα emission line displays no obvious shift with orbital phase but its profile clearly varies though the changes are not dramatic (Fig. 7). We first describe these variations. In this context, one has to keep in mind that only a single periastron passage could be observed in details, that of 2011, so that caution must be applied when generalizing the observations with φ ~ 0 to the usual behavior of the star. At phases φ > 0.1, the line displays an overall triangular shape with a basis covering about 2000 km s^-1 and a peak at  ~200 km s^-1. Superimposed on this emission are two relatively narrow absorptions with velocities of about  −175 and 45 km s^-1 (i.e., not stellar velocities). It should be noted that the RVs of these two absorptions do not seem to vary much, but neither do the stellar RVs at these phases. As time goes by, the emission peaks decrease slightly in strength, though, as do the absorptions (the blueshifted absorption disappearing only after the redshifted one). Before periastron (φ = 0.85−0.98), the line thus appears as a nearly perfect triangle. Suddenly, shortly before periastron (at about φ = 0.99 to the best of our knowledge – our data do not enable us to further pinpoint the exact time of the event), the emission peak becomes more intense and is shifted to the red. It now appears at about 300 km s^-1. In parallel, a strong absorption appears at  ~−200 km s^-1, and is also much broader than before, since it is approximately 60% wider than the previous blueshifted absorptions. A substructure of this absorption may be related to the photospheric absorption by the primary star, though this component is not obviously detected at other phases. In the WIRO data, the emission appears strongest at φ = 0.999, whereas the absorption is most intense at φ = 0.998, but this cannot be confirmed by Sophie data owing to the lack of coverage of these phases. After periastron, the emission peak decreases in strength and moves to less blueshifted velocities, while the absorption becomes less broad, and a second absorption gradually appears, coming back to the start of the cycle.

A very broad Hα emission is quite compatible with a line arising not at the photosphere, but over a wide area – maybe the wind-wind collision in our case. If the wind-wind collision zone is dense enough at periastron, it could absorb some light from the stellar photospheres, in a kind of atmospheric eclipse. In such a case, the broad absorption appearing at periastron may be related to absorption of the primary’s light by the colliding wind zone (Fig. 6). Indeed, the collision zone may well be highly distorted at periastron by Coriolis effects (e.g. Parkin & Pittard 2008), and one of its “wings” could then obscure our view of the primary, as seen from Earth. In such a configuration, the emission would come from the other zones of the collision, which are mostly redshifted and do not appear in front of one star, as seen from the observer’s point of view. The pair of narrower absorptions seen at other phases is more difficult to explain, since the density of the collision zone decreases strongly with increasing separation, but they could be explained by a similar effect if there is enough dense cold gas downstream. Their constant velocity, though, is not surprising since the stars themselves do not show much RV change at these phases. The disappearance of all structures at φ = 0.85 − 0.98 is, however, puzzling. It could be related to a physical change in the collision zone, or be simply due to viewing angle effects.

Fig. 8 Variations in the spectral parameters (left) and observed flux (middle) with phase, for the best-fit model with the first temperature and the absorptions fixed. The solid and open symbols correspond to XMM-Newton and Swift data, respectively. The right panel shows, for the same model, the flux ratio as a function of the inverse stellar separation ratio Dref/D for the XMM-Newton and Swift data. The reference data corresponds to the first observation of each dataset, i.e., data from Oct. 04 for XMM-Newton and data from Jan. 11 for Swift. The black circles corresponds to the total, 0.5–10.0 keV flux, the green squares to the soft, 0.5–2.0 keV flux, and the blue triangles to the hard, 2.0–10.0 keV flux.

4. Cyg OB2 #9, as seen in X-rays

4.1. The spectra

First, we fitted the new XMM-Newton spectra and the old ones with optically-thin thermal plasma models of the form wabs ∗ phabs ∗ (apec + apec). The fitting was performed within Xspec v12.4.0ad (which uses APEC 1.3.1): the first absorption was fixed to the interstellar one, 1.15 × 10²² cm^-2, and the second one allows for additional, in-situ absorption by the winds⁷; two thermal emission components were necessary to obtain a good fit. Because these XMM-Newton data display the higher signal-to-noise ratios, this fitting enables us to see whether temperatures and absorption change with phase. The values found for the additional absorption and the cooler temperature in the different datasets were found to remain within one sigma of each other, so that we decided to fix them. Swift data have fewer counts, hence are noisier. A fully free 2T fitting results in large error bars and sometimes erratic variations. It was thus not possible to check the absorption and cooler temperature constancy on these data. Therefore, all Swift spectra were fitted with the cooler temperature and additional absorption fixed. This seems a reasonable assumption since (1) the cooler component has a temperature typical of the intrinsic emission of hot, massive stars (see e.g. Nazé 2009, and references therein); and (2) the additional absorption is probably associated to the stellar winds, and their observed densities should not change much with phase since the two stellar components are similar (unlike what happens in WR+O systems). The contribution from a wind-wind collision should occur over a wider temperature range, and thus be the sole one responsible for the hotter component, and it should also be more variable than intrinsic emissions, especially in highly eccentric systems (see below). Table 4 yields the results of these fits, and Fig. 8 illustrates the variations of the parameters with orbital phase.

We also tried to fix the hotter temperature further, but the χ² was significantly larger for the XMM-Newton data (Table 4). However, the flux changes were similar for both XMM-Newton and Swift, so that the conclusions reported in the next section are not affected by this choice.

Two datasets were taken close to periastron, 5.5 days apart, one by XMM-Newton and one by Swift. They show a difference in flux of about 20%. The cross-calibration between Swift and XMM-Newton data is believed to be better than 20% (A. Beardmore, priv. comm., see also Tsujimoto et al. 2011). To check our results, we looked at the other sources in the Swift field-of-view. Cyg OB2 #5 and 12 are known to vary, but unfortunately, their variations are not fully explained and predictable (Linder et al. 2009; Rauw 2011). The behavior of the colliding-wind binary Cyg OB2 #8A is, however, well known (De Becker et al. 2006; Blomme et al. 2010), and the Swift data are here fully consistent with the XMM-Newton and ASCA data, when taking the errors due to a lower number of counts in the Swift data into account. The rapid brightening of Cyg OB2 #9 near periastron, on timescales of days, may thus well be real.

4.2. Archival data

The Cyg OB2 association has been observed by all major X-ray facilities since its main members were detected in X-rays in 1978 by Einstein. To see whether our conclusions on Cyg OB2 #9 are valid, we also considered these old datasets.

Einstein IPC count rates of Cyg OB2 have been reported by Waldron et al. (1998, see also the Einstein : they amount to 0.0257 ± 0.0018 cts s^-1 in mid-December 1978 (φ = 0.16), 0.0237 ± 0.0042 cts s^-1 on 1979 Nov. 22 (φ = 0.56), and 0.0322 ± 0.0014 cts s^-1 on 1980 Jun. 9 (φ = 0.79). However, Chlebowski et al. (1989,HEASARC “ipcostars” catalog) quote an IPC count rate of 0.0406 ± 0.0010 cts s^-1, which is not compatible with the previous values though they come from the same data. It is thus difficult to know exactly what happened 30 years ago.

ROSAT observed Cyg OB2 in April 1991 and 1993 with the PSPC (ObsID = 200109 and 900314, respectively) and in November 1994 and April 1995 with the HRI (ObsID = 201845 parts 1 and 2, respectively). Spectra of the source (reduced to a single bin for HRI) were extracted using the same regions as for Swift observations. In these regions, the ROSAT background-corrected count rates of Cyg OB2 #9 amount to (3.96 ± 0.36) × 10^-2 cts s^-1 in PSPC/1991, to (3.83 ± 0.15) × 10^-2 cts s^-1 in PSPC/1993, to (2.3 ± 0.2) × 10^-2 cts s^-1 in HRI/1994, and to (1.47 ± 0.06) × 10^-2 cts s^-1 in HRI/1995. The PSPC values are compatible with the PSPC count rates reported by Waldron et al. (1998)⁸. We furthermore fitted the PSPC spectra using the same model as above, with temperatures and absorption fixed (the poor signal-to-noise and low sensitivity at high energies did not permit us to free more parameters). Fitting results are reported in Table 4.

Both PSPC datasets were taken far from periastron, at phases 0.41 (in 1991) and 0.28 (in 1993), so that we expected few changes in recorded intensity. Indeed, the count rates, fluxes, and spectral parameters of the two observations agree with each other within the errors. The 1993 data was also taken at a similar phase as the sixth XMM-Newton observation, and there is also fair agreement between them. HRI observations occurred just before and just after a periastron passage at phases 0.95 for the 1994 observation (i.e., similar in phase to the third Swift observation) and 0.11 for the 1995 observation (i.e., similar in phase to the fifth Swift observation). In fact, the recorded count rates are higher in the 1994 data, which is closest to the periastron time.

Using the HRI response matrix and the same model as in the last XMM-Newton observation (i.e., same temperatures, absorptions, and normalization ratio), we found an observed flux in the 0.5–2. keV band of 8 and 5 × 10^-13 erg cm^-2 s^-1 for the 1994 and 1995 observations, respectively. These fluxes agree well with the ones observed at similar phases with Swift. Overall, the ROSAT data thus appear compatible with the variations seen in XMM-Newton and Swift observations.

ASCA observed Cyg OB2 in April 1993, at phase 0.27, nearly at the same time as the second ROSAT-PSPC observation. We extracted SIS spectra for both SIS-0 and SIS-1, and for both bright and bright2 modes (i.e., 4 spectra in total, implying that we merged data taken with different bit rates). For the source, we used the same region as in Swift and ROSAT, whereas the background region was chosen as a nearby circle of the same radius and on the same CCD. The source region is quite small for the ASCA resolution, but Cyg OB2 #9 appears close to a CCD edge. Another problem in this observation is the high background: first, straylight from the nearby extremely bright source Cyg X-3 is contaminating the field-of-view; second, Cyg OB2 #9 appears on the PSF wings of the bright Cyg OB2 #8A. Indeed, Cyg OB2 #9 is not detectable by eye in single datasets, and it appears only as a faint object in the SIS global image. The resulting spectra are thus noisy, but we nevertheless tried to fit them with the same model as above (with temperatures and absorption fixed, see Table 4). The observed flux, in the 0.8−4. keV range amounts to (8.4 ± 0.5) × 10^-13 erg cm^-2 s^-1, slightly smaller but within the errors of the value of (1.1 ± 0.4) × 10^-12 erg cm^-2 s^-1 reported by Kitamoto & Mukai (1996). Compared to the sixth XMM-Newton dataset (obtained at a similar phase) or the contemporary ROSAT observation, we again find a smaller flux in our ASCA fitting. However, considering the background and resolution problems, the agreement is not that bad.

In summary, the archival data agree well with the XMM-Newton and Swift results, providing further support to the presence of phase-locked variations in Cyg OB2 #9 since these different datasets were taken several periods apart.

4.3. A wind-wind collision in Cyg OB2 #9

The presence of bright and hard X-ray emission (with kT > 1.5 keV) is quite unusual in O-type stars. It may be linked to a magnetically-confined wind in a single star or to a wind-wind collision in a binary system. In our case, even though the nonthermal radio emission implies the presence of a magnetic field, the latter explanation seems to be the correct one. Indeed, Cyg OB2 #9 displays long-term variations, apparently compatible with the timescale of the orbital period.

Theoretical considerations show that there are two main sources of phase-locked variability amongst colliding wind binaries (Stevens et al. 1992; see also Güdel & Nazé 2009, for a review of recent observations). The first one is the absorption. Since the collision is seen alternatively through the wind of each component, changes in absorption are expected. This is especially relevant for systems with two very different winds, e.g. WR+O: two spectacular cases in this respect are γ² Vel (Willis et al. 1995) and WR22 (Gosset et al. 2009). The similarity of the two components in Cyg OB2 #9 makes such variations less prominent; indeed, no significant change in absorption is detected in the data.

The second source of phase-locked variability is the changing separation in eccentric binaries (Stevens et al. 1992). For adiabatic cases, the separation between the stars is high, and the X-ray luminosity is expected to follow a 1/D relation (i.e., it should be maximum at periastron) with a minimal impact from the velocity changes because the winds are close to their maximal speeds at all phases. For radiative cases, the X-ray luminosity and plasma temperatures should be minimum at periastron, as the winds then collide at lower speeds. Adiabatic systems should be those with long periods, like Cyg OB2 #9, but not many long-period binaries have been monitored in X-rays, and results were mixed. γ² Vel and WR22 did not show any 1/D variation, despite the changing separation of their components (Rauw et al. 2000b; Gosset et al. 2009). Cyg OB2 #8A and WR140 deviate from expectations at periastron, since the collisions then become radiative (De Becker et al. 2006; Pollock et al. 2002; De Becker et al. 2012). Hints of the influence of a 1/D relation may have been found in HD 93205 (see Fig. 2 in Antokhin et al. 2003), though the results need confirmation, in view of the huge noise, and in HD 93403, though the intrinsic emission of the stars blurred the results, yielding a luminosity change lower than that expected from a simple 1/D relation (Rauw et al. 2002). Among evolved systems, WR25 seems to display a rather “clean” 1/D variation in its luminosity (Gosset 2007; Güdel & Nazé 2009). For most of its orbit, it is also the case of η Car (Okazaki et al. 2008), though the system shows large cycle-to-cycle variations in flux, and there are strong indications that the wind collision region collapse onto the companion star at periastron (Parkin et al. 2011b, and references therein). Up to now, no comparable example among O+O binaries has been found, so Cyg OB2 #9 may well become the first one.

Indeed, while there seems to be no change in absorption, at least for our sampling, the normalization factors and fluxes clearly change with time. They show a sharp increase near periastron (see middle panel of Fig. 8). In fact, a strong correlation appears when the fluxes are compared to the inverse of the stellar separation. The correlation is not perfect, but we know that there are instrumental differences (up to 20% in flux, see above) and larger errors for the less sensitive Swift data. To reduce the impact of these instrumental effects, we performed a comparison within each dataset (see last panel of Fig. 8 and Table 5): within the errors, the fluxes vary as 1/D, as expected from an adiabatic collision in an eccentric system (Stevens et al. 1992). In fact, the linear correlation coefficient between the distance and flux ratios listed in Table 5 is 0.97 for Swift+XMM-Newton data and 0.999 for the more precise XMM-Newton data alone. In this context, it may also be noted that the increase by a factor of 1.6 in the ROSAT-HRI count rates is well matched by the decrease in separation between orbital phase φ = 0.11 and φ = 0.95 by the same factor. The only slight departure from the 1/D relation (Fig. 8) is the hard X-ray flux appearing slightly below the expected value: this may be due to the collision becoming slightly radiative because of the increased density at periastron. This would also explain the presence of cool plasma responsible of Hα absorption at these phases in the optical spectra. However, this does not affect the bulk of the X-ray emission, since the 1/D correlation is closely followed overall, with no evidence of blurring by an additional, constant intrinsic emission from the stars. The absence of dramatic changes in Hα also correlates well with an adiabatic nature for the wind-wind collision.

On the other hand, the second temperature of the spectral fitting is clearly lower near periastron in the XMM-Newton data (see left panel of Fig. 8). This may be explained by a shock occurring at a smaller separation. Indeed, a shorter distance between the two stars may imply that winds collide before they could accelerate to full speed, resulting in a cooler shocked plasma. Moreover, at small separations, radiative braking/inhibition (Stevens & Pollock 1994; Gayley et al. 1997) may take place, slowing down the winds further before they collide. Such a cooling of the hard component is predicted by models of wind-wind collisions in eccentric binaries: Pittard & Parkin (2010, see their Fig. 22) find a change by a factor of 3 in the highest temperature, whereas Parkin & Gosset (2011) find an order of magnitude variation when the shock plunges deep into the wind acceleration region. In Cyg OB2 #9, the temperature only changes by 25%, but the separation is much larger than in the case modeled by Pittard & Parkin (2010) and the wind-wind collision is closer to balance than in the case modeled by Parkin & Gosset (2011). Therefore, we expect the influence of wind acceleration in Cyg OB2 #9 to be less pronounced than found by these authors, weakening the cooling effects.

Fig. 9 Evolution with phase of the collision parameters (position of balance point, cooling parameter χ, pre-shock velocities, and average plasma temperature, kT = 0.6 × v(in      1000   km   s-1)2). Thick red lines and thin black lines correspond to the cases with and without radiative inhibition, respectively.

To put these observing results in perspective, we did a first modeling of the wind-wind collision in Cyg OB2 #9. We here use the analytical models of Stevens et al. (1992; see also a detailed example of the calculation in Parkin et al. 2011a). From the ram pressure balance, they allow us to find the position of the collision, the cone’s opening angle, plasma temperatures (from the pre-shock wind velocities⁹), and the cooling parameter (to determine the nature of the collision, adiabatic if  >1, radiative otherwise). As input, such models require wind and stellar parameters: stellar radii and separation, wind mass-loss rates, and velocities. The orbital solution has just been found, yielding the separation at each phase. However, not much information is available on the individual wind parameters, only approximate spectral types have been derived. We therefore made the following assumptions: inclination of the system of about 62 degrees (so that the true masses are compatible with the predictions of  Martins et al. 2005), typical stellar parameters of O5-5.5I and O3-4III stars taken from Martins et al. (2005), and mass-loss rates derived from the recipe of Vink et al. (2000). This yields mass-loss rates of 5 × 10^-6 M_⊙ yr^-1 and 1.2 × 10^-5 M_⊙ yr^-1 and terminal velocities (=2.6 × escape velocities) of 2540 km s^-1 and 2440 km s^-1 for the primary and secondary, respectively. We also assume the usual β-law with β = 1 for the wind velocity evolution. At maximum separation, this results in a wind momentum ratio of 0.43 (or 2.3 in favor of the secondary), which gives an opening angle of the collision zone of about 150 degrees, compatible with the radio map (see Fig. 6 of Dougherty & Pittard 2006,note that the angle estimated from that figure is approximate, because sensivity or wind opacity may affect its determination). Figure 9 shows the evolution with phase of several parameters of the collision zone. The position of the balance point (at 0.4D from the primary) does not change with phase, and the cooling parameter χ (Stevens et al. 1992) is always  ≫ 1, implying that the shocks should still be relatively adiabatic even at periastron (which nicely conforms with the X-ray observations). The agreement is not perfect, however, since the shocked plasma temperatures appear too high (about 3.6 keV at apastron and 3.1 keV at periastron vs. observed values of 2.5 keV and 1.9 keV), and the relative change in temperatures is slightly too low (13% vs. 25% observed). We therefore include the effect of radiative inhibition in the model: rather than using a simple β velocity law, we here assume the winds to be radiatively driven following the Castor et al. (1975) approximation and including the influence of the opposing star radiation field (e.g. Stevens & Pollock 1994; Gayley et al. 1997). Such a model improves the agreement with observations: the winds collide at lower speeds, even at apastron, resulting in lower plasma temperatures (about 2.7 keV at apastron and 2.0 keV at periastron), and the decrease in plasma temperature is now 25% between apastron and periastron (Fig. 9). In this case, the shock still remains mostly adiabatic throughout the orbit, as observed. Wind clumping may play a further role in decreasing plasma temperatures (Pittard 2007, and references therein). The next step will be to make full 3D hydrodynamical simulations of the wind-wind collision (cf. Pittard 2009; Parkin & Gosset 2011), to calculate synthetic spectra and lightcurves. Their comparison with observations will allow determination of the density, temperature, and velocity of the hot plasma within the wind-wind collision region, helping to constrain the stellar mass-loss rates and wind speeds in this system. Constraints on the magnetic field within the collision region will also be derived. However, this detailed modeling requires all data to be at hand (not only optical and X-ray observations – see future papers of the series); it is also clearly beyond the scope of this paper.

5. Summary and conclusions

In mid-2011 the first observable periastron passage of Cyg OB2 #9 occurred since the binarity discovery. To observe this event under all aspects, we organized a multiwavelength campaign with a simultaneous coverage at radio, infrared, optical, and X-ray energies. Here we report on the results of the optical and X-ray monitorings.

The optical spectroscopic monitoring obtained with the Sophie spectrograph enabled us to refine the orbital solution, confirming preliminary results. A spectral disentangling was also performed, revealing the two stars to be of O3−4III and O5−5.5I types for the secondary and primary, respectively. While the He ii λ 4686 Å line follows the motion of the hottest component (though with a shift), the Hα line appears extremely broad but unshifted. The line profile changes with phase, with both enhanced absorption and stronger emission at periastron.

X-ray observations by XMM-Newton and Swift have been obtained as Cyg OB2 #9 reached the 2011 periastron. The high-energy monitoring revealed a rather constant X-ray spectrum, except at periastron where the emission becomes slightly softer. The second thermal component then has a temperature reduced by about 25% at periastron, which is reproduced well by collision models including radiative inhibition. The monitoring also unveiled the luminosity variations of the system, which follow the 1/D changes predicted for adiabatic collisions remarkably closely. The only (slight) departure from that relation is a very small decrease in the hard-band flux, which could indicate that the collision has become very slightly radiative at periastron. Overall, Cyg OB2 #9 seems the best (and the first!) example of an adiabatic O+O colliding-wind system.

Online material

Acknowledgments

The authors thank Tabetha Boyajian, Michaël De Becker, and Gregor Rauw for useful comments. We also thank the Swift PI, Dr. Neil Gehrels, and the XMM-Newton project scientist, Norbert Schartel, for having made the X-ray monitoring possible. Y.N. acknowledges useful discussion with Kim Page on Swift data reduction and calibration. The Liège team acknowledges support from the European Community’s Seventh Framework Program (FP7/ 2007 − 2013) under grant agreement number RG226604 (OPTICON), the Fonds National de la Recherche Scientifique (Belgium), the Communauté Française de Belgique, the PRODEX XMM and Integral contracts, and the “Action de Recherche Concertée” (CFWB-Académie Wallonie Europe). J.M.P. would like to thank the Royal Society for funding a University Research Fellowship. We thank the OHP staff and colleagues who took the service mode observations. We thank Michael Alexander, Garrett Long, Mike Lundquist, Jessie Runnoe, Rachel Smullen, Carlos Vargas, and Earl Wood for their assistance in obtaining the data at WIRO. We thank WIRO staff James Weger and Jerry Bucher for their work that enabled this science. ADS and CDS were used for preparing this document.

References

All Tables

All Figures

	Fig. 1 The spectrum from June 2010 (black thin line) and the DIB correction (red thick line).
In the text

	Fig. 2 Spectrum observed in June 2010 (after DIB and telluric correction). The labels indicate the stellar lines (other features are interstellar). The He ii λ 6560 Å line next to Hα is strongly blended with a DIB.
In the text

	Fig. 3 Comparison of the spectral appearance of Cyg OB2 #9 in June 2010 (φ = 0.571, solid line) and in the beginning of July 2011 (φ = 0.005, thick dotted line). The spectra have been corrected for some DIBs (see text) and the telluric lines (the latter correction being imperfect around 5696 Å).
In the text

Fig. 4 RVs determined by simple fits (top) and by finding the best shift of fixed Gaussians (bottom left). Small cyan symbols correspond to WIRO data, other symbols are the same in all panels, with one colored symbol for each line. The RV measured for a narrow interstellar feature is superimposed on the He ii λ 4686 Å subpanel using asterisks, to show the quality of the wavelength calibration. The bottom right panel shows the RVs determined by disentangling. The best-fit orbital solution (Table 3) found from disentangling is superimposed on each panel.

In the text

	Fig. 5 Disentangling results. Each panel shows, from top to bottom, the spectra of the primary star, the secondary star, and the interstellar features (with continuum shifted to +1) in four different spectral domains. The normalization is imperfect, as is often the case for disentangling results when the orbital RV shifts are modest.
In the text

Fig. 6 Top: variation with phase of the RVs (solid line = primary, dashed line = secondary), separation, and position angle (defined to be zero when the primary is in front and along the line-of-sight). Bottom: the orbit of the primary star around the center of mass (see Table 3 for physical parameters). The position of periastron is shown with a big dot, of phases 0.1 to 0.9 with open squares, of phases 0.05 and 0.95 with crosses, and of phase 0.99 with a star. The line of nodes corresponds to the x-axis and the line-of-sight to the y-axis, with an Earthly observer towards the bottom. At periastron, a (arbitrary) sketch of the collision zone is shown as a pink dashed line, with the position of the secondary star shown as an open blue circle.

In the text

	Fig. 7 Evolution of the Hα line in the Sophie (top) and WIRO (bottom) datasets. Phases are quoted near each spectrum and, in the top panel, the distance between spectra is proportional to the phase difference. Note the two absorptions that disappear before periastron, the line profile then being almost triangular, and the increased absorption and emission components at periastron.
In the text

Fig. 8 Variations in the spectral parameters (left) and observed flux (middle) with phase, for the best-fit model with the first temperature and the absorptions fixed. The solid and open symbols correspond to XMM-Newton and Swift data, respectively. The right panel shows, for the same model, the flux ratio as a function of the inverse stellar separation ratio Dref/D for the XMM-Newton and Swift data. The reference data corresponds to the first observation of each dataset, i.e., data from Oct. 04 for XMM-Newton and data from Jan. 11 for Swift. The black circles corresponds to the total, 0.5–10.0 keV flux, the green squares to the soft, 0.5–2.0 keV flux, and the blue triangles to the hard, 2.0–10.0 keV flux.

In the text

	Fig. 9 Evolution with phase of the collision parameters (position of balance point, cooling parameter χ, pre-shock velocities, and average plasma temperature, kT = 0.6 × v(in      1000   km   s-1)2). Thick red lines and thin black lines correspond to the cases with and without radiative inhibition, respectively.
In the text