Analytic Dynamical Magnetosphere

An ‘Analytic Dynamical Magnetosphere’ formalism for X-ray 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, 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

1 Introduction

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 typically1 slow rotators, and so harbor DM’s. Among the magnetic B-stars, a significant fraction (about half) are also rotating slowly enough to have DM’s (Petit et al., 2013).

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.

Figure 1: Illustration of three components of material flow along a closed dipole loop line that intersects the stellar radius at a colatitude , and extends up to a maximum radius at the equatorial loop apex. Wind upflow (black arrows) from a stellar surface footpoint meets wind material from the opposite footpoint at the loop apex. This results in a pair of reverse shocks with hot post-shock gas (red) extending back from the loop apex to shock fronts at co-latitudes and , and radius . Cooled material at the loop apex (blue wedge) is pulled back by the stellar gravity, forming channels of cooled downflow (blue dashed arrows) back toward the star.

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 dipole2 magnetic field, taken here to be centered at the central radius of the underlying star.

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 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:

  1. Wind upflow

  2. Hot post-shock gas

  3. 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 non-magnetic case, the mass flux fed into a loop with footpoint scales as (Owocki & ud-Doula, 2004)

(6)

where is the radial projection cosine of the local surface field3, and the second equality applies equation  (4) for a standard dipole.

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 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).

Figure 2: Dipole magnetic field lines (blue curves) along with the retreated shock locations (red curves), labeled with cooling parameters ranging from (closest to magnetic equator) to (furthest from equator) in steps of 1 dex

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.

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 steady-state 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 co-latitude,

(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 strong-shock scaling in modeling the post-shock flow, as discussed next.

Figure 3: Mosaic of ADM hot-component properties. The left, center, and right columns are for cooling parameters , 1, and 10, respectively. The upper row shows hot post-shock temperature, scaled by maximum wind-shock temperature, . The bottom row shows log of the density in the post-shock hot gas (), scaled by the characteristic wind density .
Figure 4: Log of the density of cool material from both the cooled downflow, and the wind upflow , scaled in units of , and assuming . The 3 panels show results for apex smoothing lengths 0, 0.1, and 0.3 (left, center, right).

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),

(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 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,

(12)

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),

(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 post-shock 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 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,

(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, 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

(20)

where is obtained from equations (18) and (16).

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

(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 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. For effective4 stellar mass and radius , and so surface escape speed , conservation of gravitational + kinetic energy gives for the cold gas downflow speed,

(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 co-latitude,

(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 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.

Figure 5: Top row: For an ADM model with maximum closure radius , the top left panel shows log density (scaled by ), given by the wind outflow component for loops with apex radii , and by the cooled downflow component (with ) for loops with . The middle and right panels show the cooled downflow speed in the latitudinal () and radial ( directions, scaled by the stellar surface escape speed . Bottom row: Corresponding time averages for MHD simulations, starting in the lower left with the mean density . To reflect the scaling of line emission, we show and the density-squared-weighted velocity () in the latitudinal (middle) and radial (right) directions. The MHD simulations use parameters that give an Alfvén radius that is approximately equal to the maximum closure radius of the ADM model, . In the MHD results, the appearance of the wind speed in open regions is suppressed by truncating values outside the given colorbar ranges to white. To suppress the impact of stochastic, asymmetric north-south variations in these 2D MHD simulations, the data has been north-south symmetrized to provide a clearer correspondence to the inherent symmetry of the ADM model. In the middle and right lower panels for the MHD simulations, we have set to white any flow that is beyond the quoted colorbar ranges; this has the effect of clearly delineating between open and close field regions, thus allowing a clearer comparison with corresponding ADM results that just show closed-loop infall. The colorbar labels are in CGS units for the MHD simulation; quantitative comparison to the ADM can be done by scaling these with the associated values  g cm for characteristic wind density, and  700 km s.

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 equations (18) and (16);

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 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.

Figure 6: For the MHD model shown in figure 5, spatial variation of time-averaged rms density . Note that the overall form is very similar to the time-averaged density plot in figure 5d, but the colorbar scale here is enhanced by about a factor 7, due to the clumping associated with the complex infall of compressed material.

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

(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 emission-measure scaling law for the ADM model (see Appendix C) for a polar-view observer:

(27)

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

Figure 7: Left: Fit to observed H equivalent width (EW) light-curve (EW in units of Angstrom, with net emission counted positive) of HD191612 (red squares) as function of rotational phase, using the ADM model (black solid line). Right: Normalized flux dynamic spectra of observations (left panel) and simulations (right panel).

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

(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 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 and . 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)5. §3.1. suggests clumping factors of several factors of ten, and adopting here for example would then imply , which is in quite good agreement with predictions from radiatively driven wind theory (Vink et al., 2000). Using these parameters, the magnetic geometry is then also derived from the rotational phase variation and shape of the H equivalent-width curve, resulting here in a degenerate couple , which agrees well with the magnetic geometry derived from spectropolarimetry (Wade et al., 2011).

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

(29)

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

(30)

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 can6 be well modeled with an “XADM” analysis that is grounded in the same basic ADM scalings used here (see eqns. 32-35 of Appendix A). However, figure 8 shows that, while the idealized ADM model predicts a marked concentration of the X-ray emission in the equatorial region around the loop tops, the time-averaged X-ray emission in the MHD simulations is much more spatially extended about the magnetic equator.

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.

Figure 8: Top: For X-rays above a scaled energy , log of ADM emissivity (scaled by ) for 0.1, 1, and 10 (left, middle, right). Bottom: For MHD simulations with the parameters of the standard model in paper I, log of the time-averaged X-ray emissivity above threshold energy  keV, for cooling efficiencies tuned to give , 1, and 10; the middle () case is the same simulation used for figure 5. For the terminal speed  km/s of the associated non-magnetic wind, the terminal shock energy  keV implies the same threshold energy fraction used in the corresponding ADM models. As in figure 5, the MHD simulation data here has been north-south symmetrized to provide clearer comparison with the symmetric ADM model.

5.3 X-ray absorption

Figure 9: Spatial variation of optical depth for bound-free absorption of X-ray emission by both the cool downflow and wind outflow components of the ADM model, as well as by occultation of the opaque star. The top row shows results for a distant observer to the right, with an equator-on view, while the bottom row is for an observer at the top, with a pole-on view. The model assumes an apex smoothing length , and a terminal speed for a corresponding unmagnetized wind. The left, middle and right columns show cases with a corresponding wind optical depth 0.1, 0.3 and 1.

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.

Acknowledgments

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,

(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 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

(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) there7.

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 ,

(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 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),

(38)

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 (),

(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 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

(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 X-ray line emission that infers a power-law form for the cumulative distribution for the mass fraction undergoing a shock with temperature .

Figure 10: Top: Cumulative mass flux fraction yielding a shock temperature above , plotted vs.  on a log-log scale for scaled closure radius wind speeds from (uppermost curves) to (lowermost curves) in steps of . The dotted, full, and dashed lines correspond respectively to cooling parameters 0.1, 1, and 10. Bottom: Associated scaled differential emission measure, defined by DEM, now plotted as linear DEM vs. , for each of the same parameters cited for the top panel.

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

(42)

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

Figure 11: H emission equivalent width (counted positive) in Angstrom vs. . The red triangles show results from full radiative transfer calculations using different values of (see text). The dotted line displays the EW scaling according to eqn. C4, scaled to match the result at the lowest mass loss rate.

From the Saha-Boltzmann 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 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 ,

(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,