Accretion Discs

Suppose we have a massive object enshrouded in a gas cloud. Due to gravity, some of the gas will fall straight onto the massive object, but most of the gas particles would just move on closed Keplerian orbits. If the initial gas cloud has some angular momentum, and if some dissipation is present, then it will settle into the shape of a disk. This is called the accretion disk. Usually, it is much easier for gas to rid itself of energy than angular momentum, so most orbits in an accretion disks would be circular.

Another scenario where accretion disks occur is in a system of mutually orbiting binary stars, when one stars exudes mass (as either stellar wind or Roche lobe overflow), and the companion captures some of it.

Radial Structure

We describe an accretion disk in cylindrical coordinates. We assume that the z axis is perpendicular to the plane of the disk, that variables will only depend on the distance from the z axis (denoted by r) and the time t.

Conservation of mass takes the form

$\frac{\partial \Sigma}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma v_r \right) = 0$

Where $\Sigma = \int \rho dz$ is the surface density of the disk. The other governing equation is the conservation of angular momentum. If the flow were completely ideal, then the gas particles in the disk will continue to revolve indefinitely around the massive object. In order for the matter to eventually fall onto the surface of the massive object, some viscosity has to be involved. This can be a bit tricky, as most representations of the Navier Stokes equations assume that the viscosity coefficient $\mu$ is constant. This can easily be remedied, by recalling that strain rate is $\varepsilon = \nabla \mathbf{v}$, stress is linear in the strain $\sigma = \mu \epsilon = \mu \nabla \mathbf{v}$ and net stress is the difference between stresses acting on all sides of a fluid element $\nabla \left( \mu \nabla \mathbf{v} \right)$. The last expression is the viscosity term in the Navier Stokes equation, only the viscosity coefficient is usually assumed to be constant, so it is taken out of the brackets. As we shall later see, in our case the viscosity coefficient varies, so we can not make the same simplification. With all that being said, the conservation of angular momentum takes the form

$\frac{\partial}{\partial t} \left( \Sigma r^2 \Omega \right) + \frac{1}{r} \frac{\partial}{\partial r} \left( \Sigma r^3 \Omega v_r \right) = \frac{\partial}{\partial r} \left( \nu \Sigma r^2 \frac{\partial \Omega}{\partial r} \right)$

Where $\nu = \frac{\mu}{\rho}$ is the kinematic viscosity and $\Omega = \frac{v_{\theta}}{r}$ is the angular velocity. Usually, the angular velocity would be Keplerian $\Omega_K = \sqrt{\frac{G M}{r^3}}$, but it does not complicate things if we consider the general case. We can verify that for rigid body rotation (constant angular velocity) the dissipation term vanishes. Solving the two equations (conservation of mass and conservation of angular momentum), assuming the angular velocity is a known function of the radius and independent of time, yields

$v_r = \frac{1}{\Sigma} \frac{\partial}{\partial r} \left( \nu \Sigma r^2 \frac{\partial \Omega}{\partial r} \right) / \frac{\partial }{\partial r} \left( r^2 \Omega \right)$

$\frac{\partial \Sigma}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left[ r \frac{\partial}{\partial r} \left( \nu \Sigma r^2 \frac{\partial \Omega}{\partial r} \right) / \frac{\partial }{\partial r} \left( r^2 \Omega \right) \right] = 0$

The last equation is a non - linear diffusion equation for the mass.

We note that the viscosity mechanism is somewhat of a mystery. A simple calculation can show that simple atomic viscosity does not suffice. Typical densities and temperatures for accretion disks are $n \approx 10^{16} cm^{-3}$ and $T \approx 10^4 K$. A typical radius is of the order of magnitude of $R \approx 10^{15} cm$ (like Herbig - Haro object number 30). The diffusion time scale is

$\tau \approx \frac{R^2}{\nu} \approx \frac{R^2}{n^{-1/3} \sqrt{k T /m}} = \frac{R^2 n^{1/3} }{\sqrt{k T/m}}$

Substituting numbers yields a time scale larger than the age of the universe by a factor of about a trillion.

Currently, the leading candidate for viscosity mechanism is the magneto - rotational instability. However, this could not be the end of the story, because accretion also occurs in disks that are not ionised (like proto - planetary disks). Another promising candidates is dissipation by means of spiral shock waves.

Analytic Solutions to the Diffusion Equation

If we assume a Keplerian angular velocity and that kinematic viscosity is constant, then the non linear diffusion equation admits an analytic solution.

$\frac{\partial \Sigma}{\partial t} = \frac{3 \nu}{\sqrt{r}} \frac{\partial}{\partial r} \left( \sqrt{r} \frac{\partial \Sigma}{\partial r} \right)$

The solution can be separated into a product of a spatial and temporal component $\Sigma \left( t, r\right) = R \left(r \right) T \left(t\right)$. The temporal component is simply an exponent $T \left(t \right) \propto \exp \left( \omega t \right)$. The spatial component is a Bessel function

$R \left( r \right) = r^{1/4} \left(C_1 J_{1/4} \left(\frac{i}{\sqrt{3}} \frac{r}{\sqrt{\nu/\omega}}\right) + Y_{1/4} \left( \frac{i}{\sqrt{3}} \frac{r}{\sqrt{\nu/\omega}} \right) \right)$

Steady State Solution

If mass is supplied to the disk at a constant rate for long enough, the disk will reach steady state. The accretion rate, however, may not be equal to the mass supplied because some of the matter must become unbound and carry the excess energy and angular momentum with it.

Dropping the time derivatives from the non linear diffusion yields

$r \frac{\partial}{\partial r} \left( \nu \Sigma r^2 \frac{\partial \Omega}{\partial r} \right) / \frac{\partial}{\partial r} \left( r^2 \Omega \right) = C_1$

Where C is constant. Next we will assume a Keplerian angular velocity and solve for $\nu \Sigma$.

$\nu \Sigma = C_1 - \frac{C_2}{\sqrt{r}}$

Since the massive object is a sink for the mass, we require that the mass density vanish at its surface. This provide one boundary condition

$\nu \Sigma = C_1 \left[ 1 - \left( \frac{R_*}{r} \right)^{1/2} \right]$

The other coefficient is determined from the condition that the mass accretion rate is constant throughout

$\dot{M} = -2 \pi r \Sigma v_r = -2 \pi r \frac{\partial }{\partial r} \left[ \nu \Sigma r^2 \frac{\partial \Omega}{\partial r} \right] / \frac{\partial }{\partial r} \left( r^2 \Omega \right) = 3 \pi C_1$

Hence

$\nu \Sigma = \frac{\dot{M}}{3 \pi} \left[ 1 - \left(\frac{R_*}{r}\right)^{1/2} \right]$

Energy Dissipation Rate

So far we ignored the pressure. For that assumption to hold, all the heat from the friction must be dissipated. The work the friction does, per unit time per unit area is equal to angular velocity times the net moment exerted by the frictional force

$l = \frac{\Omega \cdot \frac{d}{dr} \left(2 \pi r^3 \nu \Sigma \frac{d \Omega}{d r} \right) }{4 \pi r}$

The reason the denominator has $4\pi$ instead of $2\pi$ is that the disk has two sides. Substituting the results from the previous section yields

$l = \frac{G M \dot{M} }{8 \pi r^3}$

Integrating over the entire disk yields

$L = 2 \int_{r_s}^{\infty} 2 \pi r dr \cdot l = \frac{1}{2} \frac{G M \dot{M} }{r_s}$

This means that at least half the potential energy of the particle is lost on its way to the surface of the massive object.

Vertical Structure

The thickness of the disk is directly related to the temperature. The hydrostatic equation in the vertical direction is

$\frac{1}{\rho} \frac{\partial p}{\partial z} = \frac{\partial}{\partial z} \left( \frac{G M}{\sqrt{r^2 + z^2}} \right)$

Assuming a thin disk $r \gg z$

$\frac{1}{\rho} \frac{\partial p}{\partial z} \approx - \frac{G M}{r^3} z$

If $H \approx z$ is the height of the disk then

$\frac{p}{\rho H} \approx \frac{c_s^2}{H} \approx \frac{G M H}{r^3}$

$H = r \frac{c_s}{v_k}$

So the height at each point is equal to the radius times the ratio between the speed of sound and the Keplerian velocity. As long as the speed of sound is much smaller than the Keplerian velocity, the disk can be considered thin.

The vertical optical depth of the disk is

$\tau = \kappa \rho H = \kappa \Sigma$

where $\kappa$ is the opacity. If the disk is opaque $\tau \gg 1$ then the temperature can be estimated using the Stefan Boltzmann law

$\sigma T^4 = l = \frac{G M \dot{M}}{r^3}$

$T = \left( \frac{G M \dot{M}}{\sigma r^3} \right)^{1/4}$

The spectrum is given by integrating the black body spectral radiance over the entire area. We will perform this calculation for a general power law dependence of the temperature on distance $T = q r^{-\beta}$

$F_{f} = \int_{r_{\min}}^{r^{\max}} B_f 2 \pi r dr = \int_{r_{\min}}^{r_{\max}} \frac{2 h f^3}{c^2} \frac{2 \pi r}{\exp \left( \frac{h f}{k T} \right) - 1} dr =$

$= \frac{2 \pi g f^3}{c^2} \left( \frac{k q}{h f} \right)^{2/\beta} \int_{\xi_{min}}^{\xi_{\max}} \frac{d\left(\xi^{2/\beta}\right)}{\exp \left( \xi \right) - 1}$

where $\xi = \frac{h f}{k T}$. At extremely low frequencies (photon at a lower energy than the thermal energy at the far end of the disk) $\xi_{\min}, \xi_{\max} \rightarrow 0$, the integral can be approximated as

$F_f \approx \frac{2 \pi h f^3}{c^2} \left( \frac{k q}{h f} \right)^{2/\beta} \left( \frac{h f}{k T_{\min}} \right)^{2/\beta-1} \frac{\frac{h f}{k T_{\min}} - \frac{h f}{k T_{\max}}}{\exp{\frac{h f}{k T_{\min}}}-1} \propto f^2$

An example for black body spectrum of an accretion disk. The code used to generated this image can be found here.

Hence, at very low frequencies, the spectrum increases as the frequency squared, independent of the temperature profile. At higher frequencies the numerical integral is approximately constant, so $F_f \propto f^{3-2/\beta}$. For $\beta = \frac{3}{4}$, $F_f \propto f^{1/3}$. At very high frequencies the dimensionless integral cannot be ignored, but it can be approximated. It will add some polynomial dependence in the frequency, but, more importantly, it will introduce a more dominant exponential dependence

$F_f \propto \exp \left(- \frac{h f}{k T_{in}} \right)$, where $T_{in}$ is the highest temperature in the disk, closest to the massive object.