Free Access
Issue
A&A
Volume 551, March 2013
Article Number A83
Number of page(s) 12
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/201219734
Published online 26 February 2013

© ESO, 2013

1. Introduction

Through their powerful stellar winds and their violent deaths as supernovae, massive stars have a tremendous impact on the ecology of their galaxies. The feedback by the stellar winds of early-type stars is indeed a major ingredient of the chemical evolution of the interstellar medium. Moreover, the mass-loss has an important impact on the evolution of the star. The quantity of material released by the wind as well as its chemical composition depend on the evolutionary stage of the star. A proper understanding of the stellar winds of massive stars is therefore fundamental for our knowledge of the stellar evolution and the impact on the interstellar medium. However, the determination of the mass-loss rate of a star is a complicated task.

Over the past decade, the importance of the structures of the stellar winds and their impact on the determination of the mass-loss rates became more and more obvious. Observationally, these structures manifest themselves through low-amplitude variability of emission lines such as Hα (Eversberg et al. 1998), and through difficulties to fit the observed spectra with model atmosphere codes that assume homogeneous winds (Bouret et al. 2005). These small-scale structures are usually attributed to hydrodynamic instabilities intrinsic to the winds (Feldmeier et al. 1997).

Accounting for the fragmented structure of the winds leads to a downward revision of the mass-loss rates by a factor between a few and ten. Owing to the different sensitivities of the spectra over various wavelength ranges to the structures, multi-wavelength analyses are required to determine consistent mass-loss rates of massive stars. For instance, Hα emission and free-free radio emission are collisional processes that scale with the wind-density squared and are thus sensitive to the clumping factor regardless of the size of the clumps. In clumped winds with a void inter-clump medium, the wind’s velocity law misses some values along a given sight line. This porosity in the velocity field is called “vorosity” (Owocki 2011). The latter modifies the Sobolev length and its impact is mainly visible on the UV resonance lines. In a fragmented wind with non-optically thin clumps, X-ray lines are sensitive to porosity, which depends on the clump continuum optical depth, and thus clump size as well as the clumping factor.

The X-ray line profiles can also give information on the hydrodynamical instabilities in the wind. In fact, considering the wind-embedded shock scenario, the emission of X-ray photons arises when two fragments with different velocities collide (Lucy & White 1980). The observable morphology of an X-ray line depends on the interplay of emission by the hot plasma and photoelectric absorption by the cool wind component. The latter effect is obviously also sensitive to the presence or absence of structures in the wind. On the other hand, the location of the X-ray emission region as well as the temperature of the X-ray emitting plasma give important information on the hydrodynamical properties of the winds.

ζ Pup (=HD 66811, O4 Ief) has been the target of numerous studies over a wide range of wavelengths. This star is often considered as a cornerstone object for the understanding of massive stars. The analysis of the far UV – UV – optical spectra by Bouret et al. (2012) confirmed a super-solar total CNO abundance, with N being highly overabundant and C being strongly depleted. On the other hand, Pauldrach et al. (2012) found quite different results in their far UV – UV spectral modeling, indicating an over-abundance of both nitrogen and carbon and a strong depletion of oxygen. They also determined a homogeneous mass-loss rate (~13 × 10-6M yr-1) about six times higher than the clumped mass-loss rates determined by Bouret et al. (2012) and Najarro et al. (2011). This situation is embarrassing as one would aim for consistent parameters obtained from studies using different wavebands and different modeling tools.

ζ Pup is also a bright X-ray source and has been observed both with Chandra and XMM-Newton at high spectral resolution. The star was indeed regularly observed with the RGS instrument onboard XMM-Newton for calibration purposes, yielding a combined data set of unprecedented quality (Nazé et al. 2012). Previous studies of ζ Pup’s high-resolution X-ray spectrum revealed new interesting results (Cassinelli et al. 2001; Kahn et al. 2001; Oskinova et al. 2006; Leutenegger et al. 2006, 2007; Cohen et al. 2010). For instance, the analyses by Cassinelli et al. (2001) and Leutenegger et al. (2006) of the forbiden,inter − combination,resonance (fir) triplets of He-like ions revealed that the inner boundary of the X-ray emission region is rather close to the star (at about 1.5 R and possibly even closer for the hottest plasma, Cassinelli et al. 2001). However, these studies were conducted on individual line profiles that were considered separately. Each specific line was in some sense assumed to be formed by a plasma at the temperature of maximum emissivity of the line. This approach thus neglects the fact that several plasma components (each one characterized by its specific temperature, emission measure and location inside the wind) may contribute significantly to a given line. To the best of our knowledge, there have been no previous attempts to fit the entire high-resolution X-ray spectrum (either RGS or HETG) of an O-type star at once by properly combining the contributions from several hot plasma components and modeling the absorption by the cool wind, accounting explicitly for the geometry of the problem.

This paper thus explores, for the first time, the possibility to fit all the lines of an X-ray spectrum between 6 and 35 Å simultaneously and coherently, assuming a discrete distribution of plasma components, characterized by their own temperature and location in the wind. To do so, we used the code developed by Hervé et al. (2012), which accounts for the intrinsic emissivity of the hot plasma as a function of temperature and chemical composition. We furthermore introduce a multi-temperature plasma spectral fitting procedure. The goal is to determine the properties of the X-ray emission region accounting for the most recent results from other wavelength domains. These properties can ultimately be compared to the predictions of existing and future hydrodynamical simulations of the instabilities of stellar winds.

This paper is organized as follows. In Sect 2, we present a short summary of the ζ Pup RGS dataset. Then we briefly recall the main assumptions in our modeling tool (Hervé et al. 2012) concerning the intrinsic emissivity and the absorption formalism (Sect. 3). We also present, in Sect. 4, the different physical processes we have included in the code to better reproduce the observational spectrum of ζ Pup. Finally, in Sect. 5, we discuss the results that we obtain, especially for the composition and the location of the different plasma components in the wind.

2. Observational data

Nazé et al. (2012) extracted and reduced 18 observations of ζ Pup with XMM-Newton, which yield more than 700 ks of useful exposure. For the details of the data reduction, we refer the reader to the paper of Nazé et al. (2012, Paper I). Using this dataset, Nazé et al. (2013, Paper II) performed an analysis of the variability of the X-ray emission, finding no significant time variability of the line intensities and profiles. Therefore, we can combine the individual RGS spectra into a single, high-quality, fluxed spectrum. The calibration of the fluxed spectrum properly accounts for the dependence of the effective area as a function of energy, but does not account for the redistribution of monochromatic response into the dispersion channels. For intrinsically narrow lines, this might introduce artefacts in spectral regions of strong instrumental sensitivity changes. However, in the case of ζ Pup, the intrinsic width of the lines is large and therefore such effects should be of very minor impact here. The instrumental width of the RGS is accounted for in our fitting process, because the synthetic spectra are broadened according to the instrumental profile.

Whilst the Chandra-HETG spectrum of ζ Pup (Cassinelli et al. 2001) has a superior spectral resolution and extends to shorter wavelengths than the RGS spectrum analysed in the present work, it has a considerably lower signal-to-noise ratio than the RGS spectrum. We therefore decided to focus our analysis on the RGS spectrum only.

3. X-ray spectral modeling code

3.1. Basic principles of the code

In order to analyse high-resolution X-ray spectra of massive stars, we have developed our own dedicated code (Hervé et al. 2012). In our model, we assume that the wind hosts several hot plasma components. Each plasma component extends between an inner and an outer boundary and has a specific kT that is constant over that region. The chemical composition is taken to be the same for all components. We then discretize the X-ray emission region of a given plasma component into a large number of cells. We compute the X-ray emissivity of each cell as a function of its temperature and the quantity of material inside the cell. For this purpose we use version 2.0.1 of the AtomDB database (Foster et al. 2012) with emissivities computed using the APEC code (Smith et al. 2001)1. This code calculates the line emission arising from many physical processes such as collisional excitation, ionization, and dielectronic recombination in a collisionally ionized, optically thin isothermal plasma in thermal equilibrium. The continuum emission is also included in the form of the bremsstrahlung emission and the two-photon decay.

We assume the plasma to be embedded in the wind and to move along with it. The velocity of the wind is determined by a β-law, i.e. where v is the terminal velocity of the wind and R is the radius of the star. For each cell, we Doppler-shift its contribution into the observer’s frame of reference accounting for the projection of the wind velocity on the line of sight (MacFarlane et al. 1991).

Since the X-ray emitting plasma is assumed to be optically thin, an X-ray photon, once emitted, either escapes or is absorbed by the cool wind material along the line of sight. Hence, the radiative transfer problem is reduced to the calculation of the optical depth of the fragmented cool wind material, τ, along the line of sight from the emitting cell to the observer. For this purpose, we consider that, at the microscopic level, the cross section for photoelectric absorption in the X-ray domain is dominated by the bound-free transitions from the K and L energy levels.

If the clumps are optically thick, the effective opacity is determined by their geometrical properties (Feldmeier et al. 2003). The effective opacity can be written as (Owocki & Cohen 2006)

where l is the scale-length of the clump, mc is its mass, κ is the mass absorption coefficient and τc the optical depth of a clump. This latter can be written as: (1)where ⟨ ρ ⟩ is the mean density of the stellar wind at position r, f the volume filling factor of the clumps, is the mass-loss rate of the star, and v is the velocity of the wind at the radius r which is determined by a β-law (see above).

With this assumption, it appears that the effective reduction of the optical depth as a result of the fragmentation depends on the ratio between the clump scale and filling factor, also called the porosity length (Owocki & Cohen 2006). This result can be generalized from the optically thick to the optically thin limit using a simplified scaled effective opacity (Owocki & Cohen 2006): (2)Using this result in a steady state (i.e. non-variable) wind, we obtain the wind optical depth (Owocki & Cohen 2006) in conventional (p,z) coordinates: (3)where , is the characteristic optical depth of the wind, and is the porosity length. In the following we adopt h(r) = h′ × r.

A slightly different formalism to describe the fragmentation of the wind is based on the concept of the fragmentation frequency (Oskinova et al. 2006). As we have shown in Hervé et al. (2012), both formalisms yield very similar X-ray line profiles and the capabilities of the current generation of high-resolution X-ray spectrographs onboard XMM-Newton or Chandra do not allow us to choose between the two formalisms. Here we adopt the porosity length formalism because it leads to an analytic expression for the integral in the equation of the optical depth. This has the advantage to allow a faster computation. In Hervé et al. (2012), we have shown that the inclusion of a few percent of the wind material in a homogeneous inter-clump medium has little impact on the emerging X-ray spectrum. In the present work, we thus consider that the inter-clump medium is void.

In either formalism, the characteristic wind optical depth τ, and thus the mass absorption coefficient κ, play a key role in the computation of the absorption of X-rays by the stellar wind. To estimate the values of κ, we used the non-LTE, line blanketed, model atmosphere code CMFGEN (Hillier & Miller 1998) to compute the ionization structure of the cool wind as a function of temperature and mass-loss rate. In these ionization calculations, CMFGEN accounts for different physical processes such as collisional excitation, photo-ionization by photospheric light, and X-rays, as well as radiative and dielectronic recombination. Once we know the ionization structure, we can compute κ.

The results are illustrated in Fig. 1 at four different wavelengths. This plot shows that the approximation of a mass absorption coefficient that is independent of the position in the wind is not a valid one in this case. Indeed, the ionization structure of the wind is not constant throughout the entire wind and the impact on the value of the mass absorption coefficient is huge. The outer regions of the wind are less ionized, and consequently, the cool material produces a stronger absorption in the outer regions.

thumbnail Fig. 1

Radial dependence of the mass absorption coefficient at different wavelengths (specified in the upper left corner of each panel). The mass absorption coefficients are calculated with the elements listed in Table 1 with abundances found in our study.

Open with DEXTER

Compared to our previous version of the code (Hervé et al. 2012) where we adopted a value of κ that was independent of r, we have decided to render the treatment of the absorption of the wind more realistic. For this purpose, we approximated the radial dependence of κ by two simple relations that hold over two ranges of values of r. A linear relation (4)for r < Rlim, and a constant (5)for r ≥ Rlim. Here κ0 and κmax are the mass absorption coefficient in the innermost part and in the outer part of the wind and γ is the slope, respectively. Rlim is the position in the wind where κ reaches its plateau value κmax (see Fig. 1). The advantage of this parametrization is to provide a more realistic description of the variation of κ as a function of the radius and, at the same time, to preserve an analytical solution for the integral of the optical depth equation. Indeed, for a linearly increasing porosity length (h(r) = h′ × r), the optical depth is now given by (6)with .

thumbnail Fig. 2

Wavelength dependence of the mass absorption coefficient at different positions inside the wind (1.5, 3, 6, and 15   R as specified in the upper left corner of each panel).

Open with DEXTER

Table 1

Chemical elements and their ionization levels used to calculate the mass absorption coefficient κ and the UV radiation field in CMFGEN.

Note that we assumed in Eq. (6) that z is included in the interval . If this last condition is not fulfilled, the boundaries of these integrals are adjusted, but the principle remains the same.

The first two integrals have analytical solutions that facilitate the computation which is true for the last integral as well, which is formally identical to the one used in the original prescription of Owocki & Cohen (2006, see our Eq. (3)).

Figure 2 illustrates the wavelength dependence of the mass absorption coefficient at four different positions in the stellar wind. The strongest absorption edge in this figure (at about 26 Å) is due to N iv. The complex shape of this relation could affect some of the frequently used diagnostic tools based on line ratios. For instance, the ratio between the flux of N vii Ly α at 24.78 Å and the flux of the N vi triplets at 28.8−29.5 Å, which is used in line-by-line analyses as a temperature indicator, could be affected by the strong N iv edge.

3.2. The treatment of the fir triplets of He-like ions

In an optically thin thermal plasma at high temperatures, some elements retain only two electrons and are thus present as ions with a He-like electron configuration (e.g. O vii and N vi). These ions produce a triplet of lines consisting of so-called forbidden 2 3S1–2 1S0 (f), inter-combination 2 3P1,2–2 1S0 (i) and resonance 2 1P1–2 1S0(r) lines (Gabriel & Jordan 1969). In low density environments or in the absence of strong UV radiation fields, the upper level of the f transition is only depopulated via the forbidden line, which is then seen as a rather strong line in the spectra (Blumenthal et al. 1972). However, in the wind of a massive star, the presence of an intense UV radiation field alters this situation. Indeed, the UV photons pump the electrons from the upper level of the f transition to the upper level of the i line (Blumenthal et al. 1972; Porquet et al. 2001). Therefore the ratio, is strongly modified, which provides an interesting tool to constrain the dilution of the UV radiation field at the position of the X-ray emission and hence to determine the location of the hot plasma in the wind.

Since the flux of UV radiation seen by the ions in the stellar wind decreases in proportion to the geometrical dilution factor , the ℛ ratio is also a function of radius. The radial modification due to the photoexcitation from the metastable upper level of the forbidden line can be expressed as follows (Blumenthal et al. 1972; Leutenegger et al. 2006): (7)where and with Hν the Eddington flux and f the sum of the oscillator strengths for 2  to all three of 2 . Ne is the free electron density. φc and Nc are the critical photoexcitation rate and the critical density (which depend on atomic parameters and the electron temperature, see Blumenthal et al. 1972). In the winds of hot, luminous massive stars, the photoexcitation rate dominates over the collisional excitation term, and in the remainder of this work, we have thus neglected the Ne/Nc term.

The treatment of these effects is now included in our code along with the above prescriptions and using the Eddington fluxes computed with our CMFGEN model. In Table 2, we indicate the physical information necessary to reproduce this modification of intensities of the forbidden and inter-combination line present in the spectra of ζ Pup.

Table 2

Parameters for the treatment of the forbidden, inter-combination, and resonance lines.

4. Multi-temperature plasma fitting procedure

The full X-ray spectrum of a massive star can usually not be fitted with a single temperature plasma model. In the present case, for instance, the most prominent ionization level of iron (Fe xvii) seen in the RGS spectrum of ζ Pup suggests the presence of a hot plasma component with a kT close to 0.50 keV. Other lines, such as the Si xiii lines around 6.7 Å, have their maximum emissivity at even higher temperatures around 0.90 keV. On the other hand, the fir triplets of He-like carbon and nitrogen are rather indicative of a lower temperature plasma (kT around 0.10 keV). Fitting the observed spectrum hence requires the combination of several models with different temperatures. At this stage, we have to stress that the HETG spectrum of ζ Pup (Cassinelli et al. 2001) as well as its EPIC spectra (Nazé et al. 2012) reveal the presence of a weak Si xiv Ly α line at 6.18 Å and of a weak S xv triplet around 5.04−5.10 Å. These ions have their strongest emission at temperatures near 1.3−1.4 keV. However, these lines are outside the sensitivity range of the RGS instrument and are therefore not included in our fit. We will return to this point in Sect. 5.

Of course, each line is not only emitted at the temperature where its emissivity reaches its maximum value, but rather exists over some range of temperatures, although with a lower intensity. Therefore most lines in the observed spectrum result from the combination of several contributions associated with plasma components at different temperatures (see Fig. 4). The main difficulty in the fitting procedure therefore consists in finding the best values of the parameters that describe our synthetic spectra, including the filling factor fhotgas (i.e. the ratio of the X-ray emitting volume and the total volume) for each plasma component.

We have thus developed a code that combines different components at different temperatures to reproduce the total spectrum of ζ Pup by a minimization of the χ2 of the fit. We based our fitting procedure on a multi-regression routine2. There are quite a number of parameters needed to describe a model. These are the general wind and stellar parameters (, v, β, R, and T), the porosity parameter h′, the abundances of the different elements that make up the hot plasma, and for each plasma component, its temperature (kT), filling factor (fhotgas), and inner and outer boundary (Rin and Rout). Not all of them are treated as free parameters in the present study: only the CNO abundances are varied with respect to solar; v, R and T are all taken from Bouret et al. (2012), whilst we adopted β = 1 to preserve the analytical solution of the optical depth integrals. Nevertheless, we are left with 4 × n + 5 parameters where n is the number of plasma components. In our case, we found that n = 4 is required to fit the observed spectrum and we are thus dealing with 21 free parameters. To obtain the best solution we used our synthetic spectra simulator code (Hervé et al. 2012) to build a database of individual plasma components (hence nine free parameters per element of the database). In doing so, each parameter is sampled over a given range that we consider reasonable (e.g. kT varies between 0.08 and 0.80 keV and varies in the range 1.75 × 10-6 to 4.0 × 10-6M yr-1).

The fitting routine then systematically explores this database, combining only plasma components that have the same chemical composition, h′ and . The synthetic multi-temperature spectrum is built as follows: (8)where Si is the observed flux at the wavelength i, fhotgas,j is the classical volume filling factor of the hot gas component j, Mj,i is the synthetic flux of each X-ray emitting plasma component at the wavelength i, and ϵ is the error of regression.

For each combination of plasma models, the best-fit values of fhotgas,j are determined in the least-squares sense and the resulting is obtained by comparing the model to the fluxed spectrum. These numbers are stored for later use, including the identification of the best-fit model based on the value of . After several tests, it appears that at least four components are needed to achieve a good fit of the fluxed RGS spectrum of ζ Pup.

We included an additional physical constraint in this fitting procedure: we imposed the global X-ray emitting region around the star to be continuous. In other words, from the lowest inner boundary radius of the emitting regions to the largest outer boundary radius, there must not exist any spatial gap where there is no hot plasma emission at all.

It is difficult to make any statement about the uniqueness of the solution. There might indeed exist other combinations of (more than four) plasma components with values of the wind parameters (, h′, etc.) outside the range explored here that could fit the spectrum. However, within the assumptions made here (four discrete plasma components,etc), we are confident that our procedure has identified the best solutions.

5. Results and discussion

To start our analysis, we used the stellar parameters of Bouret et al. (2012). Then we explored a wide range of values for all free parameters to minimize the χ2 of the fit. Our best-fit model of the full RGS spectrum of ζ Pup is illustrated in Fig. 3.

thumbnail Fig. 3

Our best-fit model (model1 of Table 3, in dashed red line) of the observed RGS spectrum of ζ Pup (represented as a solid black line). The most prominent emission lines are labelled. We note that there are many more weaker lines that contribute to the spectrum and are included in our model. The residuals (in the sense observation minus model) are shown in the panels below the spectrum.

Open with DEXTER

5.1. Temperature distribution and hot gas filling factor

The true temperature distribution of the X-ray emitting plasma in the wind of ζ Pup is almost certainly a continuum of temperatures between a minimum and a maximum temperature. Nevertheless, a discretization with four different temperatures is sufficient to correctly reproduce the observed RGS spectrum of ζ Pup. As pointed out above, these four components might not be sufficient to account for the presence of the highly ionized Si xiv and S xv lines seen in the HETG and EPIC spectra. However, although the RGS sensitivity range does not cover the S xv triplet, our model predicts an integrated line flux of 0.64 × 10-13 erg cm-2 s-1, which is not too far away from the value 0.81 × 10-13 erg cm-2 s-1 reported by Cassinelli et al. (2001), especially in view of the relative uncertainties on this flux. The latter uncertainties are likely of the order of 10%, as this triplet contains between 80 and 100 counts in the HETG spectrum, as judged from the Chandra transmission grating catalogue (Huenemoerder et al. 2011)3. Therefore, if a hotter plasma component is needed to fully account for the flux of this triplet, its emission measure is likely very low.

In this respect, it is also interesting to compare our results with the fits of the EPIC spectra reported in Paper I (Nazé et al. 2012). In that paper, we used a four-temperature thermal plasma model with an overlying “slab” of absorbing material. The best-fit plasma temperatures were found to be kT1,2,3,4 = 0.09, 0.27, 0.56, and 2.18 keV. Whilst the first three temperatures are very similar to – or mean values of – what we find here, the fourth one is clearly much higher than found in our analysis. However, the emission measure of the 2.18 keV component in the fits of Paper I is about a factor 100 below that of the 0.56 keV component. Therefore, this very hot plasma should indeed have a very limited impact on the RGS spectrum analysed here. Nevertheless, this remark implies that there exists very hot material in the wind of ζ Pup. It would thus be extremely interesting to collect an HETG spectrum of a quality comparable to the RGS spectrum, and to repeat the kind of analysis presented here, to constrain the location of this very hot plasma.

thumbnail Fig. 4

Contributions of the four plasma components to the best-fit model of the RGS spectrum of ζ Pup (model1, solid black line): kT1 = 0.10 keV (dotted red line), kT2 = 0.20 keV (dashed green line), kT3 = 0.40 keV (dotted-dashed blue line) and kT4 = 0.69 keV (in long dashed-dotted pink line).

Open with DEXTER

Table 3

Best-fit X-ray emitting plasma parameters found in our study.

In previous analyses, the line intensity ratios of hydrogen-like and helium-like ions was used as a diagnostic of the X-ray source temperature (e.g. Waldron & Cassinelli 2007). Thus ratios are expected to provide some average temperature between the various plasma components that contribute to the line formation. However, Fig. 4 indicates that the relative weights of the different plasma components in this average depend upon the element that is considered. For instance, in the case of oxygen, the O viii Lyα line is dominated by the 0.20 keV plasma, whilst the O vii triplet includes contributions from the 0.10 and 0.20 keV plasma components. Waldron & Cassinelli (2007) derived kT = 0.24 keV from the oxygen line ratio, which is consistent with the dominant role of the 0.20 keV plasma in this case. If we consider neon instead, we see that the Ne x and Ne ix lines are mostly formed in the 0.40 and 0.20 keV components, with minor contributions from the 0.69 keV component. For these ions, Waldron & Cassinelli (2007) obtained kT = 0.33 keV from the line ratio. This is consistent with the fact that, in this specific case, the two dominant components have their relative importance reverted in the lines of the H-like and He-like ions. In summary, whilst we confirm that line ratios yield some average estimate of the temperature, it is impossible to know, from a line-by-line analysis, the weights of the different plasma components in this average.

The two plasma components with the lowest temperatures (kT of 0.10 keV and 0.20 keV) extend over the widest range in radius and dominate the emission (model1 in Table 3 and Fig. 4). These emitting regions extend from a region relatively close to the surface of the star (7.5 and 1.5 R for kT of 0.10 keV and 0.20 keV) to very far out in the wind (85 and 38 R, respectively). The two hotter plasma components (kT of 0.40 keV and 0.69 keV) start at 2.7 and 3.1 R and extend over a smaller range of radii, out to 4.0 and 4.1 R, respectively.

With these results, it appears that the region between 3.0 and 4.0 R contains X-ray emitting plasma with three different temperatures. A priori, in an outwards accelerating wind, one expects the velocity jumps due to instabilities to have a lower amplitude in the inner parts of the wind. Therefore, one expects that the energy available in the shocks is only sufficient to heat the shocked matter to rather low temperatures. Farther out in the wind, the difference of velocity between two consecutive clumps of matter increases and the energy available in the shocks increases accordingly. In the outermost parts of the wind, the density and hence the emission measures become low and the same applies to the velocity jumps since the efficiency of radiative acceleration progressively drops to zero.

From this simple reasoning, one would thus expect the hottest plasma components to arise from intermediate radii in the wind, whilst the softer emission would be produced over a wider ranger of radii starting closer to the stellar surface and extending farther away from it.

With the caveat that we currently ignore where the 2.18 keV plasma, responsible for the high-energy extension of the EPIC and HETG spectra, is located, our results are in qualitative agreement with these expectations as well as with the 1D line-driven instability simulation of Feldmeier et al. (1997). These authors showed indeed that shocks can arise very close to the stellar surface (see Fig. 5, Feldmeier, priv. comm.). In our best-fit model, the temperature of the hottest plasma component that is present at a given radial position in the wind increases outwards from the innermost boundary of the emission region (1.5 R, where kT reaches 0.20 keV) to a value of kT = 0.69 keV in a region between 3 and 4 R. Farther out, the temperature of the hottest plasma component decreases again. The maximum shock temperature hence occurs in the 3−4 R range where the velocity jumps predicted by the simulations indeed reach their highest values. At positions closer to the stellar surface or farther out in the wind, the simulations predict velocity differences between two consecutive clumps that are smaller, in agreement with the lower temperatures found in these regions.

Table 4

Stellar wind parameters.

The low value of the filling factor of the kT = 0.69 keV component indicates that only a fraction of the shocks in the region between 3 and 4 R actually heat the plasma to high temperatures. To ensure that most of the kinetic energy of a fast moving clump that hits a slower clump, moving ahead of it, is indeed converted into heat, the collision must occur along the radial direction. From a statistical point of view, in a 3D wind, the probability that the velocity vectors of two colliding clumps are not collinear is most probably higher than the probability that they are collinear. Accordingly only a fraction of the collisions will produce the required heating to reach the hottest temperature.

thumbnail Fig. 5

Snapshot of the radial velocity in the stellar wind as obtained in 1D line-driven instability calculations (Feldmeier, priv. comm.).

Open with DEXTER

The overlap between the highest temperature and the lower temperature is somewhat inherent to our methodology. Indeed, whilst we worked with four different temperatures for the X-ray emitting plasma, we considered that the emission region of each plasma component is continuous (i.e. there are no gaps) between its inner and outer boundary. For instance, our results indicate the presence of a kT = 0.10 keV component relatively close to the star to fit the intensity of the lines at wavelengths above 20 Å, but also in the outer parts of the wind to fit the broad wings of the line profiles. In our code we can only assume a continuous shell that extends all the way from 7.5 to 85 R. The introduction of a discontinuous region would considerably increase the number of free parameters in the regression routine. However, the scenario of an extended, continuous shell for the coolest plasma is supported by two considerations. First, as stated above, non-collinear collisions produce cooler plasma, throughout the whole wind. Second, in a wind-embedded plasma, the hottest plasma cools progressively as it moves outwards with the wind. For instance, adopting the parametrization of the cooling function used by Feldmeier et al. (1997), we estimate a typical cooling time of ~180 s for the 0.40 keV plasma at the inner boundary of that component. In the same way, we estimate a cooling time of ~600 s for the 0.69 keV component. On the other hand, the flow times required for the hot gas to move from the inner boundary of its production site (near 2.7 and 3.1 R for the 0.40 and 0.69 keV plasma, respectively) to a location at 5 R are estimated at about 16 500 s and 13 500 s for the 0.40 and 0.69 keV gas. Hence the X-ray photons emitted at energies of kT = 0.10 keV and 0.20 keV beyond about 5 R exist through a combination of direct emission from weaker shocks and the cooling of the hot plasma initially heated to 0.40 and 0.69 keV.

5.2. Abundances

Our best-fit model of the RGS spectrum of ζ Pup yields an over-abundance of nitrogen along with a depletion of oxygen and carbon (see Table 4). Our CNO abundances agree qualitatively with those of the UV/optical studies of Bouret et al. (2012), but are at odds with the result of Pauldrach et al. (2012), at least as far as the carbon abundance is concerned. The Lyα doublet of C vi at 33.73 Å is the only isolated feature of carbon in the RGS wavelength domain. Nevertheless, it appears impossible to reproduce this line with an over-abundance of carbon as advocated by Pauldrach et al. (2012).

The other elements that are relevant for the X-ray spectrum are assumed to have solar abundances (Andres & Grevesse 1989). We cannot determine the abundances of hydrogen and helium directly with the RGS spectrum since there are no lines of these elements in the X-ray domain. However, they play a key role in the determination of the mass absorption coefficient. Therefore, by default, we used the results of Bouret et al. (2012).

The particular properties of ζ Pup can explain the strongly non-solar CNO abundances. Indeed, ζ Pup is a rapidly rotating runaway star and its kinematic properties and chemical abundances could indicate that the star was part of a binary system that underwent either an episode of mass transfer and supernova explosion of the companion or a merger event (Vanbeveren 2011).

5.3. Mass-loss rate and fragmentation

We find a slightly higher mass-loss rate ((3.5 ± 0.25)    ×    10-6   M yr-1) than the UV/optical and IR studies ((2 ± 0.2)    ×    10-6 and 2.1    ×    10-6M yr-1 respectively, Bouret et al. 2012; Najarro et al. 2011) but the difference is not important. Moreover, Bouret et al. (2012) indicate that their mass-loss rate is probably too low to be consistent with the line driving of the wind. Our determination of the mass-loss rate agrees reasonably well with the value of 2.5    ×    10-6M yr-1 obtained independently by Lamers & Leitherer (1993) from the radio fluxes of ζ Pup and by Oskinova et al. (2007) from a fit of the Hα line and the P v resonance doublet.

Our results indicate that the wind structure of ζ Pup is unlikely to be very porous. Indeed, we have tested various assumptions on the porosity parameter: h′ = 0 (no porosity), h′ = 0.07 (roughly equivalent to the fragmentation frequency n0 = 1.7    ×    10-4 s-1 derived by Oskinova et al. (2006) from a line-by-line investigation of the HETG spectrum) and an intermediate value of h′ = 0.02 (model2 in Table 3). By far the best-fit quality is achieved for the model that has no porosity (see Table 3). A similar conclusion was reached by Cohen et al. (2010), who found that porosity was not required to fit individual lines of the HETG spectrum.

From the stellar and wind parameters assumed (R, v, β) or derived () in this paper, we can estimate the threshold value of the porosity parameter h′ at which clumps become optically thick. From Fig. 1, we find that κ ~ 170 cm2 g-1 at a wavelength of 20 Å and at 4 R. With the volume filling factor of the cool gas (f = 0.05) derived by Bouret et al. (2012), we find that a clump size of l = 0.05   R, i.e. h′ = 0.25 would be needed to achieve τc = 1 at 20 Å (and at 4 R). A value of h′ < 0.02, as suggested by our fits, implies that τc = 1 at 20 Å would only be achieved in the outermost regions of the wind, where the porosity has increased significantly (assuming, as we have done here, h(r) = h′ × r). Therefore, in view of our results, it seems that porosity does not play a role in the formation of the X-ray spectrum of ζ Pup. This situation contrasts with the results of Oskinova et al. (2012) for the WN5 Wolf-Rayet star WR 6. These authors concluded that the wind of the WN5 star had to have a very porous structure to allow the observed X-ray emission to escape.

5.4. The fir triplet formation

In previous studies (Waldron et al. 2001; Leutenegger et al. 2006, 2007; Oskinova et al. 2006), the analysis of the fir triplets of the different He-like ions was made line by line to determine the hot plasma temperature and the dilution of the UV radiation field in the X-ray emission zone. For a single temperature plasma, this analysis constrains, in principle, the position of the X-ray emitting plasma in the wind.

thumbnail Fig. 6

Contribution to the observable emission for the fir components of the N vi He-like triplet. The product of the intrinsic emission per solid angle with exp(−τ) is shown as a function of z (in (p,z) coordinates with p = 0.5 R). The triplet includes contributions from the 0.10 and 0.20 keV plasma components starting at R and r = 1.5 R.

Open with DEXTER

In our analysis, we clearly showed that most spectral lines, including the He-like fir triplets, are the result of the combination of emission from several components at different temperatures (Fig. 4). This implies that the analysis of the fir triplets, assuming they are formed by a single temperature plasma and not taking into account the abundances of the elements, could be biased. We illustrate this point by the O vii and N vi triplets. In Fig. 7, one can see that the line profiles and intensities of these triplets are reasonably well fitted. However, Fig. 4 shows that both triplets are actually formed by contributions of the two lowest temperature plasmas, which produce quite different line profiles and intensities. In the case of the N vi triplet, for instance, the inter-combination line is stronger than the resonance line for the 0.10 keV plasma component (in red in Fig. 4), whilst the reverse situation applies to the 0.20 keV component. The combination of both components fits the observed triplet quite well.

Figure 6 illustrates the contributions to the observed emission of each component of the N vi triplet for a line of sight of impact parameter p = 0.5 R. Owing to the pumping by the UV photons from the upper level of the f transition, this line is strongly surpressed at the benefit of the i line. This is especially true in the innermost regions of the wind, where the 0.20 keV plasma component is the sole contributor to the N vi triplet. As a result, and the ratio between the strength of the inter-combination line and the resonance line tends towards the value of predicted by the atomic data.

thumbnail Fig. 7

Our best fit (dashed red line) of four fir triplets present in the RGS spectrum of ζ Pup (solid black line).

Open with DEXTER

thumbnail Fig. 8

Line formation region of the O vii resonance line. The left panels yield the variation of the intrinsic emission per solid angle (emissivity multiplied by the density squared) along the z axis of the (p,z) coordinates for different values of the impact parameter p. The middle panels provide the optical depth τ along the line of sight to the observer (located at z → ∞). The right panels yield the contribution to the observable emission, i.e. the product of the intrinsic emission per solid angle with exp(−τ). The morphology of this distribution results from the combination of the two plasma components at 0.10 and 0.20 keV that extend between 7.5 and 85 R for the former and 1.5 and 38 R for the latter.

Open with DEXTER

We furthermore note the presence of numerous weak lines in the AtomDB emissivities. Due to the Doppler broadening of the lines and to the limited spectral resolution of the RGS instrument, these weak lines cannot be seen individually but they contribute to the total emission. This is particularly important in the inter-combination line of N vi.

As can be seen in Figs. 3 and 7, our model has some of its strongest deficiencies in the fir triplets. For instance, the model predicts too high a flux in the red wind of the blends of the Mg xi and Si xiii triplets (Fig. 7), which likely indicates an overestimate of the inter-combination (and possibly forbidden) line with respect to the resonance line. Likewise, the resonance lines of the N vi and O vii ions are predicted with slightly too strong a blue-shift compared to the observations. There are several possible reasons for these discrepancies.

Clearly, the sampling of the emitting plasma into four discrete temperature components is a simplification that affects the results. The actual plasma in the wind of ζ Pup is more likely to have a roughly continuous temperature distribution.

Another point is the radial dependence of the hot plasma filling factor. In our model, we assumed that the filling factor of a given plasma component remains constant over the full extent of the emission region of this component. This might be an over-simplification, because the filling factor could change, especially for the 0.10 and 0.20 keV plasma components that span a wide range of radii. Although implementing a radial dependence of the filling factor into our code is in principle possible, we have refrained from doing so, because it would introduce additional free parameters, thereby rendering the search for a best-fit model even more difficult.

A third avenue to explore could be the velocity law of the hot gas. Throughout this paper, we have assumed that the X-ray emitting plasma follows the same β = 1 velocity-law as the cool gas. However, the hot plasma is highly ionized, thus lacking line opacity for the radiative acceleration by the stellar UV radiation. Although Coulomb forces help to carry the hot gas along with the cool wind, one cannot exclude the possibility that the velocity law might differ from the one of the cool wind. Testing the impact of the velocity law on our fits is beyond the scope of the present paper.

Finally, more severe discrepancies, of similar type as those observed here, motivated Leutenegger et al. (2007) to incorporate the effect of resonant scattering into their line-by-line fits of the N vi and O vii triplets4. Using the characteristic Sobolev optical depth τ0, ∗ as a free parameter, Leutenegger et al. (2007) were indeed able to improve the quality, of their fits of these two triplets, but at the expense of very high, sometimes infinite, values of τ0, ∗. Such large characteristic Sobolev optical depths are problematic though, as they require a filling factor of the hot plasma near unity, far away from the values derived here and by other studies (Hillier et al. 1993). This problem was noted by Leutenegger et al. (2007), who advocated that such high filling factors were only required in small regions where line formation takes place. However, such a configuration would then also have a strong impact on the formation of other lines, which is not seen. This is why we decided not to include resonance scattering into our model and do not consider this effect as a good candidate for solving the observed discrepancies.

Using the ℛ ratios evaluated on the HETG spectra of ζ Pup, Cassinelli et al. (2001) estimated a radial range between 1.9 and 4.0 R and between 1.7 and 2.25 R for the formation regions of the Mg xi and Si xiii triplets, respectively. Based on the same dataset, Leutenegger et al. (2006) found an inner emitting radius of the Mg xi and Si xiii triplets of (1.43 ± 0.10)   R, whilst Oskinova et al. (2006) inferred an inner radius of ≃1.5   R. At first sight, these results of Cassinelli et al. (2001) for the Mg xi triplet agree reasonably well with the location of the 0.40 keV plasma component, which dominates the formation of these lines in our model (see Fig. 4). This good agreement between the local analysis of Cassinelli et al. (2001) is due to the narrow radial domain covered by the 0.40 keV component. In contrast, Cassinelli et al. (2001) found that the Ne ix and O vii triplets were likely formed below 4 and 10 stellar radii. At first sight, this result disagrees with the location of the 0.20 keV plasma component, which extends from 1.5 to 38 R. However, the actual location of the formation region of the observable lines results from the interplay between intrinsic emission and absorption, which is rather complex (see Fig. 8). For instance, the observable emission of the O vii resonance line arises mainly from below 15 R, which agrees relatively well with the results of Cassinelli et al. (2001).

For the S xv triplet around 5.04−5.10 Å, Leutenegger et al. (2006) obtained an inner radius of only , i.e. potentially much closer to the photosphere than any of the emitting regions in our model. However, the HETG spectrum analysed by these authors is of rather poor quality at these wavelengths. Unfortunately, the RGS instrument does not cover this wavelength domain and we cannot verify whether or not the morphology of this line is well predicted by our model.

5.5. The Fe XVII λ 17.096 Å line

The forbidden lines of the He-like ions are not the only spectral features arising from metastable levels that are depopulated by the strong UV radiation field. Indeed, the transition Fe xviiλ 17.096 Å also arises from a metastable level. Mauche et al. (2001) have shown that this level can be depopulated by a strong extreme-UV continuum radiation between 190 and 410 Å. They also found that in the case of a stellar surface temperature ~55 kK, the level could be totally depopulated. Since our code does not include a full treatment of the atomic transitions of the Fe xvii ion and we have no observational constraints on the extreme UV emission of ζ Pup, we included an adjustable reduction factor to fit the intensity of this line. Our best-fit results correspond to a division of the nominal line intensity by a factor 4.5 (see Fig. 9). Therefore, it seems that the metastable level is heavily depopulated, which is not surprising because the effective temperature of ζ Pup is 40 kK. The emission seen in Fig. 9 close to 17.1 Å actually corresponds to the Fe xviiλ 17.051 Å line.

5.6. Additional unsolved problems in the fit

In addition to the discrepancies in the fits of the fir triplets discussed above, a few additional regions in the spectra are not well fitted. The main one is the region between 17.2 and 17.65 Å, where the flux is underestimated in our model (see Fig. 9). The atomic database of the AtomDB model does not contain any transitions that could easily explain this strong emission in the observation. We also checked whether there are any radiative recombination continua in that spectral range, but we could not find any. An instrumental effect can be ruled out because a similar bump is seen in the Chandra HETG spectrum of ζ Pup (Cassinelli et al. 2001, their Fig. 1) and Waldron et al. (2001, see Fig.2 in their paper) also found a similar structure in the Chandra HETG spectrum of ζ Ori.

thumbnail Fig. 9

Modification of the theoretical line profile and intensity of the Fe xvii λ 17.096 Å line due to the UV radiative flux. The observations are presented as a black line, the theoretical prediction without UV pumping as a dashed green line, and our best-fit model with an ad hoc reduction of the line due to UV pumping as a dotted red line.

Open with DEXTER

The problem could be caused an underestimate of the O vii series that affects this part of the spectrum. Although the intensity of the fir triplet is fitted correctly, an underestimate of the strength of the O vii edge near 20.8 Å in our model could actually lead to an underestimate of this series. Alternatively, the problem could point at a lack of atomic data of some elements or some peculiar ionization stages that are complicated to calculate. The failure to reproduce the flux level in this region certainly also has an impact on the determination of the attenuation factor of the Fe xvii λ 17.096 Å line as well as on the underestimated continuum beetween 18 and 19 Å. We do not think that this bump could result from a redistribution of the flux of the Fe xvii λ 17.096 Å transition because no such feature is seen in other spectra where pumping is very effective (Mauche et al. 2001).

We encounter a similar problem in the region between 12.3 and 13 Å and for two more “lines” at ~10.0 and 11.0 Å. If there exists another plasma component in the wind of ζ Pup, hotter than those included in our model, some of these features could actually be caused by numerous weak lines of higher ionization species of iron. However, with our current four-component model, these lines are not predicted.

6. Summary and conclusions

In this paper, we have presented a multi-temperature embedded plasma model that fits the RGS spectrum of ζ Pup. Figure 3 shows that the intensity of the strongest lines is well reproduced, although there remain some minor differences. As for any spectral synthesis code, our method relies on some assumptions that introduce some limitations. For instance, one of these limitations is the discretization of the temperature structure of the hot plasma in the wind, because we used “only” four temperatures. Moreover, our temperature grid itself is also discretized with a step in kT of 0.025 keV. Because the line intensities and the ℛ0 ratio for the treatment of the fir triplets are dependent on the temperature, the combination of our assumptions in the entire fitting procedure can probably explain at least some of the remaining small differences between the observed RGS spectrum of ζ Pup and our best-fit model. Additional aspects to be explored in the future are the radial dependence of the hot plasma filling factor and the velocity law of the hot gas. We furthermore assumed a simplified radial dependence of the mass absorption coefficient and used the bridging law of Owocki & Cohen (2006) to calculate the optical depth from the optically thick to optically thin regime. Finally, we stress that even the best-quality atomic data have typical uncertainties of 20%. Some of the differences between our model and the observations could thus result from the uncertainties on the atomic parameters.

Nevertheless, our code provides a major step forward in the analysis of the X-ray spectra of O-type stars. Indeed, we are now in the position to fit the full X-ray spectrum consistently. From these fits, we derived the temperature distribution in the wind and constrained the mass-loss rate as well as the abundances of a number of elements.

To the best of our knowledge, our results provide the first empirical determination of the temperature distribution of the X-ray emitting plasma in the wind of a massive star. In our analysis, the shocks producing X-ray emission arise close to the surface of the star (1.5 R) and extend over a large zone (up to about 85 R). However, the highest temperature (kT = 0.69 keV) is produced in a smaller region, at 3−4 R, in agreement with 1D hydrodynamical calculations.

Another important result of our work is our finding that porosity does not improve the quality of the fit of the global spectrum. It seems thus that the clumps that make up the wind of ζ Pup must have rather small dimensions.

For the He-like fir triplets we showed that analyses that neglect the multi-temperature nature of the emitting plasma could be biased when attempting to determine the inner radius of the emission region of the plasma. This is because several temperatures contribute to the line emissions.

We found a depletion of carbon and oxygen and an over-abundance of nitrogen. We also determined a mass-loss rate of 3.5 × 10-6   M yr-1 for ζ Pup. These results agreed with recent determinations by Bouret et al. (2012) and Najarro et al. (2011), but disagree with the analysis of Pauldrach et al. (2012).

In future studies, we will apply our code to other O-type stars with good-quality high-resolution X-ray spectra to generalize the results obtained for ζ Pup in the present paper.


1

See the section “physics” for a precise description of the calculation of the continuum and lines in a hot plasma on the AtomDB website at the url http://www.atomdb.org/

2

We used the IDL routine IMSL_Multiregress. More information are available at http://idlastro.gsfc.nasa.gov/idl_html_help/IMSL_Multiregress.html

4

Leutenegger et al. (2007) also noted that models without resonant scattering produce an inter-combination line that was not sufficiently blue-shifted compared to the observations (their Figs. 1 and 2). This discrepancy is less severe in our model thanks to the combination of several plasma components with different temperatures and locations.

Acknowledgments

This research is supported by the Communauté Française de Belgique – Action de recherche concertée (ARC) – Académie Wallonie-Europe and by an XMM+INTEGRAL PRODEX contract (Belspo). Y.N. acknowledges support by the FNRS (Belgium). Many thanks to John Hillier for making his code CMFGEN available to the community and to A. Feldmeier for providing some results of his LDI simulations.

References

All Tables

Table 1

Chemical elements and their ionization levels used to calculate the mass absorption coefficient κ and the UV radiation field in CMFGEN.

Table 2

Parameters for the treatment of the forbidden, inter-combination, and resonance lines.

Table 3

Best-fit X-ray emitting plasma parameters found in our study.

Table 4

Stellar wind parameters.

All Figures

thumbnail Fig. 1

Radial dependence of the mass absorption coefficient at different wavelengths (specified in the upper left corner of each panel). The mass absorption coefficients are calculated with the elements listed in Table 1 with abundances found in our study.

Open with DEXTER
In the text
thumbnail Fig. 2

Wavelength dependence of the mass absorption coefficient at different positions inside the wind (1.5, 3, 6, and 15   R as specified in the upper left corner of each panel).

Open with DEXTER
In the text
thumbnail Fig. 3

Our best-fit model (model1 of Table 3, in dashed red line) of the observed RGS spectrum of ζ Pup (represented as a solid black line). The most prominent emission lines are labelled. We note that there are many more weaker lines that contribute to the spectrum and are included in our model. The residuals (in the sense observation minus model) are shown in the panels below the spectrum.

Open with DEXTER
In the text
thumbnail Fig. 4

Contributions of the four plasma components to the best-fit model of the RGS spectrum of ζ Pup (model1, solid black line): kT1 = 0.10 keV (dotted red line), kT2 = 0.20 keV (dashed green line), kT3 = 0.40 keV (dotted-dashed blue line) and kT4 = 0.69 keV (in long dashed-dotted pink line).

Open with DEXTER
In the text
thumbnail Fig. 5

Snapshot of the radial velocity in the stellar wind as obtained in 1D line-driven instability calculations (Feldmeier, priv. comm.).

Open with DEXTER
In the text
thumbnail Fig. 6

Contribution to the observable emission for the fir components of the N vi He-like triplet. The product of the intrinsic emission per solid angle with exp(−τ) is shown as a function of z (in (p,z) coordinates with p = 0.5 R). The triplet includes contributions from the 0.10 and 0.20 keV plasma components starting at R and r = 1.5 R.

Open with DEXTER
In the text
thumbnail Fig. 7

Our best fit (dashed red line) of four fir triplets present in the RGS spectrum of ζ Pup (solid black line).

Open with DEXTER
In the text
thumbnail Fig. 8

Line formation region of the O vii resonance line. The left panels yield the variation of the intrinsic emission per solid angle (emissivity multiplied by the density squared) along the z axis of the (p,z) coordinates for different values of the impact parameter p. The middle panels provide the optical depth τ along the line of sight to the observer (located at z → ∞). The right panels yield the contribution to the observable emission, i.e. the product of the intrinsic emission per solid angle with exp(−τ). The morphology of this distribution results from the combination of the two plasma components at 0.10 and 0.20 keV that extend between 7.5 and 85 R for the former and 1.5 and 38 R for the latter.

Open with DEXTER
In the text
thumbnail Fig. 9

Modification of the theoretical line profile and intensity of the Fe xvii λ 17.096 Å line due to the UV radiative flux. The observations are presented as a black line, the theoretical prediction without UV pumping as a dashed green line, and our best-fit model with an ad hoc reduction of the line due to UV pumping as a dotted red line.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.