Free Access
Issue
A&A
Volume 608, December 2017
Article Number A69
Number of page(s) 9
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/201731403
Published online 08 December 2017

© ESO, 2017

1. Introduction

Massive stars have strong stellar winds. When two such stars form a binary, their stellar winds collide, leading to a number of interesting effects. The collision between the two winds heats up the material to such temperatures that it emits detectable X-ray radiation, as originally predicted by Prilutskii & Usov (1976) and recently reviewed by Rauw & Nazé (2016). In the shocks associated with the colliding-wind region (CWR), a fraction of the electrons are accelerated to high, non-thermal energies by the Fermi mechanism (Bell 1978; Eichler & Usov 1993). These electrons spiral in the magnetic field, thereby emitting synchrotron emission, which we can detect as non-thermal radio emission (Bieging et al. 1989).

In an eccentric binary, the strength of the collision (more specifically, the ram pressure) will vary with orbital phase, leading to phase-locked variations of the intrinsic emission coming from the CWR. At radio wavelengths, there is furthermore the effect of the strong free-free absorption in the stellar wind material (Wright & Barlow 1975). The changing positions of both stellar winds and the CWR lead to further phase-locked variations in the observed radio fluxes. This has been observed in a number of O+O and Wolf-Rayet+O colliding-wind binaries (e.g., White & Becker 1995; Blomme et al. 2005, 2007; Blomme & Volpi 2014). In exceptional cases, the CWR has been resolved by high-resolution radio observations (Benaglia et al. 2015, and references therein).

Table 1

Observing log of the NOEMA 3 mm observations of Cyg OB2 #8A.

Colliding winds also emit non-thermal X-rays and gamma-rays (Pollock 1987; Benaglia & Romero 2003; Pittard et al. 2006; Reimer et al. 2006; De Becker 2007; Reitberger et al. 2014), making these objects very relevant to high-energy physics (Pittard & Dougherty 2006). The regime of physical parameters (density, magnetic field, ambient radiation field, etc.) is quite different from that of other high-energy environments, and colliding-wind binaries can therefore provide important tests of our understanding of the physical processes responsible.

Colliding-wind systems can also be highly relevant for the mass-loss rate determinations in single stars. This has become a major problem in massive star research in recent years, due to the fact that winds are clumped and porous. Taking into account clumping decreases the mass loss rates by a factor 3–10 compared to what was previously thought (e.g., Puls et al. 2008). This, in turn, has major consequences for stellar and galactic evolution. Colliding-wind systems can provide an independent determination of the effect of clumping and porosity (e.g., Pittard 2007): by modelling the colliding winds, and predicting the various observational indicators, we can constrain the mass-loss rate and/or porosity in the wind.

Colliding-wind binaries have so far been relatively unexplored at millimetre wavelengths. This is mainly due to the expectation that free-free emission will start to dominate the fluxes at shorter wavelength, and that therefore the non-thermal synchrotron emission will be difficult or impossible to detect.

Furthermore, the CWR will also increase the free-free emission, as shown by Stevens (1995) and Montes et al. (2011), using simplified models for a radiative shock; and by Dougherty et al. (2003), Pittard et al. (2006) and Pittard (2010), using hydrodynamical calculations for both adiabatic and radiative shocks. The density increase leads to a considerable amount of additional free-free emission. Montes et al. (2015) tried to look for this, by comparing centimetre and millimetre fluxes for a sample of 17 Wolf-Rayet stars. For three of these, a hint from the colliding-wind region contribution was found, but their observations did not allow to clearly identify the physical mechanism.

Another approach is to look for the variability of the millimetre fluxes. An eccentric binary will show phase-locked variability of the millimetre fluxes for the same reason as for the radio fluxes. The changing distance between the two components leads to variation in the strength of the collision and hence the amount of additional free-free emission. Added to that is the variation in the free-free absorption due to the changing position of the stellar winds and the CWR.

Table 2

Measured 3 mm fluxes of Cyg OB2 #8A.

Cyg OB2 #8A (=Schulte 8A = MT 465 = BD+40 4227A) is a well-known massive colliding-wind binary, and therefore an appropriate candidate to look for variability at millimetre wavelengths. It is a member of Cyg OB2, which is an association containing a large number of massive stars (Knödlseder 2000; Wright et al. 2015). The association harbours a number of massive colliding-wind binaries that have been studied in detail (Van Loo et al. 2008; Blomme et al. 2010, 2013; Kennedy et al. 2010).

Cyg OB2 #8A was discovered to be a binary by De Becker et al. (2004). It has a period P = 21.908 ± 0.040 d, an eccentricity e = 0.24 ± 0.04, and it consists of an O6If primary and an O5.5III(f) secondary (De Becker et al. 2006). It has all the attributes of a colliding-wind binary. Its radio spectrum shows a negative spectral index (Bieging et al. 1989) and the radio flux variations are locked to the orbital phase (Blomme et al. 2010). The X-ray emission is overluminous compared to that of single O-type stars and its variability is also phase-locked (De Becker et al. 2006; Blomme et al. 2010). Among the O-type colliding-wind binaries, it is the system where the parameters of the stars, their stellar wind and their orbit are best known (De Becker 2007).

We observed Cyg OB2 #8A with the NOrthern Extended Millimeter Array (NOEMA) interferometer of the Institut de Radioastronomie Millimétrique (IRAM1), to see if we could detect the expected flux variations at millimetre wavelengths. We present the observations in Sect. 2 and the data reduction in Sect. 3. The resulting millimetre light curve is given in Sect. 4 and we model it in Sect. 5. In Sect. 6, we present the results for Cyg OB2 #8B, which is the visual companion to our main target. The conclusions are given in Sect. 7.

2. Observations

Cyg OB2 #8A was observed with the IRAM NOEMA interferometer during two different periods: end of 2014 (project: S14AW; PI: D.M. Fenech) and mid-2016 (project S16AU; PI: R. Blomme). The 2014 observations consist of three runs (a fourth run, on 2014-11-29, contains only calibrator data and no on-target observations, so it will not be considered further). The 2016 observations have 11 runs. The detailed observing log is given in Table 1.

Observations were collected at 3 mm (110 GHz; IRAM receiver no. 1 in Upper Side Band), in single-pointing mode. Dual polarization was used with a bandwidth of 2 × 3.6 GHz. An observing run consists of alternating between one or two phase calibrators and Cyg OB2 #8A, ensuring that each Cyg OB2 #8A observation is preceded and succeeded by a phase calibrator observation. In addition, for each run, the flux calibrator MWC 349 and the bandpass calibrator 3C 454.3 were observed (for the 2016-06-22 observation the bandpass calibrator 3C 273 was also observed).

During the 2014 runs, the two phase calibrators used were 2200+420 (at an angular distance Δ = 16.̊7) and 2005+403 (Δ = 4.̊9); in 2016, only 2013+370 (Δ = 5.̊4) was used. A single Cyg OB2 #8A observation consists of 30 scans (33 scans in 2016) of 45 s each. For each phase calibrator observation, three 45 s scans were done. The signals from the antennas were correlated with the WideX correlator. During the 2014 observations, the antennas were in the 6Cq configuration (with E10 missing for the 2014-12-01 run). For the 2016 observations, many changes in antenna positions were made; we therefore list the detailed configurations in Table 1.

3. Data reduction

We used the GILDAS2 software for the data reduction, specifically the packages CLIC (for calibration) and MAPPING (for the flux determination).

We split the observing runs, so that we can separately calibrate each series of 30 (or 33 scans) on Cyg OB2 #8A. The phase calibrator(s) immediately preceding and following the Cyg OB2 #8A observations were selected, as well as the bandpass and flux calibrators that are nearest in time. Based on the pipeline reduction that had run at IRAM, we flagged some data which were clearly discrepant, or were on one side of a cable phase jump, or had high system temperatures. The 2014 observations were done in a stable configuration, and therefore needed no baseline corrections. Antenna movement was frequent however during the 2016 observations, and for all these runs, baseline corrections (provided by IRAM staff) were applied.

Default automatic flagging was applied, which flags problems with antenna shadowing, data that do not have surrounding calibrator observations, and timing errors. Next, atmospheric phase corrections were applied. The receiver bandpass was then calibrated on the bandpass calibrator 3C 454.3, using the default degree of the polynomials. For the 2016-06-22 observations, the 3C 454.3 data were not usable, and the bandpass was calibrated on 3C 273. Phase calibration consists of fitting a per-antenna linear function to the phases of the phase calibrator(s). Fluxes were then calibrated on the flux calibrator MWC 349, which has a flux of 1.23 Jy. For the amplitude calibration, a per-antenna linear function was fitted to the calibrators.

As the final part of the calibration, the Data Quality Assessment procedure was run. This flags data with more than 40° root-mean-square (rms) of phase loss, and more than 20% of amplitude loss. The fraction of points that were flagged in this way are indicated in Table 2. They are a good indicator of the quality of the data, in the sense that a higher fraction of flagged points indicates worse quality. Pointing, focus and tracking errors which are too high are also flagged, but this did not occur for any of our observations. Finally, the calibrated data of Cyg OB2 #8A were written out. At this stage, all Cyg OB2 #8A data taken during a single observing run were merged again.

Next, we made a deep image by combining all the data (Fig. 1), to see if there are any other sources in the field besides Cyg OB2 #8A. We used a pixel size of 0.4′′× 0.4′′, which is about one-quarter (linear dimension) of the synthesized beam size. We took a grid size of 256 × 256 pixels, which covers somewhat more than twice the primary beam and we applied natural weighting. We then cleaned the image, stopping when we reached 0.025 mJy, which is half of the rms noise level in the image. Figure 1 shows only the inner 20′′× 20′′ part of this image.

The field is dominated by Cyg OB2 #8A, which is – as expected – a point source. A much weaker source is seen at the bottom edge of Fig. 1. We identify it as Cyg OB2 #8B, and discuss it further in Sect. 6. We checked the full image for other possible sources, but found none.

thumbnail Fig. 1

Deep 3 mm image of all Cyg OB2 #8A data combined. The colour scale values are in Jy. The contour levels are at −0.1 mJy, + 0.1 mJy and then go up in steps of 0.1 mJy to 1.6 mJy. The negative contour is indicated by the dashed, white line. The crosses indicate the positions of Cyg OB2 #8A and #8B (from SIMBAD). The synthesized beam (shown in the lower left corner) is 2.93″ × 1.76″ with a position angle of 57.̊1.

Open with DEXTER

As the data are dominated by the flux from Cyg OB2 #8A, we determined its flux by fitting a point source model to the observed visibilities for each observing run. This is a better approach than measuring the flux on images made from these visibilities, as the cleaning procedure might introduce artefacts. The fluxes we obtained are listed in Table 2.

A number of variant reductions were also tried: switching off the atmospheric phase correction, using weights in the phase or amplitude calibration, or averaging the polarizations in these calibrations, or flagging frequency parasites (i.e., sharp peaks in frequency caused by the instrument). These all gave values within the listed error bars. As the Cyg OB2 #8A flux is relatively high, we also tried to include it in the calibration process, but this did not improve results.

Table 3

Overview of models used in Sect. 5.

thumbnail Fig. 2

Observed 3 mm fluxes of Cyg OB2 #8A, plotted as a function of orbital phase in the 21.908-day binary period. The blue symbols show the 2014 data, the red ones the 2016 data. To better show the behaviour of the light curve, the phase range is extended by 0.2 on each side. Open symbols indicate duplications in that extended range. The grey data show the 6 cm observations from Blomme et al. (2010). Note that the 3 mm flux scale (left) is different from the 6 cm one (right).

Open with DEXTER

4. Millimetre light curve

The Cyg OB2 #8A fluxes listed in Table 2 clearly indicate variability. To check that this is not some artefact of the data reduction, we also list the fluxes of the phase calibrators. The phase calibrators are intrinsically variable, and their measured fluxes are clearly not correlated with the Cyg OB2 #8A ones. We can therefore exclude data reduction artefacts. The absolute flux calibration at 3 mm is good to about 10% (J. M. Winters, pers. comm.), which corresponds to about ±0.2 mJy. The flux variations are therefore also not due to changes in the levels of the absolute flux calibration.

We next plot the observed fluxes in the orbital phase diagram (Fig. 2). The fluxes are clearly correlated with the orbital phase, indicating that they are at least partly formed in the CWR. Around phase 0.7, there is a large range in the fluxes, which seems to violate the assumption that, in a colliding-wind binary, the fluxes repeat nearly perfectly from one orbit to another. We note, however, that the two extreme fluxes are of lower quality: the amplitude quality indicator in Table 2 shows that 11 and 13% of the data for these two observations have been removed by the Data Quality Assessment procedure, because they show more than 20% amplitude loss. Note that the standard error bar does not include this systematic effect. On the basis of this quality information, we consider the flux variations around phase 0.7 to be not significant.

Figure 2 also shows the 6 cm radio fluxes from Blomme et al. (2010). Note that different flux scales are used for the 3 mm and the 6 cm data. The shape of the 3 mm light curve follows that of the 6 cm one very well, but differs from it in two aspects. First, the millimetre data are above the median (1.66 mJy) about 30% of the time, while the 6 cm data are above their median about 50% of the time. This is due to the non-uniform distribution of the phases for which we have flux determinations. There are a large number of observations around the phase of maximum flux (0.7–0.8), but none around the suspected minimum (phase 0.2–0.3). Additional data around the minimum would lower the median and lead to a more equal distribution of the phase range above and below the median.

Secondly, the 3 mm light curve is phase-shifted with respect to the 6 cm one. One could consider that the uncertainty in the period is responsible for this shift. The 21.908 ± 0.040 d period was derived by De Becker et al. (2004) on the basis of spectroscopic data from the years 2000–2003. However, both X-ray data covering the years 1991–2004 and radio data covering the years 1980–2005 are consistent with this period, suggesting that the 0.040 d error bar is too conservative (De Becker et al. 2006; Blomme et al. 2010). We therefore consider it unlikely that the phase shift is due to uncertainty in the period.

Table 4

Cyg OB2 #8A parameters used in modelling.

Alternatively, we note that a phase shift in the radio light curves at different wavelengths is quite common for colliding-wind binaries. The most detailed example is WR 140, as shown by White & Becker (1995). Specifically for Cyg OB2 #8A, a hint of such a phase shift is present in Fig. 1 of Blomme et al. (2010). It can be seen when comparing the position of minimum for the 6 cm and 3.6 cm radio observations.

5. Modelling

In the following sections a number of different models will be applied to the observed data. Table 3 summarizes these models, allowing the reader to keep track of which model is used in which section.

5.1. Stellar wind contribution

We first estimate how much the two stellar winds contribute to the observed millimetre fluxes. As a first approximation, we apply the Wright & Barlow (1975) formalism. The expected flux Sν (in mJy) at frequency ν (in Hz) of a spherically symmetric, steady-state wind, with a mass-loss rate (in M yr-1) flowing out at a constant velocity ν (in km s-1) is: (1)where μ is the mean atomic weight (1.27 for solar composition), D is the distance (in kpc), is the mean squared ion charge and the Gaunt factor gff at frequency ν and electron temperature Te is given by (Leitherer & Robert 1991): (2)We apply these equations using the stellar wind parameters which were derived by De Becker et al. (2006), and which are also listed in Table 4. For Te, we take half the effective temperature of the star. The distance to the Cyg OB2 association is not well known. We take D = 1.4 kpc, following Morford et al. (2016), who give a detailed overview of the various distance determinations. Using these values, we find a 3 mm flux of 1.62 mJy for the primary, and 0.75 mJy for the secondary.

A basic assumption of the Wright & Barlow (1975) model is that the wind is flowing out at a constant velocity. Because of the wavelength-squared dependence of the free-free absorption, the millimetre fluxes are formed closer to the star than the centimetre fluxes. In this formation region, the velocity will not necessarily have reached its terminal value, and this will lead to a different value for the flux (see, e.g., Pittard 2010; Daley-Yates et al. 2016). We therefore also made a numerical integration of the specific intensities and the flux in a spherically symmetric wind, assuming a β = 0.8 velocity law. The effect for our specific set of stellar wind parameters is minor however: we now obtain 1.63 mJy for the primary, and 0.76 mJy for the secondary. Larger differences occur at shorter wavelengths where the wind is probed still deeper (Pittard 2010).

The summed flux of both stars gives 2.39 mJy, which is comparable to the highest flux value observed. This may seem surprising, but it should be realized that we are no longer dealing with spherically symmetric winds. In the overlap region between the two winds, the wind of the other component is “missing” in the sense that the material has been compressed into the CWR. Models presented in the following sections will take this effect into account.

For use in the following sections, we also derive the radius at which the optical depth of 1 is reached. Using the same assumptions as Wright & Barlow (1975), it can be shown that: (3)where (4)Filling in the numbers we find Rτ = 1 = 220 R and 140 R for the primary and secondary, respectively.

5.2. CWR synchrotron emission

We first estimate if synchrotron emission can still be a contributing factor at 3 mm. A relativistic electron with energy E will emit synchrotron radiation over a range of frequencies. The frequency νmax (in MHz) at which the maximum emission occurs is given by (Ginzburg & Syrovatskii 1965, their Eq. (2.23)): (5)where B is the local magnetic field (in G), θ is the angle between the velocity vector of the electron and the magnetic field vector, me is the mass of the electron, and c is the speed of light. This equation neglects the Razin effect (Ginzburg & Syrovatskii 1965, their Sect. 4), which decreases the flux at longer wavelengths.

Taking an order-of-magnitude value of 0.1 G for B, and neglecting the sinθ factor, we find that 3 mm (110 GHz) emission is mainly due to electrons with an energy of E ≈ 450 MeV. To know if such electrons exist in the CWR, we need to balance the Fermi acceleration against inverse Compton cooling, which is the main cooling mechanism. Inverse Compton cooling is due to the interaction of the relativistic electrons with the high-energy photons from both stars. Van Loo et al. (2005) derive an expression for the highest energy that electrons can attain under these conditions: (6)where e is the charge of the electron, χs is the shock strength, Δu is the velocity jump, σT the Thomson electron scattering cross section, and L the luminosity of the star. For an estimate, we take values which are roughly applicable to Cyg OB2 #8A at periastron: r = 100 R, log L/L = 5.7, B = 0.1 G, and for the shock: χs = 4 and Δu = 1500 km s-1). This gives E ≈ 1000 MeV. There are therefore electrons in the CWR that have sufficient energy to emit synchrotron radiation at 3 mm.

To investigate the synchrotron emission in more detail, we used the model from Blomme et al. (2010). This model follows in detail the momentum distribution of the relativistic electrons as they move away from the shock, taking into account inverse Compton and adiabatic cooling. A simplifying assumption of the model is that it does not solve the hydrodynamical equations, but it uses an analytic approach to determine the position of the shocks. From the momentum distribution, the synchrotron emissivity is calculated. This is then mapped into a three-dimensional grid (assuming rotational symmetry around the line connecting the two stars). In this grid, the radiative transfer equation is solved, taking into account the synchrotron emission as well as free-free absorption and emission due to the stellar winds. No free-free contribution from the CWR is included. The model is applied at a number of orbital phases, and the theoretical flux is determined for each of these phases. Further details of the model are given in Blomme et al. (2010).

thumbnail Fig. 3

Theoretical fluxes compared to the observed 3 mm fluxes of Cyg OB2 #8A, plotted as a function of orbital phase in the 21.908-day binary period. The black symbols with error bars show the observed data. Open symbols indicate duplications in the extended phase range (as in Fig. 2). The dark blue curve is the synchrotron model, the red curve the adiabatic thermal emission model, the green curve the radiative thermal emission model, and the light blue curve the combined synchrotron and adiabatic thermal emission model.

Open with DEXTER

We applied the Blomme et al. (2010) model to Cyg OB2 #8A, using the stellar, stellar wind, and orbital parameters from De Becker et al. (2006, see also Table 4. The simulation cube consists of 10243 cells, has a size of 4000 R, and is centred on the primary. The stellar wind material is assumed to be at 19 000 K (i.e., about half the effective temperature of the stars).

The resulting 3 mm fluxes are shown in Fig. 3 (dark blue curve). We note that the general flux level is too high, as is the flux range. One should take into consideration that no parameters were adjusted in the modelling, and some quantities are not well known (e.g., the fraction of the shock energy that is transferred to the relativistic electrons). Because the predicted flux levels coincide within an order of magnitude with the observations, the curve can be viewed as indicative of the physical processes in play but it clearly does not tell the whole story. Note that when Blomme et al. (2010) applied this same model to the 6 cm observations of Cyg OB2 #8A, they obtained theoretical fluxes that are systematically too low. The maximum of the theoretical light curve is approximately as broad as the observed one. The main discrepancy is that this theoretical maximum is phase shifted with respect to the observations.

Blomme et al. (2010) also presented a model with stellar wind parameters for the primary that are different from the De Becker et al. (2006) ones ( = 1.0 × 10-6M yr-1, ν = 2500 km s-1). While this helped to explain some features of the 6 cm observations, using these values for the 3 mm observations only shifts the maximum even further away from the observed one.

A possible explanation for the phase shift of the maximum is our neglect of the orbital motion on the shape of the CWR. On a large scale, the contact discontinuity and its associated shocks take on the shape of a spiral, as shown by Parkin & Pittard (2008) using a simplified dynamical model. This spiral shape is confirmed by the full hydrodynamical models of e.g., η Car (Parkin et al. 2011; Madura et al. 2013). For the purposes of understanding the effect on the 3 mm fluxes of Cyg OB2 #8A, however, one needs to focus on a smaller scale, which extends only somewhat beyond the Rτ = 1 distances of 220 and 140 R. At this scale the orbital motion pushes the leading edge of the CWR closer to the star with the weaker wind (see, e.g., Fig. 13 of Parkin et al. 2011). The position of the leading CWR edge in a model that does include orbital motion is therefore approximately the same as the position at an earlier phase in a model that does not include orbital motion. A first-order correction of our theoretical light curve is therefore a shift of the curve towards the right in the orbital phase diagram. This is indeed the correct direction to improve the agreement of the theoretical and observed phases of maximum. Of course, we cannot estimate the size of the shift from this qualitative argument. For this we would need detailed hydrodynamical calculations, which are outside the scope of this paper.

5.3. CWR thermal emission – adiabatic model

In the previous section, we have only considered the contribution of the synchrotron emission (and the stellar wind free-free emission and absorption) to the observed fluxes. It is expected that the highly compressed material in the CWR will also contribute free-free emission and absorption, which could be detectable in the observations (Stevens 1995; Pittard 2010; Montes et al. 2011).

To explore the effect of the thermal emission, we use another radiative transfer model, developed by Blomme & Volpi (2014). This model uses an adaptive grid scheme allowing us to efficiently determine the emergent intensities and flux from the stellar winds and the CWR. The CWR is assumed to have a shape that is centred around a cone. The position and opening angle of this cone are determined analytically (Eichler & Usov 1993, their Eqs. (1) and (3)). The CWR has a half flaring angle (α), so it is thin at the apex and expands as we move away from the line connecting the two stars (see Fig. 2 of Blomme & Volpi). The size of the CWR scales with the distance between the two stars. We assume that the density inside the CWR is four times higher than what the wind density would be at that distance (this assumption will be revised in another model we discuss – see Sect. 5.4). We refer to the Blomme & Volpi paper for further details of the model.

For each phase of the orbit, we put the two stars in a three-dimensional grid, and assign the densities and temperatures of the stellar winds to each cell. The temperature of the CWR material is a free parameter; the wind material is assumed to be at half of the effective temperatures of the stars. The grid is 16 000 R on each side, and we start with 2563 cells. We then solve the radiative transfer equation taking into account the free-free emission and absorption, and refining the cells where needed (adaptive grid scheme).

thumbnail Fig. 4

CWR of the thin radiative shock model projected on to the sky, for phases where the flux is high (left column) and where the flux is low (right column). The orbital phase is listed in the top left corner of each panel. For plotting purposes, the CWR has been cut off at a distance of 600 R.

Open with DEXTER

We explored a range of values for the temperature of the CWR, its flaring angle, and its size. The best fit we found (as judged by eye) is shown in Fig. 3 (red curve). Its parameters are a CWR temperature of 1.2 × 105 K, a half-flaring angle of 30° and a size of 1.8 times the separation between the two stars. These parameters were chosen from a large set of experiments, because they get the average flux level and the range of fluxes approximately correct. Similarly as for the synchrotron model, however, the maximum is phase shifted with respect to the observations. It is also broader than observed. The use of the alternative Blomme et al. (2010) model (with a lower mass-loss rate for the primary) does not solve the problem of the phase-shifted maximum. Similarly as in Sect. 5.2, the phase shift can probably be attributed to our neglect of orbital motion.

To find out how much the CWR is contributing to the total flux, we also ran a model without a CWR, by setting the flaring angle to zero. We then find that the 3 mm flux varies between 1.32 and 1.39 mJy, which is just marginally below the minimum of the red curve in Fig. 3. The two stellar winds therefore provide the “baseline” flux, and the CWR is responsible for the variability on top of that baseline flux. The CWR region is undetectable around phase 0.0–0.2 because the primary with the strongest wind is in front of it (conjunction is at phase 0.08). The free-free absorption in the stellar wind blocks most of the flux from the CWR. Around phase 0.4, the two stars are well separated on the sky, and that part of the CWR region which is beyond the free-free absorption of the stellar winds becomes detectable. At later phases, the secondary with the weaker wind moves in front (conjunction at phase 0.71), which reduces the detectable flux again, but not so much as when the primary was in front.

5.4. CWR thermal emission – radiative model

The model in the previous section limits the density contrast in the CWR to a factor of four. This is appropriate for an adiabatic shock. But, at least in part of the orbit the shock will be radiative, and the density contrast could be higher. This can be seen from the values of the cooling parameter χ, as defined by Stevens et al. (1992). De Becker et al. (2006) derived χ values for Cyg OB2 #8A which are around one (0.19–1.65, depending on component and orbital phase; see their Table 11). This indicates that the shocks are radiative at least part of the time (near periastron, where χ ≪ 1), and can become adiabatic near apastron.

We therefore adapted the Blomme & Volpi (2014) code to include the semi-analytical model for radiative shocks in colliding stellar winds developed by Montes et al. (2011). This model assumes that the post-shock cooling is very efficient, leading to a thin CWR with a high density contrast, and a temperature that is equal to the wind temperature. Montes et al. determine the position of the contact discontinuity semi-analytically, and from the conservation equations they derive the surface density and emission measure of the CWR material. To avoid substantial changes to our code, we do not directly use the emission measure, but we convert it into a density in a region around the contact discontinuity. As a standard value, we take 20 R for the width of that region.

Applying this model we get the results presented by the green curve in Fig. 3. The theoretical light curve shows a clear double-wave pattern, which is not seen in the data. The double wave is explained in Fig. 4, where we plot the CWR projected against the sky for a few different phases. When we see the CWR turned towards us (phases 0.125 and 0.65), there is a larger emitting area than when we see it on its side (phases 0.325 and 0.925). This explains the minima and maxima of the theoretical curve.

An important element in this explanation is that the CWR still has a significant flux contribution beyond the free-free absorption region of the stellar winds. E.g., at phase 0.325, the separation of the two stars on the sky is about 160 R, which should be compared to the Rτ = 1 radii of 220 R and 140 R (Table 4). Although these Rτ = 1 radii are based on an assumed spherically symmetric wind, they still indicate that the part of the CWR closest its apex will not be detectable. In contrast to this, for the adiabatic model from Sect. 5.3 (red curve on Fig. 3) we have limited the size of the CWR to 1.8 times the separation between the two stars, so the CWR is only just beyond the Rτ = 1 radius. Hence, it does not show the double wave effect. Of course, also in the adiabatic model we could take a much larger limit on the size of the CWR, and in that case we would also get a double wave pattern. As this does not fit the observations, however, we did not present such models. Applying the Blomme et al. (2010) alternative stellar wind parameters (lower mass-loss rate for the primary) in the radiative model results in a similar double-wave pattern.

5.5. CWR combined synchrotron and adiabatic model

We can now extend the synchrotron model of Sect. 5.2 to also include the free-free emission and absorption of the CWR. In the Blomme et al. (2010) computer code, we find those cells where the synchrotron emission is non-zero, and we add the free-free opacity and emissivity, using a density four times higher than the local stellar wind density. We consider the temperature of the CWR wind to be a free parameter.

The resulting 3 mm fluxes are shown in Fig. 3 (light blue curve). We used a CWR temperature of 7.0 × 104 K, which is lower than the value used in Sect. 5.3 to avoid having fluxes that are too high. The combined model is comparable in quality of fit to the pure synchrotron one. The phase of maximum is shifted slightly toward the observed maximum, but is still quite some distance from it.

5.6. Comparison to Pittard (2010)

Pittard (2010) used his 3D hydrodynamical models of colliding-wind binaries to predict the thermal emission from submillimetre to centimetre radio wavelengths. Three of the four models he calculated have circular orbits, and are therefore not applicable to the eccentric binary Cyg OB2 #8A. This is further confirmed by the fact that all three predicted 3 mm light curves show a double peak, which is not seen in the Cyg OB2 #8A data.

The Pittard (2010) cwb4 model is an eccentric binary, though with parameters different from Cyg OB2 #8A: it is an O6V + O6V system with a period of 6.1 days, and an eccentricity of . Each of the components has a mass loss rate of = 2 × 10-7M yr-1. and a terminal velocity of ν = 2500 km s-1. Nevertheless, the hydrodynamics are somewhat similar to what we expect for Cyg OB2 #8A, in that the CWR is radiative at periastron and adiabatic at apastron. A qualitative comparison is therefore possible.

The cwb4 model shows a single peak in the 3 mm light curve, that is close to periastron. The exact orbital phase of the peak is only slightly dependent on the viewing angle. The observed data however, peak around phase 0.8. The difference could be due to the fact that the cwb4 model is for two equal winds, while the Cyg OB2 #8A winds are not equal. This will shift the phase of maximum away from periastron, but detailed hydrodynamical calculations of the Cyg OB2 #8A system would be required to see if the maximum moves to the observed one. The flux level of the cwb4 peak is about three times the baseline flux. This compares well to the factor of about two seen in the Cyg OB2 #8A data.

6. Cyg OB2 #8B

On the map of our combined data (Fig. 1) we also detected Cyg OB2 #8B (=Schulte 8B = MT 462). We measured its flux on the combined visibility data, by fitting a model with two point sources (one for #8A, and one for #8B). We found a flux for Cyg OB2 #8B of 0.21 ± 0.033 mJy, with a position that is offset by only 0.̋4 from its SIMBAD position.

In order to convert the observed flux into a mass-loss rate, we need to know additional information about this star. Its spectral type is listed as O6 II(f) (Sota et al. 2011), O6.5 III(f) (Massey & Thompson 1991), or O7 III-II (Kiminki et al. 2007). There are no significant radial velocity changes (Kobulnicky et al. 2014), making it unlikely that it is a binary.

There are no stellar parameter determinations for this star, so we have to rely on calibrations to determine them. We follow Morford et al. (2016), who use the Martins et al. (2005) calibration to derive an effective temperature Teff = 35 644 K, a luminosity of log L/L = 5.49, and a (spectroscopic) mass of M/M = 33.68, based on the O6.5 III(f) spectral type. From the Prinja et al. (1990) calibration of terminal velocity versus spectral type, ν = 2545 km s-1 follows. Using this in Eq. (1), we find = 1.42 ± 0.2 × 10-6M yr-1.

We can compare this to the upper limits that have been derived from radio observations. Morford et al. (2016) derive an upper limit from their 21 cm observations of 0.078 mJy, corresponding to < 4.3 × 10-6M yr-1. The 0.2 mJy upper limit at 6 cm found by Bieging et al. (1989) leads to < 5.2 × 10-6M yr-1. There are no determinations from other mass loss indicators (such as Hα). It should be noted that all mass-loss rates listed here have not been clumping corrected.

The predicted mass-loss rate according to the Vink et al. (2001) formulae is 0.7 × 10-6M yr-1. The difference between this value and the observed one can be interpreted as due to clumping, with a clumping factor of four. This value is compatible with current ideas about clumping in massive-star winds (e.g., Puls et al. 2008).

7. Conclusions

We monitored the massive colliding-wind binary Cyg OB2 #8A at 3 mm with the NOEMA interferometer, with good phase coverage of its orbit. For 12 of the 14 observations, we could determine the flux.

The 3 mm light curve shows clear phase-locked variability, indicating that a substantial part of the flux comes from the colliding-wind region (CWR)3. We modelled the data using models that include synchrotron radiation due to the CWR, or free-free emission from the CWR, or both. All models give fluxes and flux ranges higher than those observed, though comparable in magnitude to the observations. With the exception of the radiative shock model, all have a single peak – as observed. The slightly better agreement of the shape of the maximum would favour the synchrotron interpretation. From a modelling point of view, however, some contribution of free-free emission from the CWR is also expected. In summary, stronger arguments will be required for a decision on how much each of the two emission mechanisms contribute to the observed fluxes.

The main problem with the models presented here is that the theoretical maxima are phase shifted with respect to the observed one. It is clear that, at least qualitatively, this is due to our neglect of the orbital motion on the shape of the CWR. To see if this is also quantitatively the correct explanation, more sophisticated modelling will be required. This should be based on solving the hydrodynamical equations, in order to derive the density and temperature distribution in the CWR. These can then be used to calculate a better theoretical light curve and will allow us to determine the relative contributions of the synchrotron and free-free radiation emission.

On the deep 3 mm image made by combining all our data, we also detected the visual companion Cyg OB2 #8B. We measured a flux of 0.21 ± 0.033 mJy, which leads to a (non-clumped) mass-loss rate of = 1.42 ± 0.2 × 10-6M yr-1. From a comparison with predicted mass-loss rates, we derive a clumping factor of four, which is compatible with current ideas about clumping in massive-star winds.


3

One could consider alternative explanations, such as phase-locked variability of the mass-loss rate, possibly caused by changing radiative inhibition in this eccentric binary (Stevens & Pollock 1994). However, also in that case a colliding-wind region would of necessity exist. The dominant effect of the CWR is shown by the fact that our models of the CWR can clearly get the magnitude of the flux variations correct – though many important details still remain to be solved.

Acknowledgments

We are especially grateful to Jan Martin Winters and Charlène Lefèvre (IRAM) for their assistance with the data reduction. We thank the referee for his/her comments, which helped to improve the paper. We have made use of the GILDAS software for the reduction of the data. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France and NASA’s Astrophysics Data System Abstract Service. J. Morford and D. Fenech wish to acknowledge funding from an STFC studentship and STFC consolidated grant (ST/M001334/1) respectively.

References

  1. Bell, A. R. 1978, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
  2. Benaglia, P., & Romero, G. E. 2003, A&A, 399, 1121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Benaglia, P., Marcote, B., Moldón, J., et al. 2015, A&A, 579, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Bieging, J. H., Abbott, D. C., & Churchwell, E. B. 1989, ApJ, 340, 518 [NASA ADS] [CrossRef] [Google Scholar]
  5. Blomme, R., & Volpi, D. 2014, A&A, 561, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Blomme, R., Van Loo, S., De Becker, M., et al. 2005, A&A, 436, 1033 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Blomme, R., De Becker, M., Runacres, M. C., Van Loo, S., & Setia Gunawan, D. Y. A. 2007, A&A, 464, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Blomme, R., De Becker, M., Volpi, D., & Rauw, G. 2010, A&A, 519, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Blomme, R., Nazé, Y., Volpi, D., et al. 2013, A&A, 550, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Daley-Yates, S., Stevens, I. R., & Crossland, T. D. 2016, MNRAS, 463, 2735 [NASA ADS] [CrossRef] [Google Scholar]
  11. De Becker, M. 2007, A&ARv, 14, 171 [NASA ADS] [CrossRef] [Google Scholar]
  12. De Becker, M., Rauw, G., & Manfroid, J. 2004, A&A, 424, L39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. De Becker, M., Rauw, G., Sana, H., et al. 2006, MNRAS, 371, 1280 [NASA ADS] [CrossRef] [Google Scholar]
  14. Dougherty, S. M., Pittard, J. M., Kasian, L., et al. 2003, A&A, 409, 217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Eichler, D., & Usov, V. 1993, ApJ, 402, 271 [NASA ADS] [CrossRef] [Google Scholar]
  16. Ginzburg, V. L., & Syrovatskii, S. I. 1965, ARA&A, 3, 297 [NASA ADS] [CrossRef] [Google Scholar]
  17. Kennedy, M., Dougherty, S. M., Fink, A., & Williams, P. M. 2010, ApJ, 709, 632 [NASA ADS] [CrossRef] [Google Scholar]
  18. Kiminki, D. C., Kobulnicky, H. A., Kinemuchi, K., et al. 2007, ApJ, 664, 1102 [NASA ADS] [CrossRef] [Google Scholar]
  19. Knödlseder, J. 2000, A&A, 360, 539 [NASA ADS] [Google Scholar]
  20. Kobulnicky, H. A., Kiminki, D. C., Lundquist, M. J., et al. 2014, ApJS, 213, 34 [NASA ADS] [CrossRef] [Google Scholar]
  21. Leitherer, C., & Robert, C. 1991, ApJ, 377, 629 [NASA ADS] [CrossRef] [Google Scholar]
  22. Madura, T. I., Gull, T. R., Okazaki, A. T., et al. 2013, MNRAS, 436, 3820 [NASA ADS] [CrossRef] [Google Scholar]
  23. Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Massey, P., & Thompson, A. B. 1991, AJ, 101, 1408 [NASA ADS] [CrossRef] [Google Scholar]
  25. Montes, G., González, R. F., Cantó, J., Pérez-Torres, M. A., & Alberdi, A. 2011, A&A, 531, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Montes, G., Alberdi, A., Pérez-Torres, M. A., & González, R. F. 2015, Rev. Mexicana Astron. Astrofis., 51, 209 [NASA ADS] [Google Scholar]
  27. Morford, J. C., Fenech, D. M., Prinja, R. K., Blomme, R., & Yates, J. A. 2016, MNRAS, 463, 763 [NASA ADS] [CrossRef] [Google Scholar]
  28. Parkin, E. R., & Pittard, J. M. 2008, MNRAS, 388, 1047 [NASA ADS] [CrossRef] [Google Scholar]
  29. Parkin, E. R., Pittard, J. M., Corcoran, M. F., & Hamaguchi, K. 2011, ApJ, 726, 105 [NASA ADS] [CrossRef] [Google Scholar]
  30. Pittard, J. M. 2007, ApJ, 660, L141 [NASA ADS] [CrossRef] [Google Scholar]
  31. Pittard, J. M. 2010, MNRAS, 403, 1633 [NASA ADS] [CrossRef] [Google Scholar]
  32. Pittard, J. M., & Dougherty, S. M. 2006, MNRAS, 372, 801 [NASA ADS] [CrossRef] [Google Scholar]
  33. Pittard, J. M., Dougherty, S. M., Coker, R. F., O’Connor, E., & Bolingbroke, N. J. 2006, A&A, 446, 1001 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Pollock, A. M. T. 1987, A&A, 171, 135 [NASA ADS] [Google Scholar]
  35. Prilutskii, O. F., & Usov, V. V. 1976, Sov. Astron., 20, 2 [NASA ADS] [Google Scholar]
  36. Prinja, R. K., Barlow, M. J., & Howarth, I. D. 1990, ApJ, 361, 607 [NASA ADS] [CrossRef] [Google Scholar]
  37. Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Rauw, G., & Nazé, Y. 2016, Adv. Space Res., 58, 761 [NASA ADS] [CrossRef] [Google Scholar]
  39. Reimer, A., Pohl, M., & Reimer, O. 2006, ApJ, 644, 1118 [NASA ADS] [CrossRef] [Google Scholar]
  40. Reitberger, K., Kissmann, R., Reimer, A., & Reimer, O. 2014, ApJ, 789, 87 [NASA ADS] [CrossRef] [Google Scholar]
  41. Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011, ApJS, 193, 24 [NASA ADS] [CrossRef] [Google Scholar]
  42. Stevens, I. R. 1995, MNRAS, 277, 163 [NASA ADS] [Google Scholar]
  43. Stevens, I. R., & Pollock, A. M. T. 1994, MNRAS, 269, 226 [NASA ADS] [CrossRef] [Google Scholar]
  44. Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265 [NASA ADS] [CrossRef] [Google Scholar]
  45. Van Loo, S., Runacres, M. C., & Blomme, R. 2005, A&A, 433, 313 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Van Loo, S., Blomme, R., Dougherty, S. M., & Runacres, M. C. 2008, A&A, 483, 585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. White, R. L., & Becker, R. H. 1995, ApJ, 451, 352 [NASA ADS] [CrossRef] [Google Scholar]
  49. Wright, A. E., & Barlow, M. J. 1975, MNRAS, 170, 41 [NASA ADS] [CrossRef] [Google Scholar]
  50. Wright, N. J., Drew, J. E., & Mohr-Smith, M. 2015, MNRAS, 449, 741 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Observing log of the NOEMA 3 mm observations of Cyg OB2 #8A.

Table 2

Measured 3 mm fluxes of Cyg OB2 #8A.

Table 3

Overview of models used in Sect. 5.

Table 4

Cyg OB2 #8A parameters used in modelling.

All Figures

thumbnail Fig. 1

Deep 3 mm image of all Cyg OB2 #8A data combined. The colour scale values are in Jy. The contour levels are at −0.1 mJy, + 0.1 mJy and then go up in steps of 0.1 mJy to 1.6 mJy. The negative contour is indicated by the dashed, white line. The crosses indicate the positions of Cyg OB2 #8A and #8B (from SIMBAD). The synthesized beam (shown in the lower left corner) is 2.93″ × 1.76″ with a position angle of 57.̊1.

Open with DEXTER
In the text
thumbnail Fig. 2

Observed 3 mm fluxes of Cyg OB2 #8A, plotted as a function of orbital phase in the 21.908-day binary period. The blue symbols show the 2014 data, the red ones the 2016 data. To better show the behaviour of the light curve, the phase range is extended by 0.2 on each side. Open symbols indicate duplications in that extended range. The grey data show the 6 cm observations from Blomme et al. (2010). Note that the 3 mm flux scale (left) is different from the 6 cm one (right).

Open with DEXTER
In the text
thumbnail Fig. 3

Theoretical fluxes compared to the observed 3 mm fluxes of Cyg OB2 #8A, plotted as a function of orbital phase in the 21.908-day binary period. The black symbols with error bars show the observed data. Open symbols indicate duplications in the extended phase range (as in Fig. 2). The dark blue curve is the synchrotron model, the red curve the adiabatic thermal emission model, the green curve the radiative thermal emission model, and the light blue curve the combined synchrotron and adiabatic thermal emission model.

Open with DEXTER
In the text
thumbnail Fig. 4

CWR of the thin radiative shock model projected on to the sky, for phases where the flux is high (left column) and where the flux is low (right column). The orbital phase is listed in the top left corner of each panel. For plotting purposes, the CWR has been cut off at a distance of 600 R.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.