Past exam of the mathematics course of the University of Cambridge / 2017 / iii / Paper 321 (Dynamics of Astrophysical Discs)

Take $0<r_1<r_2$ and positive mean kinematic viscosity $\bar{\nu }$. The steady state of the supplied viscous evolution of an accretion disk equation satisfies ${\mathcal{F}}^{\prime }=-2\pi rS(r)$. The outer no-flux condition therefore givesPositive $S$ produces positive inward accretion rate in this sign convention. Set $Y=r^{1/2}\bar{\nu }\Sigma$. In a Keplerian disk, the viscous torque in an accretion disk is proportional to $Y$, so the zero-torque inner boundary condition gives $Y(r_1)=0$. Since $Y^{\prime }=\mathcal{F}/(6\pi r^{1/2})$, a second integration yieldsThis is the steady viscous disk supplied by a radial source. It does not require $\bar{\nu }$ to be independent of radius: the equation directly determines $\bar{\nu }\Sigma$. If viscosity depends on the local thermodynamic state, the displayed formula is an implicit relation completed by a separate closure. For nonnegative integrable $S$ it gives nonnegative surface density of a disk, and the inner accretion rate equals the total supplied mass per unit time, $2\pi {\int }_{r_1}^{r_2}rS(r)dr$. Differentiating both nested integrals recovers the flux and the differential equation.

Interpret the outer Dirac delta function as a boundary source carrying its full mass, or as the limit of a narrow source lying just inside $r_2$. Then ${\int }_u^{r_2}vS(v)dv=r_2{S}_0$ for $r_1<u<r_2$, andThus the interior is the usual steady Keplerian accretion disk with a constant inward accretion rate. The product $\bar{\nu }\Sigma$ vanishes at the stress-free inner edge, and the surface density of a disk does too if the viscosity has a finite nonzero boundary value; far from that edge $\bar{\nu }\Sigma$ approaches $\dot{M}/(3\pi )$. Its density need not itself be constant because the viscosity may vary.

The outer-edge source of a viscous disk requires care about the boundary trace. The flux just inside the injecting edge is $\mathcal{F}(r_2^{-})=\dot{M}$, while the no-flux value outside the source is $\mathcal{F}(r_2^{+})=0$; the source produces the intervening jump. A simultaneously continuous interior condition $\mathcal{F}(r_2^{-})=0$ would contradict nonzero steady injection. A symmetrically regularized delta integrated only up to its center can instead give half its full weight; that would mean a different supplied mass, not a different interior solution for the same $\dot{M}$. The physically normalized one-sided limiting source removes this convention ambiguity.

In the vertical equations, $\mu$ denotes dynamic viscosity, while ${\mu }_m$ is the mean molecular weight in units of the proton mass $m_p$. The symbol $k$ denotes the Boltzmann constant. For the intended hydrostatic equilibrium calculation, put $P=K{\rho }^{\gamma }$ with $K>0$, $\gamma >1$, and set $\mathcal{R}=k/({\mu }_m{m}_p)$. Dividing the hydrostatic equation by $\rho$ and integrating from the symmetric midplane givesThe specific enthalpy vanishes at the free surface. Hence the polytropic vertical structure in stellar gravity isOn $|z|<H$, the ideal gas relation gives $T=T_m(1-z^2/H^2)$ and $P=P_m(1-z^2/H^2{)}^{\gamma /(\gamma -1)}$; the vacuum extension sets $\rho =P=0$ outside. The formula needs $\gamma >1$ for a finite positive semi-thickness and assumes the harmonic vertical gravity remains valid through this height.

There is a genuine inconsistency if all four printed structure equations are imposed literally with $F=0$ and finite positive opacity: radiative diffusion then requires $dT/dz=0$, whereas the polytropic temperature above has a nonzero gradient away from the midplane. The intended result is valid for the hydrostatic polytropic equation of state with the radiative closure replaced or neglected, for example an adiabatic opaque limit. It is not an exact simultaneous solution of the stated finite-opacity radiative system. With that radiative system and a nonvacuum ideal gas, $\gamma =1$ instead gives an isothermal solutionwhich has no finite free surface. For $\gamma \ne 1$, constant $T$ and constant $K$ would force constant density, contradicting hydrostatic balance on a nontrivial vertical interval. This is the zero-radiative-flux obstruction for a polytropic ideal gas.

Use a gas-pressure-supported, optically thick alpha disk, with $\mathcal{R}=k/({\mu }_m{m}_p)$ and positive parameters. Here $\sigma$ is the Stefan–Boltzmann constant in the Stefan–Boltzmann law, and $F_{+}$ and $F_{-}$ below are heating and cooling flux estimates per disk face. Hydrostatic equilibrium and the ideal gas relation give $\mathcal{R}T\sim {\Omega }^2{H}^2$, while $\rho \sim \Sigma /H$. Integrating viscous heating over one semi-thickness and estimating radiative diffusion giveHere the stipulated negative hydrogen ion opacity exponent is $\kappa ={\kappa }_0{\rho }^{1/3}{T}^{10}$, which is a local approximation distinct from other empirical opacity fits. Substituting $T\sim {\Omega }^2{H}^2/\mathcal{R}$ and equating heating and cooling givesDropping numerical constants of order one is essential here. This is the negative-hydrogen-opacity alpha-disk thickness scaling. Its angular thickness estimate isFor a central mass $M$, Keplerian rotation gives $\Omega =(GM/r^3{)}^{1/2}$, so $H/r\propto r^{53/82}\Sigma (r{)}^{-7/41}$. The isothermal sound speed $c_s^2=\mathcal{R}T$ isAn adiabatic sound speed differs by the constant factor $\sqrt{\gamma }$ if the adiabatic index is fixed. For $\Sigma \propto r^{-p}$, constant $\alpha$, and fixed composition, the disk aspect ratio grows as $r^{(53+14p)/82}$, so there is disk flaring when $p>-53/14$; $c_s\propto r^{(6+7p)/41}$. The problem supplies no numerical parameters or specified $\Sigma (r)$, so no unique numerical $H/r$ or unconditional radial power can be inferred.

If the additional standard closure of a constant inward accretion rate well outside a zero-torque edge is intended, $\bar{\nu }\sim \alpha \Omega H^2$ and $\bar{\nu }\Sigma \simeq \dot{M}/(3\pi )$ imply $\Sigma \propto (\alpha \Omega H^2{)}^{-1}$. Eliminating it gives $H^9\propto {\Omega }^{-38/3}$ at fixed $\dot{M},\alpha$, hence $H/r\propto r^{10/9}$ and $c_s\propto r^{11/18}$. These stronger powers use that extra steady-flow assumption. All these local scalings cease to apply if the thin disk or the stated opacity regime fails.

Assume the vertical dynamical timescale of a disk is shorter than the thermal timescale of an accretion disk, which is itself shorter than the viscous timescale. This ordering is controlled in a thin disk with $\alpha \ll 1$: $t_{\mathrm{dyn}}\sim {\Omega }^{-1}$, $t_{\mathrm{th}}\sim (\alpha \Omega {)}^{-1}$ and $t_{\nu }\sim (r/H{)}^2{t}_{\mathrm{th}}$. Therefore compare heating and cooling at fixed surface density of a disk, radius, $\alpha$ and composition. Vertical hydrostatic equilibrium adjusts on the vertical dynamical timescale of a disk, giving $H\propto T^{1/2}$ and $\rho \propto T^{-1/2}$. Consequently the stipulated negative hydrogen ion opacity varies as $\kappa \propto T^{10-1/6}=T^{59/6}$. The rates per disk face then scale asThus heating increases while cooling decreases after a positive temperature perturbation. At equilibrium $Q^{+}=Q^{-}=Q_0$, a small change gives $\delta (Q^{+}-Q^{-})=(41/6)(Q_0/T)\delta T$, which drives the perturbation further from equilibrium for positive effective heat capacity. A negative temperature perturbation likewise increases the net cooling. This proves thermal instability from negative hydrogen ion opacity within the local closure; it is the temperature-slope test of thermal stability of an accretion disk. Holding volume density rather than surface density of a disk fixed would give a different cooling exponent, but still the same qualitative instability; the fixed-column comparison is the relevant thin-disk thermal test.

Let $K={\rho }_0|\mathbf{u}{|}^2/2$ be the kinetic energy per unit volume. Dot the momentum equation with ${\rho }_0\mathbf{u}$. The Coriolis acceleration does no work because $\mathbf{u}\cdot ({\mathbf{e}}_z\times \mathbf{u})=0$. Incompressibility turns the advection and pressure terms into divergences, yieldingThus $\mathbf{F}=(K+P)\mathbf{u}$ is the energy flux and $S={\rho }_0{\Omega }^2(3x-2qz)u_x$ is the source in the requested kinetic-energy form.

To decide about conservation of energy, absorb the ordinary radial shearing-sheet tidal potential ${\Phi }_t=-3{\Omega }^2{x}^2/2$ into the density. Its material derivative cancels the $3x$ work term:For $q=0$ this is the usual conservative Jacobi energy in a shearing sheet, subject to vanishing boundary energy flux. For $q\ne 0$, the extra body force has nonzero curl, $\nabla \times [-2{\Omega }^{2}qz{\mathbf{e}}_x]=-2{\Omega }^{2}q{\mathbf{e}}_y$, so no scalar potential can absorb it. A potential $2{\Omega }^{2}qxz$ would necessarily introduce a vertical force $-2{\Omega }^{2}qx{\mathbf{e}}_z$, absent from the model. Therefore this vertical-shear model does not generically conserve the standard total mechanical energy: its maintained shear/irradiation background can supply or remove energy. This is the energy balance of the incompressible vertical-shear model.

Write the proposed background as $\mathbf{U}=U(x,z){\mathbf{e}}_y$, $U=-3\Omega x/2+q\Omega z$. It satisfies incompressibility, and $\mathbf{U}\cdot \nabla \mathbf{U}=U{\partial }_y\mathbf{U}=0$. A constant pressure has zero gradient. The Coriolis acceleration is $-2\Omega {\mathbf{e}}_z\times \mathbf{U}=2\Omega U{\mathbf{e}}_x=-(3x-2qz){\Omega }^2{\mathbf{e}}_x$, which exactly cancels the prescribed body force. Thus it is a steady state.

The background is ordinary Keplerian radial shear together with vertical shear: ${\partial }_{x}U=-3\Omega /2$ and ${\partial }_{z}U=q\Omega$. At a given radius, different heights have different azimuthal velocities. The small $q$ in this question measures vertical shear; it is distinct from the orbital shear parameter, whose radial value here remains $3/2$. The vertical shear instability draws on this height-dependent differential rotation.

Set $\mathbf{k}=(k_x,0,k_z)\ne 0$. Incompressibility gives $(\mathbf{k}\cdot \tilde{\mathbf{u}})f^{\prime }(\xi )e^{st}=0$. Since $f$ is nonconstant and differentiable, it follows that $\mathbf{k}\cdot \tilde{\mathbf{u}}=0$; thereforeThe background has no advection of a $y$-independent perturbation. Its only remaining cross-advection term is $({\mathbf{u}}^{\prime }\cdot \nabla )\mathbf{U}=(-3\Omega u_x^{\prime }/2+q\Omega u_z^{\prime }){\mathbf{e}}_y$. The velocity/time terms have spatial factor $f$, while $\nabla P^{\prime }=\tilde{P}\mathbf{k}{g}^{\prime }(\xi )e^{st}$. For nonzero pressure amplitude, $g^{\prime }$ must therefore be a constant multiple of $f$; absorbing that multiple into $\tilde{P}$ gives $g^{\prime }=f$. An additive spatial constant in $g$ is a pressure gauge. If $\tilde{P}=0$, there is no such restriction on the unused pressure profile.

Matching coefficients givesBecause the quadratic perturbation term vanishes identically, these amplitudes give an exact single-phase incompressible perturbation of a shear flow, not merely a first-order truncation. Any sufficiently differentiable profile with $g^{\prime }=f$ works locally; boundary conditions may restrict the profiles in a particular domain.

Eliminate pressure by multiplying the radial equation by $k_z$ and subtracting $k_x$ times the vertical equation. For $k_z\ne 0$, incompressibility gives ${\tilde{u}}_z=-(k_x/k_z){\tilde{u}}_x$, and the pressure-free relation is $s(k_x^2+k_z^2){\tilde{u}}_x/k_z=2\Omega k_z{\tilde{u}}_y$. Substitution in the azimuthal equation yieldsThus the vertical shear instability occurs when $2q{k}_x{k}_z>k_z^2$, or $2q(k_x/k_z)>1$ for $k_z\ne 0$. Otherwise the nonzero modes oscillate; equality is marginal in the exponential-growth sense. With $q>0$, growing disturbances have $k_x/k_z>1/(2q)$ and are strongly inclined in wavevector space.

Put $R=k_x/k_z$. Maximizing $(2qR-1)/(1+R^2)$ gives $q{R}^2-R-q=0$. For $q>0$, the unstable maximizing root isMore precisely its bracket is $1-q^2/2+O(q^4)$, and $R_{\max }\sim q^{-1}$. This is the maximum growth rate of the vertical shear instability. For negative $q$, the maximum uses the corresponding negative root and is asymptotically $\Omega |q|$. At $q=0$ there is no exponentially growing mode. The undivided formula handles $k_z=0$: $s^2=0$, not an exponential instability; an allowed vertical velocity can instead force a secular azimuthal change. The excluded $\mathbf{k}=0$ is a spatially uniform disturbance and does not belong to the stipulated nonconstant one-phase family.

The razor-thin disk approximation uses vertically integrated pressure and surface density of a disk. To close the printed equations, interpret the barotropic equation of state as $P=P(\Sigma )$, or supply an equivalent fixed vertical conversion from the printed $\rho$ to $\Sigma$. Without that conversion, the literal $P=P(\rho )$ does not close this two-dimensional system. Set $c_s^2=(dP/d\Sigma {)}_0>0$, and write $\Sigma ={\Sigma }_0+\sigma$, $\mathbf{u}={\mathbf{u}}_0+{\mathbf{u}}^{\prime }$, $h^{\prime }=P^{\prime }/{\Sigma }_0=c_s^2\sigma /{\Sigma }_0$. Axisymmetric flow here means ${\partial }_y=0$, so background advection vanishes. The background shear enters through ${\mathbf{u}}^{\prime }\cdot \nabla {\mathbf{u}}_0=-3\Omega u_x^{\prime }{\mathbf{e}}_y/2$. The complete linearized equations areThese are the compressible Keplerian shearing sheet equations for an inertial-acoustic wave. No perturbed self-gravity is present in the given force law, so no gravitational potential perturbation is to be added. This closure is the barotropic closure of a razor-thin disk.

For a normal mode proportional to $e^{ikx-i\omega t}$, the linear equations giveFor $\omega \ne 0$, eliminate $U_y=-i\Omega U_x/(2\omega )$ and $\hat{\sigma }=k{\Sigma }_0{U}_x/\omega$. The density-wave dispersion relation in a non-self-gravitating disk isThe ${\Omega }^2$ term is the radial epicyclic frequency squared for Keplerian rotation: Coriolis force and background shear restore an orbital displacement. The $c_s^2{k}^2$ term is acoustic restoration by a pressure gradient. Long wavelengths are mainly epicyclic motion; short wavelengths are mainly acoustic waves. The given model contains no disk self-gravity term.

For completeness, the full three-variable determinant is proportional to $\omega ({\omega }^2-{\Omega }^2-c_s^2{k}^2)$. Its additional $\omega =0$ branch has $U_x=0$ and $2\Omega U_y=ik{c}_s^2\hat{\sigma }/{\Sigma }_0$. This is a stationary pressure/Coriolis balance, the zero-frequency balanced mode of an axisymmetric disk, rather than a propagating density wave. At $k=0$ it reduces to a uniform surface-density change. Dividing by $\omega$ is legitimate only for the wave branches just derived.

Choose the positive-frequency branch and real radial amplitude $u$. The azimuthal amplitude of the inertial-acoustic wave is $U_y=-i\Omega u/(2\omega )$, so taking real parts givesFor $|k|c_s/\Omega \ll 1$, $\omega =\Omega [1+O((k{c}_s/\Omega {)}^2)]$. Thus, to leading long-wavelength order,The quarter-cycle phase shift and the factor $1/2$ are the epicyclic motion relation in a Keplerian shearing sheet. This approximation drops the small frequency detuning; over times $|\omega -\Omega |t\gtrsim 1$ that phase approximation is not uniform. The associated fractional surface density of a disk perturbation is $\sigma /{\Sigma }_0=(ku/\omega )\cos (kx-\omega t)$, so linearization also requires $|ku/\omega |\ll 1$. The source's arbitrary real $u$ sets the wave phase/amplitude within that small-perturbation regime.

The pressureless dust fluid in a shearing sheet has no pressure restoration and is driven by drag force toward the prescribed gas velocity; gas backreaction is omitted in the given model. After linearizing about the common background shear, axisymmetry givesApply $D_{ϵ}={\partial }_t+ϵ\Omega$ to the first equation and use the second to eliminate $v_y^{\prime }$. The exact linear scalar equation isFor the supplied leading-order epicycle, set $\theta =kx-\Omega t$. Since ${\partial }_t{u}_x^{\prime }=\Omega u\sin \theta$ and $u_y^{\prime }=u\sin \theta /2$, this becomesThe final cosine term is missing from the printed equation. Dropping it is justified at leading order in $0<ϵ\ll 1$, recovering the intended forcing, but the printed equation retains an order-${ϵ}^2$ frequency term on the left and is not an exact derivation at that order. A direct countercheck is ${\mathbf{v}}^{\prime }={\mathbf{u}}^{\prime }$: the drag vanishes and the pair solves both dust equations exactly, while $v_x^{\prime }=u\cos \theta$ leaves residual ${ϵ}^2{\Omega }^{2}u\cos \theta$ in the printed scalar equation. This is the forced dust epicycle with aerodynamic drag; its damping time is $(ϵ\Omega {)}^{-1}$.

For either scalar equation, the complementary solution isso initial velocity transients decay when $t\gg (ϵ\Omega {)}^{-1}$. In the exact corrected forced dust epicycle with aerodynamic drag, $v_{x,p}^{\prime }=u\cos \theta$ is a particular solution, and the original radial equation then gives $v_{y,p}^{\prime }=u\sin \theta /2$. ConsequentlyThis is resonant dust entrainment by an epicyclic gas wave: in the adopted long-wavelength gas approximation, the forcing is a free epicycle of the dust too, so matching the gas makes the drag zero. The convergence holds for every $ϵ>0$ in this linear model; at $ϵ=0$ transients never decay, so the zero-drag and late-time limits do not commute.

If instead the printed scalar equation is treated literally for finite $ϵ$, its particular solution isIt has amplitude $|u|/\sqrt{1+{ϵ}^2/4}$ and phase offset $\arctan (ϵ/2)$, not exact phase/amplitude agreement. Agreement in that truncated equation means only leading order as $ϵ\to 0$. This distinguishes the intended approximation from a false exact reading of the printed assertion.

There is a further order-of-limits condition: a real finite-wavelength gas wave has $\omega -\Omega \ne 0$. Leading resonant entrainment over the damping time requires $|\omega -\Omega |\ll ϵ\Omega$, or $(k{c}_s/\Omega {)}^2\ll ϵ$. Otherwise solve the dust transfer equations at the actual $\omega$; detuning can reduce the dust response even though $|k|c_s/\Omega \ll 1$.

Let ${\Sigma }_{d0}$ be the uniform dust background and put ${\delta }_d={\Sigma }_d^{\prime }/{\Sigma }_{d0}$. The linearized continuity equation is ${\partial }_t{\delta }_d=-{\partial }_x{v}_x^{\prime }$ because background shear has no divergence or axisymmetric advection. Using the late entrained velocity $v_x^{\prime }=u\cos \theta$ givesThus the coherent travelling wave has maximum fractional excess $|ku|/\Omega$, and its crest density is ${\Sigma }_{d0}[1+|ku|/\Omega ]$ when the stationary part is absent. This is the density contrast of an entrained dust wave, and requires $|ku|/\Omega \ll 1$ in the present linear approximation. Its fractional oscillating density amplitude is the same as the gas's long-wavelength amplitude, so the coherent wave produces no additional dust-to-gas concentration ratio at this order.

Velocity damping alone does not erase the time-independent $C(x)$; initial density and the history of the forcing are needed to fix it. The printed data therefore determine the travelling-wave contrast, rather than a unique total late-time concentration for arbitrary initial histories. As an explicit example, switch on the prescribed gas epicycle at $t=0$ with initially uniform dust and ${\mathbf{v}}^{\prime }=0$. The corrected velocity pair gives $v_x^{\prime }=u(1-e^{-ϵ\Omega t})\cos \theta$. Integrating continuity then leavesat late times. It is smaller by order $ϵ$ in the stipulated weak-coupling regime, and confirms the stated leading wave-crest estimate while showing why arbitrary stationary density structure is not removed by the velocity transient.