A SYSTEMATIC SEARCH FOR COROTATING INTERACTION REGIONS IN APPARENTLY SINGLE GALACTIC WOLF–RAYET STARS. I. CHARACTERIZING THE VARIABILITY

Stellar rotation has now been fully incorporated into evolutionary models of massive stars (e.g., Meynet & Maeder 2000; Heger et al. 2000), prompted by several observational facts that were unexplained by previous models without rotation (e.g., Meynet et al. 1994; Woosley et al. 1993) and that pointed to a solution related to mixing. For example, Herrero et al. (1992) had found that many O stars, and in particular fast rotators, show an enhanced He abundance; they called this the helium discrepancy. The models including rotation have been very successful in reproducing the observed abundances, predicting higher He and N abundances for higher masses and higher rotation velocities. Evolutionary models with rotation have also provided a solution for many other previously unexplained observations such as the well known blue-to-red supergiant ratio (B/R) problem (e.g., Langer & Maeder 1995). Indeed, when neglecting rotation, evolutionary models were unable to reproduce the rapid increase in B/R with increasing metallicity. As the new models predict many more red supergiants at low metallicity coupled with a shorter blue supergiant phase, they provided a natural solution.

It is now clear that rotation is an important factor in our understanding of massive stars and their evolution. It is therefore interesting to compare theoretical predictions with the available observational data. For O stars, the latest homogenous, large-scale surveys of projected rotational velocities have been carried out by Penny (1996) and Howarth et al. (1997), who respectively studied 177 and 373 stars, the latter also including early-type B stars. These studies were both based on high-resolution International Ultraviolet Explorer (IUE) spectra and used a cross-correlation technique with a narrow-line template spectrum to measure the rotational speeds. Both studies note the absence of very narrow lines for supergiants while some smaller velocities are found for main-sequence stars. Instead of interpreting this as a possible increase in rotational velocities as the star evolves, these authors favor an additional broadening mechanism such as macroturbulence. In fact, Howarth (2004) strengthened this assertion by pointing out that the shape of the absorption lines expected from pure rotation is a poor match to the observed profiles and that a line dominated by transverse turbulence is a much better match. Ignoring this possible effect (most likely small) of macroturbulence on the projected rotational velocities, their distribution peaks at ∼80 km s⁻¹ and has a very extended high-velocity tail reaching to 400 km s⁻¹. These values are not incompatible with model predictions; on the main sequence, the models show that even if the star starts out with a relatively flat internal rotation profile, the rotation becomes increasingly differential (but never more than a factor of 2 in velocity) because the convective core is contracting and therefore spinning up, while the outer layers are expanding and therefore slowing down. In spite of this, there is a global decrease in the rotation rate as a function of time as angular momentum is lost through mass loss. So, depending on their age and their initial rotation rate, the observed velocities are quite acceptable. Interestingly, Penny (1996) noted that, in fact, the higher rotational velocities are generally associated with lower mass stars. This is an important result. Indeed Meynet & Maeder (2000) predicted that the rotation rates decrease more rapidly for higher initial rotation rates and for higher mass and therefore higher mass-loss rate stars, in view of the increased loss of angular momentum. This agrees very nicely with the observations.

For WR stars, very slow average surface rotation velocities in the range 20–70 km s⁻¹ are predicted for the He-burning phase, depending on the initial rotation velocity on the main sequence and on the evolutionary state of the star (Meynet & Maeder 2003). The reason for these low values and the narrowness of this range is simple; the amount of mass that must be removed from massive O stars to turn them into WR stars is so big that most of the angular momentum is lost in the process. However, the differential rotation that steadily increased on the main sequence suddenly becomes very pronounced, as the helium core contracts at the end of core hydrogen burning. The star ends up with a core spinning faster than the initial zero-age main sequence (ZAMS) rotation velocity.

A fast-spinning core, at least in some WR stars, would support the presently favored model to explain long-soft gamma-ray bursts (GRBs). Indeed, the so-called collapsar model (Woosley 1993; MacFadyen & Woosley 1999; MacFadyen et al. 2001) involves the collapse of a rapidly rotating massive star in its presupernova phase. In this model, the central core becomes a black hole surrounded by a disk. The accretion of matter through the disk creates powerful jets in which the matter travels at very high speeds. Heger & Langer (2000) and Hirschi et al. (2004) have followed the evolution of a rotating He core to determine if the very rapid core velocities required by the collapsar model could be reached in spite of the many modes of transport and redistribution of the angular momentum taking place (mass loss, shear, turbulence, convection, etc). They found that indeed this is possible for initial rotation speeds on the main sequence that are typical of that of O stars.

However, adding magnetic fields to the models has brought new elements to the problem. Indeed, Heger et al. (2005) have shown that including magnetic torques caused by a Spruit–Tayler dynamo (Spruit 2002) dramatically reduces the angular momentum of the core leading to periods in the resulting pulsar of 15 ms at birth, which is compatible with the average value for common pulsars but much too small to make a long-soft GRB.

Nevertheless, for a rotating star, efficient chemical mixing can occur and if the mixing timescale is shorter than the thermonuclear timescale, the star can remain quasi-chemically homogeneous (Maeder 1987; Langer 1992). In that case, the star lacks the stratified layers which are required for the development of the Spruit–Tayler dynamo (Spruit 2002). This lack of a core–envelope structure prevents the magnetic core–envelope coupling described above and the spin-down of the star's core can be avoided. Woosley & Heger (2006) have shown that for a rapidly rotating star on the main sequence, if the mass loss in the WR phase is reduced by a factor of 10, massive-star evolution can indeed lead to a core with sufficient angular momentum to produce a long-soft GRB. The necessary reduced mass loss in the WR phase favors low-metallicity environments for GRBs. Finally, Yoon et al. (2006) calculated evolutionary grids for massive stars including rotation and the effects of a Spruit–Tayler dynamo for different metallicities. Still adopting WR mass-loss rates reduced by a factor of 10, they predict that the production of long GRBs should be limited to metallicities smaller than Z∼ 0.004. This translates into a peak at a redshift of z = 4.

No reliable direct observation of WR rotation rates has ever been obtained. Attempts have been made to estimate them based on the width of (what the authors had assumed to be) photospheric absorption lines in WR spectra (Massey 1980; Massey & Conti 1981), but these have been seriously questioned, as wind expansion and turbulence have not been taken into account.

It may be virtually impossible to obtain the rotational velocity of WR stars directly through photospheric line profiles but it could still be possible to deduce it from periodic variability in photometry or of spectral lines formed in the winds. Indeed, the winds of massive stars have been demonstrated to be highly variable. In particular, one type of structure in the wind that accounts for some of these variations, corotating interaction regions (CIRs), are thought to be closely linked to the rotation of the star. Indeed, these structures are caused by perturbations at the base of the wind that propagate through it and are carried by rotation. This generates a spiral-like structure in the density distribution, analogous to those commonly observed in the solar corona (Hundhausen 1972), which very likely leads to the ubiquitous discrete absorption components (DACs) in the spectra of O stars (Cranmer & Owocki 1996). These are rarely observed in WR P Cygni profiles as their absorption components are usually saturated (one exception is WR 24; Prinja & Smith 1992). However, it is believed that they also lead to a very characteristic, large-scale, periodic variability pattern in WR-wind emission lines. Dessart & Chesneau (2002) showed the unambiguous S-shape pattern in time series of model spectra of WR emission lines produced as a consequence of these variations. Observational support has been obtained for this type of variability in two clear cases through repeated spectroscopic observations. WR 6 (St-Louis et al. 1995; Morel et al. 1997; P = 3.76 days) and WR 134 (McCandliss et al. 1994; Morel et al. 1999; P = 2.25 days) show periodic variations without any indication of a companion. The correspondence between the data and models is striking when one considers that optical depth effects, which are certainly occurring in the observed data, have not been included in the models. This brings strong support to the idea that CIRs do in fact trace the rotation rate of the star at the base of the wind. According to Hamann et al. (1988, 1995), these stars both have R_* ≃ 3 R_☉. This radius corresponds by definition to a Rosseland optical depth of 20. Assuming that CIRs form at this radius, the equatorial rotation speed (at R_*) is ∼40 km s⁻¹ for WR 6 and ∼70 km s⁻¹ for WR 134, not too different from values predicted by evolutionary models. The periods and variability patterns were found to be the same for lines formed at various distances in the wind, indicating that the CIRs are governed by angular momentum conservation which results in a coherent structure in the wind. The period therefore provides a direct measurement of the rotation rate of the underlying star. To calculate the rotation velocity we need to know at what radius the perturbation leading to the CIR originates. If this radius is known, we can then calculate the rotation velocity at that radius. Much more observations and analysis of variability associated with CIRs are required to reach this goal. Therefore, the above velocities are only given as rough estimates for illustrative purposes.

In an attempt to identify more WR stars displaying CIR-type variability, we have carried out a survey of all Galactic WR stars brighter than v ∼ 12.5. Our complete sample includes 39 southern stars, which will be discussed in A.-N. Chené, & N. St-Louis (2009a, in preparation), and 25 northern stars, which we discuss here. The goal of this paper is to present the different types of spectroscopic variability found in our magnitude-limited sample of stars in a qualitative but also quantitative manner. Our aim is to use these characteristics to establish some guidelines to help identify the spectroscopic variations associated with CIRs. Once identified, these stars can then be observed more intensively and, if present, periods can be determined. In our second paper, we will build on these basis and use the complete sample to determine the frequency of the various types of line-profile variability. Follow-up papers on individual stars that have been more intensively monitored are also planned.

Our northern sample includes all Galactic WR stars visible from the Mont Mégantic Observatory that are not confirmed as SB1 or SB2 binaries. Our targets are listed in Table 1 where the stars' name, spectral type, R.A., decl., and v are provided. When available, we also provide terminal velocities measured either from IUE absorption troughs (Prinja et al. 1990), from near-infrared line profiles (Howarth & Schmutz 1992), or from optical line profiles (Hamann et al. 2006). The last column indicates each star's variability status (see below). Note that although some stars in our sample have published spectral types that indicate a possible companion (not indicated in the table), none have confirmed orbits. In the case of WR 4 and WR 128, only photometric variations have been found without radial velocity confirmation of binarity. In the case of WR 131, WR 143, WR 156, and WR 158, the claim of binarity is based solely on diluted emission lines (d.e.l.) and/or the presence of absorption lines in the spectrum (+a). However, Hamann & Gräfener (2007) argued that binary status established from the d.e.l. or +a criteria is often misleading as WR stars can produce a wide variety of spectra depending on their physical parameters. Those authors have demonstrated that in many cases these types of spectra can be reproduced entirely with a single-star model. In particular, WR 108, WR 156, and WR 158 were among the stars that they reproduced with a single-star model. Our final sample consists of 16 WN stars and nine WC stars. We have included WR 134 for comparison even though it is already known to show variability thought to be associated with CIRs. However, we have excluded WR 157 because it is too close on the sky to another star that is only ∼1 magnitude fainter. The separation between the two stars is ∼22, i.e., a value comparable with the average seeing at the Mont Mégantic Observatory, where our data were obtained. It was therefore impossible to isolate the WR star in our 15 spectrograph slit.

Our observations were obtained at the 1.6 m telescope of the Mont Mégantic Observatory in Québec, Canada, during six runs in 2001 and 2002. Details are given in Table 2 in which the run number, the dates of observations, and the Julian Dates (JD) are listed. We used a 830.8 l mm⁻¹ grating giving a resolution of Δλ = 1.6 Å (3 pixel) combined with a 2048 × 4096 EEV CCD.

Our goal is to identify new large-scale variability within our sample of stars. The extra features we are looking for on top of the main emission lines are quite large and therefore easy to find. Based on our previous experience with WR 6 and WR 134, we have made the assumption that a small number of good quality spectra, well sampled in time, are sufficient to identify new candidates showing variability possibly associated with CIRs. Of course, to quantify the changes and determine the period, a more thorough follow-up will subsequently be required. Therefore, we set out to obtain 4–5 spectra per star in our sample.

The spectra were reduced using the IRAF software package in the standard way. First a bias was subtracted, then each image was divided by a flat field. After extracting the spectra, the wavelength calibration was carried out using a He–Ar calibration lamp. Although our original spectral range was 1500 Å wide, we were forced to reject the reddest part of the spectrum due to a lack of usable lines in the calibration spectra. The final wavelength interval is Δλ = 4500–5200 Å. The signal-to-noise ratio of the spectra ranges from 100 to 120.

As our spectra are not photometrically calibrated, the final step of our data reduction procedure was to rectify the spectra. To do so, we have fitted a low-order Legendre polynomial to wavelength regions with no strong spectral lines and divided our spectra by the fitted curves.

The final spectra are shown in the top panels of Figures 1–13 for the two strongest lines present in our wavelength interval. For WN stars these are He ii λ4686 and He ii λ4860 while for WC stars they are C iii λ4650 and C iv λ5016. For each star we plot the differences between individual spectra and the global mean (which is shown in the second panel from the top). The differences have been shifted vertically for clarity. The Heliocentric Julian Date (HJD) of each observation is indicated on the left-hand side of the plot and the scale factor of the ordinate is indicated in the top right-hand corner of each panel. The graphs are organized by spectral type starting with WN stars from early to late types (Figures 1–8) followed by WC stars from early to late types (Figures 9–13). Note that for WR 4, WR 111, WR 5, and WR 154, we do not show the plot of the C iii λ4650 line because in those cases the line was on the extreme edge of the wavelength interval and therefore there was not enough continuum on one side of the line to allow for an accurate enough spectrum rectification. This could lead to a false detection of line variability. We prefer to adopt a conservative approach in this respect.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to search for significant spectroscopic variability among our targets, we have adopted the formalism of Fullerton et al. (1996). First, we have calculated the temporal variance spectrum (TVS), which is a statistically rigorous means of identifying the variable parts of our spectra. For N spectra, the TVS at wavelength j is defined by

$$\begin{eqnarray} &&({\rm TVS})_j={\frac{1}{(N-1)}}{\sum _{i=1}^N{\frac{w_i}{\alpha _{ij}}(S_{ij}-\bar{S}_j)^2}}, \end{eqnarray} \tag{ 1 }$$

where the weights of the individual spectra, w_i, are given by $w_i=(\frac{\sigma _0}{\sigma _{ic}})^2$ with

$$\begin{equation*} {\sigma _{0}}^2=\left(\frac{1}{N}\sum ^N_{i=1}{{\sigma _{ic}}^{-2}}\right)^{-1}, \end{equation*}$$

a standardized dispersion. Here, $\bar{S}_j$ is the weighted mean spectrum, σ_ic is the noise in the continuum and α_ij are wavelength correction factors given by $\alpha _{ij}=(\frac{\sigma _{ij}}{\sigma _{ic}})^2$, σ_ij being the noise associated with the flux S_ij.

Then, we have established the level of significant variability assuming that our data are governed by a reduced χ²-squared distribution with N−1 degrees of freedom. Using the standardized dispersion, σ₀, for scaling, we therefore adopt a threshold of σ²₀χ²_N−1(99%). Finally, we used $\Sigma _j(99\%)=\sqrt{\frac{({\rm TVS})_j}{\sigma _0^2 \chi ^2_{N-1}(99\%)}}$ to quantify the level of variability at each wavelength: a spectrum that reaches a level of 1, is variable at a 1σ level, if it reaches 3 it is variable at a 3σ level, etc. We have chosen a confidence level of 99%, which means that if the value of Σ_j(99%) reaches 3, we are 99% confident that the spectrum is variable at a 3σ level. We note that as the standardized dispersion, σ₀, is calculated using the noise level in the continuum regions and as the spectra are rectified, there is no real information on the variability of the continuum regions. In the bottom panels of Figures 1–13, we have plotted the Σ(99%) spectrum for each star for easy comparison with the spectra plotted above.

Inspection of Figures 1–13 reveals that the majority of the stars in our sample show significant variability in at least one of the two main spectral lines in their spectrum. One notable exception is the WN2 star WR 2. For WR 4, WR 5, WR 111, and WR 154, we did not detect variability in the weak C iv line, but we only have limited information at our disposal since we do not have usable data in the strong C iii line. It is clear that the level of variability is quite different from one star to another. The variations range from a few σ's to close to 60σ's in WR 115. As the quality of our spectra is relatively uniform, we can, in principle, compare the σ from star to star, but its value only provides information on the significance of the changes with respect to the noise level of our data.

As the line strengths are very different from star to star, comparing the significance level of the variations does not provide useful physical information. A quantity that is potentially more interesting is the fraction of the line flux that is variable. To characterize this, it is best to use a quantity that is representative of the variance in pixel j compared to the noise at that wavelength (σ_ij) rather than use a weight that is linked to the noise in the associated continuum, as the TVS does. To do this, Chené (2007) has defined a modified TVS as follows:

$$\begin{eqnarray} &&({\rm TVS}_{\rm mod})_j=\frac{1}{(N-1)} \sum _{i=1}^N{\left(\frac{\sigma _{0j}}{\sigma _{ij}}(S_{ij}-\bar{S}_j)\right)^2}, \end{eqnarray} \tag{ 2 }$$

where σ_0j is the reciprocal of the rms of the signal-to-noise ratio at pixel j of the spectral line series.

In order to quantify what fraction of the line flux (without the continuum) is variable, we divide this value by the line flux, $(\bar{S}_j -1)$. Of course this quantity has no meaning outside emission lines and in particular in continuum regions. To avoid any problems, we therefore limited our calculations to a width corresponding to a relative intensity level of 1.3 of each line. Also, for very weak lines (<1.3 in relative intensity), our ratio becomes meaningless because the line flux, the denominator, is too close to zero and the variability level, the numerator, is dominated by the noise. In such cases, we have left the panel blank. Finally, absorption components of P Cygni profiles are also excluded since in that case the value of $(\bar{S}_j -1)$ will be negative and cannot be compared directly with what is observed in emission lines.

The values of $\sigma _j=({\rm TVS}_{\rm mod})^{1/2}_j/(\bar{S}_j -1)$ are plotted in the middle panels of Figures 1–13. Several interesting conclusions can be drawn from these plots. First, for all types of stars, both lines show roughly the same level of variability. However, one has to note that in the case of WN stars, we compared lines of the same ion and therefore it is likely that in a given star, these two lines are formed in similar regions of the wind. The line widths are indeed, in most cases, very similar. Another interesting feature of the variability is that the calculated values of Σ are much stronger in P Cygni absorption components than in the corresponding emission components, i.e., the variations are more significant. This is perhaps not surprising since that part of the profile comes from a much smaller volume of the wind than the emission and therefore suffers less from canceling effects from large-scale changes in density or ionization arising in different parts of the wind. This effect is clearly present in the P Cygni profile of WR 123 (He ii λ4686) or WR119 (C iii λ4650), for example.

After a careful examination of Figures 1–13, we divide our sample into three main categories: stars showing no profile variability (NV), stars showing small-scale profile variability (SSV – σ < 5%), and stars showing large-scale variability (LSV – σ>5%). Note that the stars designated as NV are not necessarily stars that do not show spectral variability, but rather that no significant variability can be detected at a level of confidence of 99% on spectra having signal-to-noise ratio ∼100. Our conclusions are listed in the last column of Table 1. The detected subpeaks have widths that range from ∼140 km s⁻¹ for the smallest isolated structures to ∼350 km s⁻¹ for the largest. Unfortunately, our limited data set does not allow us to search for correlations between these widths and, for example, the level of variability or the width of the global line. In the following sections, we will discuss in turn the variability we found in WN and WC stars.

3.1. WN Stars

3.2. WC Stars

The results of our search for optical spectroscopic variability among presumably single WR stars can be summarized as follows. For eight stars, WR 3, WR 152, WR 128, WR 136, WR 131, WR 158, WR 108, and WR 135, variability was found but at a small level compared to the flux of the line. This type of variability is most likely associated with inhomogeneities in the dense WR wind, which move with it and reveal themselves as small-scale structures superposed on the broad emission lines. In more detailed studies, they have been found to move from line center to the edges of the lines as described for example by Moffat et al. (1988) and Lépine et al. (1996). As the driving mechanism for these radiatively driven winds has been found to be inherently unstable (Owocki et al. 1988), such variability is expected for all WR stars.

Nevertheless, for five other stars in our sample, namely WR 2, WR 4, WR 111, WR 5, and WR 154 we found no significant line-profile variability. Note however that in the last four cases, our conclusion is based solely on a single relatively weak line. In fact, WR 111, is a star that is known to display small-scale spectral variability due to the presence of inhomogeneities in its wind and has been studied by Lépine & Moffat (1999). In any case, identifying a star as nonvariable always depends on the quality of the data in hand. Down to what level of accuracy these stars are constant is, indeed, an important question. One would like to verify if all WR stars have inhomogeneities in their winds, if some stars present a different level of variability and if so, what is the cause of this difference. This kind of information can help us better understand the physical mechanisms responsible for these unstable outflows.

Nine stars in our sample presented large-scale profile variability. Of these, we have identified three new candidates for CIR-type variability: WR 1, WR 115, and WR 120. Follow-up observations with a high temporal resolution will be required to determine if these changes are periodic and if so identify the period. We have done so already for WR 1 and present our results in an upcoming paper (A.-N. Chené, & N. St-Louis, 2009b, in preparation).

Although the results presented here are not yet for our complete sample, it is interesting to compare our findings to what is known for OB stars. Although the presence of DACs in O-star P Cygni profiles is found to be common (e.g., Howarth & Prinja 1989), cyclical variability in emission lines is not found to show the same ubiquitous behavior. Morel et al. (2004) presented a long-term monitoring campaign of the Hα line of 22 bright OB supergiants. Variability was found in all cases with two possible origins: either the underlying photospheric profile and/or a wind component. Cyclical changes which very likely originate from large-scale structures in the wind were found for only two stars in their sample (∼14% of the stars). In this study we have identified four CIR-type variables (three of which are new) out of a sample of 25, or ∼16%, a very similar fraction. Alternatively, we can exclude the WN 8 and WC 9 stars from this estimate; if we do so we also obtain a similar fraction, 4/19 ∼ 21%. This striking difference between variability in P Cygni absorption components and emission lines is most likely due to the fact that in the former we are sampling a much more physically restricted volume of the wind (i.e., that in front of the star) compared to the entire wind which can be subject to canceling effects in parts of the wind which are physically distinct but that are traveling at a similar projected velocity.

In addition to the stars showing variability thought to be associated with CIRs, we have defined two distinct classes of LSV: WN8 stars which are already known to be a peculiar class of WR stars and WC9d stars. In the latter case, it remains to be checked if the changes originate in the wind of a single star because of the presence of CIRs or if these stars are in fact colliding-wind binaries. If of binary origin, a very similar variability pattern is expected with some differences that can be identified with a detailed monitoring campaign. For example, strong epoch-dependence in spite of the presence of periodic changes would suggest a single-star origin. Also, shock cones translate into a very distinct and easily recognizable variability pattern on spectral lines, translating to at most, two large excess emissions moving around periodically on the emission line.

Although our survey has revealed a relatively small fraction (∼12%) of stars showing line-profile variability thought to be associated with CIRs among WR stars, it still provides some new candidates for potential estimates of WR rotation velocities. It also provides us with the opportunity to study in more detail the CIR phenomenon for which still very little is known by providing a larger pool of candidates to observe. It is clear that better observational constraints of this phenomenon are crucial to reach a greater understanding of the physical mechanisms from which they originate.

We thank G. Skalkowsky, C. Foellmi, and G. Caron for participating in observing campaigns in 2001/2002. We are also grateful to the (anonymous) referee for comments that led to an improvement of the quality of this paper. Finally, we thank the National Sciences and Engineering Research Council (NSERC) of Canada for financial support.