© ESO, 2012

1. Introduction

With its very early spectral type (O4Infp, Walborn 1972) and a distance of only 335 pc (van Leeuwen 2007; Maíz Apellániz et al. 2008), the star Naos, better known as ζ Puppis (or HD 66811), is one of the closest and brightest massive stars. It is therefore one of the most studied objects in the O-star population. However, despite the intense work, many open questions remain on its nature.

Indeed, ζ Puppis displays several intriguing properties. First, its visible spectrum shows clear signs of helium overabundance and chemical enrichment by CNO-processed material (e.g. Pauldrach et al. 2001), as well as fast rotation (more than 200 km s^-1 for vsin(i), Penny 1996; Howarth et al. 1997). Second, it is a known runaway (e.g. from Hipparcos data, Moffat et al. 1998). These properties have led to speculations on its evolutionary status. On the one hand, the chemical enrichment and fast rotation could result from mass and angular momentum exchange through Roche lobe overflow in a binary. The star ζ Puppis could therefore have been the secondary component of such a system; and then the supernova explosion of its companion would have ejected it from its birthplace a few million years ago (van Rensbergen et al. 1996). On the other hand, ζ Puppis displays a similar Hipparcos parallax to stars in the Vela R2 association¹ (Schaerer et al. 1997), and dynamical interactions within this association could have led to the ejection of the (single) O-star (van Rensbergen et al. 1996). In this scenario, the chemical enrichment of ζ Puppis would be explained by the intense rotational mixing occurring in the fast-rotating, main-sequence progenitor (Meynet & Maeder 2000). In addition, ζ Puppis displays double-peaked emission lines, supposedly arising in a rotating wind (Conti & Leep 1974; Petrenz & Puls 1996), and a compression of the wind in the equatorial plane was detected by Harries & Howarth (1996).

Due to its brightness, ζ Puppis was one of the first massive stars observed with high resolution in X-rays (Kahn et al. 2001; Cassinelli et al. 2001). At first, its X-ray lines appeared to match expectations as they did show the broad, blueward-skewed profiles expected for the wind-embedded shock model (Owocki & Cohen 2001). However, the devil was in the details. When quantitatively fitting the line profiles, Kramer et al. (2003) found a much lower wind attenuation than expected on the basis of the mass-loss rate determined from optical and UV observations (see also Oskinova et al. 2006). They also found that the typical optical depths τ_∗ used in the wind-shock models seemed independent of wavelength, which can only be explained by invoking porosity (Feldmeier et al. 2003; Oskinova et al. 2006). To improve the fitting of the X-ray line profiles, Leutenegger et al. (2007) included the effect of resonance scattering: better fits were indeed obtained, without the need of a large reduction in the mass-loss rate, but they also show that some unexplained discrepancies remain. Reanalyzing the Chandra data of ζ Puppis, Cohen et al. (2010) argue for a reduced mass-loss rate, without the need of any porosity because their new derivation of the optical depths implies an increase with wavelength, as expected from the bound-free absorption opacity of the (cool) wind. Except for Leutenegger et al. (2007), all these studies relied on a single 68 ks Chandra observation or a 57 ks XMM-Newton exposure taken in 2000. Both facilities have their advantages: while XMM-Newton has higher overall sensitivity, Chandra has lower background and higher spectral resolution and sensivity at short wavelengths for its grating spectra. Today, however, much more data are available (see below).

Considering the uncertainties in the line profile results and the lack of new variability studies, we decided to reinvestigate ζ Puppis using the best dataset available at the present time: 18 XMM-Newton exposures, corresponding to an exposure of 1 Ms totaling  >700 ks of useful time (i.e., an improvement by an order of magnitude compared to most previous studies). This dataset thus provides the most detailed X-ray view of an O star to date. The results that we obtained will be presented in a series of papers. The first one will present the data, their reduction, and a few first results; the second will focus on the X-ray variations of ζ Puppis, using EPIC and RGS data; the last one will present a global analysis of the high-resolution X-ray spectrum, using the merged high-resolution data.

This first paper is organized as follows. The dataset and its reduction are presented in Sect. 2, the spectral fits are presented in Sect. 3, the individual line profile fitting in Sect. 4, and the results are summarized and discussed in Sect. 5.

2. XMM-Newton observations

In the past decade of XMM-Newton observations, the star ζ Puppis was observed 18 times, mostly for calibration purposes. These datasets are excellent for studying the variability of ζ Puppis since (1) the scheduled exposure times were often long (up to  ~60 ks) and (2) the observing dates probe weekly, monthly, and yearly timescales. Unfortunately, many observations were affected by soft proton background flares, resulting in total exposure times reduced by about 30%. Total net exposure times for EPIC-MOS, EPIC-pn, and RGS amount to 579 ks, 477 ks,  and 751 ks. A summary of the observations is given in Table 1. The successive columns provide the dataset ID (obsID and revolution number); the date at mid-exposure (in the format dd/mm/yy + UT time and JD–2 450 000, calculated using the start/end times of the observations listed in the on-line “Observation lokator”); the mode as well as the scheduled, performed, and effective (i.e., after cleaning flares) exposure time for both EPIC-MOS and EPIC-pn instruments; the scheduled, performed, and effective exposure time for RGS. An empty column indicates a discarded or unavailable dataset (see below). The target was placed off-axis in two observations (5.95′ off-axis in Rev. 0731 and 1.1′ off-axis in Rev. 0903).

2.1. EPIC data

The EPIC data were reduced with SAS v10.0.0 using calibration files available on January 1, 2011 and following the recommendations of the XMM-Newton team². Different modes (small-window, large-window, full frame, timing), as well as different filters (thick, medium) and position angles were used for these observations, resulting in a somewhat heterogeneous dataset. To ensure the most homogeneous analysis, hence a meaningful comparison between datasets, two decisions were taken. First, a few observations were discarded: those when the instruments were not “on” (aka CAL   CLOSED), those totally affected by flares, those with the source appearing totally or mostly in a CCD gap, those with very short exposure times (<10 ks), those taken in timing mode, and those using a unique combination of mode + filter. This trimming process results in a final dataset composed of nine observations taken with large-window + thick filter and six observations with small-window + thick filter for EPIC-MOS. For EPIC-pn, the final dataset is composed of ten observations taken with small-window + thick filter, five observations taken with large-window + medium filter mode and four observations taken with large-window + thick filter (see Table 1). Second, the extraction regions were chosen to be as constant as possible. A single source region was used for EPIC-MOS, regardless of the mode, but two background regions were defined, one for each mode since it was not possible to extract the background on the same CCD chip as the source for the small-window mode. For EPIC-pn, a single source region was used for all modes; a single background region was used for the large-window mode, but four different background regions were necessary for the small-window mode. Table 2 gives the position and shape of each of these regions.

2.1.1. Pile-up

The star ζ Puppis is rather bright in X-rays, with EPIC count rates of  ~2 cts s^-1 (for MOS) and  ~6.5 cts s^-1 (for pn with thick filter). These count rates are at the pile-up limit from the XMM Users’ handbook for the large-window modes (i.e., 1.8 and 6 cts s^-1 for MOS and pn, respectively) but well below the limits for the small-window mode (i.e., 5 and 50 cts s^-1 for MOS and pn, respectively). Some pile-up may thus affect our EPIC data taken in the large-window mode. To see how severe the pile-up is, we performed several checks. First, we inspected the event files: no event with PATTERN = 26–29 was found. The pile-up is thus moderate. Second, we run the SAS task epatplot: some small but significant deviation from the “no pile-up” configuration is detected for the large-window data, especially those taken with the medium filter.

To get rid of the pile-up, we could extract the data in an annulus centered on the source. This would require the annulus to be perfectly centered on the source, especially since the PSF is far from symmetric. Getting the exact position of ζ Puppis in the datasets is, however, an impossible task. Indeed, the detection algorithm is also disturbed by pile-up, so that the position found in this way for ζ Puppis is not accurate. For example, the pipeline-processed data from Rev. 1620 yields a position some 2.5′′ away from the Hipparcos position of ζ Puppis (this separation is the maximum found in the pipeline-processed data), while a dedicated run of the detection algorithm using only the PATTERN = 0 events gives a position that is only 1.2′′ away. One would immediately think of using nearby X-ray sources associated with well-known stars – sources that are less bright in X-rays (thus unaffected by pile-up) but still bright enough to get an accurate position in each dataset. However, such nearby sources do not exist in the neighbourhood of ζ Puppis. A perfect centroiding of annular regions is thus impossible, and always using the Hipparcos position for annular regions may alter the source’s properties, hence the results of the variability study that we wish to perform.

On the other hand, there is an alternative way in cases of mild pile-up: using only PATTERN = 0 events. After such a filtering, we compared the PATTERN = 0 spectra extracted in a circular region with the spectra extracted in annular regions (hence free of pile-up) using the usual PATTERN = 0–12 for MOS and PATTERN = 0–4 for pn. This check was done by fitting simple two-temperatures models on LW+medium data from Rev. 0156, LW+thick data from Rev. 0731 and SW+thick data from Rev. 1814 for EPIC-pn, and on LW+thick data from Rev. 0156 and SW+thick data from Rev. 1814 for EPIC-MOS. The comparison is excellent for EPIC-MOS data: fluxes and count rates differ by  <1% and best-fit spectral parameters are within the errors. The remaining difference can be attributed to the slightly larger noise in the spectra extracted in annuli and from calibration differences. The comparison is less perfect for EPIC-pn, especially for the data taken with the medium filter: the flux differences reaches 6% in this case, and best-fit spectral parameters are at 2-σ from each other. The pile-up thus still has a small influence on the EPIC-pn data taken in the large-window mode, and those data should thus be considered with caution.

2.1.2. The final files

A final check was made on the data from Rev. 1620, which yields the most discrepant (2.5′′) position for ζ Puppis, if we trust the pipeline processing. We first derived the position of ζ Puppis from the PATTERN = 0 data using the SAS task edetectchain, and then extracted the spectra using a circular region centered on that position. We compared these spectra to those extracted on PATTERN = 0 data using a circular region centered on the Hipparcos position of ζ Puppis. Both sets of spectra appear identical in Xspec: the small centroiding errors thus have no impact on the spectra as long as a circular region is used.

We therefore extracted lightcurves and spectra of ζ Puppis in a circular region centered on the Hipparcos position of the target. We used only the PATTERN = 0 event files. While this is not necessary for the small-window mode, it ensures a homogeneous data reduction. Xspec v12.6.0 was used to fit spectra, and our own software to analyze EPIC lightcurves.

Spectra were grouped using the new SAS task specgroup. It enabled us to reduce the oversampling, which may “cause problems during spectral fitting because the spectral bins are then not completely independent” (excerpt from SAS 10.0.0, online documentation). We chose an oversampling factor of five, ensuring that no spectral bin is narrower than one fifth of the full width half maximum resolution at the central photon energy of the bin. While providing more statistically correct data, this process dramatically reduces the number of spectral bins. The data were also grouped to ensure that a minimum signal-to-noise of three was reached in each spectral bin of the background-corrected spectra.

2.2. RGS data

The RGS datasets were also reduced in a standard way with SAS v10.0.0. Many new, important RGS features were modified in that version (e.g. the spectral binning in wavelength rather than in dispersion angle units). This ensures better calibration of our datasets. It also solved the calibration problems (wavelength shift and reduced flux) found when using earlier versions of the SAS for the two observations where ζ Puppis was placed off-axis. The data were extracted using the proposal position of the source, which is the same in all observations but the first two (Revs. 0091 and 0156, shift of 0.0002° in both RA and Dec) – this small position shift has no impact on the derived RGS spectra.

When detected, flares were discarded using rgsfilter³. The tasks rgsspectrum and rgsrmfgen then provided unbinned source and background spectra, as well as response matrices for each order (1, 2) and each instrument (1, 2). A final, combined spectrum was also calculated using all 18 RGS datasets and the task rgscombine. The background files and matrix responses were attached to the source spectra using the new SAS task specgroup, which we also use to ensure an oversampling factor of maximum five (see above).

Fluxed spectra combining both RGS instruments and both orders were obtained using the task rgsfluxer. A correction for off-axis angles is applied to ensure that the fluxes are real photon fluxes and not simply recorded count rates (i.e., the arf response matrix is fully taken into account). The spectra of one revolution were sampled to get 1500 spectral bins, while the spectrum combining the 18 datasets was also calculated to get 3000 spectral bins. This ensures an oversampling factor of about three and six for the former and latter cases, respectively. Two caveats should be noted. First, the rmf matrix is not fully taken into account by rgsfluxer, and there is no instrumental width correction. This needs to be accounted for when modeling the spectra (see Paper III, in prep.). Second, there are known small wavelength shifts in RGS spectra, apparently depending on the Sun’s aspect angle. However, no sign of such an effect is detected in our dataset when we use SAS v10. They would appear as small spectral variations, and there are none, see Nazé et al. (2011), hereafter Paper II.

3. EPIC spectra

With data of such high quality, the error bars on the spectra are very small, making it very difficult to get a formally acceptable fit. We thus avoided trying to get a perfect fit (i.e., χ² ~ 1), which is actually impossible to get without going into unrealistic, overcomplicated models (e.g., 10 components fits with independent, free abundances) when all instruments agree. Rather, we tried to get a fit that is as simple and realistic as possible and as close to the spectral data as possible.

Zhekov & Palla (2007) show that the high-resolution Chandra spectrum can be fitted by a shock model where the dominant temperature is 0.1–0.2 keV (ζ Puppis was actually the star with the coolest dominant plasma), though some contribution from plasma with 0.3–0.7 keV was needed to achieve a good fit. This modeling clearly shows that the plasma in ζ Puppis is fairly cool. We decided to use a simpler formalism, which does not try to reproduce shocks, but simply considers the addition of optically thin thermal plasma (without any assumption on their origins). In Xspec, the models using a distribution of thermal plasma (e.g. c6pvmkl) fail to provide a fit close to the data, so we had to fit the data using a sum of individual optically thin plasmas. As could be expected, one, two, or three temperature fits do not provide fits close to the data, and similar conclusions are reached for four-temperature fits with a single absorption (which is unsurprising in view of the c6pvmkl result). Our model of choice is thus a four-temperature model with individual absorptions: tbabs ×  ∑ vphabs × vapec, where the first component represents the interstellar absorption, fixed to 8.9 × 10¹⁹ cm^-2 (Diplas & Savage 1994), and the abundances of the absorption and emission components are assumed to be equal. Adding a fifth thermal component (e.g. at 1 keV) does not significantly improve the fit, so we stuck to the decision of using four components.

After the best-fit temperatures and absorptions were found, we investigated the impact of using nonsolar abundances. We began by keeping the helium abundance to two times solar (Repolust et al. 2004), and let the nitrogen and oxygen abundances vary. We perfomed a simultaneous fit to all available EPIC spectra for each revolution. A few conclusions could be drawn from this trial. First, the pn spectra obtained in large-window mode with the medium filter (Revs 0156 and 0552, and second pn observation of Revs 535, 538, and 542) clearly deviate from other data, showing the impact of pile-up. Since all other data (pn or MOS, LW, or SW modes with thick filter) overall agree with one another, we only discarded the pn data taken with medium filter from further analysis. Second, the absorptions and temperatures of the fits do not vary much, so we fixed them to kT_i of 0.09, 0.27, 0.56, and 2.18 keV and N_H,i of 0.1, 0.1, 0.71, 0.  × 10²² cm^-2. Finally, the origins of the high χ² can be pinpointed better. On the one hand, MOS and pn data do not always agree, especially at 0.4 keV (see Fig. 1). This explains why some fits deviate from the mean behaviour. Revs 0731 and 1071 only provided pn data, while Revs 0156 and 1164 only consist of MOS data. This difference between MOS and pn probably come from the remaining cross-calibration problems, not pile-up, since ζ Puppis is far from the pile-up limit for data taken in SW mode with the thick filter. On the other hand, all three instruments (MOS1, MOS2, and pn) sometimes display δχ² of the same sign. This is often the case near strong lines, which points toward two explanations: (1) atomic parameters are imperfect (even for APEC – Fig. 1 shows two examples: at 1.24 keV, there is flux in the data but not the model – a line is probably missing; at 1.5 keV, there is a line in the model but not the data) and (2) the asymmetric line shape of ζ Puppis influence even EPIC data. We cannot do much for the former, and the latter will actually be studied in detail in the third paper of this series, so it is beyond the scope of the present contribution.

Fig. 1 The EPIC pn (top, in green) and MOS (bottom, in red and black) data of Rev. 1983 superimposed on the best-fit model.

The next step is to free more abundances, but this could yield erratic, problematic, or even just not better results. For example, the carbon abundance clearly becomes unrealistic (whether freed last or first). Indeed, low-resolution, broad-band spectra, such as those taken by EPIC, yield few constraints on the carbon abundance and, as often happens with such data, the fitting procedure favors high values of carbon enrichment, regardless of its actual value. From previous studies, it is well known that nitrogen is overabundant in ζ Puppis, and carbon may be subsolar (e.g. 0.35 times solar in Pauldrach et al. 2001; and 0.6 times solar in Oskinova et al. 2007). However, reported abundances vary quite a lot in the literature (see other examples below) and, following Zhekov & Palla (2007), we decided to keep it to solar as for other unconstrained elements. For other elements, freeing the abundance may not improve the fitting quality or may not yield abundances significantly non-solar: this was the case of magnesium, sulphur, iron and neon when released one after another. The silicon abundance stays close to solar but does improve the χ², it was thus allowed to vary together with He, N, and O.

The best-fit results are shown in Table 3. The parameter errors are taken from the raw fitting results of Xspec (i.e., these 1-σ errors are “calculated from the second derivatives of the fit statistic with respect to the model parameters at the best-fit”, and are indicative, see Xspec manual). No error is provided for the fluxes. We could indeed use the relative error on count rates, which would yield relative errors in the 0.1–0.3% range due solely to the Poisson noise. However, relative flux differences between MOS and pn calibrations amount to about 1% in the soft energy band (where they are maximum), and the use of other similar models also yield 1% relative errors. Fluxes should therefore be considered as determined to 1/100, not 1/1000, uncertainties.

The listed parameters should not be overinterpretated because they simply represent a convenient way of fitting the data well, no more no less. The abundances, in particular, are indicative, and high-resolution data provide much more stringent constraints (see Paper III, in prep.). It is, for example, interesting to note that freeing all abundances at the same time could lead to a better χ², but to totally erratic and unrealistic abundances (carbon becoming largely overabundant, nitrogen being solar). With this in mind, it is quite remarkable that our results, which should be considered as only indicative, agree fairly well with previous abundance determinations. The helium abundance of ζ Puppis was found to be 1.2 and 3.4 times solar⁴ by Pauldrach et al. (2001) and Oskinova et al. (2006), respectively, and we found an average value of  ~2 (though this abundance is only weakly constrained in EPIC spectra). The nitrogen abundance was determined to be 1.7, 6, and 8 times solar by Zhekov & Palla (2007), Oskinova et al. (2006), and Pauldrach et al. (2001), respectively, and we again found an average value of  ~4. The oxygen abundance is the least constrained of all, with values of 0.75, 1.6, and even 0.16 times solar reported by Pauldrach et al. (2001), Oskinova et al. (2006), and Zhekov & Palla (2007), respectively. The last agrees well with our value, but it is formally unconstrained since its error was 0.23. Our silicon determination also agrees well with the value found by Zhekov & Palla (2007).

A few general conclusions can be drawn. The dominant components are those with temperatures of 0.09, 0.27 keV, and 0.56 keV (providing 20%, 43%, and 35% of the total flux, respectively), corroborating the conclusion found by Zhekov & Palla (2007), with a different formalism, that cool plasma dominates in ζ Puppis. The flux and spectral parameters do not vary much, with the largest variations obtained when only one instrument is available (e.g. Rev. 0156, where the pn data are discarded because of the piled-up associated with the medium filter). Fluxes (and their associated dispersions) in the total (0.3–4. keV), soft (0.3–0.6 keV), medium (0.6–1.2 keV), hard (1.2–4. keV), and Berghöfer’s (0.9–2. keV) bands⁵ are 15.3  ±  0.4, 4.15  ±  0.14, 8.07  ±  0.26, 3.06  ±  0.09, 5.67  ±  0.17  ×  10^-12 erg cm^-2 s^-1, respectively. Dispersions typically amount to 3%, slightly more than the typical 1% error, and a shallow decreasing trend is detected, with a decrease of only  ~5% in flux since Rev. 0156. This trend is reminiscent of aging detector sensitivity problems, and we cannot exclude this possibility on the sole basis of the XMM-Newton dataset.

Fig. 2 Spectral parameters as a function of Revolution number.

4. Wind profiles

Cohen et al. (2010) report small variations in the line profiles with wavelength, due to the energy-dependent opacity of the cool wind. However, the Chandra observation that they used was relatively short, hence subject to a much higher noise on the spectrum than for our data, and it is also rather insensitive to wavelengths  >20 Å, where the effect is expected to be the strongest. We therefore investigate here the issue again using the combined RGS spectrum, which has a better signal-to-noise ratio and extends beyond 20 Å.

The left-hand panel of Fig. 3 shows the observed Lyman α lines in velocity space. For the figure, the lines were approximately continuum-subtracted and normalized to have a peak amplitude unity, to highlight their differences. Neighboring lines can be seen for some of these Lyman α lines as bumps in the blue or red wings. We do not show the N vii Lyα line because it is blended with an N vi line. The variations with wavelength are obvious, since the peak velocity clearly appears less blueshifted for the short-wavelength lines. The comparison of width and skewness is more difficult by eye, because the RGS resolution broadens the short-wavelength lines in velocity space, blurring the trends.

To quantify the wavelength variations, we fit the lines with the same models⁶ as in Cohen et al. (2010), to ensure homogeneity. Results are provided in Table 4: the first two columns identify the considered line, the next three columns define the line shape (characteristic continuum optical depth τ_∗, radius R_∗ for the onset of the X-ray emission, and the strength of the line), and the last two columns provide details on the line ratios in the He-like fir triplets.

Several things must be noted. First, in two cases (the He-like triplets of N and O), resonance scattering was needed to achieve a good fit. Second, the Lyman α line of nitrogen is blended with a line from N vi. We fit these two lines together, assuming that the line profile parameters are identical: keeping them independent yields unrealistic results (τ_∗ ~ 0) for the weak N vi line, and the lines are so blended that little independent information is available, explaining the apparently strange results for the weakest line. The achieved fit of the nitrogen blend is far from perfect, however (χ² ~ 2). Third, the iron line at 15 Å is not very well fitted (χ² ~ 2), despite our efforts. It seems that line blends (there are numerous Fe xviii lines in the neighborhood) affect the profile. Although these lines are weak, the very low noise of our data reveals their impact, which was not obvious in the Chandra data. Finally, the fitting was done using a single terminal velocity for the wind (2250 km s^-1, Puls et al. 2006), a β = 1 exponent for the velocity law, and a power law of zero slope to represent the continuum.

Fig. 3 Left: line profiles in velocity space, of the observed Lyman α lines. Right: variations in the mean optical depth τ∗ with wavelength. Results from this work are shown with filled circles and hexagons – the latter for fits including resonance scattering, while Cohen et al. (2010) results are displayed with empty triangles. The Cohen et al. theoretical absorption is shown by a dotted line.

Because it is always compatible with zero and yields no improvement of the fits, no global line profile shift was applied, except for the iron lines near 17 Å where the improvement is significant with a shift of only –15  ±  2 km s^-1. Wavelength shifts between instruments/orders were considered, but they yielded (i) erratic results generally without any improvement of the fit, (ii) sometimes high values (incompatible with wavelength calibration reports), for example, for the iron line at 15 Å, (iii) values compatible with zero within 3 sigma, and (iv) values for one instrument/order compatible with those of another instrument/order within 3 sigma. Cross-correlation also suggests zero shifts between exposures and between different instrument/order combinations of the same exposure. Due to the broadness of its lines, ζ Puppis is indeed not the best source to find such shifts: to do so, the XMM-Newton calibration team uses sources with unresolved X-ray lines.

As shown in the right-hand panel of Fig. 3, the optical depth varies with wavelength. Our fitting thus confirms the Cohen et al. (2010) preliminary results by extending them to longer wavelengths and by decreasing the noise on most lines (except the bluest ones). The wind of ζ Puppis therefore is unlikely to be composed of clumps that are fully opaque at all wavelengths, as suggested for some porosity models. Indeed, optically thick clumps should produce a gray opacity that would be determined solely by the geometry of the clumps. In principle, these optical depth variations may be used to constrain the mass-loss rate, for a given star+wind model. This was attempted in Cohen et al. (2010), and they find that a reduced mass-loss rate with nonsolar abundances provides the best-fit results. For comparison purposes, the same theoretical opacity is shown in Fig. 3: below 20 Å, it indeed provides a good fit, but for longer wavelengths, the agreement is less good because of the nitrogen edge at 26 Å – an edge which appears very strong as the nitrogen abundance is enhanced for ζ Puppis. This discrepancy may be solved by adapting the value of the abundances and mass-loss rate of ζ Puppis, but this is beyond the scope of this paper. Clearly, additional modeling is needed to reproduce the full behavior of ζ Puppis in X-rays (to be done in Paper III, in prep.).

In contrast, the onset radius appears remarkably stable and confined in the 1–1.5 R_∗ range, which is quite usual in massive stars (Güdel & Nazé 2009, and references therein). The sole exceptions are the Nitrogen triplet and Lyα lines, but it would not be surprising that, as these lines have emissivities that peak at rather low temperatures, it simply indicates a formation farther out in the wind.

5. Conclusion

In the past decade, about 1 Ms of data were obtained by XMM-Newton on ζ Puppis. Of these, about 30% is strongly affected by flares, reducing the useful exposures to 579 ks for EPIC-MOS, 477 ks for EPIC-pn, and 751 ks for RGS. A variety of modes were used, the most reliable being the SW+Thick Filter mode. The use of the medium filter yielded piled-up data, while the data taken in LW+Thick Filter mode do not appear significantly different from those obtained with the SW+Thick Filter mode. Attention was paid to this problem, notably by choosing similar extraction regions for all datasets and using only single events.

Broad-band EPIC data taken with the thick filter were analyzed using absorbed, optically thin thermal emissions. Four temperatures are needed to reproduce the data in a reasonable way. The fits are not formally acceptable, but (i) instruments do not always agree with one another, and (ii) the reduced noise amplifies the limitations due to the imperfect atomic parameters and standard line profiles. Nevertheless, the EPIC spectra appear remarkably stable over the decade of observations, with only 3% dispersions around the average fluxes. A detailed variability study, based on lightcurves, will be presented in Paper II.

The combined high-resolution RGS spectrum confirms that the X-ray line profiles vary with wavelength. Fitting individual line profiles using a wind model yields similar onset radius for

the X-ray emission, but wind continuum opacities depending on wavelength. This is simply because the cool absorbing clumps in the wind are not fully optically thick at all wavelengths, though further modeling is needed to adequately reproduce the opacity variations.

Online material

Acknowledgments

Y.N. acknowledges L. Mahy for his help with CMFGEN, and the XMM-Newton helpdesk for interesting discussions about the data. Y.N. and G.R. acknowledge support from the Fonds National de la Recherche Scientifique (Belgium), the Communauté Française de Belgique, the PRODEX XMM and Integral contracts, and the “Action de Recherche Concertée” (CFWB-Académie Wallonie Europe). C.A.F. acknowledges support from the FNRS-Conacyt agreement, as well as from the “Programa de becas para la formacíon de jóvenes investigatores de la DAIP-GTO”. ADS and CDS were used for preparing this document.

References

All Tables

All Figures

	Fig. 1 The EPIC pn (top, in green) and MOS (bottom, in red and black) data of Rev. 1983 superimposed on the best-fit model.
In the text

	Fig. 2 Spectral parameters as a function of Revolution number.
In the text

	Fig. 3 Left: line profiles in velocity space, of the observed Lyman α lines. Right: variations in the mean optical depth τ∗ with wavelength. Results from this work are shown with filled circles and hexagons – the latter for fits including resonance scattering, while Cohen et al. (2010) results are displayed with empty triangles. The Cohen et al. theoretical absorption is shown by a dotted line.
In the text