text:
# Pendulum limit, chaos and phase-locking in the dynamics of ac-driven semiconductor superlattices

## 1 Introduction

Strongly nonlinear dynamical effects in miniband semiconductor superlattices (SSLs) driven by an intense high-frequency electric field attract growing interest in the last few years . This theoretical activity, which includes investigation of dissipative chaos  and a spontaneous generation of dc voltage in a purely ac-driven SSL , has been stimulated by a recent progress in the experimental studies of nonlinear electron transport in the SSLs driven by THz field . On other hand, recent experiments  reveal an important role of coherent plasma oscillations in nonlinear dynamics of carriers in SSLs.

The counterparts of the effects of chaos and spontaneously dc voltage generation in the THz-driven SSL are known in the rf-driven Josephson junctions for a rather long time  (for a review, see ). So-called an inverse ac Josephson effect (a quantized dc voltage generation)  has already found applications in the design of a zero-biased Josephson voltage standard , which offers several advantages in comparison with the conventional voltage standards based on an ac Josephson effect . A transition to dissipative chaos often limits the stable operation of such kind of voltage standard . Nonlinear dynamics of driven plasma oscillations in a junction is well-described by the fairly popular Resistively Shunted Junction (RSJ) model , which is equivalent to a periodically forced and damped pendulum – one of the most important paradigms in the modern nonlinear science .

Within the semiclassical approach , ac-driven plasma oscillations in a single miniband SSL are described by the set of balance equations for the average electron velocity (current), the average electron energy and the electric field (voltage) generated by an electron motion in SSL . It has been shown in ref.  that this set of three nonautonomous nonlinear first-order differential equations is isomorphic to the periodically forced Lorenz model . In comparison with the RSJ model, this system is more complex from a viewpoint of nonlinear dynamics and until recently it has been analysed mainly numerically .

However, it should be noticed that a pendulum-like behavior of miniband electrons in SSL has been established earlier for some limiting cases. First, plasma oscillations in SSL have been described within a collisionless approximation by a pendulum without damping in . Second, describing the frequency-to-voltage conversion in ac-driven SSL, Dunlap et al.  have found that their model of SSL demonstrates a behavior similar to the phase-locked solutions of damped pendulum  for $`\omega \omega _{\mathrm{pl}}`$ ($`\omega `$ and $`\omega _{\mathrm{pl}}`$ are the ac electric field frequency and the miniband plasma frequency, correspondently). These earlier findings indicate that dynamics governed by a pendulum can be underlying for nonlinear dynamics of ac-driven SSL.

In this work, we show that the driven damped pendulum could describe a stationary nonlinear dynamics of an ac-driven SSL in the limit when a characteristic scattering constant of electrons with impurities and phonons, $`\gamma `$, is less than an an ac frequency $`\omega `$. In contrast to the RSJ model, where the voltage across junction is proportional to the velocity of pendulum, the voltage across SSL depends on both the coordinate and the velocity of pendulum. These findings possess us to apply several exact and approximate results obtained earlier for a pendulum to an analytic investigation of regular and chaotic dynamics in SSL. Namely, we will show the following:
(1) Chaotic dynamics is natural for a relatively low damping, i.e. for $`\gamma /\omega _{\mathrm{pl}}1`$.
(2) The SSLs with wide minibands are more stable against perturbations leading to a transition to chaos in comparison to the case of narrow miniband SSLs.
(3) The Melnikov criterion of transition to chaos in the pendulum qualitatively explains the position of low frequency boundary of transition to chaos in SSL.
(4) For a relatively low damping, $`\gamma /\omega _{\mathrm{pl}}1`$, phase-locking in driven pendulum results in a generation of dc bias across of SSL; the dc voltage per period of SSL is approximately proportional to a half-integer multiplier of $`\mathrm{}\omega `$.
(5) The difference between generated dc voltage per superlattice’s period and $`\mathrm{}\omega `$ is proportional to $`(\gamma /\omega _{\mathrm{pl}})^2`$ at $`\gamma /\omega _{\mathrm{pl}}1`$.
(6) For a fairly strong scattering, phase-locking in driven pendulum may not result in a generation of quantized (or an approximately quantized) dc bias across SSL.
Our present analytic research explains and further develops the results obtained numerically in the previous publications .

## 2 Pendulum limit in superlattice balance equations

We consider electron transport through a single miniband of a spatially homogeneous SSL under the influence of ac current. We assume the tight-binding dependence $`\epsilon (k)=\mathrm{\Delta }/2[1\mathrm{cos}(ka)]`$ of the electron energy $`\epsilon `$ on its quasimomentum $`k`$ along SSL’s axis, where $`\mathrm{\Delta }`$ is the miniband width, and $`a`$ is the period of SSL. The electron dynamics could be described by the following balance equations  (we use notations of ref. )

$`\dot{v}`$ $`=`$ $`uw\gamma v,`$
$`\dot{w}`$ $`=`$ $`uv\gamma (ww_{\mathrm{eq}}),`$ (1)
$`\dot{u}`$ $`=`$ $`\omega _0^2v+I_0\mathrm{sin}\omega t,`$
Here the scaled variable $`u=eaE(t)/\mathrm{}`$ describes the electric field $`E`$ (or the voltage per SSL’s period), generated by an electron motion along the SSL axis; $`v=m_0\overline{V}a/\mathrm{}`$ is the scaled electron velocity ($`\overline{V}`$ is the electron velocity averaged over the time-dependent distribution function satisfying the Boltzmann equation, $`m_0=\frac{2\mathrm{}^2}{\mathrm{\Delta }a^2}`$ is the effective electron mass at the bottom of miniband) and $`w=(\overline{\epsilon }\mathrm{\Delta }/2)(\mathrm{\Delta }/2)^1`$, where $`\overline{\epsilon }`$ is the average electron energy within the miniband. Following this scaling for the energy, the lower (upper) edge of the miniband corresponds to $`w=1`$ ($`w=+1`$). Next, $`\omega `$ is the frequency and $`I_0`$ is the scaled amplitude of driving ac current; $`\gamma `$ is the phenomenological relaxation constant describing different channels of electron scattering in SSL; $`\omega _0=\omega _{\mathrm{pl}}\left(C_{\mathrm{sl}}/C\right)^{1/2}`$ is the generalized plasma frequency, where $`\omega _{\mathrm{pl}}=\left(4\pi e^2N/m_0ϵ_0\right)^{1/2}`$ is the miniband plasma frequency, $`N`$ is the number of electrons per unit volume, $`ϵ_0`$ is an average dielectric constant for the SSL material, $`C_{\mathrm{sl}}=\frac{ϵ_0S}{4\pi l}`$ is the self-capacitance of SSL with the length $`l`$ and the cross-section area $`S`$, and $`C`$ is a parasitic capacitance, which is in parallel to the SSL’s capacitance. The electron energy within miniband relaxes to the thermal equilibrium energy value, $`w_{\mathrm{eq}}`$ ($`w_{\mathrm{eq}}0`$); in the case of a nondegenerate electron gas, $`w_{\mathrm{eq}}`$ has the following temperature dependence

$$w_{\mathrm{eq}}=\frac{I_1(\mathrm{\Delta }/2k_BT)}{I_0(\mathrm{\Delta }/2k_BT)},$$
(2)
where $`I_{0,1}`$ are the modified Bessel functions, $`T`$ is the lattice temperature, and $`k_B`$ is the Boltzmann’s constant. Note that the dependence (2) well describes temperature-induced modifications of the Drude conductivity in the narrow miniband SSLs as have been observed in the experiments .

First two equations of the set (1) are the well-known balance equations of Ignatov and Romanov , which could be derived from the Boltzmann transport equation with a constant relaxation time $`\gamma `$. Thus, instead of finding the time-dependence of distribution function $`f(p,t)`$ $`(p=\mathrm{}k)`$ from the Boltzmann equation and further calculation of the average electron velocity $`\overline{V}=𝑑pf(p,t)\epsilon (p)/p`$ and the average electron energy $`\overline{\epsilon }=𝑑pf(p,t)\epsilon (p)`$, it is possible just to solve the corresponding balance equations for the average variables. Both these approaches give the same result for the tight-binding dispersion relation.

The third equation in (1) could be obtained from the Maxwell or Kirchoff equations . To the best our knowledge, three coupled balance equations in the form (1) have appeared first in the work of Tetervov , who used these equations for a description of decay of undriven ($`I_0=0`$) nonlinear plasma oscillations in SSL.

Now we want to show how the pendulum limit arises in the superlattice balance equations (1). We introduce new variables $`A`$ and $`\theta `$ as $`w=A\mathrm{cos}\theta `$, $`v=A\mathrm{sin}\theta `$. We can re-write Eqs. (1) in the form (see Appendix A)

$`\ddot{\theta }`$ $`+`$ $`A\omega _0^2\mathrm{sin}\theta +\gamma G(\theta ,A)\dot{\theta }=`$ (3)
$`I_0\mathrm{sin}\omega t+\gamma ^2w_{\mathrm{eq}}A^1F(\theta ,A),`$
$`\dot{A}`$ $`=`$ $`\gamma A\gamma w_{\mathrm{eq}}\mathrm{cos}\theta ,`$ (4)
where

$`G(\theta ,A)`$ $`=`$ $`w_{\mathrm{eq}}A^1\mathrm{cos}\theta ,`$ (5)
$`F(\theta ,A)`$ $`=`$ $`\mathrm{sin}\theta +w_{\mathrm{eq}}A^1\mathrm{sin}\theta \mathrm{cos}\theta ,`$ (6)
$`\dot{\theta }(t)`$ $`=`$ $`u(t)+\gamma w_{\mathrm{eq}}A^1(t)\mathrm{sin}\theta (t).`$ (7)
Formally, the set of equations (3)-(7) is not simpler than the original balance equations (1). However, this form is more suitable for the consideration of pendulum limits.

We start our analysis with the case of collisionless SSL, $`\gamma =0`$. As follows from Eq. (4), $`A^2=v^2(t)+w^2(t)`$ is the integral of motion in this case. For natural initial conditions $`v(0)=0`$, $`w(0)=w_{\mathrm{eq}}`$, we have $`A=w_{\mathrm{eq}}`$. The balance equations are reduced to single equation of undamped and periodically forced pendulum

$$\ddot{\theta }+(w_{\mathrm{eq}})\omega _0^2\mathrm{sin}\theta =I_0\mathrm{sin}\omega t,$$
(8)
where $`v(t)=w_{\mathrm{eq}}\mathrm{sin}\theta (t)`$, $`w(t)=w_{\mathrm{eq}}\mathrm{cos}\theta (t)`$, $`u(t)=\dot{\theta }(t)`$. Earlier the same equation has been obtained within the semiclassical Hamiltonian approach by several authors . It was used in the investigations of nonlinear plasma oscillations , the instabilities and the Hamiltonian chaos  in ideal ac-driven SSLs.

Eq. (8) is valid only for a short time $`t\gamma ^1`$. Now we want to show that a pendulum equation can also describe electron’s dynamics in the ac-driven SSL at long time ($`t\gamma ^1`$) providing $`\omega \gamma `$. Eq. (3) represents a pendulum, parameters of which are varied self-consistently due to their dependences on $`A`$ and $`\theta `$. To begin with we limit our attention to the time interval $`T=2\pi /\omega `$. During this characteristic period of oscillations (or rotations), $`T`$, the variable $`A`$ only lightly decays for $`\omega \gamma `$ (see Eq. (4)). Now we compare characteristic damping terms in Eqs. (3) and (4). Following Eqs. (3) and (5), the averaged over period of ac field dissipative term is

$$\gamma G(\theta ,A)=\gamma w_{\mathrm{eq}}A^1\mathrm{cos}\theta $$

Thus, the ratio of characteristic dissipative terms in Eqs. (3) and in Eq. (4), calculated during the time interval $`T`$, is following

$$\frac{\gamma G(\theta ,A)}{\gamma }=w_{\mathrm{eq}}A^1\mathrm{cos}\theta 1$$

This ratio is small because the term $`\mathrm{cos}\theta `$ changes its signs (oscillates), while $`A`$ stays almost constant. Next, we consider more long time interval of more number of periods $`T`$ and come to a conclusion that a dissipation in Eq. (4) is, in effect, much more strong than a dissipation in Eq. (3); the variable $`A(t)`$ approach its stationary value faster than the variable $`\theta (t)`$ approaches its attractor. Therefore, if additionally the amplitude of driving current is not large, $`I_0/\omega _0^2<1`$, the variable $`A`$ can be adiabatically eliminated . That is, we can substitute the stationary value $`A(t)A_{\mathrm{st}}(t)=w_{\mathrm{eq}}\mathrm{cos}\theta `$ in Eq. (3). Observing that $`G(\theta ,A_{\mathrm{st}})=1`$ and $`F(\theta ,A_{\mathrm{st}})=0`$, we find that Eqs. (3), (7) take the form

$$\ddot{\psi }+\gamma \dot{\psi }+(w_{\mathrm{eq}})\omega _0^2\mathrm{sin}\psi =2I_0\mathrm{sin}\omega t,$$
(9)
$$u=\dot{\psi }/2+\gamma \mathrm{tan}(\psi /2),$$
(10)
where we introduced $`\psi (t)=2\theta (t)`$. This pendulum representation of superlattice balance equations constitutes the main result of this paper.

It is useful to compare some well-known solution of Boltzmann equation for the tight-binding lattice and the corresponding solution of pendulum equation (9). This test can be simplified employing an analogy with the RSJ model of Josephson junctions, which is also governed by a periodically driven and damped pendulum. Consider a SSL driven by a given ac electric field $`E=E_0\mathrm{cos}\omega t`$ with the frequency being much greater than the miniband plasma frequency, $`\omega \omega _0`$. In this case we have $`I_0=\omega _s\omega `$ with $`\omega _s=eE_0a/\mathrm{}`$ . For such $`I_0`$ and $`\omega `$, we find that the solution of pendulum equation for the variable $`v(t)=\mathrm{sin}\theta \mathrm{cos}\theta `$ is same as the well-known expression for the time-dependent electron’s velocity $`v(t)`$ derived from the exact solution of the Boltzmann equation , if $`\omega \gamma `$ and $`I_0/\omega _0^2<1`$ (see Appendix B). Therefore, we see that the pendulum representation, Eqs. (9)-(10), gives same result as the Boltzmann equation within the range of its validity. Moreover, this range of parameters, $`I_0/\omega _0^2<1`$ and $`\omega \omega _0`$, corresponds to $`z\omega _s/\omega 1`$ for $`\omega \omega _0`$, what are quite typical values of the field strengths and the frequencies in the modern experiments on harmonics generation and detection of THz radiation in SSLs . Finally, it is worth to notice that from the viewpoint of analogy with the Josephson junctions, a described generation of harmonics of an ac electric field in SSL is completely equivalent to the standard ac Josephson effect , i.e. to a generation of current’s harmonics under the action of a sinosuidal voltage.

In this paper we are mainly focused on nonlinear dynamical effects arising within the pendulum limit. With the reference to the pendulum representation of balance equations, we can make the important conclusion: The thresholds of instabilities in an ac-driven semiconductor superlattices, including the transition to chaos, are same as in the RSJ model of an ac-driven Josephson junction providing $`\omega \gamma `$.

Before we proceed with a consideration of instabilities and chaos, we should note that the SSL balance equations (3)-(7) can be also simplified in the opposite limiting case of frequent collisions, $`\omega \gamma `$. Really, assuming $`|\dot{\theta }|\omega \gamma `$ we see from Eq. (4) that $`A(t)`$ follows $`\theta (t)`$ adiabatically. Therefore Eqs. (3), (5), (6), (7) take the form

$$\gamma \dot{\psi }+\omega _0^2\mathrm{sin}\psi =2I_0\mathrm{sin}\omega t$$
(11)
with $`u=\gamma \mathrm{tan}(\psi /2)`$. Such kind of the first order equation is well-known in the theory of Josephson junctions . Importantly, it also arises in the description of an interaction of electromagnetic radiation with the lateral SSLs . Chaos is impossible in the first order differential equation . However, the overdamped pendulum (11) driven by a two-frequency quasiperiodic force still can demonstrate very complex dynamics known as a strange nonchaotic attractor .

Now we turn to the consideration of conditions for a transition to chaos and an influence of temperature effects using the pendulum representation of balance equations (9), (10).

## 3 Periodicity and chaos

The theorem of Levi  states that the stationary dynamics of a periodically forced pendulum is never chaotic in the overdamped limit, which in our case reads

$$C\frac{\gamma }{|w_{\mathrm{eq}}|^{1/2}\omega _0}>2.$$
(12)
The numerical simulations of the driven damped pendulum performed for the wide ranges of $`I_0`$ and $`\omega `$ reveal the absence of chaos already for $`C1`$ .

For a narrow miniband at room temperature, the thermal equilibrium energy becomes close to the center of miniband, i.e. $`|w_{\mathrm{eq}}|`$ is small (see Eq. (2)). Therefore, the criterion for absence of chaos (12) can be easier satisfied in the SSLs with narrow miniband in a comparison to the case of wide miniband SSLs. In particular, for $`k_BT\mathrm{\Delta }`$ we have $`w_{\mathrm{eq}}\mathrm{\Delta }/2k_BT`$ and therefore $`C\gamma /\mathrm{\Delta }`$ (really, $`\omega _0m_0^{1/2}\mathrm{\Delta }^{1/2}`$ providing $`|w_{\mathrm{eq}}|^{1/2}\omega _0\mathrm{\Delta }`$). Thus, with a decrease of miniband width $`\mathrm{\Delta }`$, the value of $`C\gamma /\mathrm{\Delta }`$ increases making the condition (12) to be more robust.

As an example we consider the narrow miniband SSL with $`\mathrm{\Delta }=3`$ meV, $`a=22`$ nm, $`N=2\times 10^{15}`$ cm<sup>-3</sup>  and $`\gamma ^10.5`$ ps , the miniband plasma frequency is $`\omega _{\mathrm{pl}}=2.1\times 10^{12}`$ rad/sec and $`\gamma /\omega _{\mathrm{pl}}1`$. For the temperature $`T=50`$ K, $`w_{\mathrm{eq}}0.33`$ and the factor $`|w_{\mathrm{eq}}^{1/2}|`$ is 1.7. For the higher temperature, $`T=200`$ K, $`|w_{\mathrm{eq}}^{1/2}|3.5`$ and chaos is always impossible. However, chaos is still possible for a wider miniband and a longer relaxation time (see corresponding estimates in Refs. ).

The criterion of absence of strange attractor, Eq. (12), has been obtained within the pendulum limit and thus is formally valid only for $`\omega \gamma `$. However, summing up the results of numerical simulations of the superlattice balance equations, we speculate that strange attractor does not exist at any frequency $`\omega `$, if the inequality (12) is satisfied.

To find an analytic criterion of a transition to chaos in the pendulum limit, we apply the method of Melnikov . This method has been used by several authors for a determination of the conditions of transition to chaos in the periodically driven pendulum (or for the RSJ model) . In our case, the Melnikov’s criterion for a transition to chaos is

$$I_0>I_0^{\mathrm{cr}}=\frac{2\gamma \omega _0|w_{\mathrm{eq}}|^{1/2}}{\pi }\mathrm{cosh}\left(\frac{\pi }{2}\frac{\omega }{|w_{\mathrm{eq}}|^{1/2}\omega _0}\right)$$
(13)
(cf. Eq. (60) in Ref. ). This formula is derived with the assumption that conditions

$$\frac{\gamma }{|w_{\mathrm{eq}}|^{1/2}\omega _0}1,\frac{2I_0}{|w_{\mathrm{eq}}|\omega _0^2}1$$

are satisfied. For $`\omega <\omega _0`$ Eq. (13) takes the form

$$I_0^{\mathrm{cr}}=\frac{2}{\pi }\gamma \omega _0|w_{\mathrm{eq}}|^{1/2},$$
(14)
In the case of SSL driven by the ac electric field $`E=E_0\mathrm{cos}\omega t`$, one can get $`I_0=\omega _s\omega `$ with $`\omega _s=eE_0a/\mathrm{}`$ . Then, the formula (14) gives for the boundary of chaotic region in the $`\omega _s`$-$`\omega `$ plane in the form

$$\omega _s\omega =\frac{2}{\pi }\gamma \omega _0|w_{\mathrm{eq}}|^{1/2}.$$
(15)
Earlier, the low-frequency boundary of chaos in the simple form $`\omega _s\omega =\text{const}`$ has been found numerically in the works , but has been unexplained until now. Eq. (15) can qualitatively explain this dependence.

## 4 Phase-locking in pendulum and dc voltage generation in superlattice

Phase-locking is another nonlinear dynamic phenomenon which an ac-driven pendulum can demonstrate along with chaos . In the physics of Josephson junctions, phase-locking at zero bias is known as the inverse ac Josephson effect . Here we are discussing the influence of phase-locking in the pendulum on the physical properties of SSL. In the driven pendulum phase-locking means that

$$\dot{\psi }=\frac{n}{l}\omega ,$$
(16)
where $`n`$ and $`l`$ are integer numbers and $`\mathrm{}`$ stands for the time averaging over the period of ac drive $`2\pi /\omega `$ . Majority of stable phase-locked states are integer ($`n0`$, $`l=1`$); however, a fractional phase-locking ($`l>1`$) can also exist. In the case of Josephson junction, the voltage $`U`$ across the junction is proportional to $`\dot{\psi }`$; therefore phase-locking results in a generation of quantized dc voltage $`U(n/l)\omega `$ .

For the superlattice problem, the voltage is a function of both the velocity and the phase of the pendulum (see Eq. (10)). Therefore, phase-locking in the pendulum (16) determines a generation of dc voltage in the SSL as

$`u`$ $`=`$ $`{\displaystyle \frac{n}{l}}{\displaystyle \frac{\omega }{2}}\left(1+{\displaystyle \frac{\gamma }{\omega _0}}R_{nl}\right),`$ (17)
$`R_{nl}`$ $`=`$ $`{\displaystyle \frac{2\omega _0}{\omega }}{\displaystyle \frac{l}{n}}\mathrm{tan}\left({\displaystyle \frac{\psi }{2}}\right).`$ (18)
For a low damping $`\gamma /\omega _{\mathrm{pl}}1`$, we immedeatly get from (17) the dc voltage $`u(n/2l)\omega `$. Such kind of almost half-integer dc voltage states in a pure ac-driven SSL have been recently observed in the numerical simulations of balance equations .

The formula (17) tells us nothing about a stability of dc voltage states in SSL that correspond to the phase-locked states of the pendulum with different $`n`$ and $`l`$. However, numerical simulations demonstrate that the dc voltage states, which are close to the integer states (i.e. to the states with $`l=1`$ and $`n`$ being even integer), are more typical then fractional states .

What is important, the dc voltage per superlattice period is not exactly $`(n/2l)\omega `$ even for a weak damping. Moreover, a weak dependences of $`u`$ on the ac amplitude  and the frequency  have been found numerically at $`\gamma /\omega _{\mathrm{pl}}1`$. The appearance of such kind of dependences within the pendulum limit can be understood, if one takes into an account that $`u`$ is a function of both the time-average velocity and of the time-average coordinate of pendulum (see Eq. (17)). Phase-locked rotational state of the pendulum can be represented in the form

$$\psi (t)=\psi _0+n\omega +\underset{p=1}{\overset{\mathrm{}}{}}\alpha _p\mathrm{sin}(p\omega t+\mu _p),$$
(19)
where $`\psi _0`$, $`\alpha _p`$ and $`\mu _p`$ are constants and we restrict our consideration to the leading phase-locked states, $`\dot{\psi }=n\omega `$. The amplitudes $`\alpha _p`$ and the phases $`\mu _p`$ in Eq. (19) are functions of the main physical parameters of our problem, i.e. $`I_0`$, $`\omega `$ and $`\gamma `$. Substituting the expression (19) in Eq. (18) one can see that the dependence of dc voltage $`u`$ on the ac amplitude and frequency arise in Eq. 17 via the term $`R_{nl}`$.

It is important to know how this “dissipative correction” to the quantized dc voltage, $`un\omega /2=(\gamma /\omega _0)R_n`$, scales with an increase of dissipation. To estimate the dependence of $`R_n`$ on $`\gamma /\omega _0`$ for several lowest phase-locked states of the pendulum, we consider the limit of fast rotations $`\alpha 1`$. In this case formula (19) takes the form

$$\psi (t)=\psi _0+n\omega +\alpha \mathrm{sin}(\omega t+\mu ).$$
(20)
Substituting Eq. (20) in Eq. (9) and equating zero harmonics, we have the formula (cf. )

$$\omega _0^2J_n(\alpha )\mathrm{sin}(n\mu )=n\gamma \omega ,$$
(21)
which gives us an opportunity to understand the characteristic dependences of $`\alpha `$ and $`\psi `$ on $`\gamma `$ and $`\omega `$. Next, substituting Eq. (20) in Eq. (18) and calculating integrals in the limit $`\alpha 1`$, we obtain (see Appendix C)

$`\mathrm{tan}\left(\psi /2\right)_{n=1}`$ $`=`$ $`\alpha \mathrm{sin}(\mu )+𝒪(\alpha ^2),`$ (22)
$`\mathrm{tan}\left(\psi /2\right)_{n=2}`$ $`=`$ $`0.25\alpha ^2\mathrm{sin}(2\mu )+𝒪(\alpha ^4).`$
Now combining Eqs. (18), (21), (22) and using the asymptotics of the Bessel functions $`J_n(z)z^n/(2^nn!)`$ ($`z1`$), we find

$$R_n\left(1\right)^n\frac{4\gamma }{\omega _0},\text{(n=1,2)}.$$
(23)
Thus, for a weak dissipation ($`\gamma /\omega _{\mathrm{pl}}1`$), the dc voltage spontaneously generated in a pure ac-driven SSL is

$$u=n\omega /2+𝒪(\gamma ^2/\omega _{\mathrm{pl}}^2)$$
(24)
and it has a quadratic dependence on $`\gamma `$.

With a further increase of $`\gamma `$, the dissipative term in Eq. (17) becomes more important. Moreover, for a large dissipation the phase-locked rotational states (19), (20) do not exist in the pendulum anymore . However, the unquantized dc voltage still can be generated in SSL due to different mechanisms, as has been demonstrated numerically in . The discussion of relationship between a pendulum dynamics and a dc voltage generation in a superlattice for an arbitrary damping will be presented elsewhere.

## 5 Discussion and conclusion

It is instructive to compare the pendulum equation (9), derived from the superlattice balance equations in the limit of rare collisions ($`\gamma \omega `$) and corresponding motion equation that can be obtained employing the Newton law. We introduce simplest form of phenomenological dissipative term to the motion equation for electron’s momentum as

$$\dot{p}+\gamma p=eE(t).$$

Further, combining this equation and the standard definition of electron’s velocity $`V=\epsilon (p)/p`$, the tight-binding dispersion relation, and the Maxwell equation (third equation in the set (1)), we get the pendulum equation in the form

$$\ddot{\varphi }+\gamma \dot{\varphi }+\omega _0^2\mathrm{sin}\varphi =I_0\mathrm{sin}\omega t,$$

where $`\varphi =pa/\mathrm{}`$, $`|w_{\mathrm{eq}}|1`$ for simplicity. In this approach the scaled voltage is $`u=\dot{\varphi }`$. Certainly, this equation is similar to Eq. (9). The difference arises mainly in the dependence of voltage on coordinate and momentum of the pendulum (cf. Eq. (10)). In particular, the voltage is dependent not only on the momentum of pendulum but also on its coordinate in the Boltzmann equation approach.

The crude approach based on the Newton law uses a trajectory of an individual electron. Therefore, this approach can not take into an account the changes in a distribution function of electrons induced by an electric field. In contrast, the approach based on the Boltzmann equation does include such kind of corrections in the term proportional to $`\gamma `$. In the limit of infrequent collisions, $`\gamma \omega `$, this term can be often ignored because $`\dot{\theta }\omega \gamma `$. However, the dissipative contribution to voltage can become important if one considers a dc voltage generation, $`u0`$. As we have demonstrated in section 4, this dissipative correction can be important in the consideration of phase locking, $`\dot{\theta }=\omega `$. The contribution of dissipative term to the dc voltage, $`u`$, can be even more important in the case of lateral superlattices, for which $`\dot{\theta }=0`$ .

In summary, we show that pendulum equation can well describe the dynamics of electrons in a miniband of ac-driven semiconductor superlattice in the limit of weak scattering. We demonstrate that several earlier numerical findings concerning chaos and spontaneous dc voltage generation in the superlattices can be explained within such approach. We establish a link between the theory of diffusive transport in ac-driven superlattices and the theory rf-driven Josephson junctions. We point out that dissipative terms arising in the dependence of voltage across superlattice on the momentum and the coordinate of the pendulum should play an important role in many situations.

We believe that described in this work simple and instructive theory can be useful in an explanation of the results of experiments  as well as that it can stimulate other experiments on the THz-field induced nonlinear phenomena in superlattices.

## Acknowledgements

We are grateful to David Campbell, Ethan Cannon, Alexander Zharov, Pekka Pietiläinen for collaboration. We thank Laurence Eaves, Martin Holthaus, Anatoly Ignatov, Karl Renk, Mikko Saarela, Erkki Thuneberg for discussions. This research was partially supported by the Academy of Finland and the Royal Society.

## Appendix A Appendix.

In this appendix we will derive the equations (3), (4), (5), (6), (7). To begin with we introduce a complex variable $`Z=w+iv`$ and represent the superlattice balance equations (1) in the form

$$\dot{Z}=iuZ\gamma Z+\gamma w_{\mathrm{eq}}.$$
(25)
Substituting the parametrization $`Z=A\mathrm{exp}(i\theta )`$ (i.e. $`w=A\mathrm{cos}\theta `$ and $`v=A\mathrm{sin}\theta `$) in (25) and equating real and imaginary parts, we get

$`\dot{A}`$ $`=`$ $`\gamma A\gamma w_{\mathrm{eq}}\mathrm{cos}\theta ,`$ (26)
$`\dot{\theta }`$ $`=`$ $`u+\gamma w_{\mathrm{eq}}A^1\mathrm{sin}\theta .`$ (27)
From the last equation we have

$$\ddot{\theta }=\dot{u}+\gamma w_{\mathrm{eq}}\frac{d}{dt}\left[A^1\mathrm{sin}\theta \right].$$
(28)
Taking into account that

$`{\displaystyle \frac{d}{dt}}[A^1\mathrm{sin}\theta ]=A^2\dot{A}\mathrm{sin}\theta +A^1\dot{\theta }\mathrm{cos}\theta `$ (29)
$`=\gamma A^1\mathrm{sin}\theta +\gamma w_{\mathrm{eq}}A^2\mathrm{sin}\theta \mathrm{cos}\theta +\dot{\theta }A^1\mathrm{cos}\theta `$
and using the third equation of set (1), we immediately get from Eq. (28) the Eqs. (3), (5), (6). Finally, it is easy to see that Eqs. (26) and (27) are same as Eqs. (4) and (7), correspondently.

## Appendix B Appendix.

In this appendix we will show that well-known result of Tsu-Esaki-Ignatov-Romanov from the theory of ac-driven SSLs  can be reproduced within the pendulum representation of the SSL balance equations (9), (10).

Solving the Boltzmann equation for the tight-binding lattice driven by the electric field $`E=E_0\mathrm{cos}\omega t`$ with $`\omega \omega _0`$, Ignatov and Romanov  have found for $`\gamma \omega `$ the following time-dependence of the averaged electron’s velocity

$$v(t)=\underset{k=1}{\overset{\mathrm{}}{}}a_{2k1}\mathrm{sin}(2k1)\omega t,a_{2k1}=2J_0(z)J_{2k1}(z),$$
(30)
where $`z=\omega _s/\omega =(eE_0a)/\mathrm{}\omega `$. Same result has been obtained earlier by Tsu and Esaki for the limiting case $`z1`$ .

Now turn to the pendulum equation. Hence we have the ac field of a very high frequency, we can neglect the last term ($`\gamma `$) in Eq. (10) in comparison with the first term ($`\dot{\theta }\omega `$). Further, we can employ an analogy with the RSJ model of Josephson junctions. It is known for the RSJ model that the sinusoidal voltage is the solution of pendulum equation if $`\omega \omega _0`$ and $`\gamma \omega `$ . In our case it means that $`\theta (t)=z\mathrm{sin}\omega t`$ is the solution of Eq. (9). Further, substituting this time dependence of $`\theta (t)`$ in the expression for averaged electron’s velocity $`v=\mathrm{sin}\theta \mathrm{cos}\theta `$, we can determine the current generated in the SSL under the action of the ac field. (Of course, it is the analog of the ac Josephson current ). Using the Bessel expansions

$$\mathrm{cos}(z\mathrm{sin}t)=J_0(z)+2\underset{k=1}{\overset{\mathrm{}}{}}J_{2k}(z)\mathrm{cos}2kt,$$

$$\mathrm{sin}(z\mathrm{sin}t=2\underset{k=1}{\overset{\mathrm{}}{}}J_{2k1}(z)\mathrm{cos}(2k1)t,$$

we find the averaged electron’s velocity in the form

$`v(t)=2J_0(z){\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}J_{2k1}(z)\mathrm{sin}(2k1)\omega t`$ (31)
$`+2{\displaystyle \underset{k_1,k_21}{}}J_{2k_1}(z)J_{2k_21}(z)\times `$
$`\left[\mathrm{sin}\left(2(k_2k_1)1\right)\omega t+\mathrm{sin}\left(2(k_2+k_1)1\right)\omega t\right].`$
Supposing that

$$\omega _s<\frac{\omega _0^2}{\omega },$$
(32)
we have

$$z\frac{\omega _s}{\omega }<\frac{\omega _0^2}{\omega ^2}1$$

because $`\omega /\omega _01`$. Now it is easy to see that for $`z1`$ Eq. (31) is same as Eq. (30).

Thus, we have showed that for $`\gamma \omega `$ and, additionally, if the condition (32) is satisfied, the pendulum equation can reproduce the result of Tsu-Esaki-Ignatov-Romanov. Finally, combining the definition $`I_0=\omega _s\omega `$ and Eq. (32), we see that the inequality (32) just gives another condition of the validity of the pendulum approach: $`I_0/\omega _0^2<1`$.

## Appendix C Appendix.

In this appendix we will derive the equations (22). We need to calculate the integrals

$$\mathrm{tan}(\psi /2)\frac{1}{2\pi }_0^{2\pi }𝑑x\mathrm{tan}\left(\frac{nx}{2}+\frac{\alpha }{2}\mathrm{sin}(x+\mu )\right),$$
(33)
where $`n=1,2`$. In the limit $`\alpha 1`$ we can use the expansions

$`\mathrm{tan}\left({\displaystyle \frac{x}{2}}+{\displaystyle \frac{\alpha }{2}}\mathrm{sin}(x+\mu )\right)=\mathrm{tan}\left({\displaystyle \frac{x}{2}}\right)+`$
$`{\displaystyle \frac{\mathrm{sin}(x+\mu )}{\mathrm{cos}(x)+1}}\alpha 𝒪\left(\alpha ^2\right),`$ (34)
$`\mathrm{tan}\left(x+{\displaystyle \frac{\alpha }{2}}\mathrm{sin}(x+\mu )\right)=\mathrm{tan}(x)+{\displaystyle \frac{\alpha }{2}}{\displaystyle \frac{\mathrm{sin}(x+\mu )}{\mathrm{cos}^2(x)}}`$
$`{\displaystyle \frac{\alpha ^2}{4}}{\displaystyle \frac{\left(1+\mathrm{cos}^2(x+\mu )\right)\mathrm{sin}(x)}{\mathrm{cos}^3x}}+`$
$`{\displaystyle \frac{\alpha ^3}{24}}{\displaystyle \frac{\left(1+\mathrm{cos}^2(x+\mu )\right)\mathrm{sin}(x+\mu )\left(2\mathrm{cos}^2x3\right)}{\mathrm{cos}^4x}}`$
$`+𝒪\left(\alpha ^4\right).`$ (35)
After integration of the expression (C) and (35) we obtain Eqs. (22).


format: text

text:
# Radiative Heat Conduction and the Magnetorotational Instability

## 1 Preliminaries

Understanding the magnetorotational instability is key to the development of realistic models of accretion. As a physical process, the MRI justifies two sine qua non ingredients of accretion disk theory: entropy generation from the differential shear flow and turbulent angular momentum transport. Consequently, clarifying the dynamics that leads to MRI initiated turbulence in a few relevant astrophysical regimes is paramount in uncovering how the accretion process takes place.

In the context of accretion onto compact objects, the stress associated with a radiative particle species often dominates the dynamics at high temperatures but our understanding of its effect on MHD processes is incomplete at best. Agol & Krolik (1998) pointed out that compressive radiation-MHD wave modes will be strongly suppressed by diffusive loss of pressure support in a range of wave numbers: $`c_{\mathrm{ph}}/c<\stackrel{~}{k}^1<1`$, where $`\stackrel{~}{k}\mathrm{}_{\mathrm{mfp}}k`$ represents the wave number normalized to the mean-free-path and $`c_{\mathrm{ph}}`$ is the phase speed of compressive waves in the single fluid view. Such analysis is applicable to accretion disks at large wave numbers where the effect of inertial forces and stratification are minor. For the MRI, however, both inertial and stratification effects play key roles in dictating the dynamics and polarization of the modes.

Blaes & Socrates (2001) have reported on a radiation-MHD dispersion relation for axisymmetric modes that comprises the MRI under rather general circumstances. Our analysis is different from theirs mainly in the consideration of a purely toroidal field where the compressibility of the non-axisymmetric modes leads to the largest impact from radiative diffusive damping. Furthermore, aiming to understand ultra-relativistic flows very near black holes, we emphasize toroidal mode dynamics because these modes may become predominant when shear effects begin to overwhelm other inertial accelerations leading to fastest growing modes at large azimuthal scales (Araya-Góchez 2002). The scale and degree of intermittency in the turbulent flow bear a direct impact on gravitational wave emission models from hyper-accreting black holes (Araya-Góchez 2003). For a relativistically hot accretion disk with comparable vertical pressure support from diffusive and non-diffusive species (see §2), we solve for an approximate dispersion relation involving Lagrangian displacements in the horizontal regime ($`\xi _\theta /\xi _rk_r/k_\theta \mathrm{}`$). This simplification affords a lowering of the rank in the dispersion relation to second order (in $`\omega ^2`$) which, in turn, enables a relatively straightforward approximate treatment of the effects of radiative heat conduction (in contrast with Blaes and Socrates’ more complex but exact, relation for axisymmetric modes).

## 2 Radiative Viscosity vs Diffusion

The half optical depth to scattering of a standard radiation-pressure–dominated $`\alpha `$-disk is related to the vertical scale height, $``$, and to the angular rotation profile, $`\mathrm{\Omega }(r)`$, through

$$\tau _{\mathrm{disk}}=\frac{c}{2\alpha \widehat{\mathrm{A}}\mathrm{\Omega }}$$
(1)
where Oort’s A parameter is normalized $`\widehat{\mathrm{A}}\frac{1}{2}\mathrm{d}_{\mathrm{ln}r}\mathrm{ln}\mathrm{\Omega }`$ ($`3/4`$ for a non-relativistic, Keplerian flow). This relation ensues from balancing the local heating rate, $`\dot{M}\mathrm{\Omega }^2\widehat{\mathrm{A}}`$, with the local cooling rate from radiative diffusion, $`_\tau p`$; while eliminating the accretion rate in favor of the vertically integrated component of the viscous stress responsible for angular momentum transport $`\alpha p`$. Eq 1 is thus explicitly sensitive only to the local nature of the cooling (although it is implicitly subject to a suitable vertical gradient of heat deposition; see, e.g., Krolik 1999).

In the neutron-rich, neutrino cooled gas associated with a hyper-accreting black hole (see, e.g. Popham et al. 1999), the pressure contributions from radiation, pairs and neutrinos, all have the same temperature dependence, $`a/6T^4`$, with relative contributions varying only by internal degrees of freedom times particle statistics factors: 2 $`\times `$ 1, 4 $`\times `$ 7/8, and 6 $`\times `$ 7/8, respectively (with only one helicity state per neutrino). Thus, the neutrino pressure is never greater than $`p_{\mathrm{thermal}}`$ where $`p_{\mathrm{thermal}}=11/12aT^4`$ includes the pressure from photons and pairs. Yet, only the former species has radiative diffusion properties on a dynamical time scale. Thus, in this setting the neutrinos are the lone diffusive species and their stress is identified with $`p_{\mathrm{rad}}`$ below (e.g., radiative = diffusive species throughout this paper).

Partial pressure support by the non-diffusive components (baryons, toroidal field and the thermal component) does not modify Eq 1 as long as the contribution to the disk’s flare from this extra pressure is inversely proportional to the mass density decrease from the same, e.g., 1D expansion. Consequently, for a neutrino cooled $`\alpha `$-disk, a replacement in the source of opacity by neutrino scattering in the non-advective accretion regime: $`0.1\mathrm{M}_{}\mathrm{sec}^1\dot{\mathrm{M}}1\mathrm{M}_{}\mathrm{sec}^1`$ (Di Matteo et al. 2002), is all that is needed for Eq 1 to hold.

On the other hand, the “parallel” size the of fastest growing MRI eddies is $`l_{\mathrm{eddy}}=k_{\mathrm{MRI}}^1v_{_{\mathrm{Alf}}}/\mathrm{\Omega }`$. For magnetic angular momentum transport such that $`\alpha (v_{_{\mathrm{Alf}}}/c_{_{}})^2`$, the isotropic diffusion time through these eddies, $`t_{\mathrm{diff}}\mathrm{\Omega }^1`$, turns out to be similar to the timescale for fastest MRI modes to develop $`t_{\mathrm{MRI}}\mathrm{A}_{\mathrm{Oort}}^1`$ (see §5). Moreover, the optical depth through these eddies is (with $`c_{_{}}^2p_{\mathrm{tot}}/\rho `$)

$$\tau _{\mathrm{eddy}}\frac{v_{_{\mathrm{Alf}}}}{c_{_{}}}\times \tau _{\mathrm{disk}}\frac{1}{2\widehat{\mathrm{A}}\sqrt{\alpha }}\frac{c}{c_{_{}}}1.$$
(2)
Since Compton/neutrino drag go as $`𝐟_{\mathrm{drag}}\stackrel{~}{k}^2\delta 𝜷`$ and $`\stackrel{~}{k}^2𝒪(\alpha c_{_{}}^2/c^2)`$, these estimates indicate that the effect of photon/neutrino viscosity is negligible for the optically thick eddies of a radiatively cooled disk.

## 3 A Lagrangian formulation of the MRI in soft media

We take a semi-local approach by computing instantaneous (co-moving) growth rates in the limit $`k_r/k_\theta 0`$ and referring to these as the fastest growing modes on the implicit understanding of their transient nature (since $`k_r(t)=k_{0r}[\mathrm{d}_{\mathrm{ln}r}\mathrm{\Omega }]k_\phi t`$). In terms of a local, co-moving observer, azimuthal wave numbers are no longer discrete and consequently, neither are the co-moving frequencies (Ogilvie & Pringle 1996). Further, we discard field curvature and radial gradient terms. A simple meridional stratification profile sets the physical scale-length of the problem: d$`{}_{z}{}^{}\mathrm{ln}\rho =1/`$, with gas, radiative and magnetic pressures tracking the unperturbed density profile.

The first order relation between Lagrangian, $`\mathrm{\Delta }`$, and Eulerian, $`\delta `$, variations (see, e.g., Chandrasekhar & Lebovitz 1964, Lynden-Bell & Ostriker 1967) is

$$\mathrm{\Delta }=\delta +𝝃\mathbf{}.$$
(3)
The Euler velocity perturbation, $`𝐯\delta 𝐕`$, is thus related to the Lagrangian displacement, $`𝝃`$, through

$$𝐯=\{_t+𝐕\mathbf{}\}𝝃(𝝃\mathbf{})𝐕i\sigma 𝝃\xi ^r\mathrm{\Omega }_{,r}\mathrm{𝟏}_\phi .$$
(4)
The algebraic relation follows from the assumption of differential rotation and from writing exp$`i(\omega t+m\phi +k_zz)`$ dependencies for $`𝝃`$, while $`\sigma \omega +m\mathrm{\Omega }`$ denotes the co-moving frequency of the perturbations.

These geometrical equations have their more traditional equivalents in the so-called shearing sheet approximation where a co-moving, “locally Cartesian frame” $`(\widehat{r},\widehat{\phi },\widehat{\theta })(x,y,z)`$, is used along with the linearized shear velocity field, $`𝐕(x)=[\mathrm{d}_{\mathrm{ln}r}\mathrm{\Omega }]x\mathrm{𝟏}_y`$, to treat the problem locally while introducing the coriolis terms by hand. (The reader is invited to examine a companion paper: Araya-Góchez (2002), where the inertial terms are induced in a fully covariant, geometric way).

The linearized equations of motion for the displacement vector correspond to the Hill equations (Chandrasekhar 1961, Balbus & Hawley 1992)

$$\ddot{\xi }_r2\mathrm{\Omega }\dot{\xi }_\phi =4\mathrm{A}_{\mathrm{Oort}}\mathrm{\Omega }\xi _r+\frac{f_r}{\rho },\ddot{\xi }_\phi +2\mathrm{\Omega }\dot{\xi }_r=\frac{f_y}{\rho },\mathrm{and}\ddot{\xi }_\theta =\frac{f_\theta }{\rho },$$
(5)
where each over-dot stands for a Lagrangian time derivative and f, for the sum of body forces.

The Lagrangian perturbation of mass density and the Eulerian perturbation of the field follow straightforwardly from Eq , and from mass and magnetic flux conservations (under ideal MHD conditions, notwithstanding the peculiar nature of the fluid under study)

$$\frac{\mathrm{\Delta }\rho }{\rho }=\mathbf{}𝝃,\delta 𝐁=\mathbf{}\times (𝝃\times 𝐁).$$
(6)
Using the latter relation, the Lorentz force in the co-moving frame is readily laid out

$`\delta {\displaystyle \frac{1}{4\pi \rho }}𝐉\times 𝐁`$ $`=`$ $`v_{_{\mathrm{Alf}}}^2\times \left\{\mathbf{}(\mathbf{}𝝃)+_{\widehat{𝐁}}^2𝝃_{\widehat{𝐁}}\mathbf{}(\mathrm{𝟏}_{\widehat{𝐁}}𝝃)\mathrm{𝟏}_{\widehat{𝐁}}_{\widehat{𝐁}}(\mathbf{}𝝃)\right\}`$ (9)
$`\stackrel{i𝐤}{}v_{_{\mathrm{Alf}}}^2\times \{(k_i\xi _ik_{\widehat{𝐁}}\xi _{\widehat{𝐁}})𝐤+k_{\widehat{𝐁}}^2𝝃\mathrm{𝟏}_{\widehat{𝐁}}k_{\widehat{𝐁}}k_i\xi _i.\}`$
where $`\mathrm{𝟏}_{\widehat{𝐁}}`$ is a unit vector in the direction of the unperturbed field and where the scalar operator $`_{\widehat{𝐁}}\mathrm{𝟏}_{\widehat{𝐁}}\mathbf{}`$. Note that the term $`(k_i\xi _ik_{\widehat{𝐁}}\xi _{\widehat{𝐁}})k_{}\xi _{}`$ may be interpreted as a restoring force due to the compression of field lines (Foglizzo & Tagger 1995).

Next, we look in detail at the non-magnetic stress terms.

In the single fluid view, the Lagrangian variation of the specific pressure gradient, $`\rho ^1\mathbf{}p_{\mathrm{r}+\mathrm{g}}`$, contains two terms (see, e.g., Lynden-Bell & Ostriker 1967): one $`\mathrm{\Delta }\rho ^1`$, and another one $`\mathrm{\Delta }\mathbf{}p_{\mathrm{r}+\mathrm{g}}`$. In terms of the displacement vector, the first term is proportional to the equilibrium value of $`\mathbf{}p_{\mathrm{r}+\mathrm{g}}`$ which is negligible in the local treatment ($``$ radial gradient when $`\xi _\theta /\xi _r\mathrm{}`$). For the same reason, the Eulerian and Lagrangian variations of the specific pressure gradient are nearly identical.

The “thermodynamic” pressure term is thus given by

$$\delta \left(\frac{1}{\rho }\mathbf{}p_{\mathrm{r}+\mathrm{g}}\right)=\mathrm{\Gamma }\frac{p_{\mathrm{r}+\mathrm{g}}}{\rho }\mathbf{}(\mathbf{}𝝃)\stackrel{i𝐤}{}c_s^2𝐤k_i\xi _i,$$
(10)
where, for a heterogeneous fluid, $`\mathrm{\Gamma }(d_{[\mathrm{ln}\rho ]}\mathrm{ln}p)`$ represents a generalized adiabatic index (see, e.g., Chandrasekhar 1939, Mihalas & Mihalas 1984).

Assembling f from Eqs \[9 & 10\], and plugging this form into Eqs  yields

$$\ddot{𝝃}+2𝛀\times \dot{𝝃}=\{c_s^2𝐤+v_{_{\mathrm{Alf}}}^2(𝐤k_{\widehat{𝐁}}\mathrm{𝟏}_{\widehat{𝐁}})\}(𝐤𝝃)+v_{_{\mathrm{Alf}}}^2(k_{\widehat{𝐁}}^2𝝃k_{\widehat{𝐁}}\xi _{\widehat{𝐁}}𝐤)4\mathrm{A}\mathrm{\Omega }\xi _r\mathrm{𝟏}_r.$$
(11)
This equation agrees with the matrix de-composition of Foglizzo & Tagger (1995) for $`\mathrm{\Gamma }=1`$. Moreover, these authors also pointed out that because of vertical stratification and rotation, the polarization of the slow MHD (e.g., MRI) modes obey

$$\frac{\xi _\theta }{\xi _r}𝒪\left(\frac{k_r}{k_\theta }\right)$$
(12)
when $`k_r/k_\theta \mathrm{}`$. Eqs \[11 & 12\] impose an anisotropy constraint on the components of the Lagrangian displacement

$$k_{}\xi _{}=\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }+2\mathrm{\Theta }}k_{}\xi _{}\mathrm{\Lambda }k_{}\xi _{},$$
(13)
where $`\mathrm{\Theta }p__𝐁/p_{\mathrm{r}+\mathrm{g}}`$.

In this horizontal regime, with $`k_z\xi _z`$ finite, the dispersion relation derived from Eq  and from the above constraint reads

$$\widehat{\omega }^4\{(\mathrm{\Lambda }+1)\widehat{q}_{_{\widehat{𝐁}}}^2+\widehat{\chi }^2\}\widehat{\omega }^2+\mathrm{\Lambda }\widehat{q}_{_{\widehat{𝐁}}}^2\{\widehat{q}_{_{\widehat{𝐁}}}^2+4\widehat{\mathrm{A}}\}=\mathrm{},$$
(14)
where all frequencies are normalized to the local rotation rate, $`\widehat{\chi }^24(1+\widehat{\mathrm{A}})`$ is the squared of the epicyclic frequency, and $`\widehat{q}_{_{\widehat{𝐁}}}(𝐤𝐯_{\mathrm{Alf}})/\mathrm{\Omega }`$ is a frequency related to the component of the wave vector along the field (in velocity units). This expression matches the form given by T. Foglizzo (1995).

In terms of the function, $`\mathrm{\Lambda }(\mathrm{\Gamma },\mathrm{\Theta })`$, we find the fastest growing wave numbers to be given by

$$\widehat{q}_{_{\widehat{𝐁}}}^2=2\widehat{\mathrm{A}}+(\frac{1+\mathrm{\Lambda }}{2\mathrm{\Lambda }})\times \left\{\frac{2\mathrm{\Lambda }\widehat{\mathrm{A}}^2}{𝒟}\right\},\mathrm{where}𝒟1+(\frac{1\mathrm{\Lambda }}{2})\widehat{\mathrm{A}}+\sqrt{1+(1\mathrm{\Lambda })\widehat{\mathrm{A}}},$$
(15)
and where the expression in the curly brackets corresponds to the growth rate of the modes.

Notably, the compressibility of non-axisymmetric, MRI modes is imprint on deviations of $`\mathrm{\Lambda }`$ from unity

$$\frac{\mathrm{\Delta }\rho }{\rho }=(1\mathrm{\Lambda })(ik_{})\xi _{},$$
(16)
e.g., the degree of compression gets stronger with field strength and, naturally, with a softer equation of state.

## 4 A Quasi-Adiabatic Index for Wave Phenomena

In a heterogeneous fluid composed of a neutral baryonic plasma (with ideal gas properties) plus a neutrino and/or photon component, ideal MHD dictates that the field lines remain frozen to the charged massive component while pressure perturbations associated with such material gas must track density perturbations on all scales. On the other hand, collisional coupling between the gas and the diffusive species translates into a scale dependent containment of the pressure perturbation associated with the radiative species. Consequently, acoustic wave modes partially supported by radiative stress will be damped from a non-adiabatic loss of pressure support when the diffusive and wave time scales are comparable (Agol & Krolik 1998). Such phenomenology can be quantitatively implemented by considering the mathematical equivalent to the equation of state for the diffusive component of the fluid.

We start with the co-moving frame, frequency-integrated, radiative transfer equation–corrected to first order in the fluid’s motion, $`𝜷`$, and with isotropic, elastic scattering as the lone source of opacity:

$$\stackrel{~}{}_tI(𝐧)+𝐧\stackrel{~}{\mathbf{}}I(𝐧)=(1+3𝜷𝐧)J2𝜷𝐇(1𝐧𝜷)I(𝐧).$$
(17)
A tilde means normalization to the mean free path, $`\mathrm{}_{\mathrm{mfp}}`$, or to the mean crossing time, $`\mathrm{}_{\mathrm{mfp}}/c`$, while $`I_\nu (𝐧)`$ is the radiative intensity in the direction of n, and $`J`$ and H are the first two moments of the intensity.

Agol and Krolik (1998) computed the first three radiative moment equations of Eq  and, by assuming that multipoles of $`I(𝐧)`$ higher than quadrupole vanish, found an otherwise exact form for the mean intensity perturbation, $`\delta J\frac{1}{\stackrel{~}{\omega }}\stackrel{~}{𝐤}\delta 𝜷`$, in a background with $`J=I,𝐇=\mathrm{}`$ (e.g., such that $`\delta =\mathrm{\Delta }`$). Tacitly conforming with the premise of isotropy, we re-write Eq  of Agol & Krolik (1998) by taking each of the perturbed diagonal components of the radiative field stress tensor, $`\delta \underset{¯}{K}=(1/4\pi )\delta 𝑑\nu 𝑑\mathrm{\Omega }\mathrm{𝐧𝐧}I_\nu (𝐧)`$, to be identical and proportional to the perturbed radiative pressure $`\delta p_{\mathrm{rad}}`$. Thus $`\delta K_{ii}=\frac{1}{3}\delta J`$ and, upon insertion of $`\delta \mathrm{ln}V`$ in lieu of $`\frac{1}{\stackrel{~}{\omega }}\stackrel{~}{𝐤}\delta 𝜷`$, this yields the following relation for quasi-adiabatic perturbations of the radiative component of the fluid

$$\frac{\delta p_{\mathrm{rad}}}{p_{\mathrm{rad}}}=\frac{4}{3}\left\{(1i\stackrel{~}{\omega })+\frac{i}{15}\frac{\stackrel{~}{k}^2}{\stackrel{~}{\omega }}\frac{59i\stackrel{~}{\omega }}{1i\stackrel{~}{\omega }}\right\}^1\frac{\delta V}{V}.$$
(18)
Let us inspect the natural parameterization of this expression. When the non-diffusive pressure is negligible, this index represents a truly adiabatic processes for wave modes only if $`\stackrel{~}{\omega },\stackrel{~}{k}^2`$ & $`\stackrel{~}{k}^2/\stackrel{~}{\omega }1`$. Individually, both $`\stackrel{~}{\omega }=\omega (\mathrm{}_{\mathrm{mfp}}/c)|\widehat{\mathrm{A}}|(c_{_{}}/c)\tau _{\mathrm{disk}}^12\alpha \widehat{\mathrm{A}}^2(c_{_{}}/c)^2`$, and $`\stackrel{~}{k}k\tau _{\mathrm{disk}}^12\sqrt{\alpha }|\widehat{\mathrm{A}}|(c_{_{}}/c)`$, are expected to be small for the scales of interest in the magnetorotational instability problem. However, the ratio $`\stackrel{~}{k}^2/\stackrel{~}{\omega }`$ is not small and this has a non-trivial interpretation:

$$\left\{\frac{k_{}^2}{𝐤^2}\right\}\frac{\stackrel{~}{k}^2}{\stackrel{~}{\omega }}i\frac{\mathrm{}_{\mathrm{mfp}}^2/l_{\mathrm{eddy}}^2}{\mathrm{}_{\mathrm{mfp}}|\omega _{_{\mathrm{MRI}}}|/c}\frac{i}{|\mathrm{A}|}\frac{(c/\tau _{\mathrm{eddy}})}{l_{\mathrm{eddy}}}𝒪\left(2i\right)$$
(19)
which identifies with the ratio of eddy diffusion rate to MRI growth rate. $`\alpha _4\phi \upsilon `$

To zeroth order, from Eq  one thus has

$$\stackrel{~}{\mathrm{\Gamma }}\frac{4}{3}(1+\frac{i}{3}\stackrel{~}{k}^2/\stackrel{~}{\omega })^1$$
(20)
which indicates that acoustic wave perturbations (i.e., with $`\omega `$ real) on the scale of the MRI are strongly damped if radiative pressure is the predominant source of stress. Indeed, Eq  may be interpreted as a form of energy equation characterizing radiative “heat conduction” out of compressive wave modes (Agol & Krolik 1998, compare this with Eq\[51.12\] of Mihalas & Mihalas 1984).

If there is more parity between the diffusive and non-diffusive sources of stress, we find a generalized quasi-adiabatic index by defining a thermodynamic “quasi-volume”, $`\stackrel{~}{V}V^{\stackrel{~}{\eta }}`$, through Eq  as follows

$$\frac{dp_{\mathrm{rad}}}{p_{\mathrm{rad}}}=\frac{4}{3}\stackrel{~}{\eta }\frac{dV}{V}\frac{4}{3}\frac{d\stackrel{~}{V}}{\stackrel{~}{V}},$$
(21)
while computing the diffusive stress contribution to thermodynamic processes in terms of the logarithmic differential of such quasi-volume: $`d\mathrm{ln}\stackrel{~}{V}=\stackrel{~}{\eta }d\mathrm{ln}V`$. (This yields the effective volume filled by the radiative, “leaky” component of the fluid; the scale-dependent function $`\stackrel{~}{\eta }`$ is nothing but the logarithmic differential ratio between the effective and actual volumes occupied by the gas.)

We consider two cases in turn: 1- an admixture of photon radiation plus an ideal gas (with $`p_{\mathrm{gas}}p_{\mathrm{rad}}`$) and 2- an admixture of neutrinos ($`p_{\mathrm{rad}}=7/8aT^4`$) plus gas/radiation/pairs ($`p_{\mathrm{gas}}=11/12aT^4`$) and where the latter species are collisionally coupled at all scales. Furthermore, for neutrino cooled flows we make the simplification that the relatively cool neucleons ($`T\text{ }<\text{ }20`$ MeV) make a negligible contribution to the total pressure (see, e.g., Di Matteo et al. 2002)

The standard lore (see, e.g., Mihalas & Mihalas 1984, Chandrasekhar 1939) computation of the generalized adiabatic exponent, $`\mathrm{\Gamma }_1`$, for an ideal gas plus a radiative component involves setting $`dQ=dU+dW\mathrm{}`$ while working out expressions for $`dU`$ & $`dW`$ in terms of two logarithmic differentials, say, $`d\mathrm{ln}T`$ and $`d\mathrm{ln}V`$. Using $`dQ\mathrm{}`$ is clearly artificial but this is an adequate artifact when the main concern is with the non-elastic properties of the fluid and not with the amount of heat loss. Thus, we set $`dQ_{\mathrm{tot}}=dQ_{\mathrm{gas}}+d\stackrel{~}{Q}_{\mathrm{rad}}=\mathrm{}`$ with both $`dQ_{\mathrm{gas}}=\mathrm{}`$, and $`d\stackrel{~}{Q}_{\mathrm{rad}}=\mathrm{}`$ (recall that the interactions between the two species are assumed to be entirely elastic). Computation of the specific heat contribution from the gas component is straightforward $`dQ_{\mathrm{gas}}/V=p_{\mathrm{gas}}/(\mathrm{\Gamma }_{\mathrm{gas}}1)d\mathrm{ln}T+p_{\mathrm{gas}}d\mathrm{ln}V`$; while use of the quasi-volume for the radiative component yields

$$d\stackrel{~}{Q}_{\mathrm{rad}}/\stackrel{~}{V}=12p_{\mathrm{rad}}d\mathrm{ln}T+4p_{\mathrm{rad}}d\mathrm{ln}\stackrel{~}{V}=12p_{\mathrm{rad}}d\mathrm{ln}T+4\stackrel{~}{\eta }p_{\mathrm{rad}}d\mathrm{ln}V.$$
(22)
Elimination of $`d\mathrm{ln}T`$ in favor of $`d\mathrm{ln}p`$ is achieved by computing $`dp(d\mathrm{ln}T,d\mathrm{ln}V)`$ for the combined gas. Furthermore, utilizing $`\mathrm{\Gamma }_{\mathrm{gas}}=5/3`$ and the standard definition of $`\mathrm{\Gamma }_1d_{[\mathrm{ln}V]}\mathrm{ln}p(=d_{[\mathrm{ln}\rho ]}\mathrm{ln}p)`$ yields

$$\mathrm{\Gamma }_1=\frac{(4\stackrel{~}{\eta }\beta +1)(4\beta +1)+12\beta +3/2}{(1+\beta )(12\beta +3/2)},$$
(23)
where $`\beta p_{\mathrm{rad}}/p_{\mathrm{gas}}`$. Naturally, the $`\mathrm{}(\mathrm{\Gamma }_1)`$ agrees with the standard result (Chandrasekhar 1939, Mihalas & Mihalas 1984). Moreover, it is relatively straightforward to use the $`\mathrm{}(\mathrm{\Gamma }_1`$) to find the decay rate of acoustic waves partially supported by radiative stress (aside from viscosity).

For neutrino cooled accretion flows, Eq  is still valid under the understanding that $`p_{\mathrm{rad}}=7/8aT^4`$ corresponds to the neutrino pressure, while $`p_{\mathrm{gas}}=11/12aT^4`$ corresponds to photons and pairs and $`dQ_{\mathrm{rad}}/V=12p_{\mathrm{gas}}d\mathrm{ln}T+4p_{\mathrm{gas}}d\mathrm{ln}V`$. With these substitutions, an identical procedure to the one above yields

$$\mathrm{\Gamma }_1=\frac{4}{3}\frac{(1+\stackrel{~}{\eta }\beta )}{(1+\beta )}.$$
(24)
## 5 Impact of Radiative Diffusion at MRI Scales

Thus far, we have aimed to maintain the mathematical simplicity afforded by analytical solutions in order to gain physical insight on the phenomenology of radiative heat conduction and on its impact on the MRI.

We have thus worked out solutions of a simplified version of the dispersion relation for non-axisymmetric MRI modes in the limit of purely horizontal fluid displacements. Since horizontal displacements maximize the efficiency of free-energy tapping from a differential shear flow, this regime generally encompasses the unstable modes of fastest growth. The equation of motion, Eq , is applicable to standard heterogeneous fluids with arbitrary generalized index $`\mathrm{\Gamma }_1`$. The ensuing dispersion relation, Eq , shows that the MRI is sensitive to a combination of two fluid properties–the “softness” of the fluid and the strength of the field–through a single parameter: $`\mathrm{\Lambda }`$. Such parametrization, in turn, bears a direct connection to the compressibility of the modes through Eq . Some physical insight may be gained by re-writing the compressibility parameter, $`1\mathrm{\Lambda }(\mathrm{\Delta }\rho /\rho )`$ as follows (c.f. Eq )

$$1\mathrm{\Lambda }=\left(1+\frac{\mathrm{\Gamma }_1}{2}\frac{p_{\mathrm{r}+\mathrm{g}}}{p__𝐁}\right)^1,$$

while recapitulating that “2” represents the effective index of a magnetic field to perpendicular compression (see, e.g., Shu 1992). MRI modes thus involve a large degree of compression when the effective sound speed of the material medium is much smaller than the Alfvén speed.

On the other hand, in §4 (c.f. Eq ) we have demonstrated that when the disk fluid is very hot such that a radiative, diffusive species (photons or neutrinos) constitutes a significant source of stress, radiative heat conduction precludes a formal treatment of the fluid as a single fluid at the scales of relevance to the MRI; that is, compression at these scales becomes non-adiabatic. A key conceptual point of this paper is that when the diffusive species dominates the stress tensor (e.g., radiation-pressure–dominated disks), the horizontal regime becomes inaccessible because of the brisk loss of pressure support from meridional diffusion<sup>1</sup><sup>1</sup>1 Note that the ordering of wavenumber components in the horizontal regime, $`k_\phi ,k_rk_\theta `$, implies that the smallest length scale of the eddies is meridional; i.e., that the eddies are very nearly flat. . Nevertheless, we are looking to read the gross properties of the fastest growth modes from the 2D dispersion relation Eq  by requiring parity among diffusive and non-diffusive sources of stress while finding an expression for the generalized quasi-adiabatic index of the combined fluid, Eq \[23 & 24\]. It is noteworthy that such parity among sources of stress is to be expected in neutrino cooled accretion flows such as those invoked in gamma-ray burst progenitor models.

To place a zeroth order (recall $`|\stackrel{~}{\omega }|\stackrel{~}{k}^21,\stackrel{~}{\omega }/\stackrel{~}{k}^2𝒪[2i]`$) upper limit on the quasi-adiabatic index, recall that for wave phase $`i(𝐤𝐱\omega t)`$, the MRI growth rate corresponds to a purely imaginary wave frequency (i.e., since the wave phase is stationary at the co-rotation radius): $`\omega _{_{\mathrm{MRI}}}=+i|\omega _{_{\mathrm{MRI}}}|`$. Thus, combining Eqs \[19 & 20\] yields

$$\stackrel{~}{\mathrm{\Gamma }}=\frac{4}{3}\frac{1}{1+\frac{1}{3}\stackrel{~}{k}^2/|\stackrel{~}{\omega }_{_{\mathrm{MRI}}}|}\frac{4}{3}\left(1+\frac{2}{3}\left\{\frac{𝐤^2}{k_{}^2}\right\}\times \frac{\widehat{q}_{_{\widehat{𝐁}}}^2}{|\widehat{\omega }_{_{\mathrm{MRI}}}|}\right)^1.$$
(25)
An accurate estimate of the ratio $`\widehat{q}_{_{\widehat{𝐁}}}^2/|\widehat{\omega }_{_{\mathrm{MRI}}}|`$ for arbitrary values of $`\mathrm{\Lambda }(\mathrm{\Gamma }_1,\mathrm{\Theta })`$ may be found through the iterative use of Eq  (since $`\stackrel{~}{\eta }`$ and thus $`\mathrm{\Gamma }_1`$ itself depends on $`\widehat{q}_{_{\widehat{𝐁}}}^2/|\widehat{\omega }_{_{\mathrm{MRI}}}|`$), but the value of the generalized index is rather insensitive to small changes in this ratio. For Keplerian rotation and in the weak field (i.e., incompressible) limit, we have

$$\widehat{q}_{_{\widehat{𝐁}}}^2/|\widehat{\omega }_{_{\mathrm{MRI}}}|\stackrel{\mathrm{\Theta }0}{}5/4.$$

On the other hand, the ratio $`𝐤^2/k_{}^21`$ for 2D motions so it should be clear that by setting $`𝐤^22k_{}^2`$ we are merely imposing a rough upper limit to the quasi-adiabatic index:

$$\stackrel{~}{\mathrm{\Gamma }}\frac{4}{3}\stackrel{~}{\eta }_{\mathrm{max}}\stackrel{\mathrm{\Theta }0}{}\frac{1}{2}$$
(26)
(again, the result we seek is rather insensitive to the precise value of $`\stackrel{~}{\eta }_{\mathrm{max}}`$).

Using $`\stackrel{~}{\eta }_{\mathrm{max}}=3/8`$, and $`p_{\mathrm{rad}}=p_{\mathrm{gas}}`$, in Eq , the generalized quasi-adiabatic index for an ideal gas plus photon radiation works out to be

$$\mathrm{\Gamma }_11^+.$$
(27)
On the other hand, using $`\stackrel{~}{\eta }_{\mathrm{max}}=3/8`$, $`p_{\mathrm{rad}}=7/8aT^4`$ and $`p_{\mathrm{gas}}=11/12aT^4`$ in Eq , yields

$$\mathrm{\Gamma }_1=\frac{956}{1032}1^{}$$
(28)
for a neutrtino cooled accretion flow. In either case, the fluid behavior is close to isothermal (in a quasi-adiabatic sense) when the quasi-adiabatic index associated with the diffusive species drops below the isothermal value $`\stackrel{~}{\mathrm{\Gamma }}1`$.

## 6 Discussion

Our phenomenological discussion is focused on two broad dynamical properties of the instability on the presumption that these help clarify the physical picture in the non-linear stage. Thus, the growth rate is identified in order of magnitude with the inverse of the correlation time of the turbulence and the geometrical regime of fastest growth, with the shape of the turbulent eddies. We also pay particular attention to the energy deposition process at the outer scale since a major portion of this paper is centered on the issue of entropy generation through non-adiabatic compression when a diffusive species contributes non-trivially to the stress tensor.

On the basis of the discussion in the previous section, we find that the root of a sluggish growth rate for the MRI (discovered by Blaes & Socrates 2001 for axisymmetric modes) when radiative pressure dominates the non-diffusive component is the adiabatic inaccessibility of the horizontal regime associated with fastest growth in standard fluids. For non-axisymmetric modes, this regime becomes intrinsically non-adiabatic because of the brisk rate of heat transport out of anisotropic compressive perturbations. Because we begin with a toroidal field configuration where the non-axisymmetric modes are strongly compressive from the onset, the impact of radiative heat conduction is most dramatic and evidenced first hand from the dispersion relation.

While in an ideal fluid the dominant eddies tend to have $`k_\theta `$ several times larger than $`k_r`$ and $`k_\phi `$, adding a radiative diffusive species will induce conductive losses that move mostly along $`\mathrm{𝟏}_\theta `$. This, in turn, means that the modes with high $`k_\theta `$, i.e., horizontal modes, are preferentially damped. The magnitude of $`k_\theta `$ for the fastest growing modes is then be set by a marginal damping condition. Since the length scale for marginal damping is $`k_{}`$ (c.f. Eq 19), it follows that MRI eddies should be more nearly isotropic in such an environment, which goes along with the reduction in growth rate. For the same reason, energy deposition on to the thermal bath (directly associated with the radiative species, and, secondarily, to the non-radiative species through scattering) occurs on the instability time scale. On the other hand, when the source of stress is equally divided among diffusive and non-diffusive particle species, the net effect is to soften the effective index of the fluid toward isothermality $`\mathrm{\Gamma }_11`$ (c.f. Eq \[27 & 28\]) along with similar but more modest dynamical effects.

These conclusions are broadly supported by the numerical simulations of Turner et al. (2001, 2003). Their report of a standard radiation-pressure–dominated $`\alpha `$-disk, with $`p_{\mathrm{rad}}\text{ }>\text{ }p_{\mathrm{gas}}p__𝐁`$ shows that the non-linear outcome of the MRI is a porous medium with drastic density contrasts. Under nearly constant total pressure and temperature, the non-linear regime shows that density enhancements anti-correlate with azimuthal field domains–just as expected from the linear theory–and that turbulent eddies live for about a dynamical time scale while matter clumps are destroyed through collisions or by running through localized regions of shear on a similar time scale. Notably, the non-linear density contrasts may be quite large $`<\rho _{\mathrm{max}}/\rho _{\mathrm{min}}>\text{ }>\text{ }10`$ when $`p_{\mathrm{rad}}p_{\mathrm{gas}}`$.

At a very fundamental level, clump formation is intimately connected to the effects of radiative heat conduction out of compressive perturbations. This can be understood by computing the polarization properties of the fastest growing modes in the horizontal regime: $`\xi _r=\sqrt{\mathrm{\Lambda }}\xi _\phi \&|\xi _\theta ||\xi _r|`$, and reading from Eq  that at the linear stage of the instability there exists a converging flow toward the Lagrangian displacement node of the modes (see Fig. 2 of Foglizzo & Tagger 1995). In a fluid with entirely elastic (adiabatic) properties, the pressure perturbation associated with such compression will act as a restoring force to de-compress the fluid in the non-linear stage. On the other hand, when the fluid has radiative heat conduction properties on the scale of the density perturbations, no such restoring force persists on time scales longer than the inverse of the growth rate (which, following Eq , is of the same order as the diffusion time). Material clumps thus formed will maintain their integrity and survive for longer times. When the diffusive and non-diffusive pressures are comparable, these arguments still hold but only in part. In this case, compressive perturbations in the non-linear stage will decompress only in proportion to the remaining pressure support from the non-diffusive species.

Blaes & Socrates (2001) predicted qualitatively similar results for axisymmetric modes. Our analytical solutions agree with their result, $`|\widehat{\omega }_{_{\mathrm{MRI}}}|c_g/v_{_{\mathrm{Alf}}}`$, in the limit $`k_r=c_r0`$ and $`c_g0`$. This paper complements these findings and the non-linear numerical explorations of Turner et al. (2001, 2003) by presenting a simplified approximate dispersion relation for non-axisymmetric modes while interpreting the effects of radiative heat conduction in terms of an “ultra-soft” quasi-adiabatic index (i.e., $`\stackrel{~}{\mathrm{\Gamma }}<1`$) for the radiative contribution to the stress tensor. The quasi-adiabatic index and all derived thermodynamic quantities (including “quasi-specific” quantities that ensue from division by the quasi-volume, §4) are mere mathematical artifacts to manipulate thermodynamical relations in the usual manner. Furthermore, we have demonstrated that when there is parity among the sources of stress, the combined gas behaves isothermally. Our results are generally applicable to optically thick, neutrino-cooled disks since the pressure contributions from photons, pairs and neutrinos, all have the same temperature dependence whereas only the neutrino component has radiative heat conduction properties on the time and length scales of the instability.

## Acknowledgments

R.A.G. would like to acknowledge discussions with Omer Blaes and the hospitality of Caltech’s Theoretical Astrophysics and Relativity Group.