© ESO, 2018

1. Introduction

Massive stars are well known for their ability to launch dense and fast stellar winds with important consequences for their evolution and their feedback on the environment. However, numerous uncertainties remain on their properties, notably regarding their structure. The first question that arose regarded clumping. The winds of massive stars are very different from those from solar-like stars since the winds of massive stars are propelled by intense UV radiation (Lucy & Solomon 1970). Photons can be absorbed by the many transitions from metallic ions; since photons initially flow radially while scattering occurs in all directions, this leads to a net radial acceleration of those ions and, subsequently, of the entire atmospheric material after ion interactions. However, this line driving is an unstable process. If we take Doppler shifting into account, a slightly faster ion can absorb more photons and accelerate further, while slightly slower ions are deprived of photons and decelerate further since the needed photons were already absorbed by the layers just below the ion position. The line-driven winds are thus intrinsically unstable (see e.g. Lucy & Solomon 1970; Carlberg 1980) with numerous small-scale stochastic perturbations pervading the wind flow (e.g. Eversberg et al. 1998). The presence of such perturbations alters the interpretation of mass-loss diagnostics, and this impact varies with the chosen diagnostic since some of the diagnostics depend on the square of the density while others vary linearly with density. After considering those effects, actual mass-loss rates are now believed to be lower than previous used unclumped values by a factor ~ 3 (e.g. Bouret et al. 2005). These wind instabilities are also responsible for the generation of X-rays (Lucy 1982; Owocki et al. 1988; Feldmeier et al. 1997; Oskinova et al. 2004, and references therein). Stochastic variability should then arise on hourly timescales (e.g. wind flow time R_∗/v_∞ is ~1.6 h for ζ Pup) with an amplitude linked to the number of involved clumps. However, dedicated analyses of sensitive X-ray observations failed to detect such changes, implying a very large number of clumps (>10⁵ for ζ Pup; Nazé et al. 2013).

The structure of stellar winds is not restricted to small-scale clumping, however, because large-scale features are also inferred to be ubiquitous in massive hot-star winds. Indeed, variable discrete absorption components (DACs) observed in unsaturated UV P Cygni profiles of O stars (Massa et al. 1995; Howarth et al. 1995) have been interpreted as a manifestation of co-rotating interaction regions (CIRs; see e.g. Cranmer & Owocki 1996; Lobel & Blomme 2008). Related variability was detected for line profiles in the optical range as well (e.g. Prinja et al. 1996; Kaper et al. 1997). These large-scale (>R_∗) wind disturbances are believed to originate at the base of the wind from bright spots on the photosphere. The presence of such spots has also been inferred from optical variations, for example modelling of the photometry of ξ Per and ζ Pup (Ramiaramanantsoa et al. 2014, 2018) and modelling of the He ii λ4686 profile of λ Cep (Sudnik & Henrichs 2016). But, until recently, CIRs were not considered to be a potential source of X-rays since models predicted only a velocity kink and not a large jump (e.g. Lobel & Blomme 2008).

The situation changed with the detection of trends in the X-ray light curves of several massive stars. The analysis of an extensive XMM-Newton dataset of ζ Pup, while failing to unveil the expected short-term stochastic variability associated with clumping, clearly revealed shallow monotonic increases, monotonic decreases, or modulations of the X-ray light curves with a relative peak-to-peak amplitude of ~20% (Nazé et al. 2013). No full cycle could be observed over the exposure lengths, which are usually less than a day, but time analysis tools point towards potential timescales between a day and a few days. A few other (relatively bright) single massive stars hinted at a similar variability. Variations with 15% amplitude were detected in ξ Per (Massa et al. 2014). The Chandra light curve could be fitted using a sinusoidal function with the period detected in UV P Cygni profiles, although less than one cycle was observed and a linear trend was able to fit the data as well as this sinusoidal modulation. Changes by 15% were also detected in the four XMM-Newton observations of λ Cep (Rauw et al. 2015) and they appear to be compatible with the 4.1 d variation of the star observed in Hα at the same epoch, although each exposure covers only about 1 d and no single cycle was fully covered. In the past, a variation of 20% amplitude was also reported for ζ Oph using the less sensitive ASCA observatory, and it appeared possibly related to the rotation period detected in the UV (Oskinova et al. 2001). As these X-ray changes presented strong similarities (in pattern or timescale) to the optical/UV variability associated with CIRs, these features were naturally suggested to explain the X-ray variations. A definite link is not proven yet, however, as there is a clear misfit between observing cadence at X-ray energies and the possible periodicities: no massive star has ever been monitored over a full – or to be more conclusive, several consecutive – variability cycle(s) in X-rays.

To better characterize these X-ray trends, additional data are thus needed. This paper aims at such a goal, thanks to new observations of ζ Pup and revisiting old data. The choice of the target comes from its brightness: the signal to noise achieved for the other massive stars with such variability is many times below that which can be achieved for ζ Pup. Re-examining this star requires a good knowledge of its cycle length, however. In the past, several variability timescales were reported for this star. In particular, a value of ~5.1 d was proposed as the rotation period for this star by Moffat & Michaud (1981) from Hα observations, but this period was not recovered in the last 20 yr. The presence of a 15–19 h period was found in UV (Prinja et al. 1992; Howarth et al. 1995), Hα line (Reid & Howarth 1996) and even X-ray observations (ROSAT; Berghoefer et al. 1996) although no sign of this period was later found in ASCA or XMM-Newton data (Oskinova et al. 2001; Nazé et al. 2013). Recent extensive photometry now shows that timescale of 15–19 h remains coherent for a few days at best, with very frequent complete disappearances (Ramiaramanantsoa et al. 2018). Recently, Howarth & Stevens (2014) unveiled a much more stable periodicity in ζ Pup of P = 1.780938 ± 0.000093 d, from long, 40 to 236 d, Coriolis/SMEI (Solar-Mass Ejection Imager) observing runs spread over three years. The presence of this signal is confirmed in BRITE (BRIght Target Explorer) observations of ζ Pup taken in 2014–2015 as well as in observations of the He ii λ4686 line profile taken simultaneously (Ramiaramanantsoa et al. 2018). While Howarth & Stevens (2014) proposed this signal to be pulsational in nature, Ramiaramanantsoa et al. (2018) interpreted this behaviour as a rotational modulation associated with bright surface spots. Both cases could however launch CIRs in the wind, which should thus present the same recurrence timescale as detected for the variations of the He ii λ4686 wind emission line. Armed with this precise period, this paper revisits the X-ray variability of ζ Pup. Section 2 presents the X-ray data and their reduction, while Sect. 3 presents and discusses the observational results, and Sect. 4 summarizes these results.

2. Observations and data reduction

2.1. Swift

At our request, ζ Pup (O4I) was observed by Swift in December 2016 and January 2017. Each time, the X-ray emission of ζ Pup was followed with Swift-XRT (X-Ray Telescope) over one entire 1.78 d cycle with a set of 1 ks observations taken every 3–6 h. Three older datasets and a set of calibration exposures taken simultaneously with the last XMM-Newton observation were also available and are thus also analysed here; see Table 1 for the full list of observations.

Corrected (full PSF) count rates in the same energy bands as previous XMM-Newton analyses (Nazé et al. 2012, see also next subsection) were obtained, considering only grade 0, for each observation from the UK online tool¹. Heliocentric corrections were added to mid-observation times. To get spectra, XRT data were also locally reprocessed using xrtpipeline of HEASOFT v6.18 with calibrations available in winter 2016 and an explicit input of the Simbad coordinates of ζ Pup.

ζ Pup is a bright target both in the optical (V = 2.25) and at X-ray energies, hence we could not use the photon-counting mode of Swift-XRT, which would be heavily affected by optical loading for such a bright source. However, even the alternative (windowed timing) mode has some problems. Indeed, below 0.5 keV, the count rates sometimes increased erratically (up to a factor of 3) and the spectra showed a varying, unrealistic bump. Therefore, to avoid contamination by optical/UV photons as much as possible, only the best quality data (grade = 0, E> 0.5 keV) were considered. Furthermore, as recommended by the XRT team, the source spectra were extracted within a circular aperture² centred on the Simbad coordinates of the target³ with radius 10 px to minimize the background contribution, while an annular background region with radii 70 and 130 px was used⁴. Spectra were then binned (using grppha) to reach 10 counts per bin (or a signal-to-noise ratio of 3 per bin). The redistribution matrix file (RMF), to be used for calibrating the energy axis of spectra at the observing date, was retrieved from the calibration database, after checking that the general RMF matrix or the position-dependent RMF matrices⁵ yielded the same results in our case. Finally, specific ancillary response files (ARFs), used to calibrate the flux axes of spectra, were calculated for each dataset (using xrtmkarf) considering the associated exposure map.

2.2. XMM-Newton

ζ Pup was used as a calibration target for XMM-Newton, resulting in an exceptional number of 26 exposures taken since 2000. These observations were taken in a variety of modes and with various exposure lengths (see Nazé et al. 2012, for details). For this new analysis, these datasets were re-reduced using the most recent calibration available (SAS v15.0.0 with calibration files available in winter 2016) and following the recommendations of the XMM-Newton team⁶.

Fig. 1 Count rates in several energy bands of ζ Pup, and associated hardness ratios for XMM-Newton data from 2003–2006. The phases were calculated using P = 1.780938 d and T0 = 2 450 000 (Howarth & Stevens 2014). To facilitate the comparisons, the ratio between upper and lower limits of the y-axes are all equal.

This paper primarily deals with variability issues, hence it is focussed on light curves. The most sensitive data in this regard are those taken by the pn camera, but they are affected by pile-up and/or optical loading in all but one mode (small window + thick filter). This mode has been used since 2003, but the first datasets were rather short and very affected by soft proton flares, mostly due to solar activity; the pn camera was even not switched on in one case (Rev. 1164 in 2006). The pn data have the additional advantage that they are the least affected by sensivity degradation; this detector is in fact considered stable by the XMM-Newton calibration team, to within a few percents across the mission lifetime, which facilitates inter-pointing comparisons.

After the initial pipeline (epproc) processing, the pn data were filtered to keep only those with the best quality, using pattern = 0 + flag = 0 and eliminating the flares whenever needed: time intervals with count rates above 10 keV beyond 0.045 cts s^-1 were discarded. The barycentric correction was calculated and applied to the recorded times (SAS task barycen). Light curves were then extracted in the same regions as in Nazé et al. (2012): the source region is a circle centred on the Simbad coordinates of ζ Pup and with a radius of 850 px (42.5′′); the background region is also a circle centred on a nearby position (an identification number given in Table 1 refers to the background region positions listed in Table 2 of Nazé et al. 2012); these circular regions have a radius of 700 px except for Revs. 2363 and 2817, where it had to be reduced to 600 px and 450 px, respectively. We considered time bins of 5 ks and four energy bands, i.e. soft band S = 0.3–0.6 keV, medium band M = 0.6–1.2 keV, hard band H = 1.2–4.0 keV, and total band T = 0.3–4.0 keV (see Nazé et al. 2012). The light curves were further processed by the task epiclccorr, which corrects for loss of photons due to vignetting, off-axis angle, or other problems such as bad pixels. In addition, to avoid very large errors and bad estimates of the count rates, we discarded bins displaying effective exposure time lower than 50% of the time bin length. For these light curves, the background contribution is less than half the 1σ error bars, hence background plays no significant role in the variability studied here.

Fig. 2 Same as Fig. 1 but for XMM-Newton data from 2007–2012.

Fig. 3 Same as Fig. 1 but for XMM-Newton data from 2013–2017.

Fig. 4 Same as Fig. 1 but for Swift data; those from October 2016, December 2016, January 2017, and April 2017 are shown as red triangles, black dots, blue crosses, and cyan squares, respectively.

Fig. 5 Hardness ratios as a function of the source brightness for all XMM-Newton data and for all Swift data.

As explained in the next section, XMM-Newton spectra of two extreme cases were specifically analysed. In those cases, pn spectra were extracted for the chosen time intervals via the task especget in the same regions as the light curves. Dedicated ARF and RMF response matrices were also derived for each observation. The pn spectra were finally grouped, with specgroup, to obtain an oversampling factor of five and to ensure that a minimum signal-to-noise ratio of 3 (i.e. a minimum of 10 counts) was reached in each spectral bin of the background-corrected spectra. After the initial pipeline (rgsproc) processing, Reflection Grating Spectrometer (RGS) data were also filtered to keep only events in the chosen time intervals and the associated spectra and their response matrices were then calculated via the appropriate SAS tasks (rgsspectrum, rgsrmfgen). This was carried out for the default regions since ζ Pup has no neighbour of similar X-ray brightness. Finally, the RGS spectra were individually grouped with specgroup to obtain an oversampling factor of five; calibrated combinations of both RGS and both orders were also derived via rgsfluxer.

3. Results

3.1. Variations in light curve amplitudes

Figures 1–4 show the observed light curves, in several energy bands, phased with the ephemerides of Howarth & Stevens (2014, P = 1.780938 d, T₀ = 2 so that several light curves can be represented in the same graph using the same x-axis, facilitating the (direct) comparisons. The amplitude of the variations is immediately seen to change from one observation to the next. To quantify it, a variability index VI was defined as (max−min)/(max + min), i.e. the ratio between half the peak-to-peak amplitude and the mean. If the variations were sinusoidal in shape, this would yield the fractional semi-amplitude of the sine-wave. The values of this index are provided in Table 1. In the XMM-Newton datasets, the exposure lengths are similar, but VI varies between 1 and 6% in the total band, and the most recent datasets are the least variable. Indeed, the count rate in the total band varies from slightly below 6 cts ^-1 to slightly more than 7 cts ^-1 in 2007–2011 (Fig. 2), but only between 6.2 and 6.7 cts ^-1 in 2003–2006 and 2012–2016 (Fig. 3). Moreover, some light curves clearly appear above or below others⁷. Analysing all XMM-Newton datasets together therefore yields larger overall VI: 0.083 ± 0.005, 0.079 ± 0.009, 0.099 ± 0.007, and 0.095 ± 0.013 for the T, S, M, and H energy bands, respectively, confirming the value of 20% peak-to-peak variations mentioned in Nazé et al. (2013). This 20% variation amplitude is also much higher than the possibly remaining calibration uncertainties (XMM-Newton help desk, private communication). The VI for various energy bands agree, usually well within 3σ, but the maximum variability seems to be reached in the medium band.

To correct these estimates from the variability due to noise, we also determined the fractional variability amplitude F_var of the data (see Appendix A of Edelson et al. 2002). F_var is the square root of the normalized excess variance which is calculated by subtracting the average variance due to noise from the observed variance of the data and then normalizing the result by the mean. This gives an idea of the intrinsic variability of the source. The fractional variability amplitude (Table 1) varies from 0 to 4.2% in the T band (0.3–4.0 keV), reaching the highest values (>1%) for data in 2007–2013. There is negligible variability for the other years (i.e. F_var ~ 0 within 3σ), confirming the VI results. Considering all data together, we found F_var of 0.034 ± 0.002, indicating that about half of the observed variability is due to measurement noise. The fractional variability amplitude amounts to 0.033 ± 0.002, 0.040 ± 0.002, and 0.027 ± 0.002 for the S, M, and H energy bands, respectively; again an agreement within 3σ is found between F_var(S) and F_var(M) or between F_var(S) and F_var(H), but the larger value of the variability in the medium band appears here more clearly.

Fig. 6 Comparison of the pn (left) and RGS (middle and right) spectra observed at minimum brightness (first part of Rev. 1814; solid black line) and maximum brightness (first part of Rev. 2183; dotted red line). Blue segments above lines of RGS spectra indicate the intervals chosen for the moment calculations.

Fig. 7 Left: same as Fig. 4 but for the fluxes (Table 4). The green stars correspond to data from February 2013. Middle and right: comparison between the XMM-Newton light curve of April 2017 (Rev. 3172, black dots and solid lines) and the (arbitrarily scaled) simultaneous Swift light curve (red crosses and dotted lines) derived from count rates (middle) or fluxes (right).

The hardness ratios indicate only slightly smaller variability indices than the band-by-band values, with VI(M/S) = 0.058 ± 0.011 and VI(H/M) = 0.081 ± 0.015 or F_var(M/S) = 0.019 ± 0.002 and F_var(H/M) = 0.003 ± 0.002. Therefore, we examined correlations between energy bands. The correlation between soft and hard count rates is poor (r = 0.24), but those between soft and medium count rates and between medium and hard count rates are better (r of 0.79 and 0.49, respectively). Furthermore, there seems to be a correlation between hardness and brightness (Fig. 5): as the source brightens, M/S increases (correlation coefficient r of 0.35) while H/M decreases (r = −0.59), i.e. the flux of ζ Pup principally changes in the medium band.

The Swift light curves display variations of larger amplitudes, although the variability indicators (VI and F_var) can be reconciled with those measured for XMM-Newton within 3σ considering their larger errors. While the instrumental sensitivities of both observatories are different, which explains part of these differences, Swift light curves may also remain somewhat affected by optical loading. It is however encouraging to find the same brightness–hardness correlation in these data (except for two outliers, see rightmost panel of Fig. 5), demonstrating that the X-ray signal dominates over contamination. We come back to Swift data in the next section.

3.2. Spectral differences

In the XMM-Newton light curves, the minimum count rate is reached for the beginning of Rev. 1814⁸, while the maximum is reached during the first part of Rev. 2183⁹. The difference in total band count rate between these two cases amounts to at least 0.6 cts s^-1, which corresponds to >12 times the error bars on the 5 ks bins. Thus, the beginning of Rev. 1814 and Rev. 2183 are well differentiated in terms of brightness. As these two cases represent the most extreme brightness variations of ζ Pup, their pn and RGS spectra were examined in more detail to determine the nature of the changes. Figure 6 compares the spectra taken at minimum and maximum brightness. While a change in intensity is readily detected in both low- and high-resolution spectra, no other variation, such as shift of spectrum towards lower or higher energies or variation in line profile shapes, appears obvious by eye.

Fig. 8 Subgroup of XMM-Newton light curves can be phased with a period of 1.78065 d to form a relatively coherent behaviour (top), but considering all available light curves only lead to scatter plots, whatever the value of the 1.78 d period considered (bottom).

Fits were then performed for the pn spectra with Xspec 12.9.0i, using Asplund et al. (2009) for the reference solar abundances. A model wabs ∗ (vphabs ∗ vapec + vphabs ∗ vapec + vphabs ∗ vapec) was used, where the first absorption was fixed to the interstellar value of 8.9 × 10¹⁹ cm^-2 (Diplas & Savage 1994) and CNO abundances of the other components were fixed to those derived in Hervé et al. (2013)¹⁰. This model is similar to that used in Nazé et al. (2012), although that model had a fourth thermal component. In our case, adding this hottest component did not significantly improve the fit, so it was discarded for simplicity. Fitting results are provided in Table 2. In a first trial, the absorptions, temperatures, and normalizations of all thermal components were free to vary. Higher values of these components were derived for the maximum brightness case. However, a second trial was then made, fixing the temperatures to their average values. The absorptions were then found to be consistent within errors, letting only normalization factors significantly vary between minimum and maximum brightness cases. In particular, the smallest variation is detected for the hottest component (0.644 keV) and the largest variation is detected for the intermediate temperature (0.273 keV). The largest flux change is recorded in the medium band, while the flux ratios in the soft and hard bands are similar, which confirms the light curve results. The variations in ζ Pup thus seem to result from a global increase in the strength of the emission, particularly that of the warm plasma (kT ~ 0.3 keV).

To get a more detailed look at variations, X-ray lines recorded by RGS can be studied. However, the X-ray lines of ζ Pup are well known to be non-Gaussian (e.g. Kahn et al. 2001). Therefore, to best compare the RGS spectra, line moments (orders 0 to 3) were calculated in the velocity space for strong and relatively isolated lines (a method used e.g. by Cohen et al. 2006): Ne x λ12.132, Fe xvii λ15.014, 15.261, 16.780, 17.051/96, O viii λ18.967, and N vii λ24.779. Table 3 provides the results, showing that there are no significant differences between position, width, and skewness of lines at minimum or maximum brightness. Comparing the line profile shapes rather than their overall properties does not reveal any obvious feature either, considering the errors (Fig. 6). The changes thus seem limited to the line intensities, as found before for low-resolution pn spectra.

As a last investigation, the Swift spectra were also analysed. The chosen model is the same as used for XMM-Newton, but with some parameters fixed. The XMM-Newton fits showed that, if temperatures are fixed, absorptions are similar whatever the brightness, so both parameters were fixed here. Furthermore, since the Swift spectra are only examined above 0.5 keV and the first thermal component dominates at energies below 0.5 keV, the norm(1) /norm(2) ratio was also fixed. Following recommendations of the XRT team, the spectral fitting was performed considering the possibility of energy scale offsets (command “gain fit” under Xspec, with slope fixed to 1). Results are provided in Table 4 and the left panel of Fig. 7. The fluxes present lower variability indicators than the light curves, leading to a better agreement with XMM-Newton data. A direct comparison of the simultaneous XMM-Newton/Swift observations of April 2017 is shown in the middle and right panels of Fig. 7, and it further confirms that fluxes derived from Swift spectra better correlate with the XMM-Newton data than the Swift count rates. Problems due to optical loading seem negligible; the Swift spectra truly reflect the actual brightness variations of ζ Pup and can thus be trusted. But we need to be aware of their larger Poisson noise, due to shorter exposure times and lower effective area of Swift, and of the potentially larger systematic uncertainties (up to ~ 20% in the worst cases).

3.3. Shapes of the variations

There is a second striking feature in Figs. 1 to 4: the shape of the light curves. Indeed, one can observe near constancy (e.g. Rev. 2540), linear trends (e.g. Rev. 1343), or parabolic trends (e.g. Revs. 1983 and 2183), and there is even one case in which the X-ray luminosity experiences a slow increase followed by a steeper increase (Rev. 1814). Not only are the mean levels and amplitudes different, well beyond the calibration uncertainties, but the light curve shapes are sometimes incompatible with each other. Indeed, the steep increase at the end of Rev. 1814 cannot be reconciled with the shallower increases observed in, for example Revs. 1620, 1983, or 2183 (Fig. 2).

A last interesting feature is the phasing of these X-ray light curves, performed using the very precise value for the periodic variability in optical light of ζ Pup derived by Howarth & Stevens (2014) from SMEI data, i.e. P = 1.780938 ± 0.000093 d. Assuming the period to remain stable and considering its error bar, a 1σ error on the period would yield a phase shift of only 0.02 after one year, or 0.1 after five years. There should thus be little phase shift between individual X-ray light curves in each of the Figs. 1 to 4, although phase shifts as large as ~ 0.3 (for a 1σ error) can be expected between the beginning and the end of the 14 yr observing campaign in X-rays. However, the observed light curves do not form a coherent behaviour (Figs. 1 to 4 and 8, bottom). Even data taken the same year – an interval for which phasing errors and any instrumental sensitivity changes should be negligible – cannot be reconciled with each other. This is the case of Revs. 0795 and 0903, both taken in 2004, of Revs. 0980, 1071, and 1096, taken in 2005, (see Fig. 1) and, to a lesser extent (only some bins do not match within 3σ for these, see Fig. 3), of Revs. 2533 and 2540, both taken in 2013, and of Revs. 2817 and 2911, both taken in 2015. The light curves of a subgroup of observations, however, can be combined to form a coherent, sinusoid-like overall shape. This is best achieved by slightly modifying the period to, for example 1.78065 d (Fig. 8, top), a value similar to that recently found by Ramiaramanantsoa et al. (2018) – these authors have found P = 1.78063 ± 0.00025 d for the weighted average of the period values derived from the analysis of the individual datasets (SMEI 2003−2004, SMEI 2004–2005, SMEI 2005−2006, and BRITE 2014−2015). This provides marginal support for the presence of this period in the X-ray data, but only at some epochs; it should be noted that some exposures taken in between those of this subgroup (e.g. Rev. 1814) do not show the same variations (Fig. 2). The X-ray variations of ζ Pup definitely appear to differ widely from a perfectly stable modulation, both in terms of phasing and overall light curve shape. Retrospectively, this irregular behaviour may explain the variability detection obtained with ROSAT (Berghoefer et al. 1996) and the subsequent non-detection in ASCA data (Oskinova et al. 2001). The amplitude of variation found by Berghoefer is similar to what we detect with XMM-Newton; simply, as in 2007–2011, such strong variability was transient.

Such a behaviour is not totally unprecedented for this star. Indeed, Howarth & Stevens (2014) mentioned that the amplitude of the photometric cycle varied between 5 and 12 mmag. Ramiaramanantsoa et al. (2018) further showed that the light curve shape changes as well. This could perhaps be understood if the photometric changes arise in a non-constant phenomenon, for example stellar spots. The behaviour remains coherent but only for a while during the lifetime of the spot, and as new spots arise, a phase shift and a change in the light curve shape can be observed. In any case, this somewhat erratic behaviour is usual for Oef stars such as ζ Pup, where cycles with similar but not perfectly identical periodicities are often found in spectroscopic time series (Rauw et al. 2003, 2015; De Becker & Rauw 2004; Sudnik & Henrichs 2016). However, at a given date, the star should display one single behaviour. If a CIR scenario in which structures at the stellar surface launch strong wind pertubations is correct, a tight correlation is then expected between optical and X-ray light curves. In fact, in our case, this correlation between the behaviour of ζ Pup at optical wavelengths and that at high energies can be directly evaluated as some photometric campaigns were obtained at the same time as the X-ray observations of ζ Pup, for example SMEI data in 2004–2005 (Howarth & Stevens 2014) and BRITE data in 2015–2017 (Ramiaramanantsoa et al. 2018, and in prep.).

For BRITE data, the noise on individual photometric points is small, so that two comparisons are possible: (1) simultaneous X-ray and optical data (Fig. 9) and (2) the X-ray light curve and the average 1.78 d optical cycle around the same date (Fig. 10). The latter solution compares the X-ray light curves only to the 1.78 d cycle while the former compares these light curves to the full optical photometry, which includes both the 1.78 d cycle and a stochastic component of similar amplitude (Ramiaramanantsoa et al. 2018). Because the origin of the X-ray variations is uncertain, both comparisons should be considered. For SMEI data, the noise on individual photometric points is rather large, so that only averages could be used, erasing the short-term stochastic component. In fact, the optical data taken within ±15 d of the X-ray observing dates (Table 1) were combined and those averages are then compared with the associated X-ray data in Fig. 11.

Fig. 9 Comparison of the X-ray light curves with the simultaneous BRITE photometry in April 2015 (combined light curve from the two BRITE filters, for details see Ramiaramanantsoa et al. 2018), and December 2016–April 2017 (only one satellite, BRITE-Lem, in the blue filter, see Ramiaramanantsoa et al. 2018, andin prep.). For Swift, fluxes from the spectral fits (Table 4) are used since they are more reliable.

Fig. 10 Same as Fig. 9 but using epochally averaged BRITE light curves valid for the date of the X-ray observations (cycles associated with “Part V” data for April 2015, “Part VI” for Dec. 2016, “Part VII” for Jan. 2017, and “Part IX” for April 2017, for details see Ramiaramanantsoa et al. 2018). The ephemerides are those adopted by Ramiaramanantsoa et al. (2018, and in prep.).

Fig. 11 Comparison of the X-ray light curves with the average 1.78 d photometric cycle recorded with SMEI ± 15 d around the X-ray observing date. The ephemerides are those of Howarth & Stevens (2014): P = 1.780938 d and T0 = 2 450 000.

In 2004, X-ray observations were taken in April and November (Fig. 11). The former data were taken during an optical photometric minimum and the X-ray light curve is relatively stable, and hence there may be a weak correlation between the two; the latter data were taken on the descending branch between the maximum and minimum optical brightnesses, but X-rays show an increase in strength (anti-correlation). In 2005, there were three X-ray observations in April, October, and December (Fig. 11). All of these observations were taken at and after the minimum optical brightness phase, but the first and third X-ray dataset are stable so there is no correlation between X-ray and optical for these datasets; the second X-ray light curve reveals an increase in X-ray brightness so a correlation is detected. In April 2015, the X-ray observations were taken as ζ Pup displayed a marked increase in optical brightness (Figs. 9 and 10) but the X-ray light curve remained stable: there is no correlation between the two. In December 2016, ζ Pup was observed by Swift as it slowly brightened at optical wavelengths but the X-ray light curve displays a dip and then remained relatively stable (Fig. 9): the last five points of the X-ray light curve are at a similar level, while their simultaneous BRITE data show a clear increase. In January 2017, ζ Pup remained stable both in optical and X-rays. Finally, in April 2017, the X-ray light curve shows a marked decrease while the optical photometry shows both an increase and a decrease over the same interval (Figs. 9 and 10). There is thus no obvious correlation between X-ray and optical emissions in the currently available data¹¹.

Thus, in general, there does not seem to be a clear, unique correlation between the behaviour of ζ Pup in optical and X-rays, even when we allow for short, fraction of day, time lags. Furthermore, the only potential evidence for a 1.78 d cycle in the X-ray domain comes from a subset of observations (mostly from 2007–2013, see top of Fig. 8), but no optical monitoring is available these years. It is thus difficult to assess whether the brightness of ζ Pup in the optical showed a similar variability at that time, or if on the contrary the optical changes had completely disappeared, in a scenario in which the variability source had become more energetic. In summary, we cannot directly and systematically relate the X-ray variability and optical changes, and therefore we find no support for the hypothesis of a simple, direct connection. The infrequent presence of a 1.78 d timescale at X-ray energies is possible but requires further checking, using additional monitoring, to be fully demonstrated.

4. Conclusion

Because ζ Pup is a bright and nearby object, this star has been the target of numerous investigations, including at high energies. In particular, it is known to be variable in X-rays since at least Berghoefer et al. (1996). Using X-ray datasets obtained by XMM-Newton before 2010, Nazé et al. (2012, 2013) found a lack of the short-term stochastic changes expected to be associated with the clumpy nature of the massive-star winds, but they detected longer term trends. Similar features were then found in other massive stars and they are currently thought to be associated with large-scale structures in the wind, i.e. wind perturbations somehow triggered at the stellar surface (perhaps from spots or pulsations). As the latter can be detected at optical wavelengths, a correlation is expected between the X-ray and optical behaviours.

Half a decade later, new X-ray observations of ζ Pup were obtained by XMM-Newton and Swift, and a clear, stable periodicity of 1.78 d was identified in long optical photometric runs of the star (Howarth & Stevens 2014; Ramiaramanantsoa et al. 2018). We thus revisited the X-ray variability of ζ Pup.

The full dataset of X-ray light curves now reveals their variety: both the amplitude and the shape of the variations change from one exposure to the next. A subgroup of observations displays very large variations that can be combined to form a sinusoid-like shape for a period of 1.78 d. However, the other X-ray datasets behave in a different way, often with rather constant light curves of different mean levels that are difficult to reconcile with one another or with the subgroup light curves.

Hoping to gather more detailed information on the nature of the X-ray changes, we then compared XMM-Newton spectra (both EPIC at low resolution and RGS at high resolution) taken at extreme X-ray brightnesses since changing velocities or temperatures routinely help constrain the origin of high-energy variations. No significant changes were found for the temperatures, absorptions, and line shapes (centroid, width, and skewness). This notably rejects the scenario in which variability arises from changing absorption as denser wind structures come and go into the line of sight. Only changes in flux, slightly stronger in the medium (0.6–1.2 keV) energy band than at softer or harder energies, are detected. The impossibility to constrain the hot plasma properties more precisely precludes us from clarifying its nature.

The ephemerides of Howarth & Stevens (2014) are precise enough that datasets separated by a few years can be compared in phase. However, apart maybe from a subgroup, the X-ray light curves cannot be phased coherently because of incompatibilities arising even for datasets taken the same year where calibration and ephemeris uncertainties play little role. Since the 1.78 d cycle of ζ Pup varies in amplitude as well as shape, and since additional optical variations (stochastic changes of similar amplitude) are also present, a more meaningful comparison is achieved if X-ray data and optical photometry obtained simultaneously are examined. This comparison can be made thanks to SMEI in 2004–2005 and BRITE in 2015–2017 but again, no clear, single result is obtained. Sometimes a correlation between optical and X-ray data is detected, sometimes there is an anti-correlation, and there is a total absence of correlation at still other times. Therefore, there is no clear X-ray/optical link and the current data do not validate the scenario of X-rays arising in large-scale structures triggered at the surface of ζ Pup, traced by the optical data.

To solve the riddle of ζ Pup X-ray variations in the future, performing repeated, long-term, high-cadence, and multiwavelength monitoring will be needed. Only a multiwavelength study, which is best performed by combining X-rays to UV and optical data, can unveil the wind properties of the star at a given time. Such a comparison also needs to be long enough (several days/weeks, not a single snapshot) so that a clear view of the long-term (weeks) behaviour is acquired. Finally, since the variability properties of ζ Pup are not constant with epoch (>weeks), repeated observations are required, so that conclusions drawn at one time can be checked and fully validated.

Acknowledgments

Y.N. acknowledges support from the Fonds National de la Recherche Scientifique (Belgium), the Communauté Française de Belgique, the PRODEX XMM-Newton and Integral contracts, and an ARC grant for concerted research actions financed by the French community of Belgium (Wallonia-Brussels Federation). T.R. acknowledges support from the Canadian Space Agency grant FAST. A.F.J.M. is grateful for financial support from NSERC (Canada) and FQRNT (Quebec). We thank Kim Page, from the Swift team, for her kind assistance, the XMM-Newton team for their help, the referees for valuable comments that improved the paper, and the whole XRT team for their dedication towards making these Swift observations possible and usable. ADS and CDS were used for preparing this document. Based on data collected with Swift, the ESA science mission XMM-Newton (an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA), the SMEI/Coriolis instrument, and the BRITE-Constellation satellite mission (designed, built, launched, operated, and supported by the Austrian Research Promotion Agency – FFG, the University of Vienna, the Technical University of Graz, the Canadian Space Agency – CSA, the University of Toronto Institute for Aerospace Studies – UTIAS, the Foundation for Polish Science & Technology – FNiTP MNiSW, and National Science Centre – NCN).

References

All Tables

All Figures

	Fig. 1 Count rates in several energy bands of ζ Pup, and associated hardness ratios for XMM-Newton data from 2003–2006. The phases were calculated using P = 1.780938 d and T0 = 2 450 000 (Howarth & Stevens 2014). To facilitate the comparisons, the ratio between upper and lower limits of the y-axes are all equal.
In the text

	Fig. 2 Same as Fig. 1 but for XMM-Newton data from 2007–2012.
In the text

	Fig. 3 Same as Fig. 1 but for XMM-Newton data from 2013–2017.
In the text

	Fig. 4 Same as Fig. 1 but for Swift data; those from October 2016, December 2016, January 2017, and April 2017 are shown as red triangles, black dots, blue crosses, and cyan squares, respectively.
In the text

	Fig. 5 Hardness ratios as a function of the source brightness for all XMM-Newton data and for all Swift data.
In the text

	Fig. 6 Comparison of the pn (left) and RGS (middle and right) spectra observed at minimum brightness (first part of Rev. 1814; solid black line) and maximum brightness (first part of Rev. 2183; dotted red line). Blue segments above lines of RGS spectra indicate the intervals chosen for the moment calculations.
In the text

	Fig. 7 Left: same as Fig. 4 but for the fluxes (Table 4). The green stars correspond to data from February 2013. Middle and right: comparison between the XMM-Newton light curve of April 2017 (Rev. 3172, black dots and solid lines) and the (arbitrarily scaled) simultaneous Swift light curve (red crosses and dotted lines) derived from count rates (middle) or fluxes (right).
In the text

	Fig. 8 Subgroup of XMM-Newton light curves can be phased with a period of 1.78065 d to form a relatively coherent behaviour (top), but considering all available light curves only lead to scatter plots, whatever the value of the 1.78 d period considered (bottom).
In the text

Fig. 9 Comparison of the X-ray light curves with the simultaneous BRITE photometry in April 2015 (combined light curve from the two BRITE filters, for details see Ramiaramanantsoa et al. 2018), and December 2016–April 2017 (only one satellite, BRITE-Lem, in the blue filter, see Ramiaramanantsoa et al. 2018, andin prep.). For Swift, fluxes from the spectral fits (Table 4) are used since they are more reliable.

In the text

	Fig. 10 Same as Fig. 9 but using epochally averaged BRITE light curves valid for the date of the X-ray observations (cycles associated with “Part V” data for April 2015, “Part VI” for Dec. 2016, “Part VII” for Jan. 2017, and “Part IX” for April 2017, for details see Ramiaramanantsoa et al. 2018). The ephemerides are those adopted by Ramiaramanantsoa et al. (2018, and in prep.).
In the text

	Fig. 11 Comparison of the X-ray light curves with the average 1.78 d photometric cycle recorded with SMEI ± 15 d around the X-ray observing date. The ephemerides are those of Howarth & Stevens (2014): P = 1.780938 d and T0 = 2 450 000.
In the text