Mon. Not. R. Astron. Soc. 000, 000–000 (0000)  Printed 27 February 2009  (MN L^{A}TEX style file v2.2)

A mass-loss rate determination for ζ Puppis from the quantitative analysis of X-ray emission line profiles

David H. Cohen,^{1}^{⋆} Maurice A. Leutenegger,^{2} Emma E. Wollman,^{1} Janos Zsargó,^{3} D. John Hillier,^{3} Richard H. D. Townsend,^{4,5} Stanley P. Owocki,^{4}

1Swarthmore College, Department of Physics and Astronomy, Swarthmore, Pennsylvania 19081, USA 2NASA/Goddard Space Flight Center, Laboratory for High Energy Astrophysics, Code 622, Greenbelt, Maryland 20771, USA 3University of Pittsburgh, Department of Physics and Astronomy, 3941 O’Hara St., Pittsburgh, Pennsylvania 15260, USA 4University of Delaware, Bartol Research Institute, Newark, Delaware 19716, USA 5University of Wisconsin, Department of Astronomy, Madison, 475 N. Charter St., Madison, Wisconsin 53706, USA

27 February 2009

ABSTRACT

We fit every emission line in the high-resolution Chandra grating spectrum of ζ Pup with an empirical line profile model that accounts for Doppler broadening and the effects of attenuation by the bulk wind. For each of sixteen lines or line complexes we determine a best-fitting fiducial optical depth, τ_{∗} ≡ κṀ/4πR_{∗}v_{∞}, and place confidence limits on this parameter. The trend in the optical depth as a function of wavelength is completely consistent with the wavelength dependence of the atomic opacity, κ, indicating that porosity due to large-scale clumping plays no significant role in explaining the only-moderately asymmetric profile shapes observed in ζ Pup. Rather, the modest signatures of attenuation in the observed line profiles are due to the low mass-loss rate, which we determine to be 3.5 ± 0.3 × 10^{−6} M⊙ yr^{−1} from fitting the ensemble of optical depth values derived from each line. The quoted uncertainty is statistical, but the largest source of uncertainty in the derived mass-loss rate is due to the uncertainty in the elemental abundances of ζ Pup, which affects the continuum opacity of the wind, and which we estimate to be a factor of two. Even so, the mass-loss rate we find is significantly below the traditional mass-loss rate determinations for ζ Pup and in line with newer determinations that account for small-scale wind clumping. These results will have important implications for stellar and galactic evolution if ζ Pup is representative of other massive stars.

Key words: stars: early-type – stars: mass-loss – stars: winds, outflows – stars: individual: ζ Pup – X-rays: stars

1 INTRODUCTION

Massive stars can lose a significant fraction of their original mass during their short lifetimes due to their strong, radiation-driven stellar winds. Accurate determinations of these stars’ mass-loss rates are therefore important from an evolutionary point of view, as well as for understanding the radiative driving process itself. Massive star winds are also an important source of energy, momentum, and (chemically enriched) matter deposition into the interstellar medium, making accurate mass-loss rate determinations important from a galactic perspective.

A consensus appeared to be reached by the late 1990s

⋆ E-mail: cohen@astro.swarthmore.edu

that the mass-loss rates of O stars were accurately known observationally and theoretically, using the modified (Pauldrach et al. 1986) CAK (Castor et al. 1975) theory of linedriven stellar winds. This understanding was thought to be good enough that UV observations of spectral signatures of their winds could be used to determine their luminosities with sufficient accuracy to make extragalactic O stars standard candles (Puls et al. 1996).

This consensus has unraveled in the last few years, mostly from the observational side, where a growing appreciation of wind clumping – an effect whose importance has long been recognized (Eversberg et al. 1998; Hillier & Miller 1999; Hamann & Koesterke 1999) – has led to a re-evaluation of mass-loss rate diagnostics, including Hα emission, radio and IR free-free emission, and UV absorption (Bouret et al. 2005; Puls et al. 2006; Fullerton et al. 2006). Accounting for

\begin{matrix}c \\ ◯ \end{matrix} 0000 RAS

2 D. Cohen et al.

small-scale clumping that affects density squared emission diagnostics – and also ionization balance and thus ionic column density diagnostics like UV resonance lines – leads to a downward revision of mass-loss rates by a factor of several, with a fair amount of controversy over the actual factor (Hamann et al. 2008; Puls et al. 2008).

X-ray emission line profile analysis provides a good and independent way to measure the mass-loss rates of O stars. Like the UV absorption line diagnostics, X-ray emission profile diagnostics are sensitive to the wind column density and thus are not directly affected by clumping in the way density-squared diagnostics are. Unlike the UV absorption line diagnostics, however, X-ray profile analysis is not very sensitive to the ionization balance; moreover, as it relies on continuum opacity rather than line opacity, it is not subject to the uncertainty associated with saturated absorption lines that hamper the interpretation of the UV diagnostics.

In this paper, we apply a quantitative line profile analysis to the Chandra grating spectrum of the early O supergiant, ζ Pup, one of the nearest O stars to the Earth and a star that has long been used as a canonical example of an early O star with a strong radiation-driven wind. Previous analysis of the same Chandra data has established that the kinematics of the X-ray emitting plasma, as diagnosed by the line widths, are in good agreement with wind-shock theory, and that there are modest signatures of attenuation of the X-rays by the dominant cold wind component in which the shock-heated X-ray emitting plasma is embedded (Kramer et al. 2003).

The work presented here goes beyond the profile analysis reported in that paper in several respects. We analyze many lines left out of the original study that are weak, but which carry a significant amount of information. We better account for line blends and are more careful to exclude those lines where blending cannot be adequately modelled. We model the continuum emission underlying each line separately from the line itself. We use a realistic model of the spectrometers’ responses and the telescope and detector effective area. And we include the High Energy Grating (HEG) spectral data, where appropriate, to augment the higher signal-to-noise Medium Energy Grating (MEG) data that Kramer et al. (2003) reported on.

Implementing all of these improvements enables us to derive highly reliable values of the fiducial wind optical depth parameter, τ_{∗} ≡ κṀ/4πR_{∗}v_{∞}, for each of sixteen emission lines or line complexes in the Chandra grating spectrum of ζ Pup. Using a model of the wavelength-dependent wind opacity, κ, and values for the star’s radius, R_{∗}, and wind terminal velocity, v_{∞}, derived from UV and optical observations, we can fit a value of the mass-loss rate, Ṁ, to the ensemble of τ_{∗} values, and thereby determine the massloss rate of ζ Pup based on the observed X-ray emission line profiles.

In doing this, we also can verify that the wavelengthdependence of the optical depth values – derived separately for each individual line – is consistent with that of the atomic opacity of the bulk wind, and inconsistent with a gray effective opacity. This provides an independent line of evidence that the source of opacity in the soft X-ray is indeed the atomic opacity due to photoelectric absorption and stands in contrast to the predictions of highly porous wind models (Oskinova et al. 2006; Owocki & Cohen 2006). A porositydominated wind would have an opacity governed by the geometric cross section of optically thick clumps. Porosity’s effects can also be seen in the shapes of individual line profiles (Oskinova et al. 2006; Owocki & Cohen 2006). Our preliminary results indicate that profile models that include porosity are not favored over ones that do not and that very large interclump spacing is required to affect the profile shapes at all (Owocki & Cohen 2006; Cohen et al. 2008). We will extend this result in a forthcoming paper but do not address the effect of porosity on individual line profile shapes directly in the current paper.

The paper is organized as follows: We begin by describing the Chandra data set and defining a sample of well behaved emission lines for our analysis in §2. We briefly evaluate the stellar and wind properties of ζ Pup in §3. In §4 we describe the empirical profile model for X-ray emission lines and report on the fits to the sixteen usable lines and line complexes in the spectrum. We discuss the implications of the profile model fitting results in §5, and summarize our conclusions in §6.

2 THE Chandra GRATING SPECTRUM

All the data we use in this paper were taken on 28-29 March 2000 in a single, 68 ks observation using the Chandra High-Energy Transmission Grating Spectrometer (HETGS) in conjunction with the Advanced CCD Imaging Spectrometer (ACIS) detector in spectroscopy mode (Canizares et al. 2005). This is a photon counting instrument with an extremely low background and high spatial resolution (≈ 1′′). The first-order grating spectra we analyzed have a total of 21,684 counts, the vast majority of which are in emission lines, as can be seen in Fig. 1. We modelled every line in the two spectra, as we describe in §4, and indicate in this figure which of the lines we deemed to be reliable. We only include lines in our analysis that are not so weak or severely blended that interesting parameters of the line-profile model cannot be reliably constrained.

The HETGS assembly has two grating arrays - the Medium Energy Grating (MEG) and the High Energy Grating (HEG) - with full-width half maximum (FWHM) spectral resolutions of 0.0023 Å and 0.0012 Å, respectively. This corresponds to a resolving power of R ≈ 1000, or a velocity of 300 km s^{−1}, at the longer wavelength end of each grating. The wind-broadened X-ray lines of ζ Pup are observed to have v_{fwhm} ≈ 2000 km s^{−1}, and so are very well resolved by Chandra. The wavelength calibration of the HETGS is accurate to 50 km s^{−1} (Marshall et al. 2004).

The two gratings, detector, and telescope assembly have significant response from roughly 2 Å to 30 Å, with typical effective areas of tens of cm^{2}, which are a strong function of wavelength. In practice, the shortest wavelength line with significant flux in the relatively soft X-ray spectra of O stars like ζ Pup is the S xv line complex near 5 Å, and the longest wavelength line is the N vii Lyα line at 24.781 Å. The HEG response is negligible for lines with wavelengths longer than about 16 Å.

The X-ray spectrum of ζ Pup consists of emission lines from H-like and He-like ionization stages of N, O, Ne, Mg, Si, and S, and numerous L-shell lines of iron, primarily Fe xvii. The Lyα lines and often the β and even γ lines of

ζ Pup X-ray line profile mass-loss rate 3

Figure 1. The entire usable portions of the MEG (top) and HEG (bottom) first order (negative and positive orders coadded) spectra of ζ Pup. The binning is native (5 mÅ for the MEG and 2.5 mÅ for the HEG). Vertical dashed lines in the data panels themselves represent the laboratory rest wavelengths of all the lines we could reliably fit with the profile model. Solid (red) vertical lines between the two spectral plots indicate the lines we successfully fit with profile models and lines we attempted to fit but which were too blended to extract meaningful model parameters are indicated by dashed (green) lines. For all blended emission lines we show only one of these solid or dashed lines between the panels, and align it with the bluest line in the blend.

the Lyman series are seen for the H-like ions. There is a weak bremsstrahlung continuum beneath these lines. Overall, the spectrum is consistent with an optically thin, thermal plasma in ionization equilibrium with a range of temperatures from one to several million degrees present. It is possible that there are deviations from equilibrium, although the spectrum is not of high enough quality to show this. There is some evidence from the XMM-Newton RGS spectrum that a few of the emission lines are optically thick (Leutenegger et al. 2007); a possibility we will take into account when discussing the results for those lines.

3 THE STAR

ζ Puppis is a relatively nearby, bright, and well-studied early O supergiant (O4 If) that shows the enhanced nitrogen and deficient carbon and oxygen that is indicative of CNO cycle processed material. Helium is also overabundant (Puls et al. 2006). The star’s rapid rotation may explain the photospheric abundances, though they may instead have resulted from the supernova explosion that is invoked to explain its high space velocity (Vanbeveren et al. 1998). On the other hand, no special mechanism may need to be invoked if the lifetime of mass-loss of ζ Pup has removed enough of the star’s envelope to expose nuclear processed material.

There is some uncertainty regarding the distance to ζ Pup. The spectroscopic parallax (Markova et al. 2004) and trigonometric parallax (Perryman et al. 1997) are in good agreement (∼ 460 pc and 429_{−77}^{+120} pc, respectively). But it has also been suggested that ζ Pup lies farther away, at d ≈ 730 pc, where its space motion and age are consistent with an origin in the Vela Molecular Ridge (Sahu 1992; Sahu & Blaauw 1993). This larger distance implies a larger radius and an Hα mass-loss rate that is larger by a factor of two. On the other hand, the Hipparcos data has recently been reanalyzed and a smaller distance – 332 pc – has been found (van Leeuwen 2007).

We stress that the adopted distance does not affect the X-ray line profile fitting results directly. But it does affect the mass-loss rate we derive from our fits via the dependence of Ṁ on R_{∗} and τ_{∗}, and it affects the fiducial mass-loss rate to which we compare the value we derive from the X-ray profiles in this paper. The Hα and radio mass-loss rates scale as Ṁ ∝ d^{1.5} and the mass-loss rate we derive from the profile fitting results scales as Ṁ ∝ d, so the ratio scales only as the square root of the distance. Thus, any change in the distance will not strongly affect the discrepancy we find between the fiducial mass-loss rate and the one we derive from the Xray line profiles. The radius we use for our mass-loss rate calculation in this paper assumes the spectroscopic parallax distance of 460 pc, which is also assumed for the fiducial Hα mass-loss rate determination.

Detailed spectral synthesis has been carried out from the UV to the IR to determine the stellar and wind properties of ζ Pup, which we list in Table 1. Most of these parameters are taken from Puls et al. (2006). There is a range of wind property determinations in the extensive literature on ζ Pup. The terminal velocity of the wind may be as low as 2200 km s^{−1} (Lamers & Leitherer 1993), and as high as 2485 km s^{−1} (Prinja et al. 1990), though we adopt the determination by the Munich group (Puls et al. 2006), of 2250 km s^{−1}, as our standard.

Mass-loss rate determinations vary as well. This is

Table 1. Stellar and wind parameters adopted from Puls et al. (2006)

parameter  value
Distance  460 pc
Massa  53.9 M⊙
T_{eff}  39000 K
R_{∗}  18.6 R⊙
v_{rot}sinib  230 km s^{−1}
v_{∞}  2250 km s^{−1}
β  0.9
Ṁc  8.3 × 10^{−6} M⊙ yr^{−1}
Ṁd  4.2 × 10^{−6} M⊙ yr^{−1}
Ṁe  1.5 × 10^{−6} M⊙ yr^{−1}

a From Repolust et al. (2004). b From Glebocki et al. (2000). c Unclumped value from Puls et al. (2006). d Also from Puls et al. (2006), but the minimum clumping model, in which the far wind, where the radio emission arises, is unclumped, but the inner wind, where the Hα is produced is clumped. Note that the methodology of Puls et al. (2006) only enables a determination to be made of the relative clumping in different regions of the wind. This value of the mass-loss rate, therefore, represents an upper limit. e From Bouret et al. (2008), from detailed UV spectral modelling, assuming clumping and using a filling-factor approach.

partly because of the uncertainty in the distance to ζ Pup. But, it is also the case that each mass-loss rate diagnostic is subject to uncertainty: density-squared diagnostics like Hα and free-free emission are affected by clumping, no matter the size scale or optical depth of the clumps. Mass-loss rates from UV absorption lines are subject to uncertain ionization corrections. In the last few years there have been attempts to account for clumping when deriving mass-loss rates from both density-squared diagnostics and UV absorption diagnostics. We list several mass-loss rate determinations in the table, with notes about the assumptions behind each determination. The X-ray line profile diagnostics of mass-loss rate that we employ in this paper are not directly affected by clumping; although very large scale porosity (associated with optically thick clumps) can affect the profiles, as we have already discussed.

The star shows periodic variability in various UV wind lines (Howarth et al. 1995) as well as Hα (Berghoefer et al. 1996). Its broad-band X-ray properties are normal for an O star, with L_{x} ≈ 10^{−7}L_{Bol} and a soft spectrum (Hillier et al. 1993), dominated by optically thin thermal line and free-free emission from plasma with a temperature of a few million degrees. The emission measure filling factor of the wind is small, roughly one part in 10^{3}. Weak soft X-ray variability, with an amplitude of 6 percent, and a period of 16.7 hr, was detected with ROSAT (Berghoefer et al. 1996). This low-level variability appears not to affect the Chandra data.

4 EMISSION LINE PROFILE MODEL FITTING

4.1 The Model

The X-ray emission line profile model we fit to each line was first described by Owocki & Cohen (2001), building on work by MacFarlane et al. (1991) and Ignace (2001). It is a simple, spherically symmetric model that assumes that the local emission scales as the ambient density squared and that the many sites of hot, X-ray emitting plasma are smoothly distributed throughout the wind above some onset radius, R_{o}, which is expected to be several tenths of a stellar radius above the photosphere in the line-driven instability scenario (Owocki et al. 1988; Feldmeier et al. 1997; Runacres & Owocki 2002). Attenuation of the emitted X-rays occurs in the bulk, cool (T ≈ T_{eff}) wind component via photoelectric absorption, mainly out of the inner shell of elements N through Si and also out of the L-shell (n = 2) of Fe. Singly ionized helium can also make a contribution at long wavelengths. We assume that the atomic opacity of the cool wind, while a function of wavelength, does not vary significantly with radius. This is confirmed by our non-LTE wind ionization modelling, discussed in §5.1. We further assume a beta-velocity law, v = v_{∞}(1 − R_{∗}/r)^{β}, for both wind components, with β = 1 and v_{∞} = 2250 km s^{−1} as given by UV observations (Puls et al. 2006). The local velocity controls the wavelength dependence of the emissivity, the local optical depth governs the wavelength-dependent attenuation, and the density affects the overall level of emission. The first two of these effects can be visualized in Fig. 2.

We cast the expression for the line profile first in spherical coordinates, but evaluate some of the quantities explicitly in terms of ray coordinates, with the origin at the center of the star and the observer at z = ∞. We integrate the specific intensity along rays of given impact parameter, p, and then integrate over rays. Integrating over the volume of the wind, we have:

L_{λ} = 8π^{2} ∫_{−1}^{+1} dμ ∫_{R_{o}}^{∞} η(μ, r)r^{2}e^{−τ(μ,r)}dr, (1)

where L_{λ} is the luminosity per unit wavelength – it is the X-ray line profile. The angular coordinate μ ≡ cos θ, and η is the wavelength-dependent emissivity that accounts for the Doppler shift of the emitting parcel of wind material (which is completely determined, under the assumptions of spherical symmetry and the velocity law, according to its location, (μ, r)). The emissivity has an additional radial dependence due to the fact that it is proportional to the square of the ambient plasma density. The optical depth, τ, is computed along a ray, z = μr, for each value of the impact parameter, p = \sqrt{1 − μ^{2}}r, as

τ(μ, r) = t(p, z) = ∫_{z}^{∞} κρ(r′)dz′, (2)

where the dummy radial coordinate is given by r′ ≡ \sqrt{z^{′2} + p^{′2}}. The opacity, κ, does not vary significantly across a line (recall it is due to continuum processes – the strong wavelength dependence across a line profile arises purely from the geometry indicated in Fig. 2). Using the continuity equation and the beta-velocity law of the wind, we have:

Figure 2. A visualization of the wind Doppler shift and optical depth – two effects that govern the observed, broadened and asymmetric line shapes. The observer is on the left, and the light solid contours represent the line-of-sight velocity in increments of 0.2v_{∞}, with the blue shifts arising in the left hemisphere and the red shifts in the right. The star is the gray circle at the center, and the inner radius of the wind X-ray emission, R_{o}, is indicated at 1.5 R_{∗} by the solid black circle. The solid heavy contour represents the locus of points with optical depth τ = 0.33, and the dashed and dotted contours represent τ = 1 and 3, respectively. The model parameters visualized here are nearly identical to those of the best-fitting model for the Ne x Lyα line shown in Fig. 8 – R_{o} = 1.5; τ_{∗} = 2.

t(p, z) = τ_{∗} ∫_{z}^{∞} \frac{R_{∗}dz′}{r^{′2}(1 − R_{∗}/r′)^{β}} . (3)

We account for occultation of the back of the wind by the star by setting this optical depth integral to ∞ when p < R_{∗} and z < \sqrt{R_{∗}^{2} − p^{2}}. The constant at the front of eq. 3, τ_{∗} ≡ κṀ/4πR_{∗}v_{∞}, is the fiducial optical depth and is equivalent to the optical depth value along the central ray, integrated down to the stellar surface, in the case where v = v_{∞}. This quantity, τ_{∗}, is the key parameter that describes the X-ray attenuation and governs the shifted and asymmetric form of the line profiles.

We note that the optical depth integral, while generally requiring numerical integration, can be done analytically for integer values of β. We use β = 1 throughout this paper (though we report on tests we did for non-integer β values in §4.3), and for that value of the parameter, the optical depth integral along a ray with impact parameter, p, is given by:

t(p > R_{∗}, z) = \frac{R_{∗}τ_{∗}}{z_{t}} (arctan \frac{R_{∗}}{z_{t}} + \frac{π}{2}

− arctan \frac{R_{∗}μ}{z_{t}} − arctan \frac{z}{z_{t}} ) , (4)

t(p < R_{∗}, z) = \frac{R_{∗}τ_{∗}}{2z_{∗}} log (\frac{R_{∗} − z_{∗}}{R_{∗} + z_{∗}} \frac{R_{∗}μ + z_{∗}}{R_{∗}μ − z_{∗}} \frac{z + z_{∗}}{z − z_{∗}} ) , (5)

where z_{t} ≡ \sqrt{p^{2} − R_{∗}^{2}} and z_{∗} ≡ \sqrt{R_{∗}^{2} − p^{2}}, and the integral has been evaluated at z and ∞.

The intrinsic line profile function we assume for the emissivity at each location is a delta function that picks out the Doppler shift line resonance,

η ∝ δ{λ − λ_{o}[1 − μv(r)/c]}. (6)

This assumption is justified because the intrinsic line width is dominated by thermal broadening, which is very small compared to the Doppler shift caused by the highly supersonic wind flow.

Calculating a line profile model, then, amounts to solving equations 1 and 3 for a given set of parameters, R_{o}, τ_{∗}, the normalization (which determines the overall level of η), and an assumed wind velocity law, described by β and v_{∞}. This last parameter, v_{∞}, influences the emissivity term through its effect on the Doppler shift as a function of radius and spherical polar angle. And for our choice of β = 1, eqs. 4 and 5 replace eq. 3.

The model produces broad emission lines where the overall width (in the sense of the second moment of the profile), for an assumed wind velocity law, is governed primarily by the parameter R_{o}. The closer to the star’s surface R_{o} is, the more emission there is from low-velocity wind material, which contributes to the line profile only near line center. The value of τ_{∗} affects the line’s blue shift and asymmetry. The higher its value, the more blue shifted and asymmetric the profile. Large values of τ_{∗} also reduce the profile width by dramatically attenuating the red-shifted emission component of the line. The interplay of the two parameters can be seen in fig. 2 of Owocki & Cohen (2001).

4.2 Fitting the data

For each line in the spectrum, our goal is to extract values for the two parameters of interest – τ_{∗} and R_{o} – and to place formal confidence limits on these values. We begin the analysis procedure for each line by fitting the weak continuum simultaneously in two regions, one to the blue side of the line and one on the red side (but excluding the wavelength range of the line itself). We assume the continuum is flat over this restricted wavelength region. We then fit the emission line over a wavelength range that is no broader than the line itself (and sometimes even narrower, due to blends with nearby lines, which can induce us to exclude contaminated portions of the line in question). The model we fit to each line is the sum of the empirical line profile model – described by equations 1, 4, and 5 – and the continuum model determined from the fit to the two spectral regions near the line. Note that the inclusion of the continuum does not introduce any new free parameters. The overall model thus has only three free parameters: the fiducial optical depth, τ_{∗}, the minimum radius of X-ray emission, R_{o}, and the normalization of the line. In some cases, where lines are blended, we fit more than one profile model simultaneously, as we describe below, but we generally keep the two main parameters of each profile model tied together, and so the only new free parameter introduced is an additional line normalization.

We fit the wind profile plus continuum model to both the MEG and HEG data (positive and negative first orders) simultaneously, if the HEG data are of good enough quality to warrant their inclusion, and to the MEG data only if they are not. We use the C statistic (Cash 1979) as the goodness-of-fit statistic. This is the maximum likelihood statistic for data with Poisson distributed errors, which these photon-counting X-ray spectra are. Note that the maximum likelihood statistic for Gaussian distributed data is the wellknown χ^{2} statistic, but it is not valid for these data, which have many bins with only a few counts, especially in the diagnostically powerful wings of the profiles.

We determine the best-fitting model by minimization of the C statistic using the fit task in xspec. Once the bestfitting model is found, the uncertainties on each model parameter are assessed using the Δχ^{2} formalism^{1} outlined in chapter 15 of Press et al. (2007), which is also valid for ΔC. We test each parameter one at a time, stepping through a grid of values and, at each step, refit the data while letting the other model parameters be free to vary. The 68 percent confidence limits determined in this manner are what we report as the formal uncertainties in the table of fitting results, below. We also examine the confidence regions in two-dimensional sub-spaces of the whole parameter space in order to look for correlations among the interesting parameters.

We use the relatively strong and unblended Fe xvii line at 15.014 Å to demonstrate this fitting process. We show the MEG and HEG data for this line, along with the best-fitting model (the set of model parameters, τ_{∗}, R_{o}, and normalization that minimizes the C statistic) in Fig. 3. The bestfitting model parameters are: τ_{∗} = 1.97, R_{o} = 1.53 R_{∗}, and a normalization of 5.24 × 10^{−4} photons s^{−1} cm^{−2}. Using the ΔC criterion and testing each of these parameters one at a time (while allowing each of the other parameters to vary), we find that the 68 percent confidence limits on the fit parameters are 1.63 < τ_{∗} < 2.35, 1.38 < R_{o}/R_{∗} < 1.65, and 5.04 × 10^{−4} < norm < 5.51 × 10^{−4}. The confidence limits should be thought of as probabilistic statements about the chance that the true parameter values lies withing the given range, given the physical assumptions of the model.

In Fig. 4 we show 68, 90, and 95 percent confidence limits in two-dimensional τ_{∗}, R_{o} parameter space. We calculate a grid of models (typically 36 by 36), optimizing the other free parameters (just the normalization, in this case) at each point in the grid, and use values of ΔC = 2.30, 4.61, 6.17 (Press et al. 2007) to define the extent of the confidence limits. Plots such as this one are a good means of examining correlations between model parameters, in terms of their abilities to produce similar features in the line profiles. We can see what the trade offs are between parameters in a quantitative way. For example, there is a modest anticorrelation between R_{o} and τ_{∗} evident in the figure. Low values of R_{o} (shock onset close to the photosphere) reduce emission on the line wing relative to the core (more emitting material at low velocity). So although low values of R_{o} (hot plasma as close as 1.15 R_{∗}) are allowed at the 95 percent confidence limit, they require a large wind optical depth, τ_{∗} ≈ 3, to compensate. High τ_{∗} values make lines narrower, as small values of R_{o} do, but they also cause lines to be more blue-shifted and asymmetric. So, there is some degen-

1 This criterion is a specific numerical value of ΔC ≡ C_{i} − C_{min} for model realization i, where C_{min} is the C statistic value for the best-fitting model.

Figure 3. The Fe xvii line at 15.014 Å in the MEG (top) and HEG (bottom), with the best-fitting model superimposed. We have not done any rebinning of the data. The error bars represent Poisson, root-N, statistics. The dashed vertical lines indicate the laboratory rest wavelength of the emission line, and the two dotted vertical lines in each panel indicate the wavelengths associated with the Doppler shift due to the stellar wind terminal velocity of 2250 km s^{−1}. The model is shown as the thick (red) histogram, while the data are shown as (black) solid squares with error bars. The fit residuals are shown in the horizontal windows below the data, with the same one sigma error bars shown with the data.

eracy between these two parameters, but it can be broken for good quality data. We note that the confidence limits listed in the table of model fitting results, which are for individual parameters considered one at a time, will tend to differ somewhat from those inferred from these plots of joint confidence limits.

The value of τ_{∗} expected from the traditional mass-loss rate and a model of the wind opacity at 15 Å is τ_{∗} = 5.30, using the opacity model described in §5.1 (which gives a value of κ = 37 cm^{2} g^{−1} at 15 Å). The best-fitting model with fixed τ_{∗} = 5.30 is shown in Fig. 5. This model does not provide a good fit, having ΔC = 64, implying rejection probabilities well above 99.99 percent. This is the quantitative basis for claims that the X-ray emission lines of O stars in general, and ζ Pup in particular, are too symmetric and unshifted to be explained by the standard wind-shock scenario (Kahn et al. 2001; Cassinelli et al. 2001; Kramer et al. 2003). However, the primary goal of this paper is to quantify the

Figure 4. Confidence contours (68, 90, and 95 percent) for the model fitting of the the Fe xvii line at 15.014 Å. The best fit, shown in Fig. 3, is represented by the filled circle.

Figure 5. The Fe xvii line at 15.014 Å in the MEG (top) and HEG (bottom), with the best-fitting model having τ_{∗} = 5.30 superimposed. This is the value implied by the literature mass-loss rate and our wind opacity model. This fit is statistically unacceptable.

Figure 6. The MEG (top) and HEG (middle) measurements of the Si xiii helium-like complex near 6.7 Å, along with the bestfitting model. This line complex shows a relatively small degree of blue shift and asymmetry, indicative of a low τ_{∗} value, as is expected at short wavelengths, where the wind opacity is smaller. Note that there is a separate set of vertical lines – denoting the rest wavelength and the Doppler shifts associated with the wind terminal velocity – for each of the three components of the line complex (resonance, intercombination, and forbidden lines, from short to long wavelength). In this and the following three figures, we also show the 68, 90, and 95 percent confidence limits in τ_{∗}, R_{o} parameter space (bottom).

Figure 7. The derived value of τ_{∗} is significantly higher for the Mg xii Lyα line at 8.421 Å than for the shorter wavelength Si xiii complex shown in the previous figure.

mass-loss rate reduction compared to the older values from the literature by modelling the wind opacity and the effects of wind attenuation on all the line profiles simultaneously. To enable us to do this, we repeat the fitting procedure described here for the line at 15.014 Å for all of the lines and line complexes in the spectrum.

For the helium-like complexes – O vii, Ne ix, Mg xi, Si xiii, and S xv – we fit a modified version of the wind profile model in xspec that simultaneously fits three separate profiles with the basic parameters (τ_{∗} and R_{o}) tied together and accounts for the altered forbidden-to-intercombination line strength ratios due to the effects of photoexcitation out

Figure 8. The Ne x Lyα line at 12.134 Å shows an intermediate degree of blue shift and asymmetry, indicative of an intermediate τ_{∗} value, as is expected at its wavelength, where the wind opacity is larger than at the wavelength of the Mg xii Lyα line, but not as large as at longer wavelengths. Part of the red wing of this line has been excluded from the fitting because of a possible blend with an iron line.

of the 2^{3}S state, which is the upper level of the forbidden line. This model, which was first described in Leutenegger et al. (2006), assumes a spatial distribution of X-ray emitting plasma, just as the basic wind profile model does, but alters the radius-dependent line ratio according to the ultraviolet mean intensity computed from an assumed model

Figure 9. The O viii Lyα line at 18.969 Å shows a relatively large degree of blue shift and asymmetry, indicative of a higher τ_{∗} value, as is expected at longer wavelengths, where the wind opacity is larger. We did not include the very weak HEG data in the analysis of this line.

atmosphere^{2}. This model thus self-consistently accounts for the effects of the radial dependence of the individual line emissivities on both the line ratio and the profile shapes. Although the components of these complexes are blended, we can extract useful model parameters and confidence limits on those parameters by fitting each complex as a single entity.

We handle other line blends similarly, by fitting two (or more) separate line profile models (with an underlying continuum model, fit separately) with parameters forced to be the same for each component. In some cases, like the Fe xvii lines at 17.051 and 17.096 Å, where the relative intensities of the components are well constrained by atomic physics, we obtain reliable results. In other cases, like the N vii Lyα line at 24.781 Å, which is blended with the N vi Heβ line, it is impossible to accurately model the relative contributions of the two lines, and no reliable information can be obtained from fitting the line blend^{3}. In still other

2 TLUSTY O star model (Lanz & Hubeny 2003) with T = 40000 K and log g interpolated between 3.50 and 3.75. 3 For the nitrogen blend we fit a series of models with two components – one for the Lyα line and one for the Heβ line – trying different values of their relative normalizations, all within a plausible range (of 0.1 to 0.4) as implied by the Astrophysical Plasma Figure 10. Values of τ_{∗} (top) and R_{o} (bottom) derived from the model fits, shown with their 68 percent confidence limits. Line complexes and blends that were fit with multiple model components are represented by only one point.

cases, the blending is mild – through a combination of the second line being weak and the overlap region being small – and we can fit the stronger of the components reliably by simply excluding some of the data. This was the case for the Ne x Lyα line at 12.134 Å, where the extreme red wing is mildly blended with a weak iron line. The line blends that could not be reliably fit are indicated in Fig. 1 by the dashed vertical lines between the panels, and include the helium-like neon complex, which is blended with several iron lines.

After eliminating the complexes too blended to be reliable, we are left with sixteen lines and line complexes that could be fit with the wind profile model as described in the previous subsection and as demonstrated on the Fe xvii line at 15.014 Å, above. The results of these fits are summarized in Table 2. And we show four more representative line fits – spanning a wide range of wavelengths and derived values of τ_{∗} – in Figs. 6, 7, 8, and 9. Note the progression in

Emission Database (APED) (Smith et al. 2001). We found values for the fiducial optical depth ranging from less than 1 to more than 4. The Lyα line is also subject to resonance scattering, the neglect of which will cause an underestimation of the optical depth parameter, τ_{∗} (Leutenegger et al. 2007). Thus, information cannot be reliably obtained from this blended line.

Table 2. Wind profile model fit results

ion  wavelengtha  τ_{∗}  R_{o}  normalizationb
(Å)  (R_{∗})  (10^{−5} ph cm^{−2} s^{−1})
S xv  5.0387, 5.0648, 5.1015 0.01_{−.01}^{+.36} 1.41_{−.11}^{+.15}  2.56_{−.36}^{+.24}
Si xiv  6.1822  0.49_{−.35}^{+.61} 1.46_{−.14}^{+.20}  0.77_{−.14}^{+.11}
Si xiii  6.6479, 6.6866, 6.7403 0.42_{−.13}^{+.14} 1.50_{−.04}^{+.06}  11.2_{−.4}^{+.4}
Mg xi  7.8503  0.65_{−.32}^{+.19} 1.33_{−.13}^{+.12}  1.33_{−.13}^{+.17}
Mg xii  8.4210  1.22_{−.45}^{+.53} 1.34_{−.21}^{+.18}  2.95_{−.24}^{+.24}
Mg xi  9.1687, 9.2297, 9.3143 0.92_{−.16}^{+.19} 1.55_{−.06}^{+.06}  17.8_{−.5}^{+.8}
Ne x  9.7082  0.62_{−.52}^{+1.05} 1.48_{−.19}^{+.27}  0.95_{−.15}^{+.15}
Ne x  10.2388  1.95_{−.87}^{+.28} 1.01_{−.00}^{+.45}  2.99_{−.29}^{+.31}
Ne ix  11.5440  0.83_{−.44}^{+.65} 2.08_{−.36}^{+.54}  5.00_{−.50}^{+.40}
Ne x  12.1339  2.03_{−.28}^{+.24} 1.47_{−.10}^{+.11}  26.9_{−.7}^{+1.1}
Fe xvii  15.014  1.94_{−.33}^{+.32} 1.55_{−.12}^{+.13}  52.4_{−1.6}^{+2.5}
Fe xvii  16.780  2.86_{−.71}^{+.38} 1.01_{−.00}^{+.61}  23.1_{−1.2}^{+1.9}
Fe xviic  17.051, 17.096  2.52_{−.64}^{+.70} 1.47_{−.46}^{+.35}  32.7_{−1.1}^{+0.9}
O viii  18.969  3.02_{−.57}^{+.52} 1.18_{−.17}^{+.41}  37.0_{−2.6}^{+2.8}
N vii  20.9099  4.26_{−1.71}^{+2.28} 1.88_{−.87}^{+.87}  14.8_{−1.9}^{+2.3}
O vii  21.602, 21.804  1.62_{−.79}^{+1.33} 2.53_{−.50}^{+.85}  59.9_{−5.4}^{+4.9}

a Closely spaced doublets in the Lyman series lines and He-like intercombination lines are fit with a single profile model centered at the emissivity-weighted wavelength of the two components. b For the blended lines fit simultaneously, including the He-like complexes, the total normalization of all the lines in the complex is indicated. c We fit these two blended lines simultaneously, with a fixed normalization ratio of 0.9. Both line profile components were forced to have the same τ_{∗} and R_{o} values.

these profiles from fiducial optical depths, τ_{∗}, close to zero at the shortest wavelengths to significantly larger values (up to τ_{∗} = 3) at the longest wavelengths. We summarize the sixteen derived τ_{∗} and R_{o} values, along with their confidence limits, in Fig. 10.

4.3 Sensitivity of fitting results to modelling assumptions

We have made various assumptions and choices in carrying out the line-profile modelling described in the previous subsection. And we therefore have investigated many of these, again using the Fe xvii line at 15.014 Å as a test case. In this subsection, we report on the sensitivity of our results to the following assumptions and choices: background subtraction; determination of the continuum level; exclusion of portions of the line due to possible blending; inclusion of the weak HEG data; the adopted values of β and v_{∞} for the wind; and whether to allow the X-ray volume filling factor to vary with radius (as parameterized by q in f_{X} ∝ r^{−q}, where the filling factor, f_{X}, contributes to the emissivity, η – see Owocki & Cohen (2001)). We will very briefly describe those factors that we found to be unimportant, and discuss in more detail those that did make a difference. The baseline model fitting we describe here is the modelling described in the previous subsection for the 15.014 Å line, except that we fit only the MEG data (so that we may evaluate the effect of including the HEG data).

We examined the default background spectra, which were very weak, and also experimented with fitting the 15.014 Å line with and without the background spectrum subtracted and found almost no difference in the fit quality or fit parameters. We therefore opt to neglect the background when fitting each of the lines in the spectrum. The sensitivity to the continuum fit is a little greater, but still nearly negligible. When we changed the continuum level by a factor of two, none of the parameter values changed by more than ten percent. Some lines in the spectrum are blended with weaker lines. The cleanest way to handle this situation is to exclude the contaminated bins from the modelling. To test the effects of this, we eliminated 0.03 Å from the red wing of the 15.014 Å line and refit the data. We then repeated this experiment eliminating 0.07 Å - leaving only about two-thirds of the data. Even in this second, extreme case, the fit parameters varied by less than ten percent and the confidence regions only expanded slightly.

For most lines, the HEG data is significantly weaker than the MEG data. We find for the 15.014 Å line that including the HEG data changes the best-fitting model parameters by, at most, a few percent, but it does tighten the confidence limits somewhat. The effect of including the HEG data is more significant for the shorter wavelength lines, where the effective area of the HEG is larger relative to the MEG. There is very little penalty for including the HEG data, so we do so for all lines shortward of 16 Å. We also fit the MEG and HEG data separately for the 15.014 line to verify that there are not systematic effects between these two spectra; and there are not. The separate fits give results that are very similar to each other, with significantly overlapping 68 percent confidence limits for all parameters.

The original Owocki & Cohen (2001) line profile model allows for a radially varying filling factor of X-ray emitting plasma, parameterized as a power law function of radius. Values of the power-law index, q, that differ significantly from zero (no radial variation) can cause changes in the line profiles that are not insignificant, effectively weighting the emission from parts of the wind according to their velocity (via the beta-velocity law relationship between velocity and radius). However, we find that when we allow q to be a free parameter the best-fitting value is generally very close to zero. For the representative 15.014 Å line, it is q = −0.09, and q = 0 is included in the 68 percent confidence range. The general result is consistent with that found for this and other stars (Kramer et al. 2003; Cohen et al. 2006). Thus, to keep the number of free parameters manageable, we fix q = 0.

The factors discussed above have a very minor influence on the results of the line fitting. However, the remaining factors can have a significant effect.

The velocity-law exponent, β, affects line profiles for two reasons: (1) the velocity law affects the mapping between radius and Doppler shifted wavelength, and so affects the emission profile; and (2) via the continuity equation, it affects the density and so affects the level of both the emission and the absorption. Indeed, for our representative emission line, when we change the value of β from 1 to 0.8, both τ_{∗} and R_{o} change by 10 to 20 percent. The determinations of β for ζ Pup vary from at least 0.9 to 1.15, and so using a value of β = 1 seems reasonable, especially as it speeds the calculation of the line profile model by allowing the optical depth integral to be done analytically, so we use that value for all the model fitting results reported here. If, in the future, a new and more accurate determination of β is made, and it differs significantly from β = 1, then the results reported in this paper can be scaled accordingly^{4}. We also note that the X-ray emitting plasma and the bulk wind that attenuates the X-rays may not necessarily be described by the same beta velocity law. However, there is no independent evidence for this, and with the short post-shock cooling lengths expected in the relatively dense wind of ζ Pup, the X-ray emitting plasma in the wind is more likely to have a velocity close to the ambient wind velocity^{5}. And furthermore, the observed X-ray emission line widths in ζ Pup and other early O supergiants are completely consistent with the β and v_{∞} values inferred from UV and optical spectroscopy of these stars.

The terminal velocity of ζ Pup is relatively well established, with reasonable estimates from several different groups that vary by about ±10 percent about our adopted

4 Lowering β from 1 to 0.8 causes the best-fitting optical depth of the Fe xvii line at 15.014 Å to go from τ_{∗} = 1.98 to τ_{∗} = 1.66. If the value of β were to be revised upward by a similar amount, the values we derive for τ_{∗} from the line profile fitting would have to be revised upward by about 15 percent. The quality of the fits with the different values of β do not differ significantly. 5 X-ray emitting plasma is too highly ionized to be effectively driven by the photospheric UV radiation field. However, for small enough parcels, the ram pressure of the surrounding wind should keep the post-shock, hot plasma moving at the ambient velocity.

Figure 11. Values of the terminal velocity derived from fitting five strong lines with a wind profile model for which v_{∞} was allowed to be a free parameter (along with τ_{∗}, R_{o}, and the normalization). The bulk wind terminal velocity adopted from the analysis of UV profiles is indicated by the solid horizontal line. The cross-hatched area represents the 68 percent confidence region for the value of the terminal velocity derived from fitting these five points.

value of 2250 km s^{−1}. However, when we explored the effect of varying the terminal velocity in our fitting of wind profile models to the 15.014 Å line, we found that the value of τ_{∗} was quite sensitive to the assumed wind terminal velocity, even within this relatively narrow range. This is because the blue shift of the line centroid in the dimensionless, scaled wavelength parameter, x ≡ (λ/λ_{o}−1)c/v_{∞}, depends directly on the degree of wind absorption. The same observed profile appears more blue shifted in scaled wavelength units if the terminal velocity is (assumed to be) smaller. Our tests with the 15.014 Å line show that the best-fitting value for τ_{∗} ranges from 2.16 to 1.35 when we use terminal velocities between 2200 km s^{−1} and 2485 km s^{−1}. This variation is larger than that caused by every other parameter uncertainty and assumption we have explored. Thus, while we consider the value of v_{∞} = 2250 km s^{−1} to be quite reliable, future reassessments of this parameter will necessitate a rescaling of the optical depth – and mass-loss rate – results we report in this paper.

As a final test, we can treat the terminal velocity as a free parameter of the model. This enables us to see what value of the terminal velocity is preferred by the X-ray spectral data themselves. In general, the constraints on v_{∞}, while letting the other model parameters be free to vary, were not strong. But for the highest signal-to-noise lines in the spectrum, relatively tight constraints could be derived. We show the results for fitting the five most useful lines in Fig. 11. As the figure shows, these lines are all consistent with our adopted value of v_{∞} = 2250 km s^{−1}. This, of course, gives us added confidence that the value we use for the model fitting is reasonable. And, in fact, the small error bars on most of these determinations also show that significantly smaller and larger values are ruled out. The kinematics of the hot, X-ray emitting plasma seem to be the same as that of the bulk wind.

5 DISCUSSION

The most obvious new and significant result of the profile model fitting is the wavelength trend in the derived values of the fiducial optical depth, τ_{∗}, shown in the top panel of Fig. 10. The value of this parameter, which is proportional to both the mass-loss rate and the opacity of the bulk wind, increases with wavelength, which is exactly what is expected from the form of the atomic opacity. The null hypothesis of a constant value of τ_{∗} is rejected with greater than 99.9 percent confidence (χ_{ν}^{2} = 5.4 for 15 degrees of freedom). We therefore fit a model of wavelength-dependent τ_{∗}, in which the wavelength dependence derives entirely from the atomic opacity, κ(λ).

While it may seem obvious that there should be a trend in the fiducial optical depth with wavelength, this result is quite significant, in that a presumed lack of such a trend is the basis for claims that large-scale clumping and the associated wind porosity are the cause of the smaller than expected profile blue shifts and asymmetry (Oskinova et al. 2006). In the following subsections, we show how a realistic wind opacity model naturally explains the observed wavelength trend, and then how such a model can be used to make a quantitative determination of the mass-loss rate of ζ Pup.

5.1 The opacity model and the mass-loss rate determination

The opacity model depends on the abundances and, to a lesser extent, the ionization balance of the bulk stellar wind (i.e. the cooler, unshocked component). The dominant source of opacity is photoelectric absorption from the K-shell of abundant elements between N and Si, and also the L-shell of Fe. We have computed a wind opacity model using cmfgen (Hillier & Miller 1998; Zsargó et al. 2009) with atomic cross sections from Verner & Yakovlev (1995). The model is constrained by UV and optical spectra, so the abundances are derived directly from observations. Details are provided in Bouret et al. (2009). Specifically, it is found that Y_{He} = 0.16 ((Z/Z_{⊙})_{He} = 1.88 expressed as a fraction of the solar abundance), (Z/Z_{⊙})_{C} = 0.03, (Z/Z_{⊙})_{N} = 5.0, (Z/Z_{⊙})_{O} = 0.20, and (Z/Z_{⊙})_{Fe} = 1.0, where the reference solar abundances are taken from Asplund, Grevesse, & Sauval (2005). We show the wind opacity model, using these abundances and the ionization balance from the cmfgen modelling, at a single radius (r = 1.8 R_{∗})^{6} in Fig. 12, along with a solar-abundance model. The opacity is lower at most wavelengths in the cmfgen model primarily because the total abundance of metals (and most crucially the sum of carbon, nitrogen, and oxygen) is subsolar (0.53 of the Asplund, Grevesse, & Sauval (2005) value). We refer to

6 We note that there is very little variation in the opacity with radius between 1.1 R_{∗} and roughly 4 R_{∗} (at least at wavelengths where we analyze lines, below the nitrogen K-shell edge near 26 Å). By 5 R_{∗} the overall opacity is about twenty percent higher, and by 11 R_{∗} it is about a factor of two higher. The increasing opacity with radius is due to the larger fraction of singly ionized helium in the outer wind. But the wind density is so low at these distances that the outer wind does not contribute significantly to the X-ray optical depth.

Figure 12. The wavelength dependent opacity of the wind of ζ Pup computed using cmfgen modelling of the ionization balance and sub-solar abundances derived from UV and optical spectra (solid), along with a solar-abundance opacity model (dotted). Note the prominent K-shell edge of oxygen near 20 Å in the solar abundance model. In the sub-solar metallicity opacity model, this decrement is much more modest, due to the underabundance of O and overabundance of N. The overall reduction in the opacity at most wavelengths in the cmfgen model is the result of its overall subsolar metallicity and not its altered CNO abundances.

this opacity model, based on the cmfgen modelling and the abundances derived from the UV and optical spectra, as the “sub-solar metallicity” model in the remainder of the paper. Using either of these models of the opacity, and values for the stellar radius and wind terminal velocity from Table 1, we can construct a wavelength-dependent model of τ_{∗}, for which the mass-loss rate is the only free parameter. Fits with both the sub-solar metallicity wind opacity model and the solar abundance model are good (χ_{ν}^{2} ≈ 0.6 for the subsolar metallicity model and χ_{ν}^{2} ≈ 0.8 for the solar abundance model), although a higher mass-loss rate of Ṁ = 3.50×10^{−6} M⊙ yr^{−1} is found with the sub-solar metallicity model, due to its lower overall opacity. The solar abundance opacity model, which should provide a lower limiting case, gives Ṁ = 1.90 × 10^{−6} M⊙ yr^{−1}. The formal uncertainties on these derived mass-loss rates, due solely to the finite error bars on the individual τ_{∗} determinations, are about 10 percent.

The best-fitting τ_{∗} model, using the sub-solar metallicity opacities and the best-fitting mass-loss rate, is shown in Fig. 13, along with the τ_{∗} model computed using the standard unclumped mass-loss rate from the literature, Ṁ = 8.3 × 10^{−6} M⊙ yr^{−1}. The best-fitting mass-loss rate is more than a factor of two lower. If solar abundances are assumed for the opacities, the factor is more than three. The bestfitting versions of these two models are compared in Fig. 14, and have a very similar shape, implying that even with better quality Chandra data it would be difficult to distinguish them based on the X-ray data alone. We stress, though, that the abundances of ζ Pup are certainly not solar. We present this model only for comparison with the sub-solar metallicity opacity model, and as a limiting high opacity case.

Taking a closer look at the atomic opacity, we can see in the preceding three figures that the most leverage regarding the wavelength dependence of the opacity, and hence of τ_{∗}, comes at the shortest wavelengths, below the Ne K-shell

Figure 13. Values of τ_{∗} derived from the line-profile model fits, shown as points with error bars (same as the top panel of Fig. 10). The value of τ_{∗} expected from the literature mass-loss rate of 8.3× 10^{−6} M⊙ yr^{−1} is shown as the upper dashed curve. Treating the mass-loss rate as a free parameter, the best fit value of 3.50×10^{−6} M⊙ yr^{−1} is shown as the lower, varying solid curve. This model provides a formally good fit. And both of these models of the wavelength-dependent τ_{∗} use the sub-solar metallicity (cmfgen) opacity model. The horizontal dash-dot line is the best-fitting constant τ_{∗} model, as would be expected for a porosity-dominated wind. It does not provide a good fit to the data.

Figure 14. The best-fitting model, with Ṁ = 3.50 × 10^{−6} M⊙ yr^{−1}, shown in Fig. 13, is shown here again, but this time it is compared to the best-fitting solar abundance τ_{∗} model (dotted curve). The fits are of similar quality, while the solar abundance model has a lower mass-loss rate (Ṁ = 1.90 × 10^{−6} M⊙ yr^{−1}) to compensate for its higher overall opacity.

edges near 13 Å. The Fe and Ne edges and the low O abundance conspire to make the opacity rather flatter than the generally expected κ ∝ λ^{3} relationship seen from individual elements’ photoionization cross sections. Most of the strong lines in the MEG spectra of O stars are between 12 and 18 Å, where the opacity is relatively constant. This demonstrates the need for the use of realistic wind opacity models when interpreting trends in grating spectra of O stars.

Furthermore, the paucity of useful emission lines longward of the O K-shell edge makes it difficult to discriminate among various wind opacity models, although in principle, lines longward of this edge could enable us to diagnose the altered CNO-processed abundances with some certainty. And emission lines longward of the N K-shell edge near 26 Å would be especially useful, but there are none in the Chandra spectrum. The N vii Lyβ line at 20.910 Å is quite weak and does not provide a strong constraint on τ_{∗}, although it does favor the sub-solar metallicity opacity model. The longest wavelength line which we are able to reliably fit is the helium-like O vii complex near 21.8 Å. We fit the resonance and intercombination lines simultaneously (the forbidden line is not present due to 2^{3}S − 2^{3}P photoexcitation by the photospheric UV field), with the profile parameters τ_{∗} and R_{o} tied together for the two lines. However, the resonance line in this complex may be subject to resonance scattering (Leutenegger et al. 2007) – it may be optically thick to its own radiation (as distinct from the effects of continuum opacity of the overlying wind that leads to the observed skewness and blue shifts in all of the line profiles). Resonance scattering tends to make broadened, asymmetric, and blue shifted lines more symmetric, and thus the τ_{∗} value we derive from fitting this complex may be somewhat underestimated. If this is the case, then this line complex too would favor the sub-solar metallicity wind opacity model, as shown in Fig. 14. We also note that the only other line of the sixteen we analyze that is likely to be optically thick to resonance scattering is the O vii Lyα line at 18.969 Å, so the τ_{∗} determination for that line may also be somewhat underestimated.

We also can see from a careful inspection of the opacity model that the mass-loss rate determination from fitting a set of τ_{∗} values is mostly sensitive to the cross section contributions from N, O, and Fe. Alterations of O and N abundances due to CNO processing will have only a modest effect on the results, however. The sum of the contributions of C, N, and O (as well as He) is what affects the overall opacity level between about 15 Å and 20 Å, with Fe – and to a lesser extent, Ne – making a significant contribution at shorter wavelengths. This demonstrates that accurate determinations of abundances for O stars are perhaps the biggest factor in enabling the determination of clumping-independent mass-loss rates from high-resolution X-ray spectra. But when fitting a large ensemble of lines that span a relatively wide range of wavelengths, knowing the overall metallicity is probably sufficient, although including a realistic mixture of elements (and thus absorption edges) is important too.

5.2 Rejection of Gray Opacity and the Porosity Model

The rejection of the hypothesis of a gray effective opacity (wavelength-independent τ_{∗}) argues strongly against the dominance of large-scale clumping and its associated wind porosity in setting the observed properties of O star X-ray spectra. It has been shown that very large-scale clumping, in which individual clumps are optically thick in the soft X-ray continuum, can lead to a reduction in the effective opacity that might explain the only modestly shifted and asymmetric profiles observed in many O star X-ray spectra (Oskinova et al. 2006; Owocki & Cohen 2006). In this scenario, it is the physical cross section of the clumps, and not the atomic cross sections, that governs the wind opacity.

Thus the expectation is that a porosity dominated wind would lead to X-ray emission line profiles with similar morphologies and no wavelength dependence across the observed spectral range. This scenario requires a very large clump scale (Oskinova et al. 2006), however, with porosity lengths in excess of a stellar radius (Owocki & Cohen 2006), where the porosity length represents the interclump mean free path in the limit of very optically thick clumps. Furthermore, preliminary fits of X-ray line profile models that explicitly include porosity as a free parameter do not, in fact, provide better fits to the observed profiles than models that simply have reduced mass-loss rates (Cohen et al. 2008). In an upcoming paper, we will show that this conclusion holds for all strong lines in the Chandra spectrum of ζ Pup, and furthermore, that to reproduce the observed profiles with the standard, unclumped mass-loss rate of Ṁ = 8.3×10^{−6} M⊙ yr^{−1}, requires porosity lengths of more than 2 R_{∗} (which we have already shown is the case for the high signal-to-noise Fe xvii line at 15.014 Å in Cohen et al. (2008)).

5.3 Sources of Uncertainty in the Mass-Loss Rate Determination

The uncertainty in the mass-loss rate determination we have found from the fits to the ensemble of τ_{∗} values, derived from fitting the individual line profiles, come from three sources. The first is the formal uncertainty on the mass-loss rate model that stems from the uncertainties on the individual line profile fits (represented by the error bars on the τ_{∗} points in Fig. 13, for example). For the sub-solar metallicity opacity model the 68 percent confidence limit range on the fitted mass-loss rate extends from 3.25 to 3.73 × 10^{−6} M⊙ yr^{−1}, representing an uncertainty of a little less than 10 percent on the best-fitting value of 3.50 × 10^{−6} M⊙ yr^{−1}.

The second source of uncertainty arises from our imperfect knowledge of the wind terminal velocity (and, most importantly, the terminal velocity of the X-ray emitting plasma itself). However, as we have shown (see Fig. 11), the data themselves indicate that our adopted terminal velocity of v_{∞} = 2250 km s^{−1} is well supported. Three of the lines we show in that figure have best-fitting terminal velocity values near 2350 km s^{−1}, which is also the terminal velocity derived from a careful analysis of the UV line profiles by Haser (1995). When we refit the representative Fe xvii line at 15.014 Å using this higher terminal velocity, we found a reduction in our derived τ_{∗} value of 15 percent. If this scaling holds for all lines, then using this slightly higher value of the terminal velocity will lead to a downward revision of our derived mass-loss rate of about 15 percent. (Note that the terminal velocity enters into the denominator of the expression for τ_{∗}, and that will mitigate this adjustment slightly.) Similar considerations pertain to our assumption about the wind velocity parameter, β.

Finally, we estimate that the abundances derived for ζ Pup from the analysis of UV and optical data have a precision of only about a factor of two (Bouret et al. 2009). We can see from the comparison of the sub-solar metallicity model to the solar abundance model that the mass-loss rate varies by about a factor of two between these two assumed opacity models, although the solar abundance model is included in our analysis not so much as a realistic alternate model, but simply as a plausible upper bound to the atomic opacity; the sub-solar metallicity model is more realistic due to the constraints on it provided by observations in other wavelength bands, and of course, the evolved nature of ζ Pup implies that we should not expect to find solar abundances in its wind. However, the overall CNO metallicity of 0.53 solar is lower than expected, given the uniformly solar abundances in nearby massive stars (Przybilla et al. 2008). Thus, a conservative estimate of the allowed range of the mass-loss rate of ζ Pup derived from the X-ray line profile fitting is roughly 2 to 4 × 10^{−6} M⊙ yr^{−1}, with our best estimate being 3.50 × 10^{−6} M⊙ yr^{−1}.

This mass-loss rate is only a little lower than the maximum mass-loss rate of 4.2 × 10^{−6} M⊙ yr^{−1} (Puls et al. 2006), implying a small amount of clumping in the outer wind, and a small adjustment to the clumping factor in the inner wind determined by Puls et al. (2006). However, we should point out that the recent mass-loss rate determination based on a clumping analysis of the UV spectrum (Bouret et al. 2008, 2009) of only Ṁ = 1.5 × 10^{−6} M⊙ yr^{−1} could be marginally consistent with our results if the overall metallicity of ζ Pup is close to solar. This would increase the opacity in our wind model and lower the mass-loss rate by a corresponding amount, even if the individual abundances are affected by CNO processing.

5.4 Location of the X-ray Emitting Plasma

The analysis of the sixteen lines and line complexes in the Chandra spectrum of ζ Pup also enables us to derive values of the onset radius of the wind-shock X-ray emission from the profiles. These results are shown in the lower panel of Fig. 10, and are completely consistent with the expectations of the wind-shock structure induced by the line-driven instability (Feldmeier et al. 1997; Runacres & Owocki 2002). That is, an onset radius of R_{o} ≈ 1.5 R_{∗}. We have searched for a trend with wavelength in these values and found none (an unweighted fit of a linear trend shows a modest increase with wavelength, but that result is significant at only the one sigma level, and when we perform a weighted fit – with the weights inversely proportional to the uncertainties on the individual measurements – the significance is less than one sigma). Thus, the simplest interpretation is that there is a universal radius of the onset of X-ray emission and it occurs near 1.5 R_{∗} (half a stellar radius above the photosphere). This result had already been noted by Kramer et al. (2003), though we show it more robustly here. This same result can also be seen in the late O supergiant ζ Ori (Cohen et al. 2006). And this result is also consistent with the joint analysis of X-ray line profile shapes and helium-like forbidden-to-intercombination line ratios for four O stars as described by Leutenegger et al. (2006).

5.5 Comparison with Previous Analyses

Finally, let us consider why we have found a trend in wavelength for the fiducial optical depth values, τ_{∗}, derived from the same Chandra data that led Kramer et al. (2003) to report that there was no obvious trend. The two biggest factors that have led to this new result are our more careful assessment of line blends and our inclusion of many weak, but important, lines at short wavelength. Kramer et al. (2003) included only one line shortward of the Ne x Lyα line at 12.134 Å, whereas we report on nine lines or line complexes in this range (including two helium-like complexes, which Kramer et al. (2003) excluded from their analysis). While many of these lines are weak and do not provide very strong constraints when considered individually, taken together, they do provide strong constraints. As far as line blends are concerned, Kramer et al. (2003) included the N vii Lyα line at 24.781 Å and the Fe xvii complex near 15.26 Å, both of which we have determined are too blended to enable the extraction of reliable information about their intrinsic profile shapes. Furthermore, we properly account for the blended Fe xvii lines at 17.051 and 17.096 Å, fitting them simultaneously, while Kramer et al. (2003) fit them as a single line. Our use of detailed models of the effective area and resolution also may play a small role in making our results more reliable than those in the initial paper.

An additional factor of some importance is our use of a detailed model of the wind opacity. The detailed wind opacity is relatively flat over much of the wavelength range where the strong lines in the Chandra spectrum are. Specifically, from about 12 Å to about 18 Å, the presence of successive ionization edges makes the overall opacity roughly flat. Thus, for a trend to be apparent, short wavelength lines have to be included in the analysis, as we have already pointed out. Finally, the mass-loss rate reduction we have determined here is only a little more than a factor of two, while earlier analyses suggested that, without porosity, much larger mass-loss rate reductions would be required to explain the only modestly shifted and asymmetric profiles (Kramer et al. 2003; Oskinova et al. 2006). Here too, the wind opacity is key. The overall opacity of the wind is significantly lower than had been previously assumed, implying that the mass-loss rate reduction is not as great than had been previously assumed. Again, this is primarily due to the significantly sub-solar abundances (especially of oxygen) in ζ Pup.

6 CONCLUSIONS

By quantitatively analyzing all the X-ray line profiles in the Chandra spectrum, we have determined a mass-loss rate of 3.5×10^{−6} M⊙ yr^{−1}, with a confidence range of 2 to 4×10^{−6} M⊙ yr^{−1}. The largest uncertainty arises from the abundances used in the atomic opacity model. This method of mass-loss rate determination from X-ray profiles is a potentially powerful tool for addressing the important issue of the actual mass-loss rates of O stars. Care must be taken in the profile analysis, however, as well as in the interpretation of the trends found in the derived τ_{∗} values. It is especially important to use a realistic model of the wind opacity. And for O stars with weaker winds, especially, it will be important to verify that the X-ray profiles are consistent with the overall paradigm of embedded wind shocks. Here, an independent determination of the terminal velocity of the X-ray emitting plasma by analyzing the widths and profiles of the observed X-ray lines themselves will be crucial. In the case of ζ Pup, we have shown that the X-ray profiles are in fact consistent with the same wind kinematics seen in UV absorption line spectra of the bulk wind. And the profile analysis also strongly constrains the onset radius of X-ray production to be about r = 1.5 R_{∗}.

A major conclusion from the profile analysis is that there is no need to invoke large scale porosity to explain individual line profiles, as the overall wavelength trend argues strongly against porosity being the dominant cause of the reduced effective opacity of the wind. Rather, the opacity is completely consistent with the expected wavelengthdependent atomic opacity, with the lower-than-expected wind optical depths being due, simply, to a reduction in the wind mass-loss rates. This reduction (a little more than a factor of two) is consistent with other recent determinations that correct the traditional wind mass-loss rate diagnostics for small-scale clumping that affects density-squared diagnostics and ionization corrections.

ACKNOWLEDGMENTS

Support for this work was provided by the National Aeronautics and Space Administration through Chandra award number AR7-8002X to Swarthmore College and award number TM6-7003X to the University of Pittsburgh, issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics and Space Administration under contract NAS8-03060. EEW was supported by a Lotte Lazarsfeld Bailyn Summer Research Fellowship from the Provost’s Office at Swarthmore College. MAL also acknowledges support from the Provost’s office of Swarthmore College. RHDT, SPO, and DHC acknowledge support from NASA LTSA grant NNG05GC36G and JZ and DJH acknowledge support from STScI grant HST-AR-10693.02. The authors also wish to thank Marc Gagné, Alex Fullerton, Yaël Nazé, and Joachim Puls for several enlightening discussions.

REFERENCES

Asplund M., Grevesse N., Sauval A. J., 2005, in Bash F. N., Barnes T. G., eds, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis. Astron. Soc. Pac., San Francisco, p. 25 Berghoefer T. W., Baade D., Schmitt J. H. M. M., Kudritzki R.–P., Hillier D. J., Pauldrach A. W. A., 1996, A&A, 306, 899 Bouret J.–C., Hillier D. J., Lanz T., 2008, in Bresolin F., Crowther P., Puls J., eds, I.A.U. Symp. 250, Massive Stars as Cosmic Engines. Cambridge University Press, Cambridge, p. 526 Bouret J.–C., Lanz T., Hillier D. J., 2005, A&A, 438, 301 Bouret J.–C., Lanz T., Hillier D. J., Foellmo C., 2009, in preparation Canizares C. R., et al., 2005, PASP, 117, 1144 Cash W., 1979, ApJ, 228, 939 Cassinelli J. P., Miller N. A., Waldron W. L., MacFarlane J. J., Cohen D. H., 2001, ApJ, 554, L55 Castor J. I., Abbott D. C., Klein R. I., 1975, ApJ, 195, 157 Cohen D. H., Leutenegger M. A., Grizzard K. T., Reed C. L., Kramer R. H., Owocki S. P., 2006, MNRAS, 368, 1905 Cohen D. H., Leutenegger M. A., Townsend R. H. D., 2008, in Hamann W.–R., Feldmeier A., Oskinova L. M., eds, Clumping in Hot Star Winds. Universitatsverlag, Potsdam, p. 209 Dessart L., Owocki S. P., 2003, A&A, 406, L1 Eversberg T., Lepine S., Moffat A. F. J., 1998, ApJ, 494, 799 Feldmeier A., Puls J., Pauldrach A. W. A., 1997, A&A, 322, 878 Fullerton A. W., Massa D. L., Prinja R. K., 2006, ApJ, 637, 1025 Glebocki R., Gnacinski P., Stawikowski A., 2000, Acta Astron., 50, 509 Hamann W.–R., Feldmeier A., Oskinova L. M., 2008, Clumping in Hot Star Winds. Universitatsverlag, Potsdam Hamann W.–R., Koesterke L., 1999, A&A, 335, 1003 Haser S. M., 1995, Universitäts-Sternwarte der Ludwig-Maximillian Universität, München Hillier D. J., Miller D. L., 1998, ApJ, 496, 407 Hillier D. J., Miller D. L., 1999, ApJ, 519, 354 Hillier D. J., Kudritzki R.–P., Pauldrach A. W. A., Baade D., Cassinelli J. P., Puls J., Schmitt J. H. M. M., 1993, ApJ, 276, 117 Howarth I. D., Prinja R. K., Massa D., 1995, ApJ, 452, L65 Ignace R., 2001, ApJ, 549, 119 Kahn S. M., Leutenegger M. A., Cottam J., Rauw G., Vreux J.–M., den Boggende A. J. F., Mewe R., Güdel M., 2001, A&A, 365, L312 Kramer R. H., Cohen D. H., Owocki S. P., 2003, ApJ, 592, 532 Lamers H. J. G. L. M., Leitherer C., 1993, ApJ, 412, 771 Lanz T., Hubeny I., 2003, ApJ, 146, 417 Leutenegger M. A., Paerels F. B. S., Kahn S. M., Cohen D. H., 2006, ApJ, 650, 1096 Leutenegger M. A., Owocki S. P., Paerels F. B. S., Kahn S. M., 2007, ApJ, 659, 642 MacFarlane J. J., Cassinelli J. P., Welsh B. Y., Vedder P. W., Vallerga J. V., Waldron W. L., 1991, ApJ, 380, 564 Markova N., Puls J., Repolust T., Markov H., 2004, A&A, 413, 693 Marshall H. L., Dewey D., Ishibashi K., 2004, SPIE, 5165, 457 Oskinova L., Feldmeier A., Hamann W.–R., 2006, MNRAS, 372, 313 Owocki S. P., Castor J. I., Rybicki G. B., 1988, ApJ, 335, 914 Owocki S. P., Cohen D. H., 2001, ApJ, 559, 1108 Owocki S. P., Cohen D. H., 2006, ApJ, 648, 565 Pauldrach A. W. A., Puls J., Kudritzki R.–P., 1986, A&A, 164, 86 Perryman M.A.C., et al., 1997, A&A, 323, L49 Press W. H., Flannery B. P., Teukolsky S. A., Vetterling W. T., 2007, Numerical Recipes, 3^{rd} edition. Cambridge University Press, Cambridge Prinja R. K., Barlow M. J., Howarth I. D., 1990, ApJ, 361, 607 Przybilla N., Nieva M.-F., Butler K., 2008, ApJ, 688, L103 Puls J., Kudritzki R.–P., Herrero A., Pauldrach A. W. A., Haser S. M., Lennon D. J., Gabler R., Voels S. A., Vilchez J. M., Wachter S., Feldmeier A., 1996, A&A, 305, 171 Puls J., Markova N., Scuderi S., Stanghellini C., Taranova O. G., Burnley A. W., Howarth I. D., 2006, A&A, 454, 625 Puls J., Vink J. S., Najarro F., 2008, A&A Rev., 16, 209 Repolust T., Puls J., Herrero A., 2004, A&A, 415, 349 Runacres M. C., Owocki S.P., 2002, A&A, 381, 1015 Sahu M. S., 1992, University of Groningen, Groningen Sahu M., Blaauw A., 1993, in Cassinelli J. P., Churchwell E. B., eds, Massive Stars: Their Lives in the Interstellar Medium. Astron. Soc. Pac., San Francisco, p. 278 Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJ, 556, L91 van Leeuwen F., 2007, Hipparcos, the New Reduction of the Raw Data. Springer, Berlin Vanbeveren D., de Loore C., van Rensbergen V., 1998, A&A Rev., 9, 63 Verner D. A., Yakovlev D. G., 1995, A&AS, 109, 125 Zsargó J., et al., 2009, in preparation