1. Introduction
It is well known that spherically symmetric, self-gravitating collisionless equilibrium systems in which a significant fraction of the kinetic energy is stored in low-angular-momentum orbits are dynamically unstable. The associated instability is known as radial orbit instability (ROI) – see e.g. Fridman & Polyachenko (Reference Fridman and Polyachenko1984), Binney & Tremaine (Reference Binney and Tremaine2008) and Bertin (Reference Bertin2014); for a recent discussion see also Polyachenko & Shukhman (Reference Polyachenko and Shukhman2015) and references therein. Usually, the amount of radial anisotropy in a spherical system is quantified by introducing the Fridman–Polyachenko–Shukhman parameter
$$\begin{eqnarray}{\it\xi}\equiv \frac{2K_{r}}{K_{t}},\end{eqnarray}$$
where the radial and tangential kinetic energies are given respectively by
$$\begin{eqnarray}K_{r}=2{\rm\pi}\int \!{\it\rho}(r){\it\sigma}_{r}^{2}(r)r^{2}\,\text{d}r,\quad K_{t}=2{\rm\pi}\int \!{\it\rho}(r){\it\sigma}_{t}^{2}(r)r^{2}\,\text{d}r.\end{eqnarray}$$
In the expressions above,
${\it\rho}$
is the system density, and
${\it\sigma}_{r}^{2}$
and
${\it\sigma}_{t}^{2}$
are the radial and tangential phase-space averaged square velocity components (see Polyachenko & Shukhman Reference Polyachenko and Shukhman1981; Polyachenko Reference Polyachenko1992). In particular, for isotropic systems,
${\it\xi}=1$
. In the case of Newtonian force, numerical simulations show that the ROI typically occurs for
${\it\xi}\gtrsim 1.7$
, but this should be considered a fiducial value, as it is also well known that the critical value of
${\it\xi}$
above which the given system is unstable depends on the specific phase-space structure of the equilibrium configuration under study (see e.g. Merritt & Aguilar Reference Merritt and Aguilar1985; Bertin & Stiavelli Reference Bertin and Stiavelli1989; Saha Reference Saha1991; Bertin et al.
Reference Bertin, Pegoraro, Rubini and Vesperini1994; Meza & Zamorano Reference Meza and Zamorano1997; Nipoti, Londrillo & Ciotti Reference Nipoti, Londrillo and Ciotti2002; Barnes, Lanzel & Williams Reference Barnes, Lanzel and Williams2009).
Owing to its relevant astrophysical consequences, and its intrinsic physical interest for the understanding of collisionless systems, the ROI has been extensively studied in Newtonian gravity (see e.g. Palmer & Papaloizou Reference Palmer and Papaloizou1987; Gajda, Lokas & Wojtak Reference Gajda, Lokas and Wojtak2015, and references therein). For example, Nipoti et al. (Reference Nipoti, Londrillo and Ciotti2002) investigated the implications of the ROI for the Fundamental Plane of elliptical galaxies, excluding the possibility that the so-called tilt of the Fundamental Plane is just due to a systematic increase of radial orbits with galaxy mass (see also Ciotti Reference Ciotti, Nicolaci da Costa and Renzini1997; Ciotti & Lanzoni Reference Ciotti and Lanzoni1997). In a different context (Ciotti & Pellegrini Reference Ciotti and Pellegrini2004), it has been shown that the constraints on the maximum amount of radial anisotropy that can be sustained by a stable stellar system can be used to dismiss some models of mass distribution in elliptical galaxies that have been obtained by using X-ray data under the assumption of hydrostatic equilibrium of the hot interstellar medium. The ROI also plays a role in the course of violent relaxation (Lynden-Bell Reference Lynden-Bell1967), i.e. the collapse and virialization of cold self-gravitating systems (see e.g. van Albada Reference van Albada1982; Londrillo, Messina & Stiavelli Reference Londrillo, Messina and Stiavelli1991; Nipoti, Londrillo & Ciotti Reference Nipoti, Londrillo and Ciotti2006), because such collapses are dominated in their last stages by significant radial motions, and so ROI contributes to establish the dynamical and structural properties of the quasi-relaxed final states (e.g. Merritt & Aguilar Reference Merritt and Aguilar1985; Trenti, Bertin & van Albada Reference Trenti, Bertin and van Albada2005; Sylos-Labini, Benhaiem & Joyce Reference Sylos-Labini, Benhaiem and Joyce2015).
From a more general point of view, it is now known that the ROI is not a specific property of Newtonian gravity. For example, Nipoti, Ciotti & Londrillo (Reference Nipoti, Ciotti and Londrillo2011) have studied the ROI in modified Newtonian dynamics (MOND; Milgrom Reference Milgrom1983; Bekenstein & Milgrom Reference Bekenstein and Milgrom1984), comparing anisotropic MOND systems with their equivalent Newtonian systems (ENSs) with dark matter (i.e. systems in which the phase-space properties and total gravitational potential are the same as in the corresponding MOND systems). Numerical simulations showed that MOND systems are more prone to develop ROI than their ENSs. At the same time, MOND systems are able to support a larger amount of kinetic energy stored in radial orbits than one-component Newtonian systems with the same density distribution.
Some natural questions then arise. What is the dependence of the ROI on the specific force law? How much does the long- or short-range nature of the interparticle force influence the ROI? Are there any aspects of the ROI that are peculiar to the Newtonian force (which stands out among power-law forces for several unique mathematical properties)? In fact, though several features of the dynamical evolution of collisionless systems have been found to be quite independent of the specific force law, some other aspects show a clear dependence on it. For example, in the low-acceleration regime, collisionless relaxation is much less efficient in MOND than in Newtonian gravity (Ciotti, Nipoti & Londrillo Reference Ciotti, Nipoti and Londrillo2007; Nipoti, Londrillo & Ciotti Reference Nipoti, Londrillo and Ciotti2007a
,Reference Nipoti, Londrillo and Ciotti
b
), and a similar result is found for power-law forces with
${\it\alpha}<2$
(Di Cintio & Ciotti Reference Di Cintio and Ciotti2011; however, note that, in MOND, collisional relaxation seems to be faster than in Newtonian gravity; see e.g. Ciotti & Binney Reference Ciotti and Binney2004; Nipoti et al.
Reference Nipoti, Ciotti, Binney and Londrillo2008).
In particular, Di Cintio, Ciotti & Nipoti (Reference Di Cintio, Ciotti and Nipoti2013) (see also Di Cintio & Ciotti Reference Di Cintio and Ciotti2011) studied by means of direct
$N$
-body simulations the dissipationless collapses of cold and spherical systems of particles interacting via additive interparticle
$r^{-{\it\alpha}}$
forces. The main results of these works are that, almost independently of
${\it\alpha}$
, the considerably radially anisotropic final states of cold collapses have surface density profiles following the Sérsic (Reference Sérsic1968) law, and differential energy distributions well described by an exponential over a large range of energies. At fixed initial virial ratio, however, a non-monotonic trend of the Sérsic index
$m$
(the parameter determining the shape of the density profile; see e.g. Ciotti Reference Ciotti1991; Ciotti & Bertin Reference Ciotti and Bertin1999) with the force index is found. Remarkably, the
${\it\alpha}=1$
case (corresponding qualitatively to the force in the deep MOND regime) behaves very similarly to MOND, although this last theory is non-additive (the MOND field equation involves the
$p$
-Laplace operator).
Prompted by these findings and by the above questions, here we set the stage for the study of ROI in the case of power-law forces, by constructing a family of radially anisotropic Hernquist (Reference Hernquist1990) models with additive
$r^{-{\it\alpha}}$
forces, and restricting to the range
$1\leqslant {\it\alpha}<3$
. This interval contains the Newtonian
${\it\alpha}=2$
case, explores ‘short-range’ forces (
${\it\alpha}>2$
) and also long-range forces (
${\it\alpha}<2$
), down to the MOND-like case
${\it\alpha}=1$
. Smaller values of
${\it\alpha}$
are excluded by this preparatory investigation (we recall that the
${\it\alpha}=-1$
case corresponds to the harmonic force, for which no relaxation takes place; see e.g. Lynden-Bell & Lynden-Bell Reference Lynden-Bell and Lynden-Bell1982). Power-law forces with
${\it\alpha}\geqslant 3$
are also excluded by the present analysis, owing to their very peculiar properties (for example, self-gravitating systems with
${\it\alpha}=3$
can be virialized only for zero total energy). We notice that the study of the dynamics under the effect of power-law forces is by no means new: very important examples can be traced back to Newton’s Principia (see e.g. Chandrasekhar Reference Chandrasekhar1995).
This paper is structured as follows. In § 2, we construct isotropic and radially anisotropic self-gravitating spherical Hernquist models with
$r^{-{\it\alpha}}$
forces, and we study their phase-space consistency and the properties of their velocity dispersion profiles. In § 3 we present numerical results on the stability of the isotropic systems. Section 4 summarizes.
2. Radially anisotropic equilibria
The first step in the realization of the
$N$
-body simulations of ROI for systems ruled by generalized power-law forces is the construction of their equilibrium phase-space distribution function
$f$
, in a framework allowing for a tunable amount of radial anisotropy (i.e. different values of
${\it\xi}$
). It is well known that the amount of radial orbits that can be imposed onto a physically acceptable equilibrium stellar system is limited by phase-space consistency (i.e. by the positivity of the associated
$f$
): stability can be studied only for consistent systems. Therefore, before embarking on the study of the stability of a system, it is fundamental to be reassured that its equilibrium initial condition actually exists. For these reasons, in the next section we proceed to construct the phase-space distribution function for the widely used Hernquist (Reference Hernquist1990) profile (a referee pointed out to us that the Hernquist (Reference Hernquist1990) profile has been found earlier as a limiting case of a generalized isochrone model (Kuzmin & Malasidze Reference Kuzmin and Malasidze1970; Kuzmin & Veltmann Reference Kuzmin and Veltmann1973; see also Osipkov Reference Osipkov1978)), under the effect of the
$r^{-{\it\alpha}}$
force.
2.1. The phase-space distribution function
We begin by considering the problem of the construction of the phase-space distribution function
$f$
for a generic equilibrium spherical system of given density profile
${\it\rho}(r)$
, under the action of the
$r^{-{\it\alpha}}$
force with
$1\leqslant {\it\alpha}<3$
. In principle, different velocity anisotropy profiles can be imposed on a given density profile. In this work, in order to build systems with a tunable degree of radial anisotropy, we adopt the widely used Osipkov–Merritt extension of the classical Eddington (Reference Eddington1916) inversion:
$$\begin{eqnarray}f(Q)=\frac{1}{\sqrt{8}{\rm\pi}^{2}}\frac{\text{d}}{\text{d}Q}\int _{Q}^{Q_{\mathit{sup}}}\frac{\text{d}{\it\rho}_{a}}{\text{d}{\it\Phi}}\frac{\text{d}{\it\Phi}}{\sqrt{{\it\Phi}-Q}}=\frac{1}{\sqrt{8}{\rm\pi}^{2}}\int _{Q}^{Q_{\mathit{sup}}}\frac{\text{d}^{2}{\it\rho}_{a}}{\text{d}{\it\Phi}^{2}}\frac{\text{d}{\it\Phi}}{\sqrt{{\it\Phi}-Q}},\end{eqnarray}$$
where for a system of finite total mass
$Q_{\mathit{sup}}=0$
for
${\it\alpha}>1$
, and
$Q_{\mathit{sup}}=\infty$
for
${\it\alpha}\leqslant 1$
(see the following discussion). Here
${\it\Phi}$
is the (self-consistent) gravitational potential,
$Q=E+J^{2}/2r_{a}^{2}$
, where
$E$
and
$J$
are the specific (per unit mass) particle energy and angular momentum modulus, respectively,
$$\begin{eqnarray}{\it\rho}_{a}(r)\equiv \left(1+\frac{r^{2}}{r_{a}^{2}}\right){\it\rho}(r)\end{eqnarray}$$
is the so-called augmented density, and
$r_{a}$
is the anisotropy radius (Osipkov Reference Osipkov1979; Merritt Reference Merritt1985). Owing to the dependence of
$f$
on global integrals of motion, it follows from the Jeans theorem that the equilibrium phase-space distribution function
$f(Q)$
satisfies the Vlasov equation (also known as the collisionless Boltzmann equation; see e.g. Binney & Tremaine Reference Binney and Tremaine2008; Bertin Reference Bertin2014). Integration of
$f(Q)$
over the velocity space after multiplication by the appropriate velocity components shows that the anisotropy profile associated with (2.1) is
$$\begin{eqnarray}{\it\beta}(r)\equiv 1-\frac{{\it\sigma}_{t}^{2}(r)}{2{\it\sigma}_{r}^{2}(r)}=\frac{r^{2}}{r_{a}^{2}+r^{2}},\end{eqnarray}$$
where
${\it\sigma}_{t}(r)$
and
${\it\sigma}_{r}(r)$
are the tangential and radial velocity dispersion profiles of the system (see e.g. Binney & Tremaine Reference Binney and Tremaine2008). In practice, Osipkov–Merritt models are everywhere isotropic for
$r_{a}\rightarrow \infty$
and become more and more radially anisotropic for decreasing values of
$r_{a}$
; at fixed
$r_{a}$
they are isotropic for
$r\ll r_{a}$
and radially anisotropic for
$r\gg r_{a}$
. It is easy to show that the above formulae, derived in the literature for Newtonian gravity, are independent of the specific force law adopted.
From (2.1) it follows that the first step needed to recover
$f$
is knowledge of the self-consistent potential
${\it\Phi}$
of the system. For a generic density
${\it\rho}(\boldsymbol{x})$
, the
${\it\alpha}$
-potential can be written as
$$\begin{eqnarray}{\it\Phi}(\boldsymbol{x})-{\it\Phi}(\boldsymbol{x}_{0})=-\frac{G}{{\it\alpha}-1}\int _{\mathbb{R}^{3}}\left(\frac{1}{\Vert \boldsymbol{x}-\boldsymbol{y}\Vert ^{{\it\alpha}-1}}-\frac{1}{\Vert \boldsymbol{x}_{0}-\boldsymbol{y}\Vert ^{{\it\alpha}-1}}\right){\it\rho}(\boldsymbol{y})\,\text{d}^{3}\boldsymbol{y}\end{eqnarray}$$
for
${\it\alpha}\neq 1$
, and as
$$\begin{eqnarray}{\it\Phi}(\boldsymbol{x})-{\it\Phi}(\boldsymbol{x}_{0})=G\int _{\mathbb{R}^{3}}\ln \left(\frac{\Vert \boldsymbol{x}-\boldsymbol{y}\Vert }{\Vert \boldsymbol{x}_{0}-\boldsymbol{y}\Vert }\right){\it\rho}(\boldsymbol{y})\,\text{d}^{3}\boldsymbol{y}\end{eqnarray}$$
for
${\it\alpha}=1$
. In the equations above,
$G$
is a dimensional coupling constant and
$\boldsymbol{x}_{0}$
is a reference point to be chosen from convergence considerations. The value of
${\it\Phi}(\boldsymbol{x}_{0})$
is usually fixed by arguments of convenience as illustrated below. Note that, when convergence is assured, (2.5) follows from (2.4) for
${\it\alpha}\rightarrow 1$
, provided that
$\Vert \boldsymbol{x}_{0}\Vert <\infty$
.
Elementary integration shows that, for a spherical density distribution
${\it\rho}(r)$
of finite mass, the potential can be written for
$1<{\it\alpha}<3$
as
$$\begin{eqnarray}{\it\Phi}(r)=-\frac{2{\rm\pi}G}{r}\int _{0}^{\infty }\frac{(r+r^{\prime })^{3-{\it\alpha}}-|r-r^{\prime }|^{3-{\it\alpha}}}{({\it\alpha}-1)(3-{\it\alpha})}r^{\prime }{\it\rho}(r^{\prime })\,\text{d}r^{\prime },\end{eqnarray}$$
having set to zero the value of the potential for
$r\rightarrow \infty$
, and assumed
$\Vert \boldsymbol{x}_{0}\Vert \rightarrow \infty$
. For
${\it\alpha}=1$
, some care is needed in order to define
${\it\Phi}(r)$
, and it can be shown that for a system of finite mass the potential can be written as
$$\begin{eqnarray}\displaystyle {\it\Phi}(r)=\frac{{\rm\pi}G}{r}\int _{0}^{\infty }\left[(r+r^{\prime })^{2}\ln \left(\frac{r+r^{\prime }}{r_{n}}\right)-(r-r^{\prime })^{2}\ln \left(\frac{|r-r^{\prime }|}{r_{n}}\right)-2rr^{\prime }\right]r^{\prime }{\it\rho}(r^{\prime })\,\text{d}r^{\prime }, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where now
$\boldsymbol{x}_{0}=0$
, and
$r_{n}$
is an arbitrary normalization length: different choices of
$r_{n}$
just correspond to a different additive constant in the value of
${\it\Phi}$
. As expected, (2.6) and (2.7) are more complicated than their Newtonian counterparts, because Newton’s theorems do not hold for
${\it\alpha}\neq 2$
. However, in general, the potential at large radii for
$1<{\it\alpha}<3$
is dominated by the very simple asymptotic monopole term
$$\begin{eqnarray}{\it\Phi}(r)\sim -\frac{GM}{({\it\alpha}-1)r^{{\it\alpha}-1}},\end{eqnarray}$$
while in the
${\it\alpha}=1$
case the monopole term for
$r\rightarrow \infty$
diverges as
$$\begin{eqnarray}{\it\Phi}(r)\sim GM\ln \frac{r}{r_{n}}.\end{eqnarray}$$
Note that the expansions in (2.8) and (2.9) hold true also in the case of non-spherical density profiles, as can be seen from (2.4) and (2.5) for systems of finite total mass.
In the present work, we focus on systems described by the Hernquist (Reference Hernquist1990) density distribution
$$\begin{eqnarray}{\it\rho}(r)=\frac{M}{2{\rm\pi}r_{c}^{2}r(1+r/r_{c})^{3}},\end{eqnarray}$$
where
$M$
is the total mass,
$r_{c}$
is the core radius and
$r_{\ast }=(1+\sqrt{2})r_{c}$
is the half-mass radius, i.e. the radius of the sphere enclosing half of the total mass.
For the distribution (2.10) and
$1<{\it\alpha}<3$
, the potential at the centre can be evaluated analytically and reads
$$\begin{eqnarray}{\it\Phi}(0)=-\frac{GM}{r_{c}^{{\it\alpha}-1}}\frac{2\text{B}(3-{\it\alpha},{\it\alpha})}{{\it\alpha}-1},\end{eqnarray}$$
where
$\text{B}(x,y)$
is the complete Euler beta function. Curiously, the dimensionless factor in the equation above is a non-monotonic function of
${\it\alpha}$
, with maximum
$-1$
reached for
${\it\alpha}=2$
, and diverging at
$-\infty$
for
${\it\alpha}\rightarrow 1^{+}$
and
${\it\alpha}\rightarrow 3^{-}$
. For
${\it\alpha}=1$
the central potential is instead finite and reads
$$\begin{eqnarray}{\it\Phi}(0)=GM\left(1+\ln \frac{r_{c}}{r_{n}}\right)\!;\end{eqnarray}$$
its numerical value depends on the normalization length
$r_{n}$
. For Hernquist models with
${\it\alpha}=1$
, the natural choice is to adopt
$r_{n}=r_{c}$
. It is not surprising that the limit of (2.8) and (2.11) for
${\it\alpha}\rightarrow 1^{+}$
differs from (2.9) and (2.12), as two different choices of
$\boldsymbol{x}_{0}$
have been adopted for
$1<{\it\alpha}<3$
and
${\it\alpha}=1$
(see the discussion after (2.5)). In general, the Chebyshev (Reference Chebyshev1853) theorem on the integration of differential binomials assures that the integral in (2.6) can be evaluated in terms of elementary functions for the Hernquist model for all rational values of
${\it\alpha}$
. However, the formulae are cumbersome, so we computed (2.6) and (2.7) numerically; as a safety check we verified that the numerical results reproduce the asymptotic expansions at large radii ((2.8) and (2.9)) as well as the values of
${\it\Phi}(0)$
.
2.2. Consistency
In Newtonian gravity
$({\it\alpha}=2)$
, anisotropic Osipkov–Merritt systems of given
${\it\rho}(r)$
are characterized by a critical value of the anisotropy radius,
$r_{ac}$
, such that for
$r_{a}<r_{ac}$
the corresponding models are inconsistent (Ciotti & Pellegrini Reference Ciotti and Pellegrini1992; see also Carollo, de Zeeuw & van der Marel Reference Carollo, de Zeeuw and van der Marel1995; Ciotti Reference Ciotti1996, Reference Ciotti1999; Ciotti & Morganti Reference Ciotti and Morganti2010a
,Reference Ciotti and Morganti
b
; An, van Hese & Baes Reference An, van Hese and Baes2012). It is therefore important for our study to determine the minimum value
$r_{ac}({\it\alpha})$
associated with a nowhere negative (i.e. consistent)
$f(Q)$
. From (2.1) it follows that the phase-space distribution function of an Osipkov–Merritt model (independently of the force law) can be written as
$$\begin{eqnarray}f(Q)=f_{i}(Q)+\frac{f_{a}(Q)}{r_{a}^{2}},\end{eqnarray}$$
where
$f_{i}$
is the phase-space distribution function in the isotropic case (
$r_{a}\rightarrow \infty$
). For the Hernquist density profile, we found that
$f_{i}\geqslant 0$
for all the explored values of
${\it\alpha}$
, i.e. the isotropic Hernquist model is consistent under
$r^{-{\it\alpha}}$
forces, as already known for the Newtonian case. With this being assured, then the critical anisotropy radius is given by
$$\begin{eqnarray}r_{ac}^{2}=\text{sup}\left[-\frac{f_{a}(Q)}{f_{i}(Q)}\right],\quad Q\in A^{-},\end{eqnarray}$$
where the sup on the right-hand side is evaluated over the region
$A^{-}$
of
$Q$
values corresponding to
$f_{a}<0$
(see Ciotti Reference Ciotti2000).
The numerically recovered phase-space distribution functions for
${\it\alpha}=1$
, 1.5, 2 and 2.5 and different values of
$r_{a}$
(with the corresponding values of
${\it\xi}$
) are presented in figure 1. Of course,
${\it\Phi}(0)\leqslant Q<0$
for
$1<{\it\alpha}<3$
, and
${\it\Phi}(0)\leqslant Q\leqslant \infty$
for
${\it\alpha}=1$
. For
${\it\xi}\simeq {\it\xi}_{c}$
, i.e.
$r_{a}\simeq r_{ac}$
(red curves),
$f(Q)$
shows a strong non-monotonicity, owing to the almost dominant negative contribution of
$f_{a}/r_{a}^{2}$
over the isotropic and positive
$f_{i}$
. A further reduction of
$r_{a}$
would deepen even more the depression in the red curves, finally leading to a negative
$f$
, and to phase-space inconsistency. The relation between the monotonicity of
$f(Q)$
and ROI is a very important point that at this stage we cannot address, but that will be central in our forthcoming study (for example by considering the possibility to extend the Antonov laws to
$r^{-{\it\alpha}}$
forces; see Binney & Tremaine Reference Binney and Tremaine2008).
Phase-space distribution functions for Osipkov–Merritt radially anisotropic Hernquist models with
${\it\alpha}=1$
, 1.5, 2 and 2.5. The normalization constant is
$Q_{\ast }={\it\Phi}(r_{\ast })$
, where
$r_{\ast }$
is the half-mass radius. The black curves correspond to the isotropic case (
${\it\xi}=1$
), while the green ones to
${\it\xi}=1.5$
, the fiducial value for the onset of ROI in Newtonian (
${\it\alpha}=2$
) systems. The red curves show
$f(Q)$
for values of
${\it\xi}$
near the consistency limit
${\it\xi}_{c}={\it\xi}(r_{ac})$
(see § 2.2).
In figure 2 we show for
$1\leqslant {\it\alpha}<3$
the values of
$r_{ac}$
, and their associated
${\it\xi}_{c}={\it\xi}(r_{ac})$
. In the Newtonian case (
${\it\alpha}=2$
),
$r_{ac}/r_{\ast }\simeq 0.08$
, in accordance with analytic results (Ciotti Reference Ciotti1996). From figure 2 it is also apparent that systems with lower values of
${\it\alpha}$
(i.e. forces with a longer range than Newtonian) can sustain smaller values of the anisotropy radius (associated with larger values of
${\it\xi}_{c}$
), implying that they are able to store a larger fraction of the kinetic energy in low-angular-momentum orbits maintaining a nowhere negative phase-space distribution function.
Minimum value
$r_{ac}$
of the anisotropy radius for phase-space consistency (in units of the half-mass radius
$r_{\ast }$
) of the Hernquist model with Osipkov–Merritt radial anisotropy and the corresponding value of the Fridman–Polyachenko–Shukhman index
${\it\xi}_{c}={\it\xi}(r_{ac})$
, as functions of the force exponent
${\it\alpha}$
.
Figure 3(a) shows the parameter
${\it\xi}$
as a function of the anisotropy radius
$r_{a}$
for
${\it\alpha}=1$
, 1.5, 2 and 2.5. As expected,
${\it\xi}$
decreases monotonically for increasing
$r_{a}$
for all values of
${\it\alpha}$
. A complementary illustration of this behaviour is given in figure 3(b), where the trend of
${\it\xi}$
with
${\it\alpha}$
for a selection of values of
$r_{a}$
is given.
Fridman–Polyachenko–Shukhman parameter
${\it\xi}$
as a function of the normalized anisotropy radius
$r_{a}/r_{\ast }$
, and of the force exponent
${\it\alpha}$
, for Osipkov–Merritt radially anisotropic Hernquist models with
$r^{-{\it\alpha}}$
forces.
A more detailed representation of the effects of
${\it\alpha}$
and
$r_{a}$
on the internal dynamics of the models can be obtained by plotting the velocity dispersion profiles. The radial velocity dispersion
${\it\sigma}_{r}$
is given by the Jeans equation,
$$\begin{eqnarray}\frac{\text{d}{\it\rho}{\it\sigma}_{r}^{2}}{\text{d}r}+\frac{2{\it\beta}{\it\rho}{\it\sigma}_{r}^{2}}{r}=-{\it\rho}\frac{\text{d}{\it\Phi}}{\text{d}r},\end{eqnarray}$$
whose solution in the case of Osipkov–Merritt anisotropic systems is
$$\begin{eqnarray}{\it\rho}(r){\it\sigma}_{r}^{2}(r)=\frac{\displaystyle \int _{r}^{\infty }\!{\it\rho}(r^{\prime })(r^{\prime 2}+r_{a}^{2})\,\displaystyle \frac{\text{d}{\it\Phi}}{\text{d }r^{\prime }}\,\text{d}r^{\prime }}{r^{2}+r_{a}^{2}}.\end{eqnarray}$$
It is easy to show that, for a density profile decreasing at large radii as
$r^{-{\it\lambda}}$
with
${\it\lambda}>3$
, and in the presence of
$r^{-{\it\alpha}}$
forces,
$$\begin{eqnarray}{\it\sigma}_{r}^{2}(r)\simeq \frac{GM}{r^{2}+r_{a}^{2}}\left(\frac{r^{3-{\it\alpha}}}{{\it\alpha}+{\it\lambda}-3}+\frac{r_{a}^{2}r^{1-{\it\alpha}}}{{\it\alpha}+{\it\lambda}-1}\right).\end{eqnarray}$$
In particular, for
${\it\alpha}=1$
, in the limit
$r\rightarrow \infty$
,
${\it\sigma}_{r}^{2}$
flattens to the constant value
$GM/({\it\lambda}-2)$
owing to the logarithmic nature of the potential for a system of finite mass (see (2.9)). Specifically, for the Hernquist density profile given by (2.10),
${\it\lambda}=4$
, and thus
${\it\sigma}_{r}\simeq \sqrt{GM/2}$
at large radii. In figure 4 we present the radial and tangential velocity dispersion profiles
${\it\sigma}_{r}(r)$
and
${\it\sigma}_{t}(r)$
for
$r_{a}/r_{\ast }=0.5$
, 1 and 2, obtained by numerically solving (2.16). As a further test of the numerical recovery of
$f(Q)$
, we verified that the profiles of
${\it\sigma}_{r}$
and
${\it\sigma}_{t}$
extracted from
$N$
-body realizations obtained by sampling the numerically evaluated
$f(Q)$
reproduce the curves shown in figure 4. It is apparent how for increasing
${\it\alpha}$
, independently of
$r_{a}$
, the components of the velocity dispersion become more and more peaked at inner radii. As expected, for the MOND-like (
${\it\alpha}=1$
) case, the radial profile of
${\it\sigma}_{r}$
is always flat for all the explored values of
$r_{a}$
, in agreement with (2.17) and with the numerical results of Di Cintio et al. (Reference Di Cintio, Ciotti and Nipoti2013). Once the profiles of
${\it\sigma}_{r}^{2}$
and
${\it\sigma}_{t}^{2}$
are known (see (2.3)), it is immediate to evaluate
$K_{r}$
and
$K_{t}$
by numerical integration. In particular, the relatively flat profiles of
${\it\sigma}_{r}$
at large radii for small values of
${\it\alpha}$
are at the origin of the high values of
${\it\xi}$
with respect to those obtained for larger values of
${\it\alpha}$
. The behaviour of
${\it\xi}$
as a function of
$r_{a}$
and
${\it\alpha}$
can also be easily illustrated with the aid of a very simple toy model. In practice, one can integrate (2.17) for a constant density ‘atmosphere’,
${\it\rho}(r)={\it\rho}_{0}$
for
$r\leqslant R$
, under the gravitational field of a central point mass,
$GM/r^{{\it\alpha}}$
. The calculation of the radial velocity dispersion profile is elementary, while
$K_{r}$
in general is given by hypergeometric (Euler beta) functions. However, again from the Chebyshev theorem,
$K_{r}$
can be expressed in terms of elementary functions for all rational values of
${\it\alpha}$
. Then
$K_{t}$
is obtained from the virial theorem as
$K_{t}=-K_{r}+|W|/2$
, where
$W$
is given in the following equation (3.2), which for the toy model considered here gives
$W=-4{\rm\pi}G{\it\rho}_{0}MR^{4-{\it\alpha}}/(4-{\it\alpha})$
.
Radial and tangential velocity dispersion profiles of the self-consistent radially anisotropic Hernquist model with
$r^{-{\it\alpha}}$
forces, for
$r_{a}/r_{\ast }=0.5$
, 1 and 2;
$r_{\ast }$
is the half-mass radius of the distribution. The velocity dispersions are expressed in units of
$v_{\ast }=r_{\ast }/t_{\ast }$
, and
$t_{\ast }$
is given in (3.1). The flat outer parts of the green curve (
${\it\alpha}=1$
) in the top panels are in excellent agreement with the expected value from (2.17), and
${\it\lambda}=4$
.
3. Stability of isotropic systems
As a first sample of the numerical simulations that will be performed in a forthcoming work, we present here direct
$N$
-body simulations aimed at testing the stability of isotropic (
${\it\xi}=1$
) Hernquist models for various values of
${\it\alpha}$
. The initial conditions are generated as in Nipoti et al. (Reference Nipoti, Londrillo and Ciotti2002) with a Monte Carlo sampling of
$f(Q)$
. The equations of motion for
$N=25\,000$
equal-mass particles are integrated with a second-order symplectic scheme with fixed time step
${\rm\Delta}t=t_{\ast }/100$
, where the natural dynamical time scale is defined as
$$\begin{eqnarray}t_{\ast }\equiv \sqrt{\frac{2r_{\ast }^{{\it\alpha}+1}}{GM}}\end{eqnarray}$$
(see Di Cintio & Ciotti Reference Di Cintio and Ciotti2011; Di Cintio et al.
Reference Di Cintio, Ciotti and Nipoti2013). In the code, the divergence of force and potential at vanishing interparticle separation is cured by the introduction of a softening length
${\it\epsilon}$
, so that in practice the potential due to each particle is smoothed as
${\it\Phi}\propto 1/(r^{2}+{\it\epsilon}^{2})^{({\it\alpha}-1)/2}$
for
${\it\alpha}>1$
, and
${\it\Phi}\propto \ln (\sqrt{r^{2}+{\it\epsilon}^{2}})$
for
${\it\alpha}=1$
. The optimal
${\it\epsilon}$
(in units of
$r_{\ast }$
) is fixed so that the softened force on a particle placed at
$\simeq 5r_{\ast }$
from the centre of mass of the system differs less than 0.01 % from the unsoftened force. Therefore the softening length
${\it\epsilon}$
is different for different
${\it\alpha}$
:
${\it\epsilon}=0.05$
, 0.07, 0.1 and 0.125 for
${\it\alpha}=1$
, 1.5, 2 and 2.5, respectively.
The absence of the analogue of the Poisson equation (more specifically
$r^{-{\it\alpha}}$
forces do not obey a Poisson-like equation, belonging to the class of so-called Riesz potentials, associated with fractional Laplace operators (Stein Reference Stein1970)) for
$r^{-{\it\alpha}}$
forces and
${\it\alpha}\neq 2$
points at direct
$N$
-body simulations as the most natural approach for numerical studies. However, it is also important to note the well-known property that potentials associated with the
$r^{-{\it\alpha}}$
force can be expanded in terms of Gegenbauer polynomials for
${\it\alpha}\neq 1$
, and of Chebyshev polynomials for
${\it\alpha}=1$
(Arfken & Weber Reference Arfken and Weber2005), thus opening the way to the construction of tree codes (see e.g. Srinivasan, Mahawar & Sarin Reference Srinivasan, Mahawar and Sarin2005).
Evolution of the fiducial axial ratio
$c/a$
and the virial ratio
$2K/|W|$
for
$N$
-bodysimulations of isotropic Hernquist models with
$r^{-{\it\alpha}}$
forces, and
${\it\alpha}=1$
, 1.5, 2 and 2.5.
In the upper panel of figure 5, we show, for some representative values of
${\it\alpha}$
, the time evolution of the minimum-to-maximum semi-axis ratio
$c/a$
of the systems. The semi-axes are defined from the eigenvalues
$I_{1}\geqslant I_{2}\geqslant I_{3}$
of the tensor
$I_{ij}=\sum _{k}r_{i}^{(k)}r_{j}^{(k)}$
calculated within the radius containing 85 % of the total mass
$M$
(Meza & Zamorano Reference Meza and Zamorano1997; Nipoti et al.
Reference Nipoti, Londrillo and Ciotti2002). In the lower panel of that figure we also show the time evolution of the corresponding virial ratio
$2K/|W|$
, where
$K$
is the total kinetic energy and
$$\begin{eqnarray}\displaystyle W=-\int \!{\it\rho}(\boldsymbol{x})\langle \boldsymbol{x},{\rm\nabla}{\it\Phi}\rangle \,\text{d}^{3}\boldsymbol{x}=\mathop{\sum }_{j\neq i=1}^{N}m_{i}\langle \boldsymbol{x}_{i},\boldsymbol{a}_{ji}\rangle & & \displaystyle\end{eqnarray}$$
is the virial function. The second expression holds for a discrete system of particles with masses
$m_{i}$
, where
$$\begin{eqnarray}\boldsymbol{a}_{ji}=-Gm_{j}\frac{\boldsymbol{x}_{i}-\boldsymbol{x}_{j}}{\Vert \boldsymbol{x}_{i}-\boldsymbol{x}_{j}\Vert ^{{\it\alpha}+1}}\end{eqnarray}$$
is the acceleration on particle
$i$
due to particle
$j$
. We recall that, for
${\it\alpha}\neq 1$
,
$W$
is related to the system total potential energy
$U$
by the identity
$$\begin{eqnarray}W=({\it\alpha}-1)U,\end{eqnarray}$$
where
$$\begin{eqnarray}U=\frac{1}{2}\int \!{\it\rho}(\boldsymbol{x}){\it\Phi}(\boldsymbol{x})\,\text{d}^{3}\boldsymbol{x}=\frac{1}{2}\mathop{\sum }_{i\neq j=1}^{N}m_{i}{\it\Phi}_{ji},\end{eqnarray}$$
and in the case of a discrete system
${\it\Phi}_{ji}$
is the potential on particle
$i$
due to particle
$j$
. For the
${\it\alpha}=1$
case (as for systems in the deep MOND regime) the virial function is independent of time, being
$$\begin{eqnarray}W=-\frac{GM^{2}}{2}\end{eqnarray}$$
for a continuous system, and
$$\begin{eqnarray}W=-\frac{G}{2}\mathop{\sum }_{i\neq j=1}^{n}m_{i}m_{j}\end{eqnarray}$$
for a discrete system (see Di Cintio et al. Reference Di Cintio, Ciotti and Nipoti2013).
From figure 5 it is apparent that the spherical symmetry as well as the equilibrium of the initial conditions is preserved for all the considered values of
${\it\alpha}$
up to
$t/t_{\ast }=40$
, when the simulations end. The fluctuations are of the order of
$\simeq 5\,\%$
in
$c/a$
and
$\simeq 1\,\%$
in
$2K/|W|$
. The stability of the isotropic systems is confirmed also by figure 6, where we compare the final density profile with the initial profile given by (2.10). The initial density profile is preserved down to
$r/r_{\ast }\simeq 0.025$
in the best case (
${\it\alpha}=1.5$
) and to
$r/r_{\ast }\simeq 0.8$
in the worst case (
${\it\alpha}=2.5$
). In all cases the final density profiles are indistinguishable from the initial profiles for
$r\gtrsim 1.4{\it\epsilon}$
. In a forthcoming study we will perform a quantitative monitoring of the deterioration of the central profile as a function of
${\it\alpha}$
,
${\it\epsilon}$
and
$N$
, by considering the time evolution of the Lagrangian radii containing a fixed fraction of the total mass.
Final (
$t/t_{\ast }=40$
) density profiles of the same systems as in figure 5. Here
${\it\rho}_{\ast }={\it\rho}(r_{\ast })$
. The black solid line is the initial density profile given by (2.10). The vertical dashed lines indicate the softening length
${\it\epsilon}$
used in the simulations.
4. Discussion and conclusions
In this preparatory exploration of the ROI in spherically symmetric collisionless equilibrium systems with
$r^{-{\it\alpha}}$
forces, we limited ourselves to a first important step, i.e. the construction of Osipkov–Merritt radially anisotropic self-consistent equilibrium models. We have recovered the potential generated by the Hernquist density distribution under
$r^{-{\it\alpha}}$
forces and we have calculated the corresponding phase-space distribution functions for isotropic and anisotropic Osipkov–Merritt models in the range of values
$1\leqslant {\it\alpha}<3$
. We numerically determined, as a function of
${\it\alpha}$
, the minimum value of the anisotropy radius for consistency
$r_{ac}$
and the associated critical value of the Fridman–Polyachenko–Shukhman index
${\it\xi}_{c}$
, and we computed the radial and tangential velocity dispersion profiles for a few representative consistent systems. The main results of this paper can be summarized as follows:
-
(i) For
$1\leqslant {\it\alpha}<3$
the distribution functions of the isotropic
$r^{-{\it\alpha}}$
Hernquist models are monotonically decreasing and everywhere positive. Their
$N$
-body realizations are numerically stable up to
$\simeq 40t_{\ast }$
. -
(ii) For anisotropic Hernquist models, the minimum value of the anisotropy radius for consistency
$r_{ac}$
increases with increasing
${\it\alpha}$
, and the corresponding
${\it\xi}_{c}$
decreases. Therefore, systems with lower
${\it\alpha}$
can be generated with a higher amount of radial anisotropy. This appears to be in agreement with the results of Di Cintio et al. (Reference Di Cintio, Ciotti and Nipoti2013), who find increasingly anisotropic end products of dissipationless collapses, when reducing
${\it\alpha}$
from 3 to 1. -
(iii) As a general trend, for all values of
${\it\alpha}$
the Fridman–Polyachenko–Shukhman index
${\it\xi}$
decreases for increasing
$r_{a}$
, and increases with decreasing
${\it\alpha}$
at fixed
$r_{a}$
.
In the natural follow-up of this preparatory work we will determine, as a function of
${\it\alpha}$
, the critical
${\it\xi}$
for the onset of ROI, in order to understand how much this process depends on the short- or long-range nature of the interaction law, and we will study the structural and dynamical properties of the final systems.
Acknowledgements
The two referees are acknowledged for important comments that improved the presentation. L.C. and C.N. acknowledge financial support from PRIN MIUR 2010–2011, project ‘The Chemical and Dynamical Evolution of the Milky Way and the Local Group Galaxies’, prot. 2010LY5N2T.
