format: text

text:
# Local Radiative Hydrodynamic and Magnetohydrodynamic Instabilities in Optically Thick Media

## 1 Introduction

It is well known that the presence of a substantial equilibrium radiation pressure gradient can destabilize optically thick astrophysical media. Prendergast & Spiegel (1973) first speculated that compressible fluid flow subject to a large radiative flux would be unstable to the formation of buoyant rarefied regions of enhanced radiation pressure, better known as “photon bubbles.” The existence of these instabilities has never been fully demonstrated in hydrodynamic systems, but radiation pressure acting in stellar envelopes with opacity peaks can drive certain oscillation modes (now known as “strange modes”) unstable (e.g. Wood 1976; Saio, Wheeler, & Cox 1984; Gautschy 1993; Kiriakidis, Fricke, & Glatzel 1993). Strange modes have also been found to exist in hydrodynamic models of accretion disks (Glatzel & Mehren 1996, Mehren-Baehr & Glatzel 1999).

Gautschy & Glatzel (1990) showed that strange mode instabilities satisfy a “non-adiabatic reversible” (NAR) approximation, in which the modes occur in a medium with effectively zero specific heat capacity, resulting in a vanishing luminosity perturbation. A physical model for strange mode instabilities within the NAR approximation was established by Glatzel (1994). He noted that opacity peaks act as a reflection barrier for high overtone acoustic modes, sonically decoupling different radial regions in the star. The basic mechanism for the instability was the relative phase shift between the pressure and density perturbations, resulting from the constraint of having zero luminosity perturbation. Further aspects of the physics of strange modes have been elucidated by Papaloizou et al. (1997) and Saio, Baker, & Gautschy (1998). In numerical studies, Shaviv (2001) has found that both standing and propagating acoustic waves are unstable in atmospheres with sufficient radiation pressure support, even when the opacity is pure Thomson scattering.

Radiation pressure driven instabilities occur in magnetized systems as well. Arons (1992) identified and investigated photon bubble instabilities in strongly magnetized atmospheres of accreting X-ray pulsars. He showed that such instabilities are caused by the enhanced buoyancy that occurs when radiation diffuses into relatively rarefied regions in compressible perturbations. With accretion disk applications in mind, Gammie (1998) performed a linear, local stability analysis of a static, magnetized, stratified medium and found an instability over a finite range of wavenumbers that he identified as being similar in nature to photon bubbles. Using the simplifying assumption that gas and radiation were coupled together just by Thomson scattering, Blaes & Socrates (2001, hereafter BS01) found that Gammie’s (1998) instability could be extended to arbitrarily short wavelengths where it manifested itself as an overstable slow mode. In addition, the fast magnetosonic mode could also be unstable for sufficiently strong radiation fluxes, albeit with smaller growth rate.111BS01 also found that Alfvén waves could be destabilized, an effect that is almost certainly due to the fact that they considered a rotating equilibrium and these waves then have a small compressible component because rotation couples them to slow modes. Begelman (2001) has constructed a periodic shock train solution that may represent one possible nonlinear outcome of the slow mode instability.

The purpose of this paper is to investigate the local, WKB version of these compressive instabilities in a unified manner. We focus on the local driving of propagating waves, and do not address here the issue of whether and how this local driving can manifest itself as a global instability of a standing, normal mode of the object in question. We treat the thermodynamics of the medium quite generally, making no restrictive assumptions on the equality of the gas and radiation temperature in the perturbed state or using the NAR approximation. We find that hydrodynamic instabilities driven by radiative diffusion exist under a wide range of conditions, for both gas and radiation pressure dominated media. However, a medium supported purely by Thomson scattering opacity exhibits no local driving of acoustic wave instabilities, in contrast to the conclusions of Shaviv (2001). On the other hand, the anisotropic stress produced by magnetic tension widens the applicability of these instabilities to even broader classes of equilibria, even those which have opacities given by Thomson scattering only. The MHD instabilities can exist even when the equilibrium magnetic pressure is less than either the gas or radiation pressure.

All these instabilities may grow on the dynamical time scale or even faster. The instability mechanism originates from the interaction of the equilibrium radiative flux with density and opacity fluctuations in the perturbed flow, as well as changes in radiation pressure along fluid displacements. The ultimate source of free energy is provided by the stratified radiation field.

This paper is organized as follows. In section 2, we describe our basic equations and assumptions, concentrating in particular on the thermodynamics of the coupled gas and radiation fluids, and also derive the general dispersion relation for waves in a stratified, magnetized medium within the local (WKB) limit. In section 3 we discuss the solutions of this dispersion relation in the hydrodynamic limit where there is no magnetic field, generalizing and extending the work of previous authors. We discuss in some detail the physics underlying these instabilities in section 4. We then incorporate the effects of magnetic fields in section 5, and in section 6 discuss how the physics of magnetoacoustic wave instabilities is related to the hydrodynamic instabilities of the previous sections. In section 7 we briefly discuss astrophysical applications of these instabilities to accretion disks and stars, deferring more detailed applications to observed phenomena for future work. Section 8 summarizes our main conclusions. Our WKB analysis in the main body of the paper is not completely rigorous, and we provide more solid mathematical justification for our conclusions for two particular cases in an appendix. We also include an additional appendix where we examine the effects of radiative diffusion on the magnetorotational (MRI) instability in differentially rotating flows, generalizing our previous work (BS01) to include the more generic thermodynamics we use here. Readers not interested in the details of the linear perturbation theory may wish to skip ahead to sections 7 and 8 directly. Readers interested in obtaining a basic physical understanding of the causes of these instabilities may wish to focus on sections 4 and 6.

## 2 Equations and Assumptions

The basic equations of radiation magnetohydrodynamics (RMHD) have been discussed by Stone, Mihalas, & Norman (1992), and we adopt these equations here with slight changes in notation as well as some changes in the physics. The fluid equations we consider are

$$\frac{\rho }{t}+\text{r}(\rho 𝐯)=0,$$
(1)
$$\rho \left(\frac{𝐯}{t}+𝐯\text{r}𝐯\right)=\text{r}p+\rho 𝐠+\frac{1}{4\pi }(\text{r}\times 𝐁)\times 𝐁+\frac{\kappa _\mathrm{F}\rho }{c}𝐅,$$
(2)
$$\frac{u}{t}+𝐯\text{r}u+\gamma u\text{r}𝐯=\kappa _\mathrm{J}\rho cE\kappa _\mathrm{P}\rho caT_\mathrm{g}^4\kappa _\mathrm{T}\rho c\left(\frac{4k_\mathrm{B}T_\mathrm{g}}{m_\mathrm{e}c^2}\frac{h\overline{\nu }}{m_\mathrm{e}c^2}\right)E,$$
(3)
$$\frac{E}{t}+𝐯\text{r}E+\frac{4}{3}E\text{r}𝐯=\text{r}𝐅\kappa _\mathrm{J}\rho cE+\kappa _\mathrm{P}\rho caT_\mathrm{g}^4+\kappa _\mathrm{T}\rho c\left(\frac{4k_\mathrm{B}T_\mathrm{g}}{m_\mathrm{e}c^2}\frac{h\overline{\nu }}{m_\mathrm{e}c^2}\right)E,$$
(4)
$$0=\frac{1}{3}\text{r}E\frac{\kappa _\mathrm{F}\rho }{c}𝐅,$$
(5)
and

$$\frac{𝐁}{t}=\text{r}\times (𝐯\times 𝐁).$$
(6)
Here $`\rho `$, $`p`$, $`T_\mathrm{g}`$, and $`u`$ are the density, pressure, temperature, and energy density in the gas, respectively. These quantities are related to each other by

$$u=\frac{p}{\gamma 1}$$
(7)
and

$$p=\frac{\rho k_\mathrm{B}T_\mathrm{g}}{\mu },$$
(8)
where $`\gamma `$ is the ratio of specific heats in the gas and $`\mu `$ is the mean molecular mass of the gas. We assume throughout this paper that $`\gamma `$ and $`\mu `$ are constant, and thereby neglect the effects of composition gradients, ionization, and recombination.

Other symbols in equations (1)-(8) have their usual meaning: $`c`$ is the speed of light, $`a`$ is the radiation density constant, $`k_\mathrm{B}`$ is Boltzmann’s constant, and $`m_\mathrm{e}`$ is the electron mass. The vector $`𝐯`$ is the fluid velocity, $`𝐁`$ is the magnetic field, and $`𝐠`$ is the gravitational acceleration. We assume that $`𝐠`$ is time-independent throughout this paper, either because the gravitational field is due to some fixed external source (as is the case for a non-self-gravitating accretion disk), or because we choose to consider only short-wavelength perturbations that are well-described by the Cowling approximation.

The radiation energy density $`E`$ and the radiation flux $`𝐅`$ are defined in terms of frequency integrals over the radiation spectrum as measured in the local fluid rest frame:

$$E\frac{4\pi }{c}_0^{\mathrm{}}𝑑\nu J_\nu \frac{4\pi }{c}J,$$
(9)
where $`J_\nu `$ is the angle-averaged mean specific intensity in the local fluid rest frame, and

$$𝐅_0^{\mathrm{}}𝑑\nu 𝐅_\nu .$$
(10)
We have simplified the RMHD equations by assuming that the radiation stress tensor is isotropic in the local fluid rest frame, i.e. that the tensor variable Eddington factor defined by Stone et al. (1992) is one third times the identity matrix. We have therefore neglected the effects of photon viscosity. We have also neglected terms corresponding to radiation inertia in the radiation momentum equation (5), leaving us with just the radiation diffusion equation. Such terms are generally negligible for the low frequency instabilities we explore in this paper (e.g. BS01).

The two opacities $`\kappa _\mathrm{J}`$ and $`\kappa _\mathrm{F}`$ are also defined by integrals over the radiation spectrum in the local rest frame:

$$\kappa _\mathrm{J}\frac{1}{\rho J}_0^{\mathrm{}}𝑑\nu \chi _\nu ^{\mathrm{th}}(\rho ,T_\mathrm{g})J_\nu ,$$
(11)
and

$$\kappa _\mathrm{F}𝐅\frac{1}{\rho }_0^{\mathrm{}}𝑑\nu [\chi _\nu ^{\mathrm{th}}(\rho ,T_\mathrm{g})+n_\mathrm{e}\sigma _\mathrm{T}]𝐅_\nu .$$
(12)
Here $`\chi _\nu ^{\mathrm{th}}(\rho ,T_\mathrm{g})`$ is the thermal absorption coefficient at frequency $`\nu `$, $`n_\mathrm{e}`$ is the electron number density, and $`\sigma _\mathrm{T}`$ is the Thomson cross-section. The Planck mean opacity is

$$\kappa _\mathrm{P}\frac{4\pi }{\rho acT_\mathrm{g}^4}_0^{\mathrm{}}𝑑\nu \chi _\nu ^{\mathrm{th}}(\rho ,T_\mathrm{g})B_\nu (T_\mathrm{g}),$$
(13)
where $`B_\nu (T_\mathrm{g})`$ is the Planck function at the gas temperature. The Thomson opacity is

$$\kappa _\mathrm{T}\frac{n_e\sigma _\mathrm{T}}{\rho }$$
(14)
and is constant due to our neglect of ionization/recombination processes and composition gradients. Note that we have extended the RMHD equations of Stone et al. (1992) to include electron scattering. Coherent (Thomson) scattering merely allows the gas and radiation to exchange momentum, not energy. However, we have also allowed for incoherent (Compton) scattering in the last terms on the right hand side of the gas and radiation energy equations (3) and (4). We do this by treating Compton scattering in the frequency diffusion (Kompaneets) limit, neglecting anisotropies in the radiation field in the local fluid rest frame. A derivation of the Compton scattering terms can be found, e.g., in section 3.1 of Hubeny et al. (2001). The quantity $`\overline{\nu }`$ is an average frequency of the radiation field in the local rest frame of the fluid, defined by

$$\overline{\nu }\frac{1}{J}_0^{\mathrm{}}𝑑\nu \nu J_\nu \left(1+\frac{J_\nu c^2}{2h\nu ^3}\right).$$
(15)
The second term in parentheses represents the effects of stimulated scattering, and guarantees that there is no net heat exchange between the gas and the radiation when thermal equilibrium is established with the radiation field being blackbody at the local gas temperature.

In principle, evaluation of the frequency integrals in the opacities and $`\overline{\nu }`$ requires a full solution of the radiative transfer equation coupled to the fluid flow. Rather than do that here, we adopt some simplified prescriptions that we believe will not alter the essential physics we wish to explore. First, we adopt the Eddington closure relation between the zeroth and second moments of the comoving frame specific intensity, consistent with our earlier assumption that the tensor variable Eddington factor is one third times the identity matrix. Second, because we will only consider perturbations about equilibria that are in local thermodynamic equilibrium (LTE), we assume that non-LTE departures in the perturbed flow can be adequately parameterized by a blackbody radiation spectrum at a temperature $`T_\mathrm{r}`$, which may be different from the gas temperature $`T_\mathrm{g}`$.222Departures from LTE will of course generally not maintain a blackbody radiation field in the perturbed flow. Strictly speaking, we are defining $`T_\mathrm{r}`$ as the effective temperature of the radiation field, and we are assuming that the various frequency moments we need are related to each other in the perturbed flow in approximately the same way as they would be if the radiation field were blackbody.

Under these assumptions, equation (15) implies that

$$\overline{\nu }=\frac{4k_\mathrm{B}T_\mathrm{r}}{h}.$$
(16)
Equation (11) gives an expression for $`\kappa _\mathrm{J}`$ which more closely resembles the Planck mean opacity $`\kappa _\mathrm{P}`$,

$$\kappa _\mathrm{J}=\frac{4\pi }{\rho acT_\mathrm{r}^4}_0^{\mathrm{}}𝑑\nu \chi _\nu ^{\mathrm{th}}(\rho ,T_\mathrm{g})B_\nu (T_\mathrm{r}).$$
(17)
In addition, equation (12) can be replaced with an expression analogous to the Rosseland mean opacity,

$$\kappa _\mathrm{F}=\left(_0^{\mathrm{}}𝑑\nu \frac{dB_\nu }{dT_\mathrm{r}}\right)\left\{_0^{\mathrm{}}𝑑\nu \left[\frac{\chi _\nu ^{\mathrm{th}}(\rho ,T_\mathrm{g})}{\rho }+\kappa _\mathrm{T}\right]^1\frac{dB_\nu }{dT_\mathrm{r}}\right\}^1.$$
(18)
Unlike $`\kappa _\mathrm{P}`$, which is a function only of the gas density and temperature, $`\kappa _\mathrm{J}`$ and $`\kappa _\mathrm{F}`$ are also functions of the radiation temperature. Finally, the blackbody assumption implies that $`T_\mathrm{r}`$ can be related to the radiation energy density through the equation of state

$$E=aT_\mathrm{r}^4.$$
(19)
To summarize, equations (1)-(8), (13), and (16)-(19) are the basic dynamical equations we will use throughout the rest of this paper.

### 2.1 Equilibrium

For the majority of this paper we will consider the dynamical stability of static equilibria in which the equilibrium fluid velocity is zero. As mentioned in the previous section, we also assume that the equilibrium is in LTE, i.e. that the unperturbed gas and radiation temperatures are identical,

$$T_\mathrm{g}=T_\mathrm{r}T.$$
(20)
This can be somewhat problematic for applications to accretion disk models, which in some cases can be effectively thin. However, we do not presume to model the turbulent dissipation that must exist in these flows. In fact, our purpose in that application is to consider equilibria that might be useful starting points for exploring the development of this turbulence by simulation. Lacking an explicit treatment of dissipation, we are forced to consider equilibria that are in LTE and in local radiative equilibrium.

Under these assumptions, the only nontrivial partial differential equations describing our equilibrium are those expressing hydrostatic equilibrium,

$$0=\rho 𝐠\text{r}p+\frac{1}{4\pi }(\text{r}\times 𝐁)\times 𝐁+\frac{\kappa _\mathrm{F}\rho }{c}𝐅,$$
(21)
radiative equilibrium,

$$0=\text{r}𝐅,$$
(22)
and radiative diffusion,

$$0=\frac{1}{3}\text{r}E\frac{\kappa _\mathrm{F}\rho }{c}𝐅.$$
(23)
Throughout this paper we will assume that the equilibrium magnetic field is uniform, so that the Lorentz force term in equation (21) vanishes. Equation (21) then has the immediate consequence that all thermodynamic variables are constant on horizontal surfaces, i.e. those surfaces that are perpendicular to the gravitational acceleration $`𝐠`$. We define $`z`$ to be a vertical coordinate that increases upward, and $`\widehat{𝐳}`$ to be the corresponding upward unit vector. Then $`𝐠g\widehat{𝐳}`$, with $`g>0`$.

### 2.2 The Linear Perturbation in Total Pressure

When we consider the behavior of linear perturbations about the equilibrium just discussed, the problem naturally divides into two parts: the thermodynamics of the perturbed flow, and the coupling of this with the perturbed magnetic and velocity fields. We examine the former problem first by determining the perturbed total (gas plus radiation) pressure.

We begin by considering expressions for the perturbations in $`\kappa _\mathrm{J}`$ and $`\kappa _\mathrm{P}`$. The functional dependence of these two opacities on $`\rho `$, $`T_\mathrm{g}`$ and $`T_\mathrm{r}`$ immediately implies that

$$\frac{\delta \kappa _\mathrm{J}}{\kappa _\mathrm{J}}=\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}\rho }\frac{\delta \rho }{\rho }+\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}T_\mathrm{g}}\frac{\delta T_\mathrm{g}}{T_\mathrm{g}}+\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}T_\mathrm{r}}\frac{\delta T_\mathrm{r}}{T_\mathrm{r}}$$
(24)
and

$$\frac{\delta \kappa _\mathrm{P}}{\kappa _\mathrm{P}}=\frac{\mathrm{ln}\kappa _\mathrm{P}}{\mathrm{ln}\rho }\frac{\delta \rho }{\rho }+\frac{\mathrm{ln}\kappa _\mathrm{P}}{\mathrm{ln}T_\mathrm{g}}\frac{\delta T_\mathrm{g}}{T_\mathrm{g}}.$$
(25)
(We use Eulerian perturbations throughout this paper, i.e. $`\delta Q`$ is the Eulerian perturbation in the quantity $`Q`$.) Here all opacity derivatives are evaluated in the equilibrium state which is in LTE with $`T_\mathrm{r}=T_\mathrm{g}`$ and $`\kappa _\mathrm{P}=\kappa _\mathrm{J}`$. Because of this, equations (13) and (17) immediately imply

$$\frac{\mathrm{ln}\kappa _\mathrm{P}}{\mathrm{ln}\rho }=\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}\rho }$$
(26)
and

$$\frac{\mathrm{ln}\kappa _\mathrm{P}}{\mathrm{ln}T_\mathrm{g}}=\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}T_\mathrm{g}}+\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}T_\mathrm{r}}.$$
(27)
These two conditions on the absorption opacity derivatives greatly simplify the mathematical form of the thermal coupling between the gas and radiation in the perturbed flow. Linearizing the gas energy equation (3) and the radiation energy equation (4), and using equations (16), (19), and (26)-(27), we obtain

$$\frac{\delta u}{t}+\delta 𝐯\text{r}u+\gamma u\text{r}\delta 𝐯=\omega _\mathrm{a}E\frac{4(\delta T_\mathrm{r}\delta T_\mathrm{g})}{T}$$
(28)
and

$$\frac{\delta E}{t}+\delta 𝐯\text{r}E+\frac{4}{3}E\text{r}\delta 𝐯=\text{r}\delta 𝐅\omega _\mathrm{a}E\frac{4(\delta T_\mathrm{r}\delta T_\mathrm{g})}{T}.$$
(29)
Here $`\omega _\mathrm{a}`$ is an angular frequency describing the rate at which radiation and matter are thermally coupled in the perturbed flow, defined by

$$\omega _\mathrm{a}\left[\kappa _\mathrm{P}\left(1+\frac{1}{4}\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}T_\mathrm{r}}\right)+\kappa _\mathrm{T}\frac{k_\mathrm{B}T}{m_\mathrm{e}c^2}\right]\rho c\kappa _\mathrm{a}\rho c.$$
(30)
Equations (28) and (29) can be further manipulated in a physically revealing manner. We eliminate $`\text{r}\delta 𝐯`$ from them using the perturbed continuity equation (1),

$$\frac{\delta \rho }{t}+\rho \text{r}\delta 𝐯+\delta 𝐯\text{r}\rho =0.$$
(31)
In addition, we eliminate $`\delta 𝐅`$ using the perturbed radiative diffusion equation (5), which can be written as

$$\delta 𝐅=\frac{c}{3\kappa _\mathrm{F}\rho }\text{r}\delta E𝐅\left(\frac{\delta \rho }{\rho }+\frac{\delta \kappa _\mathrm{F}}{\kappa _\mathrm{F}}\right).$$
(32)
Eliminating $`\delta u`$ using equation (7), equation (28) becomes

$$\frac{\delta p}{t}c_\mathrm{g}^2\frac{\delta \rho }{t}=4(\gamma 1)\omega _\mathrm{a}E\frac{(\delta T_\mathrm{r}\delta T_\mathrm{g})}{T}(\gamma 1)\rho T\delta 𝐯\text{r}S_\mathrm{g},$$
(33)
where $`c_\mathrm{g}`$ is the adiabatic sound speed in the gas,

$$c_\mathrm{g}\left(\frac{\gamma p}{\rho }\right)^{1/2},$$
(34)
and $`S_\mathrm{g}`$ is the entropy per unit mass in the gas,

$$S_\mathrm{g}\frac{k_\mathrm{B}}{\mu (\gamma 1)}\mathrm{ln}(p\rho ^\gamma )+\mathrm{constant}.$$
(35)
Equation (29) can be written in a very similar fashion,

$$\frac{}{t}\left(\frac{1}{3}\delta E\right)c_\mathrm{r}^2\frac{\delta \rho }{t}=\frac{4}{3}\omega _\mathrm{a}E\frac{(\delta T_\mathrm{r}\delta T_\mathrm{g})}{T}\frac{1}{3}\rho T\delta 𝐯\text{r}S_\mathrm{r}+\text{r}\left(\frac{c}{3\kappa _\mathrm{F}\rho }\text{r}\frac{1}{3}\delta E\right)+\frac{1}{3}𝐅\text{r}\left(\frac{\delta \rho }{\rho }+\frac{\delta \kappa _\mathrm{F}}{\kappa _\mathrm{F}}\right),$$
(36)
where $`c_\mathrm{r}`$ is the adiabatic sound speed in the radiation,

$$c_\mathrm{r}\left(\frac{4E}{9\rho }\right)^{1/2},$$
(37)
and $`S_\mathrm{r}`$ is the entropy per unit mass in the radiation,

$$S_\mathrm{r}=\frac{4E}{3\rho T}+\mathrm{constant}.$$
(38)
Note that we have simplified the terms involving the equilibrium flux $`𝐅`$ in equation (36) by using the radiative equilibrium equation (22).

Thus far, we have not made any assumptions about the wavelengths of the perturbations we are considering, and we could proceed from this point by doing a full global perturbation analysis. We wish to focus on local instabilities in this paper, however. Hence from now on, we invoke the WKB ansatz by adopting a plane wave spacetime dependence $`\mathrm{exp}[i(𝐤𝐫\omega t)]`$ for all perturbations. Here $`𝐫`$ is the position vector of the point of interest, $`𝐤`$ is the perturbation wave vector, and $`\omega `$ is the perturbation angular frequency. As a result, the time derivative operator $`/t`$ in equations (33) and (36) is replaced by $`i\omega `$, and the spatial derivative operator r in the last two terms of equation (36) is replaced by $`i𝐤`$, i.e.

$$\text{r}\left(\frac{c}{3\kappa _\mathrm{F}\rho }\text{r}\frac{1}{3}\delta E\right)+\frac{1}{3}𝐅\text{r}\left(\frac{\delta \rho }{\rho }+\frac{\delta \kappa _\mathrm{F}}{\kappa _\mathrm{F}}\right)\frac{ck^2}{3\kappa _\mathrm{F}\rho }\left(\frac{\delta E}{3}\right)+\frac{1}{3}i𝐤𝐅\left(\frac{\delta \rho }{\rho }+\frac{\delta \kappa _\mathrm{F}}{\kappa _\mathrm{F}}\right).$$
(39)
The diffusive instabilities in which we are interested are driven by background gradients, in particular the presence of a nonzero equilibrium radiative flux $`𝐅=(c/3\kappa _\mathrm{F}\rho )\text{r}E`$. Hence this last step may at first sight appear to be dangerous: gradients in $`c/(3\kappa _\mathrm{F}\rho )`$ arising from the first term in expression (39) might be as important as the second term that depends on $`𝐅`$ as we take the short wavelength limit $`k\mathrm{}`$. As may be verified a posteriori, however, it turns out that the instabilities have $`\delta E`$ scaling as $`k^1\delta p`$ at high $`k`$, so that our ordering is in fact consistent. This makes physical sense, because at short wavelengths radiative diffusion is very fast and perturbations in the radiation temperature are smoothed out.

After some algebra, equations (8), (19), (33), and (36) can be used to derive an equation for the total (gas plus radiation) pressure perturbation in terms of the density and velocity perturbations,

$$\delta p+\frac{1}{3}\delta E=(𝒜c_\mathrm{g}^2+c_\mathrm{r}^2+𝒞)\delta \rho \frac{i\rho T}{\omega }\delta 𝐯\left[(\gamma 1)𝒜\text{r}S_\mathrm{g}+\frac{1}{3}\text{r}S_\mathrm{r}\right].$$
(40)
Here

$$𝒜\frac{\omega }{𝒟}\left[\omega +\frac{ick^2}{3\kappa _\mathrm{F}\rho }+i\omega _\mathrm{a}\left(1+\frac{4E}{3p}\right)+𝐤𝐅\left(\frac{\mathrm{\Theta }_{T\mathrm{r}}}{4E}\frac{\mathrm{\Theta }_{T\mathrm{g}}}{3p}\right)\right],$$
(41)
$$\frac{\omega }{𝒟}\left[\omega +i\omega _\mathrm{a}(\gamma 1)\left(3+\frac{4E}{p}\right)\right],$$
(42)
$`𝒞`$ $``$ $`{\displaystyle \frac{1}{𝒟}}\{{\displaystyle \frac{4i\omega _\mathrm{a}E}{\rho }}[(\gamma {\displaystyle \frac{4}{3}})\omega +{\displaystyle \frac{ick^2}{3\kappa _\mathrm{F}\rho }}(\gamma 1)]`$ (43)
$`{\displaystyle \frac{𝐤𝐅}{\rho }}[{\displaystyle \frac{\omega }{3}}(1+\mathrm{\Theta }_\rho \mathrm{\Theta }_{T\mathrm{g}})+i(\gamma 1)\omega _\mathrm{a}(1+\mathrm{\Theta }_\rho \mathrm{\Theta }_{T\mathrm{g}}\mathrm{\Theta }_{T\mathrm{r}})`$
$`+i(\gamma 1)\omega _\mathrm{a}\left({\displaystyle \frac{4E}{3p}}\right)(1+\mathrm{\Theta }_\rho )]\},`$
and

$$𝒟\left(\omega +\frac{ick^2}{3\kappa _\mathrm{F}\rho }+\frac{𝐤𝐅}{4E}\mathrm{\Theta }_{T\mathrm{r}}\right)\left[\omega +\frac{4i(\gamma 1)\omega _\mathrm{a}E}{p}\right]+i\omega _\mathrm{a}\left[\omega +(\gamma 1)\frac{𝐤𝐅}{p}\mathrm{\Theta }_{T\mathrm{g}}\right].$$
(44)
The $`\mathrm{\Theta }`$ quantities are defined as logarithmic derivatives of the flux mean opacity with respect to the variable in the subscript, i.e.

$$\mathrm{\Theta }_\rho \frac{\mathrm{ln}\kappa _\mathrm{F}}{\mathrm{ln}\rho },\mathrm{\Theta }_{T\mathrm{g}}\frac{\mathrm{ln}\kappa _\mathrm{F}}{\mathrm{ln}T_\mathrm{g}},\mathrm{and}\mathrm{\Theta }_{T\mathrm{r}}\frac{\mathrm{ln}\kappa _\mathrm{F}}{\mathrm{ln}T_\mathrm{r}}.$$
(45)
Equation (40) expresses all the thermal physics and is the one we shall use throughout the rest of the paper.

### 2.3 Coupling to Density, Velocity, and Magnetic Field Perturbations

Now that we have determined the total pressure perturbation, the only work remaining is to couple this to the perturbed continuity, gas momentum, and flux-freezing equations. Employing the WKB ansatz, the perturbed continuity equation (31) becomes

$$i\omega \delta \rho +i\rho 𝐤\delta 𝐯+\delta 𝐯\text{r}\rho =0.$$
(46)
The perturbed gas momentum equation and flux freezing equations may be written as

$$i\omega \rho \delta 𝐯=i𝐤\left(\delta p+\frac{1}{3}\delta E\right)+𝐠\delta \rho +\frac{i}{4\pi }(𝐤\times \delta 𝐁)\times 𝐁$$
(47)
and

$$i\omega \delta 𝐁=i𝐤\times (\delta 𝐯\times 𝐁),$$
(48)
respectively.

### 2.4 The Dispersion Relation

After some algebra, equations (40) and (46)-(48) may be combined to give a dispersion relation for short wavelength modes on a static, stratified and magnetized equilibrium:

$`0`$ $`=`$ $`\stackrel{~}{\omega }^2\left[\omega ^4\omega ^2k^2v_\mathrm{A}^2\stackrel{~}{\omega }^2k^2(𝒜c_\mathrm{g}^2+c_\mathrm{r}^2+𝒞)\right]`$ (49)
$`+`$ $`\stackrel{~}{\omega }^2(k^2k_z^2)\left(𝒜c_\mathrm{g}^2N_\mathrm{g}^2+c_\mathrm{r}^2N_\mathrm{r}^2\right)`$
$`+`$ $`i\stackrel{~}{\omega }^2\left[\omega ^2𝐤k^2(𝐤𝐯_\mathrm{A})𝐯_\mathrm{A}\right]\left[{\displaystyle \frac{1}{\rho }}\text{r}p\left(𝒜1\right)+{\displaystyle \frac{1}{3\rho }}\text{r}E\left(1\right)+𝒞\text{r}\mathrm{ln}\rho \right]`$
$``$ $`\left[\omega ^4\omega ^2k^2v_\mathrm{A}^22\omega ^2k_zv_{\mathrm{A}z}𝐤𝐯_\mathrm{A}+\omega ^2k_z^2v_\mathrm{A}^2+k^2(𝐤𝐯_\mathrm{A})^2v_{\mathrm{A}z}^2\right]𝐠\text{r}\mathrm{ln}\rho ,`$
where

$$\stackrel{~}{\omega }^2\omega ^2(𝐤𝐯_\mathrm{A})^2,$$
(50)
$`𝐯_\mathrm{A}𝐁/(4\pi \rho )^{1/2}`$ is the vector Alfvén speed,

$$N_\mathrm{g}^2𝐠\left(\frac{1}{\rho c_\mathrm{g}^2}\text{r}p\text{r}\mathrm{ln}\rho \right)=\frac{(\gamma 1)\rho T}{\gamma p}𝐠\text{r}S_\mathrm{g}$$
(51)
is the Brunt-Väisälä frequency in the gas,

$$N_\mathrm{r}^2𝐠\left(\frac{1}{3\rho c_\mathrm{r}^2}\text{r}E\text{r}\mathrm{ln}\rho \right)=\frac{3\rho T}{4E}𝐠\text{r}S_\mathrm{r}$$
(52)
is the Brunt-Väisälä frequency in the radiation, and the quantities $`𝒜`$, $``$, and $`𝒞`$ are defined by equations (41)-(44).

Equation (49) is a dispersion relation for eight modes. We originally had nine first order, time-dependent perturbuation equations (one continuity, three momentum, two energy, and three flux-freezing), and so we expect nine modes. However, one of these modes has zero frequency and is inconsistent with the additional constraint of Gauss’ Law that $`\text{r}\delta 𝐁=0`$. In sections 3 and 5, we provide a full discussion of the unstable waves contained in equation (49).

## 3 Hydrodynamic Instabilities

Before considering the full effects of MHD, it is useful to first explore radiative diffusion instabilities in the hydrodynamic case. Setting $`𝐯_\mathrm{A}=0`$ in the general dispersion relation (49), we obtain the following equation for hydrodynamic modes:

$`0`$ $`=`$ $`\omega ^4\omega ^2k^2(𝒜c_\mathrm{g}^2+c_\mathrm{r}^2+𝒞)+(k^2k_z^2)\left(𝒜c_\mathrm{g}^2N_\mathrm{g}^2+c_\mathrm{r}^2N_\mathrm{r}^2\right)`$ (53)
$`+`$ $`i\omega ^2𝐤\left[{\displaystyle \frac{1}{\rho }}\text{r}p\left(𝒜1\right)+{\displaystyle \frac{1}{3\rho }}\text{r}E\left(1\right)+𝒞\text{r}\mathrm{ln}\rho \right]\omega ^2𝐠\text{r}\mathrm{ln}\rho .`$
This equation can be written as a sixth order polynomial in $`\omega `$, which is to be expected given that it arises from six time-dependent hydrodynamic perturbation equations. However, one power of $`\omega `$ can be factored out, indicating that one of the modes always has zero frequency. This mode has zero total pressure perturbation $`\delta p+\delta E/3`$, at least when the wave vector is not entirely vertical.

### 3.1 Short Wavelength Limit

In the short wavelength limit, the remaining five modes described by equation (53) can be easily factored:

$$0=\left(\omega ^2k^2c_\mathrm{g}^2\right)\left(\omega +\frac{ick^2}{3\kappa _\mathrm{F}\rho }\right)\left\{\omega ^2+\frac{4i\omega _\mathrm{a}E(\gamma 1)}{\gamma p}\omega \left[1\left(\widehat{𝐤}\widehat{𝐳}\right)^2\right]N_\mathrm{g}^2\right\}.$$
(54)
From left to right, the three factors respectively correspond to adiabatic acoustic waves in the gas, a purely damped radiative diffusion mode, and gravity waves in the gas modified by damping of gas temperature fluctuations by radiative emission and absorption. Note that radiation pressure and radiative buoyancy have been lost due to the rapid radiative diffusion at short wavelengths.

The behavior and character of the gravity waves depends on the relative magnitude of the Brunt-Väisälä frequency and a characteristic gravity wave thermal coupling frequency

$$\omega _{\mathrm{thg}}\frac{2\omega _\mathrm{a}E(\gamma 1)}{\gamma p}.$$
(55)
Solving equation (54) for the gravity mode frequencies,

$$\omega =i\omega _{\mathrm{thg}}\pm \left\{\omega _{\mathrm{thg}}^2+N_\mathrm{g}^2\left[1(\widehat{𝐤}\widehat{𝐳})^2\right]\right\}^{1/2}.$$
(56)
If $`N_\mathrm{g}^2[1(\widehat{𝐤}\widehat{𝐳})^2]>\omega _{\mathrm{thg}}^2`$, then the gravity waves are damped. If $`0<N_\mathrm{g}^2[1(\widehat{𝐤}\widehat{𝐳})^2]<\omega _{\mathrm{thg}}^2`$, then this damping is so strong that the modes lose their wavelike character and have purely negative imaginary frequencies. If $`N_\mathrm{g}^2<0`$, then one of the gravity modes becomes convectively unstable, but with a growth rate that is reduced compared to ideal hydrodynamic convection. This reduction is due to the fact that emission and absorption damp gas temperature fluctuations in the wave, relative to the nearly uniform radiation temperature.

The acoustic waves are more interesting. Expanding equation (53) to first order about the high $`k`$ limit, we find

$`\omega =\pm kc_\mathrm{g}i{\displaystyle \frac{\kappa _\mathrm{F}}{2cc_\mathrm{g}}}\left\{{\displaystyle \frac{4Ec_\mathrm{g}}{3}}+{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho c_\mathrm{g}}}\left(\widehat{𝐤}𝐅\right)\left[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\right\}+𝒪(k^1).`$ (57)
The first two terms within curly brackets represent damping by radiative diffusion and emission/absorption, respectively. The third term, which depends on the equilibrium radiative flux $`𝐅`$, can give rise to instability if it dominates the first two. For a Kramers type opacity law in a gas with $`\gamma =5/3`$, $`\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}<0`$, and equation (57) then implies that downward propagating sound waves are potentially unstable, while upward propagating waves are damped. On the other hand, if there is no opacity perturbation, $`\mathrm{\Theta }_\rho =\mathrm{\Theta }_{T\mathrm{g}}=0`$, and there is no acoustic instability. Such is the case for a medium with pure Thomson scattering opacity.

Instability requires that the driving term dominate the two damping terms in equation (57), so that a rough, order of magnitude instability criterion is

$$F\stackrel{~}{\mathrm{\Theta }}\mathrm{}>E\times \mathrm{max}[c_\mathrm{g},\left(\frac{\omega _\mathrm{a}}{\kappa _\mathrm{F}\rho c}\right)\left(\frac{c^2}{c_\mathrm{g}}\right)],$$
(58)
where $`\stackrel{~}{\mathrm{\Theta }}|\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}|`$. The first part of this criterion has an intuitive interpretation: if $`\stackrel{~}{\mathrm{\Theta }}`$ is of order unity, short wavelength acoustic waves will be unstable if the radiative flux is transporting the local radiation energy density faster than the sound speed in the gas.

As we will see in section 3.3 below, the maximum growth rate of the instability, when it exists, occurs for wavenumbers such that the first term in equation (57) exceeds the second, i.e. for $`k\mathrm{}>\kappa _\mathrm{F}F\stackrel{~}{\mathrm{\Theta }}/(cc_\mathrm{g}^2)`$. In a radiation pressure dominated medium, this implies that $`k\mathrm{}>\stackrel{~}{\mathrm{\Theta }}g/c_\mathrm{g}^2`$, i.e. wavelengths shorter than the gas scale height will have maximal growth rates. This growth rate itself will be $`\stackrel{~}{\mathrm{\Theta }}(g/c_\mathrm{g})`$, faster than the reciprocal of the local free fall time by the ratio $`\stackrel{~}{\mathrm{\Theta }}c_\mathrm{r}/c_\mathrm{g}`$.

On the other hand, in a gas pressure dominated medium, the growth rate is $`\stackrel{~}{\mathrm{\Theta }}E/(H_T\rho c_\mathrm{g})`$, where $`H_T`$ is the temperature scale height. Taking this to be comparable to the pressure scale height $`c_\mathrm{g}^2/g`$, we find a much smaller growth rate $`\stackrel{~}{\mathrm{\Theta }}(g/c_\mathrm{g})(E/p)`$, slower than the reciprocal of the local free fall time by the ratio $`\stackrel{~}{\mathrm{\Theta }}(c_\mathrm{r}/c_\mathrm{g})^2`$. This growth rate occurs for wavenumbers $`\mathrm{}>\stackrel{~}{\mathrm{\Theta }}(g/c_\mathrm{g}^2)(E/p)`$, a threshold which is in violation of our WKB requirement that $`kg/c_\mathrm{g}^2`$ for gas pressure dominated media. If the gas and radiation temperatures are not tightly locked, it appears difficult to produce vigorous acoustic wave instabilities in gas pressure dominated media by this mechanism.

An expression similar to equation (57) appears to have been first derived by Hearn (1972; his equations 25, 26 and 30 and surrounding discussion), the only difference being that he does not have the radiative diffusion damping term \[the first term in curly brackets in equation (57)\]. This is because Hearn (1972) was interested in radiative amplification of acoustic waves in optically thin media, and therefore did not include the dynamical evolution of the radiation field in his analysis.333It is worth noting that Hearn’s (1972) analysis had an error in his treatment of the gas energy equation (his equation 8). In this equation the partial time derivatives should be Lagrangian derivatives. It turns out, however, that his results are correct because he also neglected gas pressure gradient contributions to the equilibrium hydrostatic balance. These two errors cancelled one another.

### 3.2 Short Wavelength Limit With $`T_\mathrm{g}=T_\mathrm{r}`$

If we first take the $`\omega _\mathrm{a}\mathrm{}`$ limit in equation (53), so that the gas and radiation are tightly thermally coupled, then a mode is eliminated from the dispersion relation which is now only fifth order, including the zero frequency mode. The fact that we have lost a mode makes sense, as this limit corresponds to replacing the two time dependent gas and radiation energy equations with a total energy equation and the time-independent condition that the gas and radiation temperatures be equal. On taking the short wavelength limit, the four modes with nonzero frequencies factor as follows:

$`0`$ $`=`$ $`\left(\omega ^2k^2c_\mathrm{i}^2\right)\left[\omega +{\displaystyle \frac{ick^2}{3\kappa _\mathrm{F}\rho }}\left({\displaystyle \frac{4(\gamma 1)E}{p+4(\gamma 1)E}}\right)\right]`$ (59)
$`\times `$ $`\left\{\omega +{\displaystyle \frac{i3\kappa _\mathrm{F}\rho }{ck^2}}\left[1\left(\widehat{𝐤}\widehat{𝐳}\right)^2\right]\left(1+{\displaystyle \frac{4E}{3p}}\right)\left[{\displaystyle \frac{\gamma p}{4(\gamma 1)E}}N_\mathrm{g}^2+{\displaystyle \frac{1}{3}}N_\mathrm{r}^2\right]\right\},`$
where

$$c_\mathrm{i}\left(\frac{p}{\rho }\right)^{1/2}=\frac{c_\mathrm{g}}{\gamma ^{1/2}}$$
(60)
is the isothermal sound speed in the gas.

These short wavelength modes resemble those of equation (54). Acoustic waves propagate at the isothermal sound speed, because the gas temperature is locked to a radiation temperature which is made nearly uniform by the rapid radiative diffusion. The two gravity modes have collapsed to a single mode in the last factor of equation (59). This mode has a purely imaginary frequency, and has therefore lost its wavelike character. It is unstable if

$$\frac{\gamma p}{4(\gamma 1)E}N_\mathrm{g}^2+\frac{1}{3}N_\mathrm{r}^2<0.$$
(61)
However, radiative diffusion strongly diminishes the growth rate of this instability in the short wavelength limit, with $`|\omega |k^2`$ as $`k\mathrm{}`$. This is because buoyancy is strongly suppressed: a perturbed parcel of gas always has the same density as its surroundings, because it has the same pressure and is forced to have the same temperature due to the rapid radiative diffusion and rapid thermal coupling with the radiation.

The acoustic waves are again unstable when first order corrections to their frequencies are made:

$$\omega =\pm kc_\mathrm{i}i\frac{\kappa _\mathrm{F}}{2cc_\mathrm{i}}\left(1+\frac{3p}{4E}\right)\left[\left(\frac{4E}{3}+p\right)c_\mathrm{i}\left(\widehat{𝐤}𝐅\right)\mathrm{\Theta }_\rho \right]+𝒪(k^1).$$
(62)
The first term inside square brackets again represents damping by radiative diffusion, while the second term will drive instability if it is larger. For, e.g., Kramers law type opacities, $`\mathrm{\Theta }_\rho >0`$ and now the upward propagating wave is unstable while the downward propagating wave is damped. Once again, a pure Thomson scattering medium exhibits no acoustic wave instability.

An order of magnitude instability criterion from equation (62) is

$$F\mathrm{\Theta }_\rho \mathrm{}>\mathrm{max}[E,p]c_\mathrm{i},$$
(63)
which should be contrasted with equation (58). Note that there is no threshold from emission and absorption because the gas and radiation temperatures are equal. Here, acoustic waves are unstable if the local radiative flux is transporting the local thermal energy density, whether it be dominated by gas or radiation, faster than the gas sound speed.

In a radiation pressure dominated equilibrium, equation (62) implies that the growth rate of the acoustic wave instability is $`\kappa _\mathrm{F}F\mathrm{\Theta }_\rho /(cc_\mathrm{i})\mathrm{\Theta }_\rho g/c_\mathrm{i}`$ from hydrostatic equilibrium. In a gas pressure dominated medium, the growth rate is $`\kappa _\mathrm{F}F\rho c_\mathrm{i}\mathrm{\Theta }_\rho /(cE)c_\mathrm{i}\mathrm{\Theta }_\rho /H_T`$. Hence the growth rate in a gas pressure dominated medium is also $`\mathrm{\Theta }_\rho g/c_\mathrm{i}`$. In contrast to the previous case where gas and radiation did not exchange heat rapidly, radiation pressure support is not required to obtain high instability growth rates when the gas and radiation temperatures are the same. In a gas pressure dominated equilibrium, the small radiation pressure fluctuations produced by the damping and driving forces produce gas pressure fluctuations that are $`(p/E)`$ times larger just by the fact that the gas and radiation temperatures are locked together by rapid absorption and emission.

Hearn (1972; equations 22 and 23) was also the first to derive an expression similar to equation (62), although again without the first damping term because he was interested in optically thin media. In addition, Hearn’s (1972) approximations resulted in a replacement of the factor $`(1+0.75p/E)`$ multiplying the damping and growth rates in equation (62) with unity. This factor becomes very important in gas pressure dominated media.

We note that our result that pure Thomson scattering media exhibit no local hydrodynamic acoustic wave instability disagrees with the conclusion of Shaviv (2001), who claimed that his “Type II” instability represents just such a local instability. In fact, this instability appears to have rather long vertical wavelengths. The growth rates actually depend on the boundary conditions , and also the location of the boundary (Glatzel 2003, private communication), indicating that it is global in nature, not local.

### 3.3 The Limit of Negligible Gas Pressure

A number of authors have studied versions of these instabilities in the limit where gas pressure is completely negligible compared to radiation pressure. The gas sound speed then vanishes, and one then loses the acoustic wave nature of the instability in the short wavelength limit.

Setting $`p=0=c_\mathrm{g}^2`$ in the hydrodynamic dispersion relation (53), and then adopting the ansatz that $`\omega k^{1/2}`$ as $`k\mathrm{}`$ (Glatzel 1994, Gammie 1998), it is straightforward to show that there are two modes with

$$\omega ^2=\frac{i\kappa _\mathrm{F}}{c}𝐤𝐅\mathrm{\Theta }_\rho =i𝐤𝐠\mathrm{\Theta }_\rho .$$
(64)
Assuming $`\mathrm{\Theta }_\rho >0`$, one of these modes corresponds to an upward propagating unstable wave, while the other corresponds to a downward propagating wave which is damped. It is interesting to note that this result agrees with the thermally locked acoustic wave frequency, equation (62), if one squares that frequency and then takes the limit of negligible gas pressure. Note that the damping terms, which set the threshold of instability, vanish in the limit of negligible gas pressure.

If we first assume negligible thermal coupling between the gas and the radiation ($`\omega _\mathrm{a}0`$), and then consider the zero gas pressure limit, we get instead two modes that reflect the form of the driving in the two temperature regime:

$$\omega ^2=\frac{i\kappa _\mathrm{F}}{c}𝐤𝐅\left[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]=i𝐤𝐠\left[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right].$$
(65)
Once again, the growth rate exhibits a $`k^{1/2}`$ dependence. Note that the damping terms again vanish in this particular limit.

The acoustic wave instabilities we have been discussing throughout this section are the local, WKB versions of the strange mode instability discussed in the stellar oscillation literature. Glatzel (1994) has discussed a physical origin of strange modes in which he presents a WKB analysis of the growth rates of purely radial modes (i.e. vertical modes in our geometry) in the limit of zero gas pressure. His equation (5.8) is (nearly444Glatzel’s equation 5.8 is $`\omega ^2=\pm i2g\mathrm{\Theta }_\rho k`$, and we actually only recover the version of this equation with a plus sign. In addition, we do not have his factor of two. The reason for both of these facts is that his definition of the wavenumber $`k`$ differs from ours. Glatzel performed a WKB analysis on a variable $`\mathrm{\Pi }`$ which was a transformation of the Eulerian pressure perturbation (his equation 5.4). Our analysis only recovers modes with short wavelengths in the pressure perturbation, and this is true of only one sign of his dispersion relation. Moreover, our wavenumber for these modes is twice his wavenumber.) identical to our equation (64). He did not recover the damping terms of equation (62) for two reasons: he assumed negligible gas pressure, and he also invoked the NAR approximation, a point that we discuss further below in section 4.1. Glatzel (1994) also presented a global analytic solution for unstable strange modes in which he showed that unstable and damped waves propagate in opposite directions (see discussion after his equation 5.17), in agreement with our conclusions.

### 3.4 Numerical Results

Table 1 summarizes the basic conclusions of the analytic work on acoustic wave instabilities presented in the previous two subsections. We now discuss numerical solutions to the full dispersion relation (53) which show how these instabilities operate in different regimes, illustrating the basic formulas in Table 1.

At sufficiently short wavelengths, photons diffuse across a wavelength faster than the wave period, and radiation temperature fluctuations therefore become small. Whether or not the gas and radiation temperatures are locked together is therefore determined by whether absorption and emission are rapid enough to drive the gas temperature fluctuations to be small as well. From the perturbed gas energy equation (33), it is clear that the gas temperature will be locked to the radiation temperature provided the wave frequency $`\omega `$ is much less than a characteristic thermal frequency

$$\omega _{\mathrm{th}}\left[\frac{4(\gamma 1)E}{p}\right]\omega _\mathrm{a}.$$
(66)
Figure 1 depicts unstable wave growth rates as a function of wavenumber for a radiation pressure dominated equilibrium, for different, small values of $`\omega _{\mathrm{th}}`$. Although such equilibria are generally dominated by Thomson scattering opacity, for purposes of illustrating the physics, we have assumed Kramers type values for the flux mean opacity derivatives: $`\mathrm{\Theta }_\rho =1`$ and $`\mathrm{\Theta }_{T\mathrm{g}}=3.5`$. In any real application to a radiation pressure dominated equilibrium, these derivatives will be reduced approximately by the ratio of the absorption opacity $`\kappa _\mathrm{a}\omega _\mathrm{a}/(\rho c)`$ to $`\kappa _\mathrm{T}`$, and this would lead to a reduction in the unstable growth rates, or possibly stabilization. We present results appropriate for more realistic applications later on below.

At the low values of $`\omega _{\mathrm{th}}`$ shown in Figure 1, gas and radiation temperatures are decoupled at high wavenumber, and the asymptotic growth rates agree with equation (57). These unstable waves propagate downward because $`\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}<0`$. As $`\omega _{\mathrm{th}}`$ increases, their asymptotic growth rates decrease and eventually damp. At the same time, a new set of modes appear at low wavenumber with growth rates that climb with increasing thermal coupling between the gas and radiation, as shown more clearly in Figure 1(b). These modes propagate upward and at higher thermal coupling turn into the one-temperature acoustic waves described by equation (62) at high wavenumber.

This is shown in Figure 2(a), which depicts the growth rates as a function of wavenumber for high values of $`\omega _{\mathrm{th}}`$, for the same equilibrium parameters used in Figure 1. In contrast to the two-temperature behavior, the asymptotic growth rates at high wavenumber must eventually fail because at short enough wavelengths the modes must return to being adiabatic in the gas alone, and the assumption of thermal locking of the gas and radiation underlying equation (62) must break down. Figure 2(b) illustrates that for wavenumbers $`k\mathrm{}>\omega _{\mathrm{th}}/c_\mathrm{g}`$, the real part of the phase velocity returns to the adiabatic gas sound speed.

The one temperature instability growth rate cuts off at even lower wavenumbers than this, however, and the reason appears to be a breakdown of thermal locking related to an interplay between the equilibrium radiation flux, absorption and emission, and fluid inertia. Numerical exploration shows that the growth rate cuts off for wavenumbers $`k\mathrm{}>(\omega _{\mathrm{th}}g|1+\mathrm{\Theta }_\rho |/c_\mathrm{i}^3)^{1/2}`$ for both gas and radiation pressure dominated equilibria, and this formula for the maximum cutoff wavenumber is included in Table 1.

For demonstration purposes, we have been considering equilibria with rather artificial opacities. Figure 3 shows the unstable wave growth rate as a function of wavenumber for a gas-pressure dominated equilibrium with a realistic Kramers type opacity law, appropriate for mid- to upper main sequence stars. For sufficiently high thermal coupling, as is indeed present in these stars (see section 7.2 below), unstable upward propagating waves exist over a broad range of high wavenumbers. Interestingly, long wavelength instabilities also appear in downward propagating waves. These downward wave instabilities exist because of the temperature dependence of the opacity ($`\mathrm{\Theta }_{T\mathrm{g}}=3.5`$). Because the gas and radiation temperatures are tightly locked, temperature fluctuations in the gas are smoothed out by radiative diffusion. Hence opacity fluctuations caused by $`\delta T`$ are smaller than opacity fluctuations caused by $`\delta \rho `$. Therefore the temperature dependence of the opacity does not affect the growth rate at short wavelengths. This agrees with equation (62). However, at long wavelengths where radiative diffusion is slow, significant gas temperature fluctuations exist in the wave, and the resulting opacity fluctuations reverse the radiative driving. Note that these long wavelength, downward unstable waves do not appear to exist in radiation pressure dominated equilibria, as shown in Figure 2.

## 4 Physics of the Acoustic Wave Instabilities

The results of the previous section are based on a brute force analytic and numerical attack on the full dispersion relation. In this section, we present a more physically-motivated derivation of some of these results. (This turns out to be a useful starting point for understanding the unstable MHD waves that we discuss below in sections 5 and 6.) It is the coupling of the equilibrium radiation flux to opacity fluctuations that drives acoustic waves unstable. To see how this works physically, consider the forces acting on a perturbed parcel of gas in the acoustic wave. \[This approach is somewhat analogous to using a zeroth order eigenfunction to compute damping or driving rates through a work integral in stellar oscillation theory (Cox 1980).\]

The gas and radiation energy equations (33) and (36) may be written in the WKB limit as

$$\delta p=c_\mathrm{g}^2\delta \rho \frac{i}{\omega }\delta 𝐯\left(\text{r}pc_\mathrm{g}^2\text{r}\rho \right)+\frac{4i(\gamma 1)\omega _\mathrm{a}E}{\omega }\left(\frac{\delta T_\mathrm{r}\delta T_\mathrm{g}}{T}\right)$$
(67)
and

$`\left(1+{\displaystyle \frac{ick^2}{3\kappa _\mathrm{F}\rho \omega }}\right){\displaystyle \frac{\delta E}{3}}`$ $`=`$ $`c_\mathrm{r}^2\delta \rho {\displaystyle \frac{i\rho T}{3\omega }}\delta 𝐯\text{r}S_\mathrm{r}{\displaystyle \frac{4i\omega _\mathrm{a}E}{3\omega }}\left({\displaystyle \frac{\delta T_\mathrm{r}\delta T_\mathrm{g}}{T}}\right)`$ (68)
$``$ $`{\displaystyle \frac{𝐤𝐅}{3\omega }}\left[(1+\mathrm{\Theta }_\rho ){\displaystyle \frac{\delta \rho }{\rho }}+\mathrm{\Theta }_{T\mathrm{r}}{\displaystyle \frac{\delta T_\mathrm{r}}{T}}+\mathrm{\Theta }_{T\mathrm{g}}{\displaystyle \frac{\delta T_\mathrm{g}}{T}}\right],`$
respectively. The radiation energy equation (68) immediately implies that for short wavelength ($`k\mathrm{}`$) acoustic waves, $`\delta E𝒪(k^1)\delta \rho `$. Physically, at high $`k`$ the radiation diffusion time across a wavelength is shorter than a wave period, so that radiation temperature fluctuations are smoothed out. On the other hand, because acoustic waves in the gas are adiabatic in this limit, the gas temperature fluctuation is given by $`\delta T_\mathrm{g}/T=(\gamma 1)\delta \rho /\rho `$. Hence the gas pressure perturbation is given by

$$\delta p=c_\mathrm{g}^2\delta \rho \frac{i}{\omega }\delta 𝐯\left(\text{r}pc_\mathrm{g}^2\text{r}\rho \right)\frac{4i(\gamma 1)^2\omega _\mathrm{a}E}{\omega }\left(\frac{\delta \rho }{\rho }\right)+𝒪(k^2).$$
(69)
The first term on the right hand side is the dominant term at short wavelengths, and represents the acoustic response of the gas to compressive perturbations. The second term arises from buoyancy. The last term is a damping term caused by radiative emission/absorption: temperature fluctuations in the gas emit and absorb photons to try and come into equilibrium with the nearly uniform radiation temperature.

To the same order in $`k^1`$, the radiation pressure perturbation is

$$\frac{\delta E}{3}=\frac{i\kappa _\mathrm{F}}{c}\left[\frac{4E}{3}\left(\frac{\omega }{k^2}\right)\right]\delta \rho +\frac{i\kappa _\mathrm{F}}{k^2c}\left(𝐤𝐅\right)\left[1+\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\delta \rho +𝒪(k^2).$$
(70)
The first term on the right hand side represents damping by radiative diffusion: a density fluctuation will try to compress photons, but these diffuse quickly and transport heat out of the fluctuation. The last term is the most interesting. It is an extra radiation pressure force arising from the equilibrium radiation flux acting on density and opacity fluctuations in the gas. As we shall see, it is this extra force that is destabilizing if it can overcome the damping forces.

Using equations (69)-(70) to eliminate the total pressure perturbation, the gas momentum equation (47) for short wavelength acoustic perturbations may be written

$`\rho {\displaystyle \frac{\delta 𝐯}{t}}`$ $`=`$ $`i𝐤c_\mathrm{g}^2\delta \rho {\displaystyle \frac{𝐤}{\omega }}\delta 𝐯\left(\text{r}pc_\mathrm{g}^2\text{r}\rho \right){\displaystyle \frac{\kappa _\mathrm{F}}{c}}𝐤\left[{\displaystyle \frac{4E}{3}}\left({\displaystyle \frac{\omega }{k^2}}\right)\right]\delta \rho {\displaystyle \frac{\kappa _\mathrm{F}}{c}}𝐤\left[{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho \omega }}\right]\delta \rho `$ (71)
$`+`$ $`{\displaystyle \frac{\kappa _\mathrm{F}}{k^2c}}𝐤\left(𝐤𝐅\right)\left[1+\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\delta \rho +𝐠\delta \rho +𝒪(k^1)\delta \rho .`$
The destabilizing force density is $`𝐤(𝐤𝐅)\delta \rho `$ times a numerical factor depending on opacity derivatives that may be positive or negative. For the purposes of discussion, let us assume for the moment that it is positive. Figure 4(a) illustrates the geometry of this force acting on a downward propagating acoustic wave. In this case $`𝐤(𝐤𝐅)\delta \rho `$ is everywhere opposite to the velocity perturbation, and the wave is therefore damped. However, for an upward propagating wave, shown in Figure 4(b), the force is everywhere in phase with the velocity perturbation, and therefore amplifies the wave. If we had chosen the numerical factor depending on opacity derivatives to be negative, then it would have been the upward propagating wave that is damped and the downward propagating wave that is amplified.

It is also possible to consider the effects of the destabilizing force density in terms of whether it induces a net lag or lead between the total pressure and density perturbations in the sound wave. Once again, only one direction of propagation is destabilizing when viewed in this way.

Examination of the growth rate in equation (57) reveals that it is just the combination of opacity derivatives $`\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}`$ that determines the action of the force, not the full factor $`[1+\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}]`$ in equation (71). Pure Thomson scattering media are not unstable by this mechanism. To see why this is so, it is necessary to examine the equation of motion (71) more closely. Using the perturbed continuity equation (46) and regrouping, equation (71) may be rewritten as

$`\rho {\displaystyle \frac{\delta 𝐯}{t}}`$ $`=`$ $`{\displaystyle \frac{i𝐤c_\mathrm{g}^2\rho }{\omega }}\left(𝐤\delta 𝐯\right){\displaystyle \frac{\kappa _\mathrm{F}}{c}}𝐤\left[{\displaystyle \frac{4E}{3}}\left({\displaystyle \frac{\omega }{k^2}}\right)\right]\delta \rho {\displaystyle \frac{\kappa _\mathrm{F}}{c}}𝐤\left[{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho \omega }}\right]\delta \rho `$ (72)
$`+`$ $`{\displaystyle \frac{\kappa _\mathrm{F}}{k^2c}}𝐤\left(𝐤𝐅\right)\left[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\delta \rho `$
$`+`$ $`\left[{\displaystyle \frac{\kappa _\mathrm{F}}{k^2c}}𝐤\left(𝐤𝐅\right)\delta \rho {\displaystyle \frac{𝐤}{\omega }}\delta 𝐯\text{r}p+𝐠\delta \rho \right]+𝒪(k^1)\delta \rho .`$
It is convenient to rewrite this equation in terms of the Lagrangian displacement ¸, where

$$\delta 𝐯=\frac{\text{¸}}{t}=i\omega \text{¸}.$$
(73)
Using the fact that to lowest order in $`k^1`$, $`\delta \rho =i\rho 𝐤\text{¸}`$ from the continuity equation (46), we obtain

$`{\displaystyle \frac{^2\text{¸}}{t^2}}`$ $`=`$ $`𝐤c_\mathrm{g}^2\left(𝐤\text{¸}\right)`$ (74)
$`+`$ $`i(𝐤\text{¸})𝐤{\displaystyle \frac{\kappa _\mathrm{F}}{c}}\left\{{\displaystyle \frac{4E}{3}}\left({\displaystyle \frac{\omega }{k^2}}\right)+{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho \omega }}{\displaystyle \frac{1}{k^2}}(𝐤𝐅)\left[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\right\}`$
$``$ $`i\left[{\displaystyle \frac{\kappa _\mathrm{F}}{k^2c}}𝐤\left(𝐤𝐅\right)\left(𝐤\text{¸}\right){\displaystyle \frac{1}{\rho }}𝐤\text{¸}\text{r}p+𝐠\left(𝐤\text{¸}\right)\right]+𝒪(k^0)|\text{¸}|.`$
The destabilizing $`𝐤𝐅`$ term that would be present in Thomson scattering atmospheres competes with the equilibrium gravitational acceleration and gas pressure gradient in the last term in square brackets in equation (74). Figure 4 illustrates this by showing that the projection of $`𝐠\delta \rho `$ along the wave acts in the opposite direction to $`𝐤\mathrm{𝐅𝐤}\delta \rho `$. Eliminating the gas pressure gradient using the hydrostatic equilibrium equation (21), the last term in square brackets in equation (74) may be written

$$\frac{\kappa _\mathrm{F}}{k^2c}𝐤\left(𝐤𝐅\right)\left(𝐤\text{¸}\right)𝐤\text{¸}\text{r}p+𝐠\left(𝐤\text{¸}\right)=\frac{\kappa _\mathrm{F}}{c}𝐤(\widehat{𝐤}\times 𝐅)(\widehat{𝐤}\times \text{¸})\text{¸}\times (𝐤\times 𝐠).$$
(75)
Because hydrodynamic acoustic waves are longitudinal to lowest order in $`k^1`$, $`\text{¸}𝐤`$, and the term involving the flux vanishes. In addition, the second term does no work as it is perpendicular to the displacement. Hence these terms contribute nothing to the growth rate - the destabilizing radiation pressure force that would be present in a Thomson scattering medium is exactly cancelled due to the hydrostatic equilibrium of the unperturbed flow, at least for longitudinal waves. In the high-$`k`$ limit, the modes with density perturbations that can couple to the equilibrium radiation force are necessarily longitudinal. This is why gravity waves do not experience radiative driving in the short wavelength limit.

Taking the scalar product of equation (74) with the wave vector $`𝐤`$, we obtain

$`{\displaystyle \frac{^2}{t^2}}(𝐤\text{¸})`$ $`=`$ $`k^2c_\mathrm{g}^2(𝐤\text{¸})+ik^2(𝐤\text{¸}){\displaystyle \frac{\kappa _\mathrm{F}}{c}}\{{\displaystyle \frac{4E}{3}}\left({\displaystyle \frac{\omega }{k^2}}\right)+{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho \omega }}`$ (76)
$``$ $`{\displaystyle \frac{1}{k^2}}(𝐤𝐅)[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}]\}+𝒪(k)|\text{¸}|.`$
Replacing $`/t`$ with $`i\omega `$, we immediately recover equation (57) for the mode frequency and growth rate.

We may repeat this procedure to better understand the instability in the limit where the gas and radiation exchange heat extremely rapidly. The gas and radiation temperatures are then locked together, so that

$$\delta T_\mathrm{g}=\delta T_\mathrm{r}\delta T=\frac{T}{4E}\delta E=T\left(\frac{\delta p}{p}\frac{\delta \rho }{\rho }\right).$$
(77)
In this regime we must replace the separate gas and radiation energy equations with a single, combined energy equation which, when perturbed, gives

$`i\omega \left({\displaystyle \frac{\delta p}{\gamma 1}}+\delta E\right)`$ $`=`$ $`{\displaystyle \frac{k^2c}{3\kappa _\mathrm{F}\rho }}\delta Ei\omega \left({\displaystyle \frac{c_\mathrm{g}^2}{\gamma 1}}+3c_\mathrm{r}^2\right)\delta \rho \rho T\delta 𝐯\text{r}(S_\mathrm{g}+S_\mathrm{r})`$ (78)
$`+`$ $`i𝐤𝐅\left[(1+\mathrm{\Theta }_\rho ){\displaystyle \frac{\delta \rho }{\rho }}+(\mathrm{\Theta }_{T\mathrm{g}}+\mathrm{\Theta }_{T\mathrm{r}}){\displaystyle \frac{\delta T}{T}}\right].`$
Once again, equation (78) immediately implies that $`\delta E𝒪(k^1)\delta \rho `$ for short wavelength acoustic waves. Rapid radiative diffusion smooths out both gas and radiation temperature fluctuations in this case. Equation (77) then implies that $`\delta p=c_\mathrm{i}^2\delta \rho `$ to lowest order, i.e. that the acoustic waves are isothermal in this limit. Equation (78) can then be solved for the radiation pressure perturbation,

$$\frac{\delta E}{3}=\frac{i\kappa _\mathrm{F}}{c}\left(p+\frac{4E}{3}\right)\left(\frac{\omega }{k^2}\right)\delta \rho +\frac{i\kappa _\mathrm{F}}{k^2c}(𝐤𝐅)(1+\mathrm{\Theta }_\rho )\delta \rho +𝒪(k^2)\delta \rho .$$
(79)
From equation (77), the gas pressure perturbation is given to the same order by

$$\delta p=c_\mathrm{i}^2\delta \rho \frac{i\kappa _\mathrm{F}}{c}\left(1+\frac{3p}{4E}\right)p\left(\frac{\omega }{k^2}\right)\delta \rho +\frac{i\kappa _\mathrm{F}}{k^2c}\left(\frac{3p}{4E}\right)(𝐤𝐅)(1+\mathrm{\Theta }_\rho )\delta \rho +𝒪(k^2)\delta \rho .$$
(80)
The first term on the right hand side of equation (79) and the second term on the right hand side of equation (80) are the damping terms due to radiative diffusion. The last terms in both equations are the potentially destabilizing terms from the equilibrium radiation flux acting on density and opacity fluctuations. Note that the opacity derivative with respect to temperature, $`\mathrm{\Theta }_{T\mathrm{g}}`$, has been lost because the gas temperature is locked to the uniform radiation temperature. There is also no buoyancy response in the pressure perturbations to this order, as expected because of the suppression of gravity waves in this limit.

Using equations (79)-(80), the gas momentum equation (47) becomes

$`\rho {\displaystyle \frac{\delta 𝐯}{t}}`$ $`=`$ $`i𝐤c_\mathrm{i}^2\delta \rho {\displaystyle \frac{\kappa _\mathrm{F}}{c}}𝐤\left(1+{\displaystyle \frac{3p}{4E}}\right)\left({\displaystyle \frac{4E}{3}}+p\right)\left({\displaystyle \frac{\omega }{k^2}}\right)\delta \rho `$ (81)
$`+`$ $`{\displaystyle \frac{\kappa _\mathrm{F}}{k^2c}}𝐤\left(𝐤𝐅\right)\left(1+{\displaystyle \frac{3p}{4E}}\right)\left(1+\mathrm{\Theta }_\rho \right)\delta \rho +𝐠\delta \rho +𝒪(k^1)\delta \rho .`$
Once again, the destabilizing force density term is proportional to $`𝐤(𝐤𝐅)\delta \rho `$, and the geometry of Figure 4 still applies. Repeating the same steps as before, we obtain

$`{\displaystyle \frac{^2\text{¸}}{t^2}}`$ $`=`$ $`𝐤c_\mathrm{i}^2(𝐤\text{¸})`$ (82)
$`+`$ $`i(𝐤\text{¸})𝐤{\displaystyle \frac{\kappa _\mathrm{F}}{c}}\left(1+{\displaystyle \frac{3p}{4E}}\right)\left[\left({\displaystyle \frac{4E}{3}}+p\right){\displaystyle \frac{\omega }{k^2}}{\displaystyle \frac{1}{k^2}}(𝐤𝐅)\mathrm{\Theta }_\rho \right]`$
$``$ $`i\left[\left(1+{\displaystyle \frac{3p}{4E}}\right){\displaystyle \frac{\kappa _\mathrm{F}}{k^2c}}𝐤(𝐤𝐅)(𝐤\text{¸})𝐤c_\mathrm{i}^2\text{¸}\text{r}\mathrm{ln}\rho +𝐠(𝐤\text{¸})\right]+𝒪(k^0)|\text{¸}|,`$
which should be compared with equation (74). The destabilizing $`𝐤𝐅`$ term that would remain in a pure Thomson scattering atmosphere again competes with other forces in the last term in square brackets. Using the equilibrium LTE condition (20), the hydrostatic equilibrium equation (21), and the radiative diffusion equation (23), the equilibrium density gradient is given by

$$c_\mathrm{i}^2\text{r}\mathrm{ln}\rho =𝐠+\left(1+\frac{3p}{4E}\right)\frac{\kappa _\mathrm{F}}{c}𝐅.$$
(83)
The last term in square brackets in equation (82) is therefore

$$\left(1+\frac{3p}{4E}\right)\frac{\kappa _\mathrm{F}}{k^2c}𝐤(𝐤𝐅)(𝐤\text{¸})𝐤c_\mathrm{i}^2\text{¸}\text{r}\mathrm{ln}\rho +𝐠(𝐤\text{¸})=\left(1+\frac{3p}{4E}\right)\frac{\kappa _\mathrm{F}}{c}𝐤(\widehat{𝐤}\times 𝐅)(\widehat{𝐤}\times \text{¸})\text{¸}\times (𝐤\times 𝐠).$$
(84)
Once again, these terms do not affect the growth rate for longitudinal waves. After taking the scalar product of equation (82) with $`𝐤`$, we finally obtain

$`{\displaystyle \frac{^2}{t^2}}\left(𝐤\text{¸}\right)`$ $`=`$ $`k^2c_\mathrm{i}^2(𝐤\text{¸})+ik^2(𝐤\text{¸}){\displaystyle \frac{\kappa _\mathrm{F}}{c}}(1+{\displaystyle \frac{3p}{4E}})[({\displaystyle \frac{4E}{3}}+p){\displaystyle \frac{\omega }{k^2}}`$ (85)
$``$ $`{\displaystyle \frac{1}{k^2}}(𝐤𝐅)\mathrm{\Theta }_\rho ]+𝒪(k)|\text{¸}|,`$
from which we recover the growth rate (62).

### 4.1 Relation to the NAR Approximation

Gautschy & Glatzel (1990) have shown that strange modes are well-described by a so-called non-adiabatic reversible (NAR) approximation. Mathematically, this approximation reduces to taking the scalar product of the perturbed radiative diffusion equation (32) with $`𝐤`$, and then setting $`𝐤\delta 𝐅=0`$. One also ignores the radiation energy equation (29). In addition, if the gas and radiation temperatures are locked together by rapid emission and absorption of photons, the gas energy equation (28) is also ignored. The perturbed radiation pressure is then given by

$$ik^2\frac{\delta E}{3}=\frac{\kappa _\mathrm{F}\rho }{c}𝐤𝐅\left(\frac{\delta \rho }{\rho }+\frac{\delta \kappa _\mathrm{F}}{\kappa _\mathrm{F}}\right).$$
(86)
This equation recovers our high $`k`$ expressions (70) and (79) for the perturbed radiation pressure, except that the NAR approximation cannot reproduce the radiation diffusion damping terms (as it neglects the radiation energy equation). Indeed, if we use equations (32), (70) and (79) to derive the perturbed radiative flux, we find

$$\delta 𝐅=\{\begin{array}{cc}3c_\mathrm{r}^2\left(\frac{\omega }{k^2}\right)𝐤\delta \rho +\left[\frac{𝐤(𝐤𝐅)}{k^2}𝐅\right]\left[1+\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\frac{\delta \rho }{\rho }+𝒪(k^1)\delta \rho \hfill & \text{if }\delta T_\mathrm{r}\delta T_\mathrm{g}\text{,}\hfill \\ (c_\mathrm{i}^2+3c_\mathrm{r}^2)\left(\frac{\omega }{k^2}\right)𝐤\delta \rho +\left[\frac{𝐤(𝐤𝐅)}{k^2}𝐅\right](1+\mathrm{\Theta }_\rho )\frac{\delta \rho }{\rho }+𝒪(k^1)\delta \rho \hfill & \text{if }\delta T_\mathrm{r}=\delta T_\mathrm{g}\text{.}\hfill \end{array}$$
(87)
When unstable driving (the second term on the right hand sides) strongly dominates radiative diffusion damping (the first term on the right hand sides), the NAR approximation that $`𝐤\delta 𝐅=0`$ is excellent.

By not invoking the NAR approximation, however, we are able to deduce the threshold criteria listed in Table 1 required for the unstable driving to exceed the damping by radiative diffusion. At the same time, we can also calculate the actual luminosity perturbation associated with the unstable wave.

## 5 MHD Instabilities in a Static Medium

We now extend the analysis of the previous section by including magnetic stresses. It turns out that magnetic fields widen the domain of the acoustic wave instabilities to include even pure Thomson scattering media. This is because magnetic tension endows acoustic waves with a mixed longitudinal/transverse character, except in the special cases of propagation purely along or perpendicular to the magnetic field.

### 5.1 Short Wavelength Limit

The dispersion relation (49) again factors easily in the short wavelength limit:

$`0`$ $`=`$ $`\left[\omega ^2(𝐤𝐯_\mathrm{A})^2\right]\left[\omega ^4\omega ^2k^2(v_\mathrm{A}^2+c_\mathrm{g}^2)+k^2c_\mathrm{g}^2(𝐤𝐯_\mathrm{A})^2\right]`$ (88)
$`\times `$ $`\left(\omega +{\displaystyle \frac{ick^2}{3\kappa _\mathrm{F}\rho }}\right)\left[\omega +{\displaystyle \frac{4i\omega _\mathrm{a}E(\gamma 1)}{\gamma p}}\right].`$
From left to right, the factors correspond to Alfvén waves, fast and slow magnetosonic waves, the damped radiative diffusion mode, and a damped $`\delta T_\mathrm{g}\delta T_\mathrm{r}`$ mode.555In our previous paper (BS01), we found only seven nonzero frequency modes in our analysis of the perturbed radiation MHD equations. This is because we did not include heat exchange between the gas and radiation through absorption and emission, a process which turns out to be important in most applications (see section 7). Setting $`\omega _\mathrm{a}=0`$ in equation (88) produces an additional zero frequency mode. The gravity waves of equation (54) have been lost because they are dominated by magnetic fields in this short wavelength limit, at least when $`𝐤𝐯_\mathrm{A}0`$, i.e. for wave vectors that are not orthogonal to the equilibrium magnetic field.

The first order corrections to the Alfvén modes vanish, so that $`\omega =\pm 𝐤𝐯_\mathrm{A}+𝒪(k^1)`$. Because Alfvén waves are purely transverse and thus do not involve density fluctuations in this short wavelength limit, they do not couple to the radiation physics.

To the same order, the magnetosonic wave frequencies are given by

$`\omega `$ $`=`$ $`\pm kv_{\mathrm{ph}}{\displaystyle \frac{i\kappa _\mathrm{F}}{2cv_{\mathrm{ph}}}}\left[{\displaystyle \frac{v_{\mathrm{ph}}^2(\widehat{𝐤}𝐯_A)^2}{2v_{\mathrm{ph}}^2v_\mathrm{A}^2c_\mathrm{g}^2}}\right]\left\{{\displaystyle \frac{4Ev_{\mathrm{ph}}}{3}}+{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho v_{\mathrm{ph}}}}\left(\widehat{𝐤}𝐅\right)\left[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\right\}`$ (89)
$`\pm {\displaystyle \frac{i\kappa _\mathrm{F}}{2cv_{\mathrm{ph}}(2v_{\mathrm{ph}}^2v_\mathrm{A}^2c_\mathrm{g}^2)}}(\widehat{𝐤}𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐅)+𝒪(k^1),`$
where $`v_{\mathrm{ph}}>0`$ is the phase speed of the fast and slow magnetosonic waves, given by

$$v_{\mathrm{ph}}^2\frac{1}{2}\left\{v_\mathrm{A}^2+c_\mathrm{g}^2\pm \left[(v_A^2+c_\mathrm{g}^2)^24(\widehat{𝐤}𝐯_\mathrm{A})^2c_\mathrm{g}^2\right]^{1/2}\right\}.$$
(90)
Equation (89) bears a strong resemblance to its hydrodynamic counterpart (57). The gas sound speed $`c_\mathrm{g}`$ has been replaced by the phase velocity of the relevant magnetosonic wave $`v_{\mathrm{ph}}`$, and the “hydrodynamic” corrections to the mode frequency have been multiplied by a factor

$$\left(\frac{v_{\mathrm{ph}}^2(\widehat{𝐤}𝐯_A)^2}{2v_{\mathrm{ph}}^2v_\mathrm{A}^2c_\mathrm{g}^2}\right)0.$$
(91)
The last term in equation (89) is a new, destabilizing term that exists only in the presence of magnetic fields. As we discuss in more detail in section 6 below, its physical origin ultimately arises from the ability of magnetic stresses to support velocity perturbations that are not purely longitudinal, i.e. that are not parallel or antiparallel to the wave vector $`𝐤`$. Note that this term vanishes for wave vectors that are either parallel or perpendicular to the equilibrium magnetic field. In each of these cases the compressive magnetosonic wave is longitudinal.

In the rapid heat exchange ($`\omega _\mathrm{a}\mathrm{}`$) limit, one of the eight modes is lost for the same reason as in the hydrodynamic case, and the remaining seven modes factor in the short wavelength limit as

$`0`$ $`=`$ $`\left[\omega ^2(𝐤𝐯_\mathrm{A})^2\right]\left[\omega ^4\omega ^2k^2(v_\mathrm{A}^2+c_\mathrm{i}^2)+k^2c_\mathrm{i}^2(𝐤𝐯_\mathrm{A})^2\right]`$ (92)
$`\times `$ $`\left[\omega +{\displaystyle \frac{ick^2}{3\kappa _\mathrm{F}\rho }}\left({\displaystyle \frac{4(\gamma 1)E}{p+4(\gamma 1)E}}\right)\right].`$
Once again, instabilities only occur in the fast and slow magnetosonic waves, whose frequencies are given to first order by

$`\omega `$ $`=`$ $`\pm kv_{\mathrm{ph}}{\displaystyle \frac{i\kappa _\mathrm{F}}{2cv_{\mathrm{ph}}}}\left[{\displaystyle \frac{v_{\mathrm{ph}}^2(\widehat{𝐤}𝐯_A)^2}{2v_{\mathrm{ph}}^2v_\mathrm{A}^2c_\mathrm{i}^2}}\right]\left(1+{\displaystyle \frac{3p}{4E}}\right)\left\{\left({\displaystyle \frac{4E}{3}}+p\right)v_{\mathrm{ph}}\left(\widehat{𝐤}𝐅\right)\mathrm{\Theta }_\rho \right\}`$ (93)
$`\pm {\displaystyle \frac{i\kappa _\mathrm{F}}{2cv_{\mathrm{ph}}(2v_{\mathrm{ph}}^2v_\mathrm{A}^2c_\mathrm{i}^2)}}\left(1+{\displaystyle \frac{3p}{4E}}\right)(\widehat{𝐤}𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐅)+𝒪(k^1),`$
where $`v_{\mathrm{ph}}`$ is given by equation (90) with $`c_\mathrm{g}`$ replaced by the isothermal sound speed $`c_\mathrm{i}`$.

Equation (90) implies that $`(2v_{\mathrm{ph}}^2v_\mathrm{A}^2c_\mathrm{g}^2)`$ is positive for fast modes and negative for slow modes. Hence if the last “magnetic” destabilizing term in equation (89), or (93), dominates the hydrodynamic terms, the unstable fast and slow waves will have opposite vector phase velocities. It is also easy to see that in this case the slow mode will always grow faster than the fast mode. However, if hydrodynamic driving (the opacity derivative terms in eqs.  and ) dominates magnetic driving, then whichever of the two modes has the larger density perturbation for a given velocity perturbation will be the one that grows fastest. If $`v_\mathrm{A}^2c_\mathrm{g}^2`$, the slow mode grows faster, whereas if $`v_\mathrm{A}^2c_\mathrm{g}^2`$, the fast mode grows faster.

### 5.2 Short Wavelength Limit With $`𝐤𝐯_\mathrm{A}=0`$: Recovery of Gravity Waves

If we consider wave vectors that are orthogonal to the equilibrium magnetic field, then the short wavelength limit recovers the fast magnetosonic modes, the radiative diffusion mode, and gives us a cubic equation describing coupled gravity and $`\delta T_\mathrm{g}\delta T_\mathrm{r}`$ modes:

$`0`$ $`=`$ $`(v_\mathrm{A}^2+c_\mathrm{g}^2)\omega ^3+{\displaystyle \frac{4i(\gamma 1)\omega _\mathrm{a}E}{p}}(v_\mathrm{A}^2+c_\mathrm{i}^2)\omega ^2\left[1(\widehat{𝐤}\widehat{𝐳})^2\right]\left(c_\mathrm{g}^2N_\mathrm{g}^2+v_\mathrm{A}^2𝐠\text{r}\mathrm{ln}\rho \right)\omega `$ (94)
$``$ $`\left[1(\widehat{𝐤}\widehat{𝐳})^2\right]{\displaystyle \frac{4i(\gamma 1)\omega _\mathrm{a}E}{p}}v_\mathrm{A}^2𝐠\text{r}\mathrm{ln}\rho .`$
The reason we recover gravity waves in this limit is that magnetic tension does not exist when $`𝐤𝐯_\mathrm{A}=0`$, and therefore no longer dominates buoyancy at short wavelengths. For weak thermal coupling between the gas and radiation, $`\omega _\mathrm{a}0`$, the roots of this dispersion relation correspond to magnetically altered gravity waves in the gas,

$$\omega ^2=\frac{1(\widehat{𝐤}\widehat{𝐳})^2}{v_\mathrm{A}^2+c_\mathrm{g}^2}\left(c_\mathrm{g}^2N_\mathrm{g}^2+v_\mathrm{A}^2𝐠\text{r}\mathrm{ln}\rho \right),$$
(95)
and a mode

$$\omega \frac{4i(\gamma 1)\omega _\mathrm{a}E}{p}\left(\frac{v_\mathrm{A}^2𝐠\text{r}\mathrm{ln}\rho }{c_\mathrm{g}^2N_\mathrm{g}^2+v_\mathrm{A}^2𝐠\text{r}\mathrm{ln}\rho }\right).$$
(96)
In the rapid heat exchange limit, $`\omega _\mathrm{a}\mathrm{}`$, the last mode is lost and the gravity waves depend only on the vertical density profile,

$$\omega ^2=\frac{1(\widehat{𝐤}\widehat{𝐳})^2}{v_\mathrm{A}^2+c_\mathrm{i}^2}v_\mathrm{A}^2𝐠\text{r}\mathrm{ln}\rho .$$
(97)
### 5.3 The Limit of Zero Gas Pressure: Photon Bubbles

Just as in the hydrodynamic case, a number of authors (Arons 1992, Gammie 1998) have investigated the slow mode instability in the limit of zero gas pressure. Because the phase speed of the slow mode then vanishes, the instability then loses its obvious connection to this mode.

Setting $`p=0=c_\mathrm{g}^2`$ in the dispersion relation (49), and then assuming $`\omega k^{1/2}`$ as $`k\mathrm{}`$, we find two modes with

$`\omega ^2`$ $`=`$ $`{\displaystyle \frac{i\kappa _\mathrm{F}}{cv_\mathrm{A}^2}}𝐤𝐯_\mathrm{A}\left[(\widehat{𝐤}𝐯_\mathrm{A})(\widehat{𝐤}𝐅)\mathrm{\Theta }_\rho (\widehat{𝐤}\times 𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐅)\right]`$ (98)
$`=`$ $`{\displaystyle \frac{i}{v_\mathrm{A}^2}}𝐤𝐯_\mathrm{A}\left[(\widehat{𝐤}𝐯_\mathrm{A})(\widehat{𝐤}𝐠)\mathrm{\Theta }_\rho (\widehat{𝐤}\times 𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐠)\right].`$
Just as in the hydrodynamic case \[eq. (64)\], equation (98) can be obtained from the thermally locked slow mode frequency in equation (93) by first squaring the frequency and then taking the limit of zero gas pressure. Again, we lose the damping terms in this limit, which determine the instability threshold for short wavelength modes.

If we first assume negligible thermal coupling between the gas and radiation, and then take the limit of zero gas pressure, we obtain instead

$`\omega ^2`$ $`=`$ $`{\displaystyle \frac{i\kappa _\mathrm{F}}{cv_\mathrm{A}^2}}𝐤𝐯_\mathrm{A}\left\{(\widehat{𝐤}𝐯_\mathrm{A})(\widehat{𝐤}𝐅)[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}](\widehat{𝐤}\times 𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐅)\right\}`$ (99)
$`=`$ $`{\displaystyle \frac{i}{v_\mathrm{A}^2}}𝐤𝐯_\mathrm{A}\left\{(\widehat{𝐤}𝐯_\mathrm{A})(\widehat{𝐤}𝐠)[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}](\widehat{𝐤}\times 𝐯_\mathrm{A})(\widehat{𝐤}\times 𝐠)\right\}.`$
In a Thomson scattering medium, $`\mathrm{\Theta }_\rho =\mathrm{\Theta }_{T\mathrm{g}}=0`$, and equations (98)-(99) are identical to the photon bubble dispersion relation derived by Gammie \[1998, his eq. (41)\]. As noted both by Gammie (1998) and BS01, this growth rate only applies for wavelengths longer than the gas pressure scale height $`c_\mathrm{g}^2/g`$. Unless this wavelength scale is optically thin, in which case our diffusion-based analysis breaks down, then the maximal growth rate occurs at wavelengths shorter than the gas pressure scale height, and is given by equations (89) or (93).

There has been some confusion as to how Gammie’s (1998) work on photon bubbles relates to the original analysis performed by Arons (1992), who also assumed a pure Thomson scattering medium with zero gas pressure. Gammie’s (1998) result can in fact be obtained as the short wavelength limit of Arons’ (1992) dispersion relation \[his eq. (33)\], at least within the latter’s assumed vertical magnetic field geometry. Arons’ (1992) numerical growth rates all climb toward higher wavenumber, and must therefore eventually recover the $`k^{1/2}`$ dependence of equation (98), provided the medium remains optically thick at these short wavelengths. Arons (1992) was mainly interested in long wavelengths where radiation diffusion was slow compared to the (radiation) acoustic wave period. Gammie (1998), on the other hand, was working in the rapid diffusion regime at shorter wavelengths.

### 5.4 Numerical Results

Table 2 summarizes the instability thresholds, growth rates, and characteristic wavenumbers for the magnetoacoustic wave instabilities in media with a variety of ratios of magnetic, gas thermal, and radiation energy densities. We have assumed a pure Thomson scattering equilibrium throughout this table, so that the hydrodynamic driving forces that we addressed in sections 3 and 4 vanish, and only the magnetically-based driving remains.

Figure 5 illustrates growth rates of unstable slow waves for small values of the thermal coupling frequency $`\omega _{\mathrm{th}}`$ for a radiation pressure dominated equilibrium, for which the asymptotic high wavenumber growth rate is accurately given by equation (89). Just as in the hydrodynamic case, the short wavelength growth rate declines with increasing thermal coupling. This continues until $`\omega _{\mathrm{th}}`$ exceeds $`g/c_\mathrm{g}`$, beyond which the gas and radiation temperatures are effectively locked and the waves reach unstable asymptotic growth rates given by equation (93), as shown in Figure 6.

Just as in the hydrodynamic case, the unstable growth rate for high thermal coupling only extends over a finite range of wavenumbers. Above a critical cutoff wavenumber, the waves again become damped. After numerical experimentation, we deduced the physical dependencies of these cutoff wavenumbers, and they are listed in Table 2. Figure 7 illustrates the accuracy of these formulas in two specific cases, one with $`v_\mathrm{A}>c_\mathrm{g}`$ and one with $`v_\mathrm{A}<c_\mathrm{g}`$.

## 6 Physics of the MHD Wave Instabilities

The magnetoacoustic wave instabilities are much more robust than their hydrodynamic counterparts, and in particular they survive unscathed even in pure Thomson scattering media. In order to see how this arises physically, we repeat the analysis of section 4 here to re-examine the competing forces acting on a parcel of gas in the wave. There are now two relevant waves that possess density fluctuations: the fast and slow magnetosonic waves. For reference, their polarizations are shown in Figure 8.

Using equation (48) to eliminate the magnetic field perturbation, and equations (69)-(70) to eliminate the total pressure perturbation, the perturbed gas momentum equation (47) may be written in the high wavenumber limit as

$`{\displaystyle \frac{^2\text{¸}}{t^2}}`$ $`=`$ $`(𝐤𝐯_\mathrm{A})^2\text{¸}+(𝐤\text{¸})(𝐤𝐯_\mathrm{A})𝐯_\mathrm{A}+(𝐤\times 𝐯_\mathrm{A})(𝐯_\mathrm{A}\times \text{¸})𝐤𝐤c_\mathrm{g}^2(𝐤\text{¸})`$ (100)
$`+`$ $`i(𝐤\text{¸})𝐤{\displaystyle \frac{\kappa _\mathrm{F}}{c}}\left\{{\displaystyle \frac{4E}{3}}\left({\displaystyle \frac{\omega }{k^2}}\right)+{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho \omega }}{\displaystyle \frac{1}{k^2}}(𝐤𝐅)\left[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}\right]\right\}`$
$`+`$ $`i\left[{\displaystyle \frac{\kappa _\mathrm{F}}{c}}𝐤(\widehat{𝐤}\times 𝐅)(\widehat{𝐤}\times \text{¸})+\text{¸}\times (𝐤\times 𝐠)\right]+𝒪(k^0)|\text{¸}|.`$
Apart from the additional magnetic terms, this equation is of course identical to its hydrodynamic counterpart, equation (74), which we discussed above in section 4. There we found that the last term in square brackets did not affect the growth rate. This continues to be true for the piece of this term that involves the gravitational acceleration $`𝐠`$, as this piece does no work. However the first piece that involves the radiative flux now can play a role because magnetosonic waves are, in general, not longitudinal due to the effects of magnetic tension.

In the short wavelength limit, the fast and slow magnetosonic modes described by equation (100) are polarized in the plane of $`𝐤`$ and $`𝐯_\mathrm{A}`$, and their Lagrangian displacements may be written as

$$\text{¸}=[v_{\mathrm{ph}}^2𝐤(𝐤𝐯_\mathrm{A})𝐯_\mathrm{A}]\psi \text{}\psi ,$$
(101)
where $`\psi `$ is some complex scalar amplitude. Taking the scalar product of equation (100) with the polarization vector , we obtain

$`{\displaystyle \frac{^2}{t^2}}(\text{}\text{¸})`$ $`=`$ $`k^2v_{\mathrm{ph}}^2(\text{}\text{¸})`$ (102)
$`+`$ $`ik^2[v_{\mathrm{ph}}^2(\widehat{𝐤}𝐯_\mathrm{A})^2](𝐤\text{¸}){\displaystyle \frac{\kappa _\mathrm{F}}{c}}\{{\displaystyle \frac{4E}{3}}\left({\displaystyle \frac{\omega }{k^2}}\right)+{\displaystyle \frac{4Ec\omega _\mathrm{a}(\gamma 1)^2}{\kappa _\mathrm{F}\rho \omega }}`$
$``$ $`{\displaystyle \frac{1}{k^2}}(𝐤𝐅)[\mathrm{\Theta }_\rho +(\gamma 1)\mathrm{\Theta }_{T\mathrm{g}}]\}i{\displaystyle \frac{\kappa _\mathrm{F}}{c}}(𝐤𝐯_\mathrm{A})(𝐤\times 𝐯_\mathrm{A})(𝐤\times 𝐅)(𝐤\text{¸}),`$
where, in the last term, we have used the fact that $`(𝐤\text{})(𝐤\times \text{¸})=(𝐤\text{¸})(𝐤\times \text{})`$. The expression for the growth rate, equation (89), now follows immediately from equation (102) and the fact that, to lowest order in $`k^1`$,

$$\frac{𝐤\text{¸}}{\text{}\text{¸}}=\frac{1}{2v_{\mathrm{ph}}^2c_\mathrm{g}^2v_\mathrm{A}^2}.$$
(103)
As illustrated in Figure 8, the polarizations of the fast and slow waves for fixed $`𝐤`$ are orthogonal, and this normalized projection of ¸ onto $`𝐤`$ therefore always has opposite signs for the two modes. Hence the unstable fast and slow modes propagate in different directions if the magnetic driving term dominates over the hydrodynamic terms.

To understand the physics of the MHD driving further, consider again the relevant force density ($`\rho `$ times the second to last term in equation ), which can be written as

$$i\frac{\kappa _\mathrm{F}\rho }{c}𝐤(\widehat{𝐤}\times 𝐅)(\widehat{𝐤}\times \text{¸})=i\frac{\kappa _\mathrm{F}\rho }{c}𝐤\left[(\widehat{𝐤}\text{¸})(\widehat{𝐤}𝐅)(\text{¸}𝐅)\right].$$
(104)
MHD driving results fundamentally from a breakdown in cancellation of the two force densities on the right hand side of this equation, whenever the waves are not longitudinal. The first, $`(\kappa _\mathrm{F}/c)\widehat{𝐤}(\widehat{𝐤}𝐅)\delta \rho `$, is the driving force based on density, not opacity, fluctuations that we illustrated in Figure 4. The second arises from a change in the radiation pressure along a fluid displacement, and is $`i(\kappa _\mathrm{F}\rho /c)𝐤(\text{¸}𝐅)`$. If we ignore the radiation diffusion damping and opacity fluctuations in equation (70), these two force densities simply combine to give the negative gradient of the Lagrangian perturbation in radiation pressure $`\mathrm{\Delta }E/3`$,

$$i\frac{\kappa _\mathrm{F}\rho }{c}𝐤(\widehat{𝐤}\times 𝐅)(\widehat{𝐤}\times \text{¸})=\text{r}\left(\frac{1}{3}\mathrm{\Delta }E\right)=\text{r}\left(\frac{\delta E}{3}+\frac{\text{¸}\text{r}E}{3}\right).$$
(105)
For any purely longitudinal wave (¸ and $`𝐤`$ parallel), the two force densities exactly cancel, and the Lagrangian radiation pressure perturbation then vanishes (again ignoring radiation diffusion damping and opacity fluctuations). For an MHD wave with fixed density perturbation $`\delta \rho `$ and fixed $`𝐤`$, magnetic tension can alter the direction of ¸ off of $`𝐤`$, provided propagation is neither along or perpendicular to the field. Spatial gradients in the flow (which determine the radiative heat flow) are no longer parallel to fluid displacements. The second force density on the right hand side of equation (104) will now be greater or less than the first, depending on the type of MHD wave (slow or fast) and the propagation direction. Note that the two force densities also cancel for vertical propagation.

The physics is very similar in the case where gas and radiation exchange heat rapidly. The perturbed gas momentum equation is then

$`{\displaystyle \frac{^2\text{¸}}{t^2}}`$ $`=`$ $`(𝐤𝐯_\mathrm{A})^2\text{¸}+(𝐤\text{¸})(𝐤𝐯_\mathrm{A})𝐯_\mathrm{A}+(𝐤\times 𝐯_\mathrm{A})(𝐯_\mathrm{A}\times \text{¸})𝐤𝐤c_\mathrm{i}^2(𝐤\text{¸})`$ (106)
$`+`$ $`i(𝐤\text{¸})𝐤{\displaystyle \frac{\kappa _\mathrm{F}}{c}}\left(1+{\displaystyle \frac{3p}{4E}}\right)\left[\left({\displaystyle \frac{4E}{3}}+p\right){\displaystyle \frac{\omega }{k^2}}{\displaystyle \frac{1}{k^2}}(𝐤𝐅)\mathrm{\Theta }_\rho \right]`$
$`+`$ $`i\left[\left(1+{\displaystyle \frac{3p}{4E}}\right){\displaystyle \frac{\kappa _\mathrm{F}}{c}}𝐤(\widehat{𝐤}\times 𝐅)(\widehat{𝐤}\times \text{¸})+\text{¸}\times (𝐤\times 𝐠)\right]+𝒪(k^0)|\text{¸}|.`$
Taking the scalar product of this equation with , we obtain

$`{\displaystyle \frac{^2}{t^2}}(\text{}\text{¸})`$ $`=`$ $`k^2v_{\mathrm{ph}}^2(\text{}\text{¸})`$ (107)
$`+`$ $`ik^2[v_{\mathrm{ph}}^2(\widehat{𝐤}𝐯_\mathrm{A})^2](𝐤\text{¸}){\displaystyle \frac{\kappa _\mathrm{F}}{c}}(1+{\displaystyle \frac{3p}{4E}})[({\displaystyle \frac{4E}{3}}+p){\displaystyle \frac{\omega }{k^2}}`$
$``$ $`{\displaystyle \frac{1}{k^2}}(𝐤𝐅)\mathrm{\Theta }_\rho ]i(1+{\displaystyle \frac{3p}{4E}}){\displaystyle \frac{\kappa _\mathrm{F}}{c}}(𝐤𝐯_\mathrm{A})(𝐤\times 𝐯_\mathrm{A})(𝐤\times 𝐅)(𝐤\text{¸}).`$
With the help of equation (103) with $`c_\mathrm{g}`$ replaced by the isothermal sound speed $`c_\mathrm{i}`$, this now gives us the growth rate of equation (93).

We have extended Hearn’s (1972) hydrodynamic acoustic wave analysis of optically thin media to the magnetosonic wave instabilities, and have found that the extra magnetic piece of the instability vanishes completely in this regime. Physically, this is easy to understand. In the limit of vanishing optical depth, the dynamics of the radiation field can be ignored, and the equilibrium flux $`𝐅`$ merely acts to globally reduce the equilibrium gravitational acceleration $`𝐠`$ to $`𝐠_{\mathrm{eff}}=𝐠\kappa _\mathrm{F}𝐅/c`$. Hence, apart from the hydrodynamic driving terms that are proportional to logarithmic derivatives of $`\kappa _\mathrm{F}`$, the radiation flux enters the gas momentum equation only through a $`\text{¸}\times (𝐤\times 𝐠_{\mathrm{eff}})`$ term \[the last term in equations (100) and (106)\]. This term performs no work. Presumably, there will still be some level of instability in media with optical depths $`\mathrm{}<1`$, but an analysis of this situation is beyond the scope of the present paper.

## 7 Astrophysical Applications

We intend to investigate applications of the theory presented in this paper to particular astrophysical phenomena in future work. For now, we limit ourselves to a brief discussion of how these instabilities might manifest themselves in accretion disks and stars.

### 7.1 Accretion Disks

A number of authors have considered applications of strange modes (Glatzel & Mehren 1996, Mehren-Baehr & Glatzel 1999) and magnetic photon bubbles or magnetoacoustic waves (Gammie 1998; BS01; Begelman 2001, 2002) to radiation pressure dominated accretion disks around black holes.

The gas temperature is probably reasonably well-locked to the radiation temperature in standard accretion disk models. For stellar mass black holes, the most important atomic absorption opacity is bremsstrahlung, for which the Planck mean opacity is given in cgs units by

$$\kappa _\mathrm{P}2.2\times 10^{24}\overline{g}_\mathrm{P}\rho T_\mathrm{g}^{7/2},$$
(108)
where

$$\overline{g}_\mathrm{P}\frac{h}{k_\mathrm{B}T_\mathrm{g}}_0^{\mathrm{}}𝑑\nu \overline{g}_{\mathrm{ff}}e^{h\nu /k_\mathrm{B}T_\mathrm{g}}$$
(109)
is an appropriate frequency-averaged Gaunt factor (Rybicki & Lightman 1979). Equation (108) gives a lower limit to the Planck mean opacity, which should also include bound-free and bound-bound contributions. These contributions are generally more important for the cooler disks thought to exist around supermassive black holes in active galactic nuclei.

By approximating $`\overline{g}_{\mathrm{ff}}`$ to be nearly constant inside the frequency integrals, it is also straightforward to show that for bremsstrahlung,

$$1+\frac{1}{4}\frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}T_\mathrm{r}}\frac{\pi ^2}{6}1.6.$$
(110)
For concreteness, consider a standard Shakura & Sunyaev (1973) disk model in which the anomalous stress is $`\alpha `$ times the total pressure. Using equations (30) and (66), the thermal coupling frequency due to bremsstrahlung is

$`\omega _{\mathrm{th},\mathrm{bremss}}`$ $``$ $`\kappa _\mathrm{P}\rho c\left(1+{\displaystyle \frac{1}{4}}{\displaystyle \frac{\mathrm{ln}\kappa _\mathrm{J}}{\mathrm{ln}T_\mathrm{r}}}\right)\left[{\displaystyle \frac{4(\gamma 1)E}{p}}\right]`$ (111)
$``$ $`3\times 10^4\mathrm{rad}\mathrm{s}^1\overline{g}_\mathrm{P}\eta ^2\alpha ^{7/8}\left({\displaystyle \frac{M}{\mathrm{M}_{}}}\right)^{7/8}\left({\displaystyle \frac{L}{L_{\mathrm{Edd}}}}\right)^2\left({\displaystyle \frac{r}{r_\mathrm{g}}}\right)^{27/16}.`$
Here $`L`$ is the disk luminosity, $`L_{\mathrm{Edd}}`$ is the (Thomson scattering) Eddington luminosity, $`M`$ is the black hole mass, $`\eta L/(\dot{M}c^2)`$ is the disk radiative efficiency, $`r`$ is the radius of the particular point of interest in the disk, and $`r_\mathrm{g}GM/c^2`$ is the gravitational radius of the black hole. We have neglected all relativistic correction factors in the disk model. The particular numerical value of the thermal coupling frequency in equation (111) assumes vertically averaged, disk interior quantities.

The corresponding thermal coupling frequency due to Compton scattering is

$$\omega _{\mathrm{th},\mathrm{Comp}}\kappa _\mathrm{T}\rho c\left(\frac{k_\mathrm{B}T}{m_\mathrm{e}c^2}\right)\left[\frac{4(\gamma 1)E}{p}\right]2\times 10^9\mathrm{rad}\mathrm{s}^1\alpha ^1\left(\frac{M}{\mathrm{M}_{}}\right)^1\left(\frac{r}{r_\mathrm{g}}\right)^{3/2}.$$
(112)
Comparing equations (111) and (112), we see that Compton scattering generally dominates the thermal coupling of gas and radiation in the inner regions of the disk, while bremsstrahlung becomes more important further out, at least in the disk interior for this class of models. For hydrodynamic acoustic waves, the requirement that the gas and radiation temperatures be locked in the wave is that $`\omega _{\mathrm{th}}`$ exceed the sound wave frequency on scales of order the gas scale height $`c_\mathrm{g}^2/g`$, i.e. that $`\omega _{\mathrm{th}}c_\mathrm{g}/g\mathrm{}>1`$. For standard Shakura & Sunyaev (1973) models,

$$\omega _{\mathrm{th},\mathrm{bremss}}\frac{c_\mathrm{g}}{g}4\times 10^4\overline{g}_\mathrm{P}\eta ^3\alpha ^1\left(\frac{L}{L_{\mathrm{Edd}}}\right)^3\left(\frac{r}{r_\mathrm{g}}\right)^{9/2}$$
(113)
and

$$\omega _{\mathrm{th},\mathrm{Comp}}\frac{c_\mathrm{g}}{g}3\times 10^1\eta \alpha ^{9/8}\left(\frac{M}{\mathrm{M}_{}}\right)^{1/8}\left(\frac{L}{L_{\mathrm{Edd}}}\right)^1\left(\frac{r}{r_\mathrm{g}}\right)^{21/16}.$$
(114)
Hence while bremsstrahlung may be unable to ensure tight thermal coupling between the gas and the radiation in the innermost radii, Compton scattering is easily adequate. We stress, however, that these estimates are very uncertain as they do not properly take into account the vertical stratification in the disk. In particular, we have assumed a gravitational acceleration appropriate at the disk scale height but a temperature appropriate to the disk midplane. Because $`\omega _{\mathrm{th},\mathrm{Comp}}c_\mathrm{g}/gT^{3/2}`$, it could be that in regions near the disk photosphere, the gas temperature may not be able to track the radiation temperature.

If the disk is treated hydrodynamically, then unstable acoustic wave growth rates are likely to be much slower than the dynamical frequency in the innermost parts of the disk because Thomson scattering dominates the Rosseland mean opacity. From Table 1, the maximum growth rate for hydrodynamic instabilities is roughly given by

$`\mathrm{\Theta }_\rho {\displaystyle \frac{g}{c_\mathrm{g}}}`$ $``$ $`{\displaystyle \frac{g}{c_\mathrm{r}}}{\displaystyle \frac{c_\mathrm{r}}{c_\mathrm{g}}}{\displaystyle \frac{\kappa _\mathrm{a}}{\kappa _\mathrm{T}}}\mathrm{\Omega }{\displaystyle \frac{c_\mathrm{r}}{c_\mathrm{g}}}{\displaystyle \frac{\kappa _\mathrm{a}}{\kappa _\mathrm{T}}}`$ (115)
$``$ $`3\times 10^5\mathrm{\Omega }\overline{g}_\mathrm{P}\eta \left({\displaystyle \frac{L}{L_{\mathrm{Edd}}}}\right)^1\left({\displaystyle \frac{r}{r_\mathrm{g}}}\right)^{3/2}`$
for a radiation pressure supported disk with $`\kappa _a\kappa _\mathrm{P}`$.

For example, Glatzel & Mehren (1996) found that a class of global strange mode instabilities was present at 25 km radii in accretion disk models around 1 $`\mathrm{M}_{}`$ black holes with accretion rates of $`10^9`$ $`\mathrm{M}_{}`$ yr-1. Our estimate (115) then gives a growth rate $`5\times 10^3\mathrm{\Omega }`$, in remarkable agreement with their calculated value $`3\times 10^3\mathrm{\Omega }`$. Mehren-Baehr & Glatzel (1999) also note that this class of modes has growth rates that decline as the ratio of gas pressure to radiation pressure increases, which is also in agreement with our estimate (115).666Mehren-Baehr & Glatzel (1999) actually state that the growth rates are inversely proportional to $`\beta 3p/E`$, but they do not present direct quantitative evidence for a dependence on $`c_\mathrm{r}^2/c_\mathrm{g}^2`$ rather than $`c_\mathrm{r}/c_\mathrm{g}`$. In addition to these “hot” instabilities, Glatzel & Mehren (1996) and Mehren-Baehr & Glatzel (1999) found a class of unstable modes in colder regions of disk models where gas pressure was substantial. They found that these modes required the presence of a vertical opacity peak due to helium or hydrogen ionization, but that these modes also satisfied the NAR approximation and were therefore likely to be related to strange modes. We speculate that these modes are related to the local, WKB driving we find in the gas-pressure dominated case. Radiation pressure support is not essential for strange modes provided the gas and radiation temperatures are locked, as these authors implicitly assumed in their analysis.

It is almost certainly true, however, that the presence of magnetic fields cannot be neglected in these flows, and these allow in principle much faster growth rates through the slow mode instability (Gammie 1998, BS01). From Table 2, the maximum growth rate in this case is $`\mathrm{\Omega }(c_\mathrm{r}/c_\mathrm{g})\mathrm{min}[1,(v_\mathrm{A}/c_\mathrm{g})]`$, which can be much faster than the orbital frequency in the radiation pressure dominated inner regions. The characteristic length scale of the instability is given by the reciprocal of the turnover wavenumber, which is the gas scale height. Again in the radiation pressure dominated inner regions, this is smaller than the disk scale height by the ratio of the gas to radiation pressure. It is far from clear how this instability might, if at all, manifest itself in the presence of MRI turbulence. The MRI acts on longer time scales (the orbital time), and possibly length scales as well, and this separation of scales between the two instabilities suggests that they might both be simultaneously present. Numerical simulations of vertically stratified, radiating shear flows will be required to address this question.

### 7.2 Stellar Envelopes on the Main Sequence

Except possibly for short wavelength acoustic waves near the stellar photosphere, the gas and radiation temperatures are tightly coupled together in main sequence stellar envelopes. Thermal locking at the turnover wavenumber $`g/c_\mathrm{i}^2`$ requires $`\omega _{\mathrm{th}}\mathrm{}>g/c_\mathrm{i}`$. For gas pressure supported envelopes, for example, this translates to the requirement that

$`\left({\displaystyle \frac{L}{\mathrm{L}_{}}}\right)\left({\displaystyle \frac{M}{\mathrm{M}_{}}}\right)^1\left({\displaystyle \frac{\kappa _F\kappa _\mathrm{a}}{\kappa _T^2}}\right)\left({\displaystyle \frac{\rho }{10^7\mathrm{g}\mathrm{cm}^3}}\right)\left({\displaystyle \frac{2\pi c_\mathrm{i}/g}{3\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}}}\right)^1\mathrm{}>1,`$ (116)
which, given realistic envelope opacities, is usually satisfied.

From Table 1, if the gas and radiation temperatures are locked together, then the vertical radiative flux will destabilize short wavelength hydrodynamic acoustic waves in a radiation pressure dominated stellar envelope if $`F\mathrm{\Theta }_\rho \mathrm{}>Ec_i`$. It is convenient to express this in terms of the flux mean (Rosseland) optical depth $`\tau _\mathrm{F}`$ of the particular envelope layer in question. Radiative diffusion implies that $`FEc/\tau _\mathrm{F}`$, so that short wavelength acoustic waves will be unstable in a radiation pressure dominated stellar envelope for all layers with optical depths satisfying

$$\tau _\mathrm{F}\mathrm{}<\left(\frac{c}{c_\mathrm{i}}\right)\mathrm{\Theta }_\rho .$$
(117)
This can also be written as

$$\mathrm{\Theta }_\rho \mathrm{}>3\times 10^5\left(\frac{c_\mathrm{i}}{10\mathrm{k}\mathrm{m}\mathrm{s}^1}\right)\tau _\mathrm{F}.$$
(118)
In a gas pressure supported stellar envelope, the instability criterion in Table 1 can be written in a similar fashion,

$$\tau _\mathrm{F}\mathrm{}<\left(\frac{c}{c_\mathrm{i}}\right)\left(\frac{E}{p}\right)\mathrm{\Theta }_\rho .$$
(119)
Now, $`p/EL_{\mathrm{Edd}}/L`$ if gas pressure dominates radiation pressure, where $`L_{\mathrm{Edd}}`$ is the Eddington luminosity appropriate for the Rosseland mean opacity $`\kappa _\mathrm{F}`$. Hence the instability criterion (119) may be rewritten as

$$\mathrm{\Theta }_\rho \mathrm{}>3\times 10^5\frac{L_{\mathrm{Edd}}}{L}\left(\frac{c_i}{10\mathrm{km}\mathrm{s}^1}\right)\tau _\mathrm{F}.$$
(120)
This equation strongly suggests that even in gas pressure supported envelopes on the mid- to upper main sequence, acoustic waves may be unstable and strange modes might exist.

We stress that we have only performed a local analysis of the stability of propagating waves in this paper. A demonstration that such waves are amplified by the background radiative flux does not necessarily imply that global compressive modes of the star will be unstable. This is because such modes can in general be viewed as superpositions of WKB waves propagating in different directions and, as we have seen, while one propagation direction may be driven, the other may very well be damped. Our instability criteria are therefore necessary, but not sufficient when it comes to global standing waves of the star. Note that the radiative driving becomes most important at lower optical depths, and it is also important to point out that we have not included the effects of turbulent damping associated with outer convection zones, when they exist.

Most main sequence stellar envelopes probably contain rather weak magnetic fields, and therefore the hydrodynamic driving that we have discussed so far is probably of greatest importance. There are exceptions, however, where MHD effects might play a role. In some chemically peculiar A and B-stars, for example, highly ordered, roughly dipolar magnetic fields of order 1 kG have been measured (e.g. Landstreet 1992). The MHD driving that we discussed in sections 5 and 6 above could produce instabilities in the upper layers of such stars, and magnetic effects will almost certainly alter the hydrodynamic driving. A particular case in point is that of roAp stars, which exhibit high order $`p`$-mode oscillations with low angular momentum quantum number $`\mathrm{}`$, largely confined to the poles of a dipole magnetic field (e.g. Kurtz 1990). Consider the possibility that hydrodynamic driving due to a large opacity variation dominates magnetic driving, but in the surface layers where the magnetic energy density dominates the gas pressure. From equation (90), the phase speed of slow modes in this case is $`|\widehat{𝐤}𝐯_\mathrm{A}/v_\mathrm{A}|c_\mathrm{i}`$, assuming the gas and radiation temperatures are locked together. Equation (93) then implies a growth rate due to hydrodynamic driving which is proportional to $`|(\widehat{𝐤}𝐯_\mathrm{A})(\widehat{𝐤}𝐅)\mathrm{\Theta }_\rho |`$. In other words, slow waves will have maximal driving if they are vertical (i.e. radial) in the regions of the magnetic poles, where the field is also vertical. This might therefore be the basic driving mechanism behind roAp oscillations. We intend to explore this further in future work.

## 8 Conclusions

We have examined the conditions for local driving of a broad class of instabilities of propagating waves in optically thick media. These waves can be driven unstable by a sufficiently strong equilibrium radiative heat flux acting on density fluctuations in the wave. The central mathematical results of this paper are the local dispersion relations: equation (49) for the full MHD case and equation (53) for hydrodynamics. Short wavelength expressions for the hydrodynamic damping and driving of acoustic waves are given in equations (57) and (62). Generalization of these for the full MHD case are given in equations (89) and (93). In addition to the acoustic waves, we have also briefly discussed short wavelength gravity waves in sections 3 and 5.

There are essentially two types of local radiative driving of the acoustic wave instabilities. Hydrodynamic driving, which occurs even in the absence of magnetic fields, requires that the wave possess fluctuations in the flux mean (i.e. Rosseland mean) opacity. A medium with pure Thomson scattering opacity possesses no such local driving. The fastest growth rates occur for waves with vector phase velocities parallel (or possibly antiparallel) to the equilibrium flux $`𝐅`$. Table 1 summarizes the properties of acoustic wave instabilities subject to hydrodynamic driving. It is this driving that is responsible for the existence of global strange mode oscillations in stars, which have normally been found in radiation pressure dominated stellar envelopes. We find that hydrodynamic driving can also produce fast growth rates in propagating acoustic waves in gas pressure dominated envelopes, at least for sufficiently large radiative heat fluxes and provided that the thermal emission and absorption effectively lock the gas and radiation temperatures together.

The presence of a large scale, equilibrium magnetic field expands the ability of an equilibrium radiative flux to drive acoustic waves unstable. In particular, MHD driving requires only that the wave possess density fluctuations, and even a medium with pure Thomson scattering opacity can have unstable magnetosonic waves. If MHD driving dominates over hydrodynamic driving, then the fastest growth rates occur for slow magnetosonic waves, and it is straightforward to show that these growth rates occur for phase velocities that are in the plane of the equilibrium flux $`𝐅`$ and magnetic field $`𝐁`$. In general, the fastest growing waves have phase velocities that are at some angle to both $`𝐅`$ and $`𝐁`$, and the MHD driving in fact vanishes for propagation along $`𝐅`$, along $`𝐁`$, or perpendicular to $`𝐁`$. The characteristics of MHD driving in various limits are summarized in Table 2. MHD driving is responsible for “photon bubble” modes in accreting X-ray pulsars and accretion disks around black holes and neutron stars. It appears to require diffusive radiation transport: MHD driving vanishes in the optically thin limit, in contrast to hydrodynamic driving which survives unscathed (Hearn 1972).

It it important to note that magnetic fields can also alter the behavior of the hydrodynamic driving that exists in the presence of opacity fluctuations, and this is summarized in equations (89) and (93). We suggest that this magnetically modified hydrodynamic driving may play a role in the excitation of observed p-mode oscillations in roAp stars.

The physics of radiative driving may be more generic than just photon-matter interactions. In a coupled two-fluid system, provided one of the fluids dominates the inertia and the other provides rapid diffusive heat transport through which momentum exchange can occur, then unstable driving of density fluctuations of the sort we have discussed here may occur. One example might be diffusive neutrino transport in proto-neutron stars or hyper-Eddington accretion flows. Rapid streaming of cosmic rays in the interstellar medium can also drive acoustic waves unstable (Begelman & Zweibel 1994), though it is unclear whether and how this might be related to the instabilities we have considered here. We hope to explore more detailed applications of radiatively driven instabilities to particular astrophysical phenomena in future work.

We are grateful to J. Arons, P. Arras, L. Bildsten, P. Chang, S. Davis, C. Fryer, C. Gammie, W. Glatzel, P. Goldreich, and N. Turner for very useful discussions. We also thank the anonymous referee for constructive criticism which greatly improved this paper. This work was supported by NSF grant AST-9970827 and UCSB/LANL CARE grant SBB-001A. OB also thanks the generous hospitality of the Caltech Tapir group, where much of this work was carried out.

## Appendix A Appendix: Additional Mathematical Justification for the WKB Damping and Driving Terms

Our WKB treatment presented in sections 3-6 is not rigorous, particularly as the interesting physics (the damping and driving effects on the underlying acoustic and magnetoacoustic waves) lies in a first order correction to the infinitely short wavelength limit. Because we invoked the WKB approximation on our separate perturbation equations before attempting to combine them into one, these corrections are vulnerable to alteration by other equilibrium vertical gradients.

A more rigorous approach would be to combine all nine of our original linear perturbation equations together into one single ordinary differential equation in $`z`$, and then apply WKB techniques on that equation alone. Unfortunately, we have not discovered a way of doing this in general. However, there are two particular cases where this can be done, and we present them here. The results of these cases fully agree with the asymptotic WKB results that we obtained in sections 3-6.

### A.1 Case 1: Vertically Propagating Hydrodynamic Acoustic Waves With Equal Gas and Radiation Temperatures

Assume that the gas and radiation exchange energy quickly enough to guarantee that $`T_\mathrm{r}=T_\mathrm{g}T`$, and neglect the effects of magnetic fields. Consider perturbations that have no variation in the horizontal direction, i.e. waves which propagate (up or down) in the vertical ($`z`$) direction. The linear perturbation equations for these waves are

$$\frac{\delta \rho }{t}+\delta v\frac{d\rho }{dz}+\rho \frac{\delta v}{z}=0,$$
(A1)
$$\rho \frac{\delta v}{t}=\frac{\delta p}{z}g\delta \rho \frac{1}{3}\frac{\delta E}{z},$$
(A2)
$$\frac{}{t}(\delta u+\delta E)+\delta v\left(\frac{du}{dz}+\frac{dE}{dz}\right)+\left(\gamma u+\frac{4}{3}E\right)\frac{\delta v}{z}=\frac{\delta F}{z},$$
(A3)
$$\frac{1}{3}\frac{\delta E}{z}=\frac{\kappa _\mathrm{F}\rho }{c}\delta F\frac{\kappa _\mathrm{F}F}{c}(1+\mathrm{\Theta }_\rho )\delta \rho \frac{\kappa _\mathrm{F}F}{c}\mathrm{\Theta }_T\frac{\delta T}{T},$$
(A4)
together with $`\delta u=\delta p/(\gamma 1)`$, $`\delta p/p=\delta \rho /\rho +\delta T/T`$, and $`\delta E/E=4\delta T/T`$. Note that in this case of vertical propagation, $`\delta 𝐯`$ and $`\delta 𝐅`$ only have nonzero components in the vertical direction.

These equations are still too complicated to combine into a single wave equation, and we must approximate the thermodynamics somewhat. The perturbed energy equation (A3) can be written just in terms of $`\delta T`$, $`\delta F`$, and $`\delta v`$ as

$$\frac{1}{T}\left(\frac{p}{\gamma 1}+4E\right)\left(\frac{\delta T}{t}+\delta v\frac{dT}{dz}\right)+\left(p+\frac{4}{3}E\right)\frac{\delta v}{z}=\frac{\delta F}{z}.$$
(A5)
Now, at very short wavelengths ($`k\mathrm{}`$), we expect rapid radiative diffusion to smooth out temperature fluctuations, so that $`\delta T=𝒪(k^1)\delta v`$ is very small. Hence the dominant terms in this equation at short wavelengths are the last term on the left hand side and the term on the right hand side.777This approximation was also used by Begelman (2001) in his treatment of nonlinear, radiatively driven MHD waves, cf. his equation (7). At the same level of approximation, we may therefore rewrite this equation as

$$\frac{}{z}\left[\left(p+\frac{4}{3}E\right)\delta v+\delta F\right]=0,$$
(A6)
which can be immediately integrated to give

$$\left(p+\frac{4}{3}E\right)\delta v+\delta F=0.$$
(A7)
Similarly, we also neglect all $`\delta T`$ terms in the momentum and diffusion equations (A2) and (A4) that are not gradients to give

$$\rho \frac{\delta v}{t}=\frac{p}{\rho T}\frac{dT}{dz}\delta \rho \frac{p}{\rho }\frac{\delta \rho }{z}\frac{p}{T}\frac{\delta T}{z}g\delta \rho \frac{4E}{3T}\frac{\delta T}{z},$$
(A8)
and

$$\frac{4E}{3T}\frac{\delta T}{z}=\frac{\kappa _\mathrm{F}\rho }{c}\delta F\frac{\kappa _\mathrm{F}F}{c}(1+\mathrm{\Theta }_\rho )\delta \rho ,$$
(A9)
respectively.

With these approximations, equations (A1), (A7), (A8), and (A9) can now be combined into a single equation. Without further loss of generality, we assume all perturbation variables have a time-dependence of the form $`\mathrm{exp}(i\omega t)`$. Then after some algebra, we obtain a single ordinary differential equation in $`\delta v`$:

$`0`$ $`=`$ $`p{\displaystyle \frac{d^2\delta v}{dz^2}}+[{\displaystyle \frac{dp}{dz}}(1+{\displaystyle \frac{3p}{4E}})\mathrm{\Theta }_\rho {\displaystyle \frac{\kappa _\mathrm{F}\rho F}{c}}]{\displaystyle \frac{d\delta v}{dz}}+[\omega ^2\rho +{\displaystyle \frac{i\omega \kappa _\mathrm{F}\rho }{c}}(p+{\displaystyle \frac{4E}{3}})(1+{\displaystyle \frac{3p}{4E}})`$ (A10)
$`{\displaystyle \frac{3p}{4E}}{\displaystyle \frac{\kappa _\mathrm{F}F}{c}}{\displaystyle \frac{d\rho }{dz}}{\displaystyle \frac{\kappa _\mathrm{F}F}{c}}\mathrm{\Theta }_\rho (1+{\displaystyle \frac{3p}{4E}}){\displaystyle \frac{d\rho }{dz}}{\displaystyle \frac{p}{\rho ^2}}\left({\displaystyle \frac{d\rho }{dz}}\right)^2+{\displaystyle \frac{p}{\rho }}{\displaystyle \frac{d^2\rho }{dz^2}}]\delta v.`$
It is convenient to define a new perturbation variable $`\delta \psi `$ by

$$\delta \psi =p^{1/2}\delta v.$$
(A11)
Then equation (A10) becomes

$`0`$ $`=`$ $`p{\displaystyle \frac{d^2\delta \psi }{dz^2}}\left(1+{\displaystyle \frac{3p}{4E}}\right)\mathrm{\Theta }_\rho {\displaystyle \frac{\kappa _\mathrm{F}\rho F}{c}}{\displaystyle \frac{d\delta \psi }{dz}}`$ (A12)
$`+`$ $`[\omega ^2\rho +{\displaystyle \frac{i\omega \kappa _\mathrm{F}\rho }{c}}(p+{\displaystyle \frac{4E}{3}})(1+{\displaystyle \frac{3p}{4E}})+{\displaystyle \frac{1}{2p}}{\displaystyle \frac{dp}{dz}}\mathrm{\Theta }_\rho {\displaystyle \frac{\kappa _\mathrm{F}\rho F}{c}}(1+{\displaystyle \frac{3p}{4E}})`$
$`+{\displaystyle \frac{1}{4p}}\left({\displaystyle \frac{dp}{dz}}\right)^2{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^2p}{dz^2}}{\displaystyle \frac{3p}{4E}}{\displaystyle \frac{\kappa _\mathrm{F}F}{c}}{\displaystyle \frac{d\rho }{dz}}{\displaystyle \frac{\kappa _\mathrm{F}F}{c}}\mathrm{\Theta }_\rho \left(1+{\displaystyle \frac{3p}{4E}}\right){\displaystyle \frac{d\rho }{dz}}`$
$`{\displaystyle \frac{p}{\rho ^2}}\left({\displaystyle \frac{d\rho }{dz}}\right)^2+{\displaystyle \frac{p}{\rho }}{\displaystyle \frac{d^2\rho }{dz^2}}]\delta \psi .`$
If we now employ the WKB ansatz that $`\delta \psi \mathrm{exp}(ikz)`$, this equation gives the dispersion relation

$$\omega ^2=k^2c_\mathrm{i}^2+i\frac{\kappa _\mathrm{F}}{c}\left(1+\frac{3p}{4E}\right)\left[\mathrm{\Theta }_\rho kF\omega \left(p+\frac{4E}{3}\right)\right]+𝒪(k^0),$$
(A13)
or

$$\omega =\pm kc_\mathrm{i}i\frac{\kappa _\mathrm{F}}{2cc_\mathrm{i}}\left(1+\frac{3p}{4E}\right)\left[\left(p+\frac{4E}{3}\right)c_\mathrm{i}\mathrm{\Theta }_\rho F\right]+𝒪(k^1),$$
(A14)
in exact agreement with equation (62) for vertically propagating waves.

An astute reader familiar with WKB techniques may notice that we cheated slightly. The formal WKB treatment first involves eliminating the first order derivative term, thereby forcing the differential equation to resemble the harmonic oscillator equation. We can do this by transforming from $`\delta v`$ to $`\stackrel{~}{\delta \psi }`$, defined by

$$\stackrel{~}{\delta \psi }=p^{1/2}\delta v\mathrm{exp}\left[\mathrm{\Theta }_\rho \frac{\kappa _\mathrm{F}\rho F}{2cp}\left(1+\frac{3p}{4E}\right)𝑑z\right].$$
(A15)
When we employ the WKB ansatz to the differential equation that then results, $`\stackrel{~}{\delta \psi }\mathrm{exp}(i\stackrel{~}{k}z)`$, we obtain

$$\omega =\pm \stackrel{~}{k}c_\mathrm{i}i\frac{\kappa _\mathrm{F}}{2c}\left(1+\frac{3p}{4E}\right)\left(p+\frac{4E}{3}\right)+𝒪(\stackrel{~}{k}^1).$$
(A16)
This appears to be an acoustic wave which is merely damped by radiative diffusion, and we have apparently lost our radiative driving term. However, one must be very careful in the physical interpretation of the wave vector $`\stackrel{~}{k}`$. It is related to the wave vector $`k`$ we employed above by

$$\stackrel{~}{k}=k+i\mathrm{\Theta }_\rho \frac{\kappa _\mathrm{F}\rho F}{2cp}\left(1+\frac{3p}{4E}\right).$$
(A17)
Substituting this expression into equation (A16), we recover the dispersion relation (A14), including the radiative driving term. In order to interpret this ambiguous behavior, note that the velocity perturbation is

$$\delta vp^{1/2}\mathrm{exp}(i\omega t+ikz)p^{1/2}\mathrm{exp}\left[i\omega t+i\stackrel{~}{k}z+\mathrm{\Theta }_\rho \frac{\kappa _\mathrm{F}\rho F}{2cp}\left(1+\frac{3p}{4E}\right)𝑑z\right].$$
(A18)
For concreteness, assume $`\mathrm{\Theta }_\rho >0`$ and consider upward propagating waves. We may choose $`k`$ to be real and large (the interpretation we have adopted throughout this paper), in which case equation (A14) is the correct dispersion relation. If radiative driving dominates damping, then the wave will grow in time. Alternatively, we may choose $`\stackrel{~}{k}`$ to be real and large, in which case equation (A16) is the relevant dispersion relation. The wave will tend to damp with time, but if the same quantitative requirement that the driving term dominates the damping term in equation (A14) still holds, then equation (A18) implies that it will grow exponentially as it propagates spatially. Either way, we have instability: the wave is driven to larger and larger amplitude by the background radiative flux.

Finally, we note that the original partial differential equations (A1), (A7), (A8), and (A9) can be combined into the following energy conservation equation for the waves,

$$\frac{}{t}\left(\frac{1}{2}\rho \delta v^2+\frac{p}{2\rho ^2}\delta \rho ^2\right)+\frac{}{z}\left(\frac{p}{\rho }\delta \rho \delta v\right)=\left(1+\frac{3p}{4E}\right)\left[\frac{\kappa _\mathrm{F}F}{c}\mathrm{\Theta }_\rho \delta \rho \delta v\left(p+\frac{4E}{3}\right)\frac{\kappa _\mathrm{F}\rho }{c}\delta v^2\right].$$
(A19)
This equation leaves no ambiguity in interpretation. The two terms on the right hand side represent radiative driving and damping by diffusion, respectively, and are consistent with the dispersion relation (A14). Note that for $`\mathrm{\Theta }_\rho >0`$, radiative driving requires $`\delta \rho \delta v>0`$, i.e. the density and velocity perturbations must be in phase, which occurs only for upward propagating waves. This is also fully consistent with the physics we described in section 4 and illustrated in figure 4.

### A.2 Case 2: Infinitely Strong Horizontal Field, Equal Gas and Radiation Temperatures

Assume again that the gas and radiation exchange energy quickly so that they have the same temperature, and consider a case where there is a horizontal equilibrium magnetic field that is so strong $`(v_\mathrm{A}/c_\mathrm{g}\mathrm{})`$ that fluid elements are constrained to move horizontally. In a Cartesian coordinate system $`(x,y,z)`$, with $`x`$ along the field and $`z`$ in the upward vertical direction, the linear perturbation equations in this case may be written as

$$\frac{\delta \rho }{t}+\rho \frac{\delta v}{x}=0,$$
(A20)
$$\rho \frac{\delta v}{t}=\frac{\delta p}{x}+\frac{\kappa _\mathrm{F}\rho }{c}\delta F_x,$$
(A21)
$$\frac{}{t}(\delta u+\delta E)+\left(\gamma u+\frac{4}{3}E\right)\frac{\delta v}{x}=\text{r}\delta 𝐅,$$
(A22)
$$\frac{1}{3}\frac{\delta E}{x}=\frac{\kappa _\mathrm{F}\rho }{c}\delta F_x,$$
(A23)
$$\frac{1}{3}\frac{\delta E}{y}=\frac{\kappa _\mathrm{F}\rho }{c}\delta F_y,$$
(A24)
$$\frac{1}{3}\frac{\delta E}{z}=\frac{\kappa _\mathrm{F}\rho }{c}\delta F_z\frac{\kappa _\mathrm{F}F}{c}(1+\mathrm{\Theta }_\rho )\delta \rho \frac{\kappa _\mathrm{F}\rho F}{c}\mathrm{\Theta }_T\frac{\delta T}{T},$$
(A25)
together with $`\delta u=\delta p/(\gamma 1)`$, $`\delta p/p=\delta \rho /\rho +\delta T/T`$, and $`\delta E/E=4\delta T/T`$. The magnetic field has vanished entirely from these equations, as it is considered to be so strong that fluid motions do not bend the field lines. As a result, Alfvén and fast magnetosonic modes have been eliminated from this system of equations. Note that, in contrast to the previous case we considered, we have not made any approximations concerning the rate of radiative diffusion.

Without loss of generality, we assume that all perturbation variables have an $`(x,y,t)`$ dependence proportional to $`\mathrm{exp}[i(k_xx+k_yy\omega t)]`$. Then with some algebraic work, all perturbation variables can be eliminated in favor of $`\delta T`$, giving us a single, second order ordinary differential equation in $`z`$, viz.

$`0`$ $`=`$ $`{\displaystyle \frac{4Ec}{3\kappa _\mathrm{F}\rho }}{\displaystyle \frac{d^2}{dz^2}}\left({\displaystyle \frac{\delta T}{T}}\right)+\left[{\displaystyle \frac{k_x^2F(1+\mathrm{\Theta }_\rho )(p+4E/3)}{\rho \omega ^2pk_x^2}}+F(\mathrm{\Theta }_T4)+{\displaystyle \frac{d}{dz}}\left({\displaystyle \frac{4Ec}{3\kappa _\mathrm{F}\rho }}\right)\right]{\displaystyle \frac{d}{dz}}\left({\displaystyle \frac{\delta T}{T}}\right)`$ (A26)
$`+`$ $`[i\omega ({\displaystyle \frac{p}{\gamma 1}}+4E){\displaystyle \frac{ik_x^2\omega (p+4E/3)^2}{(\rho \omega ^2pk_x^2)}}{\displaystyle \frac{4Ec}{3\kappa _\mathrm{F}\rho }}(k_x^2+k_y^2)`$
$`+k_x^2F(1+\mathrm{\Theta }_\rho ){\displaystyle \frac{d}{dz}}\left({\displaystyle \frac{p+4E/3}{\rho \omega ^2pk_x^2}}\right)]{\displaystyle \frac{\delta T}{T}}.`$
(We have assumed for simplicity that $`\mathrm{\Theta }_\rho `$ and $`\mathrm{\Theta }_T`$ are constants.) Now, we expect that slow magnetosonic modes will have a dispersion relation of the form $`\omega ^2=k_x^2c_\mathrm{i}^2=k_x^2p/\rho `$ plus higher order corrections. Hence adopt the WKB ansatz that $`\delta T/T\mathrm{exp}(ik_zz)`$ with large $`k_x`$, $`k_y`$, and $`k_z`$; $`\omega k`$; and $`\omega ^2k_x^2p/\rho k`$. Then equation (A26) gives the dispersion relation

$$\omega ^2=k_x^2c_\mathrm{i}^2+i\frac{k_x^2}{k^2}\frac{\kappa _\mathrm{F}\rho }{c}\left(1+\frac{3p}{4E}\right)\left[\frac{k_zF(1+\mathrm{\Theta }_\rho )}{\rho }\left(p+\frac{4E}{3}\right)\frac{\omega }{\rho }\right]+𝒪(k^0),$$
(A27)
or

$$\omega =\pm |k_x|c_\mathrm{i}=i\frac{\kappa _\mathrm{F}}{2cc_\mathrm{i}}\left(1+\frac{3p}{4E}\right)\left[\left(\frac{k_x^2}{k^2}\right)\left(p+\frac{4E}{3}\right)c_\mathrm{i}\frac{|k_x|k_z}{k^2}F(1+\mathrm{\Theta }_\rho )\right]+𝒪(k^1),$$
(A28)
which agrees completely with equation (93) in the $`v_\mathrm{A}\mathrm{}`$ limit. This confirms that our crude WKB analysis presented in the main body of the paper also appears to work in the MHD case, and is robust with regard to the radiation physics.

## Appendix B Appendix: Radiative Diffusion and Gas/Radiation Heat Exchange Effects on the Magnetorotational Instability

In our first paper (BS01), in addition to pursuing an initial exploration of radiatively driven instabilities in magnetoacoustic waves, we also examined how radiative diffusion modifies the magnetorotational instability (MRI; Balbus & Hawley 1991). Our analysis in that paper assumed pure Thomson scattering opacity, with no heat exchange between the gas and radiation. Here we wish to briefly consider the behavior of the MRI under the more general thermodynamic assumptions used throughout this paper.

We restrict consideration to axisymmetric perturbations on axisymmetric equilibria rotating with angular velocity $`\mathrm{\Omega }(R)`$, where $`R`$ is the radial distance from the rotation axis. Apart from the radial dependence of $`\mathrm{\Omega }`$, which is the source of free energy for the MRI, we completely neglect all other equilibrium gradients in the flow. In order to avoid time-dependent azimuthal fields in the equilibrium that result from the shear, we assume that the radial component of the magnetic field is zero in the equilibrium.

Using cylindrical polar coordinates $`(R,\varphi ,z)`$, the analysis proceeds very much as in BS01, with a resulting dispersion relation

$$D_{\mathrm{ms}}D_{\mathrm{BH}}+k_z^2v_{\mathrm{A}z}^2v_{\mathrm{A}\varphi }^2\left(k_z^2R\frac{d\mathrm{\Omega }^2}{dR}k^2\stackrel{~}{\omega }^2\right)\frac{k_R^2}{k_z^2}\stackrel{~}{\omega }^2\omega ^4=0.$$
(B1)
Here

$$\stackrel{~}{\omega }^2\omega ^2k_z^2v_{\mathrm{A}z}^2,$$
(B2)
$$D_{\mathrm{ms}}\omega ^2k_z^2\left(C_\mathrm{s}^2+v_{\mathrm{A}\varphi }^2\right).$$
(B3)
$$D_{\mathrm{BH}}\frac{k^2}{k_z^2}\stackrel{~}{\omega }^4\kappa ^2\stackrel{~}{\omega }^24\mathrm{\Omega }^2k_z^2v_{\mathrm{A}z}^2,$$
(B4)
and $`C_\mathrm{s}^2𝒜c_\mathrm{g}^2+c_\mathrm{r}^2+𝒞`$ in the limit of zero stratification. From equations (41)-(44), this is

$`C_\mathrm{s}^2`$ $`=`$ $`\left\{\left(\omega +{\displaystyle \frac{ick^2}{3\kappa _\mathrm{F}\rho }}\right)\left[\omega +{\displaystyle \frac{4i(\gamma 1)\omega _\mathrm{a}E}{p}}\right]+i\omega _\mathrm{a}\omega \right\}^1`$ (B5)
$`\times `$ $`\{\omega (\omega +{\displaystyle \frac{ick^2}{3\kappa _\mathrm{F}\rho }})c_\mathrm{g}^2+[\omega ^29(\gamma 1)\omega _\mathrm{a}\left({\displaystyle \frac{ck^2}{3\kappa _\mathrm{F}\rho }}\right)]c_\mathrm{r}^2`$
$`+`$ $`i\omega _\mathrm{a}\omega [1+{\displaystyle \frac{4(\gamma 1)E}{p}}]c_{\mathrm{s1}}^2\}.`$
Here $`c_{\mathrm{s1}}`$ is the total sound speed in the fluid, defined by

$$c_{\mathrm{s1}}^2\frac{16(\gamma 1)E^2+60(\gamma 1)Ep+9\gamma p^2}{9[p+4(\gamma 1)E]\rho }=\frac{\mathrm{\Gamma }_1(p+E/3)}{\rho },$$
(B6)
where $`\mathrm{\Gamma }_1`$ is the first generalized adiabatic exponent commonly used for matter and radiation (Chandrasekhar 1967).

Equation (B1) should be compared to equation (39) of BS01, which it closely resembles. Not surprisingly, altering the thermodynamics has produced only one difference: a change in the effective sound speed to $`C_\mathrm{s}`$. The basic conclusions of BS01 therefore continue to hold, with only minor quantitative changes:

If the equilibrium azimuthal field is zero ($`v_{\mathrm{A}\varphi }=0`$), then the $`k_R=0`$ MRI mode (the “channel solution”) is completely unaltered by radiation effects.

If radiative diffusion is unimportant on the scale of the MRI critical wavelength, \[$`c\mathrm{\Omega }/(3\kappa _\mathrm{F}\rho v_{\mathrm{A}z}^2)1`$\], then the MRI will be unaltered provided the azimuthal field energy density is smaller than the total thermal energy density in the fluid. To be more precise, the MRI will be unaltered provided $`v_{\mathrm{A}\varphi }^2c_\mathrm{g}^2+c_\mathrm{r}^2`$ if the gas and radiation are not thermally locked together on an orbital time, $`\omega _\mathrm{a}[1+4(\gamma 1)E/p]<\mathrm{\Omega }`$, or $`v_{\mathrm{A}\varphi }^2c_{\mathrm{s1}}^2`$ otherwise. In either case, the essential nature of the MRI is unchanged provided the magnetic field is subthermal.

If radiative diffusion is rapid on the scale of an MRI critical wavelength, then radiation pressure no longer helps enforce incompressibility in the MRI. In a radiation pressure dominated medium where the magnetic field is largely in the azimuthal direction, even subthermal (with respect to the radiation) fields can then reduce the growth rate of axisymmetric MRI modes. Unless $`v_{\mathrm{A}\varphi }^2c_\mathrm{g}^2`$ if gas and radiation are not thermally locked, or $`v_{\mathrm{A}\varphi }^2c_\mathrm{i}^2`$ otherwise, then the growth rate of the MRI will be reduced.