Modelling galactic nuclei with self-gravitating, conducting gas spheres

Getting the code

git clone git@github.com:pau-amaro-seoane/gaseous-model.git

License

* Copyright (c) 2016, Pau Amaro Seoane (pau.amaro.seoane@gmail.com)
* and Rainer Spurzem (spurzem@nao.cas.cn)
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED “AS IS” AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

Background

Some decades ago, Hachisu and collaborators (Hachisu et al 1978) and a bit later Lynden-Bell and Eggleton (1980) proposed to use transport process in a self-gravitating, conducting gas sphere as a way to mimic two-body stellar relaxation. This turned out to be an excellent way of modelling dense stellar clusters, and in this site we will briefly introduce the very basics of the subject, and a code that uses the approximation.

A number of authors(Bettwieser 1983, Bettwieser and Sugimoto 1984, Bettwieser and Spurzem 1986, Heggie 1984, Heggie and Ramamani 1989 and Louis and Spurzem 1991 implemented anisotropy and Giersz and Heggie, 1994, 1994b and Spurzem and Takahashi 1995 added a multi-mass distribution and improved the detailed form of the conductivities to have better accuracy.

Pau Amaro-Seoane, with Marc Freitag and Rainer Spurzem added to the model an accreting, central massive black hole to study its growth due to star accretion at its tidal radius. They also addressed the feedback of this process on to the core collapse as well as the post-collapse evolution of the surrounding stellar cluster (Amaro-Seoane et al 2004).

In this approach, we emulate spherically symmetric systems as a continuum; relaxation is treated as a diffusive process in phase space using the Fokker-Planck (FP) equation. We employ the local approximation to simplify the FP equation by neglecting the diffusion in position. The idea behind this is that an encounter takes place in a volume that is much smaller than the dimensions of the whole system. We model energy transfer by a local heat flux equation with an appropriately tailored conductivity.

The basis of the equations of the model is the FP equation which describes the time evolution of the probability density function. Using spherical polar coordinates, the Boltzmann equation takes the form:

$\frac{\partial f}{\partial t} + v_{\rm{r}} \frac{\partial f}{\partial r} + \dot{v}_{\rm{r}} \frac{\partial f}{\partial v_{\rm{r}}} + \dot{v}_{\rm{\phi}} \frac{\partial f}{\partial v_{\rm{\phi}}} + \dot{v}_{\rm{\theta}} \frac{\partial f}{\partial v_{\rm{\theta}}} = \left(\frac{\delta f}{\delta t}\right)_{FP}$

In the last equation the right-hand side denotes that collisions are given in terms of the FP approximation. Due to symmetry, we can define a tangential velocity $v_t^2=v_{\rm{\phi}}^2+v_{\rm{\theta}}^2$ , so that we have two velocities $v_t$ and $v_r$ to describe the system. The “centralized” moments are defined by multiplying the velocity distribution $f$ with powers of $v_t$ and $(v_r-\bar{v}_r)$ and integrating over velocity space.

The term “centralized” means that the moments are defined with respect to the mean velocity components $\bar{v}_r=\langle v_r \rangle = u$ and $\bar{v}_t=\langle v_t \rangle = 0$ , because we assume spherical symmetry. The order of a moment is defined by $n+m$ where $n$ and $m$ are the powers of velocities in the definition of moments, i.e. $\int (v_r-\bar{v}_r)^nv_t^m\,f\,\, \mathrm{d}^3v$ . The moments defined this way correspond to the density of stars, $\rho$ , the bulk velocity, $u$ , the radial and tangential pressures, $p_r$ and $p_t$ , and the radial and tangential kinetic energy fluxes, $F_r$ and $F_t$ . In order to obtain the set of differential moment equations, we multiply the last equation with powers of $v_t$ and $(v_r-\bar{v}_r)$ and integrate it in velocity space. After some recasting the integrals can be substituted by the moments. Up to second order the moment equations are the continuity equation, the Euler equation (force) and radial and tangential energy equations:

$\frac{\partial \rho}{\partial t} + \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 u \rho) = 0$
$\frac{\partial u}{\partial t} + u \,\frac{\partial u}{\partial r} + \frac{G M_r}{r^2} + \frac{1}{\rho} \frac{\partial p_r}{\partial r} + 2 \frac{p_r - p_t}{\rho r} = 0$
$\frac{\partial p_r}{\partial t} + \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 u p_r) + 2 p_r \,\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 F_r) - \frac{2 F_t}{r} = -\frac{3}{5} \frac{p_r-p_t}{\lambda_A t_{\mathrm{rx}}} + \left(\frac{\delta p_r}{\delta t}\right)_{\rm{bin}3} \frac{\partial p_t}{\partial t} + \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 u p_t ) + \frac{2 p_r u}{r} + \frac{1}{2 r^2}\frac{\partial}{\partial r}(r^2 F_t) + \frac{F_t}{r} = \frac{3}{10} \frac{p_r-p_t}{\lambda_A t_{\mathrm{rx}}} + \left(\frac{\delta p_t}{\delta t}\right)_{\rm{bin}3}$

Here $\lambda_A$ is a numerical constant related to the time-scale of collisional anisotropy decay, necessary to describe the relaxation effects on cluster evolution. It should become unity when describing GCs using higher moment models. However, this can only be confirmed by simulations; $\lambda_A=0.1$ is a physically realistic value within the half-mass radius for all numbers of particles.

The two terms on the right-hand sides of the equations for radial and tangential energy equations are the collisional terms. The fist term accounts for relaxation from uncorrelated two-body encounters and can be derived from the FP equation. The second term, which is marked with “$latex{\rm bin3}”, refers to star-binary encounters. Close 3-body or star-binary encounters generate kinetic energy. If the energy generation is high enough, this mechanism can reverse core collapse.

The radial and tangential pressure, $p_r$ and $p_t$ are related to the random velocity dispersions; $p_r= \rho \, \sigma_r^2$ and $p_t= \rho \, \sigma_t^2$ . They are linked to observable quantities in stellar clusters such as the radial velocity dispersion. The average velocity dispersion is $\sigma^2 = (\sigma_r^2 + 2\sigma_t^2)/3$ , where the factor 2 comes from the fact that there are two tangential directions. The radial energy flux of random kinetic energy is $F=(F_r + F_t)/2$ . We can see this by adding the two-moment equations for radial and tangential pressure to obtain the gas-dynamical equation for the energy density. The velocities for energy transport are defined by

$v_r = \frac{F_r}{3p_r} + u$
$v_t = \frac{F_t}{2p_t} + u$

In the case of weak isotropy, $p_r = p_t$ everywhere and hence $F_r = \frac{3}{2}F_t$ , so that $v_r = v_t$ . Therefore, the transport velocities for radial and tangential random kinetic energy are equal.

In order to close the set of moment equations, we need to set up three more equations, a mass relation and two equations accounting for heat flux. The mass relation defines the mass $M_r$ contained in a sphere of radius $r$ ,

$\frac{\partial M_r}{\partial r} = 4 \pi r^2 \rho_M$

where $\rho_M = M \cdot \rho$ is the mass density and $M$ the mass of the stellar component.

We thus obtain a set of gas dynamical equations coupled with the Poisson equation… and this is the fun part!

Since the moment equations of order $n$ obtained from the Boltzmann equation contain moments of the order $n+1$ , we need closure relations connecting the moments of order $n+1$ with lower-order moments…. right?

Ok… we can achieve it with the heat conduction closure.

This a phenomenological approach obtained in an analogous way to gas dynamics.

Think of a star… it’s just a bunch of atoms… and it somehow looks similar to a stellar cluster with a looooot of stars, right? The comparison goes beyond a simple analogy on the level of the Virial theorem… the two systems actually have similarities in heat transport, energy generation and core-halo evolution.

So… we can simply assume that heat transport is proportional to the temperature gradient:

$F = - \kappa \frac{\partial T}{\partial r} = - \Lambda \frac{\partial \sigma^2}{\partial r}$

This equation describes the heat flux in gases and liquids and… guess! : This is the reason why we usually refer to these model as the “conducting gas sphere model”.

But watch out… the use of the heat flux is based on the assumption of small mean free paths for the particles, which is certainly questionable for stellar dynamical systems.

However, we are safe… this has been tested with other models, like direct N-body or Fokker-Planck with good results (Giersz and Heggie, 1994, 1994b and Spurzem and Takahashi 1995), and there is a pretty good agreement.

In the classical approach $\Lambda \propto \rho\bar{\lambda}^2/\tau$ , where $\bar{\lambda}$ is the mean free path and $\tau$ the collisional time. Choosing the Jeans length $\lambda_J^2 = \sigma^2/(4 \pi G \rho)$ for $\bar{\lambda}^2$ and the standard Chandrasekhar local relaxation time $t_{\rm{rx}} \propto \sigma^3/\rho$ for $\tau$ , we obtain the conductivity $\Lambda \propto \rho/\sigma$ . More precisely, the conductivity takes the form

$\Lambda = \frac{3 C G m \rho N}{\sigma}$

where C is a dimensionless numerical constant of order unity. By means of the velocities of energy transport the heat flux equation can be recast to find the two closure relations in the anisotropic case

$v_r-u+ \frac{\lambda}{4 \pi G \rho t_{\rm{rx}}} \frac{\partial \sigma^2}{\partial r} = 0 \qquad v_r = v_t$ ,

where $\lambda = \frac{27 \sqrt{\pi}}{10} C$ .

We should emphasized that $\lambda$ is a free parameter that has to be determined by comparison with other models such as N-body or FP models.

In the isotropic limit, $\lambda$ is just a scaling scaling factor, but when taking into account anisotropy, it prescribes the relative speed of two processes: the decay of anisotropy and the heat flow between warm and cold regions. With increasing $\lambda$ , heat flows faster, so there is less time for gravitational encounters to “destroy” anisotropy. A larger $\lambda$ thus results in stronger anisotropy.