An ‘Analytic Dynamical Magnetosphere’ formalism for X-ray and optical emission from slowly rotating magnetic massive stars
Slowly rotating magnetic massive stars develop “dynamical magnetospheres” (DM’s), characterized by trapping of stellar wind outflow in closed magnetic loops, shock heating from collision of the upflow from opposite loop footpoints, and subsequent gravitational infall of radiatively cooled material. In 2D and 3D magnetohydrodynamic (MHD) simulations the interplay among these three components is spatially complex and temporally variable, making it difficult to derive observational signatures and discern their overall scaling trends. Within a simplified, steady-state analysis based on overall conservation principles, we present here an “analytic dynamical magnetosphere” (ADM) model that provides explicit formulae for density, temperature and flow speed in each of these three components – wind outflow, hot post-shock gas, and cooled inflow – as a function of colatitude and radius within the closed (presumed dipole) field lines of the magnetosphere. We compare these scalings with time-averaged results from MHD simulations, and provide initial examples of application of this ADM model for deriving two key observational diagnostics, namely hydrogen H- emission line profiles from the cooled infall, and X-ray emission from the hot post-shock gas. We conclude with a discussion of key issues and advantages in applying this ADM formalism toward derivation of a broader set of observational diagnostics and scaling trends for massive stars with such dynamical magnetospheres.
keywords:MHD — Stars: winds — Stars: magnetic fields — Stars: early-type — Stars: mass loss — Stars: X-rays
Hot luminous, massive stars of spectral type O and B have dense, high-speed, radiatively driven stellar winds (Castor et al., 1975). In the subset (10%) of massive stars with strong (kG), globally ordered (often significantly dipolar) magnetic fields (Petit et al., 2013), the trapping of the wind outflow by closed magnetic loops leads to the formation of a circumstellar magnetosphere. For stars with moderate to rapid rotation – such that the Keplerian corotation radius lies within the Alfvén radius that characterizes the maximum height of closed loops –, the rotational support leads to formation of a “centrifugal magnetosphere” (CM), wherein the trapped wind material accumulates into a relatively dense, stable and long-lived rigidly rotating disk (Townsend & Owocki, 2005).
In contrast, for magnetic massive stars with slow rotation, and thus , this trapped material falls back to the star on a dynamical (freefall) timescale, representing then a “dynamical magnetosphere” (DM) (Sundqvist et al., 2012; Petit et al., 2013).
Because of the rotational spindown associated with angular momentum loss through the relatively strong, magnetized wind upflow from open field regions (ud-Doula
et al., 2009), magnetic O-type stars are typically
In such DM’s the trapped wind upflow from opposite footpoints of closed loops collides near the loop apex, forming shocks that heat the gas to X-ray emitting temperatures; as this gas radiatively cools, it falls back to the star as a gravitational downflow. 2D and 3D magnetohydrodynamic (MHD) simulations of such DM’s (ud-Doula & Owocki, 2002; ud-Doula et al., 2013) show a complex and variable interplay among all three components, and this makes it difficult to derive observational signatures and discern their overall scaling trends.
Applications of these MHD models have thus been limited to a few selected O-stars, using simplified radiative transfer methods to derive synthetic spectra for H emission lines (Sundqvist et al., 2012; Grunhut et al., 2012; ud-Doula et al., 2013; Wade et al., 2015) and ultra-violet wind resonance lines (Marcolino et al., 2013; Nazé et al., 2015). These initial studies have provided strong support of the general DM concept; however, the complexity of the time-dependent 3D structure, together with computational expense of the simulations, prohibits more systematic and detailed computations of synthetic observables across the electromagnetic spectrum, as well as application to larger samples of magnetic massive stars.
Similar arguments apply to X-ray spectral synthesis. Detailed MHD simulations have been used to analyze the high spectral resolution X-ray observations available for a few selected OB stars (Petit et al., 2015; Nazé et al., 2015). But for the analysis of the much larger number of stars with low-resolution X-ray data, an analytic model can capture key observable properties and trends. Recently, ud-Doula et al. (2014, hereafter paper I) carried out an extensive MHD simulation study of radiative cooling of the hot post-shock gas in DM’s, with a focus on deriving the resulting X-ray emission as a function of the density-dependent cooling efficiency. When interpreted in terms of a simplified “X-ray analytic dynamical magnetosphere” (XADM) analysis, this led to predicted scaling laws for variation of X-ray luminosity with the wind mass loss rate. Subsequent application by Nazé et al. (2014) toward interpreting observationally inferred X-ray luminosities in a large sample of magnetic massive stars showed that, with some fixed factor adjustment to account for the limited “duty cycle” for X-ray production during the complex variations seen in MHD simulations, this XADM scaling matched quite well the observed trends for the subset of magnetic massive stars with slow enough rotation to have DM’s.
The analysis here builds on this success to develop a more general “analytic dynamical magnetosphere” (ADM) formalism that now provides explicit formulae for the spatial variation of density and flow speed (as well as temperature for the hot gas) in all three components of the closed loop magnetosphere: wind upflow, hot post-shock gas, and cooled downflow (section 2). The overall ADM analysis here is guided and tested by comparison with time-averaged results from full MHD simulations of DM’s (section 3), using the Alfvén radius in the MHD simulation to define a maximum closure radius of a dipole loop in the ADM model.
The inclusion of the cooled downflow now allows for modeling of optical emission lines like H- (section 4). Moreover, the description of the spatial distribution of both the hot and cool components allows an extension of the XADM analysis of paper I (section 5), for example by accounting for possible bound-free absorption of emitted X-rays by the cool wind and downflow (see section 5.3 and Petit et al., 2015). We conclude (section 6) with a brief summary of key issues and advantages in applying this ADM formalism toward deriving a broader range of spectral diagnostics. Appendices A and B show how the ADM model can be used to derive both the X-ray differential emission measure (DEM), as well as a shock-temperature distribution . Appendix C gives background on the ADM scalings for H emission.
2 Analytic Dynamical Magnetosphere
2.1 Basic magnetic scalings for a star-centered dipole
The ADM formalism is based on an idealized, steady-state analysis of the mass flow along closed magnetic loops that are assumed to follow the individual lines of a simple
Following figure 1, if we specify the colatitude from the north dipole axis by , then a given dipole loop line that intersects the stellar radius at extends over a band about the equator, with radius variation given by
and with the maximum radius at the loop apex on the magnetic equator,
For any position within the magnetosphere, the magnitude of the field follows the dipole scaling,
For such a dipole, the direction of the field only depends on , given by the unit vector
where and are unit vectors in the radial and latitudinal directions.
The nearly full ionization of circumstellar material around hot, OB stars implies a high conductivity and thus broad applicability of the frozen flux theorem of ideal MHD. In the context of the present model of a rigid, closed dipole field line, this means any material remains always locked onto its given field line, with flow velocity parallel to the local field direction . In terms of the local density , the mass flux thus has a divergence given by,
where the third equality uses , and the last equality defines a coordinate distance along the field line, with . This means that in any steady-state flow, wherein mass conservation requires , the local mass flux density scales with the field strength, so that constant along the flow. For a given steady input mass flux from the surface, the spatial variation of density can thus be derived from knowledge of the flow speed (or vice versa), in terms of the known spatial variation of field strength from (3).
2.2 3-Component model for mass trapped in a closed dipole loop
In full MHD simulations (see paper I), the actual flow along such closed loops is spatially structured and time-dependent, with any given loop alternating between variable intervals and regions of upflow and downflow. As demonstrated below, however, the overall conservation principles mean that, when averaged over time, and/or over many separate loops with complex structure, key characteristics from 2D and 3D MHD simulations can be relatively well characterized by an idealized ADM model that assumes each loop line can simultaneously support two oppositely directed, steady-state flows.
Specifically, as illustrated in figure 1, this ADM analysis distinguishes three distinct components of material flow within the loop:
Hot post-shock gas
Mass is fed into the loop by the wind upflow, driven by the radiative flux from the underlying star. Relative to the surface mass flux for the non-magnetic case, the mass flux fed into a loop with footpoint scales as (Owocki & ud-Doula, 2004)
where is the radial projection cosine of the local surface field
The flow speed along the loop should in principle be computed from solution of the acceleration from radiative driving, accounting for the curvature, tilt and areal divergence along the loop; but to maintain analytic tractability, our ADM analysis simply assumes this speed follows a canonical “beta” velocity law,
where the maximum speed is given here by the expected terminal speed of a corresponding unmagnetized wind, and we assume the simple case with .
As the upflow from a loop footpoint approaches the loop apex at radius , where the scaled speed reaches its maximum loop-specific value , the collision with flow from the opposite footpoint induces a pair of reverse shocks, one on each side of the loop apex. This converts the wind kinetic energy into heat, resulting then in the hot post-shock gas that extends away from the loop apex by a distance set by the post-shock cooling length . For lower luminosity stars with smaller mass loss rate and thus a lower-density wind upflow, the cooling length can become comparable to the loop apex radius, . As discussed in paper I, this forces a “shock retreat” to a lower shock radius , with thus a slower scaled shock speed and so a cooler post-shock temperature . Using the dipole shock-retreat analysis in Appendix B of paper I (recapitulated in section 2.2.2 below), figure 2 illustrates the progressive retreat of the shock with increased cooling length, as characterized by a cooling parameter , given below in equation (13).
Starting from this immediate post-shock temperature of order many MK, radiative cooling in the post-shock region causes a decrease in temperature along the loop, eventually reaching near the loop apex the same stellar photoionization equilibrium temperature of the wind upflow, typically on the order of the stellar effective temperature of a few kK. The subsonic nature of the post-shock flow means this cooling layer is almost isobaric, with the nearly constant pressure implying a strong increase in density as the temperature declines.
Much like a ball balanced at the top of a hill, the convex upward nature of the magnetic loop near its apex means such dense, cooled material is gravitationally unstable, and so begins a gravitational free-fall along the loop to form the cooled downflow component of the DM. The strong compression, and relatively slow infall speed, means the density is much higher than in the wind upflow, and the several kK equilibrium temperature means it can efficiently radiate in hydrogen emission line series, leading then to prominent optical emission in lines like H (Howarth et al., 2007; Sundqvist et al., 2012) that are much stronger than from the wind. Moreover it is cool enough that there is a significant X-ray bound-free opacity from abundant, partially ionized heavy elements like CNO and Fe, which can thus play a role in attenuating the X-rays emitted from the hot post-shock gas (Petit et al., 2015).
As noted, in 2D and 3D MHD simulation of DM’s (ud-Doula et al., 2013, paper I), these three components of wind upflow, hot post-shock gas, and cooled downflow become mixed together in complex, highly variable combinations, with 3D models showing structure down to scales of or less, such that even loop lines separated by such a small scale can be in markedly different phases for the cycle of shock-build-up, cooling, and then infall (ud-Doula et al., 2013).
But as we demonstrate in section 3 below, if one takes a suitably long time average, covering many such infall cycles, then these stochastic variations become smoothed out, leaving a distinctive and organized large-scale spatial distribution. This can be well characterized by a superposition of the density, velocity and temperature of the wind upflow, hot post-shock gas, and cooled downflow. The next subsections quantify this in terms of relatively simple ADM scalings based on quasi-steady-state conservation applied to the material from each component.
The conditions in the wind upflow are quite straightforward to specify. As noted, the speed is assumed to follow the law given in equation (7),
The density then follows from steady-state mass continuity, with base mass flux (6), and accounting for the dipole area divergence,
where is a characteristic density for an unmagnetized wind with mass loss rate and terminal speed .
Finally, the temperature of the wind upflow is expected to be of order the stellar effective temperature, but its actual value is not relevant to the ADM model, so long as the upflow is sufficiently supersonic to justify use of the strong-shock scaling in modeling the post-shock flow, as discussed next.
Hot post-shock gas
Let us next derive conditions in the hot post-shock gas in the region between the shock radius and loop apex, . Our analysis builds on and extends the X-ray analytic dynamical magnetosphere (XADM) formalism developed in paper I, deriving now explicit expressions for the spatial variation of the temperature in this post-shock cooling layer.
We begin with the equation for advective change of the specific enthalpy due to radiative cooling, characterized by a radiative cooling function , which we write in a mass-weighted form , with and respectively the mean mass per proton and per electron (see section 2 of Kee et al. (2014), and section 2.5 and Appendix B of paper I),
where is the mean molecular weight, and the temperature derivative is with respect to the field line coordinate . The nearly constant pressure of this post-shock layer implies a strong increase in density as temperaure declines.
Using steady mass flux conservation and this near constancy of the pressure , we can eliminate both the density and the flow speed in favor of the temperature . Following the analysis in Appendix B of paper I, we can then integrate (11) to obtain an implicit solution for the temperature decline from the post-shock value to the much lower (near zero) temperature at the loop apex,
here , and accounts for the mass loss weighting for a given field-line flow tube, defined as a fraction of the spherical mass loss used in the definition of the cooling parameter (taken from equation (25) of paper I),
with scaled values cm, cm/s, and /yr.
Following eqns. (B11)-(B14) of paper I, the path integral of the field strength along the loop can be evaluated as
with the absolute value operations ensuring that remains positive even in the lower hemisphere, where . Applying the boundary condition that the temperature nearly vanishes at the loop apex, and so formally taking , we find, since , that the spatial variation of temperature in this hot component of the ADM can be written as
where, for the assumed case of a highly supersonic outflow yielding a strong shock, the immediate post-shock temperature , with
In practice, to account for the effects of photoionization heating by the underlying star, we do not allow the temperature to fall below the stellar effective temperature,
where here we assume a typical hot-star value 30,000 K.
The shock radius is obtained by solving a transcendental equation for the dipole shock retreat, as derived in equation (B16) of paper I,
Figure 3 of paper I plots the variation of the scaled shock speed vs. cooling parameter , for various loop apex speeds . The red lines in figure 2 here illustrate the progressive spatial retreat of the shock away from the loop apex at the magnetic equator as the radiative cooling parameter is increased from a small () to large () values in steps of 1 dex.
Let us next obtain the spatial variation of the density for this hot, post-shock region. For a strong shock, the density of the immediate post-shock gas is simply a factor four times the incoming wind density at the shock, . Since the pressure is nearly constant in this subsonic post-shock cooling layer, we can write the spatial variation of density of the hot gas as
Using mass continuity, we can also readily derive the spatial variation of the post-shock flow speed. For a strong shock, the immediate post-shock speed is just a quarter of the incoming wind speed, , while the spatial variation is given by
where the field magnitude is given by (3).
The top two rows of figure 3 plot results for the spatial variation of (top) and (middle), for the 3 cooling parameter values 0.1, 1, and 10 (left, middle, and right columns).
Once this hot post-shock gas cools, the stellar gravity pulls the cooled material back to the star, accelerating from near zero speed at the loop apex at , into a cooled downflow along the loop.
Using mass conservation, the associated density can be computed in a completely analogous way as the wind upflow density in equation (9), but now replacing the wind upflow speed with the cooled downflow speed ,
where is a characteristic density for the downflow. To account for the fact that this infall is typically initiated from some finite length away from the exact loop apex, in the denominator here we have replaced the factor in equation (22) for the speed with
which has the effect of smoothing the density near the magnetic equator over the scale . Using this and equations (1), (3), and (6), we can rewrite (23) in a form that explicitly shows the dependence on radius and co-latitude,
The value of this apex smoothing length can be set based on results from MHD simulations, or derived from comparison to observations (see section 4).
The colorplot in figure 4 shows the spatial variation of this cooled downflow density, , for 3 selected values of the smoothing length, 0, 0.1, and 0.3 (left, center, right). The top row of figure 5 shows the density for both the wind and cooled downflow for the case (left panel), along with latitudinal and radial components of the downflow velocity (middle and right panels). The lower panel shows corresponding time-averages from full MHD simulations, as discussed further in section 3.1.
Finally, as with the wind upflow, the temperature of the cooled downflow is expected to be of order the stellar effective temperature, but its actual value and spatial variation is not set in this ADM formalism, and so can be modeled separately in the context of its relevance to each specific diagnostic.
2.3 Summary scalings for ADM components
To facilitate application of this ADM model in deriving observational diagnosts, let us collect here the scaling relations for temperature, speed, and density in each of the three model components.
For the wind upflow (section 2.2.1):
is given by equation (8);
is given by equation (10).
For the hot component (section 2.2.2):
is given by equation (21);
is given by equation (20),
For the cooled downflow (section 2.2.3):
is given by equation (22);
is given by equation (25).
In the plots here, lengths are in stellar radii , velocities are scaled by the stellar escape speed , densities are in , and temperatures are in .
The global parameters are: the maximum closed loop apex , fixed here to ; the ratio of wind terminal speed to surface escape speed, fixed here to ; the loop apex smoothing length , with standard value ; and the cooling parameter (defined in equation 13). All components also have floor temperature set to the stellar effective temperature, taken here to be K.
Finally, for all three components, the flow is along the dipole field line, with thus vectorial direction given by equation (4), with a radially positive sense for the upflow and post-shock components, and a negative sense for the cooled downflow.
3 Comparison to MHD simulations
To assess the potential for these simplified ADM scalings to provide a basis for deriving observational diagnostics, let us next compare their predictions to time-averaged results from full MHD simulations.
For consistency with the sample plots given above, we again choose as a standard the ADM model with a maximum loop closure radius . For corresponding MHD simulations, we use the standard stellar parameters of paper I, but now with a polar magnetic field G, which gives an Alfvén radius . The time averaging begins at ks, when the structure has fully relaxed from the initial condition, and extends to ks, representing about ten cycles of mass build-up and dynamical infall. To ensure the north-south symmetry of an arbitrarily long time average (or from azimuthal averaging of a full 3D model; see figure 3 of ud-Doula et al. (2013)), we also carry out a north-south averaging of all the time-averaged quantities from this 2D MHD simulation.
The remainder of this section focuses on the relatively cool ( K) material in the outflowing wind and the cooled downflow. This is then used in §4 to derive signatures in recombination-based emission lines like H. Section 5 next focuses on X-ray emission from the hot post-shock gas, and its absorption from the wind and cooled downflow.
3.1 Cooled downflow
First, for closed loop lines with apex radii , the top row of figure 5 plots ADM scalings for the density (with apex smoothing length ; left panel), and the latitudinal and radial components of velocity (middle and right panels) in the cooled downflow component. The lower panels compare corresponding time-averaged quantities from the MHD simulations. For simplicity, the lower left shows the time averaged density , but because emission profiles depend on the velocity of material weighted by its density-squared emission measure, the middle and right panels show the time averages of the speeds weighted by the density squared.
Note that in the MHD model the signatures of trapped material in closed magnetic loops is radially more extended than in the idealized, dipole form of closed loops in the ADM model, due to the dynamical interaction of the field and radiatively driven wind outflow. The dynamical variation of infall episodes means the density is less equatorially concentrated in the MHD vs. ADM models, but overall the MHD model density, as well as the radial and latitudinal velocities, show a clear boundary change linked to closed loop geometry, much as in the ADM case.
To make these comparisons semi-quantitive, the assumed stellar parameters of the MHD simulation (which are the same as the standard model of paper I) imply a surface escape speed km s, and a characteristic wind-fed loop density g cm. Using these to scale the quoted CGS values in the colorbars for the MHD simulation in the lower row of figure 5, we see that the color levels are in quite good agreement with the ADM case.
However, one key aspect of the cooled downflow not accounted for in the ADM steady-state picture is that in full MHD simulations, the material infall actually occurs in sporadic intervals of highly compressed, localized streams. This implies a significant enhancement in the mean of the density-squared, which for the two-body collision or recombination processes that underlie line emission is what sets the overall emission strength. For this MHD model, figure 6 plots now the time-averaged rms density . While the overall form is very similar to the density plot in figure 5d, note that, due to the clumped infall, the colorbar scale is now enhanced by about a factor 7. The associated “clumping factor” , which effectively sets the level of emission enhancement relative to a smooth model with the same average density, is thus typically increased by several factors of ten. In applying the ADM scalings to model such optical emission lines, one should account for this clumping by enhancing the emission by a clumping factor of this order. The next section provides a first example for the case of hydrogen H- line emission.
4 Application to Observational Diagnostics: Optical Line Emission
With this background, let us now derive scaling relations and perform a first diagnostic application of the ADM model for the hydrogen H- line. A key advantage with H is that under typical O-star wind conditions its atomic level populations remain quite close local thermodynamic equilibrium (LTE) with respect to the real population of ionized hydrogen (e.g., Puls et al., 1996; Sundqvist et al., 2011), allowing one to use relatively simple methods to analyze this line also in magnetic O-stars (Sundqvist et al., 2012).
4.1 ADM-modified scaling relation for H emission
The principal scaling of H emission in OB-star winds comes from considering the line optical depth in the Sobolev approximation. For Doppler width and directional Sobolev length , this is
where we have absorbed atomic constants of the H transition and dependencies on electron temperature and occupation number densities into the parameter , as given in Appendix C (see also Puls et al., 1996; Petrenz & Puls, 1996). For an observer viewing from above the magnetic pole, equation (26) can be used to derive an optically thin emission-measure scaling law for the ADM model (see Appendix C) for a polar-view observer:
where the function describes the dependence on the size of the closed-loop magnetosphere. Equation (27) illustrates explicitly how the standard scaling for non-magnetic wind emission is modified here by the two magnetic parameters and (essentially setting the ADM disc-thickness and size), and comparison to full radiative transfer calculations (Sect. 4.2) indicates that this simple scaling law captures quite well the principal scaling of ADM H emission under typical OB-star conditions.
4.2 A first diagnostic application
Building on the simple scaling analysis above, let us now make a first application toward using optical emission lines like H to diagnose the winds and magnetospheres of slowly rotating magnetic OB-stars. To compute synthetic spectra from the steady-state ADM density and velocity structure, we follow Sundqvist et al. (2012) and solve the formal solution of radiative transfer in a 3D cylindrical system for an observer viewing from angle with respect to the magnetic axis. Defining and as the angles that the magnetic axis and the observer line-of-sight make with the rotation axis, the variation of with rotational phase is given by
We solve the transfer equation only in the infall compoent, assuming H departure coefficients and an electron temperature structure calibrated by 1D NLTE model atmosphere calculations (e.g., Puls et al. 1996; Sundqvist et al. 2011), and using an input photospheric H line-profile as lower boundary condition (see Sundqvist et al. 2012, for more details).
The rotational phase variation of the emission can now be used to derive the magnetic geometry of the oblique rotator, and the absolute level of H emission further constrains the rate by which the magnetosphere is fed by radiatively driven wind plasma (see also Wade et al., 2015); this latter diagnostic is directly evident through the scaling in equation (27), and is analogous to how H emission from non-magnetic O-stars constrain the wind mass-loss rate. Moreover, since the smoothing length parameter affects the computation of line optical depth quite differently depending on viewing angle, the level of emission contrast between maximum and minimum phases (’high’ and ’low’ states) can, in principle, also be used to constrain this ADM parameter.
As an explicit example of the potential diagnostic power of our model, Figure 7 shows quite remarkably good fits to the observed H light-curve and dynamic spectra (Howarth et al., 2007) of the slowly rotating magnetic O-star HD191612 (Wade et al., 2011; Sundqvist et al., 2012). We use the stellar parameters derived for HD191612 by Howarth et al. (2007) (effective escape speed km/s and stellar radius ), and a maximum loop-closure radius set additionally by the observed dipole magnetic field strength (Howarth et al., 2007; Wade et al., 2011), for wind confinement-parameter (ud-Doula & Owocki, 2002). Since depends only weakly on mass loss and wind terminal speed (to the 1/4 power) and all other parameters here are known from observations, we may safely keep the closure radius fixed during the fitting.
Matching the level and variability contrast of the
observed H line profiles yields here
This is approximately a factor of three higher than the mass-loss rate
used by Sundqvist et al. (2012) to model the H rotational phase variation
in HD191612 directly from MHD simulations, and reflects the fact
that the steady-state ADM model does not account for density-squared
enhancements produced by the highly clumped streams of infalling material
(see discussion above)
In order to match the observed line widths, we have further convolved the dynamic synthetic ADM spectra presented in Figure 7 by a 150 km/s isotropic Gaussian ’macro-turbulence’. While it is not surprising the steady-state ADM models show too little velocity dispersion, we note that such extra broadening is actually required also when modeling H directly from MHD simulations (Sundqvist et al., 2012; ud-Doula et al., 2013).
While the analysis here shows very good fits of the ADM model to the hydrogen H- line in HD191612, the key aim of this first study has been to demonstrate the diagnostic potential of the model, rather than to obtain perfect estimates of all stellar, wind, and magnetic parameters. As discussed further in Appendix C, future detailed parameter-studies of different regimes will be required to fully evaluate, e.g., the accuracy of ADM-derived .
5 Application to Observational Diagnostics: X-rays
5.1 Hot post-shock gas and X-ray emission
The results in section 2.2.2 for the temperature and density of the hot post-shock gas provide a basis for computing the X-ray emission from such DM’s. Following Appendix B of paper I, for a given gas temperature let us write the spectrally integrated (mass-weighted) emission function above a specified X-ray energy as
where the approximation expresses this in terms of the total (mass-scaled) cooling function and a simple “Boltzmann” factor in the ratio . In the context of the present ADM model, let us characterize the threshold energy in terms of its ratio to the maximum shock energy . The associated spectrally integrated volume emissivity in the hot, post-shock region is then given by
For the full MHD simulations of paper I, an analogous simplified form (30) was found to reproduce quite well the spatial distribution of X-ray emission computed from a full integration of the atomic emissivity above the given threshold energy . Appendices A and B here provide a further analysis of how the ADM model can be used to derive both the differential emission measure (DEM), as well as a shock-temperature distribution .
5.2 X-ray emission in ADM vs. MHD simulations
Let us next make a direct comparison of the spatial distribution of the X-ray emission in ADM model vs. that found in MHD simulations. For the same 3 cooling parameter values ( 0.1, 1, and 10) used for the hot-gas temperature and density plots in figure 3, the top row of figure 8 plots the associated variation of the X-ray emissivity , scaled here by . The bottom row compares results for the MHD simulations, showing now the time-averaged Boltzmann corrected emission (divided by , to give results in units of a density squared) for gas above a threshold temperature of MK, corresponding to an X-ray threshold energy 0.13 keV. For the terminal speed km/s of the associated non-magnetic wind, the terminal shock energy keV then implies a threshold-energy ratio , which is thus the value used in the corresponding ADM models. As in figure 5, the MHD simulation output here has been north-south symmetrized to provide clearer comparison with the symmetric ADM model.
Paper I showed that the volume integrated X-emission from MHD simulations can
Indeed one can identify two distinct X-ray emitting regions, with distinct physical origins. The X-rays from inner loops arise from sporadic intervals of “siphon” flow between loop footpoints. This effect is still poorly understood (Bard & Townsend, 2015), but it makes only a minor overall contribution to the total X-ray emission, and is not included in the ADM model.
The second, outer-loop component arises more directly from the collisional shock and retreat that is characterized by the hot post-shock gas in the ADM analysis. The time-variable ‘sloshing’ of hot post-shock gas in the MHD simulations makes its associated time-averaged X-ray emission much more extended, but its radial onset in the MHD simulations corresponds quite well with the inner edge of the shock retreat in the ADM model.
The upshot is that, while the ADM model exaggerates the equatorial concentration of the latitudinal distribution of X-ray emission, it provides a good general description of both its overall spatially integrated value (paper I), and it radial distribution and extension away from the stellar surface.
5.3 X-ray absorption
While the shock-heated X-ray emitting plasma is expected to be mostly optically thin to its own radiation, the cool components of the wind and magnetosphere may very well be optically thick. The continuum opacity in this modestly ionized plasma is due primarily to inner-shell photoionization (bound-free opacity) of metals, and is wavelength-dependent, with generally higher opacities at longer wavelengths. The associated attenuation of the X-rays may lead to potentially observable effects that should have diagnostic power, especially in phase-dependent spectra.
In single O stars with embedded wind shock X-rays, signatures of wind absorption are seen in stars with higher mass-loss rates (yr) via the distortion of X-ray emission line shapes (Cohen et al., 2014). Broadband signatures of X-ray absorption, hardening the emergent spectrum, are also seen in both high-resolution grating spectra and lower resolution CCD spectra of single O stars (Leutenegger et al., 2010; Cohen et al., 2011). Similar X-ray absorption by the fast wind in the open field regions of magnetic O stars should be expected, while the contribution from the magnetosphere, with its steady-state upflow, shock, and downflow cycle, should be comparable to that from a spherical unmagnetized wind. In fact, because of the slower velocity and higher density of the confined wind in the magnetosphere, the degree of X-ray attenuation in the cooled downflow component of the ADM is potentially large.
The bound-free X-ray opacity in the wind and cool magnetospheric plasma is expected to more or less monotonically increase with wavelength through most of the X-ray bandpass, as shown in, e.g., figure 2 of Cohen et al. (2014). This is because for each ion that contributes to the bound-free opacity, the opacity is strongest closest to the threshold set by the ionization potential and decreases strongly toward higher energies. The opacity at any wavelength is the sum of the contributions from all abundant ions. The biggest uncertainty and potential cause of star-to-star variation in the X-ray opacity is likely due to helium, which may be singly ionized in some cases – and would thus contribute significant opacity at longer wavelengths – and may be fully ionized – and thus contribute no bound-free opacity – in other cases. The size of this helium ionization effect can be seen in figure 3 of Cohen et al. (2014). Note that the wind opacity at a fiducial photon energy of 1 keV (12 Å) corresponds to a cross section per hydrogen atom of roughly cm.
We can characterize the overall optical depth of a given ADM model by scaling the column density to the quantity . Thus a given mass-loss rate and a given opacity at a particular wavelength corresponds to a particular value, and the optical depth is degenerate in these two parameters, being proportional to their product. In figure 9 we show optical depth maps for three different values, each computed for two different viewing geometries: one in the magnetic equator and one over the magnetic pole. For the largest characteristic optical depth value shown, =1, which is expected for longer observed wavelengths in a star with a theoretical, non-magnetic mass-loss rate of less than yr, there is significant magnetospheric absorption of X-rays in even the front hemisphere in both the edge-on and pole-on views. This is in addition to occultation by the star itself, which will only be relevant for the edge-on view if the X-ray emitting plasma is in the magnetic equator. The figure shows our standard model with . A larger closed magnetosphere should produce even more attenuation. Although there will be some variation depending on the assumed location of the X-ray emitting plasma. As discussed in the previous section, the ADM models show a concentration of the weighted X-ray emissivity in the magnetic equatorial plane and near , with very little dependence on , while the MHD simulations show a somewhat more complex situation (Figure 8).
From an observational point of view, no significant X-ray absorption is seen in the phase-resolved Chandra grating observations of the prototype magnetic O star, Ori C, which has (Gagné et al., 2005). Numerical MHD simulations show magnetospheric column densities of order cm, consistent with negligible attenuation (ud-Doula et al., 2013; Petit et al., 2015, figure 10). On the other hand, significant X-ray absorption is detected in the low-resolution Chandra spectra of NGC 1624-2 (Petit et al., 2015), which is the O star with the strongest magnetic field and the largest magnetosphere ( vs. for Ori C). More X-ray absorption is seen in NGC 1624-2 when it is observed edge-on than when it is observed at a nearly pole-on phase, a viewing angle modulation that the ADM model is well suited to model. Such a large magnetosphere, with its very strong surface field, is prohibitively expensive to model using numerical MHD.
6 Summary and Future Outlook
The ADM analysis here provides readily computable formulae for the basic hydrodynamic quantities – density, temperature, and velocity – for each of the three components of a wind-fed dynamical magnetosphere – wind upflow, hot post-shock gas, and cooled downflow. Comparison with time-averaged values derived from detailed MHD simulations show, with some caveats, quite good general agreement. As such, this ADM formalism can provide a conceptually and computationally much simpler basis for synthesizing observational diagnostics, and for deriving broad scaling relations for how these depend on stellar, wind, and magnetic parameters.
For X-ray spectral bands, §5.2 compares directly the X-ray emission in ADM vs. MHD models, while §5.3 discusses how observed X-ray spectra could be affected by absorption from the cool components. Appendices A and B present a general scaling analysis of how the ADM model can be used to derive both the differential emission measure, as well as a shock-temperature distribution . This augments the XADM analysis of paper I, and so builds on the promising agreement of the derived scaling laws with observations (ud-Doula et al., 2014; Nazé et al., 2014).
To illustrate ADM spectral synthesis in optical emission lines, we exploit the relative simplicity of the H line formation process in O-type stars. § 4 explicitly demonstrates the potential diagnostic power of the ADM model by a first successful application to observations of the rotational phase variation of H emission from the O-star HD191612. The remarkably good agreement provides constraints on key physical parameters like magnetic geometry and overall mass loss rate , thus illustrating the utility of the ADM formalism even for modeling individual stars with DM’s. The analysis in Appendix C provides also general scaling relations with stellar, wind, and magnetospheric parameters.
But the simple, steady-state nature of the ADM model paves the way for future applications that require more elaborate NLTE radiative transfer. For example, in magnetic O-stars optical helium lines like HeII 4686 show clear signatures of being formed in a DM (Grunhut et al., 2012; Wade et al., 2015), and recent observations of magnetic massive stars in the infra-red (e.g., Oksala et al., 2015) suggest a strong influence of the magnetosphere also for key diagnostics in that waveband. Moreover, the physical explanation of the so-called Of?p morphological phenomena (Walborn, 1972) of magnetic O-stars is very likely related to a complex formation scenario of nitrogen and carbon spectral lines in a DM.
In summary, much as complex computer codes like cmfgen (Hillier & Miller, 1998) or fastwind (Puls et al., 2005) nowadays are routinely applied for spectroscopic analyses of non-magnetic hot stars with winds, we envision that the ADM model presented here lays the groundwork for development of NLTE radiative transfer tools for detailed spectroscopic analysis of magnetic massive stars.
This work was supported in part by SAO Chandra grant TM3-14001A and NASA Astrophysics Theory Program grant NNX11AC40G, awarded to the University of Delaware. AuD acknowledges support by NASA through Chandra Award number TM4-15001A and 16200111 issued by the Chandra X-ray Observatory Center which is operated by the Smithsonian Astrophysical Observatory for and behalf of NASA under contract NAS8- 03060. AuD also acknowledges support for Program number HST-GO-13629.008-A provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. JOS acknowledges funding from the European UnionÕs Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 656725. DHC acknowledges support of Chandra grant TM4-15001B to Swarthmore College. RHDT acknowledges NASA ATP Grant NNX08AT36H to the University of Wisconsin. VP acknowledges partial support of NNX15AG33G from NASA’s XMM-Newton Guest Observer Facility, and NASA Chandra grants TM4-15001C and GO3?14017A. We thank D. Kee and G. Wade for many fruitful discussions. We thank the anonymous referee for many constructive criticisms and comments that helped improve the final version of this paper.
Appendix A Cooling of post-shock material along fixed loop line.
For plasma of density at temperature , we can approximate the local volume emissivity of X-rays above a given energy as,
Here is the integrated plasma emission function above energy , plotted in figure A2 of paper I; as illustratred by the corresponding dashed curves in that figure, the latter approximation expresses this in terms of the total plasma emission, , times a Boltzmann factor.
Let us next integrate this over a volume trace along a fixed, closed magnetic loop line with surface footpoint at a colatitude set by , and with a surface field-line-projection . In terms of the associated local area of the magnetic flux tube along the field line coordinate , the contribution to X-ray luminosity per colatitude interval is
The second and third forms use the energy equation (11) and the field line mass flux equation (6), while the last equality uses the fraction of shock energy emitted above the threshold , as given by eqns. (35)-(37) of paper I.
This result recovers the basic XADM scaling derived in section 5 of paper I; integration over colatitude gives equation (39) there
The full ADM model here now specifies the spatial distribution of the X-ray volume emissivity, through the hot component scaling given in equation (30).
Appendix B Differential Emission Measure
Let us next consider the X-ray differential emission measure (DEM) resulting from the hot component of this ADM model. In terms of an emission measure in electron and proton number density within a volume element , the DEM contribution from a given field line can be written for a differential segment along the field line flux tube with area ,
Using the energy equation (11), this can be cast in the form,
The qualifier emphasizes that this only applies to temperatures up to the shock temperature for this field line; for , the DEM is zero. Conversely, note that for a given temperature , any field line with contributes to the global DEM at that temperature. This proves very useful for deriving the global DEM in the analysis below.
To proceed, note that, in the ADM model, both and the constant flow tube mass flux depend on the footpoint co-latitude of the field line, as set by . For a differential latitude interval , the projection of surface radial mass flux along the field implies , which scales with the surface dipole field radial projection as (Owocki & ud-Doula, 2004),
where is the standard (CAK) mass loss rate for a corresponding non-magnetic star, and the factor half accounts for the equal split of the mass loss in the two hemispheres. Let us then define the cumulative mass loss in a band about the equator (),
Note in particular that for a star with a dipole field that is strong enough to retain its dipole form at the surface for all latitudes, the total surface mass flux is a factor less than a corresponding non-magnetic, spherical wind.
For a given cooling parameter in the ADM model, the shock temperature is a monotonically increasing function of , i.e. , from which we can derive a corresponding inverse function . But for an ADM with a maximum closure radius set by the magnetic confinement parameter , there is a corresponding closure latitude set by , with corresponding maximum scaled pre-shock wind speed . This also sets a maximum post-shock temperature , written here in terms of the terminal speed shock temperature .
In these terms, we can now write the cumulative fraction of total mass flux that has a shock with temperature above a given threshold as
Now, according to (37) each field line contributing to this mass fraction contributes a corresponding amount to the total DEM for any temperature below the shock temperature, i.e., for . This implies that the total DEM from the entire ADM is given in terms of (37) by just multiplying by ,
Apart from some minor differences in notation and overall scaling, this is essentially equivalent to the result given in equation (8) of Gayley (2014), relating DEM to the shock fraction for the general case of radiatively cooled shocks. Moreover, for embedded wind shocks arising from the line deshadowing instability of radiative driving, Cohen et al. (2014) has recently presented an analysis of X-ray line emission that infers a power-law form for the cumulative distribution for the mass fraction undergoing a shock with temperature .
The top panel of figure 10 shows a log-log plot of vs. for cooling parameters 0.1, 1, and 10 (dotted, full, dashed curves), and closure radii characterized by the maximum scaled wind speed , ranging in steps of -0.2 from 1 (representing the limit of arbitrarily strong magnetic confinment) to 0.2 (with a near-surface = 1.25, representing only weak magnetic confinement, with of order unity.)
Note that, for decreasing closure speed , the upper temperature cutoff decreases as . Moreover, for a given, fixed closure speed , the temperature cutoff decreases with increasing , reflecting the stronger shock retreat from less efficient cooling, giving then slower pre-shock flow speeds, , and thus a lower maximum shock temperature, .
The bottom panel of figure 10 shows the corresponding linear-log plots of the scaled differential emission measure. For the simple ADM model here that approximates as constant fixed at a value typical of post-shock temperatures, this is defined by
This assumed constancy of was made to allow analytic solution of the shock retreat from post-shock cooling, and is justified by the fact that most of the total cooling length occurs from the initial cooling from the post-shock temperature, for which the cooling function varies only weakly with temperature. But in developing scaling predictions for an observed DEM, one could also readily apply the ADM derived to a more realistic cooling that includes its full temperature dependence, , via equation (41).
Appendix C H scaling relation
In this equation is the wind electron temperature, the helium number fraction with respect to hydrogen, and the number of free electrons per helium ion. further denotes the kinetic equilibrium (‘NLTE’) departure coefficient of atomic level , for population number density and LTE density with respect to the real population of the ground state of the next higher ion (see Mihalas, 1978).
To derive a characteristic scaling relation for the strength of H emission, let us for now neglect the influence of the photospheric absorption profile and assume an LTE source function set by the radiation temperature of the star, so that the absorption and emission in front of the stellar disc cancel. Using further the fact that the H transition is optically thin in most parts of the wind for O-stars (Puls et al., 1996), we can write down a ‘Sobolev emission measure’ scaling for a clumping-corrected ADM model in terms of integrals over impact parameter and normalized frequency ,
wherein the density and the projected velocity gradient along line of sight direction are to be evaluated at the resonance location for each frequency .
For an observer viewing from a direction along the magnetic pole, the resonance condition will then be close to the magnetic equator, with and . From (25), the cooled downflow density in this region is
Applying this and the velocity gradient scaling in equation (44) gives for the principal scaling of emission measure,