A&A 519, A111 (2010)

Non-thermal radio emission from O-type stars

IV. Cygnus OB2 No. 8A^,

R. Blomme¹ - M. De Becker^2,3 - D. Volpi¹ - G. Rauw²

1 - Royal Observatory of Belgium, Ringlaan 3, 1180 Brussel, Belgium
2 - Institut d'Astrophysique, Université de Liège, Allée du 6 Août, 17, Bât. B5c, 4000 Liège (Sart-Tilman), Belgium
3 - Observatoire de Haute-Provence, 04870 Saint-Michel l'Observatoire, France

Received 12 October 2009 / Accepted 27 May 2010

Abstract
Context. Several early-type colliding-wind binaries are known to emit synchrotron radiation due to relativistic electrons, which are most probably accelerated by the Fermi mechanism. By studying such systems we can learn more about this mechanism, which is also relevant in other astrophysical contexts. Colliding-wind binaries are furthermore important for binary frequency determination in clusters and for understanding clumping and porosity in stellar winds.
Aims. We study the non-thermal radio emission of the binary Cyg OB2 No. 8A, to see if it is variable and if that variability is locked to the orbital phase. We investigate if the synchrotron emission generated in the colliding-wind region of this binary can explain the observations and we verify that our proposed model is compatible with the X-ray data.
Methods. We use both new and archive radio data from the Very Large Array (VLA) to construct a light curve as a function of orbital phase. We also present new X-ray data that allow us to improve the X-ray light curve. We develop a numerical model for the colliding-wind region and the synchrotron emission it generates. The model also includes free-free absorption and emission due to the stellar winds of both stars. In this way we construct artificial radio light curves and compare them with the observed one.
Results. The observed radio fluxes show phase-locked variability. Our model can explain this variability because the synchrotron emitting region is not completely hidden by the free-free absorption. In order to obtain a better agreement for the phases of minimum and maximum flux we need to use stellar wind parameters for the binary components which are somewhat different from typical values for single stars. We verify that the change in stellar parameters does not influence the interpretation of the X-ray light curve. Our model has trouble explaining the observed radio spectral index. This could indicate the presence of clumping or porosity in the stellar wind, which - through its influence on both the Razin effect and the free-free absorption - can considerably influence the spectral index. Non-thermal radio emitters could therefore open a valuable pathway to investigate the difficult issue of clumping in stellar winds.

Key words: stars: individual: Cyg OB2 No. 8A - stars: early-type - stars: mass-loss - radiation mechanisms: non-thermal - acceleration of particles - radio continuum: stars

1 Introduction

All hot stars emit radio radiation due to thermal free-free emission by the ionized material in their stellar wind. A number of those stars also show non-thermal radio emission that dominates the thermal component. This non-thermal emission is believed to be due to electrons that are Fermi-accelerated around shocks (Bell 1978). As these relativistic electrons spiral around in the magnetic field, they emit synchrotron radiation, which we detect as non-thermal radio emission (Bieging et al. 1989).

Two hypotheses have been proposed for the origin of the shocks responsible for the non-thermal emission. These shocks could be due to the collision of the two winds in a binary system (Dougherty et al. 2003; Reimer et al. 2006; Pittard et al. 2006; Eichler & Usov 1993), or they could be formed by the radiative instability mechanism in the wind of a single star (Owocki & Rybicki 1984; Chen & White 1994). The binary hypothesis is supported by the non-thermal Wolf-Rayet stars (Dougherty & Williams 2000) and by a number of known O+O binaries showing non-thermal radio emission (e.g. Cyg OB2 No. 5, Contreras et al. 1997). The single-star hypothesis, on the other hand, is supported by a number of seemingly single O stars showing non-thermal radio emission.

In recent years, however, doubt has been cast on the single-star hypothesis. In previous papers of this series, we have shown that the radio emission for some of these stars shows periodic variability, suggesting binarity (HD 168 112 and Cyg OB2 No. 9, Van Loo et al. 2008; Blomme et al. 2005). Optical spectroscopy has shown that a number of these ``single'' stars are actually binaries (e.g. Cyg OB2 No. 8A and No. 9, Nazé et al. 2008; De Becker et al. 2004). Furthermore, theoretical work by Van Loo et al. (2006) showed that non-thermal radio radiation from a single star would be unlikely, as the synchrotron emission would be largely absorbed by the stellar wind material.

The study of colliding-wind systems is important for a number of reasons. First, because we can learn more about the Fermi acceleration mechanism, which is relevant to a wide range of astrophysical objects such as interplanetary shocks and supernova remnants (e.g. Jones & Ellison 1991). The physical conditions in a colliding-wind binary are quite different from those in other objects. The different range of density, magnetic field and radiation environment in these binaries provides substantial new information on the acceleration mechanism. Models, constrained by radio data, also provide predictions for the expected gamma-ray flux from these systems. These numbers are important for observations of colliding-wind binaries by current and future high-energy instruments, both ground-based and space-based ones (e.g. INTEGRAL, IXO, HESS, VERITAS).

Furthermore, the non-thermal radio emitters are useful for improving the determination of the binary frequency in clusters. The colliding-wind binaries allow us to find the longer-period systems, which are easy to detect at radio wavelengths but much harder to find with standard spectroscopic observing techniques.

Finally, colliding-wind systems are also highly relevant for 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). By determining how much of the synchrotron radiation emitted at the colliding-wind region is absorbed by the stellar wind material, we can estimate the amount of porosity in the wind.

In this paper we study the non-thermal radio emission of the colliding-wind binary Cyg OB2 No. 8A (RA = 15 0789, Dec = +4118'50 494, J2000). We also present additional X-ray observations that improve the X-ray light curve presented by De Becker et al. (2006). Cyg OB2 No. 8A is undoubtedly the O + O non-thermal radio emitter with the best constrained stellar, wind and orbital parameters so far (De Becker 2007). It has been known as a non-thermal radio emitter since the first major survey of O star radio emission by Bieging et al. (1989). It was recognized as such by its clearly negative spectral index. This system was only recently discovered to be a binary by De Becker et al. (2004). It consists of an O6If primary and an O5.5III(f) secondary. It has a period of  d and an eccentricity of . All binary parameters were determined by De Becker et al. (2006) and are listed in Table 1.

Table 1:   Parameters of Cyg OB2 No. 8A used in this work.

De Becker et al. (2006) studied the X-ray emission of this binary. They found an essentially thermal spectrum, but with an overluminosity of a factor 19-28 compared to the canonical ratio for non-interacting O-type stars. This overluminosity is due to additional X-ray emission by material heated in the colliding-wind region. The X-ray light curve shows variability of 20%. Although the curve is not well sampled, it suggests that the variability is phase-locked.

In Sect. 2 we present the radio continuum observations of Cyg OB2 No. 8A and in Sect. 3 the new X-ray data. The model for the synchrotron radio emission in a colliding-wind binary is described in Sect. 4. In Sect. 5 we apply this model both for a standard set of star and wind parameters and an alternative set chosen to obtain better agreement with the data. Conclusions are given in Sect. 6.

2 Radio data

We obtained a number of 6 cm continuum observations of Cyg OB2 No. 8A with the NRAO Very Large Array (VLA) during 2005 February 04 to March 12, covering about 1.6 orbital periods. We also retrieved a considerable number of observations from the VLA archive (also covering other wavelengths), which we used in the present analysis. To avoid introducing systematic effects, we re-reduced the archive data in the same way as our own data (see Appendix B). The resulting fluxes and their error bars are listed in Table B.1.

In Fig. 1, we plot a subset of these fluxes. We limit the figure to the 3.6 and 6 cm data as the other wavelengths have only a few observations or detections. We do not plot fluxes which were measured on images where Cyg OB2 No. 8A was offset from the centre. Though these observations follow the same general trend as the ``on-centre'' observations, they usually have large error bars and would therefore considerably confuse the figure.

Figure 1 shows clear variability in the fluxes which is phase-locked with the orbital period. This is especially obvious at the 6 cm wavelength which has the most data. Maximum flux is at phase 0.1 (where the primary is in front). The phase of minimum flux is less well determined due to a gap in the phase coverage, but is at 0.6 (close to where the secondary is in front). To check that the variability is not an instrumental effect, we also measured two other targets on the images: Cyg OB2 Nos. 63 and 9. Although these are also variable (No. 9, see Van Loo et al. 2008), their variability is not correlated with the No. 8A variability. We can therefore exclude instrumental effects.

Only a few simultaneous observations at more than one wavelength are available to study the spectral index. We find that before phase 0.5, the 3.6 and 6 cm fluxes are roughly equal, corresponding to a spectral index of about 0.0 (within the error bars). After phase 0.5, the fluxes behave differently: the 3.6 cm flux reaches higher values much more quickly than the 6 cm one, which translates into a spectral index of +0.7. Interestingly, this is close to the orbital phase where the only significantly negative spectral index (-0.7) is found, occurring between the 6 and 20 cm observations of programme AC116 (1984-12-21, phase =0.51).

Table 2:   Two new X-ray observations of Cyg OB2 No. 8A.

We next try to determine the period from the 6 cm data used in Fig. 1. To do so, we use the string-length method of Dworetsky (1983). The best period we find is 23.75 d, which is quite different from the 21.908-d value derived from optical spectroscopy (De Becker et al. 2004). However, the minimum in the string-length is not very significant. Furthermore, if we replot the De Becker et al. spectroscopic velocities on the 23.75-d period, we find that the points are scattered all over the diagram, indicating that the proposed radio period is clearly wrong. We therefore use the 21.908-d period in the further analysis.

3 X-ray data

Cyg OB2 No. 8A has been observed in X-rays two times with the XMM-Newton satellite since the publication of the De Becker et al. (2006) paper (Obs. ID 050511, PI: G. Rauw). These observations were scheduled in order to fill the gap in phase coverage of previous X-ray observations, and occurred during revolutions 1353 and 1355 (respectively at phases 0.155 and 0.326, see Table 2 for details). We adopted the same procedure as in for the data processing using the Science Analysis Software (v.6.0.0), and we refer to this study for details. We note that we rejected short time intervals of high background due to soft proton flares. We fitted the new EPIC spectra with the same three-temperature model as used by De Becker et al. (2006) for previous XMM-Newton observations. The best fit parameters are given in Table 3; they are consistent with those derived for the four previous observations, even though some variability is detected for some of them.

Table 3:   Parameters for the two new EPIC spectra of Cyg OB2 No. 8A.

In Fig. 2, we have plotted the parameters obtained from the best fit of EPIC data (MOS1 and pn) of the six XMM-Newton observations (see Table 4 in De Becker et al. 2006; and Table 3 in this paper) versus the orbital phase. The local absorption column density undergoes significant variability as a function of the orbital phase, with a maximum occurring close to periastron passage (phase 0.0). At that time the wind interaction region is behind the primary wind, which is likely to be responsible for a substantial absorption of X-rays. The plasma temperatures show only weak, or even no phase-locked variability. Only slight variations are suggested in the case of the hard emission component that might be related to the small changes in the pre-shock velocities as the separation between the stars changes along the eccentric orbit. The lower and higher temperature are indeed found, respectively, close to periastron and apastron for the hard component, but this trend should be considered with caution in view of the scatter and of the error bars of the measurements. Finally, the normalization parameters (which are related to the emission measure) undergo significant phase-locked variations, mostly if one considers the soft emission component that dominates the X-ray spectrum. The maximum coincides remarkably with orbital phases close to periastron, in agreement with the idea that the maximum X-ray emission occurs when the emission measure is expected to be maximum.

We proceeded following the same approach as in De Becker et al. (2006) to build a light curve including ROSAT-PSPC, ASCA-SIS and XMM-Newton-EPIC data (see their Fig. 6). In order to consistently compare the count rates obtained with different instruments, we used the first XMM-Newton observation as a reference (obs. 1 in De Becker et al., obtained at orbital phase 0.534). We used the 3-temperature model with the parameters obtained for the simultaneous fit of EPIC-MOS1 and EPIC-pn data of that observation, and we convolved it with the respective response matrices of ROSAT-PSPC and ASCA-SIS1 to obtain faked spectra. We then normalized the count rates of the ROSAT and ASCA observations using the count rate of the reference fake spectrum in each case, thereby deriving the normalized count rates plotted in Fig. 1. It is interesting to note that the X-ray emission level at orbital phases 0.155 and 0.326 (respectively Nos. 5 and 6 in the XMM-Newton series) is not lower than at other orbital phases. This suggests that the X-ray minimum occurs very close to periastron.

An important result of Fig. 1 is that the X-ray and radio light curves are nearly anti-correlated. X-ray flux minimum occurs close to periastron, when the radio 6 cm flux is approaching maximum. Minimum radio flux is near phase 0.6, where the X-ray flux has increased well above its minimum value (though it has not yet reached maximum). The anti-correlation will present a challenge to modelling this binary system. Although we do not model the X-ray observations in the present paper, we discuss this point further in Sect. 5.3.

4 Modelling the radio fluxes

Blomme (2005) made a very simplified model for the radio emission of the Cyg OB2 No. 8A binary system. This model consists of a point-source synchrotron emitting region orbiting a single star inside a stellar wind that includes free-free absorption. The model shows that all synchrotron emission would be absorbed by the stellar wind material and that therefore the radio fluxes should not show phase-locked variability. Only by introducing a considerable amount of porosity in the stellar wind material does the non-thermal radio emission become detectable (due to the reduction in free-free absorption).

Here we present a more sophisticated model of Cyg OB2 No. 8A that takes into account the geometric extent of the colliding-wind region, as well as the acceleration, advection and cooling of the relativistic electrons. In this section we provide an outline of the model; for details we refer to Appendix A.

The present model does not solve the hydrodynamical equations, but instead uses the Antokhin et al. (2004) equations to define the position of the contact discontinuity that separates the two stellar winds. Estimates by De Becker et al. (2006, their Table 11) show the ratio of cooling time to the flow time to be around 1. At least during part of the orbit, the shocks can therefore be expected to be radiative, resulting in considerable variability and structure. Such complications are neglected in our model.

Next we assume that the contact discontinuity and the shocks on either side of it are sufficiently close to one another that we can consider them to be in the same geometric place. We then generate new relativistic electrons at the shocks and follow their time evolution as they advect and cool, till they no longer emit significant synchrotron radiation at the wavelength we are considering.

At the shock, the new electrons have a momentum distribution that is a power law modified for the effect of inverse Compton cooling. The number of relativistic electrons is determined by , the fraction of the energy that is transferred from the shock to the relativistic electrons. We assume the typical value of (Blandford & Eichler 1987; Eichler & Usov 1993). The time evolution of the electrons is followed as they advect outward along the contact discontinuity and as they cool down due to inverse Compton and adiabatic cooling.

At each point along the tracks followed by the electrons, the synchrotron emissivity can then be calculated. We include the Razin effect, which takes into account the effect of the non-relativistic plasma on the synchrotron emission. This emissivity is then mapped from the two-dimensional orbital plane on which it was defined on to a three-dimensional volume by rotating it along the line connecting both components of the binary. Note that this assumption of cylindrical symmetry neglects any effects due to orbital motion. We refer to Lemaster et al. (2007), Parkin & Pittard (2008) and Pittard (2009) for studies of the orbital motion effects in colliding-wind systems in general. Specifically for the relatively short-period system Cyg OB2 No. 8A, not taking into account orbital motion is a major assumption, and we discuss its effect in Sect. 5.2.

The model also includes the thermal free-free opacity and emissivity due to the ionized wind material. The winds are assumed to be smooth, i.e. no clumping or porosity effects are included. The model does not include the contribution due to the increased density and temperature in the colliding-wind region (contrary to Pittard 2010) as these hydrodynamical quantities are not calculated in our model. Radiative transfer is then solved in the three-dimensional simulation volume using Adam's finite volume method (Adam 1990).

5 Results

The intention of our modelling is to see if we can reproduce qualitatively the main features of the observed Cyg OB2 No. 8A radio light curve (see Sect. 2). These are: (1) the fact that there is variability; (2) maximum flux is around phase 0.1, minimum flux around phase 0.6; and (3) the 3.6 cm and 6.0 cm fluxes are approximately equal, most of the time (i.e. spectral index 0). We do not make a formal best fit as we consider the absolute values of the fluxes to be less important, since these are sensitive to many of the assumptions in the model (e.g. magnetic field, ).

5.1 Standard model

5.1.1 Model details

The standard model uses the parameters listed in Table 1. We recall that De Becker et al. (2006) derived these from stellar parameters typical for single stars with the same spectral type as the binary components of Cyg OB2 No. 8A. We assume a surface magnetic field for each star of B=100 G, which is compatible with the fact that the magnetic field of most O-type stars is below detectable levels (e.g. Schnerr et al. 2008). The magnetic field at the shocks is directly related to this surface magnetic field (see Eq. (A.7) and Fig. 3, upper panel); this neglects any amplification of the field (such as described by, e.g. Bell 2004) or any reduction of it by magnetic reconnection. As the true rotational velocity is not known, we make the assumption that .

A comparison between the values derived from the binary orbit and typical masses of single stars corresponding to the observed spectral types indicates an inclination angle of 32 (De Becker et al. 2006). For the free-free emission and absorption we take the material to be once ionized and assign it a temperature of 28 500 K (0.75 of the effective temperature of the stars). The three-dimensional simulation volume is a cube centred on the primary that is on each side and that consists of 256³ cells. Some care must be taken to use the appropriate resolution for a given size of the simulation box; we checked that the results presented here are robust to changes in size and resolution. The size of the box is large enough to contain most, but not all, of the emission. Doubling the size and using 512³ cells shows that the flux loss is less than 10%, and that the shape of the radio light curve remains the same. With the resolution we use for the three-dimensional simulation volume, the apex of the contact discontinuity is not well resolved. This is of no consequence for our purpose, however, as the synchrotron emission coming from around the apex is completely absorbed by the free-free absorption (see also Sect. A.7).

5.1.2 Fluxes and spectral index

Figure 4 (left column) shows the resulting 3.6 cm and 6 cm fluxes as well as the spectral index between these two wavelengths. The first thing to note is that the model fluxes do show phase-locked variability, contrary to the simple point-source model (Blomme 2005). This is due to the synchrotron emission extending much further out than the point-source approximation suggests. In Fig. 5 (left column) we plot the contribution to the flux as a function of distance to the primary. To determine each flux value on that figure we artificially set the synchrotron emission to zero beyond a certain distance and re-run the radiative transfer model. The figure makes clear that the main flux contribution comes from distances of . This corresponds to approximately 10 - 30 times the separation between the two stars. The flux contribution that can be seen at small distances is due to free-free emission of the entire wind, which was not cut off in this procedure. In the figure we also indicate the distances from the star where optical depth 1 for the primary component is reached. We calculate these using the Wright & Barlow (1975) formalism and find them to be at 3.6 cm and at 6 cm. As expected, the flux we receive is formed largely outside this radius.

Spatially resolved radio observations of colliding-wind regions typically do not show such an extended non-thermal region (e.g. WR 140, Dougherty et al. 2005). The present results are not in contradiction with this. The resolved observations show the surface brightness (flux per area) while the unresolved radio data presented here show the flux (i.e. integrated over a large area). At large distances from the centre of the binary system, the surface brightness decreases substantially, hence the relatively small extent of the colliding-wind region in resolved observations - which is nicely illustrated by the simulated data for WR 140 presented by Pittard & Dougherty (2006, their Fig. 17). But, as the synchrotron emission is cylindrically symmetric around the axis connecting the two stars, the low surface brightness at larger distances gets weighted with a much larger surface area and therefore contributes substantially to the flux.

While the phase-locked variability of the flux is a positive result for the model, unfortunately we find that the phases of minimum and maximum flux do not agree with the observations. The model results are nearly in anti-phase to the observations: theoretical maximum occurs at phase 0.55 (instead of the observed 0.1) and theoretical minimum at phase 0.0 - 0.05 (instead of 0.6). Furthermore, the spectral index is not correctly predicted. The theoretical values of +0.75 to +1.0 do not agree with the observed values which are mostly close to 0.0 (Sect.  2).

5.1.3 Without Razin effect, without free-free absorption

To understand why the phases and spectral index are incorrectly modelled, we start by simplifying the model. When we leave out the free-free absorption and emission as well as the Razin effect, we find the light curve shown in the middle column on Fig. 4. The flux values are high (700 mJy at 6 cm, 400 mJy at 3.6 cm). The spectral index is approximately -1.0, clearly indicating the non-thermal nature of the spectrum.

Figure 5 (middle column) shows that in the model without Razin and without free-free absorption, the strongest contribution comes from relatively close to the apex (up to 300 , i.e. about 3 times the binary separation). This is to be expected as the collision is strongest close to the apex, so more energy is available for the acceleration of the particles. On the other hand, inverse Compton cooling should be strongest there as well, leading to a reduction of the emission.

The resulting balance between acceleration and inverse Compton cooling can be seen in Fig. 3 (lower panel), where we plot the maximum momentum () that can be reached due to inverse Compton cooling at the shock where the electrons are accelerated. We plot this value as a function of distance along the contact discontinuity ( ). For comparison we also plot the ``characteristic'' momentum for 3.6 and 6 cm emission (p₆ and p_3.6), i.e. the momentum at which the synchrotron emission is largest for a given wavelength. Its value can be derived from the fact that maximum emission occurs at frequency (see, e.g. Van Loo 2005, Eq. (2.52)), where is given by Eq. (A.18). Close to the apex, is considerably larger than p₆ (and p_3.6). Because of the (approximate) power-law distribution of the momentum, there are a large number of electrons available having their maximum emission at that wavelength. Moving towards the outer regions, becomes closer to, or even lower than, p₆ (and p_3.6), thereby explaining the decrease of the emission.

Figure 4 (middle column) shows that the maximum flux occurs at apastron (phase 0.5), the minimum at periastron (phase 0.0). At apastron, the stars are further apart and their winds have reached the highest speed before they collide. However, this effect cannot explain why maximum occurs at apastron. The energy available for the relativistic electrons is proportional to the shock energy , which scales as . The higher value for the velocity at apastron cannot compensate for the lower 1/r² value. The main reason for maximum to occur at apastron is the (approximate) scale-invariance of the colliding-wind region. As the distance between the two stars becomes larger, the size of the synchrotron emitting region grows accordingly. Although the synchrotron emissivity is lower at apastron than at periastron, the integration over the larger size of the colliding-wind region results in a higher apastron flux.

5.1.4 With Razin effect, without free-free absorption

Next, we switch back on the Razin effect but we still leave out the free-free mechanism (Figs. 4 and 5, right column). This results in a substantial reduction of the flux levels (e.g. the 6 cm maximum flux goes down from 700 mJy to 0.4 mJy). The large influence of the Razin effect on the flux can be understood by making a few simple estimates. The crucial quantity is the factor f defined in Eq. (A.16):

where is the radio frequency, the mass of the electron, c the speed of light, p the momentum of the electron and . In the definition, e is the charge of the electron and is the number density of the thermal electrons (i.e. those of the non-relativistic plasma). The f-factor takes on values between 0 and 1 with f=1 indicating no effect, and indicating a strong Razin effect. Taking the frequency at which the electron has its maximum emission, this translates into , where

(see also Van Loo 2005, Eq. (2.53)). We use Eq. (A.7) for the magnetic field B and Eq. (A.20) to relate to the mass density . Determining from the mass-loss rate and approximating the velocity by the terminal velocity results in a r^-1 scaling of . Putting in the values for the primary (Table 1) and using for the frequency corresponding to 6 cm, we find a critical radius of 1100 . The Razin effect is important below that radius, which is indeed where we see the largest effect in the model calculations (Fig. 5, comparing middle and right column).

In this model without free-free absorption, we note that the Razin effect does not change the position of maximum flux, which remains at phase 0.5. It does have an important effect on the spectral index, which becomes clearly positive (+1.9). The large change in the spectral index can be explained by the frequency dependence of the Razin effect. The frequency dependence of f (Eq. (1)) is such that , so the 6 cm flux is more affected than the 3.6 cm one. This leads to a strong change in the spectral index. Note that this change is highly sensitive to . To quantify this sensitivity, we calculate an artificial model without free-free absorption where we reduce the value in the equation by a factor 30. This leads to a spectral index that is below 0.2 during that part of the orbit which is near flux maximum. We recall that the values used for in our model (Eq. (A.20)) are not based on hydrodynamics. The true values could therefore be very different, which will certainly influence the absolute fluxes as well as the spectral index.

The variability seen in Fig. 4 (right column) is phase-locked, with minimum flux occurring at periastron. Without the Razin effect (see Sect. 5.1.3), this is mainly due to the smaller size of the synchrotron emitting region. As pointed out above, the Razin effect is more important in regions of higher density (of the non-relativistic plasma) and will therefore reduce the flux at periastron more than at apastron. This explains why the minimum still occurs at periastron (Fig. 4, right column) and why the flux contrast between apastron and periastron is higher than in the model without Razin effect (Fig. 4, middle column).

5.1.5 With Razin effect, with free-free absorption

In the final step, we again switch on the free-free absorption and emission to obtain the results in Fig. 4, left column. The absorption decreases the flux levels, but it does not introduce a shift in maximum or minimum phase (this is discussed further in Sect. 5.2). The spectral index is now somewhat less positive compared to the model without free-free absorption. This shows that both the Razin effect and free-free absorption influence the observed spectral index.

A more detailed look at the emission (Fig. 6) reveals that most of the flux we detect is formed at some distance from the shock and must therefore be due to electrons that have cooled down somewhat. The figure shows the synchrotron emissivity in a coordinate system linked to the contact discontinuity. We limit the figure to distances of 1000-3000 , which is where most of the detectable radiation comes from.

Some care must be taken in the interpretation of this figure. The tracks show the particles to be moving away from the contact discontinuity/shock (which we assume to be at the same geometric position). This differs from the true hydrodynamical situation, where the particles move from the shock towards the contact discontinuity. This difference is a direct consequence of the simplifications in our model. It has a negligible influence, however, on the resulting fluxes and spectral index because the thickness of the synchrotron emitting region is quite small.

For a given point on the contact discontinuity (i.e., for a given ), the highest values are found at from the shock (Fig. 6, dotted white line). The electrons responsible for this emission can be traced back to their acceleration point, which is some 300  earlier (in ). In addition to the strong changes orthogonal to the contact discontinuity, there is also a much weaker decreasing gradient as we move from the inside to the outside. In the flux calculation the outer regions get higher weight as this emissivity is rotated in to a three-dimensional volume.

To understand why the maximum is offset from the shock we simplify the situation by assuming that the electrons emit only at their maximum frequency and we concentrate on p₆, the characteristic momentum for 6 cm emission. At the shock at position ``1'' (Fig. 6), a certain number density, N<sub>1</sub>(p<sub>6</sub>), of relativistic electrons is responsible for the 6 cm emission. Neglecting a weak function of , this number density is proportional to (Eq. (A.6)) with the subscript indicating the position. As these electrons advect out to position ``2'', their number density decreases in proportion to . Furthermore, while advecting, these electrons cool and therefore another set of electrons becomes responsible for the 6 cm emission. These are electrons that had p > p₆ at position ``1''; their number density is less, so we introduce a factor (<1) to take that into account. At position ``2'' their number density is N₂(p₆), which is proportional to . Position ``3'' is nearly the same as ``2'', but it is on the shock. There, new relativistic electrons are generated with (where we used ). Whether position ``3'' or ``2'' has the highest number density translates in to a comparison of and . The gradient is such that (Fig. 3, upper panel), so the result depends on whether the cooling is strong enough to overcome the effect of the gradient. The model calculations show that this is not the case and therefore N₂ > N₃, resulting in a maximum emission that is offset from the shock. Only at large distances (beyond those shown in Fig. 6), does the compensating effect no longer work and there the models show the maximum to be at the shock. This, however, happens beyond the geometric region that contributes significantly to the observed radio flux.

It may seem surprising that the contribution of the outer regions of the wind is still detectable, considering that the energy going into the shock is small (as velocities are nearly tangential to the contact discontinuity). We can make a simple estimate to compare the inner and outer wind regions, based on the energy available from the shock. At a position r along the contact discontinuity, we consider a section of length . In a time interval a volume of material crosses the shock through that section. The energy input rate is then: of which a fraction is used to accelerate the particles. The velocity orthogonal to the shock ( ) is proportional to r^-1. We derive from the mass conservation equation in a spherically symmetric wind and approximate by . If we assume that some constant fraction of the energy of accelerated particles goes into the synchrotron emission, we find that the radio luminosity of section is proportional to An outer region of size therefore contributes a fraction 10^-3 of an inner region of size . If we scale this to the calculated flux value of 700 mJy for the inner region (Fig. 5, middle column), we find 0.7 mJy for the outer region. None of the inner emission is detected due to the Razin effect and free-free absorption. The resulting flux is therefore predicted to be 0.7 mJy, which is consistent with the detailed models.

In summary, the standard model shows that phase-locked radio flux variability is indeed possible in Cyg OB2 No. 8A. The model fails however in predicting the phases of maximum/minimum flux and the spectral index.

5.2 An improved model

One possible approach to solve the problem with the phases of maximum/minimum flux is to ``mirror-reverse'' the system by exchanging the wind parameters of the primary and secondary star. It is important to realize that the star and wind parameters listed in Table 1 were not derived directly from information of the binary system, but rather from typical values for single stars with the same spectral type as the binary components. It is not clear how well such a procedure works for determining stellar parameters of binary components. This leads us to a more general approach, where we also explore models with star and/or wind parameters that are different from those listed in Table 1.

In Fig. 7 we present one possible improved solution where we change the mass loss rate of the primary from to and the terminal velocity from 1873 to 2500 . We present only this model from the large number of possibilities we explored, not because it presents a best-fit solution, but rather because it requires only minimal changes to the parameters compared to our standard model. Note that these changes result in the primary now having a weaker wind (less ram pressure) than the secondary.

We see in Fig. 7 that the phase of maximum is now 0.25, which is closer to the observed value (phase 0.1). Also the position of minimum flux at phase 0.9 (observed: 0.6) is improved compared to the standard model. The main reason for this is the change in shape of the contact discontinuity as can be seen in Fig. 8. The upper part of Fig. 8 shows the position of the contact discontinuity for the standard model. The observer looks from the bottom of the figure towards the top. We neglect in this discussion that the observer is not in the orbital plane, but 58 (= ) above it. For the standard model we have maximum flux at phase 0.5 (Fig. 8, top left) because one arm of the contact discontinuity is then pointing in the direction of the observer and the star with the weaker wind (the secondary) is in front, so there is only a limited amount of absorption. At minimum (Fig. 8, top right) the contact discontinuity is well hidden by the stronger wind from the primary, which is in front.

The situation is quite different in the improved model. The secondary now has the stronger wind and therefore the contact discontinuity wraps around the primary. Maximum flux (Fig. 8, bottom left) again occurs when one arm of the contact discontinuity is pointing in the direction of the observer and the star with the weaker wind (primary) is in front, but this now happens at phase 0.25. At phase 0.5, the intrinsic synchrotron emission is even higher (see discussion Sect. 5.1.4), but the arm has rotated further along the orbit so that most of it is absorbed by the free-free absorption. At phase 0.9 (Fig. 8, bottom right) the contact discontinuity is hidden by the stronger wind from the secondary and therefore the flux is minimum.

Although the phases of minimum and maximum are now closer to their observed values, the agreement is still not very good. However, the present model does not take into account the effect of orbital motion. The importance of this effect is shown by, e.g. Pittard (2009). In the specific case of Cyg OB2 No. 8A, we have significant emission coming from 10-30 times the binary separation; orbital effects at these distances will definitely play a role. Because of the orbital motion, the contact discontinuity will get distorted resulting in a leading arm that is further from the strong-wind secondary than the trailing arm. This will shift the phases at which minimum and maximum occur. Consider, e.g. maximum flux, where the part of the contact discontinuity closest to the observer is pointing in the direction of the observer (phase 0.25 in Fig. 8, bottom left). Including orbital motion will push this further away along the orbit from the strong-wind secondary. In order to have the contact discontinuity pointing again in the direction of the observer, we need to move the secondary back somewhat in its orbit. Maximum flux will thus occur at an earlier phase. Including orbital motion will therefore result in a better agreement with the observed phases of maximum and minimum flux.

One remaining problem is the spectral index, which is still considerably different from the observed value of 0.0. The spectral index is the result of the combined influence of free-free absorption and the Razin effect. As we saw in Sect. 5.1.4, the spectral index is highly sensitive to the value of the thermal electron density , through the Razin effect. The physical basis of this effect is that the local thermal electrons result in a refraction index which changes the relativistic beaming of the synchrotron radiation thereby reducing the emitted power (Ginzburg & Syrovatskii 1965). In Sect. 5.1.5 we showed that introducing free-free absorption changes the spectral index as well.

In order to model the observed spectral index, the simplest possibility is to reduce the mass-loss rates of the stars. The reduced density will result in both less free-free absorption and a reduction of the Razin effect. As synchrotron emissivity is proportional to this suggests that the flux would go down. But this is easily compensated for by the flux increase due to the reduction of the Razin effect, as a comparison of Fig. 4 middle and right columns shows. In order to be consistent with observations of single stars, such a reduced-density wind will need to be clumped (e.g. Puls et al. 2008). This would complicate the proposed explanation: the clumping would again increase the free-free absorption, probably to levels comparable to that of a smooth high-density wind. This, in turn, could be counteracted by the destruction of clumps in colliding-wind binaries (Pittard 2007). Another possibility is the effect of porosity (Owocki & Cohen 2006; Shaviv 1998), where at least some of the clumps are optically thick thereby allowing more radiation to escape through the inter-clump medium. A further effect is that hot, low-density regions may persist even in the radiative part of the colliding-wind region (especially on the leading edge of each arm, Pittard 2009). The low density will reduce the Razin and free-free effects, resulting in a spectral index closer to the observations.

We also recall that, in our model, we left out the free-free absorption and emission due to the density and temperature changes in the colliding-wind region. Pittard (2010) presents models that study this effect. His Fig. 14 shows that the free-free contribution from the colliding-wind region has an intrinsic spectral index which can be close to 0. However, including the free-free contribution of the winds results in a positive index. This therefore does not help in explaining the observed spectral index.

5.3 X-ray emission

In Sect. 3, we noted the anti-correlated behaviour of the X-ray and radio light curves. The observed X-ray flux is caused by thermal emission due to heating in the wind-collision zone and subsequent absorption in the stellar wind material. Both the heating and the absorption change with stellar phase and thus lead to variability, just as for the radio fluxes. But the difference with the radio emission is that the optical depths are very much smaller in X-rays. For single stars, thermal X-ray emission that is formed in shocks due to the radiative driving instability can be seen down to 1-2 stellar radii (Kramer et al. 2003). Radio emission, on the other hand, can only be seen from 100 stellar radii. In a massive colliding-wind binary such as Cyg OB2 No. 8A this means that the radio emission will come from outer parts of the contact discontinuity, but the X-ray emission will be formed close to the apex of the contact discontinuity.

Minimum flux occurs at periastron in this eccentric binary, where the distance between the two components is smallest. The winds have attained only a fraction of their terminal velocities before colliding, resulting in a lower post-shock temperature. In addition, with the reduced separation, the X-ray emitting plasma is buried in denser wind material. This leads to a stronger absorption of the X-rays emitted by the plasma heated by the colliding winds, which is clearly seen in the measured local absorption column densities (Fig. 2, upper panel). This results in a strong decrease of the emerging X-ray flux at this orbital phase, even though an increase of intrinsic X-ray flux (that scales with the emission measure) is expected. At apastron, the winds have higher velocities before colliding, resulting in a higher post-shock temperature, but the absorption by the wind material surrounding the apex of the wind interaction region (where the bulk of the X-rays are produced) is significantly reduced. The global shape of the X-ray light curve results in fact from the combined effect of varying post-shock temperature as a function of the separation, changing photoelectric absorption along the line of sight, and decreasing emission measure when going from periastron to apastron.

When we replace the standard model by the improved one, the contact discontinuity now wraps around the primary rather than the secondary. This changes the absolute position of the apex, but not the relative behaviour of the X-ray emission measure as a function of orbital phase. Phases of minimum and maximum emission measure will be the same for the two models.

Figure 2 (upper panel) shows that the observed absorption column changes with orbital phase. This puts an additional constraint on the models. To see if our modelling is qualitatively consistent with this, we consider the sight-line of the observer toward the apex (Fig. 9) and we discuss the absorption column along that line. At periastron (phase 0.0) in the standard model, the sight-line passes through the stronger wind while at apastron (phase 0.5) it passes through the weaker wind. Because of the position of the apex, the sight-line also passes closer to the star at periastron than at apastron (Fig. 9, upper panel). One should also take into account that the sight-line is not in the plane of the orbit, but at an angle of 58 (= ) above it. For the standard model, maximum absorption therefore occurs near periastron, minimum absorption near apastron, which is consistent with the observed behaviour shown in Fig. 2 (upper panel). For the improved model, on the contrary, the sight-line goes through the weaker wind at periastron. This might suggest a reversal of the phases for minimum and maximum absorption. However, the dominant effect is how close the sight-line passes to the relevant star. Figure 9 (lower panel) clearly shows that this still occurs at periastron. Both the standard and improved model are therefore qualitatively in agreement with the observed behaviour of the absorption column.

Pittard & Parkin (2010) present models for the thermal X-ray emission in O+O star binaries. They discuss Cyg OB2 No. 8A, although no specific model for this system is presented. One of their models is an eccentric binary which has a light curve qualitatively similar to the lower panel in Fig. 1: it shows a roughly constant flux for a large part of the orbit, followed by a rather sharp maximum which in turn is followed by the minimum. Also their model predicts changes in the hottest component of the emission, which is indeed observed, though caution should be used in the interpretation of this result as the error bars are quite large (see Sect. 3).

5.4 Future work

The improved model presented here is qualitatively successful in explaining a number of important features of the radio emission of the colliding-wind binary Cyg OB2 No. 8A. Nevertheless, it should be realized that some important physical ingredients are still missing. The major ones are including orbital motion and solving the equations of hydrodynamics. Orbital motion will have an important effect on the phases of minimum and maximum flux. Solving the hydrodynamics equations will allow us to correctly position the shocks and take into account the instabilities and inhomogeneities, as they influence both the thermal free-free emission and absorption, and the synchrotron emission. Information from hydrodynamics will also allow us to replace the assumed shock strength by its correct value. This will influence the relativistic electron momentum distribution. In addition, clumps that are optically thick at radio wavelengths will result in porosity, letting a larger fraction of the synchrotron emission escape. Finally, future models should handle both the X-ray and radio emission simultaneously, giving us a consistent picture of these colliding-wind systems.

6 Conclusions

We determined the radio light curve of Cyg OB2 No. 8A, based on new and archive VLA observations. We find phase-locked variability in the radio fluxes. This is contrary to expectations, as the synchrotron emission due to the colliding winds should be hidden by free-free absorption in this close binary. We also present an updated version of the X-ray light curve. The radio and X-ray light curves are nearly anti-correlated, which is another unexpected result.

Our model calculation for the radio flux shows that the synchrotron emission region extends considerably beyond the apex of the colliding-wind region, explaining why we do see phase-locked variability. We achieve this result without any need for clumping or porosity in the model. Our standard model - which uses the star and wind parameters proposed by De Becker et al. (2006) - does not correctly predict the phases of minimum and maximum. When we change the wind parameters (which are not well-determined) to make the secondary have the stronger wind, we find a considerable improvement in phases of minimum/maximum. A more detailed study of Cyg OB2 No. 8A is therefore needed to better determine the star and wind parameters of its binary components. To improve the phase agreement even further we will need to include orbital motion in the model.

The nearly anti-correlated behaviour of the X-ray and radio light curves is due to their very different formation regions. X-rays are formed quite close to the apex of the contact discontinuity, where the heating is strongest and they suffer only a limited amount of absorption. Radio emission, on the other hand, is strongly absorbed, so the radio emission we can detect is formed much further out along the contact discontinuity. Because the contact discontinuity is curved, this results in X-rays and radio having their minimum/maximum flux at very different phases.

The model values for the radio spectral index (+0.75 to +1.5) differ considerably from the observed value of 0.0. However, the model values are quite sensitive to the density of the thermal electrons, which are responsible for the Razin effect and the free-free absorption. A decrease of that density would allow a better agreement with the observations. Such a decrease could be due to lower density regions caused by hydrodynamical effects in the colliding-wind region or by clumping in a wind with a reduced mass-loss rate. Although we did not require clumping or porosity to explain the phase-locked variability, it is likely that we will need it to explain the observed spectral index. Because of the high sensitivity of the spectral index to clumping conditions, it could provide a valuable pathway to investigate the well-known problem of clumping in single-star winds. Future,improved models, applied to Cyg OB2 No. 8A and other non-thermal radio emitters should therefore be capable of determining limits on this clumping.

Acknowledgements

We thank Joan Vandekerckhove for his help with the reduction of the VLA data. We are grateful to the original observers of the VLA archive data used in this paper. We also thank Patrick Müller for discussions about the numerical code. We thank the referee for comments that helped to improve this paper. D. Volpi acknowledges funding by the Belgian Federal Science Policy Office (Belspo), under contract MO/33/024. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France and NASA's Astrophysics Data System Abstract Service.

References

Online Material

Appendix A: Modelling

Most of the calculations discussed here are limited to the two-dimensional (2D) plane of the binary orbit. Only in Sect. A.7 do we rotate the results of the 2D calculations into the three-dimensional (3D) simulation volume that is used to calculate the resulting flux.

A.1 Binary orbit

Figure A.1 shows how the two stars are positioned in the XY plane of the 3D simulation volume. The primary is at the centre of the volume and the position of the secondary is determined by the orbital phase. The angle of periastron () is measured from the intersection of the plane of the sky with the plane of the orbit. The sight-line towards the observer is in the YZ plane and the inclination angle i is the angle between the plane of the sky and the orbital plane. The secondary rotates clockwise in this figure (as seen from above). Periastron corresponds to phase 0.0.

A.2 Contact discontinuity

The position of the contact discontinuity at any given orbital phase is determined by the balance of ram pressure between the velocity components orthogonal to the contact discontinuity (see equations in Antokhin et al. 2004). These equations take into account that the terminal velocity of the winds may not have been reached yet before they collide. Specifically, in the Cyg OB2 No. 8A standard model presented in Sect. 5.1, the winds collide with (primary) and (secondary) at the apex at periastron. After the stellar wind material has passed through the shock, it then moves outward parallel to the contact discontinuity (there is also a small orthogonal velocity component, see Sect. A.5). We assume that the shocks and the contact discontinuity are sufficiently close to one another that we do not need to make a distinction in their geometric positions. On each side of the contact discontinuity, we assume the stellar wind to have a velocity law and a density derived from the mass conservation law, using the assumed velocity law and mass-loss rate.

A.3 Momentum distribution relativistic electrons

We determine the momentum distribution of the relativistic electrons for a grid of points along the two shocks. We start with a population of electrons generated at a given point on the shock itself. We then follow the time evolution of this population as it advects away and cools down. The momentum distribution uses a grid that is uniform in , between p₀ = 1.0 MeV/c (for the choice of this specific value, see Van Loo et al. 2004) and the upper limit (defined below, Eq. (A.3)) These limits are sufficiently wide to cover all of the radio emission we are interested in. We stop following the electrons when their momentum falls below p₀, or when they leave the simulation volume. (In our calculations we checked that the emission beyond the simulation volume is negligible.)

In the following, we refer mainly to Van Loo (2005, hereafter VL) as the source for the equations. Citations to the original papers can be found in that reference.

A.4 New relativistic electrons

At the shock, new relativistic electrons are put into the system. Taking into account the particle acceleration mechanism and the inverse Compton cooling, it can be shown that their momenta follow a modified power-law distribution (VL, Eqs. (5.5), (5.6)):

With (VL, Eqs. (5.1)-(3)):

with and e the mass and charge of the electron, the Thomson electron scattering cross section, c the speed of light, L_* the luminosity of the star and r the distance of the point on the shock to the star. Because there are two stars, we should calculate the contribution of both:

with r₁ the distance to the primary and r₂ the distance to the secondary. Also (VL, Eq. (3.7)):

where u₁ is the upstream velocity and u₂ the downstream one and the shock strength . Throughout this work we assume strong shocks with , but we note that higher values can be expected in radiative shocks.

The normalization factor is related to , the fraction of energy that gets transferred from the shock to the relativistic particles, by (VL, Eq. (5.10)):

where is the upstream density and we assume (Blandford & Eichler 1987; Eichler & Usov 1993).

The magnetic field at distance r to the relevant star is assumed to be given by (Weber & Davis 1967; (VL, Eq. (2.48))):

where B_* is the surface magnetic field (assumed to be 100 G) and is the equatorial rotational velocity. This equation is only valid at larger distances from the star, but we use it everywhere as the region close to the star will not be visible anyway (due to free-free absorption). This equation neglects any amplification of the field (such as described by, e.g. Bell 2004) or any reduction by magnetic reconnection. For each of the two shocks, the corresponding magnetic field is taken into account.

A.5 Advection and cooling

During a time step , the relativistic electrons advect away and cool mainly due to inverse Compton and adiabatic cooling (see Chen 1992; and VL, p. 96, for a discussion on other cooling mechanisms and why they are not important in the present situation).

The equation for cooling of a single relativistic electron can be written as:

where a p² gives the inverse Compton cooling (VL, Eq. (5.2)) and b p the adiabatic cooling (which can be derived from the fact that for adiabatic cooling and using ).

Integration of Eq. (A.8) gives:

where is the momentum at the previous time step. We can invert this to:

We will also need the following derivative:

We assume that only the momentum p changes significantly with time; for the other variables, we take their value mid-way between the start and the end of the time step . The density we need in Eq. (A.10) is the one between the shock and the contact discontinuity. We do not have detailed hydrodynamical models for this, so we simply approximate it by where is the density in the smooth wind at the position of the shock. We stress once again that in these radiative shocks, the post-shock density increases rapidly downstream, maybe by several orders of magnitude, rather than the value we use.

The above equations describe what happens to a single relativistic electron. The time evolution of the momentum distribution function over the time step is:

We stop following the electrons when their momentum falls below p₀, or when they leave the simulation volume.

In our simplifying assumptions the contact discontinuity and the shock are at the same geometrical position. Taken literally, this would result in a synchrotron emitting region that has zero thickness and zero synchrotron emission. In the true hydrodynamical situation, particles move through the shock (which is at some distance from the contact discontinuity) and then move slowly towards the contact discontinuity. At the same time, they are being advected outward with a velocity parallel to the contact discontinuity that is much larger than the corresponding orthogonal velocity. The thickness of the synchrotron emitting region is set by the combination of these two velocities. In the absence of hydrodynamical information, we use the following procedure to determine these velocities. The parallel component ( ) at the contact discontinuity is determined by projecting the smooth wind velocity vector on to the contact discontinuity. At other positions in the wind, we interpolate in the grid of as a function of the y-coordinate. The orthogonal component ( ) at the contact discontinuity is expected to be of the order of the orthogonal component of the material coming into the shock, with a minimum value corresponding to the thermal expansion of the gas. As a simplification, we set it exactly equal to the incoming orthogonal component, and we use the same type of interpolation as for the parallel velocity to determine it at other positions. We reverse the orthogonal velocity direction to have the gas moving away from the contact discontinuity. Note that this situation is ``inside out'' compared to the true hydrodynamical situation where the particles move towards the contact discontinuity. The thickness of the synchrotron emitting region is so small however that this inversion has negligible influence on our results.

A.6 Emissivity

The equation for the emissivity is given by (VL, Eq. (B.1)):

with (VL, Eq. (2.51), (2.50)):

We use Eq. (A.7) for the magnetic field. Note that this is an approximation, as the shock may compress the magnetic field as well, leading to higher values than given by Eq. (A.7). The expression for is evaluated at position r and frequency . The f function corrects for the Razin effect and K_5/3 is the modified Bessel function. To calculate , the number density of the thermal electrons (i.e. those in the non-relativistic plasma), we use:

where is the smooth wind density and the mass of a hydrogen atom. The factor 1.2/1.4 takes into account a fully ionized hydrogen+helium solar composition. In these radiative shocks the post-shock density may increase by several orders of magnitude rather than by , which will influence the value of .

To reduce the computing time we set in the factors of Eq. (A.15) and we also use:

The F function is precalculated on a grid of x-values; the evaluation of the function then reduces to an interpolation.

A.7 Map into 3D grid

Using the equations from the previous sections, we can make a 2D map of the synchrotron emission (part of which is shown in Fig. 6). Note that the synchrotron emission extends on both sides of the contact discontinuity, because there is a shock on either side. Based on the cylindrical symmetry around the line connecting the two stars (at the given orbital phase), we then ``rotate'' this 2D map into the 3D simulation volume.

It is important that during this procedure we do not lose part of the emission due to interpolation or reduced resolution. To achieve this, we use volume-integrated emissivities. While doing the 2D calculation we store the volume integrated emissivity defined by:

The volume used in this definition is determined as follows. We consider the 2D surface covered by (1) the step in the coordinate along the contact discontinuity and (2) the distance travelled by the relativistic electrons between two consecutive time steps. We then rotate this surface around the line connecting the two stars over a small angle , thereby defining the emissivity volume .

When we have finished the full 2D calculation, we rotate the 2D plane into the 3D simulation volume, using steps of angle (as defined above). For each of the emissivity volumes defined above, and for a given rotation angle, we determine what fraction of the rotated volume overlaps with the cells in the 3D simulation volume. This is done by subdividing each dimension of the volume into three. For each of the resulting 3³ points we determine in which 3D simulation cell they end up after rotation. A running sum for that cell is then updated with the fraction 3^-3 of the volume integrated emissivity of the rotated volume. At the end of this procedure, the accumulated emission in each 3D simulation cell is divided by its volume, thus obtaining for that cell.

In this procedure we lose some resolution, specifically around the apex, where the 2D resolution is quite high but where the 3D resolution is much lower. However, because of the use of volume-integrated quantities, we do not lose any of the synchrotron emission. Furthermore, the apex of the wind is well hidden by the free-free absorption and will therefore not contribute to the resulting radio flux at the wavelengths we consider. If our model were to be applied to X-rays, optical emission or radio emission at much shorter wavelengths, a higher-resolution 3D grid would be required.

In the 3D model. we also include the thermal free-free opacity and emissivity, which is due to the ionized wind material. To calculate these we use the Wright & Barlow (1975) equations. The winds of both stars are assumed to be smooth. For the Gaunt factor at frequency we use the equation from Leitherer & Robert (1991):

where is the thermal electron temperature in the wind and Z the charge of the ions. The stars themselves are introduced into the simulation volume as high opacity spheres. We then solve the radiative transfer equation in Cartesian coordinates along the line of sight. We use Adam's (Adam 1990) finite volume method because of its simplicity. By repeating the whole procedure for a number of orbital phases and wavelengths, we calculate the radio light curve at these wavelengths and determine how the fluxes and spectral index change with phase.

Appendix B: Data reduction

The data reduction is accomplished using the Astronomical Image Processing System (AIPS), developed by the NRAO. We follow the standard procedures of antenna gain calibration, absolute flux calibration, imaging and deconvolution. Where possible, we apply self-calibration to the observations (see notes to Table B.1). We measure the fluxes and error bars by fitting elliptical Gaussians to the sources on the images. The error bars listed in Table B.1 include not only the rms noise in the map, but also an estimate of the systematic errors that were evaluated using a jackknife technique. This technique drops part of the observed data and redetermines the fluxes, giving some indication of systematic errors that can be present. Upper limits are three times the uncertainty as derived above. For details of the reduction, we refer to the previous papers in this series (Blomme et al. 2007; Van Loo et al. 2008; Blomme et al. 2005). We exclude from Table B.1 those observations with upper limits higher than 6 mJy (which corresponds to about 3 times the highest detection at all wavelengths).

A comparison with values previously published in the literature shows that in most cases the error bars overlap with ours. There are just a few exceptions, which we discuss further. Bieging et al. (1989) list only an upper limit for the AC116 2 cm observation, while we find a detection (though marginal at just above 3 sigma). Our detection is quite consistent in all variants we consider in the jackknife test and is at the exact position of Cyg OB2 No. 8A, so we consider it a detection. For the AS483 6 cm, Waldron et al. (1998) find a value of  mJy, which is 20% lower than our value. A possible explanation is that they did not correct for the decreasing sensitivity of the primary beam. This effect is important when the target is not in the centre of the image, and is approximately 20% for the position of No. 8A on the image. Puls et al. (2006) list an upper limit of 0.54 mJy for the AS786 6 cm observation, but we clearly detect Cyg OB2 No. 8A on the image we made. Bieging et al. (1989) find a slightly lower value of  mJy for the AC116 1984-12-21, 20 cm observation. We have no explanation for these last two discrepancies.

We also take the opportunity to correct some clerical errors in the literature: Bieging et al. (1989) date the CHUR observations as 1980-Mar-22, but it should be 1980-May-23. The Waldron et al. (1998) 1991 observation was made at 3.6 cm, not at 6 cm and the 1992 observation was made in January, not June.

Table B.1:   Reduction of the Cyg OB2 No. 8A VLA data.

Footnotes

... 8A

Partly based on observations with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).

...

Appendices are only available in electronic form at http://www.aanda.org

... NRAO

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

... spectrum

The spectral index for a pure power-law momentum distribution p^-n is -(n-1)/2, with . As we use , one would expect a -0.5 spectral index, but Pittard et al. (2006, their Fig. 4) show that the inclusion of inverse Compton cooling changes this index to a value closer to -1.0.

... Compton

The inverse Compton cooling should, in principle, use the Klein-Nishina scattering equation instead of the Thomson one. To test the importance of this, we adapted our code to use the Klein-Nishina equation, following the description in Schlickeiser & Ruppel (2010). The input photons are assumed to come from a graybody with a temperature of 40 000 K. The Klein-Nishina equation leads to a somewhat more complicated expression than Eq. (A.11) which can no longer be inverted analytically. We therefore do the time integration numerically, as well as the determination of . The resulting fluxes show fractional differences less than 10^-4. We therefore decided to continue using the simpler Thomson scattering in the model.

Copyright ESO 2010

All Tables

Table 1:   Parameters of Cyg OB2 No. 8A used in this work.

Table 2:   Two new X-ray observations of Cyg OB2 No. 8A.

Table 3:   Parameters for the two new EPIC spectra of Cyg OB2 No. 8A.

Table B.1:   Reduction of the Cyg OB2 No. 8A VLA data.

All Figures

---

Copyright ESO 2010