ASTR/ATOC-5410: Fluid Instabilities, Waves, and Turbulence August 24, 2016, Axel Brandenburg

Handout 1: Thermal Instability

The Rayleigh-Benard instability is usually characterized as a thermal instability. However, the restoring force in that case is gravity, and so it is in that sense a dynamical instability. True thermal instabilities have a restoring term in the temperature equation. The thermal instability is technically simple and can be discussed even in one dimension. An influential paper on thermal instability was written by Field (1965). Earlier work by Parker (1953) gave a less stringent stability criterion, because he assumed constant density (the isochoric case). This handout gives the basic derivation.

1 Cooling function and entropy equation

The energy equation can be formulated as an evolution equation for the specific entropy,

ρT \frac{Ds}{Dt} = −ρL + ∇ · (K∇T),  (1)

where −L = ρΛ − Γ is the difference between heating and cooling (per unit mass, and −ρL is per unit volume). Note that the cooling per unit volume is proportional to ρ^{2}, reflecting the fact the cooling is accomplished by the collision of two partners. However, in the present context, we will work with temperature and density, so we rewrite this using

DS/c_{p} = \frac{1}{γ} D ln P − D ln ρ (entropy equation)  (2)

and

D ln P = D ln T + D ln ρ (perfect gas equation).  (3)

Thus, we have

DS/c_{p} = \frac{1}{γ} D ln T − (1 − \frac{1}{γ} ) D ln ρ.  (4)

We can then write the energy equation in the form

c_{p}ρT [\frac{1}{γ} \frac{D ln T}{Dt} − (1 − \frac{1}{γ} ) \frac{D ln ρ}{Dt}] = −ρL + ∇ · (K∇T).  (5)

Using γ = c_{p}/c_{v}, we have

c_{v}ρ\frac{DT}{Dt} − c_{p}T (1 − \frac{1}{γ} ) \frac{Dρ}{Dt} = −ρL + ∇ · (K∇T).  (6)

2 Linearized cooling function

Let us linearize L (i.e., Taylor expand to first order) about ρ = ρ_{0} and T = T_{0},

L(ρ, T) ≈ L_{0} + (\frac{∂L}{∂ρ})_{T} (ρ − ρ_{0}) + (\frac{∂L}{∂T} )_{ρ} (T − T_{0}),  (7)

or, in abbreviated form using L_{ρ} = (∂L/∂ρ)_{T}, L_{T} = (∂L/∂T)_{ρ}, we have

L(ρ, T) ≈ L_{0} + L_{ρ}(ρ − ρ_{0}) + L_{T} (T − T_{0}).  (8)

If ρ_{0} and T_{0} refer to an equilibrium state, then L_{0} = 0.

3 Isochoric case (ρ = const)

If ρ = const, i.e., Dρ/Dt = 0 (isochoric case), we can write

c_{v}ρ\frac{DT}{Dt} = −ρL_{0} − ρL_{T} (T − T_{0}) + ∇ · (K∇T).  (9)

There is an equilibrium solution with T = T_{0}, ρ = ρ_{0}, and u = 0, so D/Dt = ∂/∂t provided L_{0} = 0. We can then linearize around this state, so we write T = T_{0} + T_{1}, i.e.,

c_{v}ρ_{0} \frac{∂T_{1}}{∂t} = −ρ_{0}L_{T} T_{1} + K∇^{2}T_{1}.  (10)

where we have assumed K = const. Let us now seek solutions of the form T_{1} = T̂_{1}e^{σt+ik·x}, which leads to the following characteristic equation (or eigenvalue problem),

c_{v}ρ_{0}σ = −ρ_{0}L_{T} − Kk^{2}.  (11)

where σ is the eigenvalue. All eigenvalues are positive and instability occurs when L_{T} is sufficiently negative, i.e., −L_{T} − Kk^{2} > 0, or

ρ_{0}L_{T} < −Kk^{2} (for instability).  (12)

This shows that a necessary condition for instability is

L_{T} ≡ (\frac{∂L}{∂T} )_{ρ} < 0 (for instability).  (13)

A negative L_{T} means that cooling decreases as it gets hotter, so no surprise that this leads to instability. Note, however, that thermal diffusion has a stabilizing effect.

Of course, in a compressible gas there is no reason that ρ should be constant. Another possibility is that the pressure stays constant. In that case, if it does get hotter, the density decreases, so ρL decreases further, making the system even more unstable. To treat this more general case, we have to invoke the hydrodynamic equations, which will be done in the next section.

4 The non-isochoric case

In the non-isochoric case we also consider the consider the linearized continuity equation,

\frac{∂ρ_{1}}{∂t} = −ρ_{0}u_{1},  (14)

where we have restricted ourselves to the 1-D case with u = u(x) = (u(x), 0, 0). We write the pressure gradient in the form

\frac{∂P_{1}}{∂x} = \frac{P_{0}}{T_{0}} \frac{∂T_{1}}{∂x} + \frac{P_{0}}{ρ_{0}} \frac{∂ρ_{1}}{∂x} ,  (15)

so the linearized momentum equation,

ρ_{0} \frac{∂u_{1}}{∂t} = −\frac{∂P_{1}}{∂x} = −\frac{P_{0}}{T_{0}} \frac{∂T_{1}}{∂x} − \frac{P_{0}}{ρ_{0}} \frac{∂ρ_{1}}{∂x} ,  (16)

It is convenient to introduce the (adiabatic) sound speed c_{s} via c_{s}^{2} = γP/ρ. The equilibrium value is referred to as c_{s0}. Thus, we have

\frac{∂u_{1}}{∂t} = −\frac{c_{s0}^{2}}{T_{0}} \frac{∂T_{1}}{∂x} − \frac{c_{s0}^{2}}{ρ_{0}} \frac{∂ρ_{1}}{∂x} ,  (17)

The linearized temperature equation can be written as

c_{v}ρ_{0} \frac{∂T_{1}}{∂t} − c_{p}T_{0} (1 − \frac{1}{γ} ) \frac{∂ρ_{1}}{∂t} = −ρ_{0}L_{ρ}ρ_{1} − ρ_{0}L_{T} T_{1} + K∇^{2}T_{1}.  (18)

Dividing by ρ_{0}c_{v} and T_{0}, we have

\frac{1}{T_{0}} \frac{∂T_{1}}{∂t} − (γ − 1)\frac{1}{ρ_{0}} \frac{∂ρ_{1}}{∂t} = −\frac{L_{ρ}}{c_{v}T_{0}} ρ_{1} − \frac{L_{T}}{c_{v}} T_{1}/T_{0} + γχ∇^{2}T_{1}/T_{0}.  (19)

where χ = K/(ρc_{p}) is the radiative diffusivity (as opposed to the radiative conductivity).^{1}

Collecting the three equations and making the usual ansatz ρ_{1} = ρ̂_{1}e^{σt+ik·x}, u_{1} = û_{1}e^{σt+ik·x}, T_{1} = T̂_{1}e^{σt+ik·x}, leads to the following system of algebraic equations:

σρ_{1}/ρ_{0} + iku_{1} = 0,  (20)
σu_{1} + ikc_{s0}^{2}(ρ_{1}/ρ_{0} + T_{1}/T_{0}) = 0,  (21)
σT_{1}/T_{0} − σ(γ − 1)ρ_{1} + \frac{L_{ρ}}{c_{v}T_{0}} ρ_{1} + \frac{L_{T}}{c_{v}} \frac{T_{1}}{T_{0}} = 0.  (22)

It is useful to define σ_{ρ} = ρ_{0}L_{ρ}/c_{v}T_{0} and σ_{T} = L_{T} /c_{v}, so that we have

(σ + σ_{T}) T_{1}/T_{0} − [σ(γ − 1) − σ_{ρ}] ρ_{1}/ρ_{0} = 0.  (23)

It is convenient to rewrite this in matrix form, Mq = 0, for state vector q = (ρ_{1}/ρ_{0}, u_{1}, s_{1}/T_{0})^{T},

( \begin{matrix}σ & ik & 0 \\ ikc_{s0}^{2} & σ & ikc_{s0}^{2} \\ σ_{ρ} − σ(γ − 1) & 0 & σ + σ_{T} \end{matrix} ) ( \begin{matrix}ρ_{1}/ρ_{0} \\ u_{1} \\ T_{1}/T_{0} \end{matrix} ) = 0.  (24)

which leads to the dispersion relation

σ^{2} (σ + σ_{T}) + k^{2}c_{s0}^{2} [σ(γ − 1) − σ_{ρ}] + k^{2}c_{s0}^{2} (σ + σ_{T}) = 0.  (25)

which simplifies to

σ^{2} (σ + σ_{T}) + k^{2}c_{s0}^{2} (σγ − σ_{ρ} + σ_{T}) = 0.  (26)

5 Power law cooling function

It is convenient to work with piecewise power laws, i.e.,

L = ρΛ(T) − Γ, Λ(T) = Λ_{0}T^{β}.  (27)

Let us also assume that Γ = const, so that

L_{ρ} = Λ_{0} and L_{T} = ρβΛ_{0}/T.  (28)

Since σ_{ρ} = ρ_{0}L_{ρ}/c_{v}T_{0} and σ_{T} = L_{T} /c_{v}, we have

σ_{ρ} = \frac{ρ_{0}Λ_{0}}{c_{v}T_{0}} , σ_{T} = β \frac{ρ_{0}Λ_{0}}{c_{v}T_{0}} = βσ_{ρ}.  (29)

Thus, the dispersion relation becomes

σ^{2} (σ + βσ_{ρ}) + c_{s}^{2}k^{2} [γσ + (β − 1)σ_{ρ}] = 0  (30)

Figure 1: All three roots of Equation (30).

An approximate solution can be obtained by writing

σ^{2} ≈ −c_{s}^{2}k^{2} \frac{γσ + (β − 1)σ_{ρ}}{σ + βσ_{ρ}} → σ^{2} ≈ c_{s}^{2}k^{2} \frac{1 − β}{β} (for σ ≪ σ_{ρ})  (31)

In Figure 1 we compare this approximation with the real and imaginary parts of all three roots.

We see that the instability criterion is β < 1. This covers a broader range of cases than the isochoric stability criterion, which corresponds to β < 0. The following section illustrates the instability for a parameterized (and yet realistic) cooling function. Full details can be found in Brandenburg et al. (2007).

1We’ll return later to the question why it is defined with c_{p} rather than c_{v}.

Figure 2: Net cooling vs. temperature for three values of p, given in units of [p] = 10^{−14} dyn.

Figure 3: The two stable solution branches, ρ_{C} and ρ_{W} (solid lines), and the unstable solution branch, ρ_{U} (dotted line), as a function of p. On the top the pressure is normalized by the Boltzmann constant, p/k_{B}.

Figure 4: Evolution of ln T in a space-time diagram (top) and of the mean pressure (bottom) in a one-dimensional simulation with 1024 mesh points and ν = χ = 5 × 10^{−4} Gyr k_{m}^{2} s^{−2}. During early times, the rms velocity grows exponentially at a rate of about 220 Gyr^{−1}.

Figure 5: Mean pressure vs. mean density in a onedimensional simulation (solid line), compared with the values obtained for a homogeneous system (dotted line). The dotted line (which agrees with the solid line in the stable regime and hence cannot be seen there) was obtained by plotting p(T) vs. ρ(T) using T as a parameter in Equation (33). On the right axis, the pressure is normalized by the Boltzmann constant, p/k_{B}, and at the top the number density is given. The simulation has 128 mesh points and ν = χ = 5 × 10^{−3} Gyr k_{m}^{2} s^{−2}.

6 Piecewise power law cooling function

To illustrate properties of the thermal instability in the interstellar medium, let us adopt a parameterization of the cooling function approximately equal to that given by Sánchez-Salcedo et al. (2002), which has been obtained by fitting a piecewise power law function of the form

Λ(T) = C_{i,i+1}T^{β_{i,i+1}} for T_{i} ≤ T < T_{i+1},  (32)

In Figure 2 we plot L as a function of T for constant p; three values are considered: p = 25, 35, and 50, all in units of [p] ≡ 10^{−14} dyn. This figure shows that there can be two stable states at about 10^{2} and 10^{4} K. We denote these values by T_{C} and T_{W} for the cold and warm phases. At T ≈ 10^{3} K there is an unstable equilibrium, whose temperature is denoted T_{U}. The densities of the three equilibria, obtained by solving L(T; P) = 0 for T numerically for given P and then expressing the result in terms of ρ = ρ(T, P), are plotted in Figure 3.

When the mean density is outside the range between 0.96 and 5.1 (in units of 10^{−24} g cm^{−3}), the gas is thermally stable and remains uniform. The dependence of the pressure on the density can be obtained in parametric form by calculating, using temperature as a parameter, ρ(T) and p(T), that is,

ρ(T) = \frac{Γ}{Λ(T)}, p(T) = \frac{RT}{μ} \frac{Γ}{Λ(T)},  (33)

and plotting the two against each other (see Figure 5, dotted line). The numerically obtained values for the mean pressure ⟨p⟩, for different mean densities ⟨ρ⟩, agree with those obtained under the assumption of homogeneity.

In Figure 4 we plot the evolution of ln T in a space-time diagram (top) and that of the mean pressure in a one-dimensional simulation. Here, ν = χ = 5 × 10^{−4} Gyr k_{m}^{2} s^{−2} which, together with the initial values of c_{s} = 7.5 km/s and n_{ρ} = 980 Gyr^{−1}, yields k_{F} = 720 kpc^{−1} = 23k_{1}, and hence n_{ρ}/(c_{s}k_{F}) ≈ 0.2.

In the unstable regime the pressure is, surprisingly, independent of ⟨ρ⟩, and always around ⟨p⟩ ≈ 24.2 × 10^{−14} dyn. Figure 3 shows that for ⟨p⟩ ≈ 24.2 × 10^{−14} dyn the warm and cool phases have ρ_{W} ≈ 0.19 and ρ_{C} ≈ 14.3, respectively.

References

Brandenburg, A., Korpi, M. J., & Mee, A. J., “Thermal instability in shearing and periodic turbulence,” Astrophys. J. 654, 945-954 (2007).

Field, G. B., “Thermal instability,” Astrophys. J. 142, 531-567 (1965).

Sánchez-Salcedo, F. J., Vázquez-Semadeni, E., & Gazol, A., “The nonlinear development of the thermal instability in the atomic interstellar medium and its interaction with random fluctuations,” Astrophys. J. 577, 768-788 (2002).

E. N. Parker, “Instability of thermal fields,” Astrophys. J. 117, 431-436 (1953).

$Header: /var/cvs/brandenb/tex/teach/ASTR_5410/1_Thermal_Instability/notes.tex,v 1.7 2016/08/25 02:02:35 brandenb Exp $