A DETAILED X-RAY INVESTIGATION OF ζ Puppis. II. THE VARIABILITY ON SHORT AND LONG TIMESCALES^*

Light curves of the source and background were extracted using the SAS task epiclccorr for four time bins (200 s, 500 s, 1 ks, and 5 ks) and in seven energy bands: total (0.3–4 keV), Berghöfer's band (0.9–2 keV), soft (0.3–0.6 keV), medium (0.6–1.2 keV), hard (1.2–4 keV), very hard (4–10 keV), and grand total (0.3–10 keV). The background and source regions are the same as for the spectral analysis (see Paper I).

The choice of the energy bands results from a compromise, taking into account the appearance of the ζ Puppis spectrum, count rates in each band, and previous results. Indeed, one has to test bands for which results were reported before. As explained in Section 2, a periodicity was found using ROSAT in the 0.9–2 keV band. This is why we define a "Berghöfer's band," to be able to confirm (or not) the presence of this signal in our much higher quality data. Another obvious choice is the total band, which enables us to use all available counts, and therefore get the smallest error, performing the most sensitive test of the overall level of the X-ray emission. Additional narrower bands (soft, medium, hard) are defined by taking into account the general spectral energy distribution (Figure 1). The soft and medium bands probe the brighter part of the spectrum (with count rates above ∼0.6 counts s⁻¹ keV⁻¹), while the hard band enables us to investigate the tail of the spectrum. EPIC spectra showed a dim region around 0.6 keV, hence the choice of this energy to split soft and medium energy bands. Count rates are rather similar for these "narrow" bands: 0.8, 0.4, 1, and 0.4 counts s⁻¹ for MOS (respectively 2, 2, 3.5, and 1 counts s⁻¹ for pn) in the Berghöfer, soft, medium, and hard bands, respectively. The spectra appear very noisy above 4 keV, hence the choice of that energy as an upper limit for the main energy bands. Nevertheless, a (double) check was made as we want to ensure that including all data with >4 keV does not change the conclusions. Indeed, the grand total and total bands provide indistinguishable results. The very hard band would reveal the presence of very high energy phenomena (transient or not). For ζ Puppis, however, this band always presents a very low count rate, of the order of ∼(10 ± 3) × 10⁻⁴ counts s⁻¹ for MOS (1 or 2), i.e., three orders of magnitude below the total band count rate, and similar to the background count rate in the same band for pn. This means that no bright emission at very high energies is present in this star, and we will therefore focus the analysis of ζ Puppis on the five main energy bands.

Zoom In Zoom Out Reset image size

The choice of the time bins was also made by taking several factors into consideration. First, we considered the problem of sampling different timescales. As one wished to study ζ Puppis on as many timescales as possible, both short bins and longer bins are highly desirable. However, the upper limit on the time bin length is dictated by the length of the individual exposures (13 ks for the shortest ones) since at least a few bins are required for a meaningful analysis within the exposure itself. Second, we looked at the impact of errors. Our goal is to detect the most subtle variations of ζ Puppis. To this aim, we need to combine the smallest possible errors. Considering the count rates for the different bands and instruments (see above), we can reach a 5% error on MOS data in 200 s for the total band, in 500 s for the medium and Berghöfer bands, in 1 ks for the soft and hard bands; a 3% error on pn data in 200 s for the total band, in 500 s for the soft, medium, and Berghöfer bands, in 1 ks for the hard band. A time bin of 5 ks has been added to get 1% or less error on the total band and enhance the detection of longer-term changes within one exposure.

To check the influence of the background (which is much fainter than the source in all observations), three sets of light curves were produced and analyzed individually: the raw source+background light curves, the background-subtracted light curves of the source, and the light curves of the sole background region. The results found for the raw and background-subtracted light curves of the source are indistinguishable. An example of a light curve set is given in Figure 2.

Zoom In Zoom Out Reset image size

Finally, two remarks must be made. To avoid very large errors and bad estimates of the count rates, we discarded bins displaying effective exposure time <50% of the time bin length. Our previous experience with XMM-Newton has shown us that including such bins degrades the results. It should also be noted that the SAS task epiclccorr provides equivalent on-axis, full PSF count rates, so that problems such as the presence of an offset or of bad pixels were corrected to provide comparable light curves for all exposures.

The task rgslccorr yielded raw source+background light curves, background-corrected light curves for the source, and background light curves for 10 wavelength bins: total (6–30 Å, equivalent to 0.4–2 keV), Berghöfer's band (6–14 Å, equivalent to 0.9–2 keV), soft (20–30 Å, equivalent to 0.4–0.6 keV), medium (14–20 Å, equivalent to 0.6–0.9 keV), N vi (28–30 Å), N vii (24.3–25.5 Å), O vii (21.3–22 Å), O viii (18.75–19.2 Å), Ne ix (13–14 Å), and Ne x (11.8–12.5 Å). They were calculated for each instrument (with both orders combined) and for the whole RGS (with both instruments and both orders combined). Note that the wavelength of the O vii line is at a gap of RGS2 while Ne x is in a gap of RGS1 (order 1).

As for EPIC, there were obvious choices, such as the total band or bands covering the strongest lines, but also choices dictated by previous reports (Berghöfer's band), and compromises (soft and medium bands—the soft band is similar to its "twin," the EPIC soft band, while the medium band appears in between soft and Berghöfer's bands). The limits of the wavelength intervals were carefully chosen taking into account the global aspect of the high-resolution spectra (Figure 3). Limits of the bands fall in regions as free of lines as possible; the overall lower and upper limits are dictated by the increasing noise toward short and long wavelengths. It may be noted that count rates for the soft, medium, and Berghöfer bands are similar (about 0.5 counts s⁻¹). The bands covering lines were chosen to enclose the brightest isolated lines without enclosing too much of the neighboring continuum.

Zoom In Zoom Out Reset image size

Again, as for EPIC, the choice of time bins represents a compromise between getting a sufficient number of counts to perform an analysis with a few percent relative error only and having enough bins within each exposure (of length 30 ks for the shortest ones). Bins of 500 s and 2 ks yield errors in the total band, for each RGS, of 7% and 3.5%, respectively; bins of 2 ks and 5 ks give errors in the soft, medium, and Berghöfer bands of 7% and 3.5%, respectively; bins of 5 ks and 10 ks ensure errors of 10% and 5%, respectively, for bands linked to specific lines.

Note that, contrary to epiclccorr, the SAS task rgslccorr does not correct for lower recorded fluxes if the source appears off-axis. While this does not appear to be a problem for small offsets (<15), the count rate of ζ Puppis during Rev. 0731 (where the offset reaches nearly 6') is clearly reduced, by about 25% in the total band. In addition, the count rate from Rev. 0091 at highest energies (shortest wavelengths) appears too high by about 20%. The origin of this problem is unknown (poor calibration at the earliest times of XMM-Newton operations or true brightening?—there is no EPIC data to confirm what happened). To avoid any interpretation problems, both data sets were discarded from global analyses (see Sections 4 and 5).

It is also important to remember that rgslccorr performs a randomization of the time tags, i.e., all runs of the same task, with the same parameters and input files, will produce slightly different results (the difference remaining within the error, of course). This is not the case for epiclccorr since a randomization is applied during the initial processing of the raw files (tasks epproc, emproc).

To search for the long-term variations on timescales of months to years, we computed the average count rates for each exposure and checked their constancy using χ² tests. We remind the reader of a few cautions. For RGS, the data points associated with Rev. 0731 are discarded from global analyses since the effective area changes due to the off-axis position of ζ Puppis are not taken into account in rgslccorr. The data from Rev. 0091 are also discarded since they yield suspiciously high count rates (see XMM user handbook⁵ and Section 3). For EPIC, the pn data taken with the Medium filter are most probably still slightly affected by pileup (see also Paper I). Indeed, they yield clearly erratic results: while Revs. 535, 538, and 542 provide very similar count rates for their Thick filter data, this is not the case for their Medium filter data.

Thereby considering only the best (i.e., most reliable) data, it clearly appears that the count rates are not stable over timescales of years (Figures 4 and 5). The light curves obtained by each detector decline. The rate of decline however differs among them. Over the ∼3800 days covered by the data sets, the EPIC-pn count rate decreased by about 6% in the total, soft and medium energy bands, 2% for the hard band and 4% for Berghöfer's band. The count rates of EPIC-MOS decreased by about 10% in the total, medium, and Berghöfer's bands, 6% in the hard band, and 15% in the soft band. The RGS count rates decreased by 18% in the total band, 28% in the soft band, 16% in the medium band, and 12% in Berghöfer's band. These variations are significant at the <1% level⁶ as the error bars on the average count rates are very small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fitting of EPIC spectra also showed a small decrease in flux of a few percent (Paper I). Moreover, the supernova remnant 1E0102–72, observed independently from ζ Puppis (though also for calibration purposes), shows a decrease in its MOS count rate and flux similar to what is seen for ζ Puppis (M. Guainazzi 2011, private communication). The exact origin of this systematic trend is not known (investigations are under way by the calibration teams), but such long-term effects are reminiscent of detector sensitivity ageing problems. Up to now, the pn and MOS are considered to be stable by the XMM-Newton calibration team (M. Guainazzi 2011, private communication),⁷ but there are apparently remaining imperfections in the long-term calibration. Since the decline of the registered count rates is currently not taken into account in the XMM calibration, we removed this trend by hand for further analyses and conclude that it is instrumental in nature.

Note that true variations of the source may be superimposed on these global decreasing trends. For example, ζ Puppis appears brighter during Rev. 1620 for both count rates and fluxes (derived from spectral fits, see Paper I). This is linked to the changes observed on intermediate timescales (see Section 5)—with trends on the course of days, the source may indeed appear brighter or fainter sometimes, depending on the observational sampling of these changes.

Finally, we also calculated the ratios between the RGS count rates associated with lines of the same elements: O vii/O viii, Ne ix/Ne x and N vi/N vii. Comparing globally the 16 secure RGS exposures (see above and Section 3), there is no obvious correlation between the different elements: the Ne ix/Ne x and N vi/N vii ratios do significantly vary, but only the N vi/N vii is clearly decreasing. The N lines constitute the most widely separated pair: as the long-term trend is instrumental in origin, such a calibration problem could affect some wavelengths more than others and a larger separation would then lead to larger differences in lines. The average values of these ratios, based only on count rates (i.e., not on dereddened fluxes) and excluding Revs. 0091 and 0731, are 1.44 ± 0.01 for nitrogen (both RGS, both orders), 1.28 ± 0.01 for oxygen (RGS1, both orders), and 1.86 ± 0.01 for neon (RGS2, both orders)—the choice of instruments takes into account the fact that RGS2 could not record any information on O vii, whereas the first order of RGS1 contains no data on Ne x (see Section 3).

High-resolution X-ray spectra of good quality can be derived for each exposure, and used to search for variations. The lines recorded for each observation can first be compared to the average spectrum obtained when combining all observations (Figure 6). No large, significant variation is detected. Small changes in flux (of the order of one sigma) can be spotted from time to time. They are similar to those recorded at optical wavelengths (Eversberg et al. 1998; Rauw et al. 2010), with double peaks where the blue to red ratio slightly changes. Their small amplitude and the low signal to noise of the data however prevent us from making a detailed analysis: this must await the advent of more sensitive observatories.

Zoom In Zoom Out Reset image size

To be more quantitative, a temporal variance spectrum (TVS) analysis (Fullerton et al. 1996) was performed on the combined and fluxed RGS spectra (output of rgsfluxer). The TVS computes the squared difference between individual spectra and the average one, taking the signal-to-noise ratios into account, which allows the detection of statistically significant deviations from the average. Note that the spectra are calibrated both in flux and wavelength within rgsfluxer (see Paper I for details), so that the pointing problems have no impact here. The relative weights for the different spectra in the TVS were chosen equal to the count rate errors on the total band. The TVS appears overall flat (implying constancy), with only a few distinctive features (Figure 7). For example, some peaks are found in the 20–24 Å region. This region is affected by gaps and missing CCDs, and there are therefore undefined values for wavelengths where no exposure exists. The exact position (in wavelength) where that occurs varies from one exposure to the next, notably because of pointing differences, and that results in apparent variability in the spectral set (hence peaks in the TVS). Increased noise at the shortest wavelengths (below 7 Å) and longest wavelengths (above 24 Å) is also producing a larger TVS, unrelated to intrinsic variability of the source (though the peak at 7.7 Å remains unexplained). Globally, the TVS thus indicates that the lines in the X-ray spectrum of ζ Puppis are not significantly varying from one exposure to the next. Very small changes are not excluded, however. Indeed, two small peaks of the TVS occur at 13.5 and 15 Å, i.e., at the wavelength of the strongest lines: this indicates that the line profile variations of ζ Puppis are just beyond the reach of current facilities, and they may thus be detected in the future with the better sensitivity of new X-ray observatories.

Zoom In Zoom Out Reset image size

The characteristic timescale of these intermediate-term variations is difficult to assess on individual exposures using χ² tests or BBs. These first analyses can only conclude that the timescale is longer than the typical exposure length (i.e., a day or more). Looking further at the individual light curves, no obvious oscillation is detected by eye, despite the full coverage of the putative ROSAT 17 hr period by individual exposures and the large improvement in quality over the ROSAT data. It thus seems that a sinusoidal variation of several hours and an amplitude of a few percent, such as that detected by Berghöfer et al. (1996), are transient, at best.

To give a more quantitative assessment of the variability timescale, we have performed period searches on global light curves (created by putting all individual-exposure light curves after one another, keeping their individual time tags), corrected for the long-term instrumental trend described in Section 4.1. We used the algorithms of Heck et al. (1985, see also correction by Gosset et al. 2001). As in previous sections, we considered only the best data sets (see Paper I), i.e., the 15 MOS exposures and the 10 pn exposures taken in small window + Thick filter mode, as well as the 16 RGS exposures (i.e., excluding Revs. 0091 and 0731). The small time bins were favored (200 s for EPIC, 500 s for the total RGS band, and 2 ks for the other RGS bands) since the Fourier-period-search algorithms are more affected by a small number of points than by large individual errors. Note that a more general period search technique which attempts to fit a period simultaneously with its harmonics does not yield additional insights.

The resulting periodograms are shown in Figure 10, along with the spectral window (i.e., the aliasing structure appearing naturally because of the temporal sampling). The periodograms consist of narrow peaks (since the data set covers more than 10 years) tightly disposed in much broader features (because there are only a few exposures, of max. 70 ks, within this 10 year timescale), which hamper an accurate determination of any period. However, two features are apparent for all bands except the hard one: (1) the highest peak occurs at about 0.3–0.4 day⁻¹ and (2) the right wing of this peak decreases much more slowly than that of the spectral window—some additional power is thus present around 0.7–1.3 day⁻¹.

Zoom In Zoom Out Reset image size

To assess the significance of these signals, we performed a formal iterative decomposition in frequencies—we here illustrate only the analysis for the pn total bandpass. A first extracted frequency at ν₁ = 0.365 day⁻¹ (period P = 2.74 day, semi-amplitude a = 0.25 counts s⁻¹) is significant against the null hypothesis of white noise as deduced from ad hoc simulations. The following frequencies are ν₂ ∼ 0.89 day⁻¹ (P = 1.12 day, a = 0.12 counts s⁻¹) and ν₃ ∼ 0.4 day⁻¹ (P ∼ 2.5 day, a = 0.08 counts s⁻¹; the last values are rather uncertain because they depend on the selected binning). These peaks in the Fourier periodogram are still characterized by a significance level lower than 0.001. The white noise hypothesis is thus formally rejected.

However, these frequencies do not correspond to the true content of the signal of the star. For example, when phased with the 0.365 day⁻¹ frequency, individual light curves appear one after another, with no common phase interval. As explained above, some individual runs exhibit shallow trends. The time spanned by these runs combined with their scarcity over the covered decade has an unfortunate consequence: the Fourier transform is always able to find a few frequencies that are combining the various trends in a very constructive way. However the plethoric presence of gaps in the time series and the concomitant large number of degrees of freedom is responsible for the positive combination more than the possible coherency of the signal. As a test, we detrended the individual runs with a linear function before merging them. We computed the periodogram for the new time series. All the previously reported candidate frequencies disappear from the list of frequencies. The dominant one, for the 200 s binning, is ν₄ = 1.42 day⁻¹ (P = 0.70 day, a = 0.054 counts s⁻¹) which corresponds to a peak whose significance level is around 0.001. If we detrend the X-ray light curves over the individual runs with a second degree polynomial, we obtain a time series that exhibits no outstanding peak and is thus in perfect agreement with the white noise hypothesis. Therefore, although the white noise hypothesis is rejected, the alternative hypothesis to adopt remains unclear. The star exhibits weak although significant variations at a time scale between 0.7 days and several days but the sampling does not allow us to conclude the existence of a coherent, fully deterministic signal above Poisson noise. The signal could be dominantly stochastic with some coherency at short timescale (some kind of red noise⁸). It is also possible that these changes are transient: they would then not appear as a strictly periodic signal in the global light curve. Another possibility is that any signal could be not strictly stationary, in phase or frequency, from one observing run to the next, as often seen for changes of the optical spectra in Oef stars (Rauw et al. 2003; it must be recalled that our target, ζ Puppis, belongs to this spectral type category!).

In summary, it is certain that ζ Puppis displays variations with relative amplitudes up to 10%–15% on the timescales of days—not hours (i.e., there is no trace of the signal attributed to non-radial pulsations; see Reid & Howarth 1996). However, our data set suffers from the misfit between the possible frequencies and the actual sampling. Therefore, despite its exceptional quality, it is still insufficient to detect confidently evidence for the presence of a coherent signal with daily timescales. A modulation of the X-ray flux with such daily timescale was reported for the O-type dwarf ζ Oph, and potentially associated with the corotating interacting features found through UV observations (Oskinova et al. 2001a). ζ Puppis shows a 5 day modulation also attributed to corotating regions (Moffat & Michaud 1981). New, continuous X-ray observations of ζ Puppis would be required to better characterize the timescale of the detected variability, and to establish whether these variations can be associated with such features.

To search for temporal links between the photon arrival times, the data corrected for the instrumental effect were also analyzed using autocorrelation methods (Edelson & Krolik 1988). The resulting autocorrelation functions are shown in Figure 11. Note that the binned functions were normalized by the actual number of points used and, as usual, by the dispersion around the mean (observed variance s²_ts). The former normalization is needed as the individual exposure times range from 15 to 77 ks (Paper I), so that more data points can be used for some time shifts.

Zoom In Zoom Out Reset image size

A slight positive correlation appears for the total and medium bands for time shifts below 20 ks. A slight anticorrelation is then detected for shifts of 50–60 ks, i.e., similar to the typical duration of an exposure. A similar behavior, though with even smaller amplitude, is seen for the soft, hard, and Berghöfer energy bands. This slow evolution toward lower correlation values confirms the existence of weak coherent variations with timescales of tens of ks.

The correlation amplitudes are however small, hence not highly significant, statistically speaking. However, a strong positive correlation for small shifts is expected since the wind configuration remains the same for some time (about one flow time, which is R_*/v_∞ = 5.8 ks for ζ Puppis, as v_∞ = 2250 km s⁻¹ and R_* = 18.6 R_☉, using the parameters of Oskinova et al. 2006 and references therein). For longer shifts, two behaviors are then possible. On the one hand, the time bins could be independent for shifts longer than a flow time, as in a fully stochastic wind. In this case, the correlation function would rapidly drop to zero. On the other hand, if the time bins are coherently linked by a shallow trend, the correlation function should slowly decrease. This decrease is indeed what is observed, but we would thus expect a correlation signal with larger amplitudes.

In practice, however, even a perfectly correlated signal (i.e., with a correlation of 1) will be diluted by the observational Poisson noise. The expected correlation value will then be given by

$$\begin{equation} C = 1 - \frac{s^2_{n}}{s^2_{ts}} \end{equation} \tag{ 1 }$$

where the observed variance s²_ts comprises the noise and the coherent variations of the source, whereas s²_n only corresponds to the noise. For the pn light curve in the total band with 200 s time bins (globally detrended, i.e., corrected for the instrumental effect), s²_ts is 0.089 counts² s⁻² and the variance due to Poisson noise is estimated to be around 0.049 counts² s⁻². The peak height C should then reach 0.45, which is in good agreement with the value observed in Figure 11 for small time shifts. This implies that the intrinsic correlation of the noiseless signal could be nearing one: the trends are thus real.

As already mentioned in Section 3, RGS light curves were also extracted for the brightest isolated X-ray lines, enabling us to determine whether these lines vary over shorter timescales. The observed variations are small, generally within 1σ error bars. The only significant and coherent results are found for Revs. 0156 and 1814. For the former, a linear trend with a positive slope provides a much better fit to the Ne ix line flux. For the latter, the Ne ix line appears variable while the O viii clearly increases. Note however that these changes are detected with more significance using RGS2 data than with RGS1 data.

Since the individual count rates are available, we also calculated the ratios between lines of the same elements: O vii/O viii, Ne ix/Ne x, and N vi/N vii (Figure 12). These ratios reflect the temperature (Blumenthal et al. 1972; Waldron & Cassinelli 2007), and any change in the ratios would therefore be linked to temperature variations. Formally, opacity variations may also play a role in changing these ratios. However, the latter are not as sensitive to absorption as to temperature: a change in temperature by 15%–20% results in a doubling of these ratios (note that the temperature dependence is similar for all ratios), whereas change of the absorbing column of a smooth wind by a factor of two (from 0.1 to 0.2 × 10²² cm⁻²) yields variations of 7.5%, 25%, and 35% in the Ne, O, and N ratios, respectively. As for count rates, the ratios were tested using χ² tests. For the individual exposures, only two features appear significant. For Rev. 0156, the Ne ix/Ne x is better fit by a line with negative slope. Similar but shallower possible trends are detected by eye for the two other ratios, but they are not formally significant. For Rev. 1814, the Ne ix/Ne x appears much better fit by a concave/U-shaped parabola, but no similar signal is seen for the two other ratios. Note that these two revolutions are amongst the variable ones (see above), and that the reported changes may be linked to slight changes in spectra (see Section 4.2), for which higher-sensitivity instruments are needed for a detailed characterization.

Zoom In Zoom Out Reset image size

Relative dispersions were calculated for each of the observed light curves. These light curves have two specific features which do not exist in the models (see next subsection): a temporal binning and a binning over a range of energies. In principle, this somewhat smoothes out any variability if it is present with timescales shorter than the time-bin length or if it changes strongly with energy (e.g., the light curves at 0.3 and 0.4 keV being uncorrelated). However, the Poisson noise inevitably impacts on the data and therefore prevents us from detecting low-level variability (i.e., a few percent) with time bins smaller than 200 s (see last line of Table 1)—this limit would only be changed if a more sensitive instrument was used. In addition, the emitting parcels contain plasma emitting over a range of energies, even if isothermal, and this implies some correlation over different ranges of energies. A comparison with synthetic curves can thus provide useful insights into the structure of the wind.

As mentioned below, the steps in the synthetic light curves correspond to different wind configurations. The dispersion calculations for the data were thus performed in the following way: (1) the original light curves of each exposure were first detrended using the best-fit linear trend derived from χ² calculations (see above) as the simple model (presented in next section) does not predict nor model such features; (2) the count rate values of the 200 s light curves were then sampled at 5 and 10 ks intervals (corresponding to about one and two wind flow times, respectively), i.e., only considering every 25th or 50th value; and (3) the relative dispersions of these new, reconstructed light curves were finally evaluated using

$$\begin{eqnarray*} &&{\rm rel.\ disp.}=\left[ \sqrt{ \sum (CR_i-a-b\times t_i)^2/(N-2)}\, \right] \, \right / {\rm mean,} \end{eqnarray*}$$

where a and b are found from the best-fit trend, see Section 3.2.1, and mean = ∑(CR_i/σ²_i)/∑(1/σ²_i). This enables us to try to reproduce the "independent wind configurations" of the successive "time steps" in the synthetic curves. Dispersions were also calculated for the full (i.e., considering all time bins) 200 s, 500 s, 1 ks, and 5 ks light curves (after detrending). When one compares the dispersions of the full 200 s light curves to these of the reconstructed light curves, they appear very similar. The maximum differences amount to ±4% in the worst cases, or smaller (<1% in absolute value) when the observations are longer. This is unsurprising since dispersion estimates on, e.g., 2 or 3 bins are less precise than on, e.g., 50 bins: when one samples the light curves every 25 or 50 steps, a particularly discrepant realization of the noise has more impact when there are fewer bins. The sparse sampling therefore does not change the results much, so Table 1 reports only the dispersions of the full light curves. This table gives in the first column the revolution number, in Columns 2–6 the dispersions measured for the main five energy bands and the 200 s binning, and in the 7th column the number of data points used; Columns 8–13 report similar values, but for the 5 ks binning. The last line provides for comparison the expected Poisson error, relative to the mean (i.e., $1/ \sqrt{{\rm count\,rate} \times {\rm bin\,length}}$). Note that only the EPIC data were used here, as the noise of the RGS light curves is even larger and would more easily mask subtle variations.

Looking at Table 1, one thing is obvious: despite the fact that the exposure-long trends mentioned above are often not perfectly linear, the measured relative dispersions around the best-fit lines are close to that expected on the sole basis of the Poisson noise. This means that the true short-term variability of ζ Puppis—and its potential wavelength dependence—remains hidden in the noise, implying a very small amplitude for these "intrinsic" short-term changes.

These results should now be translated into constraints on the wind structure. Considering X-rays to be emitted by many hot parcels, one may naively think that it is sufficient to notice that a 1% dispersion "naturally" corresponds to 10⁴ "emitters," using simple Poisson statistics. However, reality is not as simple: emitters suffer from different amounts of absorption depending on their location and, worse, the absorption may be clumped too. In both cases, the simple reasoning completely fails to apply, and a dedicated model is thus needed.

Embedded wind shocks are considered to be responsible for the X-ray emission of O-type stars. In these expanding and unstable winds, fast material encounters slow-moving material, giving rise to zones of dense gas, and the mutual collisions of such dense features then give rise to the X-rays. The radiation hydrodynamic (RHD) models (e.g., Feldmeier et al. 1997a, 1997b) predict that the collision, and thus the heating, occurs quickly (tens of seconds), and the cooling of the resulting hot plasma is also rather rapid, as the cooling time remains lower than the wind flow time for distances of several tens of stellar radii. A strong X-ray variability (two orders of magnitude in amplitude) was predicted. However, early ROSAT and Einstein observations failed to show such a strong variability of X-ray fluxes from massive stars (Berghöfer & Schmitt 1994a, 1994b, and references therein). To reconcile the results of the early observations with wind-shock theory, it was suggested that lateral breakup of the shells can lead to the presence of many parcels, resulting in low variability (Cassinelli et al. 1983; Feldmeier et al. 1997a): therefore two-dimensional (2D) or 3D models were needed. First attempts for 2D wind models were made by Dessart & Owocki (2003, 2005) but only isothermal winds (i.e., without X-ray generation) have so far been considered. Another example is the recent work of Cassinelli et al. (2008) on bow shocks, but the model is not yet self-consistent in terms of production and evolution for an ensemble of clumps. Unfortunately, no new, multi-dimensional hydrodynamical model has thus tackled the problem of X-ray generation in stellar winds, and we are thus left with alternative modeling paths.

In this context, Oskinova et al. (2001b) considered a spherically symmetric smooth cool wind permeated with discrete zones of hot X-ray emitting gas. It was found that X-ray variability depends on the frequency with which hot zones are generated, and on the cool wind opacity for the X-rays. It was shown that in such smooth cool wind, the variability in the soft band is expected to be smaller than in the hard band. Oskinova et al. (2004) further developed a 2D wind model where not only are hot parcels discrete, but the cool absorbing wind can be clumped too. We apply this model in the present paper to investigate the resulting X-ray variability. We only briefly recall here its basic features.

Our model was designed to correctly perform the radiative transfer of X-rays in stellar winds, and it also reproduces the features derived from RHD simulations. In this model, the X-ray emission originates from discrete zones of hot gas randomly distributed in an X-ray production zone extending from 1.2 to 100 R_*. The choice of such a large zone may at first seem surprising, since the prediction of the 1D hydrodynamical modeling places the X-ray emitting regions within a few stellar radii of the stellar surface. However, a detailed global analysis of the high-resolution spectrum of ζ Puppis (A. Hervé et al., submitted) shows that the X-ray emission zone must actually extend up to ∼85 R_* to reproduce the observed spectrum, hence our choice of a large emission zone. Note that the lower boundary of the region was chosen taking into account the analyses of individual fir triplets, which place onset radii in that range (e.g., Waldron & Cassinelli 2007; note that these analyses also consider large emitting zones, with emission formally extending to infinity). All emitting parcels of gas contain the same amount of matter. The line emission is powered by collisional excitation and therefore scales with the density squared. The density of the wind is derived based on the stellar mass-loss rate: from the mass conservation, it is $\rho (r) = \dot{M} / [4 \pi r^2 v(r)]$ where $\dot{M}$ is the mass-loss rate and the wind velocity v(r) is v_∞(1 − R₀/r)^β (with β = 1 and R₀ chosen so that the photospheric velocity v(R_*) = 0.01v_∞), as commonly used for a massive star wind. During motion, each hot fragment expands according to the continuity equation. Hence, the intrinsic unattenuated X-ray luminosity of each hot fragment scales as 1/r²v(r). The probability of finding a hot fragment in the radius interval [r, r+dr] scales with 1/v(r), i.e., emitters are concentrated at inner radii, where the wind is slow. The random radial location of fragments is determined by von Neumann's rejection method (e.g., Press et al. 1992) and their angular distribution is also random, with a uniform distribution over the sphere. Absorption and emission are decoupled: there is no self-absorption for the emitting material and no re-emission of X-rays after absorption. The model allows for further sophistications, i.e., emissivity can have a different scaling with density and density can also be a parameter. However, these parameters affect the shape of X-ray emission line profiles, not relative variability, so in the present study, we use the simplest emissivity.

Once produced, the X-ray emission propagates through the cooler stellar wind which can absorb it. The velocity of this cool wind is assumed to follow the same velocity relation as for the hot wind component (the so-called beta law, see above). The cool stellar wind is assumed to be either a smooth cool wind or a set of cool clumps (randomly distributed over the 1.5–316 R_* range), both cases having the same overall mass-loss rate and optical depth. For the fragmented wind model, we use the "cones" model of Oskinova et al. (2004), with a lateral extent of 1 degree for the spherical absorbers. A sketch of the model geometry is shown in Figure 13. Random radii are generated for each cone using the 1/v(r) probability and von Neumann's rejection method, again. The total mass of a homogeneous wind enclosed between two subsequent radii is considered to be swept up in a dense fragment with the same optical depth as the homogeneous wind material, so that the fragment location is given by r² = [∫^b_adr'/v(r')]/[∫_a^bdr'/r'²v(r')] where a and b are two subsequent random radii in the set determined by von Neumann's method. Note that, depending on the set of radii, cool clumps do not have all the same density, some being optically thick while others are not blocking much light.

Zoom In Zoom Out Reset image size

The wind parameters (v_∞, $\dot{M}$, abundances) are derived from a non-LTE atmosphere model specific to ζ Puppis (Oskinova et al. 2006), matching the optical/UV spectrum of the star. All these fix the clump location, size, density and optical depth, leaving as the only free parameter the number of hot and cool clumps. It may be noted that this model reproduces well the observed X-ray line profiles, despite its simplifications (Oskinova et al. 2004).

Our analysis of observational data did not reveal strong spectral trends in the variability. Therefore, for simplicity and clarity we simulate here only the monochromatic X-ray flux at a few selected representative wavelengths: 6 Å (∼2 keV, midpoint of the total EPIC band), 14 Å (∼0.9 keV, midpoint of the medium EPIC band), and 19 Å (∼0.65 keV, midpoint of the total RGS band).

Some comments should be made about the synthetic light curves. First, the calculated flux is monochromatic and in arbitrary units. Only relative dispersions, for example, can be calculated and compared to the data. Second, these variability light curves are not, strictly speaking, functions of time. The time-dependent radiative hydrodynamic simulations of Feldmeier et al. (1997b) show that on a timescale longer than the flow time (5.8 ks for ζ Puppis, as v_∞ = 2250 km s⁻¹ and R_* = 18.6 R_☉, see Oskinova et al. 2006 and references therein), the wind structure is renewed and is independent of the previous wind configuration. Our model reproduces this situation: each realization of our 2D stochastic wind model is independent of the previous one and each model run represents the wind configuration on a timescale shorter than the cooling time. In other words, we model the wind at some arbitrary moment of time, and compute the X-ray emergent flux at this moment; the next data point is calculated for a randomly different wind configuration for both absorbing and emitting parcels. We do not follow the wind expansion (this will be the subject of a paper by L. M. Oskinova et al., in preparation). While this approach does not allow us to model the detailed time evolution of the X-ray flux on timescales of hours, it allows us to model the relative amplitude of the X-ray variability for long stretches of randomly distributed observations, as appropriate for this XMM-Newton observing campaign. Examples of such "light curves" are shown in the right panel of Figure 13, while relative dispersions are presented in Table 2.

Several conclusions can be drawn from Table 2.

• 1.  
  As the number of emitting or absorbing parcels increases, the variability decreases, as could be expected. However, it is important to note that dispersions of 1% are only found when these numbers approach 10⁵.
• 2.  
  A smooth wind is less variable than an otherwise equivalent clumped wind.
• 3.  
  As wavelength increases (or energy decreases), the relative dispersion for a clumped wind remains stable or slightly increases whereas it decreases for a smooth wind (as already reported in Oskinova et al. 2001b). In a clumped wind, the variability is thus less energy dependent than in a smooth wind—in the extreme limiting case of a wind consisting of only opaque X-ray clumps, no energy dependence is expected (Oskinova et al. 2004, 2006).
• 4.  
  If absorbing clumps are distributed over a smaller region (R_max of 100 rather than 316 R_*), then the variability decreases. This is the effect of a smaller radial separation between clumps.

Comparing these theoretical predictions with the observational results (see previous subsection), we found that the number of emitting and absorbing parcels is huge. Indeed, even in the most favorable case of smooth cool wind absorption, which is the least variable case, more than 100,000 hot X-ray emitting zones must be present and contribute to the X-ray emission so that the relative flux variations remain below 1%. This number further increases when the cool wind fragmentation is also included in the model.

It must be underlined that no previous study has put direct constraints on the number of clumps in O-stars. Some studies exist, however, for evolved massive stars. For example, Lépine & Moffat (1999) report that "between 10³ and 10⁴ clumps in the line emission region are needed to account for the line profile variability of the WR stars" that they analyzed. The X-ray data of ζ Puppis suggest an even larger number, though this should apply to a larger zone than formation regions of optical lines. In addition, Davies et al. (2007) explained the polarization level of luminous blue variables by either a few massive, optically thick clumps (ejection rate of ≲ 0.1 clump per flow time) or many small, optically-thin clumps (ejection rate of ≳ 10³ clumps per flow time), with the latter option usually favored in the literature. Considering that the X-ray emission region and the cool wind absorption region cover several tens of stellar radii, our conclusion appears consistent with Davies' result, despite the different nature of the objects under consideration.

Thus, the sensitive XMM observations reveal that the stellar winds of O-stars are highly structured on small scales. It remains to be seen whether the theory of stellar wind instability can explain such a high degree of fragmentation. Some 2D models of the line-driven instability (LDI) in isothermal⁹ stellar winds "show that radially compressed shells that develop initially from the LDI are systematically broken up by Rayleigh–Taylor or thin-shell instabilities as these structures are accelerated outward" (Dessart & Owocki 2003), hence producing a lot of small-scale 2D structures, but the same authors later found "lateral coherence of wind structures" (Dessart & Owocki 2005) and the question of the size and number of clumps therefore remains unsettled. The hydrodynamical simulations of Feldmeier et al. (1997a; to this day the sole ones that tackled the problem of X-ray generation, though in 1D) predicted that only a few strong shocks, simultaneously present in massive star winds, are responsible for most of the X-ray emission, whereas our results strongly challenge this. The lateral break-up of the X-ray emitting shells, advocated by Feldmeier et al. (1997a) for lowering the X-ray variability, cannot provide an explanation (see Table 2). Moreover, if clumping is induced by sub-surface convection, hydrodynamical stellar evolution codes indicate that the total number of clumps in O-stars should typically amount to 6 × 10³–6 × 10⁴ (Cantiello et al. 2009). Our data indicate values an order of magnitude larger, therefore prompting further investigation on the nature of radiatively-driven stellar winds.

This data set only comprises data from the RGS. There were narrow background flares scattered all over the exposure and a broad flare affected the mid-exposure. At a significance level of 1%, no trend or deviation from a constant is detected. Count rates may be higher toward the end of the exposure, but only at the 10%–20% level. A feeble oscillation can be detected by eye in the Berghöfer band, with a recurrence time of about 10 ks, but it does not reach an amplitude larger than that of the 1σ error bars. The count rate in the medium band and with the longest time bin appears as the most variable (i.e., they show a χ² associated with the lowest significance level, but this significance level is not <1%).

A narrow background flare affected the pn data at mid-exposure. This data set shows a trend toward increasing count rates. Though less obvious in MOS2 and RGS1, this trend is clearly detected through the significant improvement of the χ² when using a linear rather than a constant fit. This trend is detected for the total, soft, and medium bands (EPIC), and for the medium band (RGS), i.e., it concerns photons with energies below 1 keV. For the EPIC-MOS1 light curve in the total band, the increase rate is 1.7 ± 0.7 × 10⁻⁶ counts s⁻² (i.e., a 3% increase over the ∼40 ks exposure). For EPIC, count rates in the total band and/or calculated using the longest time bin appear as the most variable; for RGS, the medium band data and/or the light curves with the longest time bins are the most variable.

No soft proton flare is detected for this exposure. Data are compatible with a constant count rate, without any trend detected or any obvious oscillation. The pn observation is cut into two parts for this revolution and the next two (see Paper I): here, the second pn data set (taken with Medium filter) appears more variable than the first part (taken with Thick filter), in all bands and for the smallest time bin, but this difference is not formally significant. For EPIC, the hard band light curves (especially those obtained with the smallest time bin) are the most variable; on the contrary, for the other bands, the light curves calculated using the longest time bins are the most variable. RGS data do not yield any coherent result as to which time bin and energy band is the most variable.

A large flare occurred near the end of the observation. Data are compatible with a constant, without any trend detected or any obvious oscillation. For EPIC, the hard band light curves (especially those obtained with the longest time bin) as well as light curves taken with the extreme time bins (200 s and 5 ks) are the most variable; for RGS, the shortest time bins generally yield the most variable light curves.

A large flare occurred near the end of the observation, it did not affect the MOS data. Data are compatible with a constant count rate, without any trend detected or any obvious oscillation. The first pn data set appears more variable than the second part, in all bands for the 1 ks bin, but this is not highly significant. Berghöfer's band appears as the most variable for both EPIC and RGS; the longest time bins yield the most variable RGS light curves.

Only pn data are available for EPIC. A large flare occurred during the second half of the observation. The data from pn and RGS1 indicate a slight increase (linear trend better than a constant at the <5% level) of the count rates in the second half of the observation, but this is not corroborated by the RGS2. While there is no time bin or energy band favoring the variability properties of the RGS data, the medium band data and light curves calculated using the longest time bins (except for the hard band) yield the most variable light curves in the EPIC-pn.

A large flare occurred at the end of the MOS data sets, or in the middle of pn and RGS data sets. Data are compatible with a constant, without any trend detected or any obvious oscillation. There is no time bin or energy band favoring variability.

A large flare occurred during the last third of the observation. Only pn data are available for EPIC. Data are compatible with a constant count rate, without any trend detected or any obvious oscillation. For RGS, the medium band light curve calculated with the 2 ks bin appears as the most variable; for EPIC, the hard band data are the most variable except for the longest time bins.

A large flare occurred during the second half of the observation. The count rate recorded in the hard band is not compatible with a constant—it shows a shallow decreasing trend—but only at significance levels of 1%–10% (i.e., not formally significant). For EPIC, the hard band light curve appears as the most variable, but conclusions are unclear for the time bins: MOS data show more variability in the smallest time bins, while pn data favor variability in the longest ones. For RGS, the soft band appears as the most variable.

A large flare occurred during the last third of the RGS observations. While the EPIC data (especially pn) are much improved when fitted by an increasing trend, the RGS data are more compatible with a decreasing trend. Both types of instruments agree, however, that the longest time bins generally yield the most variable light curves.

A large flare occurred near the end of the observation. Data are compatible with a constant, without any trend detected or any obvious oscillation. For EPIC, the largest variability (though not significant) is seen for the total band data while the medium data appear as the most variable in RGS.

Only pn data are available for EPIC. There is no localized flare for this exposure but the background smoothly increases in the second half of the observation. An increasing trend is detected for the pn, in the total and soft bands; it is also detected in RGS (especially RGS1) in the total band, but with a lower significance level. For RGS, the medium band light curve appears as the most variable for the longest time bins; for EPIC, the largest variations are detected for shortest time bins in the hard band and the longest time bins in the total band.

A few narrow background flares are scattered over the exposure, and a larger one occurs at the end of the observation. Data are compatible with a constant rate, without any trend detected or any obvious oscillation. For RGS, there is no time bin nor energy band favoring variability; for EPIC, the most variable data are found for the longest time bin in the hard band.

A few narrow background flares are spread over the exposure, and a larger one occurs at the beginning of the observation. Only MOS data are available for EPIC. A quadratic fit strongly improves the χ² in the total and medium bands of the MOS (which are not compatible with a constant count rate, at the 1%–10% level), and in the total band of the RGS. For both types of instruments, the longest time bins provide the most variable light curves, but the most variable energy bands are the total and soft bands for MOS and RGS, respectively.

A large flare occurred during the first third of the observations. Data in the medium, total, and Berghöfer's EPIC bands are incompatible with a constant and significantly better fitted by a linear decrease. The non-constancy in the Berghöfer's band is also detected in the RGS data, though with a lower significance level. For EPIC, the medium band data and light curves calculated using 5 ks bins appear as the most variable; for RGS, Berghöfer's band appears the most variable, especially at small time bins.

A large flare occurred near the end of the observation. Data in all EPIC bands but the hard one are incompatible with a constant and significantly better fitted by a linear increase. The trend in the total band is also detected in the RGS data. For EPIC, the total and medium band data appear as the most variable, especially for the longest time bins; for RGS, the longest time bins yield the most variable light curves for the total and Berghöfer's bands.

A small flare occurred at the beginning of the RGS observation. Data in all EPIC bands but the hard one are incompatible with a constant and are significantly better fitted by a large linear increase or, even better, a quadratic increase. The non-constancy and the improvement by an increasing trend are also detected in the RGS data for the total and medium bands. For EPIC, the hard data appear as the least variable, while the light curves calculated with the longest time bins are the most variable. For RGS, the medium band data are the most variable.

A large soft proton flare occurred near the end of the RGS observation. Data in the total and medium EPIC bands are significantly better fitted by a linear, or even quadratic, increase, while the hard and Berghöfer's bands are only better fitted by a quadratic increase. The increasing trend is also detected in the RGS data for the total band. An oscillation is visible in Berghöfer's band on both RGS and EPIC, with a recurrence time of about 50 ks. In general, the longest time bins yield the most variable light curves.