Past exam of the mathematics course of the University of Cambridge / 2005 / iii / Paper 73 (Accretion Discs)

The first equation is vertical hydrostatic equilibrium: the upward pressure gradient balances the central mass's vertical gravitational force. For $|z|\ll r$, the point-mass acceleration is $-GMz/(r^2+z^2{)}^{3/2}\simeq -{\Omega }_K^{2}z$, where Keplerian rotation gives ${\Omega }_K^2=GM/r^3$. The disc's self-gravity and vertical inertia are neglected.

The second equation is local thermal balance. The divergence of the vertical radiative flux equals viscous dissipation. In the convention of this alpha disk, $\nu =\alpha p/(\rho \Omega )$, and the dominant Keplerian shear has $rd\Omega /dr=-3\Omega /2$. ThereforeThe radiative flux increases from zero at the symmetric midplane to a positive outgoing value at the upper surface. Radial heat transport, thermal advection and rapid storage of thermal energy are excluded from this local balance.

The third equation is radiative diffusion, equivalentlyIt applies to optically thick material, with $\kappa$ the opacity per unit mass and $\sigma$ the Stefan-Boltzmann constant. Here constant electron-scattering opacity represents Thomson scattering; the formula uses a diffusive, nearly thermal radiation field. Positive upward flux requires temperature to decrease upwards.

The final equation is the ideal gas pressure law, with mean molecular weight ${\mu }_m$ in proton-mass units. Writing $\mathcal{R}=k/({\mu }_m{m}_p)$ gives $p=\mathcal{R}\rho T$. The pressure is gas pressure because radiation pressure is assumed negligible. The four equations express vertical support, local heating balance, diffusive cooling and the gas equation of state.

Work at one radius and write $\mathcal{R}=k/({\mu }_m{m}_p)$. Use the surface density to fix the mass density normalization and introduce the height scaleChoose the dimensionless variablesThe pressure and temperature scales make the hydrostatic and ideal gas equations parameter-free. The flux scale fixes the heating coefficient at $9/4$, and the choice of $h$ makes the diffusion coefficient unity. Direct substitution gives the universal ordinary differential equationsFor example, the dimensionless coefficient in the temperature equation is $3\alpha \kappa {\Sigma }^2{\mathcal{R}}^4/(16\sigma {\Omega }^5{h}^6)=1$.

Let the upper free surface be at $\zeta ={\zeta }_s$, so its physical semithickness is $H=h{\zeta }_s$. Midplane symmetry and the radiative-zero disk surface conditions giveThe last condition is precisely $\Sigma ={\int }_{-H}^H\rho dz$. The unknown surface position and the two positive midplane values $P(0),\Theta (0)$ are determined as part of the normalized boundary value problem. Symmetry makes $P,D,\Theta$ even and $Q$ odd; $P^{\prime }(0)={\Theta }^{\prime }(0)=0$ already follow from the equations and are not additional independent boundary conditions.

The zero conditions do not mean zero surface flux. Integrating the heating equation gives $Q({\zeta }_s)=(9/4){\int }_0^{{\zeta }_s}Pd\zeta >0$. Near a regular surface with finite outgoing flux, dividing the hydrostatic equation by the temperature equation yields $dP/d\Theta \simeq {\zeta }_s{\Theta }^3/Q_s$. Thus $P\propto {\Theta }^4$, $D\propto {\Theta }^3$ and $\Theta \propto {\zeta }_s-\zeta$, so mass density also vanishes. The apparent singular quotient in the diffusion equation has a finite regular limit. All boundary conditions and equations are independent of the physical parameters: this is homologous vertical structure of a constant-opacity alpha disk.

The density-weighted mean kinematic viscosity isusing $\nu =\alpha p/(\rho \Omega )$. Under the preceding rescaling, letUniqueness of the universal dimensionless solution makes $I_P$ a positive numerical constant independent of $\alpha ,\Omega ,\kappa ,\Sigma ,\mathcal{R}$. HenceFor Keplerian rotation, ${\Omega }^{-2/3}=(GM{)}^{-1/3}r$. The constant-opacity gas-pressure disk viscosity law is thereforeThe particular scaling convention changes the representation of $I_P$ but not this physical law. Its nontrivial powers come from balancing viscous heating, diffusive cooling and vertical support together, rather than assuming a fixed temperature or thickness.

Estimate characteristic bulk terms. Put $ϵ=H/r\ll 1$ and assume smooth radial variation on scale $r$, secular evolution, and $\alpha \lesssim 1$. Vertical support gives $c_s^2\sim p/\rho \sim {\Omega }^2{H}^2$, while the alpha disk prescription and slow accretion implyThese velocity estimates follow from the viscous radial evolution scale and mass conservation. They are estimates for a slowly evolving thin disc, not a claim that geometric thinness excludes independently forced rapid vertical motion.

For vertical balance, the exact central force expands asAt $z\sim H$, the neglected correction has magnitude $\rho {\Omega }^{2}H{ϵ}^2$, compared with the retained force $\rho {\Omega }^{2}H$. Radial pressure support similarly changes the rotation from its Keplerian value byVertical meridional inertia scales as $u_r{u}_z/r$ or $u_z^2/H$, both $O({\alpha }^2{ϵ}^4{\Omega }^{2}H)$. Meridional viscous accelerations can scale as $\nu u_z/H^2$ or $\nu u_r/(rH)$, both $O({\alpha }^2{ϵ}^2{\Omega }^{2}H)$. Thus the gravity correction and possible viscous correction giveSome omitted terms are smaller than this leading order.

For the heating equation, the retained dissipation is $q_0^{+}\sim \rho \nu {\Omega }^2\sim \alpha p\Omega$. Radial pressure support and the finite-height correction change the dominant shear by a relative $O({ϵ}^2)$. Smooth vertical variation of the rotation gives ${\partial }_z(r\Omega )=O(ϵ\Omega )$, whose additional squared shear contributes $O({ϵ}^2{q}_0^{+})$. Meridional shear of order $u_r/H$ gives an additional relative $O({\alpha }^2{ϵ}^2)$ contribution.

Thermal advection and compressional work are of order $\rho u_r{c}_s^2/r$ and $p{u}_r/r$, respectively. Relative to $\alpha p\Omega$ they are $O({ϵ}^2)$. Thermal storage on the viscous evolution timescale $t_{\nu }\sim r^2/\nu$ has the same order, since $p/t_{\nu }\sim \alpha p\Omega {ϵ}^2$. For radiative diffusion with the same local diffusivity in different directions, smooth gradients give $F_r/F_z\sim H/r$, soConsequentlyThis is the thin-accretion-disk truncation error for the secular, radially smooth approximation. Neglect of disc self-gravity and radiation pressure, optical thickness, the opacity model and the idealized zero surface boundary have their own validity conditions; their errors do not become small solely because $H/r$ does. Near a sharp radial boundary layer or during thermal evolution on the heating timescale, the above estimates need not apply.

The variable $\psi$ is a rotating-flow pressure potential including the conservative body-force and centrifugal potentials. For constant mass density,up to an additive function of time, where $\Phi$ is any imposed conservative potential. The last term is the centrifugal potential; its negative gradient is the centrifugal acceleration. The Coriolis force remains explicit in the equation because it is velocity-dependent and cannot be absorbed into a scalar position potential.

The further constraint is incompressible flow:For the specified shear, ${\mathbf{u}}_0\cdot \nabla {\mathbf{u}}_0=0$ and ${\nabla }^2{\mathbf{u}}_0=0$. Its Coriolis acceleration is $4\Omega Ax{\mathbf{e}}_x$, so the corresponding total basic pressure potential can be chosen as ${\psi }_0=-2\Omega A{x}^2$. This supplies the force balance required for the basic state used below.

Let the perturbation phase be $\vartheta =\mathbf{k}(t)\cdot \mathbf{r}$ and write $\psi ={\psi }_0+\mathrm{Re}[\Pi (t)e^{i\vartheta }]$. Incompressible flow imposes $\mathbf{k}\cdot \mathbf{v}=0$. Applying the background advective derivative to the phase givesCancel this position-dependent phase by choosingThe perturbation acting on the basic flow gives $({\mathbf{u}}^{\prime }\cdot \nabla ){\mathbf{u}}_0=-2A{u}_x^{\prime }{\mathbf{e}}_y$. Its self-advection vanishes exactly: the physical velocity ${\mathbf{u}}^{\prime }$ is perpendicular to $\mathbf{k}$, while every spatial derivative of ${\mathbf{u}}^{\prime }$ is proportional to a component of $\mathbf{k}$. Thus ${\mathbf{u}}^{\prime }\cdot \nabla {\mathbf{u}}^{\prime }=0$, including the products between the complex mode and its conjugate. Substituting into the full Navier-Stokes equation therefore gives an exact incompressible shearing wave, with no small-amplitude approximation:In components these equations areDifferentiating the transversality constraint closes the pressure:The coefficient $4A$ includes both the shear term in the momentum equation and $\dot{\mathbf{k}}\cdot \mathbf{v}$; treating the wavevector as fixed would miss one contribution. Initial data with $\mathbf{k}(t_0)\cdot \mathbf{v}(t_0)=0$ remain transverse under this closed evolution. These shearing waves are exact in an unbounded local shear flow or with compatible shearing-periodic boundaries; arbitrary rigid boundary conditions need not admit a single such wave.

For $k_z=v_z=0$, define the complex vertical vorticity amplitudeTake the curl of the amplitude equations, including ${\dot{k}}_x=2A{k}_y$. Incompressibility cancels the pressure and rotation terms, givingEquivalently, the uniform basic absolute vorticity has no gradient for the perturbation to advect. The general two-dimensional viscous shearing wave solution iswhere $Z_0$ is an arbitrary complex constant and $k_x(t_0)$ is an arbitrary real initial radial wavenumber. The perturbation pressure follows from part (ii), or explicitly $\Pi =(4A{k}_y^2-2\Omega k^2)Z/k^4$.

For $A\ne 0$, shift the time origin to the swing $k_x=0$ and put $T=k_x/k_y=2At$, $T_0=2A{t}_0$ and $\mathrm{Re}=A/(\nu k_y^2)$. ThenThe spatially averaged perturbation kinetic-energy mass density is $E=\rho ⟨|{\mathbf{u}}^{\prime }{|}^2⟩/2=\rho |\mathbf{v}{|}^2/4$, and $|\mathbf{v}{|}^2=|Z{|}^2/k^2$. ConsequentlyThe proportionality constant depends on the initial data; the exponential profile is not a claim of infinite physical growth from the remote past. For $A>0$, leading waves with $T<0$ can gain velocity amplitude as their wavevector magnitude decreases, before viscosity and later winding damp them. This is the Orr mechanism, and it is independent of uniform rotation for these two-dimensional disturbances.

If $A=0$, the wavevector is constant and $\mathbf{v}(t)=\mathbf{v}(t_0)e^{-\nu k^2(t-t_0)}$, so there is no shear amplification. The formulas with the physical-time integral remain valid, although $T=2At$ is no longer a useful coordinate. The high-positive-Reynolds-number case in the next part takes $A>0$.

For a fixed $\mathrm{Re}\gg 1$, optimize the gain over the initial wave tilt and a later observation time. From the energy profile,Its stationary points satisfyThere are two negative roots in the high-Reynolds-number regime: an early energy minimum and a later maximum. Write them as $T_{-}$ and $T_{+}$. Balancing the large-$|T|$ terms for the early root and the small-$|T|$ terms for the late root givesThus the optimum disturbance starts as a sufficiently tightly wound leading wave and reaches its peak slightly before $k_x=0$; finite viscosity shifts the peak away from the inviscid swing.

Normalize its initial energy at $T_{-}$. The exact maximum ratio isThe prefactor is asymptotic to $(2\mathrm{Re}{)}^{2/3}$, while $g(T_{-})/\mathrm{Re}\to -2/3$ and $g(T_{+})/\mathrm{Re}\to 0$. Therefore optimal transient amplification of a viscous shearing wave givesThe time between the optimum initial tilt and the peak is approximately $(2\mathrm{Re}{)}^{1/3}/(2A)$. Without viscosity, ever larger initial tilts would allow unbounded gain; viscosity damps them before they can swing, selecting the finite optimum. The coefficient above uses the paper's Reynolds number based on $A$, not one based on the full shear rate $2A$.

For axisymmetric disturbances, $k_y=0$, the wavevector is fixed, so there is no winding/swing mechanism. Take $k_z\ne 0$ and set $\chi =k_z^2/k^2$. Eliminating pressure with incompressibility givesThe axisymmetric incompressible epicyclic wave has inviscid frequencywhere ${\kappa }_r$ is the radial epicyclic frequency.

In a Keplerian flow $A=3\Omega /4$ and ${\kappa }_r=\Omega$. The axisymmetric disturbances therefore execute stable inertial/epicyclic oscillations, with viscous amplitude decay $e^{-\nu k^{2}t}$. Rotation opposes sustained radial displacement through angular-momentum conservation. There can be bounded exchanges of kinetic energy between components, but no large Reynolds-number-dependent amplification of this axisymmetric family. Indeed, with $w=v_x/\sqrt{\chi }$, inviscid motion conservesSince the physical velocity norm is $|w{|}^2+|v_y{|}^2$, Keplerian energy amplification is at most $\Omega /(\Omega -A)=4$; viscosity can only reduce that bound. A pure initial azimuthal disturbance can attain the inviscid factor four after a quarter epicyclic period. By contrast, the two-dimensional nonaxisymmetric Orr mechanism remains available in a Keplerian flow and has the ${\mathrm{Re}}^{2/3}$ gain derived above.

In non-rotating shear, $\Omega =0$ while $A$ remains nonzero. The axisymmetric solution isCross-stream motion transports the basic shear and generates an azimuthal streak: this is the lift-up effect. Inviscid azimuthal velocity grows linearly and its energy quadratically, rather than oscillating. With viscosity, the gain for a suitable initial cross-stream disturbance scales as $[A/(\nu k^2){]}^2$ at large Reynolds number, appreciably stronger than the two-dimensional ${\mathrm{Re}}^{2/3}$ scaling. For instance, $C_y=0$ and fixed $\chi >0$ give the leading optimized gain $4\chi A^2/(e^2{\nu }^2{k}^4)$ at $\tau \simeq 1/(\nu k^2)$. Neither mechanism is an exponentially growing hydrodynamic normal mode. Rotation bounds axisymmetric lift-up in Keplerian shear; without rotation, lift-up gives strong algebraic transient growth. If $k_z=0$ as well as $k_y=0$, incompressibility forces $v_x=0$ and removes this coupling altogether.

Consider a local shearing sheet centered at a circular orbit. Take uniform mass density, background velocity ${\mathbf{u}}_0=-2Ax{\mathbf{e}}_y$, rotation $\mathbf{\Omega }=\Omega {\mathbf{e}}_z$ and a weak uniform vertical magnetic field ${\mathbf{B}}_0=B_0{\mathbf{e}}_z$. HereUse ideal magnetohydrodynamics, neglect viscosity, magnetic diffusion and vertical stratification, and take incompressible axisymmetric disturbances proportional to $e^{ikz-i\omega t}$. They have horizontal velocity and field perturbations. Let $\mathbf{b}=\delta \mathbf{B}/\sqrt{{\mu }_0\rho }$ be the magnetic perturbation in velocity units, $v_A=B_0/\sqrt{{\mu }_0\rho }$ the Alfvén speed and ${\omega }_A=k{v}_A$ the Alfvén frequency. Magnetic pressure can be included in the perturbation's total pressure; the horizontal Lorentz force for this mode is the magnetic tension term $i{\omega }_A\mathbf{b}$.

Linearizing momentum and induction givesThe term $-2A{b}_x$ winds a radial magnetic perturbation into an azimuthal one. Omitting it would remove the shear's magnetic energy source and give the wrong stability criterion.

A useful elimination that also covers neutral modes introduces the fluid displacement $\mathbf{\xi }$. Frozen-in flux gives $\mathbf{b}=i{\omega }_A\mathbf{\xi }$, while the Eulerian perturbation velocity is $u_x={\dot{\xi }}_x$, $u_y={\dot{\xi }}_y+2A{\xi }_x$, because the displacement samples the gradient of the basic flow. The equations become the spring model of magnetorotational instability:For the stated time dependence, their determinant isThus the ideal magnetorotational dispersion relation isIt follows without dividing by $\omega$, so the neutral stability boundary is included.

To analyze it, put $a={\omega }_A^2$ and $Y={\omega }^2$. ThenFor a hydrodynamically stable rotation law, ${\kappa }_r^2>0$, the sum of the two roots is positive. Their product is negative precisely when $0<a<4\Omega A$. One squared frequency is then negative, giving an exponentially growing mode $\omega =i\gamma$ and its decaying partner. Therefore the magnetorotational instability, or MRI, has criterionAngular velocity must decrease outwards, although specific angular momentum may increase outwards and satisfy Rayleigh's circulation criterion. This is why a Keplerian disc can be stable to the axisymmetric hydrodynamic disturbances of Question 2 yet unstable once magnetic coupling is introduced. At zero magnetic frequency the unstable root becomes neutral, while sufficiently large magnetic frequency restores stability. For increasing angular velocity, $A<0$ with $\Omega >0$, this ideal vertical-field instability is absent.

The growth rate isFor $0<A<\Omega$, differentiation with respect to $a$ sets $\sqrt{{\kappa }_r^4+16{\Omega }^{2}a}=4{\Omega }^2$, yielding the maximum growth rate of ideal magnetorotational instability:For Keplerian rotation, $A=3\Omega /4$, soThus the fastest growth is on an orbital timescale. The field strength selects its wavelength, $k_{\max }\propto 1/v_A$, rather than the maximum ideal growth rate. At fixed $k$, growth does approach zero as $B_0\to 0$; retaining the finite maximum in that limit requires shorter wavelengths, which magnetic diffusion can eventually exclude. Unstable wavelengths satisfy $\lambda >2\pi v_A/\sqrt{4\Omega A}$ in this local model. A finite disc must accommodate such a wavelength, so a sufficiently strong field can suppress the available vertical modes through excessive tension.

Physically, a bent field line acts like a spring joining neighboring fluid elements. The inner element rotates faster and pulls the outer element forward, transferring angular momentum outward. The inner element loses angular momentum and moves inward; the outer element gains it and moves outward. In a rotation law with outward-decreasing angular velocity, these motions increase their angular-velocity difference and the magnetic coupling produces further separation. A weak restoring spring therefore creates a feedback instability by enabling angular-momentum exchange. Without that exchange, each displaced element keeps its angular momentum and a positive epicyclic frequency restores its orbit. Strong magnetic tension instead prevents the separation. The free energy is the differential rotation, not the initial magnetic energy.

The instability provides an important route to turbulence and angular momentum transport in sufficiently conducting accretion discs. The radial flux of azimuthal momentum includes the stressThe first term is a Reynolds stress and the second is the magnetic contribution from the Maxwell stress tensor, with the sign appropriate to momentum flux. For the growing branch, take $u_x$ real and positive. The linear equations give $u_y/u_x=({\gamma }^2+a)/(2\Omega \gamma )>0$ and $b_y/b_x=({\gamma }^2+a-4\Omega A)/(2\Omega \gamma )<0$. The last sign follows from $({\gamma }^2+a)({\gamma }^2+a-4\Omega A)+4{\Omega }^2{\gamma }^2=0$. Thus both the velocity and magnetic contributions carry angular momentum outward. The associated energy extraction is $-T_{xy}rd\Omega /dr=2A{T}_{xy}>0$. Inward accretion can then accompany outward angular-momentum transport, and dissipation of the resulting fluctuations heats the disc. This gives a physical mechanism behind an effective turbulent-viscosity description such as an alpha disk, whereas microscopic viscosity alone is generally too weak for rapid accretion.

MRI is an exponentially growing magnetic instability of outward-decreasing angular velocity; hydrodynamic shearing-wave amplification is transient and does not by itself establish sustained turbulence. Nonlinear saturation, field regeneration and transport efficiency depend on field geometry, boundaries and thermodynamics. Poor conductivity or important nonideal magnetic effects can suppress or modify the ideal instability, so the local dispersion relation is a mechanism and growth criterion rather than a universal prescription for the saturated stress.