An ‘Analytic Dynamical Magnetosphere’ formalism for Xray and optical emission from slowly rotating magnetic massive stars
Abstract
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, steadystate 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 postshock 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 timeaveraged 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 Xray emission from the hot postshock 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: earlytype — Stars: mass loss — Stars: Xrays1 Introduction
Hot luminous, massive stars of spectral type O and B have dense, highspeed, 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 longlived 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 (udDoula
et al., 2009), magnetic Otype 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 Xray 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 (udDoula & Owocki, 2002; udDoula 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 Ostars, using simplified radiative transfer methods to derive synthetic spectra for H emission lines (Sundqvist et al., 2012; Grunhut et al., 2012; udDoula et al., 2013; Wade et al., 2015) and ultraviolet 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 timedependent 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 Xray spectral synthesis. Detailed MHD simulations have been used to analyze the high spectral resolution Xray 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 lowresolution Xray data, an analytic model can capture key observable properties and trends. Recently, udDoula et al. (2014, hereafter paper I) carried out an extensive MHD simulation study of radiative cooling of the hot postshock gas in DM’s, with a focus on deriving the resulting Xray emission as a function of the densitydependent cooling efficiency. When interpreted in terms of a simplified “Xray analytic dynamical magnetosphere” (XADM) analysis, this led to predicted scaling laws for variation of Xray luminosity with the wind mass loss rate. Subsequent application by Nazé et al. (2014) toward interpreting observationally inferred Xray luminosities in a large sample of magnetic massive stars showed that, with some fixed factor adjustment to account for the limited “duty cycle” for Xray 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 postshock gas, and cooled downflow (section 2). The overall ADM analysis here is guided and tested by comparison with timeaveraged 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 boundfree absorption of emitted Xrays 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 Xray differential emission measure (DEM), as well as a shocktemperature distribution . Appendix C gives background on the ADM scalings for H emission.
2 Analytic Dynamical Magnetosphere
2.1 Basic magnetic scalings for a starcentered dipole
The ADM formalism is based on an idealized, steadystate analysis of the mass flow along closed magnetic loops that are assumed to follow the individual lines of a simple
dipole
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
(1) 
and with the maximum radius at the loop apex on the magnetic equator,
(2) 
For any position within the magnetosphere, the magnitude of the field follows the dipole scaling,
(3) 
For such a dipole, the direction of the field only depends on , given by the unit vector
(4) 
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,
(5) 
where the third equality uses , and the last equality defines a coordinate distance along the field line, with . This means that in any steadystate 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 3Component 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 timedependent, 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, steadystate flows.
Specifically, as illustrated in figure 1, this ADM analysis distinguishes three distinct components of material flow within the loop:

Wind upflow

Hot postshock gas

Cooled downflow
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 nonmagnetic case, the mass flux fed into a loop with footpoint scales as (Owocki & udDoula, 2004)
(6) 
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,
(7) 
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 loopspecific 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 postshock gas that extends away from the loop apex by a distance set by the postshock cooling length . For lower luminosity stars with smaller mass loss rate and thus a lowerdensity 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 postshock temperature . Using the dipole shockretreat 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 postshock temperature of order many MK, radiative cooling in the postshock 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 postshock 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 freefall 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 Xray boundfree opacity from abundant, partially ionized heavy elements like CNO and Fe, which can thus play a role in attenuating the Xrays emitted from the hot postshock gas (Petit et al., 2015).
As noted, in 2D and 3D MHD simulation of DM’s (udDoula et al., 2013, paper I), these three components of wind upflow, hot postshock 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 shockbuildup, cooling, and then infall (udDoula 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 largescale spatial distribution. This can be well characterized by a superposition of the density, velocity and temperature of the wind upflow, hot postshock gas, and cooled downflow. The next subsections quantify this in terms of relatively simple ADM scalings based on quasisteadystate conservation applied to the material from each component.
Wind upflow
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),
(8) 
The density then follows from steadystate mass continuity, with base mass flux (6), and accounting for the dipole area divergence,
(9) 
where is a characteristic density for an unmagnetized wind with mass loss rate and terminal speed .
Using equations (1), (3), (6), and (8), we can rewrite this in a form that explicitly shows the dependence on radius and colatitude,
(10) 
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 strongshock scaling in modeling the postshock flow, as discussed next.
Hot postshock gas
Let us next derive conditions in the hot postshock gas in the region between the shock radius and loop apex, . Our analysis builds on and extends the Xray analytic dynamical magnetosphere (XADM) formalism developed in paper I, deriving now explicit expressions for the spatial variation of the temperature in this postshock 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 massweighted 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),
(11) 
where is the mean molecular weight, and the temperature derivative is with respect to the field line coordinate . The nearly constant pressure of this postshock 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 postshock value to the much lower (near zero) temperature at the loop apex,
(12) 
here , and accounts for the mass loss weighting for a given fieldline 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),
(13) 
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
(14) 
where
(15) 
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
(16) 
where, for the assumed case of a highly supersonic outflow yielding a strong shock, the immediate postshock temperature , with
(17) 
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,
(18) 
where here we assume a typical hotstar 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,
(19) 
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, postshock region. For a strong shock, the density of the immediate postshock gas is simply a factor four times the incoming wind density at the shock, . Since the pressure is nearly constant in this subsonic postshock cooling layer, we can write the spatial variation of density of the hot gas as
(20) 
Using mass continuity, we can also readily derive the spatial variation of the postshock flow speed. For a strong shock, the immediate postshock speed is just a quarter of the incoming wind speed, , while the spatial variation is given by
(21) 
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).
Cooled downflow
Once this hot postshock 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.
For effective
(22) 
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 ,
(23) 
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
(24) 
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 colatitude,
(25) 
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 timeaverages 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),
with all 3 using the auxiliary relations (15), for the field line geometry function ), and (19), for shock location, .
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 postshock 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 timeaveraged 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 buildup and dynamical infall. To ensure the northsouth symmetry of an arbitrarily long time average (or from azimuthal averaging of a full 3D model; see figure 3 of udDoula et al. (2013)), we also carry out a northsouth averaging of all the timeaveraged 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 recombinationbased emission lines like H. Section 5 next focuses on Xray emission from the hot postshock 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 timeaveraged 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 densitysquared 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 semiquantitive, 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 windfed 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 steadystate 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 densitysquared, which for the twobody collision or recombination processes that underlie line emission is what sets the overall emission strength. For this MHD model, figure 6 plots now the timeaveraged 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 Ostar 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 Ostars (Sundqvist et al., 2012).
4.1 ADMmodified scaling relation for H emission
The principal scaling of H emission in OBstar winds comes from considering the line optical depth in the Sobolev approximation. For Doppler width and directional Sobolev length , this is
(26) 
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 emissionmeasure scaling law for the ADM model (see Appendix C) for a polarview observer:
(27) 
where the function describes the dependence on the size of the closedloop magnetosphere. Equation (27) illustrates explicitly how the standard scaling for nonmagnetic wind emission is modified here by the two magnetic parameters and (essentially setting the ADM discthickness 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 OBstar 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 OBstars. To compute synthetic spectra from the steadystate 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 lineofsight make with the rotation axis, the variation of with rotational phase is given by
(28) 
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 lineprofile 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 nonmagnetic Ostars constrain the wind massloss 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 lightcurve and dynamic spectra (Howarth et al., 2007) of the slowly rotating magnetic Ostar 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 loopclosure radius set additionally by the observed dipole magnetic field strength (Howarth et al., 2007; Wade et al., 2011), for wind confinementparameter (udDoula & 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
and .
This is approximately a factor of three higher than the massloss 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 steadystate ADM model does not account for densitysquared
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 ’macroturbulence’. While it is not surprising the steadystate 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; udDoula 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 parameterstudies of different regimes will be required to fully evaluate, e.g., the accuracy of ADMderived .
5 Application to Observational Diagnostics: Xrays
5.1 Hot postshock gas and Xray emission
The results in section 2.2.2 for the temperature and density of the hot postshock gas provide a basis for computing the Xray emission from such DM’s. Following Appendix B of paper I, for a given gas temperature let us write the spectrally integrated (massweighted) emission function above a specified Xray energy as
(29) 
where the approximation expresses this in terms of the total (massscaled) 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, postshock region is then given by
(30) 
For the full MHD simulations of paper I, an analogous simplified form (30) was found to reproduce quite well the spatial distribution of Xray 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 shocktemperature distribution .
5.2 Xray emission in ADM vs. MHD simulations
Let us next make a direct comparison of the spatial distribution of the Xray 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 hotgas temperature and density plots in figure 3, the top row of figure 8 plots the associated variation of the Xray emissivity , scaled here by . The bottom row compares results for the MHD simulations, showing now the timeaveraged 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 Xray threshold energy 0.13 keV. For the terminal speed km/s of the associated nonmagnetic wind, the terminal shock energy keV then implies a thresholdenergy ratio , which is thus the value used in the corresponding ADM models. As in figure 5, the MHD simulation output here has been northsouth symmetrized to provide clearer comparison with the symmetric ADM model.
Paper I showed that the volume integrated Xemission from MHD simulations can
Indeed one can identify two distinct Xray emitting regions, with distinct physical origins. The Xrays 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 Xray emission, and is not included in the ADM model.
The second, outerloop component arises more directly from the collisional shock and retreat that is characterized by the hot postshock gas in the ADM analysis. The timevariable ‘sloshing’ of hot postshock gas in the MHD simulations makes its associated timeaveraged Xray 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 Xray 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 Xray absorption
While the shockheated Xray 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 innershell photoionization (boundfree opacity) of metals, and is wavelengthdependent, with generally higher opacities at longer wavelengths. The associated attenuation of the Xrays may lead to potentially observable effects that should have diagnostic power, especially in phasedependent spectra.
In single O stars with embedded wind shock Xrays, signatures of wind absorption are seen in stars with higher massloss rates (yr) via the distortion of Xray emission line shapes (Cohen et al., 2014). Broadband signatures of Xray absorption, hardening the emergent spectrum, are also seen in both highresolution grating spectra and lower resolution CCD spectra of single O stars (Leutenegger et al., 2010; Cohen et al., 2011). Similar Xray 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 steadystate 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 Xray attenuation in the cooled downflow component of the ADM is potentially large.
The boundfree Xray opacity in the wind and cool magnetospheric plasma is expected to more or less monotonically increase with wavelength through most of the Xray bandpass, as shown in, e.g., figure 2 of Cohen et al. (2014). This is because for each ion that contributes to the boundfree 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 startostar variation in the Xray 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 boundfree 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 massloss 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, nonmagnetic massloss rate of less than yr, there is significant magnetospheric absorption of Xrays in even the front hemisphere in both the edgeon and poleon views. This is in addition to occultation by the star itself, which will only be relevant for the edgeon view if the Xray 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 Xray emitting plasma. As discussed in the previous section, the ADM models show a concentration of the weighted Xray 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 Xray absorption is seen in the phaseresolved 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 (udDoula et al., 2013; Petit et al., 2015, figure 10). On the other hand, significant Xray absorption is detected in the lowresolution Chandra spectra of NGC 16242 (Petit et al., 2015), which is the O star with the strongest magnetic field and the largest magnetosphere ( vs. for Ori C). More Xray absorption is seen in NGC 16242 when it is observed edgeon than when it is observed at a nearly poleon 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 windfed dynamical magnetosphere – wind upflow, hot postshock gas, and cooled downflow. Comparison with timeaveraged 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 Xray spectral bands, §5.2 compares directly the Xray emission in ADM vs. MHD models, while §5.3 discusses how observed Xray 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 shocktemperature distribution . This augments the XADM analysis of paper I, and so builds on the promising agreement of the derived scaling laws with observations (udDoula 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 Otype 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 Ostar 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, steadystate nature of the ADM model paves the way for future applications that require more elaborate NLTE radiative transfer. For example, in magnetic Ostars 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 infrared (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 socalled Of?p morphological phenomena (Walborn, 1972) of magnetic Ostars 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 nonmagnetic 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.
Acknowledgments
This work was supported in part by SAO Chandra grant TM314001A and NASA Astrophysics Theory Program grant NNX11AC40G, awarded to the University of Delaware. AuD acknowledges support by NASA through Chandra Award number TM415001A and 16200111 issued by the Chandra Xray 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 HSTGO13629.008A 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 NAS526555. JOS acknowledges funding from the European UnionÕs Horizon 2020 research and innovation programme under the Marie SklodowskaCurie grant agreement No 656725. DHC acknowledges support of Chandra grant TM415001B to Swarthmore College. RHDT acknowledges NASA ATP Grant NNX08AT36H to the University of Wisconsin. VP acknowledges partial support of NNX15AG33G from NASA’s XMMNewton Guest Observer Facility, and NASA Chandra grants TM415001C 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 postshock material along fixed loop line.
For plasma of density at temperature , we can approximate the local volume emissivity of Xrays above a given energy as,
(31) 
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 fieldlineprojection . In terms of the associated local area of the magnetic flux tube along the field line coordinate , the contribution to Xray luminosity per colatitude interval is
(32)  
(33)  
(34)  
(35) 
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 Xray volume emissivity, through the hot component scaling given in equation (30).
Appendix B Differential Emission Measure
Let us next consider the Xray 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 ,
(36) 
Using the energy equation (11), this can be cast in the form,
(37) 
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 colatitude 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 & udDoula, 2004),
(38) 
where is the standard (CAK) mass loss rate for a corresponding nonmagnetic 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 (),
(39)  
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 nonmagnetic, 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 preshock wind speed . This also sets a maximum postshock 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
(40) 
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 ,
(41) 
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 Xray line emission that infers a powerlaw form for the cumulative distribution for the mass fraction undergoing a shock with temperature .
The top panel of figure 10 shows a loglog 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 nearsurface = 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 preshock flow speeds, , and thus a lower maximum shock temperature, .
The bottom panel of figure 10 shows the corresponding linearlog plots of the scaled differential emission measure. For the simple ADM model here that approximates as constant fixed at a value typical of postshock temperatures, this is defined by
(42)  
This assumed constancy of was made to allow analytic solution of the shock retreat from postshock cooling, and is justified by the fact that most of the total cooling length occurs from the initial cooling from the postshock 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
From the SahaBoltzmann relations and the atomic constants of the hydrogen transition, the parameter in equation (26) of the main paper can be written as (Puls et al., 1996; Petrenz & Puls, 1996)
(43) 
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 Ostars (Puls et al., 1996), we can write down a ‘Sobolev emission measure’ scaling for a clumpingcorrected ADM model in terms of integrals over impact parameter and normalized frequency ,
(44) 
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
(45) 
Applying this and the velocity gradient scaling in equation (44) gives for the principal scaling of emission measure,