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

Write the radial velocity as $u(r)$, allowing $u<0$ for accretion flow, and let $c^2=\gamma p/\rho$ be the squared adiabatic sound speed. The steady continuity equation, radial Euler equations for an inviscid fluid and entropy advection equation giveOn a nonzero smooth flow branch, the polytropic equation of state is $p=K{\rho }^{\gamma }$ with constant $K$. Thus the two useful first integrals for spherically symmetric adiabatic flow areHere $\mathcal{J}$ is the signed mass flux and $\mathcal{B}$ is the Bernoulli function. The specific enthalpy is $c^2/(\gamma -1)$; the prescribed Newtonian gravitational potential has no contribution from the gas's own gravity. Eliminating the density derivative using mass conservation gives the differential formThis sonic point equation displays both the singular coefficient and the numerator that must vanish for a smooth transonic branch.

At a finite smooth sonic point $r_{\ast }$, the equation from (a) requiresSet $\delta =\gamma -1$ and $q=2/\beta -1/2$. Substitution into the Bernoulli function givesAn inflow from a warm reservoir at infinity has $\mathcal{B}=c_{\infty }^2/\delta >0$; an outflow reaching infinity has $\mathcal{B}=u_{\infty }^2/2+c_{\infty }^2/\delta \ge 0$. Therefore a nondegenerate transonic spherical flow in a power-law potential requiresFor the usual range $0<\beta <4$, the critical adiabatic exponent for spherical power-law flow is consequentlyThe printed assumption $\beta >0$ also permits $\beta \ge 4$. In that range the boxed inequality holds for every finite $\gamma >1$: there is no finite upper bound. Equivalently, take $f(\beta )=+\infty$ for $\beta \ge 4$. Extending the rational expression beyond $\beta =4$ would give an incorrect restriction.

The strict inequality also follows from local sonic-point slope discriminant analysis, including the potentially cold zero-energy endpoint. Put $X=(r_{\ast }/u_{\ast })u_{\ast }^{\prime }$. Since mass conservation implies $(c^2{)}^{\prime }/c^2=-\delta (u^{\prime }/u+2/r)$, differentiating the sonic point equation at $r_{\ast }$ yieldsIts discriminant is $8[2\beta -(4-\beta )\delta ]$. Moreover, the logarithmic Mach number derivative at the sonic point isA genuine crossing has two distinct branches and a nonzero derivative. At equality the crossing degenerates. In particular, for $\gamma =f(\beta )$ and $\mathcal{B}=0$, mass conservation and the polytropic equation of state make the Mach number constant along the scale-invariant solution; a solution that is sonic there is sonic everywhere, rather than crossing an isolated sonic point.

Use the positive mass accretion rate $\dot{M}=-\mathcal{J}$, and retain $\delta =\gamma -1$, $q=2/\beta -1/2$. Define $h=1-q\delta$, which is positive under the condition in (b). Matching the Bernoulli function to the reservoir and using the polytropic equation of state givesThe transonic spherical accretion rate in a power-law potential is thereforeThe combination $A/c_{\infty }^2$ has dimensions of length to the power $\beta$, so this expression has dimensions of mass per time. The sonic point selects the flux that connects the subsonic reservoir to the inward supersonic transonic branch.

For the endpoint limits of power-law spherical accretion, hold $A,\beta ,{\rho }_{\infty },c_{\infty }$ fixed. As $\delta \to 0$, $\ln h=-q\delta +O({\delta }^2)$, so $(q-1/\delta )\ln h\to q$. HenceThis is also obtained directly from an isothermal equation of state: the Bernoulli function becomes $u^2/2+c_{\infty }^2\ln (\rho /{\rho }_{\infty })-A/r^{\beta }=0$, so ${\rho }_{\ast }/{\rho }_{\infty }=e^q$. For $\beta =1$, $A=GM$, it recovers the isothermal Bondi accretion rate.

For $0<\beta <4$, $q>0$ and $\gamma \uparrow f(\beta )$ means $h\downarrow 0$. Since $\delta =(1-h)/q$,ThusThe finite limiting flux accompanies $r_{\ast }\to 0$ and $c_{\ast },{\rho }_{\ast }\to \infty$; it does not assert a finite-radius sonic point at the endpoint. For $\beta =1$ this is the familiar $\gamma \uparrow 5/3$ limit $\pi G^2{M}^2{\rho }_{\infty }/c_{\infty }^3$.

For completeness, the printed range $\beta \ge 4$ has the endpoint $f=+\infty$. At $\beta =4$, $q=0$ and the flux is independent of $\gamma$, namely $4\pi {\rho }_{\infty }(2A{)}^{1/2}$. For $\beta >4$, $q<0$ and $h\sim |q|\delta$, so $\dot{M}\to 0$ as $\gamma \to \infty$. These are the corresponding extended endpoint limits.

The solenoidal constraint on the magnetic field reduces to ${\partial }_x{B}_x+{\partial }_z{B}_z=0$. Locally, or globally in a simply connected cross-section, it permits a Cartesian magnetic flux function satisfyingThe transverse magnetic field lines are contours of $\psi$, because $\mathbf{B}\cdot \nabla \psi =0$. The independent component $B_y$ supplies the twist of the flux tube; it is not restricted by the solenoidal constraint.

For the Lorentz force, take the curl explicitly:Its cross product with $\mathbf{B}$ has transverse components $-({\nabla }^2\psi )\nabla \psi -B_y\nabla B_y$ and $y$ component ${\psi }_x{B}_{y,z}-{\psi }_z{B}_{y,x}$. ConsequentlyThe term $-B_y\nabla B_y/{\mu }_0$ is the transverse gradient of the magnetic pressure associated with the axial field, while the other terms include magnetic tension.

In magnetostatics, $-\nabla p-\rho \nabla \Phi +{\mathbf{f}}_L=0$. The plane-parallel Newtonian gravitational potential and all thermodynamic fields are independent of $y$, so the $y$ component of the Lorentz force must vanish:Away from a null of $\nabla \psi$, this says that $B_y$ is constant along each connected contour of the Cartesian magnetic flux function. Thus, on a regular flux region,This statement is local; disconnected contours with the same numerical flux label need not share one global function without the usual flux-tube connectivity assumption.

Since $\nabla B_y=b^{\prime }(\psi )\nabla \psi$, the transverse magnetostatic equilibrium equation becomesThis is a Cartesian magnetostatic flux-function equilibrium. Projecting it along a transverse magnetic field line gives hydrostatic balance along that line; projecting across the line balances the gas pressure and weight against magnetic pressure and magnetic tension.

The MHD induction equation can be written in material form asIts $y$ component immediately givesFor the Cartesian magnetic flux function, the $x,z$ components of the MHD induction equation, or the $y$ component of the magnetic vector potential equation, giveThe additive function of time in $\psi$ has no effect on $\mathbf{B}$. Choose this gauge so that $C=0$. Then the flux label is materially conserved:This is magnetic flux freezing in the two-dimensional geometry.

To prove the absence of axial motion, construct an invariant solution with $u_y=0$. Along each Lagrangian trajectory, $\psi$ is fixed, and the assumed divergence givesThus $B_y$ remains a function of $\psi$ and time alone. Its gradient stays parallel to $\nabla \psi$, so the $y$ component of the Lorentz force remains zero. The axial momentum equation isbecause $(B_x{\partial }_x+B_z{\partial }_z)\psi =0$. With $u_y=0$ initially, it remains zero. Hence neither axial force nor axial motion is generated. This flux-surface preservation during magnetic-tube expansion argument assumes the smooth ideal evolution for which the initial-value problem is unique; the particular in-plane rising motion need not be calculated.

Use $D={\partial }_t+w{\partial }_z$, $w=v_z$, and a prime for ${\partial }_z$. For the plane-parallel magnetohydrodynamic flow with imposed shear, $\nabla \cdot \mathbf{u}=w^{\prime }$, whileThe solenoidal magnetic field constraint gives $B_z^{\prime }=0$. The $z$ component of the MHD induction equation then gives ${\partial }_t{B}_z=0$. Thus $B_z$ is constant in both space and time.

Use the magnetic tension and magnetic pressure decompositionThe continuity equation and momentum equation now reduce toThe transverse components of the MHD induction equation areThe terms involving $a$ respectively accelerate the flow across the imposed shear flow and wind the transverse magnetic field. The closure is the isothermal equation of state $p=c_s^2\rho$, as printed in the PDF.

The logarithmic density term is the barotropic fluid energy for $p=c_s^2\rho$, rather than the thermodynamic heat content of an isolated gas. DefineFor any scalar $q$, the continuity equation implies ${\partial }_t(\rho q)+{\partial }_z(\rho wq)=\rho Dq$. Therefore $De=-c_s^2{w}^{\prime }$, $D\Phi =wg$, and dotting the momentum equations with $\rho \mathbf{v}$ givesThe MHD induction equation gives the complementary magnetic energy balanceAdding these two identities, and collecting the product derivatives, yields the barotropic magnetic energy equationHere $E=E_h+E_B$. The two stress terms in $S$ describe magnetohydrodynamic shear work: the maintained background shear flow can supply energy to, or remove energy from, the perturbation flow and magnetic field. The sign depends on the off-diagonal total stress; this is why the energy excluding the background shear flow is not generally conserved. Changing ${\rho }_0$ merely adds a constant multiple of the conserved density to $E$ and the corresponding mass flux to $F$, leaving $S$ unchanged.

For a steady flow, mass conservation gives $\rho w=\text{constant}$ and ${\rho }^{\prime }/\rho =-w^{\prime }/w$. Define the Alfvén velocity ${\mathbf{v}}_a=\mathbf{B}/\sqrt{{\mu }_0\rho }$ and $C=w^2-v_{az}^2$. Combining the horizontal momentum and MHD induction equations without dividing by $C$ givesDot these identities with $B_x/({\mu }_0\rho )$ and $B_y/({\mu }_0\rho )$, respectively:The vertical momentum equation and the isothermal equation of state giveMultiply by $C$ and eliminate the magnetic derivative to obtainNo division by $w^2-v_{az}^2$ was needed in deriving this necessary relation.

The magnetosonic critical speeds in the $z$ direction areThe plus sign gives the fast magnetosonic wave speed and the minus sign the slow magnetosonic wave speed. Thus the differential coefficient is $(w^2-c_{\mathrm{f}}^2)(w^2-c_{\mathrm{s}}^2)$. A smooth outflow proceeding from below both speeds to above both must normally pass through both magnetosonic critical speeds. At each crossing, the right-hand side must also vanish; this is the regularity at a magnetosonic point condition. The derivative coefficient changes sign at each nondegenerate crossing. For $w>0$ and $w^{\prime }>0$, the driving term must be positive below the slow point, negative between the points and positive above the fast point.

The Alfvén speed component satisfies $c_{\mathrm{s}}^2\le v_{az}^2\le c_{\mathrm{f}}^2$. With downward gravity $g>0$, the gravitational term is consequently nonnegative at the slow point and nonpositive at the fast point. The magnetohydrodynamic shear work contribution must balance it at each point, and can provide the upward driving needed to pass the fast point. For $a=0$ and strictly positive $g$, generic separated slow and fast points cannot satisfy the required zero numerator; special vanishing-gravity or coincident-speed cases need separate treatment.

There is also Alfvén-point compatibility in a plane-parallel sheared flow. At $w^2=v_{az}^2$, the original transverse equations requireThese restrictions are not generally visible as a zero of the scalar differential coefficient, which there equals $-v_{az}^2(v_{ax}^2+v_{ay}^2)$. In particular, a strictly accelerating regular solution must have $B_x=0$ at that point. The scalar relation is therefore a necessary wind equation, not a substitute for regularity of all the original ideal magnetohydrodynamic equations. Degenerate cases such as a purely longitudinal magnetic field can merge characteristic speeds and reduce the number of distinct critical conditions.

Use the fluid displacement $\mathbf{\xi }$, with velocity perturbation $\delta \mathbf{u}={\partial }_t\mathbf{\xi }$. Denote Eulerian perturbations by $\delta$ and the corresponding Lagrangian pressure perturbation by ${\Delta }_{L}p=\delta p+\mathbf{\xi }\cdot \nabla p$. For adiabatic perturbations, the linearized mass conservation, momentum, adiabatic equation of state and Poisson equation areThese self-gravitating adiabatic displacement equations retain the perturbation of the star's own Newtonian gravitational potential. In the uniform-density interior, $\delta \rho =-\rho \nabla \cdot \mathbf{\xi }$, $\nabla \Phi ={\omega }_d^2\mathbf{r}$, and $\nabla p=-\rho {\omega }_d^2\mathbf{r}$. The dynamical frequency of a uniform-density star ${\omega }_d=(GM/R^3{)}^{1/2}$ sets the natural timescale. For a normal mode with time factor $e^{-i\omega t}$, replace ${\partial }_t^2$ by $-{\omega }^2$. The requested interior analysis needs no surface or exterior matching conditions.

The spatial factor $r^l{Y}_l^m$ is a regular solid harmonic. Besides the Laplace equation ${\nabla }^{2}F=0$, homogeneity gives $\mathbf{r}\cdot \nabla F=lF$. The divergence of the fluid displacement isDefine the scalar dilation amplitude $D_r=r{U}^{\prime }+(l+3)U+(l/r)V^{\prime }$; this is the quantity denoted $\Delta$ in the question, not the notation ${\Delta }_L$ for a Lagrangian pressure perturbation. Then mass conservation gives $\hat{\rho }=-\rho D_r$.

The equilibrium pressure gradient and homogeneity identity giveConsequently the adiabatic equation of state yieldsFor the force equations, the product rule givesEquating the coefficients of $F{\mathbf{e}}_r$ and $\nabla F$ in the self-gravitating adiabatic displacement equations givesFinally, applying the Laplacian to the gravitational perturbation givesso the Poisson equation becomesThese are the required interior equations for uniform-density stellar oscillation.

There is a radial degeneracy at $l=0$: $F$ is spatially constant and $\nabla F=0$, so the displacement is independent of $V$. The second force equation then cannot be inferred by equating independent vectors. For nonzero $\omega$ it may be imposed as an auxiliary definition of $V$, but it is not an additional physical radial equation. At zero frequency the radial equations should be used directly. This distinction matters for interpreting the zero factor in (c).

For the polynomial stellar-mode coefficient reduction, write the highest coefficients asSet $x={\omega }^2/{\omega }_d^2$, $L=n(2l+n+1)$ and $k=l+n+1$. The highest powers in the equations of (b) giveThe factor $3$ follows from $4\pi G\rho =3{\omega }_d^2$; the positive $\gamma d/2$ term comes from the highest, negative $r^2$ coefficient of the equilibrium pressure.

For $l>0$, the final spectrum found below has $x\ne 0$. Eliminating $u,v$ first givesThe remaining pressure relation isA nontrivial leading coefficient has $d\ne 0$, soThusor, in dimensional form,

For the printed stability continuation, define $H=\gamma L/2-4$. The two branches for nonradial stellar oscillations areFor every $l>0$, their product is $-l(l+1)<0$. Therefore one branch oscillates and the other grows exponentially, for every positive $\gamma$. With the chosen time dependence, $x_{-}<0$ gives $\omega =i{\omega }_d\sqrt{-x_{-}}$ and growth rate ${\omega }_d\sqrt{-x_{-}}$. Increasing $\gamma$ does not remove this nonradial instability. At large $n$ or large $\gamma$, $x_{+}\sim \gamma L/2-4$ is a stiff compressive branch, while $x_{-}\sim -l(l+1)/(\gamma L/2-4)$ is a slower unstable buoyancy branch.

The buoyancy frequency supplies the physical explanation. The background density has zero gradient, but pressure decreases outward:The specific entropy of the perfect gas decreases outward, because $p/{\rho }^{\gamma }$ does. An adiabatically displaced fluid parcel therefore experiences destabilizing buoyancy. The model has convective instability of a uniform-density star, even when its radial compression modes are stable. The divergence of this expression at the zero-pressure surface also reflects the artificial equilibrium; the interior negative sign is already decisive.

For $l=0$, the formal quadratic factors asThe physical radial displacement has no $V\nabla F$ component. Using its equations directly givesIndeed $d=(n+1)u$, $j=-3u/n$, and $h=u+\gamma (n+1)u/2$ yield this result directly in the radial momentum equation. The extra identically zero factor in the two-component determinant need not represent an independent radial mode, since $V$ is redundant when $l=0$.

The lowest allowed degree, $n=2$, has ${\omega }_{\mathrm{radial}}^2=(3\gamma -4){\omega }_d^2$. Hence all these radial modes are stable for $\gamma >4/3$; the fundamental radial mode is marginal at $\gamma =4/3$ and unstable below it. Higher $n$ have their individual thresholds $\gamma =8/[n(n+1)]$. These radial stability of a uniform-density star thresholds coexist with nonradial convective instability. They are conclusions within the polynomial interior family assumed in the question; no boundary-condition quantization is being asserted.