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

Write $\Sigma ={\Sigma }_0{x}^{-5/2}g$ because the prescribed kinematic viscosity is $\bar{\nu }={\bar{\nu }}_0{x}^{5/2}$. Substituting $r=r_{0}x$ and ${\partial }_t=(3{\bar{\nu }}_0/4r_0^2){\partial }_{\tau }$ into the viscous evolution of an accretion disk givesHence the dimensionless diffusion equation isFor the useful inverse-radius coordinate $y=x^{-1/2}$, this becomes simply $g_{\tau }=y{g}_{yy}$. The positive factor in $\tau =1+3{\bar{\nu }}_{0}t/(4r_0^2)$ makes increasing $\tau$ correspond to increasing physical time.

There is a distinction between the named torque variable and the physical viscous torque in an accretion disk. For Keplerian rotation, the outward torque magnitude isHere $M_{\star }$ denotes the central mass. The zero-torque inner boundary condition therefore requires $x^{1/2}g\to 0$ as $x\to 0$. Requiring $g\to 0$ is stronger. The nontrivial finite-mass profile used below satisfies both; this distinction matters when counting integration constants.

Let $g=G(\xi )$ with $\xi ={\tau }^{-1}{x}^{-1/2}$. The derivatives of this similarity solution areThe first-derivative terms cancel in the dimensionless diffusion equation, leaving $G^{\prime \prime }+G^{\prime }=0$ for $\xi >0$. ThusAt the inner boundary $x\to 0$, $\xi \to \infty$. Interpreting the named variable as vanishing there sets $A=0$. Even if only the true viscous torque in an accretion disk $x^{1/2}g$ is required to vanish, a nonzero $A$ would give infinite total mass near $r=0$; finite mass again selects $A=0$. With nonnegative surface density, choose $B=C>0$; its value fixes the initial mass normalization. ThereforeThe surface density vanishes exponentially at the origin and has an $x^{-5/2}$ tail at infinity. Its peak occurs at $x_{\max }=1/(25{\tau }^2)$, with ${\Sigma }_{\max }/{\Sigma }_0=C(5\tau {)}^5{e}^{-5}$. The peak moves inward and grows as time increases.

The plot uses $C=1$ and shows the actual surface density, with logarithmic radius to resolve both peaks. The later profile is larger at every fixed radius. Its total mass is finite; its total angular momentum is not, as the next calculation establishes.

For the finite-mass similarity solution above, the total disk mass isSubstitute $\xi =1/(\tau \sqrt{x})$, so $x={\tau }^{-2}{\xi }^{-2}$ and $x^{-3/2}|dx|=2\tau d\xi$. This givesThe requested linear increase of mass is therefore correct. It is supplied by the infinite-radius boundary, rather than by a source-free isolated finite disk. The inward mass accretion rate iswhich tends to zero at the origin and to $3\pi C{\bar{\nu }}_0{\Sigma }_0$ at infinity.

For Keplerian rotation, the specific angular momentum is $\ell =\sqrt{G{M}_{\star }r}$. The total disk angular momentum would beThis integral diverges at $\xi =0$, corresponding to infinite radius. There is no finite constant total angular momentum for the nonzero printed similarity solution. The seemingly time-independent last expression is an infinite quantity; it cannot prove a finite conservation law.

The problem remains visible under a cutoff. With $r\le r_{0}R$, write $J_R=J_{\ast }{E}_1(1/(\tau \sqrt{R}))$, where $J_{\ast }=4\pi C{\Sigma }_0\sqrt{G{M}_{\star }}{r}_0^{5/2}$ and $E_1$ is the exponential integral. ThenEven the finite change after subtracting a common logarithmic divergence is not zero. Direct integration of the disk diffusion equation gives the boundary identityfor this profile the outer expression has a nonzero limit, giving the same finite change rate. An inner zero torque alone does not impose zero outer angular-momentum flux.

More generally, the printed family $A+B{e}^{-\xi }$ can have finite angular momentum at both ends only if $A=0$ and $A+B=0$, hence it is trivial. Thus the request for a nonzero finite constant $J_d$ is incompatible with the stated viscosity exponent, similarity ansatz and infinite radial domain. A physical finite-domain or altered similarity problem would need different boundary data or assumptions; silently assigning a finite value to the divergent integral does not repair it.

The four equations describe, respectively, vertical hydrostatic equilibrium in the central object's gravity, the vertical divergence of radiative flux balancing viscous dissipation, optically thick radiative diffusion, and the equation of state combining ideal gas gas pressure with radiation pressure. Here $\mu$ is dynamic viscosity, not the mean molecular weight ${\mu }_m$; $\kappa$ is the opacity, $\sigma$ the Stefan-Boltzmann constant, and $\Omega$ the local orbital frequency.

In the radiation-pressure limit, $p=4\sigma T^4/(3c)$, so radiative diffusion givesCompare this with vertical hydrostatic equilibrium. Wherever $\rho \ne 0$,At fixed radius, constant opacity makes $dF/dz=c{\Omega }^2/\kappa$. The viscous dissipation equation therefore yieldsindependent of height within the optically thick radiation-supported layer. The result assumes the supplied gravity, heating and diffusion equations; optically thin surface corrections are outside that approximation.

Use characteristic, vertically averaged quantities at fixed radius. Let $H$ be a disk thickness, $\rho \sim \Sigma /H$, and let the alpha disk prescription give $\mu \sim \alpha p/\Omega$, up to a convention-dependent constant. From the preceding radiation pressure and constant-opacity balance, $\mu \sim c/\kappa$ is independent of $\Sigma$; thusVertical hydrostatic equilibrium gives $p\sim \rho {\Omega }^2{H}^2\sim \Sigma {\Omega }^{2}H$, henceThe vertically averaged kinematic viscosity satisfies $\bar{\nu }\Sigma =\int \mu dz\sim \mu H$. ConsequentlyAt fixed radius the omitted factors are independent of $\Sigma$, so this is the requested inverse-density scaling of the viscous torque in an accretion disk.

This is an order-of-magnitude vertical model, as requested. Imposing both exact pointwise $\mu =\text{constant}$ and exact pointwise $\mu \propto p(z)$ would force $p$ to be constant, contrary to nontrivial vertical hydrostatic equilibrium. The characteristic-pressure prescription supplies a scaling, not an exact simultaneous vertical solution to all those idealized equations.

Write $W(\Sigma )=\bar{\nu }\Sigma$. In a small radial patch $r=r_c+\chi$ with disturbance scale much smaller than $r_c$, the leading disk diffusion equation is ${\partial }_t\Sigma \simeq 3{\partial }_{\chi }^{2}W$. Taking the linearization around a uniform surface density ${\Sigma }_c$ givesFor a local Fourier mode $\delta \Sigma \propto e^{st+ik\chi }$,because $W\propto {\Sigma }^{-1}$ has $W^{\prime }<0$. Thus the effective density diffusion coefficient is negative and the disk has a viscous instability of an accretion disk.

A localized overdensity is a maximum of $\Sigma$ but a minimum of $W$. At its center ${\partial }_{\chi }^{2}W>0$, so ${\partial }_t\Sigma >0$: it grows. The local radial mass flux is proportional to $-{\partial }_{\chi }W$, directing matter toward that minimum from both sides. This is the physical feedback: an overdense patch has a smaller viscous torque in an accretion disk and attracts further accumulation. The arbitrarily rapid $k\to \infty$ growth in this ideal negative-diffusion model is not a resolved prediction below the validity scale of the thin-disk and local viscous closure.

Let $\mathbf{u}=-(3/2)\Omega x{\mathbf{e}}_y+\mathbf{v}$ and $\rho ={\rho }_0(1+\delta )$, with $\delta$ and $\mathbf{v}$ small. Define the sound speed by $c^2=(dp/d\rho {)}_0$, assumed positive. The background Keplerian shear balances the radial Coriolis force and the shearing-sheet tidal potential. For disturbances preserving axisymmetry ${\partial }_y=0$, its advective derivative ${\mathbf{u}}_0\cdot \nabla$ vanishes, but $\mathbf{v}\cdot \nabla {\mathbf{u}}_0=-(3/2)\Omega v_x{\mathbf{e}}_y$ remains. The linearized Euler equations and continuity equation areThe azimuthal coefficient is $(2-3/2)\Omega =\Omega /2$, not $2\Omega$: it combines the background velocity gradient with the Coriolis force. The uniform background has no buoyancy or vertical stratification in this stated model.

Substitute the axisymmetric Fourier mode into the linearized equations. With hats denoting amplitudes,For a nondegenerate nonzero frequency, eliminate the azimuthal velocity and useContinuity then givesEquivalently, taking the determinant of the original four amplitude equations avoids division by zero and gives the full dispersion relationThe two squared frequencies areThe discriminant is $({\Omega }^2-c^2{k}^2{)}^2+4{\Omega }^2{c}^2{k}_x^2\ge 0$. For $c^2>0$, both squared frequencies are nonnegative. The high-frequency branch is predominantly acoustic and the low-frequency branch is an inertial wave. Zero-frequency limiting modes, including $k_z=0$, are retained by the determinant calculation.

Take $ck\gg \Omega$, retaining a finite low-frequency branch and a nonzero wavevector. In the dispersion relation, neglect ${\omega }^4$ against $c^2{k}^2{\omega }^2$ and neglect ${\Omega }^2$ against $c^2{k}^2$ in their coefficient. The inertial wave limit isChoose the signed branches $\omega =ϵ\Omega k_z/k$, $ϵ=\pm 1$. The phase velocity normal to the phase planes and the group velocity areTherefore $\mathbf{k}\cdot {\mathbf{c}}_g=ϵ\Omega (k_z/k-k_z{k}^2/k^3)=0$, soThe geometric reason is that this limiting dispersion relation is homogeneous of degree zero in $\mathbf{k}$: its frequency depends on direction, not magnitude. Euler homogeneous function theorem gives $\mathbf{k}\cdot {\nabla }_{\mathbf{k}}\omega =0$. Exact orthogonality holds for this inertial wave limit; finite compressibility introduces corrections. At a zero-frequency degeneracy the usual propagating-wave interpretation needs a limiting treatment.

For the polytropic equation of state $p=K{\rho }^{1+1/n}$, take $K>0$ and $n>0$ for the physical model. The specific enthalpy relative to vacuum isIts gradient obeys $\nabla Q={\rho }^{-1}\nabla p$, so the Euler equations becomeFor the material derivative $D_t={\partial }_t+\mathbf{u}\cdot \nabla$, $D_{t}Q=(Q/n\rho )D_t\rho$. The continuity equation thus becomesThe choice $Q=0$ at zero mass density fixes the additive convention in this equation. Although the momentum equation alone is unchanged by adding a constant to $Q$, the polytropic relation and enthalpy evolution must use the same vacuum reference.

For $\mathbf{u}=\alpha y{\mathbf{e}}_x-\beta x{\mathbf{e}}_y$, the velocity divergence is zero and $\mathbf{u}\cdot \nabla \mathbf{u}=-\alpha \beta (x{\mathbf{e}}_x+y{\mathbf{e}}_y)$. Substituting into the steady Euler equations yieldsIntegration givesThe steady enthalpy equation also requires $\mathbf{u}\cdot \nabla Q=0$. Its $xy$ coefficient is $\alpha (3{\Omega }^2+\alpha \beta -{\beta }^2)$, so, since $\alpha >0$,For a bounded patch with positive mass density and a vacuum fluid free boundary, choose $Q_0>0$ and require both quadratic coefficients to confine the fluid. With $\Omega >0$ and $\alpha ,\beta >0$, this meansThe lower bound is imposed by $\alpha >0$; the upper bound by bounded positive $Q$. At $\beta =2\Omega$, $Q$ is constant and has no nontrivial finite vacuum boundary; $\beta \le \sqrt{3}\Omega$ violates or degenerates the stipulated elliptical flow.

For a polytropic elliptical patch in a shearing sheet, the physical fluid free boundary is $p=\rho =Q=0$. It is the ellipsewhere $a_x^2=2Q_0/[\beta (2\Omega -\beta )]$ and $a_y^2=2Q_0/[\alpha (2\Omega -\beta )]$. In the interior $\rho =[Q/((n+1)K){]}^n$. The boundary is a material streamline because $\mathbf{u}\cdot \nabla (\beta x^2+\alpha y^2)=0$. Since $\alpha <\beta$, it is elongated in the azimuthal $y$ direction.

This is a pressure-supported fluid patch in the supplied tidal model; no additional planetary self-gravity potential was assumed.

At $x_n=0$ and $y_n=hn$, every radial pair displacement vanishes, so $f_x^n=0$. The particles with indices $n+l$ and $n-l$ give opposite azimuthal accelerations,which cancel pairwise. The acceleration sum has absolute convergence because ${\sum }_{l\ge 1}{l}^{-2}<\infty$. With zero velocities and $x_n=0$, both the Coriolis force and radial tidal force vanish. The equally spaced row is therefore an equilibrium of the infinite gravitating particle lattice. Its absolute gravitational potential would require an additional reference or regularization; the acceleration used here is well defined.

Let ${\mathbf{d}}_{jn}=h(j-n){\mathbf{e}}_y$ and $\Delta \mathbf{\xi }={\mathbf{\xi }}_j-{\mathbf{\xi }}_n$, with ${\mathbf{\xi }}_n=(x_n^{\prime },y_n^{\prime })$. Linearizing the Newton's law of universal gravitation usesFor an azimuthal separation, the matrix is $\mathrm{diag}(1,-2)$. After canceling the equilibrium forces, the linearized shearing sheet equations areThe transverse gravitational response tends to restore a radial displacement, while the longitudinal response reinforces an azimuthal density concentration. Both sums converge for the bounded Fourier modes considered next.

Set $\xi =hk$. For the discrete Fourier mode, the displacement difference in the $l$th neighbor pair carries $e^{il\xi }-1$, andThus the infinite system reduces toThe gravitating lattice Fourier kernel $F$ is nonnegative, even and $2\pi$-periodic, with $F(0)=0$. The uniform convergence and absolute convergence of its differentiated series givesTo prove the maximum, it is not sufficient to maximize each cosine term separately: the even terms do not reach their individual maxima at $\pi$. Instead use the convergent integral representation $l^{-2}={\int }_0^{\infty }t{e}^{-lt}dt$. Summing the geometric series for sines givesFor $0<\xi <\pi$, the denominator is positive and $\sin \xi >0$, so $F^{\prime }>0$. For $\pi <\xi <2\pi$, $F^{\prime }<0$. The interchange is justified by ${\sum }_l{\int }_0^{\infty }t{e}^{-lt}|\sin (l\xi )|dt\le {\sum }_l{l}^{-2}<\infty$. HenceTheir value iswhere $\zeta$ is the Riemann zeta function. This maximizing Fourier mode alternates the signs of neighboring particle displacements.

Put $A=(Gm/h^3)F(hk)$. The two amplitude equations have coefficient matrixIts determinant gives the dispersion relationWriting $q=A/{\Omega }^2\ge 0$, its squared growth rates areand the four growth rates are both square roots of each displayed value. DefineFor $0<q<q_{-}$, the discriminant is positive, the sum of the squared rates is $q-1<0$, and their product is $2q(3-q)>0$; both squared rates are negative. These modes have purely imaginary $s$. For $q_{-}<q<q_{+}$, the squared rates form a nonreal complex conjugate pair; their square roots include roots with positive real part. This is overstability, meaning growing oscillations. For $q\ge q_{+}$, at least the plus squared rate is positive because $q>1$, so a real growing root exists. Consequently every $q>q_{-}$ is exponentially unstable.

Take the maximizing mode $hk=\pi$, so $q=(Gm/(h^3{\Omega }^2))F(\pi )$. It is unstable wheneverThis is also the threshold for the onset of exponential growth among all lattice modes, since $0\le F\le F(\pi )$. At equality the maximizing mode has repeated imaginary roots and can have algebraic growth; it is a marginal collision of oscillatory modes, not a strict-decay assertion. The $k=0$ mode has $F=0$, an epicyclic motion sector and a zero-frequency epicyclic guiding center sector with possible secular shear drift. These neutral or algebraic behaviors are distinct from the strict exponential instability established by the requested inequality.