A numerical anisotropic model

In this section we introduce the fundamentals of the numerical method we use to model our system. We give a brief description of the mathematical basis of it and the physical idea behind it. The system is treated as a continuum, which is only adequate for a large number of stars and in well populated regions of the phase space. We consider here spherical symmetry and single-mass stars. We handle relaxation in the Fokker-Planck approximation, i.e. like a diffusive process determined by local conditions. We make also use of the hydrodynamical approximation; that is to say, only local moments of the velocity dispersion are considered, not the full orbital structure. In particular, the effect of the two-body relaxation can be modelled by a local heat flux equation with an appropriately tailored conductivity. Neither binaries nor stellar evolution are included at the presented work. As for the hypothesis concerning the BH, see section ().

For our description we use polar coordinates, $r$ $\theta$, $\phi$. The vector ${\bf v} = (v_i), i=r,\theta,\phi$ denotes the velocity in a local Cartesian coordinate system at the spatial point $r,\theta,\phi$. For succinctness, we shall employ the notation $u=v_{\rm r}$, $v=v_\theta$, $w=v_\phi$. The distribution function $f$, is a function of $r$, $t$, $u$, $v^2+w^2$ only due to spherical symmetry, and is normalised according to

(r,t) = f(r,u,v^2+w^2,t) dudvdw.

Here $\rho(r,t)$ is the mass density; if $m_{\star}$ denotes the stellar mass, we get the particle density $n=\rho/m_{\star}$. The Euler-Lagrange equations of motion corresponding to the Lagrange function

L = 12(r^2 + r^2.^2 + r^2 ^2 .^2) - (r,t)

are the following

$${\dot u} = - \frac{\partial \Phi}{\partial r} + {v^2\!+\!w^2\over r}$$

$${\dot v} = - {uv\over r} + {w^2\over r\tan\theta}$$ (15)

$${\dot w} = - {uw\over r} - {vw\over r\tan\theta}$$

And so we get a complete local Fokker-Planck equation,

ft+v_rfr+v_r fv_r+v_f v_+v_fv_= ( ft )_FP

In our model we do not solve the equation directly; we use a so-called momenta process. The momenta of the velocity distribution function $f$ are defined as follows

<i,j,k>:=^+_- v^i_r v^j_ v^k_ f(r, v_r, v_,v_,t)dv_rdv_dv_;

We define now the following moments of the velocity distribution function,

$$<0,0,0> := \rho = \int f dudvdw$$

$$<1,0,0> := {u} = \int uf dudvdw$$

$$<2,0,0> := p_{\rm r} + \rho {u}^2 = \int u^2 f dudvdw$$

$$<0,2,0> := p_\theta = \int v^2 f dudvdw$$ (16)

$$<0,0,2> := p_\phi = \int w^2 f dudvdw$$

$$<3,0,0> := F_{\rm r} + 3{u}p_{\rm r} + {u}^3 = \int u^3 f dudvdw$$

$$<1,2,0> := F_\theta + {u}p_\theta = \int uv^2 f dudvdw$$

$$<1,0,2> := F_\phi + {u}p_\phi = \int uw^2 f dudvdw,$$

where $\rho$ is the density of stars, $u$ is the bulk velocity, $v_{\rm r}$ and $v_{\rm t}$ are the radial and tangential flux velocities, $p_{\rm r}$ and $p_{\rm t}$ are the radial and tangential pressures, $F_{\rm r}$ is the radial and $F_{\rm t}$ the tangential kinetic energy flux LS91. Note that the definitions of $p_i$ and $F_i$ are such that they are proportional to the random motion of the stars. Due to spherical symmetry, we have $p_{\theta} = p_{\phi }=: p_{\rm t}$ and $F_{\theta} = F_{\phi} =: F_{\rm t}/2$. By $p_{\rm r} = \rho\sigma_{\rm r}^2$ and $p{_{\rm t}} = \rho\sigma_{\rm t}^2$ the random velocity dispersions are given, which are closely related to observable properties in stellar clusters.

$F = (F_{\rm r} + F_{\rm t})/2$ is a radial flux of random kinetic energy. In the notion of gas dynamics it is just an energy flux. Whereas for the $\theta-$ and $\phi-$ components in the set of Eqs. (2.20) are equal in spherical symmetry, for the $r$ and $t$- quantities this is not true. In stellar clusters the relaxation time is larger than the dynamical time and so any possible difference between $p_{\rm r}$ and $p_{\rm t}$ may survive many dynamical times. We shall denote such differences anisotropy. Let us define the following velocities of energy transport:

$$v_{\rm r} = {F_{\rm r} \over 3 p_{\rm r}} + u,$$ (17)

$$v_{\rm t} = {F_{\rm t} \over 2 p_{\rm t}} + u.$$

In case of weak isotropy ($p_{\rm r}$=$p_{\rm t}$) $2F_{\rm r}$ = $3F_{\rm t}$, and thus $v_{\rm r}$ = $v_{\rm t}$, i.e. the (radial) transport velocities of radial and tangential random kinetic energy are equal.

The Fokker-Planck equation (2.18) is multiplicated with various powers of the velocity components $u$, $v$, $w$. We get so up to second order a set of moment equations: A mass equation, a continuity equation, an Euler equation (force) and radial and tangential energy equations. The system of equations is closed by a phenomenological heat flux equation for the flux of radial and tangential RMS ( root mean square) kinetic energy, both in radial direction. The concept is physically similar to that of LBE80. The set of equations is

$$\begin{eqnarray*} \frac{\partial{\rho}}{\partial t} + \frac{1}{r^2}\frac{\partial} {\partial r} (r^2u\,\rho)= 0 \end{eqnarray*}$$

$$\begin{eqnarray*} \frac{\partial u}{\partial t}+u\,\frac{\partial u}{\partial r}... ...r}}{\partial r} + 2\,\frac{p_{\rm r} - p_{\rm t}}{\rho\, r} = 0 \end{eqnarray*}$$

$${ \frac{\partial{p_{\rm r}}}{\partial {t}} + \frac{1}{r^2} \frac{... ...l u}{\partial r} + \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 F_{\rm r}){}}$$

$${} - \frac{2F_{\rm t}}{r} = -\frac{4}{5} \frac{(2p_{\rm r}-p_{\rm t})} {\lambda_A t_{\rm relax}}$$

$$\begin{eqnarray*} \lefteqn{ \frac{\partial{p_{\rm t}}}{\partial {t}} + \frac{1}{... ...2}{5} \frac{(2p_{\rm r}-p_{\rm t})}{\lambda_A t_{\rm relax}}, \end{eqnarray*}$$

where $\lambda_A$ is a numerical constant related to the time-scale of collisional anisotropy decay. The value chosen for it has been discussed in comparison with direct simulations performed with the $N$-body code GS94. The authors find that $\lambda_A=0.1$ is the physically realistic value inside the half-mass radius for all cases of $N$, provided that close encounters and binary activity do not carry out an important role in the system, what is, on the other hand, inherent to systems with a big number of particles, as this is.

With the definition of the mass $M_{\rm r}$ contained in a sphere of radius $r$

M_rr = 4 r^2 ,

the set of Eqs.(2.22) is equivalent to gas-dynamical equations coupled with the equation of Poisson. To close it we need an independent relation, for moment equations of order $n$ contain moments of order $n\,+\,1$. For this intent we use the heat conduction closure, a phenomenological approach obtained in an analogous way to gas dynamics. It was used for the first time by [Lynden-Bell and Eggleton, 1980] but restricted to isotropy. In this approximation one assumes that heat transport is proportional to the temperature gradient,

F = -Tr = - ^2r

That is the reason why such models are usually also called conducting gas sphere models.

It has been argued that for the classical approach $\Lambda\propto\bar{\lambda}^2/\tau$, one has to choose the Jeans' length $\lambda_J^2 = \sigma ^2/(4\pi G\rho)$ and the standard Chandrasekhar local relaxation time $t_{\rm relax}\propto \sigma ^3/\rho$ LBE80, where $\bar{\lambda}$ is the mean free path and $\tau$ the collisional time. In this context we obtain a conductivity $\Lambda\propto \rho/ \sigma$. We shall consider this as a working hypothesis. For the anisotropic model we use a mean velocity dispersion $\sigma^2 = (\sr^2 + 2\st^2)/3$ for the temperature gradient and assume $v_{\rm r} = v_{\rm t}$ BS86. Forasmuch as, the equations we need to close our model are

$${v_{\rm r} - u + \frac{\lambda}{4\pi \,G\rho\,t_{\rm relax}} \frac{\partial \sigma^2}{\partial r} = 0}$$ (19)

$$v_{\rm r} = v_{\rm t}.$$