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

Several decades ago, Hachisu and collaborators (Hachisu et al. 1978) followed by Lynden-Bell and Eggleton (1980) proposed using transport processes in a self-gravitating, conducting gas sphere to simulate two-body stellar relaxation. This proved remarkably effective for modelling dense stellar clusters, and this site introduces the fundamentals alongside an implementation of this approximation.

Subsequent work by various authors (Bettwieser 1983, Bettwieser and Sugimoto 1984, Bettwieser and Spurzem 1986, Heggie 1984, Heggie and Ramamani 1989, and Louis and Spurzem 1991) incorporated anisotropy, while Giersz and Heggie (1994), 1994b and Spurzem and Takahashi (1995) introduced multi-mass distributions and refined conductivity formulations for improved accuracy.

Pau Amaro-Seoane, working with Marc Freitag and Rainer Spurzem, extended the model to include a central accreting massive black hole, studying its growth through stellar accretion at the tidal radius and the feedback effects on core collapse and post-collapse cluster evolution (Amaro-Seoane et al. 2004).

This approach treats spherically symmetric systems as continua, modelling relaxation as a diffusive process in phase space via the Fokker-Planck (FP) equation. The local approximation simplifies the FP equation by neglecting positional diffusion – justified as encounters occur in volumes much smaller than the system’s overall dimensions. Energy transfer follows a local heat flux equation with carefully tailored conductivity.

The model’s foundation is the FP equation describing the time evolution of the probability density function. In spherical coordinates, the Boltzmann equation becomes:

$\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}$

Here the right-hand side represents collisional terms in the FP approximation. Symmetry permits definition of a tangential velocity $v_t^2=v_{\rm{\phi}}^2+v_{\rm{\theta}}^2$ , reducing the system to two velocity components ( $v_t$ and $v_r$ ). The “centralised” moments derive from multiplying the velocity distribution $f$ by powers of $v_t$ and $(v_r-\bar{v}_r)$ and integrating over velocity space.

“Centralised” indicates moments defined relative to mean velocity components $\bar{v}_r=\langle v_r \rangle = u$ and $\bar{v}_t=\langle v_t \rangle = 0$ (from spherical symmetry). Moment order equals $n+m$ where $n$ and $m$ are velocity powers in $\int (v_r-\bar{v}_r)^nv_t^m\,f\,\, \mathrm{d}^3v$ . These correspond to stellar density ( $\rho$ ), bulk velocity ( $u$ ), radial/tangential pressures ( $p_r$ , $p_t$ ), and kinetic energy fluxes ( $F_r$ , $F_t$ ).

Multiplying the equation by velocity powers and integrating yields differential moment equations. To second order, these comprise:

Continuity: $\frac{\partial \rho}{\partial t} + \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 u \rho) = 0$
Euler (force): $\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$
Energy equations:
- Radial: $\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}$
- Tangential: $\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 governing collisional anisotropy decay timescales, crucial for relaxation effects. While theoretically $\lambda_A=1$ for globular clusters in higher-moment models, simulations suggest $\lambda_A=0.1$ provides physically realistic values within the half-mass radius across particle numbers.

The right-hand terms in the energy equations represent collisional effects: the first term describes relaxation from uncorrelated two-body encounters (derivable from FP theory), while the “$latex{\rm bin3}$” terms account for star-binary encounters. Sufficiently energetic three-body interactions can reverse core collapse.

Pressures relate to velocity dispersions: $p_r= \rho \, \sigma_r^2$ and $p_t= \rho \, \sigma_t^2$ , connecting to observable cluster properties. The average dispersion is $\sigma^2 = (\sigma_r^2 + 2\sigma_t^2)/3$ (factor 2 accounting for two tangential dimensions), with radial energy flux $F=(F_r + F_t)/2$ evident from combining moment equations.

Energy transport velocities are defined by:

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

Under weak isotropy ( $p_r = p_t$ everywhere), $F_r = \frac{3}{2}F_t$ ensures $v_r = v_t$ , equating transport velocities for radial and tangential energy.

Closure requires three additional equations: mass relation and two heat flux equations. The mass relation defines $M_r$ within radius $r$ :

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

where $\rho_M = M \cdot \rho$ is mass density for stellar mass $M$ .

This yields gas dynamical equations coupled with Poisson’s equation… where the real fun begins!

Since nth-order moment equations contain (n+1)th-order moments, we require closure relations… correct?

Enter the heat conduction closure – a phenomenological approach analogous to gas dynamics.

Consider a star… fundamentally an atomic ensemble bearing intriguing similarities to stellar clusters. Beyond Virial theorem analogies, both systems share heat transport, energy generation and core-halo evolution characteristics.

We therefore assume heat transport proportional to temperature gradient:

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

This gas/liquid heat flux equation explains why we term this the “conducting gas sphere model”.

Caution: this assumes small mean free paths – questionable for stellar systems. However, validation exists through agreement with N-body and FP models (Giersz and Heggie 1994, 1994b; Spurzem and Takahashi 1995).

Classically, $\Lambda \propto \rho\bar{\lambda}^2/\tau$ , with mean free path $\bar{\lambda}$ and collisional time $\tau$ . Choosing Jeans length $\lambda_J^2 = \sigma^2/(4 \pi G \rho)$ for $\bar{\lambda}^2$ and Chandrasekhar’s relaxation time $t_{\rm{rx}} \propto \sigma^3/\rho$ for $\tau$ yields conductivity $\Lambda \propto \rho/\sigma$ , specifically:

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

where dimensionless constant $C \approx 1$ . Transforming via energy transport velocities gives anisotropic closure relations:

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

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

Note $\lambda$ is a free parameter requiring calibration against N-body/FP models. In isotropic cases it merely scales solutions, but with anisotropy it balances two processes: anisotropy decay versus heat flow between regions. Larger $\lambda$ accelerates heat flow, preserving anisotropy – thus stronger $\lambda$ produces greater anisotropy.