Past exam of the mathematics course of the University of Cambridge / 2007 / iii / Paper 65 (Astrophysical Fluid Dynamics)

In the energy-conserving blast regime, neglect ambient pressure, radiative losses, ejecta mass and gravity compared with the explosion. The cold ambient ideal gas has no sound speed scale, and its power-law mass density has no preferred radius. The point explosion introduces only its energy $E$ and elapsed time $t$. This scale-free setting supports a self-similar blast wave; extra physics would limit the time interval over which that description applies.

The mass swept up to radius $R$ isIt is finite because $\beta <3$. With no energy loss, the kinetic energy and internal energy both scale as $C{R}^{3-\beta }{\dot{R}}^2$. Therefore $R^{3-\beta }{\dot{R}}^2$ is constant and integration givesThe dimensionless coefficient $A$ depends on the adiabatic exponent and ambient exponent through the similarity profiles. Indeed, $[E/C]=L^{5-\beta }{T}^{-2}$, which independently verifies the scaling. The shock wave decelerates for $0<\beta <3$, since its time exponent lies between $2/5$ and $1$.

Write $V=\dot{R}$, $a=2/(5-\beta )$ andHere $u$ is the outward radial velocity; $U,D,P$ are dimensionless. The relationsare sufficient for the reduction. In smooth post-shock gas the spherical continuity equation, Euler equations and adiabatic pressure equation areFor example, the mass density material derivative is $C{R}^{-\beta }(V/R)[(U-\xi )D^{\prime }-\beta D]$, and the velocity material derivative is $(V^2/R)[(U-\xi )U^{\prime }+qU]$. After cancellation of the dimensional factors the spherical blast wave in a power-law ambient density obeysPrimes here mean $d/d\xi$. In the pressure equation the scale derivative $2q-\beta =-3$ explains the constant term. Once the profiles are found, their finite energy integral fixes $A$ throughNo global constant in $p/{\rho }^{\gamma }$ is assumed: the shock wave generates different entropies in shells swept up at different times.

Assume $\gamma >1$. The upstream ideal gas is at rest in the laboratory frame and has mass density ${\rho }_1=C{R}^{-\beta }$. In the shock frame its incident speed is $V$. The strong-shock limit of the Rankine-Hugoniot conditions for a perfect gas gives compression $\chi =(\gamma +1)/(\gamma -1)$ and downstream shock-frame speed $V/\chi$. Transforming back to the laboratory frame gives $u_2=V(1-1/\chi )=2V/(\gamma +1)$.

The normal momentum jump yieldswhere the cold upstream pressure has been neglected. Thus the boundary values in the normalization of part (b) areThe velocity jump uses the relative speed $V-u_2$, not the laboratory speed $u_2$. This distinction is essential to all three shock wave values.

Put $U=b\xi$ and $D=d\xi$. The shock wave conditions fix $b=2/(\gamma +1)$ and $d=(\gamma +1)/(\gamma -1)$. The continuity equation then reduces toso a linear mass density profile requiresMomentum requires $P^{\prime }=-db(b-1+q){\xi }^2$. At this exponent $q=-2(\gamma -1)/(\gamma +1)$, hence integration gives $P=2{\xi }^3/(\gamma +1)+P_0$. Substitution in the pressure equation makes the cubic coefficient vanish because $(\gamma +1)b=2$, while its constant term requires $(3\gamma b-3)P_0=0$. For $\gamma >1$ this implies $P_0=0$. Therefore all three equations and all shock wave conditions hold for the linear-profile spherical blast wave,The prescribed range $0<\beta <3$ restricts this family to $1<\gamma <7$.

For $\gamma =5/3$ one has $\beta =2$, $a=2/3$, $U=3\xi /4$, $D=4\xi$ and $P=3{\xi }^3/4$. The energy integral isConsequently $E=(3\pi /2)CR{\dot{R}}^2=(2\pi /3)C{R}^3/t^2$, which fixes the radius and dimensional fields completely:The kinetic energy and internal energy integrands are each $9{\xi }^5/8$, so each accounts for half the explosion energy. Integrating the mass density also gives the correct swept-up mass $4\pi CR$.

Let the equilibrium ideal gas be at rest, with uniform ${\rho }_0,p_0$ and ${\mathbf{B}}_0=B_0{\mathbf{e}}_z$. This local wave calculation neglects gravity, viscosity and magnetic diffusion. Define the adiabatic sound speed $c\equiv v_s=(\gamma p_0/{\rho }_0{)}^{1/2}$, the Alfvén speed $v_a=B_0/({\mu }_0{\rho }_0{)}^{1/2}$, and a wavevector $\mathbf{k}=k(\sin \theta ,0,\cos \theta )$, where $k=|\mathbf{k}|$ and $\theta$ is its angle to the equilibrium magnetic field. Perturbations have dependence $e^{i\mathbf{k}\cdot \mathbf{x}-i\omega t}$ and phase speed $v=\omega /k$.

The linearized ideal magnetohydrodynamic equations for mass conservation, pressure, induction and momentum giveThe relation between pressure and mass density follows from eliminating $\mathbf{k}\cdot {\mathbf{u}}^{\prime }$ between the two material evolution equations; it applies to the nonzero-frequency wave sector. The induction result automatically has $\mathbf{k}\cdot \mathbf{b}=0$.

The $y$ MHD wave polarization is perpendicular to both the field and wavevector. It has $\delta p=\delta \rho =0$ and obeysThus the transverse Alfvén wave hasIn the $x$-$z$ plane the two remaining velocity components instead obey the symmetric matrix eigenproblemTaking its determinant gives $v^4-(c^2+v_a^2)v^2+c^2{v}_a^2{\cos }^2\theta =0$. Solving this quadratic equation in $v^2$ derives both magnetoacoustic wave phase speeds,The plus sign is the fast magnetosonic wave, the minus sign the slow magnetosonic wave. The radicand is nonnegative, since it equals $[(c^2-v_a^2{)}^2+4c^2{v}_a^2{\sin }^2\theta ]/4$. Both squared speeds are nonnegative. Each nonzero mode has forward and backward propagation; an advected entropy disturbance is a separate zero-frequency mode in this static equilibrium.

Set $ϵ=v_a^2/c^2\ll 1$. Expanding the square root in part (a) givesHence, to first absolute order in $ϵ$,The high-pressure limit of magnetosonic waves can also be expressed as nonnegative speeds. Retaining the next squared-speed term gives the first relative correction for the slow mode:The exact Alfvén speed along the propagation direction remains $v_a|\cos \theta |$.

The fast magnetosonic wave is predominantly an acoustic wave, with velocity almost parallel to $\mathbf{k}$ and appreciable mass density and gas pressure variations. Its magnetic pressure correction is weak. The Alfvén wave is exactly incompressible, polarized perpendicular to the $\mathbf{k}$-${\mathbf{B}}_0$ plane and restored by magnetic tension.

The slow magnetosonic wave is polarized inside that plane and almost perpendicular to $\mathbf{k}$. Indeed, the leading sound speed part of the matrix in (a) annihilates a vector perpendicular to $\mathbf{k}$, so only the smaller magnetic restoring force acts on that MHD wave polarization. Its compression is smaller by order $ϵ$, with gas pressure changes almost cancelling magnetic pressure changes. It is often called the pseudo-Alfvén wave MHD wave polarization in this high pressure limit. Thus two distinct almost-incompressible magnetic MHD wave polarizations have the same leading speed, while the third is predominantly sound. For exactly perpendicular propagation the Alfvén wave and slow magnetosonic wave frequencies vanish; for parallel propagation the two transverse magnetic MHD wave polarizations are degenerate.

Take $B_0=B_z>0$, $v_a=B_0/\sqrt{{\mu }_0{\rho }_0}$ and $\zeta =k(z-v_{a}t)$. The transverse ideal magnetohydrodynamic induction equation and magnetohydrodynamic momentum equation are satisfied by $u_x=-B_x/\sqrt{{\mu }_0{\rho }_0}$, $u_y=u_z=0$, with constant mass density. However, the longitudinal magnetohydrodynamic momentum equation also requiresHere $B_x^2=a^2{B}_0^2{\cos }^2\zeta$ is not constant. Balancing its gradient would requireThe adiabatic pressure equation instead gives ${\partial }_{t}p=0$, since this velocity is divergence-free and has no component along the spatial variation. No choice of $p_{\ast }(t)$ removes the spatially varying time derivative. Thus a nontrivial finite-amplitude linearly polarized sinusoidal Alfvén wave is not an exact pure traveling wave in a compressible ideal fluid. The magnetic-pressure obstruction to a linearly polarized Alfvén wave is quadratic in $a$, so the wave is a valid linear approximation and drives a compressive response at higher order. In an incompressible model pressure is a constraint variable; it is not governed by the same compressible adiabatic evolution law.

Allowing longitudinal flow and varying mass density does not rescue a regular one-dimensional periodic traveling solution with precisely this transverse magnetic field. To see this, let $U=u_z-v_a$ and use the traveling coordinate. Integrated continuity, induction and transverse momentum givewhere $m,A,D$ are constants. Eliminating $u_x$ givesAt any zero of $B_x$, finite mass density and velocity require the right side to vanish. Wherever $B_x\ne 0$, the coefficient then vanishes and fixes $\rho ={\mu }_0{m}^2/B_0^2$, constant; continuity fixes $U$, and adiabatic entropy transport fixes $p$. The longitudinal magnetohydrodynamic momentum equation again contradicts the varying $B_x^2$. The alternative $m=0$ makes transverse momentum force $B_x$ constant. Trivial cases such as $a=0$ or $k=0$ are excluded from the nontrivial wave conclusion.

For the circular MHD wave polarization, $B_x^2+B_y^2=a^2{B}_0^2$ is constant. Choose any constant ${\rho }_0>0,p_0>0$, put $v_a=B_0/\sqrt{{\mu }_0{\rho }_0}$ and takeThis gives a circularly polarized nonlinear Alfvén wave of arbitrary amplitude. All fields depend only on $z,t$, and $u_z=0$, so $\mathbf{u}\cdot \nabla \mathbf{u}=0$ and $\nabla \cdot \mathbf{u}=0$. Constant mass density and pressure then satisfy continuity and the adiabatic pressure equation, and $\nabla \cdot \mathbf{B}=0$ because $B_z$ is constant.

The remaining ideal magnetohydrodynamic induction equation and transverse magnetohydrodynamic momentum equation areSince ${\partial }_t=-v_a{\partial }_z$ on the traveling profiles and ${\mathbf{u}}_{\perp }=-(v_a/B_0){\mathbf{B}}_{\perp }$, both identities hold exactly when $v_a^2=B_0^2/({\mu }_0{\rho }_0)$. The longitudinal magnetohydrodynamic momentum equation holds because $p+|{\mathbf{B}}_{\perp }{|}^2/(2{\mu }_0)$ is constant. The circularly polarized wave is therefore an exact nonlinear solution even though the fluid is compressible. Compressibility permits mass density changes; it does not require every possible motion to compress the gas.

Use a steady axisymmetric flow governed by ideal magnetohydrodynamics and take the Newtonian gravitational potential to be axisymmetric. Work locally on connected regular magnetic flux surfaces with nonzero poloidal flow. As usual for a wind without an imposed toroidal loop voltage, $E_{\varphi }=0$; regularity on an included rotation axis guarantees this condition. It is needed to exclude steady cross-field flow in an annular domain.

Since $\nabla \varphi ={\mathbf{e}}_{\varphi }/R$, the poloidal magnetic flux function givesThe poloidal magnetic field is divergence-free. The zero toroidal electric field, with $\mathbf{E}=-\mathbf{u}\times \mathbf{B}$, gives ${\mathbf{u}}_p\parallel {\mathbf{B}}_p$. Define $k$ by $\rho {\mathbf{u}}_p=k{\mathbf{B}}_p$. Steady mass conservation then givesThus $k=k(\psi )$ is the magnetohydrodynamic mass loading, the mass flux per unit poloidal magnetic flux. Subtracting $(k/\rho )\mathbf{B}$ from the full velocity leaves only an azimuthal component. The cross product isThe steady ideal magnetohydrodynamic induction equation requires $0=\nabla \times (\omega \nabla \psi )=\nabla \omega \times \nabla \psi$. Consequently the field-line angular velocity is $\omega =\omega (\psi )$ and

Neither gas pressure nor axisymmetric gravity exerts an azimuthal torque. Multiplying the azimuthal momentum equation by $R$ combines the cylindrical curvature terms into angular-momentum transport:Insert $\rho {\mathbf{u}}_p=k{\mathbf{B}}_p$ and use ${\mathbf{B}}_p\cdot \nabla k=0$. It follows thatso the magnetohydrodynamic angular-momentum invariant isThe two terms account respectively for matter and Maxwell stress tensor transport of angular momentum.

Finally, subtracting $\gamma$ times the logarithmic mass density equation from the logarithmic pressure equation gives $D\ln (p/{\rho }^{\gamma })/Dt=0$. The specific entropy of a fixed-composition ideal gas is a function of $p/{\rho }^{\gamma }$. Steadiness and axisymmetry reduce its advection to ${\mathbf{u}}_p\cdot \nabla s=0$. HenceThese flux-label functions are local on connected regular surfaces; disconnected components with the same numerical label need not share constants without an additional matching condition.

On one poloidal flux surface, the field-line angular velocity $\omega (\psi )$ is constant. In a frame rotating at that rate the relative velocity is $\mathbf{v}=\mathbf{u}-R\omega {\mathbf{e}}_{\varphi }=k\mathbf{B}/\rho$, parallel to the magnetic field. The Lorentz force is perpendicular to this motion, and the Coriolis force also does no work. Centrifugal acceleration comes from the effective potential $-R^2{\omega }^2/2$. Since entropy is constant along the surface, pressure work integrates to specific enthalpy $w$, with $dw=dp/\rho$ along the flow.

Therefore kinetic energy relative to the rotating poloidal flux surface, enthalpy and effective gravitational-plus-centrifugal potential have a constant sum:This is a corotating energy integral, not the laboratory energy alone: Maxwell stress tensor stresses can transfer energy and angular momentum in the laboratory frame. In fact the expression equals $|\mathbf{u}{|}^2/2+w+\Phi -\omega R{u}_{\varphi }$, namely the total magnetohydrodynamic Bernoulli invariant minus $\omega \ell$. Differential rotation between surfaces does not invalidate the argument, because the chosen angular velocity is fixed along the particular surface being followed.

All invariants in this part refer to the selected magnetic field line. Its squared poloidal Alfvén number isAt the Alfvén surface crossing it equals unity, so ${\rho }_a={\mu }_0{k}^2$ and $M_{Ap}^2=1/y$. Subtract the two azimuthal integrals in part (a) to findFinite azimuthal fields and velocity at $y=1$ require the Alfvén-surface regularity condition for an axisymmetric wind $\ell =\omega R_a^2$. HenceHere $v_{\varphi }$ denotes the azimuthal velocity relative to field-line rotation. For gravity of a point mass at $z=0$, $\Phi =-GM/R=-GM/(R_{a}x)$.

Insert these quantities in the corotating energy integral with $w=0$, and divide by the gravitational energy scale $GM/R_a$. This gives the cold radial magnetohydrodynamic wind integral,withBoth are dimensionless. The apparent $0/0$ at $x=y=1$ must be evaluated on the smooth solution branch, not assigned an arbitrary value.

A constant-energy curve obeys $f_x+f_{y}dy/dx=0$. At a nondegenerate magnetosonic point, the vanishing coefficient of the mass density derivative requires simultaneous vanishing of its numerator:The curve must admit a real crossing slope. Expanding about such a point gives $f_{yy}{m}^2+2f_{xy}m+f_{xx}=0$, $m=dy/dx$, so a regular saddle point of a scalar function crossing has $f_{xy}^2-f_{xx}{f}_{yy}>0$. These are the form of the slow and fast regularity conditions in a warm wind, with the enthalpy contribution retained in $f$.

In the strictly cold approximation used to obtain the displayed $f$, the slow speed is zero. There is a cold-limit degeneracy of a slow magnetosonic point: a finite-speed interior slow crossing need not survive. Indeed,so this cold equation has no ordinary sub-Alfvénic slow critical point. A small retained enthalpy is needed to locate that point, which can approach the launching boundary as the sound speed tends to zero. The fast critical point still obeys the simultaneous conditions above wherever the reduced expression is nonsingular.

Let ${\rho }_0$ denote the midplane mass density and choose the symmetry plane at $z=0$. Isothermal hydrostatic equilibrium and the Poisson equation for Newtonian gravity giveso $c_s^2(d^2/d{z}^2)\ln \rho =-4\pi G\rho$. For $\rho ={\rho }_0{\mathrm{sech}}^2(z/H)$, direct differentiation yieldsThus the profile is a self-gravitating isothermal slab when $H^2=c_s^2/(2\pi G{\rho }_0)$. Integrating over both sides gives $\Sigma =2{\rho }_{0}H$, henceThe corresponding potential is ${\Phi }_0(z)-{\Phi }_0(0)=2c_s^2\ln \cosh (z/H)$, whose gradient vanishes at the midplane and tends to $\pm 2\pi G\Sigma$ far from the slab. The profile parameter $H$ is not its asymptotic exponential mass density scale height: that is $H/2$.

To distinguish Eulerian and Lagrangian fluid perturbations, use $\mathbf{\xi }=({\xi }_x,{\xi }_y,{\xi }_z)$ for the fluid displacement, $\varphi ={\Phi }^{\prime }$ for the Eulerian potential perturbation, and $\delta \rho$ for the Eulerian mass density perturbation. The equilibrium mass density $\rho (z)$ is the profile from part (a); subscript $z$ below denotes differentiation. Linearized mass conservation and the isothermal equation of state giveThe velocity perturbation is $-i\omega \mathbf{\xi }$. Linearized momentum initially has the formSince $\nabla {\Phi }_0=-c_s^2\nabla \ln \rho$, the pressure and perturbed mass density gravity terms combine into $-\rho \nabla (c_s^2\delta \rho /\rho )$. Define $W=c_s^2\delta \rho /\rho +\varphi$. The closed linearized system isThe $y$ component is a neutral transverse sector; it vanishes for nonzero-frequency modes.

To prove reality rather than assume it, multiply momentum by $\rho {\mathbf{\xi }}^{\ast }$ and integrate over $z$. With $I=\int \rho |\mathbf{\xi }{|}^{2}dz>0$, integration by parts and the conjugate continuity equation giveUsing the conjugate Poisson equation for Newtonian gravity in the last integral givesTake disturbances with finite $I$, finite pressure and gravitational integrals, and vanishing pressure work and gravitational boundary terms. For $k\ne 0$, require the potential to decay at both infinities and the fluid displacement to make $\rho {\xi }_z^{\ast }W\to 0$ as well. In a source-free outer region the decaying potential is proportional to $e^{-|k||z|}$. The isothermal displacement energy of a self-gravitating slab then givesSince $I$ is strictly positive for a nonzero displacement, ${\omega }^2$ is real. Positive values describe oscillations; negative values give pure growth and decay. The result is the energy identity for the pressure-gravity displacement operator, a self-adjoint operator in the mass density-weighted inner product. The identity excludes an oscillatory growing eigenmode under these boundary conditions. At $k=0$ the same identity applies when its boundary terms vanish and its integrals remain finite, with the gravitational integral involving only $|{\varphi }_z{|}^2$.

Consider the nontrivial gravitational mode with $k\ne 0$. At ${\omega }^2=0$, horizontal momentum requires $ikW=0$, hence $W=0$ andSubstituting into the Poisson equation for Newtonian gravity and using $4\pi G{\rho }_0/c_s^2=2/H^2$ givesSet $\tau =\tanh (z/H)$ and $\nu =kH$. Since ${\partial }_z=(1-{\tau }^2){\partial }_{\tau }/H$, this becomes the associated Legendre equation of degree one and order $\nu$:

To verify the two explicit solutions, put $\eta =z/H$ and $T=\tanh \eta$. They arebecause $e^{2\eta }=(1+\tau )/(1-\tau )$. Differentiating $F_{+}$ twice giveswhich cancels identically against $[2(1-T^2)-{\nu }^2]F_{+}$. Reflection $\eta \mapsto -\eta$ proves the same for $F_{-}$. Their Wronskian is $2\nu (1-{\nu }^2)$, evaluated at $\eta =0$; thus they are generally independent but coincide or become dependent at the exceptional orders $0,\pm 1$.

Use $n=|k|H>0$ to analyze decay, which is unaffected by the sign of $k$. The solution decaying at $+\infty$ is proportional to $e^{-n\eta }(n+\tanh \eta )$. At $-\infty$ its leading term is $(n-1)e^{n|\eta |}$, so it also decays there only if $n=1$. At this valueTherefore the marginal fragmentation mode of an isothermal slab hasAt $n=1$ the second independent solution obtained by reduction of order grows at infinity, so the coincidence of the two explicit formulas does not supply an extra admissible mode. At $k=0$ a uniform vertical translation gives another neutral disturbance, and mass-preserving relabeling displacements can have $\delta \rho =\varphi =0$. Neither is the nonzero-wavenumber fragmentation threshold just calculated.

The unstable fragmentation modes have wavelengths greater than $2\pi H$. Gas-pressure support becomes less effective relative to self-gravity at long wavelength. The side of the threshold can also be checked from the energy identity: at $|k|=1/H$ the marginal mass density disturbance has zero quadratic energy. Keep that mass density disturbance and decrease $|k|$; the pressure integral stays fixed, while the gravitational term becomes more negative because the positive inverse operator $(-{\partial }_z^2+k^2{)}^{-1}$ increases as $k^2$ decreases. Thus the quadratic energy becomes negative, giving an unstable mode on the longer-wavelength side. This argument concerns the gravitational sector and its stated decay conditions.