# Johnson-Nyquist Noise

## Contents

### I. Fluid Basics

#### Fluid Approximation; Fluid Element

Fluids are ‘things that flow,’ and the equations of hydrodynamics treat a fluid as a continuous medium with well-defined macroscopic properties (pressure, density, etc.) at each point. This can be applied to astrophysical systems such as the dynamics of stars in galaxies, the gas inside stars, stellar winds, accretion discs, etc. Some key astrophysical numbers are the following:

\begin{align} \rho _{gas,MW} & \sim 10^{6} \ particle/m^{3} \\ \rho _{gas,Earth} & \sim 2.7 \times 10^{25} \ particles/m^{3} \\ T,\rho _{cold \ molecular \ gas} & \sim 10 \ K, 10^{8} \ particles/m^{3} \\ T,\rho _{warm \ atomic \ gas} & \sim 10^{4} \ K, 10^{5-6} \ particles/m^{3} \\ T,\rho _{hot \ gas \ (ex: \ corona)} & \sim 10^{6} \ K, 10^{3} \ particles/m^{3} \\ M_{MW} & \sim 10^{12} M_{\odot } \\ M_{globular \ cluster} & \sim 10^{5}M_{\odot } \\ M_{galaxy \ cluster} & \sim 10^{14}M_{\odot } \end{align}\,\!

The fluid equations are based on the concept of a ‘fluid element,’ or a region over which we can define local variables. There are some conditions on the size of this region:

1) Small enough so that properties Q are constant: $L \ll \frac{Q}{|\nabla Q|}$
2) Large enough to contain many particles: $nL^{3} \gg 1$
3) Large enough that particles ‘know’ about local conditions through collisions: $L \gg \lambda _{mfp}$

A collisional fluid has a well-defined distribution of particle speeds and therefore a particular equation of state. It is difficult to write down equation of states for collisionless systems and therefore the focus henceforth will be on collisional systems.

#### Derivatives; Eulerian vs. Lagrangian

1) Eulerian: Consider a small volume at a fixed spatial position through which a fluid flows with physical variables specified as functions of time (ex: $\rho = \rho (\vec{r},t)$)
2) Lagrangian: Examining the change in variables of a particular fluid element $\vec{a}$; the reference system is comoving with the fluid (ex: $\rho = \rho (\vec{a},t)$)

The Eulerian approach is usually more useful, though the Lagrangian approach is important if the path of an individual element is important (ex: tracing a parcel of gas after it’s enriched by a supernova). For steady flows (quantities at a given position do not change), the Eulerian approach is particularly useful. A difficulty of Eulerian codes is choosing how to distribute your grid (where there should be fine gridding), and a difficulty of Lagrangian codes is computing properties at each point. Techniques such as Smoothed Particle Hydrodynamics (SPH) and Adaptive Mesh Refinement (AMR) are evolving to address these issues.

The mathematical relationship between the two fluid descriptions is the following, where the left-hand side is the rate of change in the fluid particle’s frame, and the right-hand side accounts for the rate of change at a fixed position and the fluid particle’s motion. The second term is known as the ‘convective derivative’.

\begin{align} \frac{DQ}{Dt} & = \frac{\partial Q}{\partial t} + (\vec{v} \cdot \nabla ) Q \end{align}\,\!

Fluid velocity, or a particle path, is defined as:

\begin{align} \vec{v} & = \frac{d\vec{r}}{dt} \end{align}\,\!

#### Conservation Laws; Continuity Equation

For a fixed volume V with surface element $d\vec{S}$, the rate of change of mass can be set equal to the outward mass flow across the volume:

\begin{align} -\frac{\partial }{\partial t}\int \rho dV & = -\int \rho \vec{v} d\vec{S} \\ & = -\int \nabla (\rho \vec{v})dV \\ \frac{\partial \rho }{\partial t} + \nabla (\rho \vec{v}) & = 0 \\ \frac{D\rho }{Dt} +\rho (\nabla \cdot \vec{v}) & = 0 \end{align}\,\!

For an incompressible flow, Dρ / Dt = 0, and fluid elements preserve their density along their paths. Therefore, $\nabla \cdot \vec{v} = 0$ and incompressible flows are divergence free.

#### Pressure; Stress Tensor; Momentum/Force/Euler Equation

In order to conserve momentum, the forces between fluid elements need to be considered. The force across a surface is related to the stress tensor (force per area in i acting on a surface whose normal is j) by Fi = TijSj. The net force acting over the surface due to pressure is:

\begin{align} \vec{F} & = -\int P d\vec{S} \\ & = -\int \nabla P dV \end{align}\,\!

Therefore, the acceleration (force per mass) contribution from pressure is $-\frac{\nabla P}{\rho }$. Note that the force and pressure gradient act in opposite directions.

There is an additional acceleration term due to gravity as well. Putting things together, the Euler/momentum equation is:

\begin{align} \frac{\vec{F}}{m} & = \frac{D\vec{v}}{Dt} \\ & = -\frac{\nabla P}{\rho } + \vec{g} \\ \frac{D\vec{v}}{Dt} & = -\frac{\nabla P}{\rho } - \nabla \Phi \\ \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla )\vec{v} & = -\frac{\nabla P}{\rho } - \nabla \Phi \end{align}\,\!

The above equations state that the momentum of a fluid element changes as a result of pressure gradients and gravitational forces. Remember that force is a change in momentum.

Switching to tensor notation can illuminate the concept of ram pressure. The rate of change of momentum density can be expanded and its terms can be replaced using the Continuity and Euler equations:

\begin{align} \frac{\partial }{\partial t}(\rho v_{i}) & = \rho \frac{\partial v_{i}}{\partial t} + v_{i}\frac{\partial \rho }{\partial t} \\ & = -\rho v_{j}\partial _{j}v_{i} - \partial _{i}P - \rho \partial _{i}\Phi - v_{i}\partial _{j}(\rho v_{j}) \\ & = -\partial _{j}(\rho v_{j}v_{i}+P\delta _{ij})-\rho \partial _{i}\Phi \\ \frac{\partial }{\partial t}(\rho v_{i}) & = -\partial _{j}T_{ij}-\rho \partial _{i}\Phi \end{align}\,\!

This is the conservative form of the momentum equation, where Tij = ρvivj + Pδij (units of momentum flux, or pressure). The first term in this stress tensor represents ram pressure due to the bulk motion of a fluid (pressure on a body moving through a fluid) and the second term represents thermal pressure associated with random motions that are isotropic. Thermal pressure is independent of i and j and acts perpendicular to any surface in the fluid.

The gravitational force on a fluid in astronomical situations is very important, and gravity can be defined by $\vec{g} = -\nabla \Phi = -\frac{GM}{r^{2}}\hat{r}$.

\begin{align} \int \nabla \Phi d\vec{S} & = \int \frac{GM}{r^{2}}\hat{r}d\vec{S} \\ \int \nabla ^{2} \Phi dV & = \int \frac{GM}{r^{2}}\hat{r}\frac{r^{2}d\Omega }{\hat{r}} \\ & = \int GM4\pi \\ & = 4\pi G \int \rho dV \\ \nabla ^{2}\Phi & = 4\pi G\rho \end{align}\,\!

The equation above is Poisson’s equation, which relates any mass density distribution to a potential.

#### Equation of State; Energy Equation; Bernoulli’s Principle

So far, we have four unknowns, $\Phi , \rho , \vec{v},$ and P, and three equations (Continuity, Euler, Poisson). We need an equation of state relating P with other variables in order to close the system.

Most of the fluids in the Universe can be approximated as ideal gases (warm and dilute). The first law of thermodynamics states that $dQ = TdS = d\varepsilon + P dV$, where dQ is the heat absorbed by the fluid, PdV is the work done by the fluid, and $d\varepsilon$ is the internal energy content of the fluid. The specific heat equations are:

\begin{align} C_{V} & = \frac{\partial Q}{\partial T}\Big|_{V} = \frac{\partial \varepsilon }{\partial T}\Big|_{V} \\ C_{P} & = \frac{\partial Q}{\partial T}\Big|_{P} \end{align}\,\!

Therefore, dQ = CVdT + PdV. To find how the specific heats are related to each other, the equation of state for an ideal gas is substituted in.

\begin{align} PV & = Nk_{B}T \\ P dV + VdP & =& Nk_{B}T \\ dQ & = C_{V}dT + Nk_{B}dT-VdP \\ C_{P} & = C_{V} + Nk_{B} \end{align}\,\!

For an adiabatic ideal gas, there is no heating or cooling due to dissipative processes, so dQ = 0.

\begin{align} C_{V}dT + PdV & = 0 \\ C_{V}dT + \frac{Nk_{B}T}{V}dV & = 0 \\ T & \propto V^{1-\gamma } \\ P & \propto V^{-\gamma } \\ P & \propto T^{\frac{\gamma }{\gamma -1}} \end{align}\,\!

These are adiabatic relationships between V, T, and P, for $\gamma = \frac{C_{P}}{C_{V}}$. Also, since $\rho \propto V^{-1}$, P = kργ. This is a barotropic (function of density) equation of state.

For a monatomic gas, $\varepsilon = \frac{3}{2}Nk_{B}T$, $C_{V} = \frac{3}{2}Nk_{B}$, $C_{P} = \frac{5}{2}Nk_{B}$, and γ = 5 / 3. For a diatomic gas, $\varepsilon = \frac{5}{2}Nk_{B}T$, $C_{V} = \frac{5}{2}Nk_{B}$, $C_{P} = \frac{7}{2}Nk_{B}$, and γ = 7 / 5. For an isothermal gas, temperature remains constant and γ = 1.

The total energy density (energy per volume) of a fluid is $E = \rho (\frac{1}{2}v^{2}+\Phi + \epsilon )$, where ε is the internal energy per mass. The first law of thermodynamics can be tweaked (per mass) and the Continuity and Euler equations can be substituted in to obtain the energy conservation equation.

\begin{align} d\epsilon & = dQ - PdV \\ \frac{D\epsilon }{Dt} & = -P\frac{DV}{Dt} - \dot{Q}_{cool} \\ & = \frac{P}{\rho ^{2}}\frac{D\rho }{Dt} - \dot{Q}_{cool} \\ \frac{DE}{Dt} & = \frac{E}{\rho }\frac{D\rho }{Dt} + \rho \vec{v}\frac{D\vec{v}}{Dt} + \rho \frac{D\Phi }{Dt} + \rho \frac{D\epsilon }{Dt} \\ & = \frac{E}{\rho }\frac{D\rho }{Dt} + \rho \vec{v}\frac{D\vec{v}}{Dt} + \rho \frac{D\Phi }{Dt} + \frac{P}{\rho }\frac{D\rho }{Dt} - \rho \dot{Q}_{cool} \\ \frac{\partial E}{\partial t} + \nabla [(E+P)\vec{v}] & = -\rho \dot{Q}_{cool} + \rho \frac{\partial \Phi }{\partial t} \end{align}\,\!

The cooling function is positive if the medium is cooled and negative for external heating. Possible energy transport mechanisms include conduction (transfer of heat from warm areas to cooler ones as a consequence of random particle motions), convection (transfer of energy due to temperature gradients), and radiation (energy carried by photons).

Recalling that enthalpy is $H = \varepsilon + PV$, h can be defined as enthalpy per unit mass: $h = \epsilon + \frac{P}{\rho }$. Therefore, the energy flux is $(E+P)\vec{v} = \rho \vec{v} (\frac{1}{2}v^{2}+\Phi +h)$.

To recap, the three main quantities that are conserved are mass m, momentum $m\vec{v}$, and energy $\varepsilon$. Each of their conservation laws involve the rate of change of their densities (ρ, $\rho \vec{v}$, and E, respectively). Their fluxes are given by their densities times a velocity ($\rho \vec{v}$, Tij, and $(E+P)\vec{v}$, respectively).

#### Simple Solutions; Hydrostatic Equilibrium; Polytropic Star; Lane-Emden Equation

For the case of hydrostatic equilibrium, $\vec{v} = 0$ and $\frac{\partial }{\partial t} = 0$. Therefore, the simplified momentum and Poisson equations are:

\begin{align} \frac{\nabla P}{\rho } & = -\nabla \Phi \\ \nabla ^{2}\Phi & = 4\pi G\rho \end{align}\,\!

An equation of state is needed to close the system. For an isothermal atmosphere affected by only gravity:

\begin{align} P & = \frac{k_{B}T}{\mu m_{p}}\rho = A\rho \\ \frac{\nabla P}{\rho } & = \vec{g} = -g \hat{z} \\ \frac{A}{\rho }\nabla \rho & = -g \\ \rho (z) & = \rho _{0}e^{-z/z_{s}} \\ z_{s} & = \frac{k_{B}T}{\mu m_{p}g} \end{align}\,\!

Another case of hydrostatic equilibrium is that of a star. Stars can be modeled as self-gravitating polytropes where $\frac{1}{\rho }\frac{dP}{dr} = -\frac{d\Phi }{dr}$. Since ρ = ρ(Φ) and P = P(Φ), it follows that P = P(ρ), and therefore stars are barotropes. The barotrope equation of state can be parametrized as a polytrope $P = k\rho ^{1+\frac{1}{n}}$, where n is the polytropic index. To solve for this case, the polytropic equation of state should be substituted into the equation of hydrostatic equilibrium.

\begin{align} \frac{\nabla P}{\rho } & = \frac{1}{\rho }\nabla (k\rho ^{1+\frac{1}{n}}) \\ (n+1)\nabla (k \rho ^{1/n}) & = -\nabla \Phi \\ \rho & = \Big[\frac{\Phi _{0}-\Phi }{(n+1)k}\Big]^{n} \\ \rho _ c & = \Big[\frac{\Phi _{0}-\Phi _{c}}{(n+1)k}\Big]^{n} \\ \rho & = \rho _{c}\Big[\frac{\Phi _{0}-\Phi }{\Phi _{0}-\Phi _{c}}\Big]^{n} \end{align}\,\!

In the formulas above, the subscript 0’s refer to the surface of the star, the subscripts c’s refer to inside the star, and no subscript refers to outside the star. Defining $\theta = \frac{\Phi _{0}-\Phi }{\Phi _{0}-\Phi _{c}} = \Big(\frac{\rho }{\rho _{c}}\Big)^{1/n}$ and using Poisson’s equation:

\begin{align} \nabla ^{2} \Phi & = 4\pi G \rho \\ \nabla ^{2} \theta & = -\frac{4\pi G\rho _{c}}{\Phi _{0}-\Phi _{c}}\theta ^{n} \end{align}\,\!

Because of spherical symmetry, a new variable can be defined: $\xi = \Big(\frac{4\pi G\rho _{c}}{\Phi _{0}-\Phi _{c}}\Big)^{1/2}r = \Big(\frac{4\pi G}{(n+1)k}\rho _{c}^{1-\frac{1}{n}}\Big)^{1/2}r$ and the equation becomes the Lane-Emden equation.

\begin{align} \frac{1}{\xi ^{2}} \frac{d}{d\xi }\Big(\xi ^{2}\frac{d\theta }{d\xi }\Big) & = -\theta ^{n} \end{align}\,\!

The boundary conditions for the Lane-Emden equation are that $\theta =1, \frac{d\theta }{d\xi }=0, \rho =\rho _{c}$ at ξ = 0 (center) and ρ = 0,Φ = Φ0 at θ = 0 (surface).

There are analytic solutions to the Lane-Emden equation for n = 0,1,5 that involve some mathematical simplifications and use of the boundary conditions. The key result is that the radius of the star (where θ = 0) increases with increasing n.

A useful mass-radius scaling relation comes out of the Lane-Emden equation, since all stars with the same n have the same θ(ξ).The central density ρc sets the physical scale of ρ and r.

\begin{align} M & = \int 4\pi \rho r^{2} dr \\ & = 4\pi \Big(\frac{4\pi G\rho _{c}^{1-\frac{1}{n}}}{k(n+1)}\Big)^{-3/2} \int \theta ^{n} \xi ^{2} d\xi \\ M & \propto \rho _{c}^{\frac{1}{2}(\frac{3}{n}-1)} \end{align}\,\!

From the definition of ξ, $r \propto \rho _{c}^{\frac{1}{2}(\frac{1}{n}-1)}$, so the mass-radius relation for polytropic stars is found to be $M \propto R^{\frac{3-n}{1-n}}$. For an adiabatic star, $1+\frac{1}{n}=\gamma$ and γ = 5 / 3 for a monatomic gas, so $M \propto R^{-3}$. This is characteristic of non-relativistic degenerate stars, like white dwarfs.

Going back to the Euler equation, it can be re-written using the definition of vorticity $\vec{w} = \nabla \times \vec{v}$.

\begin{align} \frac{\partial \vec{v}}{\partial t} + \nabla (\frac{1}{2}\vec{v}^{2}+\Phi ) + \frac{\nabla P}{\rho }-\vec{v} \times \vec{w} & = 0 \end{align}\,\!

There are two special cases that can be applied here. The first is the case of steady flows of ideal fluids. Recalling enthalpy, $h = \epsilon + \frac{P}{\rho }$, so $dh = d\epsilon + Pd(\frac{1}{\rho })+\frac{1}{\rho }dP = Tds + \frac{1}{\rho }dP$, where everything is per mass. Therefore, for constant S (ideal fluids), $\nabla h = \frac{1}{\rho } \nabla P$. Re-writting the Euler equation, dotting it with $\vec{v}$, and using this substitution yields Bernoulli’s function, which is in parentheses below. It is comprised of three terms: kinetic energy per mass, gravitational energy per mass, and enthalpy per mass.

\begin{align} \vec{v} \cdot \nabla (\frac{1}{2}\vec{v}^{2} + \Phi + h) & = 0 \end{align}\,\!

For steady flows, $\frac{\partial \mathcal{B}}{\partial t} = 0$, and since $\vec{v} \cdot \nabla \mathcal{B} = 0$, then $\frac{D\mathcal{B}}{Dt} = 0$ and Bernoulli’s function is constant along fluid elements.

The second case to consider is irrotational flows, for which $\vec{w} = 0$. In this case $\nabla \mathcal{B} = 0$ and Bernoulli’s function is constant everywhere. This means that the flow is curl-free. In general, the Euler equation can be written as the following:

\begin{align} \frac{\partial \vec{v}}{\partial t} + \nabla (\frac{1}{2}\vec{v}^{2}+\Phi + h) - \vec{v} \times \vec{w} & = 0 \\ \frac{\partial \vec{v}}{\partial t} + \nabla \mathcal{B} - \vec{v} \times \vec{w} & = 0 \end{align}\,\!

Taking the cross product of $\nabla$ with the above yields $\frac{\partial \vec{w}}{\partial t} = \nabla \times (\vec{v} \times \vec{w})$, which is Helmholtz’s equation. It implies that vorticity flux is conserved.

### II. Waves

#### Sound Waves

Sound waves are important in astrophysics because they provide the principal mechanism by which disturbances propagate in fluids. Sound-crossing timescales are often used to determine whether given regions have had time to respond dynamically to disturbances. In order to derive the familiar wave dispersion relation, the fluid equations are employed with perturbations. For steady-state, no gravity, and $\vec{v}_{0} = 0$, the perturbations are:

\begin{align} \vec{v} & = \vec{v}_{1} \\ P & = P_{0} + P_{1} \\ \rho & = \rho _{0} + \rho _{1} \end{align}\,\!

It’s noted that these perturbations are Lagrangian, but they will be substituted into the Eulerian fluid equations. In this case, the unperturbed quantities are uniform so $\frac{DQ}{Dt} = \frac{\partial Q}{\partial t}$. The first order fluid equations are:

\begin{align} \frac{\partial \rho _{1}}{\partial t} + \rho _{0}\nabla \vec{v}_{1} & = 0 \\ \frac{\partial \vec{v}_{1}}{\partial t} & = -\frac{1}{\rho _{0}}\nabla P_{1} \end{align}\,\!

Taking $\frac{\partial }{\partial t}$ of the continuity equation and substituting the momentum equation yields $\frac{\partial ^{2}\rho _{1}}{\partial t^{2}} - \nabla ^{2}P_{1} = 0$. Assuming the fluid is barotropic and a change in density maps to a change in pressure, the sound speed can be defined as $v_{s}^{2} = \frac{\partial P}{\partial \rho } = \frac{P_{1}}{\rho _{1}}$. Therefore, the equation simplifies to the free wave equation:

\begin{align} \frac{\partial ^{2}\rho _{1}}{\partial t^{2}} - v_{s}^{2}\nabla ^{2} \rho _{1} & = 0 \end{align}\,\!

The solution to the free wave equation is $\rho _{1}(\vec{r},t) = \int d^{3}k e^{i(\vec{k}\cdot \vec{r}-\omega t)} \cdot c(\vec{k})$. Plugging this into the equation yields the dispersion relation $v_{s}^{2} = \frac{\omega ^{2}}{k^{2}}$. The group velocity is defined to be $v_{g} = \frac{\partial \omega }{\partial k} = v_{s}$ and the phase velocity is defined to be $v_{\phi } = \frac{\omega }{k} = v_{s}$, so in this case vg = vφ, and these velocities are independent of k (non-dispersive).

A similar process can be followed to obtain the wave equation for $\vec{v}_{1}$, which has a similar solution. A relation can be found between the two perturbations and the solutions for $\vec{v}_{1}$ and ρ1 can be plugged in.

\begin{align} \frac{\partial \vec{v}_{1}}{\partial t} & = -\frac{v_{s}^{2}}{\rho _{0}}\nabla \rho _{1} \\ -i\omega \vec{v}_{1} & = -\frac{v_{s}^{2}}{\rho _{0}} (ik)\rho _{1} \\ \vec{v}_{1} & = \frac{\omega }{k}\frac{\rho _{1}}{\rho _{0}} \hat{k} \end{align}\,\!

From this equation, the fluid velocity and density perturbations are in phase. Conceptually, sound waves are density perturbations which create pressure gradients which create a force and therefore acceleration, or propagation of the disturbance. For a stiff equation of state (high $\frac{dP}{d\rho }$), there is a large restoring force for small density perturbations and implies rapid propagation.

Some practical sound speed examples include the isothermal sound speed, $v_{s}^{2} = \frac{\partial P}{\partial \rho } |_{T} = \frac{k_{B}T}{\mu m_{p}}$ and adiabatic sound speed, $v_{s}^{2} = \frac{\partial P}{\partial \rho } |_{S} = \gamma \frac{k_{B}T}{\mu m_{p}}$. Therefore, $v_{s,adiabatic} = \sqrt {\gamma }v_{s,isothermal}$. For air (diatomic), γ is greater than 1 and therefore the adiabatic sound speed is greater than the isothermal sound speed.

#### Gravity Waves; Surface Water Waves; Capillary Waves

If there are external forces present, such as sound waves propagating in an isothermal atmosphere with constant gravity, the z-direction equations of motions will be affected. The perturbations are the following:

\begin{align} \vec{v} & = \vec{v}_{1} \\ P & = P_{0} + P_{1} \\ \rho & = \rho _{0} \\ \nabla \Phi _{0} & = -\vec{g} \\ \nabla \Phi _{1} & = 0 \end{align}\,\!

The 0th order momentum equation is:

\begin{align} \frac{1}{\rho _{0}}\nabla P_{0} & = -\nabla \Phi _{0} = \vec{g} = -g \hat{z} \\ P_{0} & = -\rho _{0}gz + const \end{align}\,\!

The 1st order fluid equations are::

\begin{align} \rho _{0} \nabla \vec{v}_{1} & = 0 \\ \frac{\partial \vec{v}_{1}}{\partial t} & = -\frac{1}{\rho _{0}}\nabla P_{1} \end{align}\,\!

For waves traveling along $\hat{x}$, the components of the equations are:

\begin{align} \frac{\partial }{\partial x}v_{1x} + \frac{\partial }{\partial z}v_{1z} & = 0 \\ \frac{\partial }{\partial t}v_{1x} & = -\frac{1}{\rho _{0}} \frac{\partial }{\partial x}P_{1} \\ \frac{\partial }{\partial t}v_{1z} & = -\frac{1}{\rho _{0}} \frac{\partial }{\partial z}P_{1} \end{align}\,\!

Putting the equations together yields Laplace’s equation:

\begin{align} \Big(\frac{\partial ^{2}}{\partial x^{2}}+\frac{\partial ^{2}}{\partial z^{2}}\Big) P_{1} & = 0 \end{align}\,\!

The solution to this is $P_{1}(x,z,t) \propto e^{i(kx-\omega t)}f(z)$. Plugging this into Laplace’s equation yields solutions for f(z) and v1z and v1x.

\begin{align} -k^{2}f+\frac{\partial ^{2}f}{\partial z^{2}} & = 0 \\ f(z) & \propto e^{\pm kz} \\ P_{1}(x,z,t) & = Ae^{i(kx-\omega t)}(e^{kz} + Be^{-kz}) \\ v_{1z}(x,z,t) & = \frac{Ak}{\rho _{0}i\omega }e^{i(kx-\omega t)}(e^{kz} - Be^{-kz}) \\ v_{1x}(x,z,t) & = \frac{Ak}{\rho _{0}\omega }e^{i(kx-\omega t)}(e^{kz} + Be^{-kz}) \end{align}\,\!

There’s a boundary condition that v1z = 0 at the bottom of the sea (z = − H). This implies that B = e − 2kH. Recalling 2cosh(x) = ex + ex and 2sinh(x) = exex, plugging this into the previous solutions gives:

\begin{align} P_{1}(x,z,t) & = Ae^{i(kx-\omega t)}e^{-kH}2cosh[k(z+H)] \\ v_{1z}(x,z,t) & = \frac{Ak}{\rho _{0}i\omega }e^{i(kx-\omega t)}e^{-kH}2sinh[k(z+H)] \\ v_{1x}(x,z,t) & = \frac{Ak}{\rho _{0}\omega }e^{i(kx-\omega t)}e^{-kH}2cosh[k(z+H)] \end{align}\,\!

Before deriving the dispersion relation, the concept of a fluid displacement $\vec{\xi }(\vec{r},t)$ must be introduced. The perturbed fluid velocity is related to this displacement by $\vec{v}_{1} = \frac{d\vec{\xi }}{dt}$, so $\xi _{z} = \frac{v_{1z}}{-i\omega }$. Also, at the air-water interface (z = 0), a boundary condition is that there is constant pressure along a fluid element. Mathematically, $P_{1} + \vec{\xi } \cdot \nabla P_{0} = 0$ at z = 0. The proof of this comes from the relationship between Eulerian and Lagrangian perturbations. A Eulerian perturbation follows a fixed point at different times: $P(\vec{r},t) = P_{0}(\vec{r})+P_{1}(\vec{r},t)$, where the second term on the right is the Eulerian perturbation. A Lagrangian perturbation follows a fluid element undergoing a displacement: $P(\vec{r}+\vec{\xi },t) = P_{0}(\vec{r})+P_{1L}$. Solving for P1L gives $P_{1L} = P(\vec{r},t)+\vec{\xi }\nabla P - P_{0}(\vec{r}) = P_{1}(\vec{r},t) + \vec{\xi } \nabla P$. Since there is constant pressure along a fluid element at the interface, P1L = 0, and the surface boundary condition follows:

\begin{align} P_{1} + \vec{\xi } \nabla P_{0} & = 0 \end{align}\,\!

The dispersion relation can be found for these gravity waves by applying this condition.

\begin{align} P_{1} + \vec{\xi }_{z}\frac{dP_{0}}{dz} & = 0 \\ cosh[k(z+H)] + \frac{k}{\rho _{0}\omega ^{2}}sinh[k(z+H)](-\rho _{0}g) |_{z=0} & = 0 \\ cosh[kH] = \frac{gk}{\omega ^{2}}sinh[kH] \\ \omega ^{2} = gktanh[kH] \end{align}\,\!

There are two limits worth discussing. The first is the deep water/short wavelength limit: $kH \gg 1$ or $H \gg 1/k \sim \lambda$. In this case, tanh[kH]˜1 and $\omega = \sqrt {gk}$. The phase velocity in this regime is $v_{\phi } = \sqrt {\frac{g}{k}}$, which implies that the wave is dispersive and shorter wavelengths travel more slowly. There is also no dependence on H, which means that small disturbances don’t know about the sea floor. The second regime is the shallow water/long wavelength limit: $kH \ll 1$ or $H \ll \lambda$. In this case, tanh[kHkH and $\omega = \sqrt {gH}k$. The phase and group velocities are equal to each other: $v_{\phi } = v_{g} = \sqrt {gH}$. This wave is non-dispersive. For example, tsunamis behave as shallow water gravity waves and all tsunami waves travel at the same speed.

#### Shock Waves; Jump Conditions

If a fluid is subject to a non-linear disturbance (ex: compression by a large factor or a large acceleration), the result is a shock. Astrophysics shock phenomena include gas free-falling onto the surface of stars (gravity produces high speeds and then this fluid has a sudden deceleration when it meets other fluid, which produces a shock) or colliding interstellar clouds.

For a fluid that flows faster than the disturbance speed vs, the Mach number is $M = \frac{v}{v_{s}} = \frac{1}{sin \alpha }$, where α is the angle of the Mach cone. The flow is undisturbed until it reaches the disturbance and then its properties change discontinuously in a shock. Typically, a cold fast flow becomes a hot, high pressure, high density, slow flow. The conservation of mass, momentum, and energy still apply for shocks and can be derived from the fluid equations. First is mass conversation (take the integral over dx of the continuity equation), noting that the mass flux in equals the mass flux out:

\begin{align} \frac{\partial \rho }{\partial t} + \frac{\partial }{\partial x} (\rho v_{x}) & = 0 \\ \frac{\partial }{\partial t} \int \rho dx + \rho v_{x} |_{\Delta {x}/2} - \rho v_{x} |_{-\Delta {x}/2} & = 0 \\ \rho _{1}v_{1} = \rho _{2}v_{2} \end{align}\,\!

Applying the same method to the conservative form of the momentum equation (and noting that Φ is continuous across shocks) yields:

\begin{align} \rho _{1}v_{1}^{2}+P_{1} & = \rho _{2}v_{2}^{2}+P_{2} \end{align}\,\!

The momentum conservation equation represents how a shock converts ram pressure to thermal pressure. And finally, using the energy conservation equation and noting that $\dot{Q}_{cool} = 0$:

\begin{align} \frac{1}{2}v_{1}^{2}+\epsilon _{1}+\frac{P_{1}}{\rho _{1}} & = \frac{1}{2}v_{2}^{2}+\epsilon _{2}+\frac{P_{2}}{\rho _{2}} \end{align}\,\!

The energy conservation equation represents the conversion of kinetic energy into enthalpy. For an adiabatic ideal gas, dQ = 0 so $d\epsilon = -Pd(\frac{1}{\rho })$. Since P = kργ, the differential equation can be solved to get $\epsilon = {\frac{1}{\gamma -1}}\frac{P}{\rho }$. Recall that ε is energy per mass, while $\varepsilon = \frac{1}{\gamma -1}P$ is energy per volume.

The expression for ε can be substituted into the energy conservation equation, and $v_{s}^{2} = \gamma \frac{P_{1}}{\rho _{1}}$ for an adiabatic gas. Re-arranging and substituting amongst the three conservation equations (the Rankine-Hugoniot relations) produces the following relationships:

\begin{align} M_{1}^{2} & = \frac{v_{1}^{2}\rho _{1}}{\gamma P_{1}} \\ \frac{v_{1}}{v_{2}} & = \frac{\rho _{2}}{\rho _{1}} = \frac{\gamma +1}{\gamma -1+2/M_{1}^{2}} \\ \frac{P_{2}}{P_{1}} & = \frac{2\gamma M_{1}^{2}-(\gamma -1)}{\gamma +1} \\ M_{2}^{2} & = \frac{(\gamma -1)M_{1}^{2}+2}{2\gamma M_{1}^{2}-(\gamma -1)} \end{align}\,\!

In the strong shock limit, $M_{1} \gg 1$, $\frac{v_{1}}{v_{2}} = \frac{\gamma +1}{\gamma -1}$, $\frac{P_{2}}{P_{1}} \gg 1$, $\frac{T_{2}}{T_{1}} \gg 1$, and $T_{2} \sim \frac{2(\gamma -1)\mu m_{p}}{(\gamma +1)^{2}k_{B}}v_{1}^{2}$ (really high temperature).

For isothermal shocks, the energy conservation equation doesn’t apply since heat is lost. Instead, T1 = T2 and the sound speed is the same before and after the shock. Re-arranging and substituting the first two conservation equations shows that the density ratio can have very high values for high Mach numbers, unlike the adiabatic case in which the ratio was capped.

#### Blast Waves; Sedov-Taylor Similarity Solutions

Supernovae release around $10^{44} \ J$ of thermal and kinetic energy on minute timescales. Their shocks expand and sweep up gas, creating large bubbles in the ISM. These strong explosions can be characterized as blast waves of energy E (it goes from adiabatic expansion to a blast wave after the shock has swept up a mass similar to the initial ejecta mass). There is a pre-shock density ρ1 before the explosion and a post-shock density ρ2 inside the blast wave shell. In the pre-shock’s rest frame, the post-shock velocity is $v_{2}’ = v_{1}-v_{2} = v_{1}(1-\frac{\gamma -1}{\gamma +1}) = \frac{2}{\gamma +1}v_{1}$. As the shock moves outward, most mass within a radius R gets swept up into a shell of thickness D. For a monatomic gas with γ = 5 / 3:

\begin{align} {\frac{4}{3}}\pi R^{3} \rho _{1} & = 4\pi R^{2}D\rho _{2} \\ D & = \frac{R}{3}\frac{\rho _{1}}{\rho _{2}} \\ & = \frac{R}{3}\frac{v_{2}}{v_{1}} \\ & = \frac{R}{3}\frac{\gamma -1}{\gamma +1} \\ & \sim \frac{R}{12} \\ D & \ll & R \end{align}\,\!

The shell is very thin compared to the size of the explosion. The explosion energy E is $\frac{1}{2}mv_{1}^{2} = \frac{1}{2}\frac{4}{3}\pi R^{3}\rho _{1}\dot{R}^{2}$. Therefore, $E \propto \frac{\rho _{1}R^{5}}{t^{2}}$. Solving for R shows how the shock evolves with time:

\begin{align} R & = A\Big(\frac{E}{\rho _{1}}\Big)^{1/5}t^{2/5} \end{align}\,\!

This scaling holds if $M_{1} \gg 1$ and if the swept up mass is greater than the mass of ejecta. In other words, it applies to a shock that has been around long enough to sweep up a lot of mass, but not long enough so that energy loss becomes important. Radiative losses become important when the age of the system is higher than the cooling timescale of the post-shock gas.

The blast wave is expected to undergo a spherically-symmetric evolution in space and time, and the only parameters in the problem are the explosion energy and external density. Therefore, the fluid equations can be written in non-dimensional variables using scale coordinates and the fluid evolution is understood as a similarity solution (the flow at any location and any time should look the same as it did at some other location at an earlier time). A dimensionless distance parameter can be defined: $\xi = \frac{r}{R(t)}$. Therefore, any parameter X can be plotted as a function of ξ instead of r and it will have the same shape but be scaled up or down by a time-dependent factor: $X = X_{2}(t) \tilde{X}(\xi )$. In order to find these self-similar solutions, variables are changed from (r,t) to (ξ,R).

\begin{align} \frac{\partial X}{\partial r} & = X_{2}\frac{d\tilde{X}}{d\xi }\frac{\partial \xi }{\partial r} |_{t} \\ \frac{\partial X}{\partial t} & = \tilde{X}(\xi )\frac{dX_{2}}{dt} + X_{2} \frac{d\tilde{X}}{d\xi }\frac{\partial \xi }{\partial t} |_{r} \end{align}\,\!

For example, $r = \xi \Big(\frac{Et^{2}}{\rho _{1}}\Big)^{1/5}$, so the velocity of the shock is $v = \frac{dr}{dt} = {\frac{2}{5}}\xi \Big(\frac{E}{\rho _{1}t^{3}}\Big)^{1/5}$. If the velocity is known/calculated, then t can be solved for as a function of ξ, and therefore r is known as a function of ξ, as well as M and other parameters depending on r.

### III. Fluid Instabilities

#### Gravitational Instability in 3D; Jean’s Length; Static vs. Expanding Medium

For a fluid in a steady state, small perturbations can grow in time (unstable system), diminish (stable system), or oscillate. Fluid instabilities are responsible for convection in stars, the formation of galaxies and stars in the Universe, etc.

In a static, 3D medium, suppose the perturbations are:

\begin{align} \vec{v} & = \vec{v}_{1} \\ P & = P_{0} + P_{1} \\ \rho & = \rho _{0} + \rho _{1} \\ \Phi & = \Phi _{0} + \Phi _{1} \end{align}\,\!

The 0th order fluid equations are:

\begin{align} \frac{\partial \rho _{0}}{\partial t} + \nabla (\rho _{0}\vec{v}_{0}) & = 0 \\ \frac{\partial \vec{v}_{0}}{\partial t} + (\vec{v}_{0} \cdot \nabla )\vec{v}_{0} & = -\frac{1}{\rho _{0}}\nabla P_{0} - \nabla \Phi _{0} \\ \nabla ^{2} \Phi _{0} & = 4\pi G \rho _{0} \end{align}\,\!

The first equation is trivial, since $\vec{v}_{0}=0$ and the system is in steady-state. The left-hand side of the momentum equation is also zero, as well as the first term on the right-hand side since P0 is constant in time and space. Therefore, that leaves $\nabla \Phi _{0} = 0$ which contradicts Poisson’s equation that states $\nabla \Phi _{0} = \frac{4}{3}\pi G\rho _{0}$. This discrepancy is called Jean’s swindle, and normally the background density is ignored (the two equations imply that a static universe is empty). The first order fluid equations are:

\begin{align} \frac{\partial \rho _{1}}{\partial t} + \rho _{0} \nabla \vec{v}_{1} & = 0 \\ \frac{\partial \vec{v}_{1}}{\partial t} & = -\frac{1}{\rho _{0}}\nabla P_{1} - \nabla \Phi _{1} \\ \nabla ^{2} \Phi _{1} & = 4\pi G \rho _{1} \end{align}\,\!

Taking the partial time derivative of the continuity equation and substituting in the momentum equation, and using the adiabatic sound speed and Poisson’s equation yields a wave equation with a driving term:

\begin{align} \frac{\partial ^{2}\rho _{1}}{\partial t^{2}} - v_{s}^{2}\nabla ^{2}\rho _{1} - 4\pi G\rho _{0}\rho _{1} & = 0 \end{align}\,\!

The solution to this equation is $\rho _{1}(\vec{r},t) = \int c(\vec{k})e^{i(\vec{k}\cdot \vec{r}-\omega t)} d^{3}k$. Plugging this in results in the dispersion relation for Jean’s instability:

\begin{align} \omega ^{2} & = v_{s}^{2}k^{2}-4\pi G \rho _{0} \end{align}\,\!

Jean’s wavenumber can be defined as $k_{J} = \sqrt {\frac{4\pi G\rho _{0}}{v_{s}^{2}}}$ so that $\omega ^{2} = v_{s}^{2}(k^{2} - k_{J}^{2})$. Therefore, if k > kJ, then ω is real and ρ1 oscillates (stable solution). If k < kJ, then ω is imaginary and $\rho _{1} \propto e^{\pm \omega t}$ (Jean’s instability).

Also note that Jean’s wavelength is $\lambda _{J} = \frac{2\pi }{k_{J}}$ and Jean’s mass is $M_{J} = \frac{4}{3}\pi \rho _{0}\lambda _{J}^{3}$. For λ > λJ, the system is unstable and will collapse.

There are two dimensionless analysis methods to obtain the relationship between the size of a collapsing cloud and other parameters. The first is by using the condition for collapse that Fgravity > Fpressure:

\begin{align} F_{g} & >& F_{P} \\ \frac{GMm}{R^{2}} & \sim P \cdot A \\ \frac{G\rho ^{2}}{(\frac{4}{3}\pi R^{3})^{2}R^{2}} & \sim v_{s}^{2}\rho 4\pi R^{2} \\ R & >& \frac{v_{s}}{\sqrt {G\rho }} \propto \lambda _{J} \end{align}\,\!

The second method is by comparing the gravity free-fall time scale with the sound crossing time (determined by pressure):

\begin{align} t_{g} & <& t_{P} \\ \frac{R}{v} & \sim \frac{R}{v_{s}} \\ \frac{R}{\sqrt {\frac{GM}{R}}} & \sim \frac{R}{v_{s}} \\ \frac{R^{3/2}}{\sqrt {G\rho \frac{4}{3}\pi R^{3}}} & \sim \frac{R}{v_{s}} \\ R & >& \frac{v_{s}}{\sqrt {G\rho }} \propto \lambda _{J} \end{align}\,\!

Jean’s instability is the basic reason why the Universe is not uniform. Structures have arisen from local over-densities of mass exceeding the local Jeans mass, and as the regions collapsed, sound waves did not have time to run ahead of the collapse and set up the pressure gradients required to offset it.

The case of perturbations in an expanding, 3D medium is relevant for cosmology and the expansion of the Universe. The perturbations are the following (the density perturbation is a fractional deviation from a uniform background, and the velocity perturbation is a peculiar velocity):

\begin{align} \vec{v} & = \vec{v}_{0} + \vec{v}_{1} \\ & = H\vec{r} + \vec{v}_{1} \\ & = \frac{\dot{a}}{a}\vec{r} + \vec{v}_{1} \\ P & = P_{0} + P_{1} \\ \rho & = \rho _{0} + \rho _{1} \\ & = \rho _{0}(t) + \delta \rho _{0}(t) \\ \Phi & = \Phi _{0} + \Phi _{1} \end{align}\,\!

The 0th order fluid equations are:

\begin{align} \frac{\partial \rho _{0}}{\partial t} + \nabla (\rho _{0}\vec{v}_{0}) & = 0 \\ \frac{\partial \vec{v}_{0}}{\partial t} + (\vec{v}_{0} \cdot \nabla )\vec{v}_{0} & = - \nabla \Phi _{0} \\ \nabla ^{2} \Phi _{0} & = 4\pi G \rho _{0} \end{align}\,\!

From the continuity equation, plugging in $\nabla \cdot \vec{v}_{0} = 3\frac{\dot{a}}{a}$ results in $\rho _{0} \propto a^{-3}$. The first order fluid equations are:

\begin{align} \frac{\partial \rho _{1}}{\partial t} + \rho _{0} \nabla \vec{v}_{1} +\nabla (\rho _{1}\vec{v}_{0}) & = 0 \\ \frac{\partial \vec{v}_{1}}{\partial t} + (\vec{v}_{0} \cdot \nabla )\vec{v}_{1} + (\vec{v}_{1} \cdot \nabla )\vec{v}_{0} & = -\frac{1}{\rho _{0}}\nabla P_{1} - \nabla \Phi _{1} \\ \nabla ^{2} \Phi _{1} & = 4\pi G \rho _{1} \end{align}\,\!

The continuity equation can be simplified by substituting ρ1 = δρ0 and the relationship found before that $\dot{\rho }_{0} = -3 \frac{\dot{a}}{a}\rho _{0}$. The momentum equation can also be simplified using sound speed and similar substitutions. The two equations become:

\begin{align} \dot{\delta } + \frac{\dot{a}}{a}\vec{r}\nabla \delta + \nabla \vec{v}_{1} & = 0 \\ \frac{\partial \vec{v}_{1}}{\partial t} + \frac{\dot{a}}{a}(\vec{r} \cdot \nabla )\vec{v}_{1} + \frac{\dot{a}}{a}\vec{v}_{1} & = -v_{s}^{2}\nabla \delta - \nabla \Phi _{1} \end{align}\,\!

It is useful to go into Fourier space to manipulate these equations. Also recall that the physical coordinate $\vec{r}$ is related to the comoving coordinate $\vec{x}$ by $\vec{r} = a(t) \vec{x}$.

\begin{align} [\delta (\vec{r},t), \vec{v}_{1}(\vec{r},t),\Phi _{1}(\vec{r},t)] & = \int d^{3}k e^{i\vec{k}\cdot \vec{r}/a(t)} [\delta _{k}(t),\vec{v}_{k}(t),\Phi _{k}(t)] \end{align}\,\!

In this Fourier space, the derivatives can be taken for δ, $\vec{v}_{1}$, and Φ1 with respect to t or r and be substituted into the three simplified fluid equations, which become:

\begin{align} \dot{\delta }_{k} + \frac{i\vec{k}}{a}\vec{v}_{k} & = 0 \\ \frac{\partial }{\partial t}(a\vec{v}_{k}) + iv_{s}^{2}\vec{k}\delta _{k} + i\vec{k}\Phi _{k} & = 0 \\ \frac{\vec{k}^{2}}{a^{2}}\Phi _{k} & = -4\pi G \rho _{0} \delta _{k} \end{align}\,\!

Finally, combining these three equations yields a wave equation with a damping term:

\begin{align} \ddot{\delta }_{k} + 2 \frac{\dot{a}}{a}\dot{\delta }_{k} + \frac{v_{s}^{2}}{a^{2}}(k^{2}-k_{J}^{2})\delta _{k} & = 0 \\ k_{J} & = \sqrt {\frac{4\pi G\rho _{0}a^{2}}{v_{s}^{2}}} \end{align}\,\!

Again, instability occurs when k < kJ.

Taking the case where the k term is ignored, so $\ddot{\delta }_{k} + 2 \frac{\dot{a}}{a}\dot{\delta }_{k} - 4\pi G\rho _{0}\delta _{k} = 0$, the perturbation solutions can be solved for in different types of universes. For a matter dominated Einstein-de-Sitter Universe, Ω = Ωm = 1, $4\pi G\rho = \frac{3}{2}H^{2}$, and $a \propto t^{2/3}$. Therefore, $\ddot{\delta }_{k} + \frac{4}{3t}\dot{\delta }_{k} - \frac{2}{3t^{2}}\delta _{k} = 0$. Taking an ansatz of $\delta _{k} \propto t^{\alpha }$ results in two solutions (growing and decaying):

\begin{align} \delta _{+} \propto t^{2/3} \propto a \\ \delta _{-} \propto t^{-1} \end{align}\,\!

#### Gravitational Instability in Rotating Disks: Uniform vs. Differential Rotation; Epicycle Frequency

Consider a uniform rotating thin sheet with constant angular velocity $\vec{\Omega } = \Omega \hat{z}$ with the following perturbations in the x-y plane (no 0th order velocity in the rotating frame):

\begin{align} \vec{v} & = \vec{v}_{1}(x,y,t) \\ P & = P_{0} + P_{1} \\ \Phi & = \Phi _{0} + \Phi _{1} \\ \Sigma (x,y,t) & = \Sigma _{0} + \Sigma _{1}(x,y,t) \end{align}\,\!

In this case, pressure has units of force/length. The density is ρ = Σδ(z), so $v_{s}^{2} = \frac{P_{1}}{\Sigma _{1}}$. The 0th order fluid equations are:

\begin{align} \frac{\partial \Sigma _{0}}{\partial t} + \nabla (\Sigma _{0}\vec{v}_{0}) & = 0 \\ \frac{\partial \vec{v}_{0}}{\partial t} + (\vec{v}_{0} \cdot \nabla )\vec{v}_{0} & = -\frac{1}{\Sigma _{0}}\nabla P_{0} - \nabla \Phi _{0} - 2\vec{\Omega } \times \vec{v}_{0} + \Omega ^{2}(x\hat{x}+y\hat{y}) \\ \nabla ^{2} \Phi _{0} & = 4\pi G \Sigma _{0}\delta (z) \end{align}\,\!

The last two terms in the momentum equation are contributions from the coriolis force and centrifugal force. The first order fluid equations are:

\begin{align} \frac{\partial \Sigma _{1}}{\partial t} + \Sigma _{0} \nabla \vec{v}_{1} & = 0 \\ \frac{\partial \vec{v}_{1}}{\partial t} & = - \frac{v_{s}^{2}}{\Sigma _{0}}\nabla \Sigma _{1} - \nabla \Phi _{1} - 2\vec{\Omega } \times \vec{v}_{1} \\ \nabla ^{2} \Phi _{1} & = 4\pi G \Sigma _{1}\delta (z) \end{align}\,\!

As for the expanding medium instability case, going into Fourier space is useful (only considering waves traveling in the x direction here).

\begin{align} [\Sigma _{1}(x,y,t),v_{1x}(x,y,t),\Phi _{1}(x,y,t)] & = \int dk e^{i(kx-\omega t)} [\Sigma _{1}(k), v_{1x}(k), \Phi _{1}(k)] \end{align}\,\!

The above transformation of Φ1 is valid for z = 0, but otherwise $\Phi _{1} \propto e^{i(kx-\omega t)}f(z)\Phi _{1}(k)$ and $\nabla ^{2}\Phi _{1} = 0$, so $f(z) \propto e^{-k|z|}$ (f needs to be finite for large z). Therefore, $\Phi _{1}(x,y,t) \propto e^{i(kx-\omega t)}e^{-k|z|}\Phi _{1}(k)$.

Integrating the first order Poisson equation from z = − ε to z = ε gives:

\begin{align} \frac{\partial \Phi _{1}}{\partial z}|_{z=\epsilon } - \frac{\partial \Phi _{1}}{\partial z}|_{z=-\epsilon } & = 4\pi G\Sigma _{1}(1) \\ -ke^{-k\epsilon }\Phi _{1} - ke^{-k\epsilon }\Phi _{1} & = 4\pi G\Sigma _{1} \\ \Phi _{1} = \frac{-2\pi G}{k}\Sigma _{1} \end{align}\,\!

Using Fourier Space, the first order continuity equation yields $v_{1x} = \frac{\omega }{k}\frac{\Sigma _{1}}{\Sigma _{0}}$. The momentum equation can be broken into x and y components, and it is found that $v_{1y} = -\frac{2\Omega i}{k}\frac{\Sigma _{1}}{\Sigma _{0}}$. Substituting relationships into the x-component yields the dispersion relation:

\begin{align} \omega ^{2} & = k^{2} v_{s}^{2}-2\pi Gk\Sigma _{0} + 4\Omega ^{2} \end{align}\,\!

There are three terms in this expression. The first is a pressure term that serves to stabilize the disk. The second is a gravity term that causes the disk to be unstable. The third is a rotation term that stabilizes.

Several comments can be made about this result. If there’s no rotation and negligible sound speed, the disk is always unstable (ω2 < 0). If there is no rotation (only gravity and pressure), ω2 < 0 for a Jeans wavenumber $k < k_{J} = \sqrt {\frac{2\pi G\Sigma _{0}}{v_{s}^{2}}}$ (long modes are unstable). If there is negligible sound speed (no pressure), then ω2 < 0 for $k > \frac{2\Omega ^{2}}{\pi G \Sigma _{0}}$ (short modes are unstable). This implies that on larger scales, rotation stabilizes gravity.

In general, the stability criterion is:

\begin{align} \frac{v_{s}\Omega }{G\Sigma _{0}} > \frac{\pi }{2} \end{align}\,\!

For a thin disk of collisional stars, vs becomes σ, the 1-dim velocity dispersion. Therefore, $\frac{\sigma \Omega }{G\Sigma _{0}} > 1.68$ for stability.

For a differentially rotating disk of gas, Ω = Ω(R) is not a constant so it needs to be replaced by the epicycle frequency.

\begin{align} \kappa ^{2} & = R \frac{\partial \Omega ^{2}}{\partial R} + 4\Omega ^{2} = \frac{2\Omega }{R}\frac{d}{dR}(R^{2}\Omega ) \end{align}\,\!

The form of κ depends on vcirc(R) = RΩ(R) or the mass distribution. For Keplerian orbits, κ = Ω. For flat rotation curves, $\kappa = \sqrt {2}\Omega$. And for solid-body rotation, κ = 2Ω. In the case of solid-body rotation, the disk is stable if $Q = \frac{v_{s}\kappa }{\pi G\Sigma _{0}} > 1$. This is called the Toomre Q parameter.

#### Rayleigh-Taylor Instability

When there are two layers of fluids under the influence of gravity, Rayleigh-Taylor instabilities result. This phenomenon is especially relevant for atmospheric and ocean sciences. Recalling the 0th order fluid equation solutions for water waves, the expressions for pressure for the top and bottom fluids in a two-layer system are:

\begin{align} P_{+,0} & = -\rho _{+,0}gz + C_{+} \\ P_{-,0} & = -\rho _{-,0}gz + C_{-} \end{align}\,\!

Also recalling the first order pressure solutions from the water waves case, $P_{+,1} \propto e^{-kz}$ for z > 0 and $P_{-,1} \propto e^{kz}$ for z < 0. The first order momentum equation is $\frac{\partial \vec{v}_{1}}{\partial t} = \frac{-1}{\rho _{0}}\nabla P_{1}$, and the solution for v1z has been derived previously. Therefore:

\begin{align} -i\omega v_{+,1z} & = \frac{kP_{+,1}}{\rho _{+,0}} \\ -i\omega v_{-,1z} & = \frac{-kP_{-,1}}{\rho _{-,0}} \end{align}\,\!

There are two boundary conditions to apply here. The first is that fluid displacements are continuous at z = 0. Remember that $\vec{v}_{1} = \frac{d\vec{\xi }}{dt}$, then $\xi _{z} = \frac{v_{1z}}{-i\omega }$ for both layers of fluids at z = 0.

The second boundary condition is that the force of the top fluid on the bottom fluid must equal the force of the bottom fluid on the top fluid. Without surface tension, this means that P + ,1L = P − ,1L. With surface tension, this means that $P_{+,1L} - \mathcal{T}\frac{\partial ^{2} \xi _{z}}{\partial x^{2}} = P_{-,1L}$. These both apply at z = 0. The left-hand side of the second equation is the total downward force per area. Recall also that $P_{1L} = P_{1} + \vec{\xi }\nabla P_{0}$, where P1L was set to zero for the water waves boundary condition (constant pressure along a fluid element). Applying this to the top fluid:

\begin{align} P_{+,1L} & = P_{+,1}+\xi _{z}\frac{\partial P_{+,0}}{\partial z} \\ & = \frac{-i \omega }{k}v_{+,1z}\rho _{+,0} + \frac{v_{1z}}{-i\omega }(-\rho _{+,0}g) \end{align}\,\!

Doing the same thing for the bottom fluid:

\begin{align} P_{-,1L} & = P_{-,1}+\xi _{z}\frac{\partial P_{-,0}}{\partial z} \\ & = \frac{i \omega }{k}v_{-,1z}\rho _{-,0} + \frac{v_{1z}}{-i\omega }(-\rho _{-,0}g) \end{align}\,\!

The two equations are related by the surface tension term.

\begin{align} \frac{-i \omega }{k}v_{+,1z}\rho _{+,0} + \frac{v_{1z}}{-i\omega }(-\rho _{+,0}g) - \mathcal{T}\frac{\partial ^{2} \xi _{z}}{\partial x^{2}} & = \frac{i \omega }{k}v_{-,1z}\rho _{-,0} + \frac{v_{1z}}{-i\omega }(-\rho _{-,0}g) \end{align}\,\!

Simplifying this expression and noting that v + ,1z = v − ,1z because of the first boundary condition produces the dispersion relation:

\begin{align} \omega ^{2}(\rho _{+,0}+\rho _{-,0}) & = gk(\rho _{-,0}-\rho _{+,0})+k^{3}\mathcal{T} \end{align}\,\!

For ρ − ,0 > ρ + ,0 (denser fluid on the bottom), ω2 > 0 and the system is stable (perturbations oscillate and travel as ei(kx − ωt) waves. If $\rho _{-,0} \gg \rho _{+,0}$ (like for water waves), then $\omega ^{2} \sim gk + \frac{T}{\rho _{-,0}}k^{3}$ and the solution is always stable.

For ρ + ,0 > ρ − ,0 (denser fluid on the top), ω2 < 0 when $gk(\rho _{-,0}-\rho _{+,0}) + k^{3}\mathcal{T} < 0$, or $k^{2} < k_{crit}^{2} = \frac{g(\rho _{+,0}-\rho _{-,0})}{\mathcal{T}}$. Therefore, instability arises for long modes, and increased surface tension helps to stabilize the perturbations.

For the case of negligible surface tension, $\omega ^{2} = \frac{gk(\rho _{-,0}-\rho _{+,0})}{\rho _{-,0}+\rho _{+,0}}$, and ω2 < 0 when ρ − ,0 < ρ + ,0.

And finally, note that the downward gravitational acceleration $\vec{g}$ can be mapped into an upward $\vec{a}$ acceleration. Therefore, a light fluid accelerating into a heavy fluid also gives rise to the same instability as a denser fluid on top of a less dense fluid.

An application of a Rayleigh-Taylor instability is a supernova explosion. For a blast wave, a dense shell of gas is on top of less dense gas. This gives rise to an instability that leads to the production of filaments in the post-shock gas.

#### Kelvin-Helmholtz Instability

Kelvin-Helmholtz instabilities are similar to Rayleigh-Taylor instabilities (two layers of different densities), except that there is an additional velocity shear between the layers. This type of instability arises in clouds, the Sun’s corona, and with wind blowing over water.

The derivation of the dispersion relation for this type of instability is similar to that for Rayleigh-Taylor, except that $\vec{v}_{+,0} \ne 0$ and $\vec{v}_{-,0} \ne 0$. The resulting dispersion relation is:

\begin{align} (\omega - kv_{+,0})^{2}\rho _{+,0} + (\omega - kv_{-,0})^{2}\rho _{-,0} & = gk(\rho _{-,0}-\rho _{+,0})+k^{3} \mathcal{T} \end{align}\,\!

This can be rewritten as a quadratic equation for ω / k. When surface tension is negligible, instability arises for short wavelength modes, or:

\begin{align} k & >& \frac{(\rho _{-,0}^{2} - \rho _{+,0}^{2})g}{\rho _{-,0}\rho _{+,0}(v_{+,0}-v_{-,0})^{2}} = k_{crit} \end{align}\,\!

Note that if surface tension is present, the short modes can be stabilized. Also, if g = 0, any wavenumber is Kelvin-Helmholtz unstable. If kcrit is small, Kelvin-Helmholtz instabilities are common, and if kcrit is large, only very small wavelengths are unstable.

### IV. Sticky Stuff

#### Viscous Stress Tensor and Force; Heuristic and Mathematical Derivations

Recall the conservative form of the momentum equation, $\frac{\partial }{\partial t}(\rho v_{i}) = -\partial _{j}T_{ij}-\rho \partial _{i}\Phi$. This equation was based on the premise that in the fluid element’s frame, its momentum does not change as a result of the motion of neighboring fluid elements. Rather, it is the result of momentum flowing with the fluid and thermal pressure gradients. This equation, however, needs to be modified if there’s a process that can transfer momentum associated with velocity differences between elements. This type of process is called a viscous process (a measure of fluid ‘friction’ or resistance to being deformed by shear stress), and an extra term representing the viscous stress tensor must be added onto the expression for the momentum flux stress tensor.

\begin{align} T_{ij} & = \rho v_{i}v_{j}+P\delta _{ij} - \sigma _{ij} \end{align}\,\!

Consider a shear viscosity in a linear shear flow (velocity gradient in the x-direction is perpendicular to the y-direction of the streamlines). If different streamlines can communicate, then they transfer momentum to each other and this net momentum flux makes up the new term above. The momentum transferred per particle across a boundary at x (from particles at x + l / 2 and xl / 2) is $p_{y}(x) \sim m \frac{\partial v_{y}}{\partial x} l$. The number of particles per unit time passing the boundary at x is N˜nAvth, where $v_{th} \sim \sqrt {\frac{kT}{m}}$. Therefore, $\dot{p}_{y}(x) \sim (m \frac{\partial v_{y}}{\partial x} l)(nAv_{th})$. The viscous force on a volume element AΔx is $\frac{\partial \dot{p}_{y}}{\partial x}\Delta x$, and therefore the viscous force per unit volume is $\frac{F_{y}}{A\Delta x} \sim \frac{\partial }{\partial x}(\rho \nu \frac{\partial v_{y}}{\partial x})$, where ρ = mn and ν = lvth. ν is called the kinematic viscosity. The dynamic viscosity is η = ρν, and the viscous shear stress tensor is $\sigma _{yx} = \eta \frac{\partial v_{y}}{\partial x}$. Finally, the force per unit volume is the spatial gradient of the viscous stress tensor. Note that viscosity is the first spatial gradient of velocity.

A mathematical derivation of the stress tensor has the general form $\sigma _{ij} = \eta (\partial _{j}v_{i}+\alpha \partial _{i}v_{j}+\beta \delta _{ij}\partial _{k}v_{k})$. Alpha and beta can be determined from the conditions that σij = 0 for uniform solid-body rotation (α = 1) and σij = 0 for uniform isotropic expansion (such as the Hubble Flow) (β = − 2 / 3). Putting things together:

\begin{align} \sigma _{ij} & = \eta (\partial _{j}v_{i}+\partial _{i}v_{j}-\frac{2}{3} \delta _{ij}\partial _{k}v_{k}) \end{align}\,\!

#### Navier-Stokes Equation; Reynolds Number

Taking the derivative of the viscous tensor yields $\eta (\nabla ^{2}v_{i}+\frac{1}{3}\partial _{i}(\nabla \cdot \vec{v}))$. Putting things into the momentum equation produces the Navier-Stokes equation:

\begin{align} \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla )\vec{v} & = -\frac{\nabla P}{\rho } - \nabla \Phi + \nu (\nabla ^{2}\vec{v} + \frac{1}{3}\nabla \cdot (\nabla \cdot \vec{v})) \end{align}\,\!

The last parenthetical term in the Navier-Stokes equation is the viscous force per unit mass. Note that viscosity shows up as a diffusion term in the momentum equation, and this kinematic viscosity has units of [cm2 / s] (0.01 for water, 0.15 for air, 1 for oil, etc.).

The viscous stress tensor is symmetric and traceless. If there is bulk viscosity (ex: shocks), then diagonal terms are added: $\sigma _{ij} \rightarrow \sigma _{ij} + \gamma \delta _{ij}\partial _{k}v_{k}$. The extra term is associated with isotropic motion (compression or expansion).

The Reynolds number is a dimensionless quantity and is the ratio of inertial forces ($\vec{v} \cdot \nabla \vec{v}$) to viscous forces ($\nu \nabla ^{2} \vec{v}$).

\begin{align} \mathcal{R} & = \frac{vL}{\nu } \end{align}\,\!

For $\mathcal{R} \ll 1$, viscosity is important, and for $\mathcal{R} \gg 1$, the flow is inviscid. For very high Reynolds number ($\mathcal{R} > 3000$), flows tend to be turbulent, while low Reynolds number flows are laminar.

#### Applications: Poiseuille Flow; Stokes Flow

An application of viscous flows is the Poiseuille flow, which is a laminar, constant flow of an incompressible uniform viscous fluid through a tube of constant cross section (no gravity). The Navier-Stokes equation becomes:

\begin{align} (\vec{v} \cdot \nabla )\vec{v} & = -\frac{\nabla P}{\rho } + \nu \nabla ^{2}\vec{v} \end{align}\,\!

If the velocity is taken to be constant along the flow in the z-direction, then $\frac{\partial v}{\partial z} = 0$ but $\frac{\partial v}{\partial x} \ne 0$ and $\frac{\partial v}{\partial y} \ne 0$. In the z-direction (for a pressure difference over L):

\begin{align} 0 & = -\frac{1}{\eta }\frac{\partial P}{\partial z} + (\partial _{x}^{2}+\partial _{y}^{2})v \\ (\partial _{x}^{2}+\partial _{y}^{2})v & = -\frac{1}{\eta }\frac{\Delta P}{\Delta L} \end{align}\,\!

Using cylindrical coordinates yields a velocity solution $v(R) = -\frac{1}{\eta }\frac{\Delta P}{\Delta L}(\frac{R^{2}}{4} + a ln R + b)$, where a and b are integration constants. The boundary condition of v(R = 0) = finite gives a = 0 and the boundary condition v(R = R0) = 0 gives the Poiseuille flow profile:

\begin{align} v(R) & = \frac{1}{4\eta }\frac{\Delta P}{\Delta L}(R_{0}^{2} - R^{2}) \end{align}\,\!

Some comments can be made about this flow profile. The maximal flow occurs at the tube’s center, and the wider the tube, the faster the flow. The fluid at the boundary of the tube communicates with the fluid at the center through the viscous force. Finally, the mass flux through the tube can be found:

\begin{align} \dot{M} & = \int v \rho dA \\ & = \int v \rho 2\pi R dR \\ & = \frac{\pi \Delta P}{8\nu \Delta L}R_{0}^{4} \end{align}\,\!

Another example of a viscous flow is Stoke’s flow. This is a ‘creeping flow’ around a sphere of radius a moving through a viscous fluid. For a constant, low-Reynolds number viscous flow of an incompressible fluid, the Navier-Stokes equation becomes $0 = -\nabla P - \eta \nabla ^{2} \vec{v}$. Taking the divergence and curl of this expression gives $\nabla ^{2} P = 0$ and $\nabla ^{2} \vec{w} = 0$, respectively. The boundary condition of this problem is $\vec{v} = 0$ at r = a and the general solution for velocity is $\vec{v} = V(A(r)cos \theta \hat{r} - B(r) sin \theta \hat{\theta })$. Doing a lot of math (using the incompressibility condition to solve for B in terms of A, solving for the φ-component of vorticity, taking $\nabla ^{2} w = 0$ to find the form of A and using the boundary condition to find B and the r and θ components of velocity, solving for pressure using the Navier-Stokes equation, etc.) finally yields an expression for the z-component of the force exerted by the fluid on the sphere. Integrating this over the sphere’s surface ($\int F_{z} a^{2} sin \theta d\theta d\phi$) gives the drag force on the sphere.

\begin{align} F_{drag} & = 6\pi \eta v a \end{align}\,\!

The terminal speed of the sphere is reached when the drag force and buoyancy force are equal to the gravitational force.

\begin{align} F_{d} + F_{b} & = F_{g} \\ 6\pi \eta v_{terminal} a + \rho _{f}g\frac{4}{3}\pi a^{3} & = \rho _{s}g\frac{4}{3}\pi a^{3} \\ v_{terminal} & = \frac{2}{9}\frac{ga^{2}}{\nu }(\frac{\rho _{s}}{\rho _{f}}-1) \end{align}\,\!

#### Accretion Disk: Molecular vs. Turbulent Viscosity

Consider an axisymmetric thin disk ($\frac{\partial }{\partial \phi } = 0$) that is rotating with angular velocity $\vec{\Omega } = \Omega \hat{z}$. For accretion disks, the quantities of interest are the orbital motion (vφ) and radial motion due to viscosity (vr). The continuity equation in cylindrical coordinates becomes:

\begin{align} \frac{\partial \Sigma }{\partial t} + \nabla (\Sigma \vec{v}) & = 0 \\ \frac{\partial \Sigma }{\partial t} + \frac{1}{R}\frac{\partial }{\partial R}(\Sigma R v_{R}) & = 0 \end{align}\,\!

Across some ΔR of the disk, the change in mass is ΔM = ρ2πRΔRdz. Therefore, the inward mass accretion rate is $\dot{M} = -\int \frac{\Delta M}{\Delta t}$. Noting that $\Sigma = \int \rho dz$, the mass accretion rate becomes $\dot{M} = -2\pi R \Sigma v_{R}$.

From the momentum (Navier-Stokes) equation, where $\vec{F}$ is any other relevant forces, the φ component can be looked at.

\begin{align} \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla )\vec{v} & = -\frac{\nabla P}{\rho } + (\nabla \nu \nabla ) \vec{v} + \vec{F} \\ \frac{\partial v_{\phi }}{\partial t} + (\vec{v} \cdot \nabla )v_{\phi } & = \frac{\partial v_{\phi }}{\partial t} + v_{R}\frac{\partial v_{\phi }}{\partial R} + \frac{v_{R}v_{\phi }}{R} \\ & = \frac{\partial v_{\phi }}{\partial t} + \frac{v_{R}}{R}\frac{\partial }{\partial R}(R^{2} \Omega ) \\ \Sigma \Big[ \frac{\partial v_{\phi }}{\partial t} + \frac{v_{R}}{R}\frac{\partial }{\partial R}(R^{2} \Omega )\Big] & = (\nabla \Sigma \nu \nabla )v_{\phi } \\ \Sigma \Big[ \frac{\partial v_{\phi }}{\partial t} + \frac{v_{R}}{R}\frac{\partial }{\partial R}(R^{2} \Omega )\Big] & = \frac{\partial }{\partial R}(\nu \Sigma \frac{\partial v_{\phi }}{\partial R}) + \frac{1}{R}\frac{\partial }{\partial R}(\nu \Sigma v_{\phi }) - \frac{\nu \Sigma v_{\phi }}{R^{2}} \end{align}\,\!

The relationship vφ = RΩ was used above, as well as $\frac{\partial v_{\phi }}{\partial \phi } = v_{R}$. All derivatives with respect to z are zero. Manipulating the continuity and momentum equations results in:

\begin{align} \frac{\partial }{\partial t}(\Sigma Rv_{\phi }) + \frac{1}{R}\frac{\partial }{\partial R}(\Sigma R^{3}\Omega v_{R}) = \frac{1}{R}\frac{\partial }{\partial R}(\nu \Sigma R^{3}\frac{\partial \Omega }{\partial R}) \end{align}\,\!

This equation represents the net viscous torque exerted on the annulus. The first term represents the Eulerian rate of change of angular momentum (per unit area), and the second term represents the net rate of angular momentum loss from the area due to the flow of angular momentum in the radial direction. The angular momentum content of this region is evolving as a result of the torque. Defining the torque as $T = 2\pi \nu \Sigma R^{3} \frac{\partial \Omega }{\partial R}$, then the right hand side of the equation is $\frac{1}{2\pi R}\frac{\partial T}{\partial R}$. Recall also that the torque is the viscous force times R, where the viscous force is the product of the area and the viscous stress.

Assuming $\dot{\Omega } = 0$ and a Keplerian disk ($v_{\phi } = R\Omega = \sqrt {\frac{GM}{R}}$), combining the continuity and angular momentum conservation equations gives:

\begin{align} \frac{\partial \Sigma }{\partial t} & = \frac{3}{R}\frac{\partial }{\partial R}(R^{1/2}\frac{\partial }{\partial R}(\nu \Sigma R^{1/2})) \end{align}\,\!

This is the basic equation for the time evolution of surface density in a Keplerian accretion disk. One solution to this is the case of constant ν. Changing variables from $s = 2\sqrt {R}$ yields an expression that can be solved for with separation of variables. The solution is a combination of an exponential function and a bessel function. Using dimensionless variables $x = \frac{R}{R_{0}}$ and $\tau = \frac{12\nu }{R_{0}^{2}}t$, then the surface density can go from being a function of R and t to x and τ.

\begin{align} \Sigma (x,\tau ) & = \frac{m}{\pi R_{0}^{2}}\frac{1}{\tau x^{1/4}}e^{\frac{-1+x^{2}}{\tau }}I_{1/4}\Big(\frac{2x}{\tau }\Big) \end{align}\,\!

The main take-away here is that viscosity spreads out a ring. Greater values of τ lowers the surface density. Most of the mass moves inwards, losing energy and angular momentum, but a small amount moves outwards to conserve angular momentum. From the continuity equation and surface density evolution equation:

\begin{align} v_{R} & = -3 \nu \frac{\partial }{\partial R}ln(\Sigma R^{1/2}) \end{align}\,\!

And using the dimensionless variables, this becomes:

\begin{align} v_{R} & = -\frac{3\nu }{R_{0}}\frac{\partial }{\partial x}(\frac{1}{4}ln x - \frac{1+x^{2}}{\tau }+ln I_{1/4}\Big(\frac{2x}{\tau }\Big)) \end{align}\,\!

The behavior of I1 / 4(z) is flat for low z’s and steeper for high z’s. Therefore, for 2x > τ:

\begin{align} v_{R} & \sim \frac{3\nu }{R_{0}}\Big(\frac{1}{4x} + \frac{2x}{\tau } - \frac{2}{\tau }\Big) \end{align}\,\!

And for 2x < τ:

\begin{align} v_{R} & \sim -\frac{3\nu }{R_{0}}\Big(\frac{1}{2x}-\frac{2x}{\tau }\Big) \end{align}\,\!

For x > 1, R > R0, and vR > 0. Therefore, the outer parts move outward. For τ > 4x2, vR < 0, and the inner part moves inward. At very long times, most of the mass will have moved inward, but the mass moving outward takes the angular momentum.

What is the source of viscosity in accretion disks? The inward velocity is $v_{R} \sim \frac{\nu }{R}$, so the timescale for orbital evolution is $t_{inflow} \sim \frac{R}{v_{R}} \sim \frac{R^{2}}{\nu } \sim \frac{R}{v_{\phi }}\frac{Rv_{\phi }}{\nu } \sim t_{orbital}\mathcal{R}$. Therefore, the timescale for radial evolution due to viscosity exceeds the local orbital timescale by a factor on the order of the Reynolds number. The Reynolds number is generally very large, since $\nu \sim v_{thermal} \lambda _{mfp} \sim v_{thermal}\frac{1}{\sigma n}$, and typical numbers give $\mathcal{R} \sim 10^{10}$. Molecular viscosity is very slow!

In order to avoid such a large Reynolds number, the viscosity must be high. If turbulence causes high viscosity, then ν˜vturblturb˜vsH, where H is the disk thickness. Disks can therefore be characterized as ν = αvsH, where α parametrizes ignorance. This characterization is known as the Shakura-Sunyaev disk.

### V. Kinetic Theory

#### Phase Space; Boltzmann Equation

Kinetic theory is a way of describing a fluid as a collection of particles in constant, random motion. It deals with a statistical distribution of a ‘gas’ made up of particles that travel and collide. The basic quantity is a phase space distribution: $f(\vec{x}, \vec{p}, t)$. The number of particles is $dN = f(\vec{x}, \vec{p}, t)d^{3}xd^{3}p$. A simple example of a phase space distribution is:

\begin{align} f(p) = \frac{g}{h^{3}}\frac{1}{\rho ^{E/kT}\pm 1} \end{align}\,\!

The continuity fluid equation equivalent in kinetic theory is the collisionless Boltzmann equation.

\begin{align} \frac{\partial f}{\partial t} + \frac{\partial }{\partial x_{i}}(\dot{x}_{i}f) + \frac{\partial }{\partial \rho _{i}}(\dot{\rho }_{i}f) & = 0 \end{align}\,\!

The second term is the flux out of a spatial volume element, and the third term is the flux out of a momentum volume element. The zero on the right-hand side means that it is collisionless. Expanding the derivatives out and recalling Hamilton’s definitions ($\dot{x}_{i} = \frac{\partial H}{\partial p_{i}}$ and $\dot{p}_{i} = -\frac{\partial H}{\partial x_{i}}$), the equation becomes:

\begin{align} \frac{\partial f}{\partial t} + \dot{x}_{i}\frac{\partial }{\partial x_{i}}f + \dot{p}_{i}\frac{\partial }{\partial p_{i}}f & = 0 \\ & = \frac{Df}{Dt} \end{align}\,\!

#### Similarities to and Differences from Fluids: Infinite Hierarchy, Closure Relations

If the fluid is collisional, then $\frac{Df}{Dt} \ne 0$. To proceed, moments of f are the following:

\begin{align} n & = \int f d^{3} p \\ \langle \vec{v}\rangle & = \frac{\int \vec{v} f d^{3} p}{\int f d^{3} p} \\ \langle v_{i}v_{j}\rangle & = \frac{\int v_{i}v_{j} f d^{3} p}{\int f d^{3} p} \end{align}\,\!

Note that $\sigma _{ij}^{2} = \langle v_{i}v_{j}\rangle - \langle v_{i}\rangle \langle v_{j}\rangle$. From the moments, macroscopic properties of the fluid can be determined. The moments of the Boltzmann equation are taken using the relationships found previously (and $\dot{p} = -m \nabla \Phi$) to yield the continuity equation and Euler equation:

\begin{align} \frac{\partial n}{\partial t} + \nabla (n\langle \vec{v}\rangle ) & = 0 \\ \frac{\partial }{\partial t}\langle v_{j}\rangle + \langle v_{i} \rangle \partial _{i}\langle v_{j} \rangle & = -\partial _{j} \Phi - \frac{1}{n}\partial _{i}(n\sigma _{ij}^{2}) \end{align}\,\!

Notice that most of the momentum equation is recovered. However, there is no equation of state in kinetic theory. Instead, the 2nd order velocity moment is used ($\sigma _{ij}^{2} \rightarrow \frac{P}{\rho } \propto T$). The second moment of the Boltzmann equation is the heat equation.

\begin{align} \frac{3}{2} \frac{\partial T}{\partial t} + (\nabla \cdot \vec{v})T + \frac{1}{\rho }\nabla \cdot \vec{q} & = 0 \end{align}\,\!

The heat flux is $q_{i} = \frac{1}{2} \rho \sum \limits _{j=1}^{3} \sigma _{ijj}$.

#### Collisionless vs. Collisional Matter: CMB Physics

Recall the 3D gravitational instability in an expanding universe. There were three linearized equations (continuity, momentum, Poisson) in co-moving coordinates and in k-space.

\begin{align} \dot{\delta }_{k} + \frac{i\vec{k}}{a}\vec{v}_{k} & = 0 \\ \frac{\partial }{\partial t}(a\vec{v}_{k}) + iv_{s}^{2}\vec{k}\delta _{k} + i\vec{k}\Phi _{k} & = 0 \\ \frac{\vec{k}^{2}}{a^{2}}\Phi _{k} & = -4\pi G \rho _{0} \delta _{k} \end{align}\,\!

Defining $\theta = \nabla \cdot \vec{v}$, then θk = ikvk.

\begin{align} \dot{\delta }_{k} + \frac{\theta _{k}}{a} & = 0 \\ \dot{\theta }_{k} + \frac{\dot{a}}{a}\theta _{k} - v_{s}^{2}k^{2}\frac{1}{a}\delta _{k} - \frac{k^{2}}{a}\Phi _{k} & = 0 \\ k^{2} \Phi _{k} & = -4\pi G a^{2} \rho _{0} \delta _{k} \end{align}\,\!

The full linear perturbation theory results in similar expressions with a few extra terms. Extra moments need to keep getting taken to try to close the system. These equations can be applied to baryons, cold dark matter, photons, etc.

### VI. MagnetoHydroDynamics

#### Motivation; Maxwell’s Equations

In the presence of magnetic fields, current can be induced in electrically conducting fluids, which creates forces on the fluid. In other words, if a fluid is composed of charged particles, their behavior can be modified by electromagnetic fields. Their velocity component parallel to the magnetic field is unaffected, but the velocity component perpendicular gives rise to a force that leaves their perpendicular motion constrained. Some examples of electrically conducting fluids are liquid metals (mercury, sodium, etc.) and plasma (hot ionized gas). The equations governing MHD are a combination of the Navier-Stokes fluid equation and Maxwell’s equations.

Maxwell’s equations in CGS units are:

\begin{align} \nabla \cdot \vec{E} & = 4\pi \rho _{e} \\ \nabla \cdot \vec{B} & = 0 \\ \nabla \times \vec{E} & = -\frac{1}{c}\frac{\partial \vec{B}}{\partial t} \\ \nabla \times \vec{B} & = \frac{4\pi }{c}\vec{J} + \frac{1}{c}\frac{\partial \vec{E}}{\partial t} \end{align}\,\!

Recall that $\vec{J}$ is the charge current density with units of charge per area per time. ρ is the charge density. Combining the first and fourth equations give a conservation law similar to the continuity equation:

\begin{align} \frac{\partial \rho _{e}}{\partial t} + \nabla \cdot \vec{J} & = 0 \end{align}\,\!

For steady currents, $\frac{\partial \rho _{e}}{\partial t} = 0$, so Ampere’s law for steady currents is $\nabla \cdot \vec{J} = 0$.

In the absence of $\vec{J}$, taking the cross product of the fourth equation and substituting in the third equation gives a waves equation:

\begin{align} \Big(\nabla ^{2} - \frac{1}{c^{2}}\frac{\partial ^{2}}{\partial t^{2}}\Big)\vec{B} & = 0 \end{align}\,\!

#### Ohm’s Law; Induction Equation; Magnetic Reynolds Number

Ohm’s law is an empirical relationship between $\vec{J}$ and $\vec{E}$. Note that it is not always linear (it fails at strong voltages), but in general, in the rest frame of the conductor: $\vec{J}’ = \sigma E’$. The electric conductivity is σ. Since J˜I / A and E˜V / L, then $V \sim \frac{L}{A\sigma }I = IR$. In the lab frame, the conductor moves with some velocity, so Ohm’s law becomes $\vec{J} = \sigma (\vec{E} + \frac{\vec{v}}{c} \times \vec{B})$. In astronomical systems, σ is usually very large (perfect conductor). It’s noted that in this case, $|\vec{E}| \sim \frac{v}{c}|\vec{B}| \ll |\vec{B}|$, so the right-most term in Ampere’s law (fourth Maxwell equation) is typically small (flow is highly non-relativistic).

Ohm’s, Faraday’s, and Ampere’s equations can be combined (using the above approximation to ignore the $\vec{E}$ term in Ampere’s law) into a single equation for $\vec{B}$.

\begin{align} \frac{\partial \vec{B}}{\partial t} & = \nabla \times (\vec{v} \times \vec{B}) + \eta \nabla ^{2} \vec{B} \end{align}\,\!

This is the induction equation that describes the evolution of the magnetic field. The magnetic diffusivity is $\eta = \frac{c^{2}}{4\pi \sigma }$.

The magnetic Reynolds number is the ratio of the two terms on the right-hand side of the induction equation. The term with η is a diffusion term.

\begin{align} \mathcal{R}_{m} \sim \frac{vL}{\eta } \end{align}\,\!

At low $\mathcal{R}_{m}$, the diffusion term dominates and $\frac{\partial \vec{B}}{\partial t} \sim \eta \nabla ^{2}B$ so $\tau \sim \frac{L^{2}}{\eta }$ (field decay timescale). For high $\mathcal{R}_{m}$, the convective term dominates (remember that this case is still for high σ), and $\frac{\partial \vec{B}}{\partial t} \sim \nabla \times (\vec{v} \times \vec{B})$. In astronomical systems, L is very large, and both τ and $\mathcal{R}_{m}$ are much greater than 1.

For the case of high $\mathcal{R}_{m}$, the equation is similar to the Helmholtz equation $\frac{\partial \vec{w}}{\partial t} = \nabla \times (\vec{v} \times \vec{w})$, where $\vec{w}$ is vorticity. This equation states that the flux of vorticity is conserved. Therefore, for high $\mathcal{R}_{m}$, the flux of $\vec{B}$ is conserved (the magnetic field moves with the fluid, or is ‘frozen in’ to a perfectly conducting fluid).

The Lorentz force on a charge q moving at speed $\vec{u}$ is $\vec{F} = q\vec{E} + \frac{q}{c}\vec{u} \times \vec{B}$. Therefore, the force density is $\vec{F} = \rho _{e} \vec{E} + \frac{\rho _{e}}{c}\vec{u} \times \vec{B}$. Assuming the first term is typically small and using $\vec{J} = \rho _{e}\vec{u}$:

\begin{align} \vec{F} & \sim \frac{1}{c}\vec{J} \times \vec{B} \\ & \sim \frac{1}{c}\frac{c}{4\pi }(\nabla \times \vec{B}) \times \vec{B} \\ & \sim \frac{1}{4\pi }(-\frac{1}{2}\nabla B^{2}+(\vec{B} \cdot \nabla )\vec{B}) \end{align}\,\!

This is the force density (force per volume) contribution from the magnetic field. Thus, the Navier-Stokes equation (force per mass) becomes:

\begin{align} \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla )\vec{v} & = -\frac{1}{\rho }\nabla \Big(P + \frac{B^{2}}{8\pi }\Big) - \nabla \Phi + \frac{1}{4\pi \rho }(\vec{B} \cdot \nabla )\vec{B} + \nu (\nabla ^{2}\vec{v} + \frac{1}{3}\nabla \cdot (\nabla \cdot \vec{v})) \end{align}\,\!

The $\frac{B^{2}}{8\pi }$ term is a magnetic pressure term that arises from squeezing magnetic fields together (due to a non-uniform magnetic field). The other magnetic term is a magnetic tension term arising from a gradient along $\vec{B}$ (wanting to straighten out the curved field lines).

A measure of field strength is the ratio between thermal pressure and magnetic pressure, or $\beta \sim \frac{P}{\frac{B^{2}}{8\pi }}$. If $\beta \gg 1$, it is a weak field, and if $\beta \ll 1$, it is a strong field.

#### Ideal MHD; Magnetohydrodynamic Waves; Alfven Speed

Consider a non-viscous, perfectly conducting fluid ($\sigma \rightarrow \infty$) in a $\vec{B}$ field in the absence of gravity. The linear perturbations are the following:

\begin{align} \rho (\vec{x},t) = \rho _{0} + \rho _{1}(\vec{x},t) \\ \vec{v}(\vec{x},t) = v_{1}(\vec{x},t) \\ \vec{B}(\vec{x},t) = \vec{B}_{0} + \vec{B}_{1}(\vec{x},t) \end{align}\,\!

The Euler, Navier-Stokes, and induction first order equations are (for infinite conductivity):

\begin{align} \frac{\partial \rho _{1}}{\partial t} + \rho _{0} \nabla \cdot \vec{v}_{1} & = 0 \\ \frac{\partial \vec{v}_{1}}{\partial t} & = -\frac{v_{s}^{2}}{\rho _{0}}\nabla \rho _{1} + \frac{1}{4\pi \rho _{0}}(\nabla \times \vec{B_{1}}) \times \vec{B_{0}} \\ \frac{\partial \vec{B}_{1}}{\partial t} & = \nabla \times (\vec{v}_{1} \times \vec{B}_{0}) \end{align}\,\!

Taking the partial time derivative of the second equation and substituting in the first and third gives the MHD equation:

\begin{align} \frac{\partial ^{2}\vec{v}_{1}}{\partial t^{2}} - v_{s}^{2}\nabla (\nabla \cdot \vec{v}_{1}) + \vec{v}_{A} \times (\nabla \times (\nabla \times (\vec{v}_{1} \times \vec{v}_{A}))) & = 0 \end{align}\,\!

The Alfven velocity is defined as $\vec{v}_{A} = \frac{1}{\sqrt {4\pi \rho _{0}}}\vec{B}_{0}$

Going into Fourier space, where $\vec{v}_{1} = \int d^{3}k \vec{v}_{1}(\vec{k})e^{i(\vec{k} \cdot \vec{x} -\omega t)}$, the wave dispersion equation is found:

\begin{align} -\omega ^{2}\vec{v}_{1} + (v_{s}^{2}+v_{A}^{2})(\vec{k} \cdot \vec{v}_{1})\vec{k} + \vec{v}_{A} \cdot \vec{k}[(\vec{v}_{A} \cdot \vec{k})\vec{v}_{1} - (\vec{v}_{1} \cdot \vec{k})\vec{v}_{A} - (\vec{v}_{A} \cdot \vec{v}_{1})\vec{k}] & = 0 \end{align}\,\!

For the case of $\vec{k} \perp \vec{v}_{A}$ (or $\vec{k} \perp \vec{B}_{0}$), looking at the $\vec{v}_{1}\hat{k}$ component:

\begin{align} \omega ^{2} \vec{v}_{1} & = (v_{s}^{2}+v_{A}^{2})(\vec{k} \cdot \vec{v}_{1})\vec{k} \\ & = (v_{s}^{2}+v_{A}^{2}) k^{2}(\hat{k} \cdot \vec{v}_{1})\hat{k} \\ \omega ^{2} & = (v_{s}^{2}+v_{A}^{2})k^{2} \\ v_{\phi } & = \sqrt {v_{s}^{2}+v_{A}^{2}} \end{align}\,\!

The phase velocity depends on both the hydrostatic and magnetic pressure. Any longitudinal disturbances result in rarefactions and compressions, like an ordinary sound wave. However, the magnetic field lines are bunched together in the compressions since it’s frozen to the fluid, so there’s an additional magnetic pressure resisting the compression.

For the case of $\vec{k} \parallel \vec{v_{A}} \parallel \vec{B}_{0}$, there are two types of wave motion. The first is longitudinal ($\vec{v}_{1} \parallel \vec{k}$):

\begin{align} -\omega ^{2}v_{1} + (v_{s}^{2}+v_{A}^{2})k^{2}v_{a} - v_{A}^{2}k^{2}v_{1} & = 0 \\ \omega ^{2} & = v_{s}^{2}k^{2} \end{align}\,\!

These are normal sound waves. The fluid does not ‘know’ about the compressions and rarefactions. The second type of wave motion is transverse ($\vec{v}_{1} \perp \vec{k}$):

\begin{align} -\omega ^{2}v_{1} + v_{A}k(v_{A}kv_{1}) & = 0 \\ \omega ^{2} & = v_{A}^{2}k^{2} \end{align}\,\!

These are Alfven waves (pure MHD waves). Magnetic tension provides a restoring force on displaced fluid elements. The motion of the fluid elements and perturbations are perpendicular to the wave propagation, and when the perturbation crosses the magnetic field line, it ‘plucks’ it and Alfven waves travel along the magnetic field lines.

#### Magnetorotational Instabilities

Magnetorotational instabilities are fluid instabilities that cause an accretion disk to become turbulent. Consider a non-viscous accretion disk in a magnetic field with cylindrical coordinates. The Navier-Stokes equation says that:

\begin{align} \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla )\vec{v} & = -\frac{1}{\rho }\nabla \Big(P + \frac{B^{2}}{8\pi }\Big) - \nabla \Phi + \frac{1}{4\pi \rho }(\vec{B} \cdot \nabla )\vec{B} \\ \frac{D\vec{v}}{Dt} & = -\frac{1}{\rho }\nabla \Big(P + \frac{B^{2}}{8\pi }\Big) - \nabla \Phi + \frac{1}{4\pi \rho }(\vec{B} \cdot \nabla )\vec{B} \end{align}\,\!

Going into a frame co-rotating with RΩ:

\begin{align} u_{\phi } & = v_{\phi } - R\Omega \\ u_{z} & = v_{z} \\ u_{R} & = v_{R} \end{align}\,\!

Recalling the definition of epicycle frequency, the R and φ components of the Navier-Stokes equation are:

\begin{align} \frac{D u_{R}}{Dt} - 2\Omega u_{\phi } & = -\frac{1}{\rho } \frac{\partial }{\partial R}(P + \frac{B^{2}}{8\pi }) + \frac{1}{4\pi \rho }(\vec{B}\cdot \nabla )B_{R} \\ \frac{Du_{\phi }}{Dt} + \frac{\kappa ^{2}}{2\Omega }u_{R} & = -\frac{1}{\rho }\frac{1}{R}(P + \frac{B^{2}}{8\pi }) + \frac{1}{4\pi \rho }(\vec{B}\cdot \nabla )B_{\phi } \end{align}\,\!

The induction equation (for $\sigma \rightarrow 0$ then $\eta \rightarrow 0$) can be expanded out using a vector identity (and noting no magnetic monopoles and incompressibility) to become:

\begin{align} \frac{D\vec{B}}{Dt} & = (\vec{B} \cdot \nabla )\vec{v} \end{align}\,\!

Again, the R and φ components are looked at:

\begin{align} \frac{DB_{R}}{Dt} & = (\vec{B} \cdot \nabla )v_{R} = (\vec{B} \cdot \nabla ) u_{R} \\ \frac{DB_{\phi }}{Dt} & = R(\vec{B} \cdot \nabla )(\Omega + \frac{u_{\phi }}{R}) \end{align}\,\!

If there are small perturbations $\vec{u} = \vec{u}_{1}$ and $\vec{B} = \vec{B}_{0}\hat{z} + \vec{B}_{1}$, then the perturbation modes of interest are those with $\vec{k} \parallel \vec{B}_{0}$ and $\vec{v}_{1} \perp \hat{k}$ in the disk plane. Going into Fourier space, with $\vec{u}_{1} \propto e^{i(kz-\omega t)}$ and $\vec{B}_{1} \propto e^{i(kz-\omega t)}$, the four equations for R and φ are linearized in k-space:

\begin{align} -i \omega u_{1R} - 2\Omega u_{1\phi } & = \frac{ik}{4\pi \rho }B_{0}B_{1R} \\ -i \omega u_{1\phi } + \frac{\kappa ^{2}}{2\Omega }u_{1R} & = \frac{ik}{4\pi \rho }B_{0}B_{1\phi } \\ -i \omega B_{1R} & = ikB_{0}u_{1R} \\ -i \omega B_{1\phi } & = ikB_{0}u_{1\phi } + B_{1R}\frac{\partial \Omega }{\partial ln R} \end{align}\,\!

There are four equations and four unknowns, and combining them yields the dispersion relation:

\begin{align} \omega ^{4} - \omega ^{2}(\kappa ^{2}+2(kv_{A})^{2})+(kv_{A})^{2}\Big((kv_{A})^{2}+\frac{\partial \Omega ^{2}}{\partial ln R}\Big) & = 0 \end{align}\,\!

Several comments can be made about this relation. First, if there is no rotation, then κ = 0 and Ω = 0, and Alfven waves are recovered ($\omega ^{2} = k^{2}v_{A}^{2}$). This makes sense because the waves considered had $\hat{k} \parallel \vec{B}_{0}$ and $\vec{v}_{1} \perp \hat{k}$. Second, the disk is MRI unstable if ω2 < 0, which occurs when $k^{2}v_{A}^{2} < -\frac{\partial \Omega ^{2}}{\partial ln R}$. For Keplerian orbits, $\Omega \propto \frac{v_{R}}{R} \propto \frac{R^{-1/2}}{R} \propto R^{-3/2}$, so the dispersion relation is $k^{2}v_{A}^{2} < 3\Omega ^{2}$. Third, unstable modes grow with e | ω | t, and | ω | has a maximum that occurs when $|\omega _{max}| = \frac{1}{2}|\frac{d\Omega }{d lnR}| = \frac{3}{4}\Omega$ for a Keplerian disk. This is fast! Less than one rotation gives a whole exponential growth factor.

Finally, what is the physical mechanism behind MRI? Rotation de-stabilizes an accretion disk system because of magnetic tension. The key here is that the angular velocity Ω decreases with larger R, so particles in an inner orbit move ahead of particles in an outer orbit, which increases the magnetic tension (‘massless spring’) and the inner orbit particle feels a retarding torque, losing angular momentum and falling to a smaller orbit (while the outer particle gains angular momentum and moves outward). This stretches the spring more, increasing the torques and the motion becomes unstable.

Without rotation, a particle displaced perpendicularly to a magnetic field line feels a restoring force that stabilizes. However, this is not true for a rotating system and this is what leads to MRIs.

### VII. Turbulence

#### Laminar vs. Turbulent Flows

Turbulence is common in high Reynold number flows and common in daily lives. Properties of turbulence include disorder (irregularity in time and space), a wide range of interacting scales, vorticity-eddies, large dissipation, and efficient transport and mixing of momentum, heat, chemicals, etc.

A flow in which the kinetic energy dies out due to viscosity is a laminar flow.

In 1941, Kolmogorov postulated that for very high Reynolds numbers, small scale turbulent motions are statistically isotropic. He also hypothesized that the statistics are universally and uniquely determined by the scale size and energy dissipation/cascade rate $\varepsilon$. Turbulent flows have a hierarchy of scales through which energy cascade takes place. For example, dissipation of kinetic energy takes place at scales on the order of the Kolmogorov length, while the input of energy comes from the decay of larger scales. So, there is a range of scales that have formed at the expense of the energy of the large ones. Large eddies break into smaller eddies, which give rise to even smaller eddies, etc. until a sufficiently small length scale is reached that viscosity can dissipate the kinetic energy.

The mean kinetic energy (per mass) of a flow is related to the energy spectrum function E(k) by:

\begin{align} \frac{1}{2}\langle u_{i}u_{i}\rangle & = \int E(k) dk \end{align}\,\!

Therefore, $E(k) \propto \Big[\frac{l^{3}}{t^{2}}\Big]$. The energy cascade rate is $\varepsilon \propto \Big[\frac{l^{2}}{t^{3}}\Big]$. If the energy spectrum function is written like $E(k) \propto \varepsilon ^{\alpha } k^{\beta }$, where $k \propto 1/l$, then dimensional analysis produces the characteristic Kolmogorov spectrum:

\begin{align} E(k) \propto \varepsilon ^{2/3}k^{-5/3} \end{align}\,\!

On very small scales (large k-scales), then dissipation plays a role, and on very large scales (small k-scales), the system is no longer isotropic (geometrical features of the boundaries come into play). However, between these scales, the Kolmogorov spectrum describes the distribution of energy as a function of eddy/vortex size.

### Appendix: Mathematical Formulas

\begin{align} d\vec{r} & = dx\hat{x} + dy\hat{y} +dz\hat{z} \\ & = dr\hat{r} + rd\theta \hat{\theta } + r sin\theta d\phi \hat{\phi } \\ & = dr\hat{r} + rd\phi \hat{\phi } + dz\hat{z} \\ (\vec{v} \cdot \nabla ) Q & = v_{x}\frac{\partial Q}{\partial x} + v_{y}\frac{\partial Q}{\partial y} + v_{z}\frac{\partial Q}{\partial z} \\ & = v_{r} \frac{\partial Q}{\partial r} + \frac{v_{\theta }}{r}\frac{\partial Q}{\partial \theta } + \frac{v_{\phi }}{r sin \theta }\frac{\partial Q}{\partial \phi } \\ & = v_{r}\frac{\partial Q}{\partial r} + \frac{v_{\phi }}{r}\frac{\partial Q}{\partial \phi } + v_{z}\frac{\partial Q}{\partial z} \\ \nabla ^{2} & = \frac{\partial ^{2}}{\partial x^{2}} + \frac{\partial ^{2}}{\partial y^{2}} + \frac{\partial ^{2}}{\partial z^{2}} \\ & = \frac{1}{r^{2}}\frac{\partial }{\partial r}\Big(r^{2} \frac{\partial }{\partial r}\Big) + \frac{1}{r^{2}sin\theta }\frac{\partial }{\partial \theta }\Big(sin\theta \frac{\partial }{\partial \theta }\Big) + \frac{1}{r^{2} sin^{2} \theta }\frac{\partial ^{2}}{\partial \phi ^{2}} \\ & = \frac{1}{r}\frac{\partial }{\partial r}\Big(r\frac{\partial }{\partial r}\Big) + \frac{1}{r^{2}}\frac{\partial ^{2}}{\partial \phi ^{2}} + \frac{\partial ^{2}}{\partial z^{2}} \end{align}\,\!