Nonlinear evolution of the radiation-driven magneto-acoustic instability (RMI)

Rodrigo Fernández11affiliation: Einstein Fellow and Aristotle Socrates22affiliation: John N. Bahcall Fellow Institute for Advanced Study. Einstein Drive, Princeton, NJ 08540, USA.

We examine the nonlinear development of unstable magnetosonic waves driven by a background radiative flux – the Radiation-Driven Magneto-Acoustic Instability (RMI, a.k.a. the “photon bubble” instability). The RMI may serve as a persistent source of density, radiative flux, and magnetic field fluctuations in stably-stratified, optically-thick media. The conditions for instability are present in a variety of astrophysical environments, and do not require the radiation pressure to dominate or the magnetic field to be strong. Here we numerically study the saturation properties of the RMI, covering three orders of magnitude in the relative strength of radiation, magnetic field, and gas energies. Two-dimensional, time-dependent radiation-MHD simulations of local, stably-stratified domains are conducted with Zeus-MP in the optically-thick, highly-conducting limit. Our results confirm the theoretical expectations of Blaes and Socrates (2003) in that the RMI operates even in gas pressure-dominated environments that are weakly magnetized. The saturation amplitude is a monotonically increasing function of the ratio of radiation to gas pressure. Keeping this ratio constant, we find that the saturation amplitude peaks when the magnetic pressure is comparable to the radiation pressure. We discuss the implications of our results for the dynamics of magnetized stellar envelopes, where the RMI should act as a source of sub-photospheric perturbations.

Subject headings:
radiative transfer — diffusion — MHD — instabilities — methods: numerical

1. Introduction

Many gravitationally-bound astrophysical systems remove their binding energy by the diffusion of radiation through optically thick regions. Examples are main-sequence stars and accretion disks onto compact objects. In general, these systems tend to be good conductors and therefore support magnetic fields.

Blaes & Socrates (2003; hereafter BS03) found that these environments are susceptible to magnetosonic overstability. The physical mechanism responsible for secular driving involves short wavelength compressible fluctuations that grow exponentially due to the presence of a background radiative flux. Fluid motion along the direction of the equilibrium magnetic field, resulting from magnetic tension forces, couples to the radiative flux perturbation. In the event that these two perturbations are in phase, the radiation field performs work on the fluid oscillation and increases its amplitude (Figure 1).

A surprising result of BS03 is that weakly magnetized and/or gas pressure dominated equilibria are subject to the same radiative driving mechanism responsible for overstability in radiation pressure and/or magnetic pressure dominated environments. The instability mechanism is so generic that even when the radiation is degenerate (e.g., neutrinos in core-collapse supernova environments, Socrates et al. 2005), it may still operate. Furthermore, a similar mechanism operates in the effect that the energy-transporting particles are charged and diffuse primarily along magnetic field lines (Socrates et al. 2008).

Before BS03, such phenomena were thought to be restricted to radiation-pressure-dominated media that are strongly magnetized. For these conditions, the instability is often referred to as “photon bubbles” (Prendergast & Spiegel 1973; Arons 1992; Gammie 1998; Blaes & Socrates 2001; Begelman 2001), though neither buoyancy, nucleation, nor surface tension play a role in the driving mechanism. Numerical investigation of the nonlinear development of this instability has been limited, with only a few calculations of strongly-magnetized and mostly radiation pressure dominated atmospheres reported in the literature (Klein et al. 1996; Hsu et al. 1997; Davis et al. 2004; Turner et al. 2004, 2005, 2007; Jiang et al. 2012). A major difficulty for numerical studies results from the fact that the radiation diffusion time over a wavelength smaller than the gas pressure scale height needs to be resolved in order to attain numerical convergence (Turner et al. 2005). In practice, the oscillation period itself is often much shorter than the e-folding time of the mode, requiring a large number of time steps to capture the secular amplification of the waves.

Here we take advantage of the parallel radiation-MHD code Zeus-MP (Hayes et al. 2006) to begin an exploration of the more extended parameter space for unstable modes found by BS03. To reflect the physical character of the driving mechanism, we designate it the Radiation-Driven Magneto-Acoustic Instability (RMI).

The parameter space to be explored is motivated by the radiative envelopes of massive stars, many of which are known to support up to kG-strength magnetic fields (e.g., Wade et al. 2012 and references therein). In main-sequence O stars, the radiation pressure only marginally exceeds the equilibrium gas pressure; a magnetic field near equipartition with the gas provides favorable conditions for the growth of the RMI (e.g., BS03, Turner et al. 2004). Going down the main sequence in mass and luminosity, radiation pressure becomes increasingly less important, yet the instability may persist near the photosphere.

The paper is organized as follows. Section 2 contains an overview of the physics of the RMI. Section 3 describes our assumptions, the numerical method employed, and the parameter space explored. Section 4 presents our results. Section 5 briefly discusses implications for massive stars. We conclude with a summary of our work in Section 6. Appendix A contains a brief derivation of the growth rate of the RMI and the stability criteria. Appendix B describes verification tests of Zeus-MP. Mathematical symbols have their usual meanings; the most frequently used ones are defined in Table 1.

Figure 1.— Driving of magnetosonic modes by a background radiative flux (BS03). The shaded region represents a magnetosonic plane wave, with darker and lighter regions corresponding to overdensities and rarefactions, respectively. The perturbation to the flux does work on the component of the velocity perturbation perpendicular to the wave vector , provided that the latter is not entirely parallel or perpendicular to the magnetic field . The acceleration of gravity is denoted by , and velocity vectors are drawn for a slow mode in the weak field limit.

2. Overview of the RMI

The photon bubble instability was discovered numerically by Arons (1992) and Gammie (1998) in atmospheres that are radiation pressure dominated, with either strong () or super-strong () magnetic fields. In their computations, they found the slow magnetosonic wave to be susceptible to over-stability, driven by the presence of a background radiative flux . Blaes & Socrates (2001) subsequently found that the fast magnetosonic mode could also be driven unstable, and that in the presence of rotation, the splitting between the Alfvén and slow wave allows for the – otherwise incompressible – Alfvén branch to be driven as well. The effects of magnetic stratification in the magnetic- and radiation-dominated regime have been addressed by Tao & Blaes (2011).

The stability criteria and dynamics of the RMI for a constant magnetic field were derived in BS03. They found that (i) stability is determined by a balance between driving from the background flux , damping due to radiative diffusion, and departures from thermal equilibrium, (ii) the instability can be thought of as magnetosonic waves that are secularly driven by the background flux , (iii) the instability mechanism operates independently of the degree of thermal coupling between the gas and radiation, and (iv) the instability mechanism operates even when the magnetic and/or radiation pressure is/are less than the equilibrium gas pressure.

In order to determine the relevant parameter space for numerical simulations of the RMI, we must understand the physics that determines whether a given atmosphere is unstable to RMI driving. A brief derivation is found in Appendix A. Below we provide a general overview.

Symbol Quantity
gas energy density
radiation energy density
entropy per baryon
baryon number density
gas pressure
isothermal gas sound speed
gravitational acceleration
radiation pressure
radiative flux
flux-mean opacity
absorption opacity
magnetic field
Alfvén velocity
magnetic pressure
wave vector
mode frequency
diffusion frequency (eq. [2])
heat exchange frequency (eq. [4])
rapid-diffusion parameter (eq. [13])
critical Eddington ratio (eq. [14])
speed of light
box size
adiabatic gas sound speed at
density at
Table 1Frequently-Used Symbols

2.1. Characteristic Frequencies

We consider here the evolution of plane waves of frequency and wave vector on an ideal MHD fluid with constant background magnetic field and gravitational acceleration . Furthermore, we take the various frequency-weighted opacities to be constant in order to isolate the RMI from opacity-driven instabilities.

All of the unstable RMI parameter regimes found by BS03 involve heat rapidly diffusing across a wavelength in comparison to the frequency of oscillation. That is




is the radiative diffusion frequency, and is the flux-mean opacity. Radiation diffuses rapidly in comparison to the mode oscillation time, i.e., the perturbations are highly non-adiabatic.

In this paper, we restrict ourselves to the case where the gas and radiation are thermally well coupled, so that both are described by a single temperature. This condition is satisfied in the radiative envelopes of massive stars. Quantitatively, we demand that




is a characteristic heat exchange frequency (BS03), with the thermal absorption (Planck-mean) opacity, and the speed of light. In this limit, gas-acoustic perturbations travel at the isothermal 111This is not the same as using an isothermal equation of state, since is a local quantity. sound speed, .

2.2. Stability Criteria

As shown in Appendix A, the criterion for RMI driving is given by (eqns. [A.2]-[A21])

The expression above quantifies the competition between the rate of radiative energy transfer into a fluctuation over its wavelength due to the presence of a background flux , and the rapid diffusion of energy that results from compression. Equation (2.2) is obtained as a first order expansion in the small parameters , , and , where is the gas pressure scale height. Equation (2.2) differs from equation (A.2) only by angular factors, being thus accurate to within a factor of order unity. A more accurate value than equation (A.2) can be obtained by solving the eighth-order dispersion relation of BS03.

For slow magnetosonic waves, we may write


where we have used equation (A16) and defined


The growth rate can be approximated by




is the Eddington ratio. The angular dependence of equation (A.2) has been replaced by a global factor of that accounts for the product of the sines of the relative angles between , , and .

For fast magnetosonic waves, the corresponding instability criterion reads


with a corresponding asymptotic growth rate


The extra factor of has also been included.

It is instructive to visualize the instability criteria for slow and fast waves in a plane with the Eddington ratio on the x-axis and the ratio on the y-axis, as is shown in Figure 2. The critical stability curves can be obtained by relating the Eddington ratio to the ratio of radiation to gas pressure through




is the ratio of gas pressure scale height to radiation pressure scale height .

The resulting curves have an additional free parameter,


where is the optical depth over a gas pressure scale height. The parameter is proportional to the ratio of diffusion speed to the isothermal sound speed. When , there is a critical value of the Eddington ratio


above which fast and slow waves are unstable for any magnetic field strength. If , then only slow modes can be destabilized, and for every Eddington ratio there is a maximum field strength for which the slow wave RMI operates.

Figure 2.— Unstable parameter regions in the thermally locked case (mode frequencies smaller than eq. [4]). If the “rapid-diffusion” parameter (eq. [13]) exceeds unity, then a critical value of the Eddington ratio (eq. [14]) exists. For , unstable fast-magnetosonic modes are possible, and slow-magnetosonic modes are unstable for any magnetic field. Below , unstable slow modes are possible only up to a maximum value of given by equation (5), and fast modes are stable. If , then for any Eddington ratio there exists a maximum value of for which slow modes can be destabilized.

2.3. Aspects of the Stability Criteria

The growth/damping rate and stability criteria above give insight regarding the properties of the equilibrium about which the RMI takes place. Balance between radiative driving and damping is, respectively, a balance between energy input from the background radiative flux and energy loss due to diffusive cooling that results from compression, which is .

Consider, for example, the slow magnetosonic wave in the strong field limit, so that equation (2.2) becomes


where is the wavelength of an oscillation with wave number , and is the oscillation period. The left hand side is the upper limit – if, e.g., the perturbations were nonlinear – for the rate of radiative energy transfer into a fluctuation of scale . The characteristic velocity at which the background radiative field transfers radiative energy from one layer to another is given by , and the corresponding time at which radiative diffusion transfers energy across the wavelength is given by . The right hand side of the criteria for instability results from entropy or pressure perturbations that are sourced by changes in density. Since radiation diffuses, effectively, instantaneously over a wavelength during the oscillation time , the rate of energy loss per unit volume is capped by the relative contribution to pressure or entropy changes that arise from compression, which is and the oscillation time .

By re-inserting (eq. [A22]), the instability criteria can be written as


where the phase velocity (eq. [A4]), and we have taken . The maximum value of the left hand side is achieved at the photosphere where .

Equation (16) re-enforces several facts that are important in attempting to understand where RMI phenomena take place. Instability vanishes in limit since the ratio of the radiation flux to energy density vanishes as well. The free energy for the instability is rooted in the anisotropy of the radiation field, which leads to a radiative flux that can then perform work on fluid fluctuations. Therefore, in the center or stars (or, e.g., at the surface of last scattering for the cosmic background radiation) RMI driving does not operate. Consequently, instability favors regions relatively near the photosphere.

The presence of in the denominator of the left hand side implies that the instability favors magneto-acoustic oscillations with small phase velocity. Given the oscillation wavelength, the rate at which radiation removes energy from a fluid disturbance is given by the phase velocity. That is, the rate at which radiative diffusion can remove energy from the wave is relatively small for relatively slow oscillations.

3. Numerical Setup

Our aim is to follow the nonlinear development of the RMI in a small patch of stellar interior. Below we describe the equations solved, numerical method, initial conditions, boundary conditions, and choice of parameters for simulations.

3.1. Equations

We solve the radiation-magnetohydrodynamic equations in the ideal MHD and grey diffusion limits


Our goal is to isolate the growth and nonlinear development of the RMI, which serves as justification for setting the flux-mean opacity to a constant, as in the case of the Thomson scattering opacity. Radiation-hydrodynamic instabilities such as the mechanism and the so-called strange mode instability (e.g., Glatzel 1994), which rely on opacity variations, are thus excluded.

The Planck-mean and intensity-mean opacities normally enter independently into the radiation and gas energy equations (eqns. [19]-[20]). Here they are set equal to one another, implicitly assuming that the radiation follows a Planck distribution with its own temperature (e.g., BS03). We refer to these opacities collectively as the absorption opacity . In our calculations, the ratio is fixed to a value such that the condition for thermal locking (eq. [3]) is satisfied.

Equations (17)-(22) are closed by an equation of state for the matter fluid, which we take to be an ideal gas of adiabatic index . It follows that , where and are the mean molecular weight and proton mass, respectively. In the diffusion approximation, the radiation stress tensor is related to the radiation energy density by , where is the identity tensor. We denote the scalar radiation pressure by . No flux-limiters are employed (eq. [21]).

A cartesian coordinate system is adopted, with the acceleration of gravity pointing in the negative direction. Throughout this study, we restrict ourselves to two-dimensional simulations, with denoting the coordinate direction transverse to gravity.

3.2. Numerical Method

To evolve the system of equations (17)-(22), we employ the publicly available finite-difference code Zeus-MP (Hayes et al. 2006). The MHD part is evolved using the Modified Characteristics – Constrained Transport method of Hawley & Stone (1995). Shocks are captured with the artificial viscosity prescription of VonNeumann & Richtmyer (1950). The radiation part is evolved implicitly via the Conjugate Gradient method (Hayes et al. 2006), and interacts with the MHD sector via operator splitting.

As a basic test of our implementation, we verify that the initial conditions (§3.3) do not evolve when unperturbed. We have also tested the elimination of the flux limiter in Zeus-MP by solving the diffusion equation in a unit square with periodic boundary conditions and constant initial . A more stringent test involves comparing eigenfrequencies of magnetosonic eigenmodes on a uniform background in the thermally-locked limit with the asymptotic values found by BS03. Details on both of these tests are provided in Appendix B.

aaValues of , , , and at the center of the box. bbThe parameters and are determined by requiring that is 100 times the slow magnetosonic frequency of a mode with wavelength equal to , at domain center. bbThe parameters and are determined by requiring that is 100 times the slow magnetosonic frequency of a mode with wavelength equal to , at domain center. ccWidth of the computational box in the direction transverse to gravity, in units of at . The vertical size is always . ddAngle between the imposed magnetic field and the background radiative flux. eeResolution in the directions transverse and parallel to gravity, respectively. ffRatio of the root-mean-square fluctuation over the mean (eqns. [35]-[36]) at positions , respectively.
0.01 0.3 31 1.2 2
0.03 3.5
0.1 11
0.3 28
1 60
3 91
10 ggThe diffusion solver breaks down before saturation, due to order unity density fluctuations. 110
1 0.01 5.7 9.8 2
0.1 18 34
1 40 79
10 40 79
0.1 0.3 31 11 2
0.1 0.3 31 11 4
1 0.3 31 60 2
Table 2Models Evolved and Time-Averaged Properties

3.3. Equilibrium Structure and Initial Conditions

We initialize the computational domain with an atmosphere in hydrostatic, radiative, and thermal equilibrium. The equations describing this atmosphere are obtained from equations (18)-(21) and the gas equation of state, by setting , , and dropping the time derivatives:


The equilibrium magnetic field is taken to be uniform in space and therefore does not contribute to the equilibrium structure of the atmosphere.

The problem is completely specified by nine parameters at one location. These are the ratio of box size to the gas pressure scale height ; the Eddington ratio ; the ratio of radiation to gas pressure ; the adiabatic index ; the mean molecular weight , the ratio of isothermal sound speed to light speed ; the optical depth over a gas pressure scale height ; the ratio of magnetic to gas pressure ; and the angle between magnetic field and radiative flux .

We initialize the domain in such a way that there is cancellation to within a few parts in between all the forces. This is possible because Zeus-MP computes gradients using finite differences on a staggered mesh. Evolving this initial state without additional perturbations excites vertical fast magnetosonic modes, which decay (there is no RMI driving if the wave-vector is perpendicular to the magnetic field and/or parallel to the radiative flux).

Random initial perturbations in density are applied from to for all . The amplitude is , and the random number is varied from cell to cell. The internal energy density is also perturbed following the density, , so that the gas temperature and radiation energy density remain uniform everywhere.

For simplicity, we consider atmospheres with , which are stably-stratified in the rapid diffusion limit (BS03).

3.4. Boundary Conditions

At the inner vertical boundary, we fix the density, gas-, and radiation energy density to their steady-state values, by continuing the initial solution beyond the boundary. The radial velocity is reflected at this boundary, while the tangential velocity is set to zero in the ghost cells. The magnetic field is set to its uniform initial value. This arrangement ensures that all waves are reflected at the inner boundary, and that the domain does not collapse under the action of gravity. A constant radiative flux is thus incident on the box from below.

At the outer boundary, we use an outflow boundary condition that allows waves to leave the domain but avoids runaway mass loss. To account for the fact that this is an internal patch inside a hydrostatic structure rather than a domain with a ‘vacuum’ external boundary, we set the outflow radial velocity so that the mass flux leaving the domain is that due do the excess density over the steady-state value. At a given transverse coordinate, we thus set


where the subscripts and label values in the initial and time-dependent solution at the last active cell inside the outer radial boundary. This prescription automatically includes the case where the density is much larger than the steady-state value, smoothly approaching a pure outflow boundary condition for the radial velocity. The density, gas energy density, magnetic field, and tangential velocity from the last active zone are copied into the two ghost zones. If the radiative flux between the last pair of active zones is positive, the radiation energy density in the ghost zones is set so as to keep this flux constant across cells. If instead the flux is zero or negative, then the radiation energy density is set to have zero gradient, yielding zero flux.

The horizontal boundary condition is periodic for all variables.

3.5. Models Evolved

We run two base sequences of models that vary the ratios of radiation to gas pressure and magnetic to gas pressure . These base sequences employ square computational boxes of two gas pressure scaleheights on a side , with a resolution of cells, and a time step equal to the diffusion time for a wavelength . The magnetic field points in the direction perpendicular to gravity . This choice of field geometry allows waves to cross the domain many times while their amplitude grows. We perform a few runs that explore the effects of changing the box size, resolution, and field orientation. All models are shown in Table 2, with parameters showing values at the center of the domain.

In all simulations, we take , , and . The optical depth in the box is chosen so that the diffusion frequency is times larger than the slow magnetosonic frequency for a wavelength equal to . The wave number dependence of ensures that all modes of shorter wavelength are well within the rapid-diffusion limit. This choice fixes the value of and thus .

The ratio is chosen so that is times larger than the slow magnetosonic mode of wavelength equal to . Our choice is motivated by results from tests of radiation-MHD waves on a uniform background (Appendix B).

4. Results

4.1. Linear Growth and Saturation

The general behavior of all models involves damping of the initially applied perturbations, and the subsequent emergence of unstable propagating waves. Because the background magnetic field is perpendicular to gravity in most models, waves are allowed to cross the domain many times222To lowest order, the group velocity of slow magnetosonic modes lies along the magnetic field.. This leads to exponential growth of the kinetic energy.

The evolution of the volume-integrated kinetic energy for a gas-dominated model (, , ) is shown in Figure 3a, broken up into vertical and horizontal components,


respectively. Initial density perturbations drive vertically-propagating waves that eventually decay. The horizontal kinetic energy shows no such initial transients, growing out of numerical noise from the beginning of the simulation. Both energies saturate at nearly same time, closely following each other after the RMI overtakes the decaying waves.

Figure 3.— Panel (a): Evolution of the volume-integrated vertical and horizontal kinetic energies (eqns. [27]-[28]) as a function of time, for the model with and (, resolution ). Panel (b): rate of change of the horizontal kinetic energy in the same model, obtained by differencing variables at consecutive time steps (black). Also shown is the sum of the terms that make up the kinetic energy equation (green; eq. 30). Panel (c) shows each of these terms (eqns. [31]-[34]) separately.

The equation governing the horizontal kinetic energy is


By itegrating over the computational volume, we may write




are the contributions from advection, gas pressure gradients, radiation pressure gradients, and the Lorentz force, respectively. We combine the contribution from magnetic pressure and tension, since they cancel each other out to a large degree due to their phase offset during linear growth.

Figure 3b shows the evolution of the rate of change of the horizontal kinetic energy of the same model shown in panel (a). This rate is in excellent agreement with the sum of the terms comprising equation (30), an additional verification of the accuracy of the solution method.

The terms that make up the rate of change of the horizontal kinetic energy are shown separately in Figure 3c. The gas and radiation pressure gradient terms grow at the same rate for e-folding times. The advective and Lorentz force components grow more slowly initially, but then increase their growth rate as the amplitude increases.

Saturation in this model occurs when the Lorentz force reaches the same amplitude as the radiation pressure component. These two contributions are out of phase, as expected from the form of the RMI forcing in the linear regime (Appendix A); the saturation amplitude is small in this model (; Table 2). This result is consistent with that of Turner et al. (2005), who found that the RMI stopped growing when the field ‘buckled’. The much lower radiative forcing in this model means that the amplitude of magnetic field fluctuations required are much smaller. This balance of stresses at saturation does not apply to other models with lower , however. A more careful study of this interplay between magnetic and radiation stress tensors will be addressed in future work.

The growth rates inferred from the evolution of are shown in Figure 4 for the sequence of models with varying radiation pressure and fixed magnetic field (; ; ). When radiation pressure is low, the measured values agree within a factor of a few with the growth rate of the most unstable mode from the dispersion relation of BS03. Agreement improves as radiation makes an increasingly larger contribution to the total pressure, and the degree of forcing becomes larger. Comparing the growth rates measured in the model with for resolutions of and shows a decreasing agreement (at the few percent level) with coarser resolution, pointing to numerical dissipation as the likely cause of the discrepancy with the analytic value at low radiation pressure. The fact that numerical growth rates exceed the analytic prediction at large radiation pressure can perhaps be attributed to long-wavelength effects, as the BS03 value is a correction of the order of to the magnetosonic frequency. We nevertheless take these results as a confirmation of the predictions of BS03 on the growth of the RMI when gas pressure dominates.

Figure 4.— Linear growth rates inferred from the evolution of the horizontal kinetic energy as a function of the ratio of radiation to gas pressure, for the sequence of models with fixed magnetic field (squares). The solid line shows the growth rate of the most unstable mode from the dispersion relation of BS03. The dashed lines are the approximate growth rates of slow (red) and fast (blue) modes from equations (7) and (10), respectively.
Figure 5.— Top: Instantaneous deviation of the density relative to its initial value during saturation, for models with increasing radiation pressure (as labeled) and constant magnetic field (, resolution , ). Curves show a few magnetic field lines. Note that gravity points in the negative direction. Bottom: Profiles of the time- and horizontal root-mean-square fluctuation (eq. [36]) divided by mean value (eq. [35]) of density (blue), radiative flux (red), and magnetic field (green), for the same models shown in the top panels.
Figure 6.— Same as Figure 5, but for a range of models spanning different ratios of magnetic pressure to gas pressure, as labeled. The radiation to gas pressure is fixed to , and the resolution is . When the magnetic field dominates, the usual photon bubble behavior is recovered (Panel d).

Snapshots of the fractional change in density relative to the initial value during saturation are shown in Figure 5 for models with fixed magnetic field and varying radiation to gas pressure ratios. When the radiation pressure is low, amplitudes are small and the flow is dominated by domain-filling structures. The magnetic field becomes significantly deformed when . A larger contribution from radiation pressure increases the saturation amplitude and amount of small-scale structure.

The model with and deserves special attention. The amplitude of the fluctuations becomes large enough that the diffusion solver fails to converge. Future studies should address this region of parameter space using a suitable closure for the radiation moment equations and a larger computational domain.

Another set of snapshots is shown in Figure 6 for models with fixed and increasing magnetic field. Despite the fact that the contribution from radiation is significant, very small amplitudes are obtained when magnetic fields are low. In fact, no growth is seen at all for the lowest field configuration when the resolution is , in contrast to most other models. Decreasing the radiation pressure to and keeping leads to no growth even at our baseline resolution. Increasing the magnetic field strength increases the amplitude and spatial size of the dominant structures up to the point where . Further increases in magnetic field cause the saturation amplitude to decrease slightly, giving rise to the structures seen by Turner et al. (2005).

4.2. Time-Averaged Properties

We now discuss the quantitative behavior of the simulations in the saturated state. Saturation is defined as a period during which the kinetic energy is constant when averaged over timescales longer than the oscillation period or growth time. In some cases (e.g., ), the system achieves an initial phase of saturation lasting a few hundred sound crossing times, but then mass begins to leave the box at an increasingly rapid rate, changing its global properties. We always choose the time interval for averaging such that the mass in the box remains nearly constant.

Figure 7.— Time-averaged properties as a function of the ratio of radiation to gas pressure, for models with (, resolution ). Panel (a) shows the ratio of total kinetic energy to the sum of the gas and radiation energies. Error bars show the size of the corresponding root mean square fluctuation. Panels (b), (c), and (d) show the ratio of the root mean square fluctuation to mean value at domain center of the density, vertical radiative flux, and magnetic field strength, respectively.

The time-average of the volume-integrated vertical kinetic-, horizontal kinetic-, gas-, and radiation energies is denoted by , , , and , respectively. For field variables, we compute time- and horizontal averages


where stands for a scalar such as density. The root-mean-square fluctuation is computed as


For all the models that reach saturation, we report in Table 2 the ratio of total time-average kinetic energy to time-averaged internal plus radiation energies , where , and the ratio of vertical to horizontal kinetic energies . In addition, ratios of root-mean-square fluctuation to mean value of density , vertical component of the radiative flux , and magnetic field strength are provided at three points in the computational box. The lower panels in Figures 5 and 6 show that these averaged variables change considerably with altitude, sometimes by up to an order of magnitude. To capture this variation, values at , , and of the box height are shown in Table 2.

Figure 7 shows the kinetic energy as a function , for the sequence of runs with constant and resolution . The kinetic energy is a monotonic function of the ratio of radiation to gas pressure, increasing the fastest at . A similar trend is observed for the root-mean-square fluctuation of the density, radiative flux, and magnetic field strength, whose values at domain center are also shown in the same figure. Such increase in steady-state amplitudes is clearly the consequence of the fact that the radiation field is the ultimate source of energy (§2.3). The amplitude of the fluctuations in radiative flux and density are comparable to each other, while the magnetic field perturbations are larger by an order of magnitude at box center. This discrepancy is in part a consequence of the different spatial distributions of fluctuations in the box (e.g., lower panels of Figure 5).

The dependence of the saturated quantities on the ratio of magnetic to gas pressure is shown in Figure 8 for the sequence of simulations with . As pointed out in §4.1, the kinetic energy increases with field strength when the inequality obtains. Maximal amplitude is obtained when the contributions to the pressure are all comparable to each other. Further increase in the magnetic field suppresses vertical motions, leading to a lower total kinetic energy. This is clearly shown by the sharp decrease in the ratio of vertical to horizontal kinetic energy when going from to in Figure 8b. The amplitude of the magnetic field fluctuations follow the magnitude of the kinetic energy.

Figure 8.— Time-averaged properties as a function of the ratio of magnetic to gas pressure, for models with (, resolution ). Panel (a) shows the ratio of the total kinetic energy to the sum of the gas and radiation energies. The ratio of the vertical to horizontal kinetic energies is shown in panel (b), and the normalized root mean square fluctuation at domain center of the magnetic field strength in panel (c).

Table 2 also shows that varying the inclination of the magnetic field relative to gravity changes the saturation amplitude. The optimal configuration appears to be a field perpendicular to gravity. When the magnetic field comes close to the vertical direction, the growth rates are eventually suppressed by the sine of the angle between field and radiative flux (equation A.2). This change of saturation amplitude with angle was also seen by Turner et al. (2005).

4.3. Dependence on Numerical Parameters

We have selected a gas-dominated model (, , ) for studying the dependence of the saturation properties on the choice of numerical parameters. Figure 9 shows that the time-averaged kinetic energy in the box as a function of resolution is approaching convergence when using cells per gas pressure scale height in each direction ( cells in the box). The difference with using half the number of cells per direction decreases the kinetic energy by . We thus consider our simulations are spatially well resolved. Note that using less than cells per gas pressure scale height can significantly underestimate the saturation amplitude.

Turner et al. (2005) found that the choice of time step is also relevant for correctly capturing the growth of the instability when using an implicit solver for the radiation sector. We have performed an additional set of simulations (not shown in Table 2) that vary the time-step employed on the model with resolution. No appreciable difference in the saturation amplitude is found when increasing the minimum time step from the fiducial value, which resolves diffusion over a wavelength 3.5), to the diffusion time for a wavelength . Minor differences of the order of dex in the amplitude of fluctuations of field variables are the only difference. One possible cause of this is that the nonlinear phase is dominated by large scale structures that are well resolved.

Figure 9.— Time-averaged kinetic energy as a function of spatial resolution. Parameters of the model are , , and . Error bars denote the root mean square fluctuation.

The size of the box in the horizontal direction can change the value of the saturation amplitude by a factor of order unity. Table 2 shows that the simulation with resolution has less than half the kinetic energy (relative to the sum of the gas and radiation energies) of that where the horizontal size is (). In both cases, the dominant mode fills the domain laterally (as in Figure 5b), and the mode with longer wavelength seems to extract more power from the radiation field. Further increasing the box size to causes the relative kinetic energy to decrease. The dominant mode in this case has a lateral wavelength of . The relative kinetic energy is larger than that of the model with a side of , but less than that of the model with a side of , in agreement with a trend of monotonically increasing power with dominant mode wavelength. Further work should investigate the origin of the dominant mode wavelength, when periodic boundary conditions are used.

5. Application: Massive Stars

Massive main-sequence and blue supergiant stars have radiative envelopes, with an important contribution from radiation pressure to the overall support. Many OB stars have recently been found to have magnetic fields up to kG strength or even higher (see, e.g., Petit et al. 2012 for a recent compilation). Conditions near the surface of these stars are clearly susceptible to the RMI333Following BS03, our analysis assumes a constant magnetic field on length scales comparable to the gas pressure scale height. This is a good starting point to address the prevalence of this instability in stars with radiative envelopes, though it is not necessarily a realistic approximation. .

The winds of massive stars are known to be time-variable and show clumpiness over a wide range of spatial scales. For instance, mass loss rates inferred from UV resonance lines – which scale linearly with density – yield estimates that are at least an order of magnitude lower than those from collisional emission processes (e.g., H), which depend on the density squared (Fullerton et al. 2006).

Variability phenomena are traditionally divided into large- and small-scale. Below we elaborate on the potential impact of the RMI on each of these two categories.

5.1. Large-Scale Variability

It is currently thought that cyclic variability on scales comparable to the stellar radius (e.g., time-dependent absorption superimposed on the blue edge of P Cygni profiles) originates in so-called Corotating Interacting Regions (CIRs), in analogy with the Solar Wind (Mullan 1984). In this model, azimuthal fluctuations at the base of the wind cause faster streams to overtake slower ones, generating shocks and variations in the optical depth that follow a spiral pattern due to angular-momentum conservation (e.g., Cranmer & Owocki 1996, Dessart 2004).

The source of these fluctuations in massive stars is not established yet. Proposed mechanisms include flow in a spatially inhomogeneous magnetic field (Underhill & Fahey 1984) and multi-periodic non-radial pulsations (Kaufer et al. 2006), the latter possibly excited by near surface convection zones from opacity peaks (Cantiello et al. 2009). In at least one source (HD 64760), the spectral line variability can be accounted for if the spots do not corotate with the stellar surface (Lobel & Blomme 2008), lending support to the multi-periodic non-radial pulsation hypothesis. From our results, the RMI is a clear candidate as a source of such pulsations.

5.2. Small-Scale Variability

Line-driven winds are subject to the Line Deshadowing Instability (LDI) above the photosphere (Lucy & Solomon 1970; MacGregor et al. 1979; Owocki & Rybicki 1984). When the velocity of a parcel of gas in the wind is perturbed outward, the corresponding Doppler shift moves the relevant rest-frame transitions into an energy range where the continuum flux from the star is higher, increasing the outward acceleration. The instability tends to favor small spatial scales. The nonlinear development of the LDI results in extended regions of sharp density and velocity fluctuations that lead to the compression of wind material into thin radial shells (Owocki et al. 1988). These shells are subsequently broken up by secondary instabilities, resulting in a clumped wind (Dessart & Owocki 2003). It is currently thought that the LDI can account for most of the phenomenology requiring small scale inhomogeneities, such as transient emission line substructures, black troughs of saturated UV resonance lines, discrepant mass-loss estimates, and X-ray emission (e.g., Fullerton 2011).

A remaining question is the presence of clumping at the base of the wind. Simulations of the ‘self-excited’ LDI find that damping due to scattered radiation (Lucy 1984) suppresses the onset of fluctuations below stellar radii (Runacres & Owocki 2002). A number of observational studies have found evidence for clumping much closer to the base of the wind. In particular, Najarro et al. (2011) used different line diagnostics to reconstruct the radial profile of the clumping factor (average of over the squared average of ) for the entire wind inside stellar radii in Pup, finding that the peak of this factor occurs inside stellar radii. Recently, Sundqvist & Owocki (2012) have reported LDI simulations that agree with the Najarro et al. (2011) clumping profile when including limb-darkening effects (which decrease line-drag from scattered radiation) and generic sound wave perturbations at the base of the wind.

The RMI can readily drive sub-photospheric density perturbations on small spatial scales. BS03 found that in the thermally locked limit (eq. [3]) the instability operates up to a wave number

where is given by equation (6). For a near equipartition magnetic field, this translates into a minimum unstable wavelength of


where is the density in units of  g cm, the acceleration of gravity in units of  cm s, the isothermal sound speed in units of  km s, and is the electron scattering opacity. Comparing to the gas pressure scale height at the photosphere,


shows that the RMI in its WKB version operates over a dynamic range of at least in wavelength.

6. Summary

In this paper we have studied the nonlinear development of the Radiation-Driven Magneto-Acoustic Instability (RMI) in the regime where the gas pressure is comparable or dominates over that due to the radiation and magnetic field, and in the limit of rapid heat exchange between matter and radiation (eq. [3]). Our primary results are:

1. – We have confirmed the predictions of Blaes and Socrates (2003). In particular, the RMI instability operates when gas pressure dominates over radiation pressure and when the magnetic field is weak. The growth rate of the kinetic energy is consistent with that expected for the most unstable mode (Fig. 4).

2. – The saturation amplitude is an increasing function of the ratio of radiation to gas pressure (Fig. 7) over all of the dynamic range studied here, . At the lowest end of this interval, the saturation amplitude is small but non-negligible ().

3. – The saturated state varies non-monotonically with magnetic field strength (Fig. 8). For moderate forcing by radiation (), the saturation amplitude increases with field strength as long as magnetic pressure is sub-dominant. Once the magnetic field dominates, vertical motions are suppressed and increasingly smaller amplitudes are obtained.

4. – The RMI may be a source of sub-photospheric perturbations in the envelopes of massive main-sequence and blue supergiant stars. The instability can serve as a source of large spatial scale perturbations for Corotating Interacting Regions. The RMI also operates over a wide dynamic range in wavelength, providing sufficient seed perturbations for the Line Deshadowing Instability even if a convection zone is absent.

The computational requirements to capture the RMI in realistic contexts are very demanding. The instability is strongest in regions near the photosphere, where diffusion is fast relative to the mode period and where the pressure scale height is much smaller than the stellar radius. Capturing the instability requires resolving diffusion over the pressure scale height in space and in time, introducing an additional computational cost. For example, the ratio of the sound crossing time over the diffusion time at a given wavelength is proportional to the diffusion speed over the sound speed (, eq. [13]). For massive stars this ratio can be as high as , requiring a comparable number of diffusion steps per Courant time to capture the RMI.

Future work will address the saturation of the RMI, and further application to astrophysical systems.

We thank the anonymous referee for helpful comments that improved the manuscript. RF is supported in part by NASA through Einstein Postdoctoral Fellowship grant number PF-00062, awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060, and by NSF grant number AST-0807444. AS acknowledges support from a John N. Bahcall Fellowship, awarded by the Institute for Advanced Study, Princeton. This research was also supported in part by the National Science Foundation through TeraGrid/XSEDE resources (Catlett 2007). Computations were performed at LONI Queen Bee and the IAS Aurora cluster.

Appendix A Derivation of Linear Growth Rates and Stability Criteria

a.1. Basic Magneto-Acoustic Waves

In the absence of stratification and a changing radiation field, compressive fluid perturbations are forced by a combination of their own acoustic response and the Lorentz force. In this limit, conservation of momentum to linear order can be written as


where denotes Eulerian perturbation, and when (eq. [3]). We assume that all perturbed quantities take plane-wave form and are . Using the linearized mass conservation equation, the linearized induction equation


and the component of the linearized momentum equation along , one can rewrite equation (A1) as


where . Projecting onto and using mass conservation yields the magnetosonic dispersion relation,


a.2. Radiative Correction: Driving and Diffusive Damping

We now include stratification in the background state and allow the radiation field (or equivalently the temperature) to vary. The exchange of heat between radiation and energy is treated as a correction to the mode frequency. We thus have


where is a driving or damping rate such that . The linearized momentum equation is now


where is the Lagrangian displacement, and is the radiation pressure in the optically thick limit.

Defining , using the condition of hydrostatic equilibrium , and expressing Eulerian perturbations in terms of the Lagrangian displacement


one can re-express momentum balance as




contains the driving and damping terms. The second term on the left hand side of equation (A11) includes the terms arising from acoustic and magnetic forces (equation A3). The second term on the right hand side is perpendicular to the Lagrangian displacement and hence cannot do work on the perturbation (e.g., Socrates et al. 2005). Terms of higher order in have been discarded.

To first order in , the growth rate is


Before evaluating equation (A13), it is instructive to explicitly separate the driving and damping terms in equation (A12). To do so we solve for the temperature perturbation in the linearized equation for the total energy density. Defining , we have






To first order in and in , the temperature perturbation is


The first term represents damping by radiative diffusion (eq. [A16]): entropy is lost upon compression444In cosmology, this effect is known as “Silk damping”, after Silk (1968).. Note that we have dropped the first term on the left hand side of equation (A.2); it becomes important if , contributing with a global factor of order unity to the growth rate (e.g., compare and the maximum growth rate from the dispersion relation in Figure 4; the choice of parameters is such that ).

We now evaluate equation (A13) by inserting equation (A17) into equation (A12), projecting onto , using equation (A3) to eliminate a term proportional to , and keeping the same order of approximation, we find

The second term inside the square brackets represents driving by the background radiative flux (BS03). Note that only the component of the flux perpendicular to contributes to the driving, and that driving vanishes when the magnetic field is completely parallel or perpendicular to the wave vector. By using equation (A3), the denominator of equation (A13) can be written as


Equations (2), (A13), (A.2), (A16) and (A19) then yield the growth rate to first order in , , and ,