# Fluid Dynamics Course Notes

### 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:

{\displaystyle {\begin{aligned}\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{aligned}}\,\!}

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: ${\displaystyle L\ll {\frac {Q}{|\nabla Q|}}}$ \2) Large enough to contain many particles: ${\displaystyle nL^{3}\gg 1}$ \3) Large enough that particles ‘know’ about local conditions through collisions: ${\displaystyle 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: ${\displaystyle \rho =\rho ({\vec {r}},t)}$) \2) Lagrangian: Examining the change in variables of a particular fluid element ${\displaystyle {\vec {a}}}$; the reference system is comoving with the fluid (ex: ${\displaystyle \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’.

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

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

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

#### Conservation Laws; Continuity Equation\

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

{\displaystyle {\begin{aligned}-{\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{aligned}}\,\!}

For an incompressible flow, ${\displaystyle D\rho /Dt=0}$, and fluid elements preserve their density along their paths. Therefore, ${\displaystyle \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 ${\displaystyle i}$ acting on a surface whose normal is ${\displaystyle j}$) by ${\displaystyle F_{i}=T_{ij}S_{j}}$. The net force acting over the surface due to pressure is:

{\displaystyle {\begin{aligned}{\vec {F}}&=-\int Pd{\vec {S}}\\&=-\int \nabla PdV\end{aligned}}\,\!}

Therefore, the acceleration (force per mass) contribution from pressure is ${\displaystyle -{\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:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

This is the conservative form of the momentum equation, where ${\displaystyle T_{ij}=\rho v_{i}v_{j}+P\delta _{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 ${\displaystyle i}$ and ${\displaystyle 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 ${\displaystyle {\vec {g}}=-\nabla \Phi =-{\frac {GM}{r^{2}}}{\hat {r}}}$.

{\displaystyle {\begin{aligned}\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{aligned}}\,\!}

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, ${\displaystyle \Phi ,\rho ,{\vec {v}},}$ and ${\displaystyle P}$, and three equations (Continuity, Euler, Poisson). We need an equation of state relating ${\displaystyle 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 ${\displaystyle dQ=TdS=d\varepsilon +PdV}$, where ${\displaystyle dQ}$ is the heat absorbed by the fluid, ${\displaystyle PdV}$ is the work done by the fluid, and ${\displaystyle d\varepsilon }$ is the internal energy content of the fluid. The specific heat equations are:

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

Therefore, ${\displaystyle dQ=C_{V}dT+PdV}$. To find how the specific heats are related to each other, the equation of state for an ideal gas is substituted in.

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

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

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

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

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

The total energy density (energy per volume) of a fluid is ${\displaystyle E=\rho ({\frac {1}{2}}v^{2}+\Phi +\epsilon )}$, where ${\displaystyle \epsilon }$ 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.

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

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 ${\displaystyle H=\varepsilon +PV}$, ${\displaystyle h}$ can be defined as enthalpy per unit mass: ${\displaystyle h=\epsilon +{\frac {P}{\rho }}}$. Therefore, the energy flux is ${\displaystyle (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 ${\displaystyle m}$, momentum ${\displaystyle m{\vec {v}}}$, and energy ${\displaystyle \varepsilon }$. Each of their conservation laws involve the rate of change of their densities (${\displaystyle \rho }$, ${\displaystyle \rho {\vec {v}}}$, and ${\displaystyle E}$, respectively). Their fluxes are given by their densities times a velocity (${\displaystyle \rho {\vec {v}}}$, ${\displaystyle T_{ij}}$, and ${\displaystyle (E+P){\vec {v}}}$, respectively). \

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

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

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

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

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

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

Because of spherical symmetry, a new variable can be defined: ${\displaystyle \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.

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

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

There are analytic solutions to the Lane-Emden equation for ${\displaystyle 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 ${\displaystyle \theta =0}$) increases with increasing ${\displaystyle n}$. \

A useful mass-radius scaling relation comes out of the Lane-Emden equation, since all stars with the same ${\displaystyle n}$ have the same ${\displaystyle \theta (\xi )}$.The central density ${\displaystyle \rho _{c}}$ sets the physical scale of ${\displaystyle \rho }$ and ${\displaystyle r}$.

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

From the definition of ${\displaystyle \xi }$, ${\displaystyle r\propto \rho _{c}^{{\frac {1}{2}}({\frac {1}{n}}-1)}}$, so the mass-radius relation for polytropic stars is found to be ${\displaystyle M\propto R^{\frac {3-n}{1-n}}}$. For an adiabatic star, ${\displaystyle 1+{\frac {1}{n}}=\gamma }$ and ${\displaystyle \gamma =5/3}$ for a monatomic gas, so ${\displaystyle 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 ${\displaystyle {\vec {w}}=\nabla \times {\vec {v}}}$.

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

There are two special cases that can be applied here. The first is the case of steady flows of ideal fluids. Recalling enthalpy, ${\displaystyle h=\epsilon +{\frac {P}{\rho }}}$, so ${\displaystyle dh=d\epsilon +Pd({\frac {1}{\rho }})+{\frac {1}{\rho }}dP=Tds+{\frac {1}{\rho }}dP}$, where everything is per mass. Therefore, for constant ${\displaystyle S}$ (ideal fluids), ${\displaystyle \nabla h={\frac {1}{\rho }}\nabla P}$. Re-writting the Euler equation, dotting it with ${\displaystyle {\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.

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

For steady flows, ${\displaystyle {\frac {\partial {\mathcal {B}}}{\partial t}}=0}$, and since ${\displaystyle {\vec {v}}\cdot \nabla {\mathcal {B}}=0}$, then ${\displaystyle {\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 ${\displaystyle {\vec {w}}=0}$. In this case ${\displaystyle \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:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

Taking the cross product of ${\displaystyle \nabla }$ with the above yields ${\displaystyle {\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 ${\displaystyle {\vec {v}}_{0}=0}$, the perturbations are:

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

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 ${\displaystyle {\frac {DQ}{Dt}}={\frac {\partial Q}{\partial t}}}$. The first order fluid equations are:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

Taking ${\displaystyle {\frac {\partial }{\partial t}}}$ of the continuity equation and substituting the momentum equation yields ${\displaystyle {\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 ${\displaystyle v_{s}^{2}={\frac {\partial P}{\partial \rho }}={\frac {P_{1}}{\rho _{1}}}}$. Therefore, the equation simplifies to the free wave equation:

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

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

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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 ${\displaystyle {\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, ${\displaystyle v_{s}^{2}={\frac {\partial P}{\partial \rho }}|_{T}={\frac {k_{B}T}{\mu m_{p}}}}$ and adiabatic sound speed, ${\displaystyle v_{s}^{2}={\frac {\partial P}{\partial \rho }}|_{S}=\gamma {\frac {k_{B}T}{\mu m_{p}}}}$. Therefore, ${\displaystyle v_{s,adiabatic}={\sqrt {\gamma }}v_{s,isothermal}}$. For air (diatomic), ${\displaystyle \gamma }$ 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:

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

The 0th order momentum equation is:

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

The 1st order fluid equations are::

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

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

Putting the equations together yields Laplace’s equation:

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

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

{\displaystyle {\begin{aligned}-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{aligned}}\,\!}

There’s a boundary condition that ${\displaystyle v_{1z}=0}$ at the bottom of the sea (${\displaystyle z=-H}$). This implies that ${\displaystyle B=e^{-2kH}}$. Recalling ${\displaystyle 2cosh(x)=e^{x}+e^{-x}}$ and ${\displaystyle 2sinh(x)=e^{x}-e^{-x}}$, plugging this into the previous solutions gives:

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

Before deriving the dispersion relation, the concept of a fluid displacement ${\displaystyle {\vec {\xi }}({\vec {r}},t)}$ must be introduced. The perturbed fluid velocity is related to this displacement by ${\displaystyle {\vec {v}}_{1}={\frac {d{\vec {\xi }}}{dt}}}$, so ${\displaystyle \xi _{z}={\frac {v_{1z}}{-i\omega }}}$. Also, at the air-water interface (${\displaystyle z=0}$), a boundary condition is that there is constant pressure along a fluid element. Mathematically, ${\displaystyle P_{1}+{\vec {\xi }}\cdot \nabla P_{0}=0}$ at ${\displaystyle 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: ${\displaystyle 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: ${\displaystyle P({\vec {r}}+{\vec {\xi }},t)=P_{0}({\vec {r}})+P_{1L}}$. Solving for ${\displaystyle P_{1L}}$ gives ${\displaystyle 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, ${\displaystyle P_{1L}=0}$, and the surface boundary condition follows:

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

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

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

There are two limits worth discussing. The first is the deep water/short wavelength limit: ${\displaystyle kH\gg 1}$ or ${\displaystyle H\gg 1/k\sim \lambda }$. In this case, ${\displaystyle tanh[kH]\sim 1}$ and ${\displaystyle \omega ={\sqrt {gk}}}$. The phase velocity in this regime is ${\displaystyle 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 ${\displaystyle H}$, which means that small disturbances don’t know about the sea floor. The second regime is the shallow water/long wavelength limit: ${\displaystyle kH\ll 1}$ or ${\displaystyle H\ll \lambda }$. In this case, ${\displaystyle tanh[kH]\sim kH}$ and ${\displaystyle \omega ={\sqrt {gH}}k}$. The phase and group velocities are equal to each other: ${\displaystyle 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 ${\displaystyle v_{s}}$, the Mach number is ${\displaystyle M={\frac {v}{v_{s}}}={\frac {1}{sin\alpha }}}$, where ${\displaystyle \alpha }$ 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 ${\displaystyle dx}$ of the continuity equation), noting that the mass flux in equals the mass flux out:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

Applying the same method to the conservative form of the momentum equation (and noting that ${\displaystyle \Phi }$ is continuous across shocks) yields:

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

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

The expression for ${\displaystyle \epsilon }$ can be substituted into the energy conservation equation, and ${\displaystyle 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:

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

In the strong shock limit, ${\displaystyle M_{1}\gg 1}$, ${\displaystyle {\frac {v_{1}}{v_{2}}}={\frac {\gamma +1}{\gamma -1}}}$, ${\displaystyle {\frac {P_{2}}{P_{1}}}\gg 1}$, ${\displaystyle {\frac {T_{2}}{T_{1}}}\gg 1}$, and ${\displaystyle 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, ${\displaystyle T_{1}=T_{2}}$ 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 ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle \rho _{1}}$ before the explosion and a post-shock density ${\displaystyle \rho _{2}}$ inside the blast wave shell. In the pre-shock’s rest frame, the post-shock velocity is $\displaystyle 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 ${\displaystyle R}$ gets swept up into a shell of thickness ${\displaystyle D}$. For a monatomic gas with ${\displaystyle \gamma =5/3}$:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

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

This scaling holds if ${\displaystyle 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: ${\displaystyle \xi ={\frac {r}{R(t)}}}$. Therefore, any parameter ${\displaystyle X}$ can be plotted as a function of ${\displaystyle \xi }$ instead of ${\displaystyle r}$ and it will have the same shape but be scaled up or down by a time-dependent factor: ${\displaystyle X=X_{2}(t){\tilde {X}}(\xi )}$. In order to find these self-similar solutions, variables are changed from ${\displaystyle (r,t)}$ to ${\displaystyle (\xi ,R)}$.

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

For example, ${\displaystyle r=\xi {\Big (}{\frac {Et^{2}}{\rho _{1}}}{\Big )}^{1/5}}$, so the velocity of the shock is ${\displaystyle 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 ${\displaystyle t}$ can be solved for as a function of ${\displaystyle \xi }$, and therefore ${\displaystyle r}$ is known as a function of ${\displaystyle \xi }$, as well as ${\displaystyle M}$ and other parameters depending on ${\displaystyle 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:

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

The 0th order fluid equations are:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

The first equation is trivial, since ${\displaystyle {\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 ${\displaystyle P_{0}}$ is constant in time and space. Therefore, that leaves ${\displaystyle \nabla \Phi _{0}=0}$ which contradicts Poisson’s equation that states ${\displaystyle \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:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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:

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

The solution to this equation is ${\displaystyle \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:

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

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

Also note that Jean’s wavelength is ${\displaystyle \lambda _{J}={\frac {2\pi }{k_{J}}}}$ and Jean’s mass is ${\displaystyle M_{J}={\frac {4}{3}}\pi \rho _{0}\lambda _{J}^{3}}$. For ${\displaystyle \lambda >\lambda _{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 ${\displaystyle F_{gravity}>F_{pressure}}$:

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

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):

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

The 0th order fluid equations are:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

The continuity equation can be simplified by substituting ${\displaystyle \rho _{1}=\delta \rho _{0}}$ and the relationship found before that ${\displaystyle {\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:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

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

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

Again, instability occurs when ${\displaystyle k. \

Taking the case where the ${\displaystyle k}$ term is ignored, so ${\displaystyle {\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, ${\displaystyle \Omega =\Omega _{m}=1}$, ${\displaystyle 4\pi G\rho ={\frac {3}{2}}H^{2}}$, and ${\displaystyle a\propto t^{2/3}}$. Therefore, ${\displaystyle {\ddot {\delta }}_{k}+{\frac {4}{3t}}{\dot {\delta }}_{k}-{\frac {2}{3t^{2}}}\delta _{k}=0}$. Taking an ansatz of ${\displaystyle \delta _{k}\propto t^{\alpha }}$ results in two solutions (growing and decaying):

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

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

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

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

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

Integrating the first order Poisson equation from ${\displaystyle z=-\epsilon }$ to ${\displaystyle z=\epsilon }$ gives:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

Using Fourier Space, the first order continuity equation yields ${\displaystyle 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 ${\displaystyle v_{1y}=-{\frac {2\Omega i}{k}}{\frac {\Sigma _{1}}{\Sigma _{0}}}}$. Substituting relationships into the x-component yields the dispersion relation:

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

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 (${\displaystyle \omega ^{2}<0}$). If there is no rotation (only gravity and pressure), ${\displaystyle \omega ^{2}<0}$ for a Jeans wavenumber ${\displaystyle k (long modes are unstable). If there is negligible sound speed (no pressure), then ${\displaystyle \omega ^{2}<0}$ for ${\displaystyle 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:

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

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

For a differentially rotating disk of gas, ${\displaystyle \Omega =\Omega (R)}$ is not a constant so it needs to be replaced by the epicycle frequency.

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

The form of ${\displaystyle \kappa }$ depends on ${\displaystyle v_{circ}(R)=R\Omega (R)}$ or the mass distribution. For Keplerian orbits, ${\displaystyle \kappa =\Omega }$. For flat rotation curves, ${\displaystyle \kappa ={\sqrt {2}}\Omega }$. And for solid-body rotation, ${\displaystyle \kappa =2\Omega }$. In the case of solid-body rotation, the disk is stable if ${\displaystyle 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:

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

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

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

There are two boundary conditions to apply here. The first is that fluid displacements are continuous at ${\displaystyle z=0}$. Remember that ${\displaystyle {\vec {v}}_{1}={\frac {d{\vec {\xi }}}{dt}}}$, then ${\displaystyle \xi _{z}={\frac {v_{1z}}{-i\omega }}}$ for both layers of fluids at ${\displaystyle 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 ${\displaystyle P_{+,1L}=P_{-,1L}}$. With surface tension, this means that ${\displaystyle P_{+,1L}-{\mathcal {T}}{\frac {\partial ^{2}\xi _{z}}{\partial x^{2}}}=P_{-,1L}}$. These both apply at ${\displaystyle z=0}$. The left-hand side of the second equation is the total downward force per area. Recall also that ${\displaystyle P_{1L}=P_{1}+{\vec {\xi }}\nabla P_{0}}$, where ${\displaystyle P_{1L}}$ was set to zero for the water waves boundary condition (constant pressure along a fluid element). Applying this to the top fluid:

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

Doing the same thing for the bottom fluid:

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

The two equations are related by the surface tension term.

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

Simplifying this expression and noting that ${\displaystyle v_{+,1z}=v_{-,1z}}$ because of the first boundary condition produces the dispersion relation:

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

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

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

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

And finally, note that the downward gravitational acceleration ${\displaystyle {\vec {g}}}$ can be mapped into an upward ${\displaystyle {\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 ${\displaystyle {\vec {v}}_{+,0}\neq 0}$ and ${\displaystyle {\vec {v}}_{-,0}\neq 0}$. The resulting dispersion relation is:

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

This can be rewritten as a quadratic equation for ${\displaystyle \omega /k}$. When surface tension is negligible, instability arises for short wavelength modes, or:

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

Note that if surface tension is present, the short modes can be stabilized. Also, if ${\displaystyle g=0}$, any wavenumber is Kelvin-Helmholtz unstable. If ${\displaystyle k_{crit}}$ is small, Kelvin-Helmholtz instabilities are common, and if ${\displaystyle k_{crit}}$ 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, ${\displaystyle {\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.

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

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 ${\displaystyle x}$ (from particles at ${\displaystyle x+l/2}$ and ${\displaystyle x-l/2}$) is ${\displaystyle p_{y}(x)\sim m{\frac {\partial v_{y}}{\partial x}}l}$. The number of particles per unit time passing the boundary at ${\displaystyle x}$ is ${\displaystyle N\sim nAv_{th}}$, where ${\displaystyle v_{th}\sim {\sqrt {\frac {kT}{m}}}}$. Therefore, ${\displaystyle {\dot {p}}_{y}(x)\sim (m{\frac {\partial v_{y}}{\partial x}}l)(nAv_{th})}$. The viscous force on a volume element ${\displaystyle A\Delta x}$ is ${\displaystyle {\frac {\partial {\dot {p}}_{y}}{\partial x}}\Delta x}$, and therefore the viscous force per unit volume is ${\displaystyle {\frac {F_{y}}{A\Delta x}}\sim {\frac {\partial }{\partial x}}(\rho \nu {\frac {\partial v_{y}}{\partial x}})}$, where ${\displaystyle \rho =mn}$ and ${\displaystyle \nu =lv_{th}}$. ${\displaystyle \nu }$ is called the kinematic viscosity. The dynamic viscosity is ${\displaystyle \eta =\rho \nu }$, and the viscous shear stress tensor is ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle \sigma _{ij}=0}$ for uniform solid-body rotation (${\displaystyle \alpha =1}$) and ${\displaystyle \sigma _{ij}=0}$ for uniform isotropic expansion (such as the Hubble Flow) (${\displaystyle \beta =-2/3}$). Putting things together:

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

#### Navier-Stokes Equation; Reynolds Number\

Taking the derivative of the viscous tensor yields ${\displaystyle \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:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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 ${\displaystyle [cm^{2}/s]}$ (${\displaystyle 0.01}$ for water, ${\displaystyle 0.15}$ for air, ${\displaystyle 1}$ for oil, etc.). \

The viscous stress tensor is symmetric and traceless. If there is bulk viscosity (ex: shocks), then diagonal terms are added: ${\displaystyle \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 (${\displaystyle {\vec {v}}\cdot \nabla {\vec {v}}}$) to viscous forces (${\displaystyle \nu \nabla ^{2}{\vec {v}}}$).

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

For ${\displaystyle {\mathcal {R}}\ll 1}$, viscosity is important, and for ${\displaystyle {\mathcal {R}}\gg 1}$, the flow is inviscid. For very high Reynolds number (${\displaystyle {\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:

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

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

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

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

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

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:

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

Another example of a viscous flow is Stoke’s flow. This is a ‘creeping flow’ around a sphere of radius ${\displaystyle a}$ moving through a viscous fluid. For a constant, low-Reynolds number viscous flow of an incompressible fluid, the Navier-Stokes equation becomes ${\displaystyle 0=-\nabla P-\eta \nabla ^{2}{\vec {v}}}$. Taking the divergence and curl of this expression gives ${\displaystyle \nabla ^{2}P=0}$ and ${\displaystyle \nabla ^{2}{\vec {w}}=0}$, respectively. The boundary condition of this problem is ${\displaystyle {\vec {v}}=0}$ at ${\displaystyle r=a}$ and the general solution for velocity is ${\displaystyle {\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 ${\displaystyle B}$ in terms of ${\displaystyle A}$, solving for the ${\displaystyle \phi }$-component of vorticity, taking ${\displaystyle \nabla ^{2}w=0}$ to find the form of ${\displaystyle A}$ and using the boundary condition to find ${\displaystyle B}$ and the ${\displaystyle r}$ and ${\displaystyle \theta }$ 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 (${\displaystyle \int F_{z}a^{2}sin\theta d\theta d\phi }$) gives the drag force on the sphere.

{\displaystyle {\begin{aligned}F_{drag}&=6\pi \eta va\end{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}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{aligned}}\,\!}

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

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

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

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

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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 ${\displaystyle T=2\pi \nu \Sigma R^{3}{\frac {\partial \Omega }{\partial R}}}$, then the right hand side of the equation is ${\displaystyle {\frac {1}{2\pi R}}{\frac {\partial T}{\partial R}}}$. Recall also that the torque is the viscous force times ${\displaystyle R}$, where the viscous force is the product of the area and the viscous stress. \

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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 ${\displaystyle \nu }$. Changing variables from ${\displaystyle 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 ${\displaystyle x={\frac {R}{R_{0}}}}$ and ${\displaystyle \tau ={\frac {12\nu }{R_{0}^{2}}}t}$, then the surface density can go from being a function of ${\displaystyle R}$ and ${\displaystyle t}$ to ${\displaystyle x}$ and ${\displaystyle \tau }$.

{\displaystyle {\begin{aligned}\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{aligned}}\,\!}

The main take-away here is that viscosity spreads out a ring. Greater values of ${\displaystyle \tau }$ 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:

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

And using the dimensionless variables, this becomes:

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

The behavior of ${\displaystyle I_{1/4}(z)}$ is flat for low z’s and steeper for high z’s. Therefore, for ${\displaystyle 2x>\tau }$:

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

And for ${\displaystyle 2x<\tau }$:

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

For ${\displaystyle x>1}$, ${\displaystyle R>R_{0}}$, and ${\displaystyle v_{R}>0}$. Therefore, the outer parts move outward. For ${\displaystyle \tau >4x^{2}}$, ${\displaystyle v_{R}<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 ${\displaystyle v_{R}\sim {\frac {\nu }{R}}}$, so the timescale for orbital evolution is ${\displaystyle 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 ${\displaystyle \nu \sim v_{thermal}\lambda _{mfp}\sim v_{thermal}{\frac {1}{\sigma n}}}$, and typical numbers give ${\displaystyle {\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 ${\displaystyle \nu \sim v_{turb}l_{turb}\sim v_{s}H}$, where ${\displaystyle H}$ is the disk thickness. Disks can therefore be characterized as ${\displaystyle \nu =\alpha v_{s}H}$, where ${\displaystyle \alpha }$ 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: ${\displaystyle f({\vec {x}},{\vec {p}},t)}$. The number of particles is ${\displaystyle dN=f({\vec {x}},{\vec {p}},t)d^{3}xd^{3}p}$. A simple example of a phase space distribution is:

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

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

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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 (${\displaystyle {\dot {x}}_{i}={\frac {\partial H}{\partial p_{i}}}}$ and ${\displaystyle {\dot {p}}_{i}=-{\frac {\partial H}{\partial x_{i}}}}$), the equation becomes:

{\displaystyle {\begin{aligned}{\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{aligned}}\,\!}

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

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