Radio emission models of colliding wind binaries

# 3D models of radiatively driven colliding winds in massive O+O star binaries - II. Thermal radio to sub-mm emission

## Abstract

In this work the thermal emission over cm to sub-mm wavelengths from the winds in short-period O+O-star binaries is investigated (potential non-thermal emission is presently ignored). The calculations are based on three-dimensional hydrodynamical models which incorporate gravity, the driving of the winds, orbital motion of the stars, and radiative cooling of the shocked plasma. The thermal emission arises from the stellar winds and the region where they collide. We investigate the flux and spectrum from a variety of models as a function of orbital phase and orientation of the observer, and compare to the single star case. The emission from the wind-wind collision region (WCR) is strongly dependent on its density and temperature, being optically thick in radiative systems, and optically thin in adiabatic systems. The flux from systems where the WCR is highly radiative, as investigated for the first time in this work, can be over an order of magnitude greater than the combined flux from identically typed single stars. This excess occurs over a broad range of wavelengths from cm to sub-mm. In contrast, when the WCR is largely adiabatic, a significant excess in the thermal flux occurs only below 100 GHz.

In circular systems with (near) identical stars the observed variability in synthetic lightcurves is typically less than a factor of 2. Eccentric systems may show order of magnitude or greater flux variability, especially if the plasma in the WCR transitions from an adiabatic to a radiative state and vice-versa - in such cases the flux can display significant hysteresis with stellar separation. We further demonstrate that clumping can affect the variability of radio lightcurves.

We investigate the spectral index of the emission, and often find indices steeper than . Synthetic images display a variety of morphologies, with the emission sometimes resembling an intertwined “double-helix”. We conclude by comparing our results to observations. The predictions made in this paper are of interest to future observations with the next generation of radio and sub-mm telescopes, including the EVLA, e-MERLIN, ALMA, and the SKA, and future upgrades to the VLBA.

###### keywords:
shock waves – stars: binaries: general – stars: early-type – stars: mass loss – stars: winds, outflows – radio continuum: stars
12

## 1 Introduction

Luminous stars of OB and Wolf-Rayet (WR) spectral types drive fast, dense, ionized winds with mass-loss rates that are high enough to significantly affect their evolution and their surroundings. The winds emit free-free radiation over a wide range of wavelengths, from cm emission in the radio band, to m emission in the near-infra-red. Since the bremsstrahlung opacity is proportional to the square of the wavelength, the observed emission originates above a characteristic radius that increases with wavelength (typical values are at 20 cm (1.5 GHz), and at 1 mm (300 GHz)). At wavelengths where the wind is optically thick, the thermal emission can be used to determine stellar mass-loss rates. This has mostly been done using radio data (e.g. Abbott et al. 1980, 1981; Bieging, Abbott & Churchwell 1989; Leitherer, Chapman & Koribalski 1995; Contreras et al. 1996; Scuderi et al. 1998; Benaglia, Cappa & Koribalski 2001; Schnerr et al. 2007), since the wind can be assumed to be at terminal velocity, and the interpretation of the flux is not strongly dependent on details of the ionization conditions (Wright & Barlow 1975; Panagia & Felli 1975). Wind asphericity typically does not alter radio derived mass-loss rates by more than a factor of two (Schmid-Burgk 1982).

The free-free flux is expected to behave like , where the spectral index at cm wavelengths (Wright & Barlow 1975). However, the spectral index of “thermal” sources3 at shorter (mm and far-infrared) wavelengths increases to around to (e.g. Williams et al. 1990; Leitherer & Robert 1991; Altenhoff, Thum & Wendker 1994; Nugis, Crowther & Willis 1998). Several possible causes for this steepening have been suggested, including a decrease of the wind temperature or the ionization state with radius or an increase of the wind speed with radius (Schmid-Burgk 1982; Leitherer & Robert 1991). However, since temperature gradients only weakly affect the flux through the gaunt factor, while the wind speed changes relatively little between the cm and mm emitting regions, these effects are unable to explain the difference between theoretical and observed spectral indices between cm and mm wavelengths. Schmid-Burgk (1982) also noted that deviations from the expected behaviour “should not easily be blamed” on non-spherical winds. This led both Williams et al. (1990) and Leitherer & Robert (1991) to suggest a gradient in ionization as the most likely cause.

Another possibility involves a radial variation of the degree of clumping within the wind. For instance, an excess of mm emission over that expected from scaling the radio emission to shorter wavelengths can be caused by a decline in the degree of clumping between the regions of the wind emitting at mm and cm wavelengths (Runacres & Blomme 1996; Blomme et al. 2002). Nugis et al. (1998) have also argued for models which combine clumping and a highly ionized outer wind region. A recent review of small-scale structure in massive star winds can be found in Puls, Vink & Najarro (2008).

At very high frequencies, the photospheric emission, which is rapidly rising, begins to dominate the emission from the wind. This transition occurs at roughly the same frequency as the whole wind becomes optically thin. The intrinsic thermal spectrum from the wind only then flattens, as the emission then arises from essentially a blackbody of constant projected linear radius at all frequencies greater than the critical value given by Eq. 12 of Wright & Barlow (1975). However, the dominance of the photospheric emission means that observationally the emission retains a strong wavelength dependence.

The presence of a companion star may also change the observed thermal flux, its spectral index, and the inferred mass-loss rates. First, there are now two stars, rather than one, and the winds of both will contribute to the observed emission. Moreover, the interaction of the winds may also produce significant free-free emission. Stevens (1995) investigated the total thermal flux from a wide colliding winds binary (CWB) as a function of the wind momentum ratio, and discovered that it could be about 50 per cent higher than that expected from the more massive wind alone. In more recent work, Pittard et al. (2006) investigated the emission from an adiabatic WCR, where the plasma in the WCR remains at high temperature as it flows out of the system, as a function of the stellar separation, . If the WCR remains optically thin, the thermal flux scales as . Moreover, since optically thin emission scales as , the flux from the WCR may dominate the thermal flux from the system, particularly at cm wavelengths4. However, there has been no investigation to date of the thermal emission from systems where the WCR is highly radiative. This situation clearly needs to be addressed.

In this paper we investigate the thermal radio to sub-mm emission from short period O+O-star binaries computed from the three-dimensional hydrodynamical models described in Pittard (2009a, hereafter Paper I). We aim to illuminate the relative contribution of the WCR to the total free-free flux, the effect that this may have on the observed spectral index, and the size and nature of orbital related variations, for a variety of models where the behaviour of the plasma in the WCR is either highly radiative or largely adiabatic.

The layout of this paper is as follows. We briefly describe in Section 2 the hydrodynamical models and in Section 3 our method of calculating the thermal cm to sub-mm emission. In Section 4 we present our results, focussing first on calculations of the emission from the wind of a single star, and then on calculations based on our binary models. In Section 5 we compare our findings against observations, and in Section 6 we summarize and conclude this work.

## 2 Summary of the hydrodynamical models

The radio calculations in this paper are based on the three-dimensional hydrodynamical models described in Paper I. The models incorporate the radiative driving of the stellar winds, gravity, orbital effects, and cooling, and are summarized in Tables 1 and 2. The models were not designed to simulate particular systems, but were aimed at mimicking much of the interesting phenomena that can occur.

Model cwb1 is of an O6V+O6V system with a circular orbit and a period of 3 days. The WCR is highly radiative, and significantly distorted by orbital effects, showing strong aberration and downstream curvature. Model cwb1 is similar to DH Cep (see Linder et al. 2007, and references therein), HD 165052 (Arias et al. 2002; Linder et al. 2007), and HD 159176 (De Becker et al. 2004b; Linder et al. 2007), all of which have near identical main-sequence stars of spectral type O6O7, and circular or near-circular orbits with periods near 3 days. Hence the hydrodynamics of, and emission from, the WCR in model cwb1 should be a reasonable approximation to these systems.

Model cwb2 has identical parameters to model cwb1 except for an increase in the orbital period to 10 days. The winds now have more room to accelerate before they collide, and the postshock gas remains largely adiabatic as it flows out of the system. Both the aberration angle and the downstream curvature of the WCR are lessened relative to model cwb1. Model cwb2 is similar to HD 93161A, an O8V + O9V system with a circular orbit and an orbital period of 8.566 days (Nazé et al. 2005), albeit with slightly more massive stars and powerful winds. Another system which is not too dissimilar is Plaskett’s star (HD 47129, Linder et al. 2006, 2008), though this object contains stars which have evolved off the main sequence.

Model cwb3 examines the interaction of unequal winds in a hypothetical O6V + O8V binary. The stars have the same separation as those in model cwb2, but a slightly longer orbital period. The stronger wind from the primary star collides at higher speeds than the slower secondary wind, resulting in postshock plasma which is hotter and more adiabatic.

Model cwb4 investigates the effect of an eccentric orbit, which takes the stars through a separation of (i.e. the separations of the stars in the circular orbits of models cwb1 and cwb2). The WCR is radiative at periastron and adiabatic at apastron, and its aberration and downstream curvature are phase dependent. A surprising finding from Paper I is that dense cold clumps formed in the WCR at periastron still persist near the apex of the WCR at apastron. This is because the clumps have relatively high inertia, and require a significant fraction of the orbtial period to be accelerated out of the system. Some well-known O+O binaries with eccentric orbits include (in order of increasing orbital period) HD 152248 (; Sana et al. 2004), HD 93205 (; Morrell et al. 2001), HD 93403 (; Rauw et al. 2002), Cyg OB2#8A (; De Becker, Rauw & Manfroid 2004; De Becker et al. 2006) and  Orionis (; Bagnuolo et al. 2001).

## 3 Modelling the radio emission from CWBs

We use the model developed by Dougherty et al. (2003) and Pittard et al. (2006) to determine the thermal radio to sub-mm emission from each of our colliding winds simulations. Detailed notes on our method can be found in these papers, but the salient points are summarized below.

To calculate the radio emission we read our hydrodynamical models into a radiative transfer ray-tracing code, and calculate appropriate emission and absorption coefficients for each cell. A synthetic image on the plane of the sky is then generated by solving the radiative transfer equation along suitable lines of sight through the grid. The thermal emission () and absorption () coefficients at frequency are given by Rybicki & Lightman (1979) as

 εffν = 6.8×10−38Z2neniT−1/2e−hν/kTgff, (1) αffν = 3.7×108T−1/2Z2neniν−3(1−e−hν/kT)gff, (2)

where is the mean ionic charge, and are electron and ion number densities, is the temperature of the gas, and is a velocity averaged Gaunt factor (Hummer 1988). All quantities are in cgs units. The densities and are determined from the cell density, composition and ionization, where the ionization is specified as a function of cell temperature. The appropriate ionization for a particular temperature regime is determined and the absorption and emission coefficients are then evaluated. The unshocked winds are assumed to be isothermal at  K, and each model is assumed to be at a distance of 1 kpc.

In the majority of the models in this work we assume that the winds are smooth (i.e. not clumped). This enables the characteristic radius of emission and the surface to better fit within our hydrodynamical grids, at the expense of underestimating the emission compared to real systems where the winds are often strongly clumped. Our assumption of smooth winds also removes the need to specify the radial stratification of the clump volume filling factor, , which remains poorly known (see e.g. Puls et al. 2008). Table 3 lists the characteristic radius, , and flux density, , of free-free emission from our model O6V star at various frequencies between 5 and 2000 GHz (corresponding to wavelengths between 6 cm and 150 m). The radius of optical depth unity, .

All of the hydrodynamical grids in our models are large enough to contain the characteristic radius for free-free emission at GHz, while models cwb2 and cwb3 are also able to contain the characteristic radius for emission at GHz (compare Table 3 with Table 3 in Paper I). This enables us to explore the thermal radio emission from these systems, in addition to the X-ray emission explored in a previous paper (Pittard 2009b, hereafter Paper III).

Clumping increases the emission and absorption coefficients in Eqs. 1 and 2 by a factor of , where is the volume filling factor of the clumps, and increases the characteristic radius by a factor , and the flux density by a factor . This means that the flux from a model with a given and , if interpreted as arising from a wind with no clumping, gives (incorrectly) the mass-loss rate (note that one can also take an alternative approach: a particular observed flux can be translated into a mass-loss rate assuming that there is no clumping, or a mass-loss rate , assuming that there is clumping). However, even when the stellar winds are clumpy, the WCR may be considerably less so. This is expected to be the case when the postshock gas behaves largely adiabatically. In such situations the destruction time of the clumps can be significantly shorter than the flow time of the shocked gas out of the system, and the WCR, though highly turbulent, is noticeably smoother than the stellar winds (see Pittard 2007). Therefore, in previous work (Dougherty et al. 2003; Pittard et al. 2006; Pittard & Dougherty 2006) we have assumed that the winds are clumpy and the WCR is smooth(er).

In this work we also investigate the emission assuming clumpy winds with , but a WCR which is smoother to varying degrees i.e. . This provides more realistic fluxes from the unshocked winds and insight into how clumping in different regions of the system affects the observed lightcurves and spectra. In these models the coefficients in Eqs. 1 and 2 are multiplied by a factor of () for cells in the unshocked winds (WCR), respectively.

## 4 Results

### 4.1 Emission from the wind of a single star

We first investigate the cm to infrared free-free emission from the wind of a single star, namely from the O6V star used in our binary simulations. The wind density as a function of radial distance is obtained from our solution to the modified CAK equations (Castor, Abbott & Klein 1975; Pauldrach et al. 1986). At high frequencies the emission arises from within the acceleration zone of the wind, and it is no longer possible to use, e.g., Eq. 8 in Wright & Barlow (1975) to evaluate the flux. We therefore numerically integrate the flux as a function of impact parameter, following the method of Wright & Barlow (1975). The wind is assumed to begin at the point where the outflow velocity is equal to the isothermal sound speed, which is in our models. As in the binary models which follow we assume that the wind is smooth, isothermal and maintains a constant ionization: these assumptions allow us to isolate the effect of the wind acceleration on the resulting emission. We also include the photospheric flux, which is modelled as a black body for simplicity. The photospheric flux is attenuated at wavelengths where the wind is optically thick.

Fig. 1(a) shows the resulting spectrum obtained from the wind and photosphere of our model O6V star. At GHz frequencies (cm wavelengths) the total emission is dominated by the wind and arises from large distances where the wind has very nearly reached its terminal speed. The spectral index (Fig. 1b) is near , as expected. The emission starts to probe the acceleration region of the wind at  GHz, where a steepening of the spectrum becomes apparent. It is from this point that the flux increasingly exceeds that from a calculation where the wind is forced to emanate from the star at its terminal speed (i.e. instantaneous acceleration). As shown in Fig. 1(a), such a terminal speed wind becomes optically thin at frequencies above 700 GHz (Eq. 12 of Wright & Barlow 1975), with declining flux thereafter (see also Bertout et al. 1985). In contrast, the sharp increase in the wind density of the CAK model as the photosphere is approached results in a steep rise in the thermal flux. At frequencies around  Hz (), the spectral index of the emission from the wind steepens to reach a maximum of . This steepening is due purely to the observed emission occuring from regions where the wind is rapidly accelerating.

The wind also starts to become optically thin at  Hz (), and there is a steep rise at higher frequencies in the spectral index of the total emission as the photospheric emission becomes dominant at  Hz (; for stars with denser winds, the emission from the wind may dominate the photospheric emission until m wavelengths). The maximum spectral index for thermal free-free emission is , the value for a blackbody source where the flux , as given by the Rayleight-Jeans approximation to the Planck function. This is indeed the case in a calculation where the wind and photospheric temperatures are equal. However, when there is a sharp change between the wind and photospheric temperatures, such as in our simple model, spectral indices greater than +2.0 can be obtained. This is apparent in Fig. 1(b), where the spectral index actually reaches at  Hz, due to the steeply rising photospheric flux as the absorption of the overlying wind declines.5 The photospheric emission from our model O6V star peaks at  Hz (), and subsequently declines.

For a smooth, spherical, isothermal wind of constant ionization, the spectral index, , is related to the acceleration index, (where ), by (e.g. Leitherer & Robert 1991)

 α=4ϵ+1.82ϵ+3. (3)

Note that is a local approximation (assumed valid over the formation region corresponding to a particular wavelength). The value of is not valid over the whole distance range, in the same way that the value of in the “-law” (see below) is.

To obtain requires that . This is easily achieved in accelerating winds. Fig. 1(c) shows as a function of radius for the CAK wind profile and standard “-law” parameterizations (where

 v(r)=v0+(v∞−v0)(1−R∗/r)β, (4)

and is the initial wind velocity at the stellar radius, assumed here to be ). The CAK velocity profile is reasonably well approximated by a velocity law with . Winds with slower acceleration have larger values of , and vice-versa. Winds with large have smaller values of close to the star, but larger values further from the star where their winds continue to accelerate. Values of greater than 100 are achieved in the innermost regions of strongly accelerating winds. In the CAK model at .

The maximum frequency which can be explored in the 3D hydrodynamical models is set by the necessity of sufficiently resolving the acceleration region where the characteristic radius of emission occurs. Models cwb1 and cwb4 have a numerical resolution of . In models cwb2 and cwb3 the numerical resolution is . Therefore, it is only possible to examine frequencies up to  THz, since the optical depth for an observer at infinity to is less than unity at frequencies greater than this. In the following we restrict our investigation to  THz (m).

### 4.2 Emission from binary systems

#### Model cwb1

The hydrodynamical grid used to compute model cwb1 contains within its volume the characteristic radius of emission from the O6V wind at  GHz. The optical depth unity surface, which is slightly smaller than , fits completely inside the grid for  GHz. Therefore, we limit our study of the free-free intensity images and lightcurves from model cwb1 to  GHz, but show spectra down to 8 GHz for ease of comparison with our other models.

Fig. 2 shows intensity images of the free-free emission from model cwb1 at 43, 100, 250, and 1000 GHz, for an observer located directly above the orbital plane (). It is immediately apparent that the WCR dominates the emission, even at a frequency as high as 1 THz. The morphology of the emission from the WCR reflects its underlying structure and clumpiness, caused by its high susceptibility to dynamical instabilities (see Paper I for a detailed discussion of the hydrodynamics). Near the apex of the WCR the projected emission from different inhomogeneities merges together, but further downstream individual clumps and bowshocks can be identified. It is clear that a great deal of potential emission is not captured due to the limited physical extent of the hydrodynamical grid used in the model. For obvious reasons it is impossible to know just how much flux is “lost” in this way as dense clumps in the WCR flow out through the grid boundaries, but it is clear that this problem is worse at lower frequencies where structure in the WCR remains optically thick for larger downstream distances6. Nevertheless, these images, and the lightcurves and spectra which follow, provide important new insight into the emission at cm and sub-mm wavelengths from CWBs.

The spatial scale of the images is 1.12 by 1.12 mas, which is too small to be resolved with the EVLA (angular resolution of 4 mas at 50 GHz) or ALMA (highest angular resolution expected to be about 10 mas). In contrast, the VLBA has the necessary spatial resolution, but even with all current upgrade plans may struggle with the necessary sensitivity to obtain images of this system within a reasonable integration time (e.g. 10 hrs) to avoid smearing of the essential features due to orbital motion. In 10 hrs, with 256 Mbit/s bandwidth at 43 GHz, the VLBA + GBT + phased-VLA obtains a sensitivity of 0.029 mJy/beam. Bandwidth improvements to 4 Gbit/s will increase the sensitivity by a factor of 4. A signal to noise ratio (SNR) of 5 (the minimum needed to provide imaging constraints) in 10 hrs thus requires a flux of 0.0725 mJy/beam. The spatial resolution at 43 GHz is  mas, so each arm of the WCR falls within a beam (see Fig. 2), which is roughly 50 per cent of the total flux, or 0.02 mJy. Hence a 10 hr exposure is insufficient to obtain a SNR of 5 - instead  hrs is required. Having said this, it has already been noted that there is a significant loss of flux at such frequencies from model cwb1 due to the finite size of the hydrodynamical grid, so the future detection of such a system with the VLBA is not totally beyond the realm of possibility (e.g. if the flux from model cwb1 at 43 GHz was actually higher). Systems with stronger winds (higher mass-loss rates) than assumed for model cwb1 may further improve the possibility of such imaging.

Fig. 3 shows the free-free emission lightcurves at 43 GHz (top panels) to 1000 GHz (bottom panels), as a function of the inclination angle, , of the system. The left panels show the combined emission from the unshocked winds and the WCR, while the middle and right panels display the unshocked and shocked components, respectively. The stars are aligned along the lines-of-sight at phases 0.0 and 0.5, and are at quadrature at phases 0.25 and 0.75. The variations in flux are entirely due to changes in the occultation and circumstellar absorption, as the intrinsic emission is constant. Note that the presence of the companion star breaks the spherical symmetry of the unshocked winds. For example, because of the twin radiation fields, the winds accelerate more slowly towards each other, and more quickly away from each other. There are also regions of lower velocity wind in the “shadows” behind each star (see Paper I for further details).

The lightcurves of the WCR emission (right panels in Fig. 3) would display dual symmetry about phases corresponding to both quadrature (0.25, 0.75) and conjunction (0, 0.5) if there were no orbital effects on the WCR (i.e. no aberration or coriolis induced curvature - cf. Pittard & Stevens 1997), because of the equal winds and constant stellar separation. Orbital effects break the symmetry about quadrature, so that in the case of identical winds a double periodicity in the lightcurve is expected. However, dynamical instabilities which form in the WCR can break this symmetry too, as revealed by close examination of Fig. 3.

At the lower frequencies the flux from the unshocked winds is typically only a small fraction of the expected combined flux from the winds (see the horizontal red lines across the plots in the middle column of Fig. 3). For instance, at 43 GHz the thermal flux from a single smooth (unclumped) wind of the same parameters as used in model cwb1 is 0.017 mJy, while the combined unshocked winds contribution in Fig. 3(b) is mJy. This disparity reduces at higher frequencies - at 250 GHz, the two unshocked winds contribute mJy, whereas the expected single star value is 0.046 mJy - and is eliminated by 1000 GHz.

The lack of observed flux from the unshocked winds as shown in Figs. 3(b) and (e) is likely due to the removal of a large fraction of the unshocked material from both of the winds by the presence of the WCR. In contrast, at the highest frequencies considered (250 and 1000 GHz), the characteristic radius of emission moves so close to the surface of the stars (particularly at 1000 GHz) that the emission from this region is not strongly disturbed by the presence of the WCR at larger distances. A small additional effect at the lowest frequencies arises from the limited extent of the computational grid compared to the characteristic radius of emission, which causes an underestimate of the flux. This can be seen in the spectra in Fig. 4, where a distinct downturn in the unshocked (“ff winds”) flux is visible at the lower frequencies investigated. However, at higher frequencies between  GHz the removal of unshocked wind material by the WCR is the main cause of the “missing” flux.

At 1000 GHz, the combined flux from the unshocked winds is mJy, which actually slightly exceeds twice the theoretical value of 0.095 mJy for our single star. This reflects the fact that the characteristic radius of emission is now deep within the acceleration zones of the winds, while the theoretical flux is estimated assuming that the entire wind is moving at its terminal speed. This excess over the “theoretical” value was also seen in the single-star calculations presented in the previous section.

The thermal emission from the unshocked winds is generally pretty steady with orbital phase. Slight maxima occur near to quadrature, with minima near to conjunction, as expected. The 1000 GHz lightcurve is not so smooth in appearance - its more “noisy” behaviour is a consequence of lines of sight through the WCR becoming optically thin in some directions, but not in others. The density of the thin layer of cooled postshock gas is near the apex of the WCR. The linear absorption coefficient of this gas at 1000 GHz is , while the thickness of this layer is typically , giving an optical depth of approximately unity. Hence we see the emitting regions of the WCR and the unshocked winds near the stellar photospheres through a semi-opaque, semi-porous WCR.

The free-free emission from the WCR, on the other hand, displays much greater phase variability. It has already been noted that the flux from the WCR dominates the emission from the unshocked winds. Fig. 3 shows that the WCR flux is around greater at  GHz, though the true value is likely to be somewhat higher. Future calculations with a larger grid are necessary to accurately determine this ratio. The lightcurves of the emission from the WCR typically display two minima per orbital cycle, with a phase-lead with respect to conjunction of the stars which appears to increase with frequency. Whereas the X-ray lightcurves of model cwb1 in Paper III produce reasonably sharp maxima, the lightcurves in the left column of Fig. 3 display much broader maxima.

An explanation for the behaviour of the lightcurve of the emission from the WCR now follows. Figs. 2 and 3 show that at radio to sub-mm frequencies, the most intense emission from model cwb1 occurs from the cold, dense sheet of post-shock material of cooled gas within the WCR. The intrinsic emissivity of this region is very high, since the free-free emission coefficient , and the density in the WCR is orders of magnitude greater than in the unshocked winds. The WCR is also optically thick, and, because of the downstream curvature induced by the coriolis effect, it displays a larger projected surface to the observer, and hence a larger flux, when the stars are at quadrature than at conjunction (see Fig. 2). The amplitude of the variability decreases slightly with frequency as the optical depth through the WCR declines.

Fig. 4 shows radio-to-sub-mm spectra as a function of orbital phase for an observer in the orbital plane (panels a and b), and directly above/below it (panel c). In all cases the optically thick emission from the WCR dominates the observed flux (remember that there is likely a significant loss of flux at the lowest frequencies displayed). However, the emission from the WCR is less dominant over the unshocked winds for an observer above or below the orbital plane (compare panel c to panels a or b), since the orientation of the layer of cold dense gas to the observer is less favourable (which reduces the flux observed from it), while that of the winds is more favourable (both winds are seen side-to-side on the sky).

The spectrum of the total thermal flux shows considerable curvature for an observer in the orbital plane (panels a and b), but significantly less curvature for an observer with . Fig. 4(d)-(f) shows the frequency dependence of the spectral index from these spectra. Between 100 and 250 GHz the spectral index of the unshocked winds, , is , , and , for viewing angles of (, ), (, ), and (), respectively. These are steeper than the canonical for optically thick, isothermal, terminal speed winds, reflecting the fact that the winds in our model are accelerating and asymmetrical. increases towards 1000 GHz as the emission probes deeper into the acceleration region of the winds. At 1000 GHz, , , and for the same orientations. The variation of with viewing angle clearly shows that the presence of the companion star has a non-negligible impact on this parameter, with large variations in between conjunction and quadrature.

The reduction in the spectal index at the highest frequencies of the emission from the WCR, , is an indication that parts of the WCR are becoming optically thin. Between 100 and 250 GHz, , , and for the three viewing angles noted above. At 1000 GHz, , , and , respectively. For the total free-free emission between 100 and 250 GHz, , , and , while at 1000 GHz, , , and . The less dominant contribution of the WCR at results in the least variation in with frequency (see Fig. 4(d)-(f)).

Fig. 5(a) shows the variation of the total thermal emission for an observer in the orbital plane between conjunction and quadrature. The emission at quadrature is dominant over the emission at conjunction up to  GHz. Between 500 and 1000 GHz the emission at these phases is almost identical. Note that in this frequency range the maxima and minima in the lightcurves no longer occur at quadrature and conjunction. It is not possible to state whether this behaviour continues beyond 1000 GHz, for the reasons mentioned earlier in Sec. 4.1.

#### Model cwb2

Fig. 6 shows intensity images of the free-free emission from model cwb2 at 15, 43, 100, 250, and 1000 GHz, for an observer located directly above the orbital plane (). The curved shape of the WCR is clearly seen. The WCR dominates the emission at the lower frequencies ( GHz), while the unshocked winds dominate the emission at higher frequencies. At 15 GHz the WCR is very bright in comparison to the winds. The brightest part of the WCR is its apex, between the stars. Further downstream the gas at the trailing edge of the WCR is brighter than at the leading edge, reflecting the higher densities there (see Paper I). The angular scale of all of the images in Figs. 6-8 is 2.7 by 2.7 mas, which is again too small to be resolved with EVLA/ALMA but may be resolved with the VLBA.

Intensity maps from model cwb2 for an observer in the orbital plane are shown in Fig. 7. The features in these images depend on the relative orientation and the degree of emission/absorption of the winds and the WCR, and display significant changes as a function of frequency and orbital phase. At conjunction only the foreground wind is seen, since this blocks the emission from the wind of the more distant star. The surface brightness of the foreground wind is less than that of the WCR at 15 GHz, so there is a “hole” in the emission from the WCR at conjunction. In contrast, the surface brightness of the unshocked winds is greater than that of the WCR when  GHz. The sharp, curved edges to the emission from the WCR near the borders of the images is due to limb-brightening. If the hydrodynamical grid were large enough to contain several spirals of the WCR around the stars one would expect to see an expanding series of progressively fainter, concentric, edges. The unshocked winds emission at 1000 GHz arises from a region very close to the surface of the stars.

Some very interesting features arise in the synthetic images as the phase of the orbit advances. From the observer’s point-of-view, the foreground star slides to the right, which at  GHz exposes more of the brighter WCR to the left by phase 0.125. At higher frequencies the winds are brighter than the WCR. At phase 0.25 (quadrature) the apex of the WCR is seen side on, with the unshocked winds to either side. The brightest part of the WCR is between the stars as expected, and is visible in these images even at 1 THz. At low frequencies stronger absorption of the WCR occurs on the right-hand-side of the image due to the greater depth of intervening wind material. The right-hand column of Fig. 7 shows the system viewed at phase 0.375, halfway between quadrature and conjunction. The image now displays a prominent double-helix-like structure formed through the limb-brightening of the WCR. The vertical curvature of the WCR induced by the increasing phase-lag of shocked gas with height/depth from the orbital plane is strikingly displayed (see also the warped shock surfaces of the WCR shown in Paper I). At low frequencies, this double-helix-like structure is partially obscured by the foreground wind (now to the left of the image).

Fig. 8 displays images obtained for an observer with an inclination angle . The brightest parts of the WCR are again those which are limb-brightened, with the overall morphology of the WCR being shaped like an “S”. In some of the images displayed (e.g. the middle-left image with  GHz and ) the arms of the WCR do not trace back to a central point, due to the additional emission from the winds. This figure also demonstrates that the observed shape of the WCR, particularly at the lower frequencies, owes much to the location of any foreground absorber (i.e. a wind), and is not necessarily a good indicator of the actual shape and position of the WCR (a point made already in Dougherty et al. 2003).

Lightcurves of the thermal emission from model cwb2 are shown in Fig. 9. Because of the larger size of the computational volume, lightcurves at 15 and 22 GHz are also displayed. The lightcurves of the total observed emission differ distinctly in their shape compared to those from model cwb1 (left panels of Fig. 3) - those from model cwb2 display maxima which are less broad. This reflects the differences in the optical depth through the WCR in the two models. Fig. 10 illustrates that the emission from the WCR of model cwb2 becomes optically thin at relatively low frequencies, in contrast to that of model cwb1 where it remains optically thick to  GHz.

The free-free emission from the unshocked winds at 43 GHz is still substantially below the combined emission from the winds of two single stars (panel h in Fig. 9), but compared to model cwb1 (see Fig. 3b) the flux from model cwb2 is roughly twice as high. This reflects both the larger computational domain in model cwb2, which means more of the emission at 43 GHz is captured, plus the wider separation of the stars i.e. the WCR removes a smaller segment of the 43 GHz emitting region from the unshocked winds. Once again, it is only by 250 GHz that the flux from the unshocked winds is comparable to twice the flux from a single star. And, as was also the case for model cwb1, by 1000 GHz the unshocked wind flux exceeds twice the single star value (due to the emission probing deep into the acceleration region of the winds). However, the 1000 GHz lightcurve of unshocked wind emission from model cwb2 (Fig. 9q) is almost completely constant (the main exception being a triangular shaped minimum centered on each conjunction), in stark contrast to the greater variability seen from model cwb1 (Fig. 3k). The minimum in Fig. 9(q) fails to approach 50 per cent of the unshocked wind flux. This is expected, since two perfectly symmetric winds around stars in conjunction will not completely occult, because wind emission from the background star will not be completely absorbed for impact parameters where the winds start to get optically thin. Therefore, a perfect occultation by two stars with stellar winds is not possible, and hence a dip of 50 per cent is never reached. An additional (smaller) part may be played by the fact that the winds are no longer spherically symmetric, since their acceleration is inhibited towards each other, but is boosted in other directions. A slight lag in the phase of the minima relative to conjunction is increasingly apparent at lower frequencies.

Now consider the observed emission arising from the WCR (see the right column of Fig. 9). A maximum in the emission again occurs near quadrature, just as in model cwb1, though this time for a totally different reason. In model cwb2, maxima occur when the line of sight of the observer to the apex of the WCR (where the intrinsic emissivity is greatest) is through one of the low-opacity arms of the WCR. This contrasts with model cwb1, where maxima occur when the optically thick WCR is oriented almost face on.

Although the free-free emissivity per unit volume in the WCR is much smaller in model cwb2 than in cwb1 due to the reduced densities and higher temperatures, the optical depth within this material is also reduced. The net result is that the observable flux from the WCR at  GHz is much higher in model cwb2 than in model cwb1. But above  GHz, when the opacity through the circumstellar environment is low, the reduced densities in the WCR limit the intrinsic emission compared to model cwb1, so that the observed flux from the WCR in model cwb1 exceeds that from model cwb2 (compare Figs. 4 and 10).

At 43 GHz the free-free flux from the WCR varies by a factor of 2 for observers with a inclination (Fig. 9i). At lower frequencies the variation is even more pronounced, as the relative contribution of free-free emission from the WCR becomes increasingly dominant (the free-free flux from the winds scales as , whereas it scales as from the optically thin WCR).

Interestingly, for  GHz, the cwb2 WCR lightcurves display a remarkable resemblance to the  keV X-ray lightcurves from this model shown in Paper III, including the asymmetrical slope of the minimum and the faster rise to maximum than the fall from it. Moreover, the 100 and 250 GHz WCR lightcurves (panels l and o in Fig. 9) are similar to the  keV X-ray lightcurve. On the other hand, the  keV X-ray lightcurve is most similar to the 43 and 100 GHz lightcurves from the unshocked wind emission (see panels h and k in Fig. 9).

Radio-to-sub-mm spectra for a variety of viewing orientations are shown in Fig. 10. The emission from the WCR produces a “bump” on top of the unshocked winds spectrum in a manner which is not too dissimilar to a weak synchrotron source. Note, however, that this work only considers the free-free emission from the WCR. Such behaviour was previously noted in Pittard et al. (2006). It is interesting that Contreras et al. (1996) purposefully observed a sample of O and WR stars at 43 GHz to minimize the effects of non-thermal emission on determinations of mass-loss rates. In this respect Fig. 10 shows that even at 43 GHz the contribution to the total flux due to thermal emission from the WCR may not be negligble, complicating this goal.

Nevertheless, at high frequencies the thermal emission from the unshocked winds always dominates the emission. This occurs at  GHz in model cwb2. Naively, therefore, one might expect that sub-mm observations may be preferred to radio observations in respect of determining the mass-loss rates of the winds. However, there is then the additional complication that the acceleration regions of the winds may be being probed. Another problem is that the effects of clumping in the winds may be more severe at mm wavelengths than in the radio (see e.g. Runacres & Blomme 1996; Blomme et al. 2002). Moreover, the flux from sub-mm observations is typically not as well calibrated as the flux from radio observations.

The lower panels of Fig. 10 show the spectral indices as a function of frequency. (as shown by the pink line) is consistently around +1.0, irrespective of the orientation of the system to the observer. At lower frequencies is slightly larger, though this may simply reflect some of the flux loss which occurs from the model due to the finite grid size. At high frequencies the emission from the WCR has a spectral index . increases at lower frequencies as the circumstellar environment absorbs some of the intrinsic free-free emission from the WCR. The spectral index of the total free-free emission, , lies between and . Between 43 and 1000 GHz, , +0.67, and +0.65 in Figs. 10(d)-(f). These are comparable to the values observed from some stars (e.g. Williams et al. 1990; Leitherer & Robert 1991; Altenhoff et al. 1994; Nugis et al. 1998). However, the spectral index from the unshocked winds over this frequency range is steeper (, +0.99, and +1.00), being offset by the flatter spectra from the WCR emission (, +0.00, and , in Figs. 10(d)-(f), respectively).

The variation of the total thermal emission for an observer in the orbital plane between conjunction and quadrature is shown in Fig. 5(b). The emission at quadrature is dominant over the emission at conjunction over the entire frequency range considered.

#### Model cwb3

The left panel of Fig. 11 shows intensity images of the free-free emission from model cwb3 (an unequal winds simulation) at 43 GHz for an observer directly above the orbital plane. A comparison against the 43 GHz synthetic image in Fig. 6 reveals several differences. Firstly, the position of maximum brightness (at the apex of the WCR) is different in the two images: in Fig. 11, the apex of the WCR is pushed closer to the O8V star (which is the southern star in this image). Secondly, the emission from the downstream arms of the WCR is also in different locations, reflecting the changed positions of ram-pressure balance between the two winds. Finally, there is a greater emission contrast between the leading and trailing edges of the leading arm of the WCR, which reflects the changes in the underlying density distribution (see Paper I). Conversely, there is also a reduced emission contrast between the leading and trailing edges of the trailing arm, compared to that from model cwb2 (again see Paper I for an explanation of the density distributions within the WCRs of these models).

The right panel in Fig. 11 shows a synthetic image at 250 GHz for an observer with and . The greater emission from the O6V wind is clearly visible, along with the “S”-shaped emission from the WCR.

Fig. 12 shows intensity images of the free-free emission from model cwb3 at 100 GHz as a function of phase for an observer located in the orbital plane (). Again, significant differences are apparent compared to the 100 GHz images from model cwb2 (see Fig. 7). At conjunction with the weaker wind from the O8V star in front (), there are two main differences. First, the limb brightened edge on the leading arm (to the right side of the image) is projected closer to the centre of the image and shows greater curvature, while the corresponding region of limb brightening on the trailing arm is located off the side of the image and is not visible. Second, the reduced radius of the optical depth unity surface in the O8V wind creates a smaller silhouette against the brighter WCR than occurs in model cwb2.

Differences also occur at other viewing angles. At phase 0.125 (), the WCR is seen mostly face on, with the foreground O8V wind silhouetted against it (slightly to the right of centre). The vertical curvature of the WCR is again clear at phase 0.25 (quadrature), reflecting the wrapping of the WCR around the O8V star (to the right of centre) due to the stronger wind from the O6V star (to the left of centre). The apex of the WCR, which occurs closer to the O8V star (reflecting the position of the ram pressure balance of the winds), is seen sideways on. The greater emission from the O6V wind, relative to that from the O8V wind, is also clearly visible. The double-helix-like structure seen from model cwb2 at phase 0.375 disappears in model cwb3. Instead, the brightest regions of the projected emission from the WCR primarily traces the limb-brightened, dense, O8V gas on the trailing edge of the leading arm (see Paper I).

The images from the second half of the orbit ( phase ) resemble those in the upper row of Fig. 12 which show the first half of the orbit. However, some small differences are seen, which lead also to differences in the lightcurves. For example, the conjunction at phase 0.5 now has the denser O6V wind in the foreground, whereas the weaker O8V wind was in front of the WCR at phase 0.0. This changes the nature of the silhouette.

Fig. 13 shows lightcurves at a range of frequencies from radio to sub-mm of the free-free emission from model cwb3. The lightcurves are reasonably similar to those shown in Fig. 9 from model cwb2, so we show a smaller sample of frequencies. The slight differences in the total flux and the morphology of the variation have their origin in the slighter weaker secondary wind in model cwb3. The maxima now show asymmetrical peaks due to the different wind strengths. The highest peak occurs when viewing down the trailing arm. This is straighter than the leading arm, as the curvature of the WCR opposes the coriolis-induced curvature, and allows more flux to escape through the low opacity plasma within the WCR.

Fig. 14 shows radio-to-sub-mm spectra from model cwb3 for three different viewing orientations. Strong similarities to the results from model cwb2 are again seen. Since model cwb3 consists of O6V and O8V stars, rather than two identical O6V stars as in the other models, we also plot in Fig. 14 the individual contributions of the unshocked winds to the total thermal flux. As expected, the flux from the denser O6V wind usually dominates that from the O8V wind (see Fig. 14b and c). However, for favourable orientations (e.g. when the O8V star is directly in front of the O6V star), its contribution to the total thermal flux can be comparable to, and at some frequencies exceed, the contribution from the O6V wind (Fig. 14a).

The spectral index as a function of frequency is shown in Fig. 14(d)-(f). The values of are very similar to those in Fig. 10(d)-(f) from model cwb2. Note that the high frequency spectral index from the O8V wind exceeds that from the O6V wind. This is because the emission probes deeper into the acceleration region of the O8V wind at a given (high) frequency.

Fig. 5(c) shows the variation of the total thermal emission for an observer in the orbital plane between conjunction (phase 0.0, O8V star in front) and quadrature (phase 0.75). Again, the behaviour is similar to that from model cwb2, though the variation between phases is slightly weaker. The flux, of course, is also reduced.

#### Model cwb4

Unlike the previous models which all had circular orbits, model cwb4 simulates a CWB with an eccentric orbit (). The intrinsic emission now varies with phase, whereas it was constant in the circular orbit models cwb1cwb3. Fig. 15 shows intensity images of the free-free emission from model cwb4 at 1000 GHz for an observer directly above the orbital plane. The images show striking variations in their brightness and morphology as a function of orbital phase, reflecting the dramatic changes in the WCR during the orbit (see Paper I for full details of the hydrodynamics).

Discussion of this figure begins, most easily, at orbital phase 0.5 (apastron). At this phase we see the emission regions of the two winds (located very close to the stellar surfaces due to the high frequency), plus between them emission from the hot WCR. There are also small but very bright regions of emission within and close to the WCR. These are very dense, cold, clumps of previously shocked gas which formed during the previous periastron passage when the WCR became highly radiative. Most of the dense, cold, gas which formed during this time has now been cleared out of the inner parts of the system. The dense clump and its tail visible in the top left of the image are actually located outside of the WCR, since the high inertia of the clumps causes them to move on almost ballistic trajectories, while the hot gas in the WCR responds more rapidly to the changing positions of the stars (see Paper I for more details).

The process of clearing out the remaining dense clumps continues as the stars advance in orbital phase. By phase 0.7, no clumps are left within the central regions of the system. At this point the stars are moving in a slow waltz towards each other. The density within the WCR steadily increases as the stellar separation reduces, resulting in brighter emission from the WCR, so that by phase 0.95 the apex of the WCR is brighter than the unshocked winds themselves, despite the high frequency! The spatial distribution of the emission from the WCR also reflects the WCR’s increasing curvature as periastron is approached.

By periastron the gas in the densest parts of the WCR in model cwb4 has undergone substantial radiative cooling. Cool ( K), dense gas forms in the WCR between the stars and for some distance downstream (see Figs. 11 and 15 of Paper I). This gas is optically thick to frequencies exceeding 1000 GHz. The enhancement in the density of this gas relative to that of the neighbouring unshocked winds means that its emission coefficient is significantly higher, as indicated by the images shown in Fig. 15. Further downstream, the plasma in the WCR remains hotter, with some gas at temperatures exceeding  K (see Fig. 15(e) in Paper I). This plasma is optically thin to emission above  GHz, and its emission coefficient is much smaller than that of the optically thick cold gas nearer the apex of the WCR. The wide variation in plasma density and temperature throughout the WCR demonstrates that models of the emission from such systems must be based on time-dependent hydrodynamical calculations.

Even more mass within the WCR has cooled back to the wind temperature ( K) by phase 0.05. Newly shocked gas in the WCR becomes largely adiabatic by phase 0.2, as the stars continue to separate and the pre-shock velocity of the winds increases. No new clumps form after this phase until the next periastron passage. Instead, the slow process of clearing the existing clumps out of the system begins.

It is likely that the majority of the emission between orbital phase is captured by the numerical grid, but it is clear that significant flux losses occur from the simulation between phase as optically thick clumps leave the finite boundaries of the grid. This effect is likely to be greatest at the lower frequencies since the clumps stay optically thicker for greater outflow distances, and should be kept in mind when examining the remaining figures in this section.

Fig. 16 displays lightcurves of the free-free emission from model cwb4. Since the typical temperature of the plasma in the WCR varies between “cold” near periastron, and “hot” for the rest of the orbit (excepting the dense clumps formed within the WCR), we have chosen to separately identify the emission from cold ( K) and hot ( K) gas, rather than from the unshocked winds and the WCR. At times when the WCR is hot (such as at apastron), this leads to the same identification of gas as in the previous sections. However, during the small range in phases around periastron, our analysis groups the cold gas within the WCR with the unshocked gas from the winds, rather than with the hotter parts of the WCR on the leading arms and downstream. Hence, the “cold” and “hot” components cannot be directly compared to the “winds” and “WCR” components from model cwb1, although as noted above it is still reasonable to compare to model cwb2 at other phases. Note that Fig. 16 shows observed fluxes, so, for example, the contribution of hot material to the flux includes absorption from overlying cold material in the circumstellar environment.

Fig. 16 shows that the emission typically reaches a maximum around periastron, as was already obvious from the images at  GHz and , presented in Fig. 15. The varying circumstellar absorption due to the passage of the stars in front of the WCR complicates the interpretation of the lightcurves when the observer is in the orbital plane. Therefore, we shall concentrate on the lightcurve for an observer at . At 43 GHz the emission peaks just prior to periastron, while at 1000 GHz the maximum occurs just after periastron. There are two reasons for this behaviour. One reason is related to the changes with orbital phase of the optical depth of sightlines through the WCR, and is discussed in more detail below. A minimum centered on phase 0.1 also occurs in the 43 GHz lightcurve - this is primarily due to the high optical depth of the cold dense plasma formed in the WCR during the periastron passage. Most of the emission from this material is therefore absorbed. By phase 0.3 many parts of the WCR become optically thin, and the emission from background material starts to become visible. However, the emissivity of this material is also dropping, and this becomes the dominant factor affecting the observed flux at subsequent phases. The 43 GHz emission reaches a minimum around phase , and then subsequently increases as the stars approach each other, due to the increase in density of the hot plasma in the WCR.

In contrast, the 1000 GHz lightcurve of the total emission (Fig. 16j) does not display a minimum at phase 0.1. This is partly because some parts of the WCR remain optically thin at this phase (although a small reduction in the emission from “hot” plasma relative to the trend is nevertheless visible in Fig. 16l). It is also partly because some of the emission at 1000 GHz is from the unshocked stellar winds (which is responsible for some of the emission seen in Fig. 16k), and the stars are near quadrature at phase 0.1.

Fig. 17 shows radio-to-sub-mm spectra from model cwb4 as a function of orbital phase for an observer in the orbital plane at . The spectrum from the WCR is optically thin between phase 0.5 and 0.9 over the displayed frequency range. Comparing against the spectrum from model cwb2 in Fig. 10(a), we see that both the “hot/WCR” and the “cold/winds” flux are in good agreement at high frequencies. However, there is less emission from model cwb4 than from model cwb2 at low frequencies. At 10 GHz, the “winds” flux from model cwb2 is about three times greater than the “cold” flux from model cwb4, while the “WCR” flux is more than an order of magnitude greater than the “hot” flux. This lack of emission from model cwb4 at the lower frequencies reflects the smaller grid compared to that of model cwb2.

Since a large fraction of the emission from the WCR is optically thick at periastron, it is not surprising to find that the spectral index of the “hot” emission at phase 1.0 becomes positive like the unshocked winds ( between 15 and 100 GHz - see Fig. 18b). In fact, the “hot” spectrum is remarkably straight at periastron, with very little curvature (see Figs. 17a and k, and Fig. 18b).

Due to the dramatic cooling at periastron, by phase 0.1 the majority of the gas within the WCR is now cold. This results in a much lower flux from the “hot” component of the emission, which is now mostly low density gas on the leading edges of the WCR and a narrow strip of rapidly cooling gas near the apex and trailing edges of the WCR. Conversely, the emission from the “cold” component is significantly bolstered by the substantial amount of gas which has cooled within the WCR, being more luminous at 100 GHz than the emission at periastron.

The change in viewing angle into the system also causes the low frequency emission to drop (see e.g. panel l of Fig. 17). A gradual recovery occurs in the low frequency emission up to phase 0.3, after which it declines again. These changes are largely driven by changes to the optical depth to the “hot” flux, which shows more variability at low frequencies than the “cold” flux. The level of variability is expected to be smaller in calculations where the hydrodynamical grid extends to larger radii.

Although some parts of the WCR remain hot even at periastron (for instance the leading edge of the shocks), the majority of the mass within the WCR is optically thick at phases 0.0 and 0.1. At the highest frequencies considered, the WCR is only marginally optically thick at phase 0.1, with the optical depth through its surface occasionally exceeding unity, but more typically being less than this. The WCR becomes increasingly optically thin thereafter, as the stars continue their separation.

Just as the flux from the hot part of the WCR varies with phase, so does the contribution from the unshocked winds plus cold parts of the WCR. At phase 0.5 this part of the emission has a straight slope with between 200 and 1000 GHz. The flux and spectral slope are relatively steady up to phase 0.9. At periastron, the first gas in the WCR to cool back down to K creates a noticeable “bump” on the cold gas spectrum at  GHz (see Fig. 18b). The extra emission becomes more extended in its range of frequencies as more and more gas passing through the WCR cools back to the temperature of the unshocked winds, and reaches its maximum contribution at phase 0.1, with as much flux from cold gas compared to phase 0.5 at 200 GHz. Thereafter the bump declines as the cold plasma in the WCR is ablated and/or expelled out of the system, so that the “cold” spectrum is almost entirely generated by emission from the unshocked winds by apastron.

These changes to the spectra also affect the spectral index of the emission, as shown in Fig. 18, with the most rapid changes to occuring around periastron passage. At phase 0.9, the spectral index of the total emission increases from at 55 GHz, to at 1 THz, due to the steady increase of flux from the cold unshocked winds relative to the hot WCR. At periastron, the increasing opacity of the hot WCR leads to a significant increase in to between 20 and 400 GHz. The enhancement to the emission around 100 GHz due to the formation of cold gas increases to nearly at 100 GHz. The spectral index of the total emission is in the range to between 20 and 600 GHz. By phase 0.1, monotonically decreases from low to high frequencies. However, the computed values of the spectral index at low frequencies shown in Fig. 18(c) are unlikely to be quantitatively correct, given the known underestimate of the flux due to the finite grid size.

Fig. 19 shows the variation of the total thermal flux with orbital separation for an observer with . Bearing in mind that our simulations likely underestimate the true flux from this system between phase , it nevertheless reveals that there is a strong hystersis to the emission at all frequencies. It is also interesting that at low frequencies the emission is stronger as the stars approach periastron, whereas at high frequencies the emission is stronger as the stars recede from periastron. This is because at low frequencies the winds are very opaque and the emission most easily escapes when the WCR is full of hot plasma, which is the case in the second half of the orbit. In contrast, at high frequencies the winds are always transparent, and stronger intrinsic emission from the cold, dense, plasma which exists in the WCR following periastron passage can be observed.

Fig. 19 also shows the theoretical flux from a single O6V wind which accelerates instantaneously to its terminal velocity (lower horizontal line). The upper horizontal line shows twice this value. We have already seen from many of the lightcurves in this paper that at the lower frequencies the models at times significantly underestimate the true flux due to the finite extent of the numerical grid, and we again see this here. At higher frequencies (e.g.  GHz in model cwb4) this is no longer an issue, at least over orbital phases .

The flux at 43 GHz from model cwb4 typically exceeds twice the theoretical single star flux, although it is below this value at phases 0.05 and 0.1 when the stars orbit deep into each other’s wind - at these times the optically thick contribution from the WCR is not able to offset the loss of flux from the unshocked winds due to the presence of the WCR (however, particularly at phase 0.1, we expect the true flux to be significantly underestimated due to optically thick, cool, clumps flowing through the grid boundaries). By 250 GHz, the total flux from model cwb4 always exceeds twice the expected single star value. This is again because the characteristic region of emission is now close enough to the stars that the WCR no longer “cuts” out significant parts of this region. At phase 0.5 the extra emission from the WCR, which is mostly from hot gas, is negligible (see Fig. 17f). However, as the stars approach periastron the emission from the WCR becomes significant as the plasma within it increases in density and decreases in temperature, causing the flux rise shown in Fig. 19(e). The dip in the flux at periastron is caused by the foreground wind occulting the apex of the WCR where the brightest emission arises.

Pittard et al. (2006) noted that the thermal radio flux from the WCR in an adiabatic system scales as , so that we would expect a variation of just over a factor of two from our model. However, the variations from model cwb4 are obviously much greater. For example, Fig. 16(l) shows that for an observer with , the flux increases by a factor of 3.4 between phases 0.6 and 0.9, compared to a increase of 1.7. This is because the pre-shock velocities change as the stellar separation varies, affecting the post-shock temperature, whereas the analysis in Pittard et al. (2006) is only valid in the scale-free limit where the winds collide at constant velocity.

Finally, we note that the brightness of the WCR emission is highly dependent on its density. In models cwb1 and cwb4 the cooled postshock gas becomes very dense, partly because of the assumption that there is no significant magnetic pressure which would otherwise support this gas against collapse. If magnetic pressure were to become significant, the densities in the WCR would be limited and the emissivity of the WCR would be reduced.

#### Clumpy winds

Figs. 20 and 21 show how the lightcurves and spectra from model cwb2 are modified if it is assumed that the unshocked winds (and at times also the shocked gas in the WCR) are clumpy (note that the mass-loss rates remain the same). Fig. 20(b) shows how clumping with a volume filling factor increases the flux from the unshocked winds compared to the smooth winds case. The flux increases by a factor of 2.8, which is a little below the expected increase of 4.64 from the relation , and is likely due to the greater loss of flux from the numerical grid because of its finite size.

Fig. 20(c) shows how the flux from the WCR varies with clumping. Comparing the smooth to the results, we see that the larger optical depth unity surfaces in the winds of the latter calculation lead to greater attenuation of the WCR flux. Hence the observable flux from the WCR can be affected by changes to the clumping in the unshocked winds, even though in both of these calculations the WCR is assumed to be smooth and the intrinsic emission from the WCR is identical.

Pittard (2007) showed that an adiabatic WCR rapidly smoothes out clumps if their density contrast and size is not too high, so that the destruction timescale of the clump is shorter than the dynamical timescale for material within the WCR to flow out of the system. However, this work also revealed that the plasma in the WCR contains a multitude of weak shocks, and is far from being perfectly smooth. To account for this, we have also computed calculations where there is some “residual clumping” within the WCR. Fig. 20(c) shows that the flux from the WCR gradually increases as the WCR is made more clumpy. When , the observed flux from the WCR at quadrature is roughly equal to that from the smooth winds case, with the extra absorption from the unshocked winds (with ) being offset by extra intrinsic emission from the WCR.

At the extreme limit of there being no “smoothing” effect by the WCR (i.e. ), the observed flux from a clumpy WCR can dramatically exceed its smooth counterpart, as shown in Fig. 20(c) when . However, we regard this possibility as extremely unlikely in systems with an adiabatic WCR, and more typically would expect when . Indeed, the level of the continuum in the theoretical X-ray spectra shown in Fig. 3 of Pittard (2007) shows that the intrinsic flux increases by only 20 per cent7, equivalent to the WCR smoothing the clumping from to .

Fig. 21 shows the spectra obtained from calculations where the winds are significantly clumped () and the degree of clumping within the WCR is varied between (i.e. a smooth WCR) and (i.e. the WCR has no smoothing effect on the clumps). In Fig. 21(a) we see that clumping increases the free-free flux from the winds, at the expense of increased absorption of the free-free flux from the WCR (particularly at the lower frequencies). Increasing the clumping within the WCR has little noticeable effect at first (see Fig. 21b), but large enhancements in the flux are apparent when the clumping within the WCR is extreme (see Fig. 21c). In fact, in such cases the increase in flux from the WCR exceeds the increase seen from the winds, due to the differences in optical depth within these regions.

Figs. 20 and 21 demonstrate that even with considerable clumping in the winds, free-free emission from the WCR can still escape and create variability with orbital phase.

## 5 Comparison with observations

### 5.1 O-star binaries

The parameters in model cwb1 are similar to those in HD 215835 (DH Cep), HD 165052, and HD 159176. Currently, none of these systems is detected in the radio or sub-mm bands, though Bieging et al. (1989) note 6 cm upper limits of 0.27 mJy for HD 165052, and 1.05 mJy for HD 159176. Since we do not capture all of the flux from model cwb1 (particularly at the lower frequencies investigated), a direct comparison is not possible. However, we can obtain a very conservative upper limit to the 5 GHz flux from model cwb1 by extrapolating the 1000 GHz flux to 5 GHz with a spectral index of (the true spectral index may be significantly steeper than this). This yields a 5 GHz flux of  mJy, which is still consistent with the observational upper limits noted above, even if moderate clumping was added to our models.

However, these systems, and others like them, should be detected by the next generation of interferometers. For instance, the EVLA will have a point source sensitivity better than 1 Jy between  GHz (1 rms in  hrs). e-MERLIN’s specifications are almost as good. The Square Kilometre Array (SKA), which will operate over  GHz, will probe even deeper. ALMA’s band 3 ( GHz) also has a continuum sensitivity of 0.06 mJy, and should be just sensitive enough to detect these systems. Obviously the chances for detection improve if the winds are clumped.

Model cwb2 is similar to HD 93161A, and also shares some similarities with Plaskett’s star (HD 47129). Benaglia & Koribalski (2004) failed to detect HD 93161A with ATCA at 3 and 6 cm. Observations centered on HD 93129A (about 3 arcmin from HD 93161A) reveal that at the position of HD 93161A the rms noise is  mJy/beam at 3 cm, and  mJy/beam at 6 cm, giving upper limits of 0.6 mJy and 0.75 mJy, respectively. These limits are consistent with the fluxes predicted from model cwb2. On the other hand, Plaskett’s star was detected at 5 GHz by Persi et al. (1988), who reported a flux of 0.2 mJy. This flux is higher than we see from model cwb2, and may indicate significant clumping within the winds of this system. The fact that Plaskett’s star is at a distance of 1.5 kpc (rather than the 1 kpc of the models in this work) makes this difference even worse.

Model cwb4 bears some similarities to some well-known O+O systems with eccentric orbits, including HD 152248, HD 93205, HD 93403, Cyg OB2 #8A, and  Ori. Of these systems, HD 152248, HD 93205, and HD 93403 are currently not detected in the radio. HD 152248 fell within the field of view of ATCA observations of Sco OB1 (Setia Gunawan et al. 2003), and upper limits (at both 3 cm and 6 cm) are 0.24 mJy (Stevens, priv. communication, 2009). In contrast, Cyg OB2 #8A is a bright, non-thermal, and strongly variable radio source (Bieging et al. 1989). Following the discovery of this system’s binarity (De Becker, Rauw & Manfroid 2004), Blomme (2005) showed that the radio flux density undergoes phase-locked variations8. The non-thermal emission from our models will be investigated in a future paper.

Finally,  Orionis (HD 37043) was detected by Howarth & Brown (1991) at 3.6 cm, with a flux of  mJy (see Lamers & Leitherer 1993). Since it was observed only at one frequency it is not known if this emission is predominantly thermal or non-thermal. Moving the system to a distance of 1 kpc yields a flux of  mJy, comparable to the predicted 8 GHz flux from model cwb2, and similar to the flux from model cwb4 at certain phases.

### 5.2 WR and LBV systems

Although this work has focussed on O+O systems, thermal radio emission is also seen from WR and LBV stars (e.g. Abbott et al. 1986; Leitherer et al. 1997; Dougherty & Clark 2008). However, there are few accounts in the literature on the variability of such sources9. For instance, despite a large multi-frequency monitoring campaign of  Velorum, little evidence for variation was found (Sean Dougherty, priv. communication, 2009).

Radio variability at the per cent level, has, however, been recently reported from four WN9 stars in the Arches cluster (sources AR1, AR3, AR4, and AR8: Lang et al. 2005). All but AR1 have a positive spectral index (; observational uncertainities do not allow an accurate measurement of the spectral index of AR1). Lang et al. (2005) suggested the variability might indicate the presence of a time-variable non-thermal component, and/or changes to the underlying mass-loss rate or velocity of the wind. However, as this work demonstrates, an alternative explanation is that these sources are binaries, and that the variability is caused by changes to the circumstellar opacity along sightlines into the system as the stars orbit each other, or due to a non-zero orbital eccentricity. This view is supported by the fact that AR1 and AR4 (identified as sources A1N and A1S by Wang, Dong & Lang 2006) have strong and hard X-ray emission with a prominent Fe K line at 6.7 keV, all of which are characteristics of colliding wind binaries10.

Three point-like radio sources were also detected in the Quintuplet cluster (QR6, QR8, and QR9: Lang et al. 2005). While these all have positive spectral indices indicative of a strong stellar wind, we have demonstrated in this work that the bulk of the emission could actually arise from a radiative, and optically thick, WCR. In Westerlund 1, six WR stars are detected in the radio, but only one has a thermal spectrum (WR L; Simon Clark, priv. communication, 2009). Note, however, that WR L looks like it might be a CWB from it’s X-ray properties (Clark et al. 2008).

Variable thermal emission with has also been detected from WR 89, WR 113, and WR 138 (Montes et al. 2009). For WR 89, the spectral index of the emission, , while for WR 113, , and for WR 138. WR 89 is classified as a visual binary (van der Hucht 2001), while WR 113 is a WC8+O8-9 system with an orbital period of 29.7 d (Niemela et al. 1999). WR 138 is a spectroscopic binary with a period of 1538 d (Annuk 1990), though this system may also harbour a close companion with a short period of order days (Lamontagne et al. 1982; Moffat & Shara 1986). It is likely that the observed variability is related to the binarity of these systems. Montes et al. (2009) also note that the radio emission from WR 98a and WR 104 (two of the so-called “pinwheel” systems) is variable and the spectral index is positive, but less than the value expected from a single wind (i.e. ). These systems could harbour non-thermal emission (i.e. be “composite” systems), or the emission could be purely thermal, combining optically thick emission from the winds with optically thin emission from the WCR.

It is also interesting to note that the broad maxima seen in the lightcurves of model cwb1 are similar to those seen in the radio lightcurve of the LBV  Carinae (Duncan & White 2003). The latter is usually explained in terms of time-variable ionization of dense gas within the spoked equatorial disc seen in the optical (see also Kashi & Soker 2007). The spatial extent of the brightest part of the radio emission spans a width of 1 to several arcseconds, depending on the orbital phase. At a distance of 2.3 kpc, this corresponds to thousands of AU. However, 3-dimensional simulations of the structure of the WCR reveals density enhancements on similar spatial scales (Parkin et al. 2009; Gull et al. 2009). Thus, it is possible that at least some of the variability could be due to time variations in the properties and orientation of the optically thick shocked wind of the primary LBV star. Furthermore,  Car showed a sharp peak in the 7 mm flux density in 2003.5 (and again in 2009.0, Abraham, private communication, 2009), coinciding with periastron passage (Abraham et al. 2005a). This spike has been interpreted as the shocked plasma from the secondary wind becoming optically thick at high temperatures (Abraham et al. 2005b). It is interesting that we also see a sharp maximum in the flux from model cwb4 near periastron. However, for this we have a different interpretation. At low frequencies the spike in model cwb4 occurs before periastron, when the density and thus emissivity of the WCR is increasing. However, the WCR is optically thin at this time. In contrast, at high frequencies ( GHz), the greatest flux from model cwb4 occurs after periastron. At such frequencies, the majority of the unshocked winds is optically thin, allowing the bright emission from the cold dense clumps created during periastron passage to be seen. It is clear that further modelling of  Car is necessary, in order to better understand the observed flux variations, on both short and long timescales, from this fascinating system.

## 6 Summary and conclusions

This work examines the thermal radio-to-sub-mm emission from colliding wind binaries, focussing in this first instance on short-period O+O binaries (the potential non-thermal emission will be studied in due course). We investigate how the thermal emission is influenced by the presence of a nearby companion star, and the extent to which variability with orbital phase occurs. The emission is calculated from grids of density and temperature generated from the 3D hydrodynamical models presented in Pittard (2009a). The models span a variety of different physical conditions, including radiative and adiabatic wind-wind collision regions in systems with circular orbits, and also a system with an eccentric orbit where the WCR repeatedly cycles from radiative to adiabatic and back to radiative again.

We find distinct differences in the spectra and lightcurves from systems with radiative and adiabatic WCR’s. In the former (model cwb1), the WCR is optically thick to wavelengths as short as sub-mm, and radiates as a dense sheet of emission. The WCR dominates the total emission from the system, and only begins to show signs of becoming optically thin at  GHz. The flux (which has a positive spectral index) can be over an order of magnitude greater than from the wind of a single star, and highlights potentially serious consequences for observational estimates of mass-loss rates where sources are assumed to be single stars, but are in fact binaries. The flux is highest just before quadrature, when the inner regions of the WCR are viewed almost face on (due to the large aberration of the WCR), and have the largest projected extent. The lightcurve maxima are broad, with minima occuring slightly before conjunction. For an observer in the orbital plane, the total flux (from the unshocked winds plus the WCR) varies by a factor of two or so. Although there is some loss of flux from the models, the spectral index of the emission is likely to be steeper than +0.6, and varies with orbital phase and the orientation of the observer. Synthetic images for an observer directly above the orbital plane display an “S”-shaped region of emission which traces the location of the WCR.

In contrast, the emission from adiabatic WCR’s is optically thin, and most easily escapes the system when sight-lines to the observer are themselves through the WCR. In model cwb2 the WCR dominates the emission at  GHz, with free-free absorption creating a low frequency turnover typically at  GHz. The emission from the WCR typically contributes only 10 per cent of the total flux at 1000 GHz (0.3 mm). The lightcurves show relatively narrow maxima at phases just after quadrature. The broad, and relatively flat minima inbetween the maxima indicate phases where there is considerable absorption of the emission from the WCR. There is a gradual drop in the degree of variability with frequency as the circumstellar environment becomes less opaque and increasingly transparent. At 1000 GHz the vast majority of the WCR emission escapes the system. Between 100 and 250 GHz, the lightcurves change morphology. At  GHz, the maxima become almost completely flat-topped, and the minima are now caused by occultation of a narrow emission region near the surface of each star. The spectral index is steepest at low frequencies (where absorption of emission from the WCR is important) and at high frequencies (where the emission is almost entirely from deep within the acceleration region of the individual winds). Images of the emission for an observer in the orbital plane reveal an intertwined “double-helix” showing limb-brightened parts of the WCR, against which the foreground wind is sometimes silhouetted. At other viewing angles (e.g. ) an “S”-shaped region of emission is again seen.

Adiabatic systems where one of the winds is weaker than the other generally display similar properties, the main differences being slight asymmetries to the lightcurves (e.g. differing heights to the maxima), and a change from a “double-helix” emission structure to a single limb-brightened structure emitted by the relatively dense trailing edge of the trailing arm of the WCR.

One of the more surprising findings in this work is the strength of the flux variations with orbital phase from a system with an eccentric orbit (model cwb4), which exceed an order of magnitude from the WCR), and may be almost as strong from the system as a whole. This is caused by dramatic changes to the physical properties of the WCR, which is highly radiative at periastron, but adiabatic at apastron, and changes from optically thick to optically thin, respectively. This creates a strong frequency-dependent hystersis to the emission around the orbit. At low frequencies the emission is stronger as the stars approach periastron, due to the low opacity through the hot WCR at this time. The flux variations easily exceed the scaling proposed by Pittard et al. (2006), which is invalidated by changes to the pre-shock wind speeds as the stars orbit deep within the acceleration zone of each other’s wind. In contrast, at high frequencies the emission is stronger as the stars recede from periastron, because of the high intrinsic emission from dense, cold, post-shock gas created during periastron passage. Even larger variations are seen when the observer is in the orbital plane, indicating that changes in the circumstellar absorption also play an important role.

An examination of the effects of clumping reveal that in addition to higher fluxes, clumping can also effect the variability of radio lightcurves. In the case of strong clumping in the unshocked winds but a smooth WCR, the variability decreases. However, if the clumps are not smoothed by the WCR (a case which we consider rather unlikely), the variations of the flux with orbital phase can be enhanced compared to the smooth winds case.

We find that the fluxes from our models are consistent with the mostly upper limits which exist to date from the O+O systems which they most resemble. On the other hand, there are many instances where radio variability is observed from WR sources with positive spectral indices, some with . Many of these sources are likely to be binaries, with the variation linked to changes in the absorption through the circumstellar material, as the orientation of the system to the observer changes, or perhaps due to changes in the separation of the stars. However, we also caution that systems with positive spectral indices are not necessarily thermal - instead, one could be witnessing the absorption of a non-thermal component. This possibility must be carefully considered, especially if the spectral index is calculated at relatively low frequencies (e.g. between 5 and 8 GHz). To safely distinguish between thermal and non-thermal emission in close, spatially unresolved systems where the spectral index is positive may require flux measurements over a wide frequency range, additional observations at X-ray and -ray energies, and detailed comparisons with theoretical models.

It is difficult to predict how the fluxes from the models investigated in this work will scale to other systems (e.g. WR+O systems), since there are so many key variables, such as whether the WCR is optically thick or thin, the recent thermal history of the plasma in the WCR, the relative sizes of the radio photospheres in each wind to the stellar separation, the orientation of the system, etc. In many ways this work, in revealing some of these complications, demonstrates that it may not always be valid to apply some of the scaling laws noted in Dougherty et al. (2003) and Pittard et al. (2006). In this respect there is a clear need for further numerical simulations of CWBs with a range of key parameters (mass-loss rates, stellar separations, etc.). Such simulations will aid the study and interpretation of the resulting emission from these fascinating systems.

The future commissioning of more sensitive telescopes, such as the EVLA, e-MERLIN, ALMA, and SKA, should result in detections of short-period O+O systems, thus allowing the predictions of this work to be tested. Although the synthetic images shown in this work are of an angular scale which is too small to be resolved with EVLA or ALMA, planned upgrades to the VLBA offer the exciting prospect of spatially resolving the emission from tight CWBs. The EVLA, ALMA, and other facilities may also resolve some of the wider O+O (and WR+O) CWBs.

## acknowledgements

I would like to thank the referee for a helpful and timely report, Sean Dougherty and Ross Parkin for some of the inspiration behind this research, Ian Stevens for use of his modified CAK code, and Ian, Paula Benaglia, Simon Clark, and Perry Williams for providing some information on the current radio observations of early-type stars. I would also like to thank the Royal Society for a University Research Fellowship.

### Footnotes

1. pagerange: 3D models of radiatively driven colliding winds in massive O+O star binaries - II. Thermal radio to sub-mm emissionReferences
2. pubyear: 2009
3. We repeat that in this work we concentrate solely on the thermal emission from the winds of early-type stars, though observations reveal that there is also a significant population of non-thermal emitters which appears to be linked to the presence of a companion star (Dougherty & Williams 2000; De Becker 2007) and the resulting collision of the winds (see Pittard et al. 2006, and references therein, for models of the non-thermal emission). The non-thermal emission from the hydrodynamical models presented in Pittard (2009a), which are the basis for the current work, will be investigated in future papers.
4. Since this emission has a negative spectral index it can mimic a synchrotron emission component in the system - see Pittard et al. (2006) for further details.
5. In reality, there will not be such a sharp division in temperature between the wind and photosphere, and one would expect that the maximum value of the spectral index determined from observations would be closer to +2.0. Spectral indices greater than +2.0 have in fact been previously reported from Be stars (Dougherty, Taylor & Waters 1991), but we note that these were based on 3 upper limits to the flux, and so are meaningless.
6. One can estimate a rough upper limit for the amount of flux lost by taking the flux at 1000 GHz, and extrapolating to lower frequencies with a spectral index of - this reveals that the actual flux from the WCR at 43 GHz may be up to an order of magnitude higher.
7. It was this very behaviour which led Pittard (2007) to note that “clumping-independent” measurements of the mass-loss could potentially be obtained from the X-ray flux.
8. Other O+O systems with strongly varying non-thermal radio emission are HD 168112 (Blomme et al. 2005), HD 167971 (Blomme et al. 2007), and Cyg OB2#9 (Van Loo et al. 2008).
9. Of course, the variability of many non-thermal sources is much better studied.
10. On the other hand, source AR6 clearly shows a non-thermal radio spectrum (Lang, Goss & Rodríguez 2001, 2005), and since it is associated with strong X-ray emission (Wang et al. 2006), it is also likely to be a colliding wind binary.

### References

1. Abbott D. C., Bieging J. H., Churchwell E., Cassinelli J. P., 1980, ApJ, 238, 196
2. Abbott D. C., Bieging J. H., Churchwell E., 1981, ApJ, 250, 645
3. Abbott D. C., Bieging J. H., Churchwell E., Torres A. V., 1986, ApJ, 303, 239
4. Abraham Z., et al., 2005a, A&A, 437, 977
5. Abraham Z., Falceta-Gonçalves D., Dominici T., Caproni A., Jatenco-Pereira V., 2005b, MNRAS, 364, 922
6. Altenhoff W. J., Thum C., Wendker H. J., 1994, A&A, 281, 161
7. Annuk K., 1990, Acta Astronomica, 40, 267
8. Arias J. I., Morrell N. I., Barbá R. H., Bosch G. L., Grosso M., Corcoran M., 2002, A&A, 333, 202
9. Bagnuolo W. G. Jr., Riddle R. L., Gies D. R., Barry D. J., 2001, ApJ, 554, 362
10. Benaglia P., 2010, ASP Conf. Ser., “High Energy Phenomena in Massive Stars” (arXiv:0904.0533)
11. Benaglia P., Cappa C. E., Koribalski B. S., 2001, A&A, 372, 952
12. Benaglia P., Koribalski B. S., 2004, A&A, 416, 171
13. Bertout C., Leitherer C., Stahl O., Wolf B., 1985, A&A, 144, 87
14. Bieging J. H., Abbott D. C., Churchwell E., 1989, ApJ, 340, 518
15. Blomme R., 2005, in Rauw G., Nazé Y., Blomme R., Gosset R., eds, Proc. JENAM 2005, Massive Stars and High-Energy Emission in OB Associations. IAGL, Liège, p. 45
16. Blomme R., De Becker M., Runacres M. C., Van Loo S., Setia Gunawan D. Y. A., 2007, A&A, 464, 701
17. Blomme R., Prinja R. K., Runacres M. C., Colley S., 2002, A&A, 382, 921
18. Blomme R., Van Loo S., De Becker M., Rauw G., Runacres M. C., Setia Gunawan D. Y. A., Chapman J. M., 2005, A&A, 436, 1033
19. Castor J. I, Abbott D. C., & Klein R. I., 1975, ApJ, 195, 157 (CAK)
20. Clark J. S., Muno M. P., Negueruela I., Dougherty S. M., Crowther P. A., Goodwin S. P., de Grijs R., 2008, A&A, 477, 147
21. Contreras M. E., Rodríguez L. F., Gómez Y., Velázquez A., 1996, ApJ, 469, 329
22. De Becker M., 2007, A&ARv, 14, 171
23. De Becker M., Rauw G., Manfroid J., 2004, A&A, 424, L39
24. De Becker M., Rauw G., Pittard J. M., Antokhin I. I., Stevens I. R., Gosset E., Owocki S. P., 2004b, A&A, 416, 221
25. De Becker M., et al., 2006, MNRAS, 371, 1280
26. Dougherty S. M., Clark J. S., 2008, in P. Benaglia, G. Bosch, C. Cappa, eds, “Massive Stars: Fundamental Parameters and Circumstellar Interactions”, RevMexAA, 33, 68
27. Dougherty S. M., Taylor A. R., Waters L. B. F. M., 1991, A&A, 248, 175
28. Dougherty S. M., Pittard J. M., Kasian L., Coker R. F., Williams P. M., Lloyd H. M., 2003, A&A, 409, 217
29. Dougherty S. M., Williams P. M., 2000, MNRAS, 319, 1005
30. Duncan R. A., White S. M., 2003, MNRAS, 338, 425
31. Gull T. R., et al., 2009, MNRAS, 396, 1308
32. Hummer D. G., 1988, ApJ, 327, 477
33. Kashi A., Soker N., 2007, MNRAS, 378, 1609
34. Lamers H. J. G. L. M., Leitherer C., 1993, ApJ, 412, 771
35. Lamontagne R., Koenigsberger G., Seggewiss W., Moffat A. F. J., 1982, ApJ, 253, 230
36. Lang C. C., Goss W. M., Rodríguez L. F., 2001, ApJ, 551, L143
37. Lang C. C., Johnson K. E., Goss W. M., Rodríguez L. F., 2005, AJ, 130, 2185
38. Leitherer C., Chapman J. M., Koribalski B., 1995, ApJ, 450, 289
39. Leitherer C., Chapman J. M., Koribalski B., 1997, ApJ, 481, 898
40. Leitherer C., Robert C., 1991, ApJ, 377, 629
41. Linder N., Rauw G., Pollock A. M. T., Stevens I. R., 2006, MNRAS, 370, 1623
42. Linder N., Rauw G., Sana H., De Becker M., Gosset E., 2007, A&A, 474, 193
43. Linder N., Rauw G., Martins F., Sana H., De Becker M., Gosset E., 2008, A&A, 489, 713
44. Moffat A. F. J., Shara M. M., 1986, AJ, 92, 952
45. Montes G., Pérez-Torres M. A., Alberdi A., González R. F., 2009, astro-ph 0810.5026
46. Morrell N. I., et al., 2001, MNRAS, 326, 85
47. Nazé Y., Antokhin I. I., Sana H., Gosset E., Rauw G., 2005, MNRAS, 359, 688
48. Niemela V. S., Gamen R., Morrell N. I., Giménez Benítez S., 1999, in K. A. van der Hucht, G. Koenigsberger, P. R. J. Eenens, eds, Proc. IAU Symp. No. 193, “Wolf-Rayet Phenomena in Massive Stars and Starburst Galaxies”, p. 26
49. Nugis T., Crowther P., Willis A., 1998, A&A, 333, 956
50. Panagia N., Felli M., 1975, A&A, 39, 1
51. Parkin E. R., Pittard J. M., Corcoran M. F., Hamaguchi K., Stevens I. R., 2009, MNRAS, 394, 1758
52. Pauldrach A. W. A., Puls J., Kudritzki R. P., 1986, A&A, 154, 86
53. Persi P., et al., 1988, ASSL, 142, 227
54. Pittard J. M., 2007, ApJ, 660, L141
55. Pittard J. M., 2009a, MNRAS, 396, 1743 (Paper I)
56. Pittard J. M., 2009b, MNRAS, in preparation (Paper III)
57. Pittard J. M., Dougherty S. M., Coker R. F., O’Connor E., Bolingbroke N. J., 2006, A&A, 446, 1001
58. Pittard J. M., Dougherty S. M., 2006, MNRAS, 372, 801
59. Pittard J. M., Stevens I. R., 1997, MNRAS, 292, 298
60. Puls J., Vink J. S., Najarro F., 2008, Astron. Astrophys. Rev., 16, 209
61. Rauw G., Vreux J.-M., Stevens I. R., Gosset E., Sana H., Jamar C., Mason K. O., 2002, A&A, 388, 552
62. Runacres M. C., Blomme R., 1996, A&A, 309, 544
63. Rybicki G. B., Lightman A. P., 1979, Radiative processes in astrophysics (New York, Wiley-Interscience, 1979. 393 p.)
64. Sana H., Stevens I. R., Gosset E., Rauw G., Vreux J.-M., 2004, MNRAS, 350, 809
65. Schmid-Burgk J., 1982, A&A, 108, 169
66. Schnerr R. S., Rygl K. L. J., van der Horst A. J., Oosterloo T. A., Miller-Jones J. C. A., Henrichs H. F., Spoelstra T. A. Th., Foley A. R., 2007, A&A, 470, 1105
67. Scuderi S., Panagia N., Stanghellini C., Trigilio C., Umana G., 1998, A&A, 332, 251
68. Setia Gunawan D. Y. A., Chapman J. M., Stevens I. R., Rauw G., Leitherer C., 2003, in K. A. van der Hucht, A. Herrero, C. Esteban, eds, Proc. IAU Symp. No. 212,“A Massive Star Odyssey, from Main Sequence to Supernova”, p. 230
69. Stevens I. R., 1995, MNRAS, 277, 163
70. van der Hucht, K. A., 2001, New Astronomy Review, 45, 135
71. Van Loo S., Blomme R., Dougherty S. M., Runacres M. C., 2008, A&A, 483, 585
72. Wang Q. D., Dong H., Lang C., 2006, MNRAS, 371, 38
73. Williams P. M., van der Hucht K. A., Sandell G., Thé P. S., 1990, MNRAS, 244, 101
74. Wright A. E., Barlow M. J., 1975, MNRAS, 170, 41
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters