A Deep Exposure in High Resolution X-Rays Reveals the Hottest Plasma in the ζ Puppis Wind

ζ Puppis is one of the optically brightest O-stars, and is close at an estimated distance of only 332 pc (Howarth & van Leeuwen 2019). As such, this O4If(n)p star (Sota et al. 2011) has been the target of many important studies to characterize its variability and properties. Much of this work was summarized in a recent analysis of optical spectroscopy and high-precision photometry by Ramiaramanantsoa et al. (2018), to whom we refer readers for relevant background discussion about the target star.

ζ Pup's relatively high X-ray flux of about ${10}^{-11}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ (1–$40\,\mathring{\rm A}$) allows detailed analysis of its spectrum when dispersed by the High Energy Transmission Grating Spectrometer (HETGS) on board the Chandra X-ray Observatory. An existing $68\,\mathrm{ks}$ Chandra/HETGS spectrum of ζ Pup was analyzed by Cassinelli et al. (2001) who demonstrated that the spectrum is well described by multitemperature thermal plasma models. To take advantage of the HETGS high resolution by driving down the statistical errors, we have obtained an additional 813 ks of exposure time for this star.

Below about 9 Å, the resolution of the HETGS and the low emission line density allows a relatively clean separation between lines and continuum. The thermal Bremsstrahlung (free–free) continuum shape is sensitive to the hottest plasmas expected (∼10 MK). In addition, the emission lines in this region are from highly ionized species (H- and He-like ions) of relatively abundant elements (Mg, Si, S, Ar, Ca, and Fe), and are also formed at high temperatures.

The XMM-Newton spacecraft has been used to study this star, starting with Kahn et al. (2001) and continuing with the extensive series of papers (e.g., Nazé et al. 2012, 2018). However, the sensitivity of the RGS instrument on board the XMM-Newton spacecraft is much lower below 9 Å and its lower resolution impedes the clear separation of the continuum and emission lines. Hence, HETG is the best current instrument for study of both the line and continuum emission in the short wavelength region we concentrate on here.

In this paper, our primary interest is in determination of plasma temperatures, in particular the characterization and measurement of high-temperature lines and continuum in the Chandra/HETGS spectrum of ζ Pup. If X-ray emission originates solely from embedded-wind shocks (Lucy & White 1980; Feldmeier et al. 1995), then the maximum temperature is determined by the highest relative velocities of the colliding structures, and the amount of emission is indicative of the emitting volume. Hot thermal plasmas emit strongly in bound–bound lines of highly ionized states, in thermal Bremsstrahlung radiation, and in bound-free continuum emission. Bremsstrahlung and bound-free emission drop exponentially for photon energies larger than about kT of the plasma (in which k is Boltzmann's constant and T is the plasma temperature). For some plasma temperatures, bound-free emission can exceed that from thermal Bremsstrahlung at high energies (Landi 2007; Kaastra et al. 2008). However, the continuum shape is similar to that from Bremsstrahlung. The drop in the spectrum at high energies (short wavelengths) is sensitive to the hottest temperatures present in a multithermal plasma. Thermal emission lines are also very sensitive to temperature, having temperatures of peak emissivities that increase with atomic number; they typically emit most strongly over a range in temperature of about 0.3 dex (defined by the full width at half maximum of their emissivity curves). Hence, we concentrate here on the 2–9 Å region due to the relative sparseness of emission lines, the high-temperature sensitivity, and the reduced importance of wind absorption. The temperature sensitivity occurs because for the expected plasma temperatures from embedded-wind shocks (2–10 MK), the spectrum will drop accordingly somewhere in the 1–10 Å range. In this regime, the wind absorption is relatively unimportant since the continuum opacity (due mainly to K-shell absorption by metals) drops very steeply to short wavelengths, so we do not have to be concerned with the wind density structure to interpret the spectral shape.

This new data set allows us to probe the highest temperature plasma with unprecedented precision, providing a basic test for any models of X-ray production. Our aim is to use simple but robust methods to infer this maximum temperature.

The observations analyzed here are part of a Cycle 19 Large Project (PI Waldron, Sequence Number 201176), to which we also added the first HETGS spectrum from 2000 (PI Cassinelli, Observation ID 640). Chandra Observations were made with the HETG in conjunction with the ACIS-S detector array, a configuration known as HETGS, the High Energy Transmission Grating Spectrometer (Canizares et al. 2005). An observing log is given in Table 1. We note that secular changes in the HETGS sensitivity have largely occurred at longer wavelengths than those we are concerned with here (less than 10% change below $5\,\mathring{\rm A}$), making it relatively easy to seamlessly include the ObsID 640 data with the more recent data.¹² Any remaining small changes in the instrument's properties are dealt with through calibration, as described below.

Spectra were extracted using standard CIAO procedures (version 4.11; Fruscione et al. 2006), using calibration database version 4.8.2. For each of the 22 observations, we extracted the positive and negative first diffraction-order count spectra for the two types of gratings (High and Medium Energy Gratings, or HEG and MEG). For each of these count spectra, responses were computed from CIAO programs—the effective area files, or "ARFs," and the response matrices, or "RMFs." The latter encodes the energy-dependent line-spread function and enclosed count fraction for the extraction aperture used.

The variability in the 2–9 Å spectral region is small enough to justify our analysis of the overall mean spectrum (see Section 4). Detailed analysis of the X-ray variability of this star will be described in a separate paper (J. S. Nichols et al. 2020, in preparation).

While background is usually negligible for HETGS¹³ due to the ability to do spectral order-sorting through rejection of events using the intrinsic energy resolution of the detector, background does become appreciable in our analysis in the 1.7–3.0 Å region where the source flux and the effective area both drop (the effective area at 3 Å is about $75\,{\mathrm{cm}}^{2}$ but is only $18\,{\mathrm{cm}}^{2}$ at $1.8\,\mathring{\rm A}$). The interesting high-temperature lines of Fe xxv, $\lambda 1.85\,\mathring{\rm A}$ (maximum emissivity at $\sim 60\,\mathrm{MK}$) occur in this region, along with the continuum from the highest temperature plasma. Hence, we have included background spectra, taken from regions adjacent to the dispersed spectra, in our analysis. This is always done by adding background-count to model-count spectra during fitting; the only time we subtract background counts from the data is for the purposes of visualization.

Fitting and modeling were done using the Interactive Spectral Interpretation System (ISIS;¹⁴ Houck & Denicola 2000) which provides interfaces to an atomic database ("AtomDB;" Foster et al. 2012) for atomic data and emissivities used to construct plasma models (including bound-free, free–free, two-photon, and dielectronic recombination processes); it provides access to Xspec models, from which we used "phabs" for absorption, which is based on cross-sections from Verner et al. (1996). This absorption model is adequate for the spectral region of interest since the major sources of opacity are inner-shell electrons of high Z ions. The ionization state and nonsolar abundances of C, N, and O do not matter below $\sim 10\,\mathring{\rm A}$. ISIS also has features that facilitate simultaneous modeling of the 88 spectra (22 observations with 4 orders each) and their associated responses.

The fundamental paradigm in X-ray spectral analysis is iterative forward folding, in which the instrumental response is applied to a model spectrum via an integral equation in order to generate model counts. A background (empirical or model) is optionally added to the model counts. The model parameters are adjusted (for fitting done here, usually with an amoeba-subplex optimization method) to minimize the residuals between the observed and model counts, using Poisson counting statistics. This is necessary because while the response matrix is nearly diagonal for diffraction grating instruments, it still has significant off-diagonal terms (representing the spectral resolution and instrumental profile), and cannot be inverted uniquely. The only robust method is to maintain each count spectrum and their individual responses.¹⁵ We do use an approximate "unfolding" for visualization in division of the counts by the model count rate per unit flux. This does not deconvolve line profiles, nor is it guaranteed to be accurate in regions of rapidly changing spectral response with wavelength. But it is good for visualization of the data in model space and for assessing the overall quality of the fit (residuals are still formed between modeled and observed counts). To improve efficiency in fitting, we combine data sets across observations and orders for each grating type; that is, we combine HEG positive and negative orders, and likewise for MEG. This means that models only need to be evaluated twice, once for each grid. Internally to the system, however, each data set and response are still unique: model counts are evaluated for each array, then counts are summed, models are summed, and residuals computed. Hence, we can work with any subset of the 22 observations without creating actual files (counts and responses) for each permutation of interest.

We proceeded in iterative steps in modeling the spectral region below $9\,\mathring{\rm A}$. Since emission lines are present that are formed at very different temperatures (e.g., Mg xi and S xvi), the plasma is not isothermal, so we began by fitting a two-temperature-component model to provide overall plasma characteristics and to guide further analysis. For our more definitive multitemperature thermal plasmas, each temperature component is constrained to have the same line shape, elemental abundances, and multiplicative absorption function. Because for the purposes of this paper we are primarily interested only in the line strengths (as a probe of temperature, abundance, and absorption), a Doppler shifted and broadened Gaussian is used to model the emission line shapes. A more detailed analysis of the lines would include distributed wind absorption and emission in an expanding medium (see, e.g., Owocki & Cohen 2001; Oskinova et al. 2006; Zhekov & Palla 2007), but as can be seen from the fits displayed later, this simpler model is adequate for our purposes of modeling the overall line fluxes and separating the lines from the continuum. Whatever the actual physical effects determining the line shapes, it is fortunate that the lines seem simply shifted rather than strongly asymmetric (see the fits in Cassinelli et al. 2001; Waldron & Cassinelli 2001). Likewise, the absorption is not really a foreground screen of absorption; our assumption is similar to the "exospheric" approximation (Owocki & Cohen 1999), which assumes all emission is from above the wind optical depth for X-ray emission, ${\tau }_{\lambda }=1$. This is reasonable for the shorter wavelengths, which can be seen from geometrically deep in the wind, as far down as the photosphere (see Figure 5 in Oskinova et al. 2006). We add a "foreground" screen of absorption to approximate the distributed wind absorption. This should be reasonable in understanding the shape of the spectral energy distribution. In the more detailed analysis of broadband X-ray absorption of Leutenegger et al. (2010), it was found that below $10\,\mathring{\rm A}$, the foreground absorbing screen approximation was adequate.

The free parameters per plasma temperature component are then the temperature and normalization (emission measure). The line width and shift were determined independently by fitting line-dense regions (∼9–13 Å) with an AtomDB model template, and then frozen. We used a uniform Doppler shift of −450 km s⁻¹ and a Gaussian broadening term (in addition to thermal broadening) with a full width at half maximum of about 1900  km s⁻¹, which is reasonable given the terminal wind velocity of 2250 km s⁻¹ (Puls et al. 2006). The apparent Doppler shift is a signature of wind absorption, due to line-of-sight absorption of the receding wind's X-ray emission. In general, these parameters change with wavelength, due to the region of formation in the wind. Using the mid-spectrum gives us an average value, and later during fitting we bin the spectrum substantially so that we are not so sensitive to the specific values; we are primarily interested in the continuum shape and the integrated line fluxes as temperature diagnostics. We defer detailed line profile analysis to a follow-up paper (in preparation).

Abundances of Mg, Si, S, and Ar were left free; for 2T fits, these parameters are degenerate with T. The result is a provisional model that gives a reasonable match to the observed spectrum. It required temperatures of about 6.0 and 14 MK, with emission measures of $2.6\times {10}^{55}$ and $2.3\,\times {10}^{54}\,{\mathrm{cm}}^{-3}$.

The absorption term required a column of ${N}_{{\rm{H}}}\sim 4\,\times {10}^{21}\,{\mathrm{cm}}^{-2}$. For comparison, the interstellar column is only $2\times {10}^{20}\,{\mathrm{cm}}^{-2}$, using the color excess in the blue-to-visual color magnitudes, $E(B-V)=0.04$ (Howarth & van Leeuwen 2019) and the standard dust-to-gas ratio scaling relation, ${N}_{{\rm{H}}}\,\approx 5\times {10}^{21}\times E(B-V)\,{\mathrm{cm}}^{-2}$. While this spectral region is less sensitive to wind absorption than the overall X-ray spectrum, it is not negligible. Our model is not sufficient to apply to the entire HETG region, but since we are mainly interested in the highest temperature, this is well constrained (to about 0.1 dex). Improvements to a global model will come when line profiles are modeled in detail.

Using this provisional model, we evaluated the presence of "line-free" regions by taking the ratio of the total to the continuum-only model count spectrum, and defined line-free as regions where the deviation from the model continuum is below 5%. Over the 2–9 Å range, about 45% of the resolution elements are line-free (at HEG or MEG resolution; the lines are resolved in each channel). Above 9 Å, there are no such regions. Over 2–9 Å, 20% of the photon flux is from the line-free regions, whereas the complement 80% of the region's flux is about equally divided between line and continuum photons.

We attempted fits to just the line-free region, but the lower counts made for much more uncertain parameters. The hotter plasma continuum only becomes appreciable below about 3–4 $\,\mathring{\rm A}$, where the signal drops. The lines of Mg, Si, S, Ar, Ca, and Fe are needed to constrain plasmas having temperatures from 6 to 60 $\,\mathrm{MK}$. The lines contribute significantly to the counts improving the overall statistics, and they provide substantial leverage on temperature determination.

For a more physically motivated model than the two-temperature fit, such as might be encountered with an ensemble of shocks in a clumpy wind (see, for example Cassinelli et al. 2008; Ignace et al. 2012), we have adopted a power-law differential emission measure (DEM) model of the form,

$$\begin{eqnarray}&&d\mathrm{EM}/{dT}=[{n}_{e}{n}_{{\rm{H}}}\,{dV}/{dT}]={D}_{0}(r){T}^{-\beta },\end{eqnarray} \tag{ 1 }$$

where dEM/dT is the differential emission measure, V is the volume of X-ray-emitting plasma, D₀ is a normalization factor, n_e and n_H are the electron and hydrogen number densities, respectively, T is the plasma temperature which we limit to a range T_min, T_max, and the exponent β. In the Appendix we derive an analytic solution for the continuum spectral energy distribution in the limit of ${T}_{\min }=0$ (also see Brown 1974).

We have implemented a power-law DEM in ISIS using a weighted sum of AtomDB model spectra on an evenly spaced logarithmic temperature grid with parameters T_min, T_max, the DEM exponent β, an overall normalization term (proportional to the volume emission measure, but not identical to D₀ in Equation (1) because the stellar distance is also included), the plasma model parameters shared by all components (Gaussian broadening, Doppler shift, and abundances), with multiplicative absorption.

Our free parameters are the normalization, β, ${T}_{\min }$, ${T}_{\max }$, ${N}_{{\rm{H}}}$, and abundances of the most significant species in the region, Mg, Si, S, and Ar. We have obtained fits with these free parameters, but the statistic's minimum is somewhat complicated with degenerate or poorly constrained parameters—better minima are often found when additional searches are performed, such as for confidence limits on parameters. Hence, we use a Markov Chain Monte Carlo method (an ISIS implementation of the Foreman-Mackey et al. 2013, Python code) to search parameter space and form confidence contours.¹⁶ From these runs, we could easily see a strong degeneracy between the normalization and the minimum temperature. Since the minimum temperature is not an interesting parameter (and is set to zero in the theoretical derivation in the Appendix), we froze ${T}_{\min }$ to the best-fit value of 6.0 dex K; it was always nearly a factor of 10 or more below ${T}_{\max }$. Furthermore, the shape of the spectrum, which is the primary diagnostic of the highest temperature plasma, is not very sensitive to ${T}_{\min }$. The fits then produced closed confidence contours for the normalization, and did not change other confidence limits. This does introduce an overall systematic uncertainty in the model, but not one that affects our primary results.

Figure 1 shows the data, model, and residuals. We only fit the spectrum below $9\,\mathring{\rm A}$, and in a few near-continuum bins near $13\,\mathring{\rm A}$ which are important for constraining N_H. We also ignored the Si xiii intercombination and forbidden lines near $6.7\,\mathring{\rm A}$, since these are strongly affected by the photospheric UV radiation field, as can be seen in the residuals when we evaluate the model over this region. (In our model, we also set the Fe relative abundance to 0.5, since this gave a better extrapolation of our fit to the 10–12 Å region, but it is otherwise not important to the spectral region fit.)

Zoom In Zoom Out Reset image size

Figure 2 shows parameter-pair confidence contours for the normalization, maximum temperature, absorbing column, and differential emission measure exponent. The relative abundance parameter pair confidence contours are shown in Figure 3. The reference abundances relevant to the region fit are listed in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Reference Abundances

Z	Elem.	$\mathrm{log}N$	Mass Frac.
		(dex)
1	H	12.00	5.98 × 10−1
2	He	11.21	3.85 × 10−1
6	C	7.60	2.85 × 10−4
7	N	9.10	1.05 × 10−2
8	O	8.14	1.30 × 10−3
10	Ne	8.02	1.25 × 10−3
12	Mg	7.68	7.06 × 10−4
13	Al	6.54	5.55 × 10−5
14	Si	7.60	6.63 × 10−4
16	S	7.21	3.08 × 10−4
18	Ar	6.49	7.32 × 10−5
20	Ca	6.43	6.40 × 10−5
26	Fe	7.59	1.29 × 10−3

Note. Abundances for H, He, C, N, and O are from Bouret et al. (2012), and the remainder are from the solar mass-fraction values of Asplund et al. (2009).

Download table as:  ASCIITypeset image

In fitting the cumulative spectrum, we have assumed that the mean flux is representative of the typical flux and that we do not have large temporal variability in the spectrum. To assess global variability, we take our power-law emission measure model, which gives a reasonable empirical match to the mean spectrum (as shown in Figure 1) and we fit only the normalization for each individual observation; we use the model spectrum and the calibration to remove any instrumental effects that might be present in a simple count-rate analysis. Tracking changes in the normalization over time will give us a first-order assessment of variability. It is within the realm of possibility that the spectral shape changes without changing the normalization, but that would require enough special circumstances to render it unlikely. Figure 4 shows the model normalizations as a function of the start times of the individual observations. While there are deviations from the mean of 0.226, the maximum fluctuations are about 12%, while most points are within one standard deviation (0.011) of the mean as expected for a constant source. While the very first HETGS spectrum of ζ Pup from 2000 has the lowest flux, its spectral distribution is essentially identical to the recent epoch, so we chose to include it in the analysis; it contributes about 9% of the counts below $10\,\mathring{\rm A}$. Hence, we believe fitting the cumulative spectrum is meaningful.

Zoom In Zoom Out Reset image size

The highest temperature is a diagnostic of the maximum relative velocity of the plasma undergoing strong shocks (see, for example, Stevens et al. 1992; Feldmeier et al. 1997; Cassinelli et al. 2008; Ignace et al. 2012). The maximum plasma temperature best fit for ζ Pup is $7.08\pm 0.03\,$ dex MK (or $1.04\,\mathrm{keV}$) (90% confidence). Using the relation for a strong shock, ${kT}=\tfrac{3}{16}\mu {m}_{{\rm{H}}}{v}^{2}$, where T is the post-shock temperature and v is the preshock velocity, and with the mean molecular weight μ = 0.67 for the adopted abundances of highly ionized elements (mainly H- and He-like ions), and using the approximation for the number of stripped electrons from Rutherford et al. (2013, Appendix C.7), we find a maximum relative shock velocity v ≈ 900 km s⁻¹ (with a formal statistical error of ∼50 km s⁻¹ based on the temperature uncertainty). As this velocity well below the wind terminal speed, it is within the capabilities of radiative acceleration to achieve, and would be slow for a wind–wind collision in a binary. For wind expansion with a "beta-law" exponent 0.9 (see, for example, Eversberg et al. 1998), a wind speed of 900 km s⁻¹ occurs at about 1.5 stellar radii above the photosphere and accelerates to about 1200 km s⁻¹ by 2 stellar radii. The line-based analysis of Cohen et al. (2010) found that X-rays shortward of $10\,\mathring{\rm A}$ first appear at a turn-on radius consistent with 1.5 stellar radii. Their results are consistent with our finding that the absorption is relatively low shortward of $10\,\mathring{\rm A}$, so any absence of emission below 1.5 stellar radii could not be due to absorption effects. However, it has not been shown that an absence of emission at low radii is necessary to fit the X-ray emission line profiles, and additional evidence that X-ray creation may extend much deeper will be addressed in a separate publication.

The presence of a stand-off distance for the initiation of X-ray creation could be taken as evidence of line-deshadowing-instability (LDI) activity, as it is a self-excited instability that could take some distance to appear. The absence of a stand-off distance would suggest instead that variations in the acceleration right from the surface are the crucial cause. Both of these effects are considered in Sundqvist & Owocki (2013) where they find that surface perturbations are essential for generating strong shocks at low radii, and our analysis is consistent with early appearance of strong shocks, although a more in-depth analysis is required. If surface-related effects are the crucial element, then the main explanatory factor for the hard X-ray heating must trace back to the cause of variations at the boundary. Such variations could be inherent in the line-driving critical point conditions (Sundqvist & Owocki 2015) or might require a specific physical mechanism such as subsurface convection (Cantiello et al. 2009; Cantiello & Braithwaite 2011).

Our spectral fits led to a differential emission measure exponent of 2.6 (±0.2 90% confidence). The AtomDB continuum is not purely Bremsstrahlung, but also contains radiative recombination edges and other contributions, so has additional features not included in the analytic solution. However, if we evaluate the AtomDB's true continuum for our differential emission measure plasma model and compare to Equation (A6) using Table 3, we find an exponent near 5/2. Comparison of the spectrum to a simple power-law energy distribution also shows us that the curvature is significant, appearing below $5\,\mathring{\rm A}$, and that we have indeed detected the high-temperature cutoff. Our best fit is in very good agreement with the value of 7/3 which the hydrodyamical calculations of Cassinelli et al. (2008) predicted for a bow-shock around a single clump. The models shown by Krtička et al. (2009) have a power-law index of about 2.4 (see their Equation (10)), which is also consistent with our fitted slope. We note that their results also required enforcing some type of surface variations.

Table 3.  Solutions of ${\gamma }_{m}$ for Select Values of m

m	β	${\gamma }_{m}(x)$
0	3/2	${e}^{-x}$
1	5/2	$(x+1){e}^{-x}$
2	7/2	$({x}^{2}+2x+2){e}^{-x}$
3	9/2	$({x}^{3}+3{x}^{2}+6x+6){e}^{-x}$

Download table as:  ASCIITypeset image

The best-fit value of ${N}_{{\rm{H}}}$ of about $0.7\times {10}^{22}\,{\mathrm{cm}}^{-2}$ seems fairly large, and also causes a spectral model extrapolation to longer wavelengths to lie well below the data beyond about $17\,\mathring{\rm A}$. However, we do not expect a constant ${N}_{{\rm{H}}}$ to apply to the entire range in our approximation, since at longer wavelengths, the optical depth is larger, and one cannot see as geometrically deeply into the wind. Furthermore, the phabs model is not valid (for stellar wind absorption modeling) at longer wavelengths because it does not incorporate warm absorbers, which have edges at different wavelengths than neutral species. Below $9\,\mathring{\rm A}$, in our range of interest, the equivalent optical depths rapidly drop below unity (Oskinova et al. 2006), and we maintain our assumption that the intrinsic continuum properties dominate in the model; the standard phabs model is accurate here in determining the spectral curvature. The value of ${N}_{{\rm{H}}}$ we derived is consistent with the Fe xvii, $15\,\mathring{\rm A}$ line profile fitting of Cohen et al. (2010), who found that the optical depth on the central ray (${\tau }_{* }$) was about 2.0, and this is equivalent to ${N}_{{\rm{H}}}\sim 0.5\times {10}^{22}\,\,{\mathrm{cm}}^{-2}$.

The shock heating rate was studied by Zhekov & Palla (2007), and by Cohen et al. (2014a) who placed upper limits on the heating rate above $T\sim 10\,\mathrm{MK}$, based on emission lines alone. We do firmly detect the high-temperature lines of Ar xvii and Ca xix, whose maximum emissivities occur near $20\,\mathrm{MK}$ and $30\,\mathrm{MK}$, respectively (7.3 and 7.5 dex). However, our spectral fits do not require plasmas this hot, so their emission comes from the low-temperature tail of their emissivity curves; at $10\,\mathrm{MK}$, Ar xvii is at about 25% of its maximum emissivity, while Ca xix is at about 10%. Nondetection of both Ar xviii and Ca xx is consistent with this upper temperature limit.

There is a suggestion in the data of excess counts from Fe xxv ($1.85\,\mathring{\rm A}$), whereas the model's line emission here is quite weak. According to the model and the atomic database, the features in the model are due to a blend of many weak Fe xxiii and Fe xxiv lines, whose emissivities peak near $32\,\mathrm{MK}$. The He-like Fe xxv emission lines peak near $63\,\mathrm{MK}$ and at $32\,\mathrm{MK}$ are still about 8 times stronger than the Fe xxiii–Fe xxiv blends. However, at the $\sim 10\,\mathrm{MK}$ maximum model temperature, Fe xxiii and Fe xxiv are only at 1% of their peak emissivity, while Fe xxv is nearly an order of magnitude fainter than that. We do not think there is any significantly detected signal in Fe xxv.

Instrumental background also begins to dominate below $\sim 2\,\mathring{\rm A}$, as can be seen in the large errorbars in Figure 1; the upturn is due to amplification of the background noise. The lack of significant Fe emission and the continuum below $\sim 3\,\mathring{\rm A}$ provide good constraints on the shape of the spectrum and support the power-law emission measure fit here.

Abundance determinations rely on the line-to-continuum flux ratio, the theoretical emissivities, and the emission measure distribution, assuming collisional ionization equilibrium in an optically thin plasma. Mass-loss rates derived from emission lines are proportional to the adopted abundances, since they scale with the opacity (see, for example, Oskinova et al. 2006; Cohen et al. 2014b). We have assumed the emission line fluxes are only affected by foreground absorption, which is not rigorously true, since both the absorption and emission occur in the wind. We know that lines are subject to continuum absorption with Doppler effects from the wind, since line centroids are blueward of their rest wavelengths, and the red wings are affected by continuum absorption. These effects are mitigated, however, at the shorter wavelengths where the continuum opacity is smaller.

We adopted our reference abundance set for the nuclear-processed elements from Bouret et al. (2012), and the remainder as solar from Asplund et al. (2009). We give the mass fraction and logarithmic number fractions in Table 2. We allowed the abundance parameters to range from half the reference to 1.5 times the reference values. Our results are shown in Figure 3 and show near-reference values for S and Ar, and slightly lower values for Mg and Si. They generally agree with the values found by Zhekov & Palla (2007; but with significantly smaller uncertainties). Refinement of abundances will require more detailed line-based analysis and consideration of systematic effects due to assumed temperature distributions.

This deep Chandra/HETG spectrum of ζ Pup has allowed us to rigorously determine that there is little to no plasma above $\sim 12\,\mathrm{MK}$ for a power-law differential emission measure model, and that the emission measure distribution is consistent with that expected from shocked clumps deep in the wind. Assuming efficient X-ray creation (though see Steinberg & Metzger 2018), this limits its relative preshock velocities to about 900 km s⁻¹. This is below the wind velocity in the acceleration zone, but is still difficult to reconcile with traditional LDI without some perturbations, which can then be amplified by LDI. The hydrodynamical calculations of Sander et al. (2017) further exacerbate this: they found that the ζ Pup wind acceleration cannot be described by a simple beta-law, but accelerates much more slowly, reaching only $0.1\times {v}_{\infty }$ (≈200 km s⁻¹) by 1.5 stellar radii above the photosphere. The inferred inefficient radiative acceleration deep in the wind presents additional challenges to understanding strong shocks within the purview of pure radiative acceleration, so additional physical mechanisms may be required. In addition, the 2D and 3D simulations of Steinberg & Metzger (2018) also present tension with our high temperature at small radii conclusion. Their simulations find that post-shock temperatures can be one or two orders of magnitude below the usual strong-shock prediction in the case of highly radiative shocks, though it is unclear how severe this restriction would be for much more weakly radiative shocks. Hence, there is clear motivation for additional theoretical and observational work to understand and reconcile these apparent differences.

Additional constraints will come from detailed line-profile modeling using radial wind structural and dynamic models, and will provide a better idea of the effective absorbing column and region of X-ray formation. The assumed foreground absorbing screen with a single value for ${N}_{{\rm{H}}}$ is an approximation; line profile fitting using detailed wind structure models will help with the interpretation of the spectral shape, especially at longer wavelengths, and also lead to important quantities such as the mass loss rate. At the short wavelengths of primary interest here, the absorption is mainly affecting the normalization and spectral curvature, but not the high energy cutoff—we are largely decoupled from the details of wind structure. Hence, we are confident that the maximum temperature is of the order we have found, and this will constrain plasma properties in future models.

Support for D.P.H. for this work was provided by NASA through the Smithsonian Astrophysical Observatory (SAO) Grant GO8-19011C to MIT, and by contract SV3-73016 to MIT for Support of the Chandra X-Ray Center (CXC) and Science Instruments.

Chandra General Observer Program, Cycle 19 supported W.W. by GO8-19011A, J.S.N. by GO8-19011B, N.A.M. by GO8-19011D, N.D.R. by GO8-19011E, and R.I. by GO8-19011F.

A.F.J.M. acknowledges financial support from NSERC (Canada) and FQRNT (Quebec).

N.A.M. also acknowledges support from the UWEC Office of Research and Sponsored Programs through the sabbatical and URCA programs, and from a Chandra Research Visitor award.

J.S.N. also acknowledges support by CXC.

Y.N. acknowledges support from the Fonds National de la Recherche Scientifique (Belgium), the Communauté Française de Belgique (for financial support of the OHP campaign), the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme (contract XMaS).

CXC is operated by SAO for and on behalf of NASA under contract NAS8-03060.

This research has made use of ISIS functions (ISISscripts) provided by ECAP/Remeis observatory and MIT (http://www.sternwarte.uni-erlangen.de/isis/).

Facility: CXO (HETG/ACIS). -

Software: CIAO (Fruscione et al. 2006), ISIS (Houck & Denicola 2000).

Consider an emission measure distribution that is a power law in temperature as given by

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{EM}}{{dT}}=\displaystyle \frac{{n}^{2}(r)}{{T}_{0}}{\left(\displaystyle \frac{T}{{T}_{0}}\right)}^{-\beta }\,{dV},\end{eqnarray} \tag{ A1 }$$

where EM is the emission measure, T is temperature, n(r) is the number density of the X-ray emitting plasma at radius r, T₀ is a temperature constant at that radius, β is the power-law exponent at that radius, and dV is the volume element for a shell. In what follows, both T₀ and β are taken as constants for all radii.

The cooling function from bremsstrahlung is given by

$$\begin{eqnarray}&&{\rm{\Lambda }}(E,T)={{\rm{\Lambda }}}_{0}{\left(\displaystyle \frac{{T}_{0}}{T}\right)}^{1/2}\,{e}^{-E/{kT}},\end{eqnarray} \tag{ A2 }$$

where ${{\rm{\Lambda }}}_{0}$ is a constant. Note that there are several processes that can contribute to the emissivity in the X-ray band. At the short wavelengths of interest for this application, two-photon emission is negligible. Bound-free emission can be important across the X-ray band. However, for a wind that is predominantly H and He, such as ζ Pup, and given relatively high-temperature plasmas that contribute to the short wavelengths, Landi (2007) showed that Bremsstrahlung is most important below around $10\,\mathring{\rm A}$. Even when bound-free has a nonnegligible contribution, as long as edges are not important in the wave band of interest, the cooling has the same scaling as Bremsstrahlung. In this scenario, bound-free would only impact the proportionality constant, whereas our focus is on slope information from the observed continuum. With these caveats, the following shows how the observed continuum relates to a high-temperature cutoff combined with a power-law differential emission measure.

The spectral energy distribution (SED) for the continuum is dL/dE:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dL}}{{dE}} & = & \displaystyle \int \,{\rm{\Lambda }}(E,T)\displaystyle \frac{d\mathrm{EM}}{{dT}}\,{dT},\\ & = & {{\rm{\Lambda }}}_{0}\,\displaystyle \int {\left(\displaystyle \frac{{T}_{0}}{T}\right)}^{1/2}\,{e}^{-E/{kT}}\,\displaystyle \frac{{n}^{2}(r)}{{T}_{0}}{\left(\displaystyle \frac{T}{{T}_{0}}\right)}^{-\beta }\,{dT}\,{dV}\\ & = & \displaystyle \frac{{L}_{0}}{{{kT}}_{0}}\,\displaystyle \int {\left(\displaystyle \frac{{T}_{0}}{T}\right)}^{1/2}\,{e}^{-E/{kT}}{\left(\displaystyle \frac{T}{{T}_{0}}\right)}^{-\beta }\,\displaystyle \frac{{dT}}{{T}_{0}},\end{array}\end{eqnarray} \tag{ A3 }$$

where in the final line, the integration over volume has been carried out and a constant with units of luminosity, L_0, was introduced. For the present study we assume that wind absorption can be ignored.

If the temperature ranges from zero to infinity, then formally, the integration above will produce an SED with ${dL}/{dE}\propto {E}^{-\beta +1/2}$. A high-temperature cutoff leads to a downturn of the continuum toward higher energies. We now assume that the constant T₀ represents the maximum temperature achieved in every shell of the wind. The integration limits are zero to T₀. It is convenient to introduce a change of variable for normalized energy, with dummy variable t = E/kT and normalized energy $x=E/{{kT}}_{0}$. The integral now becomes

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{{dE}}=\displaystyle \frac{{L}_{0}}{{{kT}}_{0}}\,{x}^{-\beta +1/2}\,{\int }_{x}^{\infty }\,{t}^{\beta -3/2}\,{e}^{-t}\,{dt}.\end{eqnarray} \tag{ A4 }$$

It is convenient to introduce

$$\begin{eqnarray}&&m=\beta -3/2,\end{eqnarray} \tag{ A5 }$$

in which case the specific luminosity becomes

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{{dE}}=\displaystyle \frac{{L}_{0}}{{{kT}}_{0}}\,{x}^{-m-1}\,{\gamma }_{m}(x),\end{eqnarray} \tag{ A6 }$$

where

$$\begin{eqnarray}&&{\gamma }_{m}(x)={\int }_{x}^{\infty }\,{t}^{\beta -3/2}\,{e}^{-t}\,{dt}.\end{eqnarray} \tag{ A7 }$$

Physically, ${x}^{-m-1}$ is the power-law SED when there is no maximum temperature cutoff; it is ${\gamma }_{m}(x)$ that eventually leads to an exponential decline in the SED at high energies. This integral has an analytic solution when m is a nonnegative integer, given by

$$\begin{eqnarray}&&{\gamma }_{m}(x)=\left[{x}^{m}+{{mx}}^{m-1}+m(m-1){x}^{m-2}+...+m!\right]\,{e}^{-x}.\end{eqnarray} \tag{ A8 }$$

Here the exponential decline in normalized energy is made explicit. Table 3 provides values of β and solutions of ${\gamma }_{m}$ for $m=0,1,2,$ and 3.

Note that the energy flux spectrum is proportional to ${dL}/{dE};$ the photon flux spectrum is proportional to ${E}^{-1}\,{dL}/{dE}\,\propto {x}^{-m-2}{\gamma }_{m}(x)$.