text:
# Local Dynamical Instabilities in Magnetized, Radiation Pressure Supported Accretion Disks

## 1 Introduction

The physical state of the radiation pressure dominated, innermost regions of accretion disks around black holes has been uncertain ever since the early days of accretion disk theory. Standard alpha disk models in which the viscous stress is assumed to scale with the radiation pressure are subject to thermal and viscous instabilities (Lightman & Eardley 1974, Shakura & Sunyaev 1976). These instabilities are sensitive to the assumed prescription for the anomalous viscosity (e.g. Piran 1978), so it remains unclear whether or not they are actually present in real flows. If they are, the accretion flow may adopt a radically different state from that usually envisaged in thin accretion disk theory, such as the multiphase equilibrium recently proposed by Krolik (1998).

In addition to these secular instabilities, dynamical instabilities also exist. First and foremost, a differentially rotating flow with a negative angular velocity gradient and an initially weak magnetic field is unstable to the magnetorotational instability (MRI, Balbus & Hawley 1991). The turbulence resulting from this instability is currently the most plausible mechanism known for generating the anomalous viscosity required in accretion disk models. Dissipation of magnetohydrodynamic waves excited by this turbulence by photon diffusion and photon viscosity has recently been examined by Agol & Krolik (1998). Photon diffusion might also affect the linear development of the instability itself in a laminar, radiation pressure dominated flow, but this issue has not been examined previously.

Gammie (1998) has suggested that the overstable photon bubble modes discussed by Arons (1992) in the context of X-ray pulsars also exist in radiation dominated accretion flows in general. However, his instability analysis was limited to studying a static equilibrium where the effects of differential rotation were entirely neglected. Pietrini & Krolik (2000) have recently investigated convective instabilities in unmagnetized, radiation pressure dominated, differentially rotating flows. They assumed a constant vertical density profile, which necessarily leads to an unstable entropy gradient. Whether or not such unstable gradients exist in reality depends on the vertical dissipation profile, which in turn depends on the structure of the MHD turbulence.

These dynamical instabilities may all play a role at some level in the dynamics and thermodynamics of the radiation pressure dominated portion of accretion disks. It is the purpose of this paper to provide a unified description of all three by deriving a general linear dispersion relation which incorporates them all. Such an analysis will hopefully prove to be a useful guide to numerical simulations which explore the nonlinear development of these instabilities and their effects on the resulting turbulent state of the inner accretion disk.

This paper is organized as follows. In section 2 we discuss our basic equations and assumptions. Then in section 3 we focus exclusively on the MRI by deriving its dispersion relation in the absence of vertical stratification so that the other instabilities are suppressed. We are therefore able to study just the effects of photon diffusion on what is probably the most important of these instabilities. In section 4 we then present the full dispersion relation for a vertically stratified, differentially rotating medium, and discuss its solutions. We discuss the relevance of these solutions to astrophysical accretion disks in section 5 and then summarize our conclusions in section 6.

## 2 Equations and Assumptions

In Eulerian coordinates in an inertial frame, the equations of radiation MHD which we use in this paper are

$$\frac{\rho }{t}+\text{r}(\rho 𝐯)=0,$$
(1)
$$\rho \left(\frac{𝐯}{t}+𝐯\text{r}𝐯\right)=\text{r}p\rho \text{r}\mathrm{\Phi }+\frac{1}{4\pi }(\text{r}\times 𝐁)\times 𝐁+\frac{\kappa _{\mathrm{es}}\rho }{c}𝐅,$$
(2)
$$\frac{p}{t}+𝐯\text{r}p=\frac{\gamma p}{\rho }\left(\frac{\rho }{t}+𝐯\text{r}\rho \right),$$
(3)
$$\frac{E}{t}+𝐯\text{r}E+\frac{4}{3}E\text{r}𝐯=\text{r}𝐅,$$
(4)
$$𝐅=\frac{c}{3\kappa _{\mathrm{es}}\rho }\text{r}E,$$
(5)
and

$$\frac{𝐁}{t}=\text{r}\times (𝐯\times 𝐁).$$
(6)
Here $`\rho `$ is the fluid mass density, $`𝐯`$ is the velocity, $`p`$ is the gas pressure, $`𝐁`$ is the magnetic field, $`\gamma `$ is the gas adiabatic index, $`E`$ is the radiation energy density, $`𝐅`$ is the radiative flux, and $`c`$ is the speed of light. We have followed the notation of Stone, Mihalas, & Norman (1992) with regard to variables associated with the radiation field. In order to facilitate comparison with the work of Gammie (1998), we note that his $`p_g`$ is our $`p`$, his $`J`$ is $`cE/4\pi `$, and his $`𝐇`$ is $`𝐅/4\pi `$.

We neglect the self-gravity of the fluid, and assume that the gravitational potential $`\mathrm{\Phi }`$ originates from the central mass. We completely neglect relativistic effects. We also assume a pure electron scattering opacity $`\kappa _{\mathrm{es}}`$, neglecting absorption opacity entirely. This should be a good first approximation to the inner parts of accretion disks around black holes which are electron scattering dominated. Note, however, that this has the consequence that the gas and radiation are coupled together purely by momentum exchange. They do not exchange heat, which is why the gas energy equation (3) simply reduces to the condition for adiabatic flow. This approximation still recovers the basic MRI and photon bubble instabilities, and so we believe captures the physics that we wish to explore here.

We have closed the radiation moment equations by assuming that the radiation field is close to isotropic, so that the stress tensor is diagonal with elements given by one third the radiation energy density. We are therefore neglecting the effects of photon viscosity considered for example by Agol & Krolik (1998). The radiation momentum equation (5) deserves special comment. We have assumed a simple diffusion form, and have in particular neglected time derivative and velocity-dependent terms that, e.g. Stone et al. (1992) and Gammie (1998) have retained. The full equation is

$$\frac{1}{c^2}\left(\frac{𝐅}{t}+𝐯\text{r}𝐅+𝐅\text{r}𝐯\right)=\frac{1}{3}\text{r}E\frac{\kappa _{\mathrm{es}}\rho }{c}𝐅.$$
(7)
We are interested in perturbations with wavelengths $`\lambda `$ smaller than the vertical disk scale height $`H`$, and angular frequencies $`\omega `$ or growth rates which are of order the local angular velocity of the disk $`\mathrm{\Omega }`$ (MRI, Balbus & Hawley 1991) or possibly a factor $`(2\pi H/\lambda )^{1/2}`$ larger (photon bubbles, Gammie 1998). All the terms on the left hand side of equation (7) will be negligible provided

$$\frac{\omega }{\mathrm{\Omega }}\tau _H\left(\frac{c}{c_\mathrm{s}}\right)\mathrm{and}\frac{v}{c}\frac{\tau _H}{2\pi }\left(\frac{\lambda }{H}\right),$$
(8)
where $`\tau _H\kappa _{\mathrm{es}}\rho H`$ is the scattering optical depth across a scale height, and $`c_s`$ is the total (gas plus radiation pressure) sound speed.

Well outside the innermost stable circular orbit, the scattering depth in the inner, radiation pressure dominated zone of a standard $`\alpha `$-viscosity disk (Shakura & Sunyaev 1973)111It is almost certainly the turbulence generated by the MRI which is responsible for “viscous” stresses in accretion disks, and therefore a value of $`\alpha `$. We are nevertheless using $`\alpha `$ to give us some handle on the ambient conditions in a fictitious, laminar accretion disk prior to the development of any instability. The reader should therefore view these scalings with appropriate caution. is given by

$$\tau _H1\times \alpha ^1\eta \left(\frac{L}{\mathrm{L}_{\mathrm{Edd}}}\right)^1\left(\frac{R}{r_\mathrm{g}}\right)^{3/2}.$$
(9)
Here $`\eta 0.1`$ is the radiative efficiency, $`L/\mathrm{L}_{\mathrm{Edd}}`$ is the luminosity scaled with Eddington, and $`r_\mathrm{g}`$ is the gravitational radius of the central object. The scattering depth is therefore always significantly greater than unity except perhaps in the innermost parts of disks accreting near the Eddington limit. Radiative equilibrium, hydrostatic equilibrium, and the assumption that the viscous stress is $`\alpha `$ times the total pressure leads to a ratio of the speed of light to the sound speed $`c_\mathrm{s}(4E/9\rho )^{1/2}`$ given by

$$\frac{c}{c_\mathrm{s}}\alpha \tau _H.$$
(10)
Hence conditions (8) are almost always satisfied in the interior of accretion disks, provided we do not consider wavelengths that are extremely small. The only exception is for near-Eddington disks in the innermost few gravitational radii, where our Newtonian treatment of the magnetohydrodynamics breaks down anyway. We therefore choose to neglect the left hand side of equation (7), which therefore reduces to the diffusion equation (5). The fact that Gammie (1998) found flux limited electromagnetic waves as solutions of his dispersion relation is a reflection of the fact that he retained these terms. In most cases, these terms have negligible effect on the photon bubble modes, which are at much lower frequency.222Indeed, we have numerically recomputed all the results in this paper with the terms on the left hand side of equation (7) included and found no significant differences for realistic physical parameters. Note that equation (5) allows us to write the radiation pressure very simply in terms of a gradient in the radiation energy density $`E`$, so the momentum equation (2) becomes

$$\rho \left(\frac{𝐯}{t}+𝐯\text{r}𝐯\right)=\text{r}\left(p+\frac{1}{3}E\right)\rho \text{r}\mathrm{\Phi }+\frac{1}{4\pi }(\text{r}\times 𝐁)\times 𝐁.$$
(11)
Eliminating $`𝐅`$ using equation (5), we may write the radiation energy equation (4) as

$$\frac{E}{t}+𝐯\text{r}E+\frac{4}{3}E\text{r}𝐯=\text{r}\left(\frac{c}{3\kappa _{\mathrm{es}}\rho }\text{r}E\right).$$
(12)
In contrast to the radiation momentum equation (7), the left hand side of this equation can be neglected only if

$$\frac{\omega }{\mathrm{\Omega }}\frac{4\pi ^2}{3\tau _H}\left(\frac{c}{c_\mathrm{s}}\right)\frac{4\pi ^2}{3}\alpha \left(\frac{H}{\lambda }\right)^2\mathrm{and}\frac{v}{c}\frac{2\pi }{3\tau _H}\left(\frac{H}{\lambda }\right).$$
(13)
Because $`\alpha `$ may be quite small and $`\tau _H`$ is not very small, these conditions may not be satisfied, so we must retain the left hand side of the radiation energy equation (4). In fact, it is neglecting the right hand side of this equation (photon diffusion) which returns us to the standard ideal MHD equations with gas and radiation pressure. When the conditions (13) start to become satisfied, photon diffusion starts to modify standard MHD. It is this effect which gives rise to photon bubbles, and which can modify the MRI. We therefore choose to retain the left hand side of equation (4), although the right hand side can be dominant.

### 2.1 Equilibrium

We consider a narrow annulus of the disk which, in equilibrium, is stationary and axisymmetric with a purely azimuthal flow velocity $`𝐯=R\mathrm{\Omega }\widehat{\varphi }`$. (We assume cylindrical polar coordinates $`\{R,\varphi ,z\}`$ throughout this paper.) Following the usual analysis of the MRI as well as previous analyses of photon bubbles, we assume the existence of a locally uniform magnetic field. Because we are considering a differentially rotating system, we assume this field has no radial component, which would otherwise be sheared into a time-dependent azimuthal field. Hence $`𝐁=B_\varphi \widehat{\text{Œ}}+B_z\widehat{𝐳}`$. With the exception of the external gravitational potential $`\mathrm{\Phi }`$, we assume that all radial gradients are much smaller than vertical gradients of the same quantities, i.e. that the disk is geometrically thin. The radiative diffusion equation (5) then implies that the equilibrium radiative flux is mainly in the vertical direction, and thus the radial component may be neglected. We will allow for the possibility of radial fluid gradients creating small departures from Keplerian flow by retaining a general form for the angular velocity profile $`\mathrm{\Omega }(R)`$ in our dispersion relation.

With these assumptions, the only nontrivial equilibrium equations are those expressing vertical hydrostatic balance,

$$\frac{1}{\rho }\frac{p}{z}+\frac{\kappa _{\mathrm{es}}}{c}F_zg=0,$$
(14)
radiative equilibrium,

$$\frac{F_z}{z}=0,$$
(15)
and radiation diffusion,

$$F_z=\frac{c}{3\kappa _{\mathrm{es}}\rho }\frac{E}{z}.$$
(16)
Here $`g\mathrm{\Phi }/z`$ is the magnitude of the local vertical gravitational acceleration produced by the tidal field of the central mass.

This equilibrium is characterized by a number of parameters which are important in describing wave modes. First are the gas and radiation sound speeds, defined as

$$c_\mathrm{g}=\left(\frac{\gamma p}{\rho }\right)^{1/2}\mathrm{and}c_\mathrm{r}=\left(\frac{4E}{9\rho }\right)^{1/2},$$
(17)
respectively. Then there are the Alfvén speeds associated with each component of the magnetic field,

$$v_{\mathrm{A}\varphi }=\frac{B_\varphi }{(4\pi \rho )^{1/2}}\mathrm{and}v_{\mathrm{A}z}=\frac{B_z}{(4\pi \rho )^{1/2}}.$$
(18)
Internal hydrodynamic modes are driven by equilibrium gradients with associated oscillation frequencies. The square of the epicyclic frequency is given by

$$\kappa ^2=\frac{1}{R^3}\frac{d}{dR}(R^4\mathrm{\Omega }^2).$$
(19)
It is convenient to define the square of the Brunt-Väisälä frequency separately in the gas and radiation as

$$N_\mathrm{g}^2=g\left[\frac{1}{\rho c_\mathrm{g}^2}\left(\frac{p}{z}\right)\frac{\mathrm{ln}\rho }{z}\right]$$
(20)
and

$$N_\mathrm{r}^2=g\left[\frac{1}{3\rho c_\mathrm{r}^2}\left(\frac{E}{z}\right)\frac{\mathrm{ln}\rho }{z}\right],$$
(21)
respectively. Then the square of the total Brunt-Väisälä frequency is

$$N^2=g\left[\frac{1}{\rho (c_\mathrm{g}^2+c_\mathrm{r}^2)}\frac{}{z}\left(p+\frac{1}{3}E\right)\frac{\mathrm{ln}\rho }{z}\right]=\frac{c_\mathrm{g}^2N_\mathrm{g}^2+c_\mathrm{r}^2N_\mathrm{r}^2}{c_\mathrm{g}^2+c_\mathrm{r}^2}.$$
(22)
Finally, there is a characteristic length scale associated with radiative diffusion, given by how far photons diffuse on an orbital time. This can be used to define a diffusion wavenumber as

$$k_{\mathrm{diff}}\left(\frac{3\kappa \kappa _{\mathrm{es}}\rho }{c}\right)^{1/2}.$$
(23)
### 2.2 Perturbation Equations

Linearizing our equations of motion about the equilibrium, and assuming a WKB space-time dependence $`\mathrm{exp}[i(k_R𝑑R+k_z𝑑z\omega t)]`$, we arrive at the following perturbation equations:

$$i\omega \delta \rho +i\rho 𝐤\delta 𝐯+\delta v_z\frac{\rho }{z}=0,$$
(24)
$$i\omega \rho \delta v_R2\mathrm{\Omega }\rho \delta v_\varphi =ik_R\left(\delta p+\frac{1}{3}\delta E\right)i\frac{k_R}{4\pi }[B_\varphi \delta B_\varphi +B_z\delta B_z]+i\frac{k_z}{4\pi }B_z\delta B_R,$$
(25)
$$i\omega \rho \delta v_\varphi +\frac{\kappa ^2}{2\mathrm{\Omega }}\rho \delta v_R=i\frac{k_z}{4\pi }B_z\delta B_\varphi ,$$
(26)
$$i\omega \rho \delta v_z=ik_z\left(\delta p+\frac{1}{3}\delta E\right)g\delta \rho i\frac{k_z}{4\pi }B_\varphi \delta B_\varphi ,$$
(27)
$$i\omega \delta p+\delta v_z\frac{p}{z}=i\gamma p𝐤\delta 𝐯,$$
(28)
$$i\omega \delta E+\delta v_z\frac{E}{z}+i\frac{4}{3}E𝐤\delta 𝐯=i𝐤\delta 𝐅,$$
(29)
$$\delta F_R=\frac{ic}{3\kappa _{es}\rho }k_R\delta E,$$
(30)
$$\delta F_\varphi =0,$$
(31)
$$\delta F_z=F_z\frac{\delta \rho }{\rho }\frac{ic}{3\kappa _{es}\rho }k_z\delta E,$$
(32)
$$i\omega \delta B_R=ik_zB_z\delta v_R,$$
(33)
$$i\omega \delta B_\varphi R\frac{d\mathrm{\Omega }}{dR}\delta B_R=iB_\varphi 𝐤\delta 𝐯+ik_zB_z\delta v_\varphi ,$$
(34)
and

$$i\omega \delta B_z=ik_RB_z\delta v_R.$$
(35)
The radial and vertical induction equations together guarantee that the divergence of the magnetic field perturbation vanish. Hence one or the other can be replaced with

$$ik_R\delta B_R+ik_z\delta B_z=0.$$
(36)
Equations (24) and (28)-(32) can be solved to yield a simple expression for the perturbed total pressure, viz.

$$\delta p+\frac{1}{3}\delta E=\left(c_\mathrm{g}^2+\frac{\omega }{\omega _{\mathrm{diff}}}c_\mathrm{r}^2\right)\delta \rho \frac{\rho }{g}\left(c_\mathrm{g}^2N_\mathrm{g}^2+\frac{\omega }{\omega _{\mathrm{diff}}}c_\mathrm{r}^2N_\mathrm{r}^2\right)\frac{i\delta v_z}{\omega }\frac{k_zF_z}{3\omega _{\mathrm{diff}}}\frac{\delta \rho }{\rho },$$
(37)
where

$$\omega _{\mathrm{diff}}\omega +i\frac{k^2}{k_{\mathrm{diff}}^2}\kappa .$$
(38)
Each term on the right hand side of equation (37) has a simple interpretation. The first is the acoustic pressure perturbation in response to rarefaction or compression. The second is the buoyancy force (note that $`i\delta v_z/\omega `$ is simply the vertical Lagrangian displacement of a fluid element). The last term is the radiation pressure perturbation arising from the enhanced (diminished) optical depth in compression (rarefaction) regions. The middle term can give rise to convection, while the last term is responsible for photon bubble modes.

The twelve equations (24)-(34) and (36) relate twelve perturbation variables. Eight of these equations have an explicit $`\omega `$ dependence, and so we therefore expect an eighth order dispersion relation. However, as we show in Appendix A, one of these modes is in fact a zero frequency solution which merely represents a perturbation to the equilibrium. The interesting physics therefore lies in a seventh order dispersion relation. In section 4, we will present this full dispersion relation, but we will first discuss in section 3 the simpler case of the MRI in a non-stratified equilibrium.

## 3 The Magnetorotational Instability in the Absence of Stratification

In the fictitious laminar flow we are considering, it is presumably the MRI which generates the turbulence whose dissipation, combined with cooling, sets up the thermal profile which then might drive photon bubble and convective modes. We therefore begin by considering the MRI in isolation, but in a flow which is still hot enough that radiation pressure dominates gas pressure. In the absence of equilibrium stratification, the perturbation equations can be solved to yield the following 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.$$
(39)
Here

$$\stackrel{~}{\omega }^2\omega ^2k_z^2v_{\mathrm{A}z}^2,$$
(40)
$$D_{\mathrm{ms}}\omega ^2k_z^2\left(c_\mathrm{g}^2+\frac{\omega }{\omega _{\mathrm{diff}}}c_\mathrm{r}^2+v_{\mathrm{A}\varphi }^2\right)$$
(41)
are magnetosonic dispersion terms, and

$$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$$
(42)
is the MRI dispersion relation (if set to zero) in a non-stratified, incompressible medium (Balbus & Hawley 1991).

The dispersion relation (39) describes how the MRI is coupled to magnetosonic modes when allowance is made for compressibility and radiation. Without radiation, it reduces to equation (64) of Blaes & Balbus (1994) for $`k_R=0`$, which is all they considered. The conclusion of those authors remains valid here: provided the equilibrium azimuthal magnetic field is sufficiently subthermal, with respect to either the gas or radiation energy densities, then the MRI remains essentially unchanged. Mathematically, this is because the Alfvén speed is then much smaller than the relevant sound speed. For the MRI, $`|\omega |\mathrm{\Omega }kv_{\mathrm{A}z}`$, and the dispersion relation reduces to

$$\left(c_\mathrm{g}^2+\frac{\omega }{\omega _{\mathrm{diff}}}c_\mathrm{r}^2\right)D_{\mathrm{BH}}=0,$$
(43)
so that the MRI mode frequencies remain unchanged from their incompressible form. Physically, in subthermal magnetic fields, the relevant wave speed for the MRI (the Alfvén speed) is much slower than the sound speed, so the associated fluid motions are very nearly incompressible.

The question of how subthermal the azimuthal field must be in order not to affect the MRI is determined by radiative diffusion. From the form of the dispersion relation (39), it is clear that the MRI will be unaffected provided

$$v_{\mathrm{A}\varphi }^2c_\mathrm{g}^2+c_\mathrm{r}^2\left(1+\frac{\mathrm{\Omega }^4}{v_{\mathrm{A}z}^4k_{\mathrm{diff}}^4}\right)^{1/2}.$$
(44)
A sufficient condition for this to be true is that $`v_{\mathrm{A}\varphi }c_\mathrm{g}`$. However, in a radiation pressure dominated region of an accretion disk, it is reasonable to consider azimuthal fields for which $`c_\mathrm{g}<v_{\mathrm{A}\varphi }<c_\mathrm{r}`$. In this case the MRI can be modified if the characteristic wavenumber for instability $`\mathrm{\Omega }/v_{\mathrm{A}z}`$ exceeds the diffusion wavenumber $`k_{\mathrm{diff}}`$, i.e. photons have time to diffuse across a wavelength in an orbital time. This reduces the effective sound speed, making the fluid more compressible. The criterion that the MRI be unaffected is more stringent: $`v_{\mathrm{A}\varphi }c_\mathrm{r}(k_{\mathrm{diff}}v_{\mathrm{A}z}/\mathrm{\Omega })`$.

Compressibility acts to reduce the growth rate of the MRI in the presence of azimuthal fields, but it does not remove the instability (Blaes & Balbus 1994).333We do not have a rigorous proof that the actual stability criterion remains unchanged in our case. However, we suspect that this is true based on our numerical calculations illustrated in Fig. 1. In addition, the dispersion relation (39) admits a zero frequency root for $`\omega `$ when $`k^2v_{\mathrm{A}z}^2=R(d\mathrm{\Omega }^2/dR)`$ (see section 4), which is the condition for the marginally stable wavenumber of the MRI in incompressible MHD in the absence of stratification. This is illustrated in Fig. 1, which depicts the maximum growth rate as a function of azimuthal field strength and diffusion wavenumber. We have restricted consideration to vertical ($`k_R=0`$) modes, as these are expected to have the fastest growth rates, and we consider an equilibrium with negligible gas pressure ($`c_\mathrm{g}=0`$). At high values of $`k_{\mathrm{diff}}v_{\mathrm{A}z}/\mathrm{\Omega }`$, photon diffusion is negligible, and we recover the MRI behavior for an ideal compressible fluid. Subthermal azimuthal fields then produce no effect on the growth rate of the instability, but when the field starts to become thermal in strength, there is a modest reduction in the growth rate. For $`k_{\mathrm{diff}}v_{\mathrm{A}z}/\mathrm{\Omega }\mathrm{}<1`$, photon diffusion is important. In this regime, azimuthal fields can reduce the growth rate quite substantially if $`v_{\mathrm{A}\varphi }\mathrm{}>c_\mathrm{r}(k_{\mathrm{diff}}v_{\mathrm{A}z}/\mathrm{\Omega })`$. For the parameters chosen in the figure, this corresponds to $`v_{\mathrm{A}\varphi }/v_{\mathrm{A}z}\mathrm{}>1`$ for $`k_{\mathrm{diff}}v_{\mathrm{A}z}/\mathrm{\Omega }=0.1`$. Higher values of $`k_{\mathrm{diff}}v_{\mathrm{A}z}/\mathrm{\Omega }`$ (less diffusion) require higher azimuthal field strengths to affect the instability.

We can estimate how small the MRI growth rate becomes when diffusion is important by setting $`k_R=c_\mathrm{r}=0`$ in equation (39) and then solving it in the limit $`c_\mathrm{g}0`$. The MRI then satisfies

$$\omega ^2\frac{k^2c_\mathrm{g}^2v_{\mathrm{A}z}^2\left(R\frac{d\mathrm{\Omega }^2}{dR}+k^2v_{\mathrm{A}z}^2\right)}{v_{\mathrm{A}z}^2\left(R\frac{d\mathrm{\Omega }^2}{dR}+k^2v_{\mathrm{A}z}^2\right)+v_{\mathrm{A}\varphi }^2(\kappa ^2+k^2v_{\mathrm{A}z}^2)}.$$
(45)
For $`v_{\mathrm{A}\varphi }v_{\mathrm{A}z}`$, the unstable growth rate is therefore reduced to $`\mathrm{\Omega }(c_\mathrm{g}/v_{\mathrm{A}\varphi })`$.

Note that regardless of photon diffusion, the MRI is virtually unaffected by compressibility provided the field is sufficiently vertical, i.e. $`v_{\mathrm{A}\varphi }v_{\mathrm{A}z}`$ (e.g. the $`v_{\mathrm{A}\varphi }/v_{\mathrm{A}z}=0.1`$ curve shown in Fig. 1). This is because a flow with a purely vertical equilibrium field will have vertical ($`k_R=0`$) MRI modes identical to the incompressible case, as the dispersion relation (39) permits $`D_{\mathrm{BH}}=0`$ as a solution.

To summarize, the MRI is essentially unaffected by photon diffusion, unless the azimuthal component of the field is comparable to or larger than the vertical component, and the vertical component of the field is weak enough that photon diffusion is important on the scale of an unstable wavelength $`v_{\mathrm{A}z}/\mathrm{\Omega }`$. While the instability still persists, even subthermal azimuthal field components can then reduce the growth rate to $`\mathrm{\Omega }(c_\mathrm{g}/v_{\mathrm{A}\varphi })`$ if condition (44) is violated.

### 3.1 The Relation Between the MRI, Alfvén and Slow Waves

In the next section we will examine stratified equilibria, which turn out to be subject to a number of unstable modes which can couple to each other in complicated ways. To help disentangle this complexity, it is helpful to briefly consider the relation between the MRI and classical MHD waves.

Normally, in the absence of radiative diffusion, the MRI can be viewed as a differential rotation driven destabilization of the slow magnetosonic mode in a weakly magnetized medium (Balbus & Hawley 1998). The distinction between slow and Alfvén modes is delicate, however. In the limit where the magnetic field becomes very weak, the slow magnetosonic mode becomes incompressible and has a dispersion relation which is the same as that of Alfvén waves. It is in fact identical in this limit to the “pseudo-Alfvén” mode which exists in incompressible fluids. The prefix “pseudo” is used because, although the mode is degenerate with the Alfvén mode, the eigenfunction is still different. Alfvén waves have velocity perturbations which are perpendicular to the plane of the wave vector and equilibrium magnetic field. Slow, or pseudo-Alfvén, waves have velocity perturbations in the plane of the wave vector and equilibrium magnetic field. For example, if $`B_\varphi =0`$ in our geometry, Alfvén waves would involve only $`\delta v_\varphi `$, whereas slow waves would have nonzero $`\delta v_R`$ and $`\delta v_z`$.

When the Coriolis forces associated with rotation are included, radial and azimuthal motions are coupled together, so the distinction between slow and Alfvén modes is blurred further. It is therefore perhaps accidental that the stable modes to which the MRI converts are normally the slow modes at high wavenumber. Figure 2(a) illustrates this fact in the case of negligible radiative diffusion. However, when radiative diffusion is important, the MRI can convert to the Alfvén wave \[Fig. 2(b)\] or the slow wave \[Fig. 2(c)\].

This would only be of academic interest were it not for the fact that it is the slow mode which continuously transforms at low wavenumbers into the unstable photon bubble mode in static, stratified equilibria (Gammie 1998). The relationship between this instability and the MRI in linear theory therefore depends on whether the Alfvén or slow mode is more closely tied to the MRI. Also, as we will show in the next section, the Alfvén mode is itself unstable in stratified equilibria.

## 4 The Dispersion Relation in the Presence of Vertical Stratification

With considerable algebra, the general perturbation equations from section 2 can be used to derive the following seventh order dispersion relation:

$`D_{\mathrm{ms}}D_{\mathrm{BH}}+k_z^2v_{\mathrm{A}z}^2v_{\mathrm{A}\varphi }^2\left(k_z^2R{\displaystyle \frac{d\mathrm{\Omega }^2}{dR}}k^2\stackrel{~}{\omega }^2\right){\displaystyle \frac{k_R^2}{k_z^2}}\stackrel{~}{\omega }^2\omega ^4`$
$`+\left(D_{\mathrm{BH}}{\displaystyle \frac{k_R^2\omega ^2\stackrel{~}{\omega }^2}{k_z^2}}\right)\left[{\displaystyle \frac{k_zk_R^2F_z}{3\omega _{\mathrm{diff}}\rho }}+\left(gi{\displaystyle \frac{k_z^2F_z}{3\omega _{\mathrm{diff}}\rho }}\right){\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\rho }{z}}\right]`$
$`+i2\mathrm{\Omega }\omega k_zk_Rv_{\mathrm{A}z}v_{\mathrm{A}\varphi }\left[{\displaystyle \frac{k_zk^2F_z}{3\omega _{\mathrm{diff}}\rho }}i{\displaystyle \frac{k_z^2F_z}{3\omega _{\mathrm{diff}}\rho }}{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\rho }{z}}2ik_zg\right]`$
$`+k_R^2\stackrel{~}{\omega }^2\left(c_\mathrm{g}^2N_\mathrm{g}^2+{\displaystyle \frac{\omega }{\omega _{\mathrm{diff}}}}c_\mathrm{r}^2N_\mathrm{r}^2\right)+{\displaystyle \frac{k_R^2k_z\stackrel{~}{\omega }^2F_z\omega ^2}{3\omega _{\mathrm{diff}}\rho }}{\displaystyle \frac{g}{\rho }}{\displaystyle \frac{\rho }{z}}k_R^2v_{\mathrm{A}\varphi }^2\omega ^2=0.`$ (46)
The first three terms are identical to the non-stratified dispersion relation (39) discussed in the last section. All other terms arise from the stratification. Note that in order for the WKB approximation to be consistent, we must have $`|k_z||\mathrm{ln}\rho /z|`$, so some of the terms with explicit dependence on the vertical density gradient should probably be taken to be negligible. In a radiation pressure dominated accretion disk, the WKB approximation also demands that

$$k\frac{g}{c_\mathrm{r}^2}\frac{\kappa _{\mathrm{es}}F_z}{cc_\mathrm{r}^2}.$$
(47)
The dispersion relation (46) reduces to those of previous authors under various limits. First, in the incompressible limit where $`(c_\mathrm{g}^2+\omega c_\mathrm{r}^2/\omega _{\mathrm{diff}})\mathrm{}`$, we recover the MRI dispersion relation in the presence of vertical stratification (Balbus & Hawley 1991), albeit slightly modified by radiative diffusion:

$$\frac{k^2}{k_z^2}\stackrel{~}{\omega }^4\left[\kappa ^2+\frac{k_R^2}{k_z^2}\left(\frac{c_\mathrm{g}^2N_\mathrm{g}^2+\omega c_\mathrm{r}^2N_\mathrm{r}^2/\omega _{\mathrm{diff}}}{c_\mathrm{g}^2+\omega c_\mathrm{r}^2/\omega _{\mathrm{diff}}}\right)\right]\stackrel{~}{\omega }^24\mathrm{\Omega }^2k_z^2v_{\mathrm{A}z}^2=0.$$
(48)
In addition, if we first assume that gas pressure is negligible ($`c_\mathrm{g}^20`$), adopt the ansatz that $`\omega k^{1/2}`$ and then take the short wavelength limit $`k\mathrm{}`$, we obtain two roots given by

$$\omega ^2=ig\left(\frac{k_zk_R^2}{k^2}\right)\left(\frac{B_z}{B}\right)^2.$$
(49)
This is identical to the asymptotic photon bubble instability dispersion relation obtained by Gammie (1998, eq. 41), adapted to our geometry. In spite of the presence of rotation, the photon bubble instability survives unscathed in this limit.

Finally, if we set the equilibrium magnetic field equal to zero, the equilibrium density to be constant, and the gas pressure to be negligible, then after cancelling two resulting zero modes, the dispersion relation reduces to fifth order:

$$3\omega ^5+i\frac{k^2c}{\kappa _{\mathrm{es}}\rho }\omega ^4\left(3\kappa ^2+\frac{4k^2E}{3\rho }\right)\omega ^3i\frac{k^2c\kappa ^2}{\kappa _{\mathrm{es}}\rho }\omega ^2+\left(\frac{4k_z^2\kappa ^2E}{3\rho }3k_R^2g^2\right)\omega +\frac{k_zk_R^2gc\kappa ^2}{\kappa _{\mathrm{es}}\rho }=0.$$
(50)
This equation resembles the dispersion relation presented by Pietrini & Krolik (2000) when the dissipation terms that they included are neglected. Their dispersion relation still has several additional terms compared to equation (50), however. Each of these terms are small and arise from the left hand side of the full radiation momentum equation (7), which we neglected.

Despite the complexity of the general dispersion relation, the MRI turns out to be a robust instability. An analytic indication that this is so can be found by just setting $`\omega =0`$ in equation (46), which leads to the following equation:

$$k_z^2v_{\mathrm{A}z}^2\left(k^2v_{\mathrm{A}z}^2+R\frac{d\mathrm{\Omega }^2}{dR}\right)\left[i\frac{k_zk_R^2F_zk_{\mathrm{diff}}^2}{3\rho k^2\kappa }+\left(g\frac{k_z^2F_zk_{\mathrm{diff}}^2}{3\rho k^2\kappa }\right)\frac{1}{\rho }\frac{\rho }{z}k_z^2c_\mathrm{g}^2\right]k_R^2k_z^2v_{\mathrm{A}z}^2c_\mathrm{g}^2N_\mathrm{g}^2=0.$$
(51)
If the gas Brunt-Väisälä frequency $`N_\mathrm{g}=0`$, this equation shows that there will exist a zero frequency root if $`k<k_{\mathrm{crit}}[R(d\mathrm{\Omega }^2/dR)]^{1/2}/v_{\mathrm{A}z}`$, the usual critical wavenumber for the onset of the MRI. Even if $`N_\mathrm{g}0`$, there will be a zero frequency root for $`k_R=0`$ and $`k^2v_{\mathrm{A}z}^2=R(d\mathrm{\Omega }^2/dR)`$. While not a rigorous proof, these facts suggest that the critical wavenumber for the onset of the MRI is unchanged by the effects of radiation. Hence the instability criterion for the MRI is still that the angular velocity decreases outward in magnitude. These expectations are born out by our numerical calculations. Radiative diffusion does not affect the stability criterion, although it can dramatically affect the growth rates, as discussed in section 3.

### 4.1 Photon Bubble Modes and Other MHD Wave Instabilities

In the high wavenumber limit $`k\mathrm{}`$, the dispersion relation can be factored to give

$$\left(\omega +i\frac{k^2c}{3\kappa _{\mathrm{es}}\rho }\right)(\omega ^2k_z^2v_{\mathrm{A}z}^2)[\omega ^4\omega ^2k^2(c_\mathrm{g}^2+v_{\mathrm{A}z}^2+v_{\mathrm{A}\varphi }^2)+k_z^2v_{\mathrm{A}z}^2k^2c_\mathrm{g}^2]=0.$$
(52)
The first factor is a strongly damped radiation diffusion mode, the second describes Alfvén waves, and the third corresponds to the fast and slow magnetosonic waves. Rotation is irrelevant at these high wavenumbers. We have lost the radiation effects on the fast and slow modes because in this limit the wavelength is much less than the diffusion length scale.

It is worthwhile examining the behavior of the MHD wave modes to higher order in $`k^1`$, because they all turn out to be unstable for sufficiently high vertical photon fluxes. If, as we have been generally assuming, the radiation and magnetic energy densities in the plasma are much larger than the gas thermal energy density, then the fast and slow modes have frequencies given by

$$\omega \pm k(v_{\mathrm{A}z}^2+v_{\mathrm{A}\varphi }^2)^{1/2}i\frac{3\kappa _{\mathrm{es}}\rho }{2c(v_{\mathrm{A}z}^2+v_{\mathrm{A}\varphi }^2)}\left[c_\mathrm{r}^2\left(\frac{k_R^2}{k^2}v_{\mathrm{A}z}^2+v_{\mathrm{A}\varphi }^2\right)\frac{k_zk_R^2F_zv_{\mathrm{A}z}^2}{3k^3\rho (v_{\mathrm{A}z}^2+v_{\mathrm{A}\varphi }^2)^{1/2}}\right]$$
(53)
and

$$\omega \pm k_zc_\mathrm{g}\left(\frac{B_z}{B}\right)i\frac{\kappa _{\mathrm{es}}k_R^2F_z}{2cc_\mathrm{g}k^2}\left(\frac{B_z}{B}\right),$$
(54)
respectively. The lowest order correction to the Alfvén mode is of order $`k^1`$, and this correction is real, implying that any instability or damping is at still higher order. As we shall see, instabilities in the Alfvén waves peak at finite wavenumbers.

The first imaginary term in equation (53) represents damping of the fast modes by radiative diffusion. This term agrees with the optically thin damping rate derived by Agol & Krolik (1998), although the waves need not be optically thin for this rate to apply. The last term in equation (53) shows that a sufficiently large vertical radiative flux can overcome the damping and destabilize the fast waves. The instability criterion can be obtained by comparing the last two terms of equation (53):

$$F_z>3\rho c_\mathrm{r}^2\left(1+\frac{k^2v_{\mathrm{A}\varphi }^2}{k_R^2v_{\mathrm{A}z}^2}\right)\frac{k}{k_z}(v_{\mathrm{A}z}^2+v_{\mathrm{A}\varphi }^2)^{1/2}.$$
(55)
Crudely, if the speed at which photons are diffusing upward in the equilibrium exceeds the Alfvén speed, then fast modes propagating upward at some nonzero angle to the vertical will be unstable. For large photon fluxes, the growth rate of this instability is

$$\mathrm{Im}(\omega )\frac{g}{2v_A}\left(\frac{k_zk_R^2}{k^3}\right)\left(\frac{B_z}{B}\right)^2.$$
(56)
The instability criterion is most easily satisfied, and the largest growth rates are achieved, for vertical magnetic fields. Azimuthal fields tend to be stabilizing.

Equation (54) reveals that slow modes propagating downward at some nonzero angle to the vertical are also unstable at high photon fluxes. These modes have the fastest growth rates, and at lower values of $`k`$ become the photon bubble modes described by Gammie’s (1998) dispersion relation (49). Adopting Gammie’s (1998) photon bubble ansatz of $`\omega k^{1/2}`$ and retaining only the leading order terms as $`k\mathrm{}`$ in the general dispersion relation, we find

$$\omega ^2=\left(k_z^2c_\mathrm{g}^2i\frac{k_zk_R^2\kappa _{\mathrm{es}}F_z}{ck^2}i\omega \frac{3\kappa _{\mathrm{es}}\rho c_\mathrm{r}^2k_z^2}{ck^2}\right)\left(\frac{B_z}{B}\right)^2.$$
(57)
The middle term gives the photon bubble dispersion relation (49) when use is made of the hydrostatic equilibrium condition (14) in the radiation dominated limit. The other terms compete with the photon bubble term at low and high wavenumbers, leading us to expect that the photon bubble dispersion relation (49) will only hold in the range of wavenumbers given by

$$\frac{9\rho ^2c_\mathrm{r}^2\kappa _{\mathrm{es}}^2}{c^2g}\mathrm{}<k\mathrm{}<\frac{g}{c_\mathrm{g}^2}.$$
(58)
A condition equivalent to this was also noted by Gammie (1998). When the last inequality breaks down at high wavenumbers, the slow mode dispersion relation is better described by equation (54).

Figure 3 presents a numerical solution to the full dispersion relation (46) which illustrates these high wavenumber instabilities and compares their growth rates to the analytic expressions above. In addition to the fast and slow modes, the Alfvén modes also appear to go unstable for large photon fluxes. In contrast to the fast and slow modes, the Alfvén wave growth rates tend to zero as $`k\mathrm{}`$. Rotation directly affects the unstable growth rates of the Alfvén modes, and also affects the fast and slow modes at low wavenumbers. We describe this, and the coupling to the MRI, in the next subsections.

### 4.2 Zero Azimuthal Field Case

We now return to low wavenumbers where rotation is important. We first consider the case where the equilibrium field is vertical ($`B_\varphi =0`$). The general dispersion relation (46) then reduces to

$`D_{\mathrm{ms}}D_{\mathrm{BH}}{\displaystyle \frac{k_R^2}{k_z^2}}\stackrel{~}{\omega }^2\omega ^4+\left(D_{\mathrm{BH}}{\displaystyle \frac{k_R^2\omega ^2\stackrel{~}{\omega }^2}{k_z^2}}\right)\left[{\displaystyle \frac{k_zk_R^2F_z}{3\omega _{\mathrm{diff}}\rho }}+\left(gi{\displaystyle \frac{k_z^2F_z}{3\omega _{\mathrm{diff}}\rho }}\right){\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\rho }{z}}\right]`$
$`+k_R^2\stackrel{~}{\omega }^2\left(c_\mathrm{g}^2N_\mathrm{g}^2+{\displaystyle \frac{\omega }{\omega _{\mathrm{diff}}}}c_\mathrm{r}^2N_\mathrm{r}^2\right)+{\displaystyle \frac{k_R^2k_z\stackrel{~}{\omega }^2F_z\omega ^2}{3\omega _{\mathrm{diff}}\rho }}=0.`$ (59)
Note that $`D_{\mathrm{BH}}=0`$ is a solution of this equation when $`k_R=0`$, so we immediately conclude that in the absence of an azimuthal field, the vertical MRI modes are completely unaffected by compressibility, stratification, and radiation. This is not surprising, as the vertical MRI modes involve horizontal motions only, and therefore do not involve any excess forces from the vertical radiative flux.

Equation (59) is complicated, with many parameters. For numerical computations, we restrict consideration to cases where the gas pressure is negligible, so that $`c_\mathrm{g}`$ and $`N_\mathrm{g}`$ are both zero. In addition, we neglect convection and couplings to gravity modes by considering only isentropic equilibria, so that $`N_\mathrm{r}=0`$. Under these assumptions, we may use equation (14) to write the acceleration due to gravity and the equilibrium vertical density gradient in terms of the photon flux, viz.

$$g=\frac{k_{\mathrm{diff}}^2}{3\kappa \rho }F_z\mathrm{and}\frac{1}{\rho }\frac{\rho }{z}=\frac{k_{\mathrm{diff}}^2}{3\kappa \rho c_\mathrm{r}^2}F_z.$$
(60)
Hence the only new parameter is $`F_z`$.

Figure 4 presents numerical calculations of the unstable mode growth rates in the diffusive regime at low wavenumber, for various values of photon flux. The introduction of stratification causes the slow (photon bubble) and Alfvén wave instabilities to appear and have increasing growth rate as $`F_z`$ gets larger. The fast wave instability appears only after a critical threshold in $`F_z`$ is reached \[cf. Fig. 4(b)\]. This is quantitatively consistent with our instability criterion (55), which predicts that the fast wave will be unstable for $`F_z/(3\rho v_{\mathrm{A}z}^3)>125`$ for the parameters chosen in Fig. 4. Fig. 4(c) should be compared with Fig. 3 which has similar radiation parameters. Note that for the value of $`k_{\mathrm{diff}}`$ chosen in Fig. 4, the MRI converts to Alfvén waves at $`k_{\mathrm{crit}}`$ \[cf. Fig. 2(b)\], and indeed the unstable Alfvén waves show up for $`k>k_{\mathrm{crit}}`$. The MRI itself appears to be largely unaffected by stratification, except that the maximum growth rate increases with increasing $`F_z`$. It is unclear what would happen at higher fluxes because we are already close to the WKB limit (47) at these high fluxes. Given that the vertical ($`k_R=0`$) MRI modes are unaffected by flux, they are not always the fastest growing modes in the presence of nonzero flux. Modes with a slight nonzero $`k_R`$ in fact grow faster for the equilibrium parameters shown in Fig. 4(c). The fact that $`F_z`$ can affect the MRI growth rates for $`k_R0`$, albeit modestly, suggests that geometrically thick accretion flows with significant radial flux in the equilibrium could have a more marked effect on the MRI than the geometrically thin configurations considered in this paper.

The photon bubble mode is always the most rapidly growing of the three MHD wave instabilities that exist in stratified media. It appears to survive for $`k<k_{\mathrm{crit}}`$, albeit at a smaller growth rate than expected from extrapolation of the high $`k`$ behavior. Rotation therefore acts to mildly suppress, and does not kill, this instability when $`B_\varphi =0`$.

### 4.3 Nonzero Azimuthal Field

We finally turn to the full dispersion relation (46) with nonzero azimuthal field, but continuing to neglect gas pressure and vertical entropy gradients. Figure 5 presents results for the low wavenumber instability growth rates corresponding the the $`B_\varphi =0`$ case shown in Fig. 4(b). As we stated above, the MRI is a robust instability which even in the presence of a nonzero azimuthal field and radiation, still exists for $`k<k_{\mathrm{crit}}`$. However, as expected from our results in section 3, the nonzero azimuthal field produces a marked reduction in the overall growth rate. The maximum MRI growth rates in Fig. 5 are comparable to the $`k_R=0`$ growth rates shown in Fig. 1 at the chosen diffusion wavenumber $`k_{\mathrm{diff}}v_{\mathrm{A}z}/\mathrm{\Omega }=0.1`$. The vertical radiative flux makes little difference to the MRI compared to the nonstratified case, just as we found in the previous subsection.

The azimuthal field component has a noticeable effect on the other instabilities. First, the fast and Alfvén waves are quickly stabilized. This is consistent with our findings in section 4.1, where we noted that azimuthal fields would stabilize the fast waves. For the parameters chosen in Fig. 5, equation (55) predicts fast wave stability for $`B_\varphi \mathrm{}>0.4B_z`$.

The photon bubble growth rate also decreases with increasing azimuthal field, in quantitative agreement with equation (49). At high values of $`B_\varphi /B_z`$, this mode appears to split into a high wavenumber part (the true photon bubble), and a low wavenumber part \[cf. Fig. 5(b)\]. The latter dominates the MRI growth rate except at the lowest wavenumbers, at least for the parameters chosen in Fig. 5. At high wavenumbers, this mode evolves into the Alfvén mode.

## 5 Astrophysical Relevance

It turns out that diffusion is actually very important in standard $`\alpha `$-viscosity models of accretion disks around black holes and neutron stars (Shakura & Sunyaev 1973). Neglecting relativistic correction factors, one can show that in the radiation pressure dominated inner zone,

$$k_{\mathrm{diff}}H2\alpha ^{1/2}\mathrm{radians}.$$
(61)
Surprisingly, this relation is independent of accretion rate, radius, and central object mass. Combining this with the WKB condition for the MRI, $`k\mathrm{\Omega }/v_{\mathrm{A}z}H^1`$, we find

$$\frac{k_{\mathrm{diff}}v_{\mathrm{A}z}}{\mathrm{\Omega }}\alpha ^{1/2}.$$
(62)
Depending on the assumed strength of the initial magnetic field, realistic accretion disk models ($`\alpha 0.010.1`$) will often be in the regime where photon diffusion can significantly slow the growth of the MRI in the presence of azimuthal field components (cf. Fig. 1). Our analysis shows that in this case the growth rate would be reduced to $`\mathrm{\Omega }(c_\mathrm{g}/v_{\mathrm{A}\varphi })`$. In the radiation pressure dominated inner zone, one can conceive of $`v_{\mathrm{A}\varphi }`$ being somewhat less than $`c_\mathrm{r}`$. Neglecting relativistic correction factors once again,

$$\frac{c_\mathrm{g}}{c_\mathrm{r}}4\times 10^3\alpha ^{1/8}\eta \left(\frac{M}{\mathrm{M}_{}}\right)^{1/8}\left(\frac{L}{\mathrm{L}_{\mathrm{Edd}}}\right)\left(\frac{R}{r_\mathrm{g}}\right)^{21/16},$$
(63)
where $`M`$ is the mass of the central object. The MRI growth rate can therefore become quite sluggish if the azimuthal field energy density exceeds the thermal energy density in the gas by an amount comparable to the radiation/gas thermal energy density ratio. We stress however that even in this regime, the MRI still exists as a dynamical instability, and will still presumably grow eventually to nonlinear amplitude. It might be interesting to see if the resulting turbulence depends on how the initial magnetic field energy density compares to the gas and radiation thermal energy densities, but it is difficult to predict any results from linear theory alone.

Note that all these conclusions refer to a fictitious, unstable initial state which cannot exist in nature. The equations above, in particular that for the diffusion wavenumber (62), depend on $`\alpha `$, which is presumably a consequence of the MRI. It is therefore difficult to determine how the MRI will be affected when it is responsible for the ambient conditions in the disk in the first place.

This fact is even more important when we consider the stratified calculations that we presented in section 4. Multiple instabilities are simultaneously present in this case, and the modes with the fastest growth rates presumably dominate the nonlinear evolution. One would at first conclude that the photon bubble/slow mode at high wavenumbers therefore controls the dynamics. From equation (54), the growth rate is at most

$$\mathrm{Im}(\omega )10^2\mathrm{\Omega }\alpha ^{1/8}\eta ^1\left(\frac{L}{\mathrm{L}_{\mathrm{Edd}}}\right)\left(\frac{M}{\mathrm{M}_{}}\right)^{1/8}\left(\frac{R}{r_\mathrm{g}}\right)^{21/16},$$
(64)
substantially faster than the orbital frequency $`\mathrm{\Omega }`$.

However, it is important to recognize the different sources of free energy which are driving these instabilities. The MRI feeds off the differential rotation in the flow, transporting angular momentum outward which ultimately allows material to accrete and release gravitational potential energy. We therefore view this instability as being fundamental to the dynamics of the accretion disk. Dissipation of the gravitational energy released and radiative transport processes set up vertical gradients which drive the photon bubble and other MHD wave instabilities. Convection will also presumably occur if the resulting entropy gradient has the right sign (Pietrini & Krolik 2000). It is the vertical gradients which provide the source of free energy for these extra instabilities, but these gradients are fundamentally determined by the action of the MRI.

We speculate that on the largest vertical scales, where diffusion is negligible on an orbital time, the MRI will always operate as usual, launching perhaps a turbulent cascade to smaller scales. How this energy is ultimately dissipated is still unclear. If there is sufficient vertical flux, it may not be as simple as radiative damping of compressive MHD waves (Agol & Krolik 1998), as those waves are themselves unstable for some propagation directions. Given the high growth rates at short wavelengths from equation (64) above, it seems likely to us that photon bubbles will play an important role in thermal transport on small scales.

Concrete answers to the question of how radiative diffusion and radiative flux really affect the turbulent angular momentum and heat transport in the disk will require investigation by numerical simulations. One must be careful to choose the fictitious unstable initial state with care when pursuing such simulations, keeping in mind that accretion power, and therefore the MRI, is the driver for all this complex dynamics.

## 6 Conclusions

We may summarize our conclusions as follows. In the absence of an azimuthal field, the vertical ($`k_R=0`$) MRI modes have the standard, incompressible behavior elucidated by Balbus & Hawley (1991). If the azimuthal field component dominates the vertical field component and the Alfvén speed is larger than the effective sound speed \[cf. eq. (44)\], then the growth rate of the MRI is reduced to $`\mathrm{\Omega }(c_\mathrm{g}/v_{\mathrm{A}\varphi })`$, well below the usual rapid growth rate of order the orbital frequency $`\mathrm{\Omega }`$. However, the instability still persists. Large vertical photon fluxes in the equilibrium do not significantly affect the most rapidly growing MRI modes.

At short wavelengths where rotation is unimportant, both the fast and slow MHD waves are destabilized by sufficiently high vertical photon flux in the equilibrium. These instabilities require that the waves propagate at some nonzero angle to the vertical, and affects upward propagating fast waves and downward propagating slow waves. The slow waves appear to always have the larger growth rates. At smaller wavenumbers these are well-described by Gammie’s (1998) photon bubble dispersion relation (49), but at larger wavenumbers they asymptote to a maximum value given by equation (54). At small wavenumbers where rotation is important ($`k<k_{\mathrm{crit}}`$), the photon bubble growth rate is somewhat suppressed by rotation, at least at modest fluxes. Nonzero azimuthal fields stabilize the fast wave instability and reduce the photon bubble growth rate. In addition to the fast and slow wave instabilities, the vertical radiative flux also drives unstable modes at intermediate wavenumber which are related to Alfvén waves at high wavenumber.

The ambient conditions in the radiation pressure dominated inner region of standard alpha disk models are generally in a regime where radiative diffusion can slow the growth of the MRI. In addition, photon bubble modes are expected at short wavelengths with growth rates that can significantly exceed the orbital frequency. It is not clear how these instabilities will evolve into the nonlinear, turbulent regime, but we speculate that the MRI will dominate the largest scales and that photon bubbles (and perhaps unstable fast and Alfvén waves) may affect both the dissipation and thermal transport on smaller scales. A complete understanding of the turbulent angular momentum and heat transport will require numerical simulations.

Finally, we stress that we have concentrated here on geometrically thin accretion disk configurations where the photon flux is primarily in the vertical direction and the rotation curve is near-Keplerian. Geometrically thick or quasi-spherical, radiation pressure supported flows will have radiation pressure gradient components in the radial direction and may therefore exhibit more interesting dynamical coupling between the radiative flux and the MRI.

We acknowledge useful conversations with Eric Agol, Yuri Levin, and Jim Stone. This work was supported by NSF grants AST-9970827 and PHY-9907949, and NASA grant NAG5-7075.

## Appendix A Appendix A: The Zero Frequency Mode of the General Dispersion Relation

Despite the fact that the perturbation equations (24)-(34) and (36) have eight separate $`\omega `$-dependent terms, we only end up with a dispersion relation that is seventh order. The reason is that there is a zero frequency mode, which we discuss in some detail in this appendix.

Setting $`\omega =0`$ in the perturbation equations, the radial induction equation (33) immediately implies that $`\delta v_R`$ must vanish. The continuity equation (24) and azimuthal momentum equation (26) then imply that $`\delta v_z`$ and $`\delta B_\varphi `$ must also vanish. The gas energy equation (28) is then trivially satisfied, which is the mathematical reason why a zero frequency mode is a possible solution of the perturbation equations.

Excluding the trivial azimuthal flux equation (31), the remaining seven equations may be written

$$2\mathrm{\Omega }\rho \delta v_\varphi +ik_R\left(\delta p+\frac{1}{3}\delta E\right)+i\frac{k_R}{4\pi }B_z\delta B_zi\frac{k_z}{4\pi }B_z\delta B_R=0,$$
(A1)
$$ik_z\left(\delta p+\frac{1}{3}\delta E\right)+g\delta \rho =0,$$
(A2)
$$𝐤\delta 𝐅=0,$$
(A3)
$$\delta F_R+\frac{ic}{3\kappa _{es}\rho }k_R\delta E=0,$$
(A4)
$$\delta F_z+F_z\frac{\delta \rho }{\rho }+\frac{ic}{3\kappa _{es}\rho }k_z\delta E=0,$$
(A5)
$$R\frac{d\mathrm{\Omega }}{dR}\delta B_Rik_zB_z\delta v_\varphi =0,$$
(A6)
and

$$k_r\delta B_R+k_z\delta B_z=0.$$
(A7)
These seven equations have eight unknowns, so that any one of the nonzero perturbation amplitudes determine the other seven for this zero frequency mode.

It is perhaps surprising that this mode has such a complicated eigenfunction, and it is therefore worthwhile understanding its physical meaning. Equation (A6) is particularly interesting. This arose from the azimuthal component of the induction equation which describes, for example, how a radial magnetic field component is sheared out by the differential rotation to produce a time-dependent azimuthal field component. In this case, however, this is exactly compensated by an azimuthal velocity perturbation which creates azimuthal field in the opposite direction out of vertical field.

Further insight into the physical nature of this mode may be obtained by taking various limits. In the limit $`k_z0`$, we have $`\delta B_R=\delta \rho =\delta E=\delta 𝐅=0`$, and the only nontrivial condition comes from the radial momentum equation:

$$2\mathrm{\Omega }\rho \delta v_\varphi +ik_R\delta p+i\frac{k_R}{4\pi }B_z\delta B_z=0.$$
(A8)
In other words, any random field of perturbations of $`\delta v_\varphi `$ can be balanced by radial gradients of gas and/or magnetic pressure to produce a slightly different equilibrium.

In the opposite limit $`k_R0`$, we have $`\delta v_\varphi =\delta B_z=\delta B_R=\delta F_z=0`$, and

$$ik_z\left(\delta p+\frac{1}{3}\delta E\right)+g\delta \rho =0,$$
(A9)
$$\delta F_R+\frac{ic}{3\kappa _{es}\rho }k_R\delta E=0,$$
(A10)
$$F_z\frac{\delta \rho }{\rho }+\frac{ic}{3\kappa _{es}\rho }k_z\delta E=0.$$
(A11)
In this limit, a random field of density perturbations is balanced by vertical gradients in gas and radiation pressure, again producing a slightly different equilibrium.

Finally, we note that if we neglect the vertical stratification of the equilibrium, then this zero mode has all perturbations vanishing except for the density. In this case a random field of density perturbations will merely be advected along with the equilibrium flow without producing additional net accelerations.


format: text

text:
# Evidence for novel polaronic Fermi-liquid in doped oxides

## Abstract

The modern multi-polaron theory predicts a polaronic Fermi-liquid state, but such a state has not been unambigously confirmed by experiment so far. Here we report theoretical and experimental studies of the isotope effects on the low temperature kinetic properties of the doped ferromagnetic manganites. Our results provide clear evidence for a novel polaronic Fermi-liquid state in doped oxides.

Over the last decades, the very interesting phenomena such as high-temperature superconductivity and colossal magnetoresistance have been discovered in several doped oxides (e.g., cuprates , bismuthates , and manganites ). These doped oxides are characterized by significant carrier densities ($`10^{21}`$ cm-3) and low mobility of the order or even less than the Mott-Ioffe-Regel limit ($`ea^2/\mathrm{}`$ 1 cm 2/Vs, where $`a`$ is the lattice constant) at room temperature. The very nature of the low-temperature ‘metallic’ state of these materials can be hardly understood within the framework of the canonical theory of metals. However, it has been known for many decades starting from the pioneering work by Landau  that the local lattice deformation can transform electrons into ‘self-trapped’ polarons in doped ionic insulators. Further extensive theoretical studies (for review see  have shown that polarons behave like heavy particles, and can be mobile with metallic conduction at sufficiently low temperatures. Under certain conditions  they form a polaronic Fermi liquid with some properties being very different from ordinary metals. One of them is the effective mass of polarons. It is well known that electrons can change their mass in solids due to the interactions with ions, spins, and themselves. The renormalized (effective) mass of electrons is independent of the ion mass $`M`$ in ordinary metals where the Migdal adiabatic approximation is believed to be valid. However, if the interactions between electrons and nuclear ions are strong and/or the adiabatic approximation is not applied electrons form polarons (quasiparticles dressed by lattice distortions), so their effective mass $`m^{}`$ will depend on $`M`$. This is because the polaron mass $`m^{}=m\mathrm{exp}(A/\omega )`$ , where $`m`$ is the band mass in the absence of the electron-phonon interaction, $`A`$ is a constant, and $`\omega `$ is a characteristic phonon frequency which depends on the masses of ions. Hence, there is a large isotope effect on the carrier mass in polaronic metals, in contrast to the zero isotope effect in ordinary metals.

The total exponent of the isotope effect on $`m^{}`$ is defined as $`\alpha _m^{}=d\mathrm{ln}m^{}/d\mathrm{ln}M_i`$ ($`M_i`$ is the mass of the $`i`$th atom in a unit cell). From this definition and the expression for the polaron mass $`m^{}`$ mentioned above, one readily finds

$$\alpha _m^{}=\frac{1}{2}\mathrm{ln}(m^{}/m).$$
(1)
Interestingly, the same isotope exponent is predicted for the inter-site bound pairs of polarons (i.e., bipolarons) . A large isotope effect on the effective mass of supercarriers has been observed in the high-temperature cuprate superconductors , which suggests that polarons and/or bipolarons should be involved in high-temperature superconductivity. Furthermore, in doped manganites, the polaronic nature of the carriers in the high-temperature paramagnetic insulating state has been shown by several independent experimental results such as the sign anomaly of the Hall effect , the Arrhenius behavior of the drift and Hall mobilities , the giant oxygen-isotope shift of the Curie temperature $`T_C`$ , and the electron paramagnetic resonance linewidth which is proportional to the conductivity .

On the other hand, the very nature of the charge carriers in the low-temperature metallic state of doped manganites has not been understood. In a theory of colossal magnetoresistance for ferromagnetic manganites  polarons are considered as the carriers even in the low-temperature metallic state, while others  believe that polaronic effects are not important at low temperatures. In particular, it was shown that, in low $`T_C`$ polycrystalline manganites ($`T_C<`$120K), the resistivity and $`T_C`$ are extremely sensitive to the external pressure  and chemical pressure  while the thermopower is small and nearly pressure independent . This unusual behavior in the low $`T_C`$ materials was taken as evidence for a vibronic electronic state , or may be explained in terms of a percolative conduction mechanism . In contrast, from isotope-effect studies , it was concluded that polaronic effects may be involved in the ferromagnetic metallic state of the low $`T_C`$ compounds. Moreover, in high $`T_C`$ materials, it is even less clear whether polaronic effects are important in the low temperature metallic state.

Theoretically it has been demonstrated that the coherent motion of small polarons can lead to metallic conduction at low temperatures . However, there have been no convincing and clear experiments to confirm this theoretical prediction. In order to provide decisive evidence for the ‘metallic’ polaronic liquid in doped ferromagnetic manganites, it is essential to study the isotope effects on the low-temperature kinetic and thermodynamic properties. These properties are dominated by the impurity and electron-electron scatterings, so the only source of the isotope effect might be the polaron mass dependence on the ion mass. For example, replacing 16O with 18O will increase the polaron mass, leading to the changes in the residual resistivity, and in other kinetic and thermodynamic properties. Here we provide clear evidence for the polaronic liquid in doped ferromagnetic manganites at low temperatures from theoretical and experimental studies of the isotope effects on kinetic and thermodynamic properties.

We first define and calculate three exponents: $`\beta _C=d\mathrm{ln}C_{el}/d\mathrm{ln}m^{}`$, $`\beta _\rho =d\mathrm{ln}\rho _o/d\mathrm{ln}m^{}`$, and $`\beta _S=d\mathrm{ln}S/d\mathrm{ln}m^{}`$ in the polaronic liquid. Here $`C_{el}`$ is the electronic specific-heat at low temperatures, $`\rho _o`$ is the residual resistivity, and $`S`$ is the low-temperature thermopower. The above definition is equivalent to $`C_{el}(m^{})^{\beta _C}`$, $`\rho _o(m^{})^{\beta _\rho }`$, and $`S(m^{})^{\beta _S}`$. We assume that the Fermi level (chemical potential) lies above the mobility edge, so polarons propagate in the Bloch states as heavy fermions  scattered by the short-range impurity potential. Since doped holes in manganites mainly reside on the oxygen orbitals , their bands are almost one-dimensional (1D) due to a large anisotropy in the hopping integrals (i.e., $`t_{pp\pi }=(1/4)t_{pp\sigma }`$ ). The polaronic mass enhancent factor could be larger along the directions with a much smaller hopping integral . This will increase the anisotropy further so that the polaron bands become more 1D-like. Then the polaronic propagator, $`G(k,E)=[Ek^2/2m^{}\mathrm{\Sigma }(k,E)]^1`$ is found analytically with the self-energy, which is momentum ($`k`$)-independent in the non-crossing approximation,

$`\sigma (ϵ)`$ $`=`$ $`{\displaystyle \frac{ϵ}{3}}\left({\displaystyle \frac{1i3^{1/2}}{2}}\right)\left[{\displaystyle \frac{1}{16}}+{\displaystyle \frac{ϵ^3}{27}}+\left({\displaystyle \frac{1}{256}}+{\displaystyle \frac{ϵ^3}{216}}\right)^{1/2}\right]^{1/3}`$ (2)
$``$ $`\left({\displaystyle \frac{1+i3^{1/2}}{2}}\right)\left[{\displaystyle \frac{1}{16}}+{\displaystyle \frac{ϵ^3}{27}}\left({\displaystyle \frac{1}{256}}+{\displaystyle \frac{ϵ^3}{216}}\right)^{1/2}\right]^{1/3}.`$ (3)
Here we introduce a dimensionless energy, $`ϵ=E/ϵ_0`$, and the self-energy $`\sigma (ϵ)=\mathrm{\Sigma }(E)/ϵ_0`$ using the scattering rate $`ϵ_0=(D^2m^{}/\mathrm{}^2)^{1/3}`$ as the energy unit. The constant $`D`$ is the second moment of the Gaussian white-noise potential, $`D=2(n_{im}v^2/a^2)`$, where $`n_{im}`$ is the impurity density, and $`v`$ is the coefficient of the $`\delta `$-function impurity potential. This self-energy is almost exact above the mobility edge as shown with the exact solution for a one-dimensional particle in a random white-noise potential .

The specific heat is proportional to the density of states (DOS) at the (dimensionless) Fermi level, $`\mu =E_F/ϵ_0`$. Assuming three degenerate spin-polarized oxygen bands, we find

$$\frac{C_{el}}{T}=\frac{2\pi k_B^2}{\mathrm{}a^2}\left(\frac{m^{}}{ϵ_0}\right)^{1/2}g(\mu )$$
(4)
with $`g(ϵ)=\mathrm{}\sigma (ϵ)`$. The residual conductivity $`\rho _o^1`$ and the low-temperature thermopower $`S`$ are found with the Kubo formalism  as

$$\rho _o^1=\frac{e^2}{2\pi ^2(m^{})^2}_{\mathrm{}}^{\mathrm{}}𝑑kk^2\left[\mathrm{}G(k,\mu )\right]^2,$$
(5)
and

$$S=\frac{\pi ^2k_B^2T}{3eϵ_0}\frac{d\mathrm{ln}\rho _o}{d\mu }.$$
(6)
If there coexist electron-like and hole-like polarons, Eq. 5 must be multiplied by a factor $`f_S`$ with $`|f_S|1`$. We would like to mention that Eq. 4 should be underestimated due to the fact that the bands are not strictly one-dimensional and that the transport relaxation rate is generally smaller than $`ϵ_0`$. Calculating the integral in Eq. 4 leads to

$$\rho _o^1=\frac{e^2}{2a^2\pi (m^{}ϵ_0)^{1/2}}F(\mu ),$$
(7)
where

$$F(\mu )=\frac{\stackrel{~}{\mu }[\stackrel{~}{\mu }+(\stackrel{~}{\mu }^2+g^2)^{1/2}]+g^2/2}{g[(\stackrel{~}{\mu }^2+g^2)(\stackrel{~}{\mu }+(\stackrel{~}{\mu }^2+g^2)^{1/2})]^{1/2}},$$
(8)
$`\stackrel{~}{\mu }=\mu \mathrm{}\sigma (\mu )`$ and $`g=g(\mu )`$.

These expressions allow us to calculate the exponents: $`\beta _C`$, $`\beta _\rho `$, and $`\beta _S`$. The number of carriers does not depend on the isotope substitution , so that the chemical potential itself depends on the polaron mass with the exponent

$$\frac{d\mu }{d\mathrm{ln}m^{}}=\frac{2N}{3g},$$
(9)
where $`N=_{\mathrm{}}^\mu g(ϵ)𝑑ϵ`$ is proportional to the total number of states under the Fermi level. Keeping that in mind we finally obtain

$$\beta _C=\frac{1}{3}\left(1\frac{2Ng^{}}{g^2}\right),$$
(10)
$$\beta _\rho =\frac{2}{3}\left(1+\frac{NF^{}}{gF}\right),$$
(11)
and

$$\beta _S=\frac{1}{3}\left(\frac{2N(F^{}/FF^{\prime \prime }/F^{})}{g}1\right).$$
(12)
The prime and the double primes mean the first and the second derivative with respect to $`\mu `$. The numerical results for the exponents are shown in Fig. 1 as a function of the Fermi level $`E_F/ϵ_0`$. To compare with the experimentally measured oxygen-isotope isotope exponents one should multiply $`\beta _C,\beta _\rho ,\beta _S`$ by the oxygen-isotope exponent $`\alpha _m^{}^O`$. Because of the logarithmic dependence on the mass ratio (see Eq. 1), the magnitude of the exponent $`\alpha _m^{}`$ is expected to be of the order of unity. At large filling, i.e., $`\mu >>1`$, the theory predicts sizable isotope effects on all thermodynamic and kinetic properties with the exponents:

$$\beta _C=1;\beta _\rho =2;\beta _S=1.$$
(13)
This is precisely what one expects from the simple Born approximation for the impurity scattering. The ratio remains the same for any dimension of the polaron spectrum in this doping regime. Interestingly, if the doping is not so high, $`\mu 1`$, the exponent for the thermopower is very small (see Fig. 1), while the exponent for the residual resistivity remains larger than unity. For $`\mu =0.25`$ and $`\alpha _m^{}^O\alpha _m^{}1.1`$,

$$\beta _C0.45;\beta _\rho 1.45;\beta _S0.03,$$
(14)
and

$$\alpha _C^O0.5;\alpha _\rho ^O1.6;\alpha _S^O0.0,$$
(15)
i.e., negligible isotope effect on $`S`$, and more than 3 times lower isotope effect on $`C_{el}`$ than on $`\rho _o`$. Applying the exact density of states in the random potential , we find even a somewhat lower exponent $`\beta _C`$ for $`\mu 0`$ than that in Fig. 1. All these findings are in sharp contrast with those in ordinary metals where $`\alpha _m^{}=0`$, and there are no isotope effects on $`C_{el},\rho _o`$ and $`S`$.

The choice of $`\mu =0.25`$ and $`\alpha _m^{}^O=1.1`$ is consistent with the experimental results for Nd0.67Sr0.33MnO3. Using the measured bare plasma frequency $`\mathrm{}\mathrm{\Omega }_p`$ $``$ 3.3 eV  and the total carrier density per cell of 0.33, we obtain the bare hopping integral $`t=t_{pp\sigma }=\mathrm{}^2/2a^2m`$ = 0.9 eV, which is very close to the value (0.87 eV) calculated for La2CuO4 . From Eq. 1 and $`\alpha _m^{}\alpha _m^{}^O1.1`$, we get the effective hopping integral $`t^{}=\mathrm{}^2/2a^2m^{}`$ $``$ 100 meV. By setting $`\mu `$ = 0.25, $`t^{}`$ = 100 meV, $`ϵ_0`$ = 27 meV, we calculates $`C_{el}/T`$ = 26 mJ/moleK2 and $`\rho _o`$ = 185 $`\mu \mathrm{\Omega }`$ cm by using Eq. 3 and Eq. 6, respectively. These values are in remarkably good agreement with the measured values: $`C_{el}/T`$ = 25(3) mJ/moleK2 , and $`\rho _o`$ $``$ 170 $`\mu \mathrm{\Omega }`$ cm for the best single crystals . Furthermore, the value of $`ϵ_0`$ is in good agreement with the measured transport relaxation rate ($``$24 meV) for La0.7Ca0.3MnO3 . This justifies the choice of $`\mu =0.25`$ and $`\alpha _m^{}^O=1.1`$ for Nd0.67Sr0.33MnO3. For La0.75Ca0.25MnO3, we calculate $`C_{el}/T`$ = 7.5 mJ/moleK2 and $`\rho _o`$ = 150 $`\mu \mathrm{\Omega }`$ cm by setting $`\mu =0.70`$ and $`\alpha _m^{}^O=0.7`$. The measured $`\rho _o`$ for La0.75Ca0.25MnO3 is 120-150 $`\mu \mathrm{\Omega }`$ cm (see Fig. 2a and discussion in Ref. ). From our preliminary specific-heat data for La0.75Ca0.25MnO3, we estimate $`C_{el}/T`$ $``$ 7 mJ/moleK2, which is slightly lower than the one ($``$8 mJ/moleK2) for La0.8Ca0.2MnO3 . With $`\mu =0.70`$ and $`\alpha _m^{}^O=0.7`$, we have

$$\alpha _C^O0.5;\alpha _\rho ^O1.2;\alpha _S^O0.17,$$
(16)
Experimentally, the oxygen-isotope effect on the intrinsic low-temperature resistivity has been studied  in high-quality thin films of La0.75Ca0.25MnO3 and Nd0.7Sr0.3MnO3. Our thin films were thick (2000 $`\AA `$), and were annealed at 940oC for 10 hours, so no strains could remain in the films. We did back-exchange, and we knew precisely the oxygen-isotope enrichment (more details can be found in Ref. ). In Fig. 2 we plot the low-temperature resistivity of the oxygen-isotope exchanged films of (a) La0.75Ca0.25MnO3; (b) Nd0.7Sr0.3MnO3. In both cases, the residual resistivity $`\rho _o`$ for the 18O samples is larger than for the 16O samples by about 15$`\%`$. The observed large oxygen-isotope effect on $`\rho _o`$ is consistent with the above theoretical prediction for polaronic metals. From the definition of the oxygen-isotope exponent: $`\alpha _\rho ^O`$ = $`d\mathrm{ln}\rho _o/d\mathrm{ln}M_O`$ ($`M_O`$ is the oxygen mass), we obtain $`\alpha _\rho ^O=1.4(2)`$ for Nd0.7Sr0.3MnO3, and $`\alpha _\rho ^O=1.2(2)`$ for La0.75Ca0.25MnO3. These values are in excellent agreement with Eq. 14 and Eq. 15, respectively.

The quantitative agreement between the resistivity data and the theory indicates that the low-temperature metallic state in doped manganites is a polaronic Fermi-liquid. To further confirm this, one should experimentally demonstrate that $`\alpha _S^O`$ is negligible and that $`\alpha _C^O`$ is rather small in these materials. Our preliminary results of the oxygen-isotope effect on the low-temperature specific heat in La1-xCaxMnO3 are consistent with Eq. 15 within the experimental uncertainty. The results will be published elsewhere. Here we show the data of the thermopower $`S`$ at low temperatures for the isotope-exchanged La1-xCaxMnO3. The detailed measurement technique was described in Ref. . Fig. 3a shows the thermopower $`S`$ below 60 K for the 16O and 18O samples of La1-xCaxMnO3 with $`x`$ = 0.20 and 0.25. The isotope-exchanged samples are the same as those reported in . It is apparent that the oxygen-isotope effect on $`S`$ is negligible below 20 K for both compositions ($`\alpha _S^O0.0`$), as seen more clearly from Fig. 3b. This is in sharp contrast to the large oxygen-isotope shift of the Curie temperature TC, which is 21 K for $`x`$ = 0.20 and 14.5 K for $`x`$ = 0.25 . The fact that the slope of the linear part in S (see Fig. 3) is very sensitive to the doping or the carrier concentration and that the 16O and 18O samples have the same slope implies that the carrier concentrations for two isotope samples of a fixed $`x`$ should be nearly the same, in agreement with the previous results . The negligible $`\alpha _S^O`$ and the substantial $`\alpha _\rho ^O`$ $``$ -1.3 observed in these samples are in quantitative agreement with the theoretical prediction (Eq. 14 and Eq. 15). We are not aware of any other theoretical models which can quantitatively explain the experimental values of several physical quantities, such as electronic specific heat, residual resistivity, scattering rate, and three different oxygen-isotope exponents. Therefore, the present theoretical and experimental studies of the low-temperature kinetic and thermodynamic properties provide strong evidence for the existence of novel polaronic Fermi-liquid in doped metallic oxides (e.g., magnetoresistive manganites and high-temperature superconductive cuprates).
Acknowlegdment: We would like to thank K. A. Müller for valuable comments. The work at Zürich was supported by the Swiss National Science Foundataion. The work at Houston was supported in part by the US National Science Foundation, and work at Berkeley by the US Department of Energy.

Figure captions

Fig.1. The numerical results for three exponents: $`\beta _\rho `$, $`\beta _C`$, and $`\beta _S`$ as a function of the Fermi level $`\mu `$ $`E_F/ϵ_0`$. See text for the definition of the exponents.

Fig.2. The low-temperature resistivity of the oxygen-isotope exchanged films of (a) La0.75Ca0.25MnO3; (b) Nd0.7Sr0.3MnO3. In both cases, the residual resistivity $`\rho _o`$ for the 18O samples is larger than for the 16O samples by 15$`\%`$. The oxygen-isotope exponent $`\alpha _\rho ^O`$ was calculated from the definition: $`\alpha _\rho ^O`$ = $`d\mathrm{ln}\rho _o/d\mathrm{ln}M_O`$ ($`M_O`$ is the oxygen mass).

Fig.3. The temperature dependence of the thermopower $`S`$ of the 16O and 18O samples of La1-xCaxMnO3 below 60 K (Fig. 3a) and below 20 K (Fig. 3b). The solid and dash lines in Fig. 3b are the fitted curves by $`S=S_o+BT`$ with a constraint that the absolute value of $`S_o`$ should be less than the systematic measurement error 0.1 $`\mu `$V/K (Ref. ). The magnitudes of the exponent $`\alpha _S^O`$ were estimated from the values of the fitting parameter $`B`$.


text:
# Dynamical Interpretation of Chemical Freeze-Out in Heavy Ion Collisions

## Abstract

It is demonstrated that there exists a direct correlation between chemical freeze-out point and the softest point of the equation of state where the pressure divided by the energy density, $`p(\epsilon )/\epsilon `$, has a minimum. A dynamical model is given as an example where the passage of the softest point coincides with the condition for chemical freeze-out, namely an average energy per hadron $``$ 1 GeV. The sensitivity of the result to the equation of state used is discussed.

Over the last few years the question of chemical equilibrium in heavy ion collisions has attracted much attention . Assuming thermal and chemical equilibrium within a statistical model, it has now been shown that it is indeed possible to describe the hadronic abundances produced at beam energies ranging from 1 to 200 AGeV. The observation was made that the parameters of the chemical freeze-out curve obtained at CERN/SPS, BNL/AGS and GSI/SIS all lie on a unique freeze-out curve in the temperature $`T`$ – chemical potential $`\mu _B`$ plane . Recently, a surprisingly simple interpretation of this curve has been proposed: the hadronic composition at the final state is determined solely by the average energy per hadron, $`E_{had}/N_{had}`$, of approximately 1 GeV in the rest frame of the system under consideration . Using this finding $`E_{had}/N_{had}1`$ GeV as a heuristic definition for the chemical freeze-out point, in this letter we study physical conditions which are realized around this point in evolution of the system. We show that the chemical freeze-out point is intimately related to the softest point of equation of state defined by the minimum of the pressure-to-energy density ratio, $`p(\epsilon )/\epsilon `$, as a function of $`\epsilon `$ .

Our considerations are essentially based on the recently proposed Mixed Phase (MP) model  which is consistent with the available QCD lattice data . The underlying assumption of the MP model is that unbound quarks and gluons may coexist with hadrons forming an homogeneous quark/gluon–hadron phase . Since the mean distance between hadrons and quarks/gluons in this mixed phase may be of the same order as that between hadrons, their interaction with unbound quarks/gluons plays an important role defining the order of the phase transition. Recently the importance of quark/gluon–hadron interactions was also discussed in the context of dilepton production .

Within the MP model  the effective Hamiltonian is written in the quasiparticle approximation with the density-dependent mean–field interaction. Under quite general requirements of confinement for color charges, the mean–field potential of quarks and gluons is approximated by the following form:

$$U_q(\rho )=U_g(\rho )=\frac{A}{\rho ^\gamma }$$
(1)
with the total density of quarks and gluons

$$\rho =\rho _q+\rho _g+\underset{j}{}n_j\rho _j$$

where $`\rho _q`$ and $`\rho _g`$ are the densities of unbound quarks and gluons outside hadrons, while $`\rho _j`$ is the density and $`n_j`$ is the number of valence quarks inside the hadron of type $`j`$. The presence of the total density $`\rho `$ in (1) describes the interaction between all components of the mixed phase. The approximation (1) recovers two important limiting cases of the QCD interaction. Namely, if $`\rho 0`$ the interaction potential goes to infinity, i.e. an infinite energy is required to create an isolated quark or gluon which ensures the confinement of color objects. In the other extreme of high energy density corresponding to $`\rho \mathrm{}`$, we obtain the asymptotic freedom regime.

The use of the density-dependent potential (1) for quarks and the hadronic potential described by a modified non-linear mean–field model  requires certain constraints to satisfy thermodynamic consistency . For the chosen form of the Hamiltonian these conditions show that $`U_g(\rho )`$ and $`U_q(\rho )`$ should be independent of the temperature. From these requirements one also obtains an expression for the quark–hadron potential . In the case when the quark-gluon component is neglected, the MP model is reduced to the interacting hadron gas model with a non-linear mean–field interaction .

A detailed study of the pure gluonic medium with $`SU(3)`$ color symmetry and a first order phase transition allows to fix the values of $`\gamma =0.62`$ and $`A^{1/(3\gamma +1)}=250MeV`$. These values are then generalized to the $`SU(3)`$ system with dynamical quarks. For two light flavors at zero baryon density, $`n_B=0`$, the MP model is consistent with the results from lattice QCD  with the deconfinement temperature $`T_{dec}=153MeV`$ and the crossover type of phase transition. The model is then proposed to be extended to baryon-rich systems in a parameter–free way .

A particular consequence of the MP model is that for $`n_B=0`$ the softest point  of the equation of state is located at comparatively low values of the energy density: $`\epsilon _{SP}0.45GeV/fm^3`$ which is in a good agreement with recent lattice estimates . This value of $`\epsilon _{SP}`$ is close to the energy density inside a nucleon. Thus, reaching this value signals that we are dealing with a single ’big’ hadron consisting of deconfined matter. For baryonic matter the softest point is gradually washed out at $`n_B\text{ }>0.4n_0`$. As shown in , this thermodynamic behavior differs drastically from both the interacting hadron gas model which has no softest point and the two–phase approach, based on the bag model, having by construction a first order phase transition and the softest point at $`\epsilon _{SP}>1GeV/fm^3`$ independent of $`n_B`$ . These differences should manifest themselves in the expansion dynamics.

In Fig.1 we show the evolution trajectories for central $`Au+Au`$ collisions in the $`T\mu _B`$ plane together with the freeze-out parameters obtained from hadronic abundances. Using the quark-gluon string model , the initial state was estimated as a state of hot and dense nuclear matter inside a cylinder in the center-of-mass frame with radius $`R=4fm`$ and length $`L=2R/\gamma _{c.m.}`$. The calculated temporal behavior of energy and baryon densities is close to that in RQMD or UrQMD transport models and the procedure of defining the initial state of a fireball is described in . The subsequent isentropic expansion of this fireball is treated within the Frankfurt expansion model  based on a scaled hydrodynamical prescription assuming $`tV`$. As seen in Fig.1, the turning points of these trajectories, i.e. the points where $`T/\mu _B`$ changes the sign, are in a good agreement with the extracted chemical freeze-out parameters and clearly correlate with the smooth curve for the ideal gas model with $`E_{had}/N_{had}=1`$ GeV  defined above as a condition of chemical freeze-out. Existence of the turning point is related to two limiting equilibrium regimes for temperature behavior of the chemical potential: ultrarelativistic regime dominated by light particles (quarks, pions) when $`\mu _BT^2`$ and non-relativistic one defined mainly by nucleons and deltas with $`\mu _BT^1`$
\- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Footnote: Other approach have been recently applied in  where conditions for thermodynamic equilibrium in heavy ion collisions were studied within the transport UrQMD model and the calculation results at a fixed time moment were approximated by the ideal hadron gas equation of state to extract the time dependence of $`T`$ and $`\mu _B`$. It is of interest that our results for dynamical trajectories below the turning point practically coincide with that from  for equilibrium part of $`T\mu _B`$ trajectory. However, at higher temperature, where according to  the system is in a non-equilibrium state, their results differ from ours and exhibit no turning point.
\- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
. Uncertainty in the initial $`\epsilon `$ and $`n_B`$ in 10-20$`\%`$ shifts slightly the dynamical trajectory but the turning point nevertheless stays on the $`E_{had}/N_{had}=1`$ GeV curve. Notice that there is a finite region in $`T`$ where $`\mu _B/T0`$, i.e. the chemical potential in the expanding system is kept constant around this turning point.

The observed correlation is further elucidated in Fig.2. The quantity $`p/\epsilon `$ is closely related to the square of the velocity of sound and characterizes the expansion speed
\- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Footnote: In simplified hydrodynamic models was shown that, for example, the transverse expansion of a cylindrical source is governed by the pressure-to-enthalpy ratio, $`p/(p+\epsilon )`$ , i.e. by $`p/\epsilon `$ rather than by the sound velocity.
\- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
, so the system lives for an appreciable amount of time around the softest point which facilitates to reach the chemical equilibrium. This statement is evidenced directly in the right-hand side of Fig.2 where the time evolution of $`p/\epsilon `$ is shown: The system spends about 1/3 of the total expansion time in a state near the softest point. During this time the baryon chemical potential $`\mu _B`$ remains practically constant.

In Fig.2 (left-hand side) the position of the softest point correlates with the average energy per hadron being about 1 GeV for all beam energies except for 2 AGeV case. One should note that the quantity $`\epsilon /\rho _{had}`$, where $`\epsilon `$ is the total energy density, coincides with $`E_{had}/N_{had}`$ considered in  only in the case when there are no unbound quarks/gluons in the system. In the MP model, all components are interacting with each other and therefore the quantity $`E_{had}`$ is not well defined. However the admixture of unbound quarks at the softest point $`\epsilon _{SP}`$ amounts to about $`13\%`$ and $`8\%`$ at beam energies $`E_{lab}=150`$ and $`10AGeV`$, correspondingly. Thus thermodynamical quantities and in particular the ratios of hadron abundance in the MP and resonance gas model are very close to each other at the freeze-out point.

The MP equation of state plays a decisive role for the regularity considered here, describing both the order of the phase transition and the deconfinement temperature. The two-phase (bag) model exhibits a first order phase transition with $`T_{dec}=160MeV`$ and has a spatially separated Gibbs mixed phase but the corresponding trajectories in the $`T\mu _B`$ plane are quite different from those in the MP model as shown in . The exit point from the Gibbs mixed phase at $`E_{lab}=150`$ AGeV is close to the corresponding freeze-out point in Fig.1. However it does not lead to the observed correlation with the softest point position in the whole energy range considered. The interacting hadron gas model has no softest point effect as was demonstrated in  but nevertheless it exhibits the turning point in $`T\mu _B`$ plane. These facts are seen also from Fig.1 and Fig.2 where for $`E_{lab}=2AGeV`$ quark-hadron interactions are negligible ($`1\%`$). Indeed, at this energy instead of a minimum one has a monotonic fall-off specific for hadronic models with only a small irregularity in $`p/\epsilon `$ near the point $`\epsilon /\rho _{had}1GeV`$
\- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Footnote: Note that at the SIS energies the chemical freeze-out point practically coincides with the thermal freeze-out .
\- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
. At higher energies the dynamical trajectories are very different for hadronic and mixed phase models but the positions of turning points in the $`T\mu _B`$ plane turn out to be still rather close to each other and both lying on the curve corresponding to $`\epsilon /\rho _{had}=1GeV`$. Therefore, one could expect similar results for the ratios of hadronic abundance in these two models. However the model predictions differ essentially as long as evolution of the system is concerned.

It is noteworthy that similarly to the results presented in Fig.2, the softest point of the equation of state correlates with an average energy per quark, $`\epsilon /\rho 350MeV`$ which is close to the constituent quark mass. So, at higher values of $`\epsilon /\rho `$ we are dealing with a strongly-interacting mixture of highly-excited hadrons and unbound massive quarks/gluons forming (in accordance with Landau’s idea ) an ’amorphous’ fluid suitable for hydrodynamic treatment. Below the softest point the interaction deceases, the relative fraction of unbound quarks/gluons is getting negligible, higher hadronic resonances decay into baryons and light mesons. At this stage the application of hydrodynamics becomes questionable.

In summary, we have shown in this paper an example of an equation of state where there exists a direct correlation between chemical freeze-out point and the softest point of the equation of state. This description correlates with the observation that chemical freeze-out occurs when the average energy per hadron, $`E_{had}/N_{had}`$, drops below the value of 1 GeV as found in . In the examples considered such a clear correlation is observed for the equation of state of the mixed phase model, where the energy per hadron seems to be a relevant observable defining chemical freeze-out conditions.

We thank B. Friman, Yu. Ivanov, E. Kolomeitsev and W. Nörenberg for useful discussions. E.G.N., K.R. and V.D.T. gratefully acknowledge the hospitality at the Theory Group of GSI, where some part of this work has been done. J.C. gratefully acknowledges the hospitality of the Physics Department of the University of Bielefeld. This work was supported in part by DFG under the program of scientific-technological collaboration (project 436 RUS 113/558/0) and by Russian Found of Fundamental Investigations (project 00-02-04012). K.R also acknowledges the partial support of the State Committee for Scientific Research (KBN).