Two-component Jaffe models with a central black hole. I: the spherical case
L. Ciotti A. Ziaee Lorzad
astro-ph.GA
 
Perfect diode in quantum spin chains
    Dario Poletti
    December 30, 2023
====================================



 

  Dynamical properties of spherically symmetric galaxy models where
  both the stellar and total mass density distributions are described
  by the Jaffe (1983) profile (with different scale-lenghts and
  masses), are presented.  The orbital structure of the stellar
  component is described by Osipkov–Merritt anisotropy, and a black
  hole (BH) is added at the center of the galaxy; the dark matter halo
  is isotropic.  First, the conditions required to have a nowhere
  negative and monothonically decreasing dark matter halo density
  profile, are derived. We then show that the phase-space distribution
  function can be recovered by using the Lambert-Euler W function,
  while in absence of the central BH only elementary functions appears
  in the integrand of the inversion formula. The minimum value of the
  anisotropy radius for consistency is derived in terms of the galaxy
  parameters.  The Jeans equations for the stellar component are
  solved analytically, and the projected velocity dispersion at the
  center and at large radii are also obtained analytically for generic
  values of the anisotropy radius. Finally, the relevant global
  quantities entering the Virial Theorem are computed analytically,
  and the fiducial anisotropy limit required to prevent the onset of
  Radial Orbit Instability is determined as a function of the galaxy
  parameters. The presented models, even though highly idealized,
  represent a substantial generalization of the models presentd in
  Ciotti et al. (2009), and can be useful as starting point for more
  advanced modeling the dynamics and the mass distribution of
  elliptical galaxies.




celestial mechanics – galaxies: kinematics and dynamics – galaxies:
elliptical and lenticular, cD




§ INTRODUCTION


Spherically symmetric galaxy models, despite their simplicity, are
useful tools for theoretical and observational works in Stellar
Dynamics, and for the modelization of stellar systems (e.g., Bertin
2000, Binney & Tremaine 2008). Quite obviously spherical symmetry is
an oversimplification when considering the vast majority of stellar
systems, and a useful spherical model must compensate this limitation
with other features, that make its use preferred or even recommended,
especially in preliminary investigations. Among the important features
of a useful spherical model here we list analytical simplicity,
structural and dynamical flexibility, i.e., possibility to add to the
stellar component a dark matter halo with adjustable density profile,
or alternatively to specify the total density profile, to include the
dynamical effects of a central black hole, to control orbital
anisotropy

For example, the density profile of the stellar distribution of the
model, once projected, should be similar to that of early-type
galaxies, i.e. to the de Vaucouleurs (1948) R^1/4 law, or better,
to its generalization, the so-called R^1/m law (Sersic 1963).
Unfortunately the R^1/m law doesn't allow for an explicit
deprojecton in terms of elementary functions, however the so-called
γ models (Dehnen 1993, Tremaine et al. 1994) in projection are
well fitted over a large radial range, by the R^1/m law. This is
especially true for the Jaffe (1983) and Hernquist (1990) models.

Another important feature of a useful spherical model is the
possibility to reproduce the large scale observational properties of
the total density profile of early-type galaxies.  In fact,
analysis of stellar kinematics (e.g. Bertin et al. 1994, Rix et
al. 1997, Gerhard et al. 2001), as well as several studies combining
stellar dynamics and gravitational lensing support the idea that the
dark and the stellar matter in elliptical galaxies are distributed so
that their total mass profile is described by a density distribution
proportional to r^-2 (e.g., see Treu & Koopmans 2002, 2004; Rusin
et al. 2003; Rusin & Kochanek 2005; Koopmans et al. 2006; Gavazzi et
al. 2007; Czoske et al. 2008; Dye et al. 2008, Nipoti et al. 2008, see
also Shankar et al. 2017). It is clear that simple dynamical models of
two-component galaxies can be useful as starting point of more
sophisticated investigations based on axysimmetric or triaxial galaxy
models (e.g., Cappellari et al. 2007, van den Bosch et
al. 2008). Simple models with flat rotation curve have been in fact
constructed (e.g. Kochaneck 1994, Naab & Ostriker 2007). In
particular we recall the family of two-component galaxy models whose
total mass density is proportional to r^-2, while the visible
(stellar) mass is described by the γ models (Ciotti et
al. 2009, hereafter CMZ09; see also the double power-law models of
Hiotelis 1994). These latter models have been used in
hydrodynamical simulations of accretion onto the central supermassive
black hole (hereafter, BH) in elliptical galaxies (Ciotti & Ostriker
2012, and references therein).  We notice that other models built with
the same approach have been recently applied for the interpretation of
observations (Poci et al. 2017).  We also remark that the approach
used to build these models is different from the standard one, where a
dark matter halo (herefater, DM) is added to the stellar
distribution (e.g. Ciotti & Renzini 1993; Ciotti et al. 1996,
hereafter CRL96; Ciotti 1996, 1999; Sect. 4.4 in CMZ09 ).

A third important feature of a useful spherical model, strictly
related to the previous point, is the possibility to easily compute
the dynamical properties of the stellar component in presence of a
central BH, and possibly to be proved dynamically consistent (see
Sect. 3.1). In fact, supermassive BHs with a mass of the order of
≃ 10^-3 are routinely found at the center of the
stellar sferoids of total mass  (e.g., see Magorrian et al. 1988,
Kormendy & Ho 2013).

Following the arguments above, this paper builds on the CMZ09 model,
and present an even more general (and realistic) class of models,
containing the CMZ09 model as a limit case. On one side, we maintain
the assumption of a Jaffe profile for the stellar distribution, but
now the total density profile is described by another Jaffe law
(instead of a pure r^-2 law), so that the total mass of the models
(that we call JJ models) is finite. At the same time, the scale-lenght
of the total density is a free parameter and so we can reproduce an
r^-2 profile over an arbitrary large radial range. Finally, a
central BH of arbitrary mass (missing in CMZ09 models) is considered
when solving the dynamical equations. For JJ models we show that the
Jeans equations for the stellar component with Osipkov-Merritt
(Osipkov 1979, Merritt 1985a, hereafter OM) radial anisotropy can be
solved analytically, and the projected velocity dispersion at the
center and at large radii can be expressed by means of extremely
simple formulae for generic values of the model parameters. The
positivity of the phase-space density distribution function of the
stellar component (hereafter DF), the so-called consistency, is easily
investigated by using a remarkable property of JJ models, i.e. the
fact that the radial coordinate can be written in terms of the total
potential in terms of the so-called Lambert-Euler W function. By
using this property, we determine the maximum amount of radial
anisotropy allowable for consistency as a function of the galaxy
parameters. These results add to the large amount of phase-space
information already available about one and two-component γ
models (e.g., Carollo et al. 1995, Ciotti 1996, 1999; Baes et
al. 2005, Buyle et al. 2007, Ciotti & Morganti 2009).  As a byproduct
of our analysis we also found that the one-component Jaffe model, at
variance with statements in the literature, cannot be supported
by purely radial orbits.  We note that the W function also appears
in the recenty discovered analytical solution of the isothermal Bondi
accretion problem in Jaffe galaxies with central BH (Ciotti &
Pellegrini 2017), and this fact suggests a first natural application
of JJ models outside the field of Stellar Dynamics, namely in the
field of BH accretion and AGN feedback.

The paper is organized as follows. In Section 2 the main structural
properties of the models are presented. In Section 3 an investigation
of the phase-space properties of the models is carried out both from
the point of view of necessary and sufficient conditions for
consistency, and from direct inspection of the DF. In Section 4 the
solution of the Jeans equation with OM radial anisotropy is presented,
together with the projection of the velocity dispersion profile at
small and large radii.  In Section 5 the important properties related
to the Virial Theorem and global energetic are explicitly calculated,
and the maximum amount of radial anisotropy that can be sustained by
the model without developing Radial Orbit Instability is
estimated. The main results are summarized in Section 6, while more
technical details are given in the Appendix.





§ THE MODELS



As anticipated in the Introduction, the present models are
characterized by a total density distribution (stars plus DM)  described by a
Jaffe (1983) profile; the stellar density distribution  is
also described by a Jaffe profile, in general with a different scale
radius. For future use we recall that the Jaffe density of total mass 
 and scale length  is given by

    (r)= 4π r^2( + r)^2.

The
cumulative mass contained within the sphere of radius r,
and the associated gravitational potential (with the natural condition
of vanishing at infinity, pertintent to systems of finite mass), are
given by 

    (r)= r+r,  (r)=Glnr+r.


Moreover the Jaffe model belongs to the family of γ-models

    (r)=(3-γ) 4π
      r^γ(+r)^4-γ,  
    0≤γ <3,

where  is the total mass,  is a scale-length, and
eq. (1) is obtained for γ=2. The cumulative mass within the
sphere of radius r is given by

    (r)=×(r +r)^3-γ,

so that the half-mass (spatial) radius is =/(2^1
  3-γ-1), and = for the Jaffe model.  For generic
values of γ the projected density at radius R in the
projection plane is given by

    (R)=2∫_R^∞ (r)rdr√(r^2-R^2),

(e.g., Binney & Tremaine 2008), but unfortunately it cannot be
expressed in terms of elementary functions.  However, for γ=2

    (R)=^2×1 4η +
    √(1-η^2)-(2-η^2)  arcsech(η) 2π (1-η^2)^3/2,
                                                            0<η <1 ;
             1 4 - 2 3π,
                                         η=1;
             1 4η-
    √(η^2-1)+(η^2-2) arcsec(η) 2π(η^2-1)^3/2,
                            η >1;

where η≡ R/. In the central and in the very external
regions the projected density profile behave like a power law, with

    (R)∼^2×1 4η,
                                                            R→ 0;
             1 8η^3,
                              R→∞,

respectively.  Finally, an important structural property that we will
consider in the following is the projected mass (R) contained
within the cylinder of radius R. It can be proved that for spherical
systems of finite total mass

    (R)≡2π∫_0^R Σ (R)  R  dR =M-4π∫_R^∞ρ(r)r√(r^2-R^2)dr.

It follows that the projected mass of the Jaffe model is given by (R)=× g(η), where 

    g(η)=η×π 2-η  arcsech(η)√(1-η^2),
                                                            0<η <1 ;
             π 2-1,
                                         η=1;
             π 2-η  arcsec (η)√(η^2-1),
                            η >1.

In particular, the effective radius  of the Jaffe profile
(i.e., the radius in the projection plane encircling half of the total
mass), where g(η_ e)=1/2, is ≃ 0.7447 (in the
Jaffe original paper the slightly erroneous value of 0.763 is
reported).




 §.§ Stellar and total mass distribution



We denote our family of models as “JJ” models, to indicate that it
is a two-component Jaffe model, even though constructed in a different
way with respect to other two-components Jaffe models in the
literature (CLR96, Ciotti 1996, 1999). The properties of the stellar
component are obtained with = and = in eqs. (1)-(9),
while for the galaxy total density distribution (stars plus DM)
= and =.  We adopt  and  as the natural mass
and length scales, and we define

    s≡r,   ξ≡,  ≡ = +1.

From the request that the DM component has a non-negative total mass
 it follows that ≡/≥ 0, and so ≥
1. It is important to note that the request of a non-negative 
does not prevent the possibility of an unphysical, locally negative DM
density. This case will be excluded with the introduction of an
additional constraint, determined in Sect. 2.2.
We also define

    ≡ 4π^3,  ≡G,

as the natural density and potential scales.
With these conventions,
eqs. (1) and (2) for the galaxy model become

    (r)= ξ s^2(ξ+s)^2,

and

    (r)= sξ+s,   (r)=ξlnsξ+s.

We note here an important connection of the JJ models with the models
in CMZ09. In fact, the total galaxy density profile in CMZ09
(eq. [6] therein) can be written as

    ^CMZ(r)=^2 4π G r^2= s^2,

where  is the constant circular velocity. As the total mass
associated with eq. (14) diverges, the parameter =^2/
is not the ratio of the total-to-stellar mass as in JJ
models. An elementary integration shows that =()/,
i.e., it is the total mass contained within the half mass radius
of the Jaffe stellar density profile, normalized to the total stellar
mass.  From eqs. (12) and (14) it follows that the total density
distribution (and the associated quantities, such as the cumulative
and the projected mass profiles, and the force field) of CMZ09 models
can be obtained from JJ models with the substitution

    =ξ

in the corresponding quantities, and then considering the limit for
ξ→∞. Some care is needed for the case of the
potential. In fact JJ models have finite mass and vanishing potential
at infinity, while the logarithimic potential of CMZ09 models

    ^CMZ(r) =^2ln s,

diverges for s→∞.  The proper way to reobtain
^CMZ from eq. (13) is to apply the substitution (15) to
JJ scaled potential +  (lnξ)/ξ, and then to take
the limit ξ→∞.




 §.§ The dark matter distribution: positivity and monotonicity



Before studying the dynamical properties of the models, it is
important to determine the conditions for the positivity and radial
monotonicity of the density distribution of the DM halo. While as
anticipated in Sect. 2.1 the request of positivity is natural, a brief
comment is in order to justify the requirement of monotonicity. In fact,
it can be shown that monotonicity of the density as a function of the
potential is necessary for the positivity of the phase-space
distribution function. From the second Newton's Theorem, the
gravitational potential of a spherical system is necessarily radially
monotone, so that the density profile must be a monotone function of
radius (Ciotti & Pellegrini 1992, hereafter CP92, see also Sect. 3).

As already found in the simpler two component model of CMZ09, also in
JJ models not all values of  and ξ are compatible with a
nowhere negative DM distribution . Curiously, it is
possible to obtain analytically the positivity condition for the more
general family of two-component γ models, built with the same
approach of JJ models. For γγ models the DM
distribution can be written as:

    (r)=(3-γ) s^γ[ξ (ξ+s)^4-γ - 1 (1+s)^4-γ]:

note that  of γγ models is not a γ
  model, unless the stellar and total length scale are equal, so that
in general the local DM-to-stellar mass ratio (r)/ (r)
depends on r. It is easy to verify that the total DM mass
  associated with  is =4π^3 ( -1).





In Appendix A we determine, for given 0≤γ <3, the
conditions on  and ξ to have ≥ 0 for r≥ 0. In the
case of the JJ models (i.e. γγ models with γ=2),
the positivity condition (A2) reduces to

    ≥(ξ)= max(1ξ,ξ);

a DM halo of a model with = is called a minimum
  halo. From equation above it follows that more and more DM is
needed for a total density distribution  more and more
contracted, or more and more expanded than the stellar distribution
. The minimum value =1 can be only adopted when ξ=1,
i.e. when the stellar and total density are proportional, and so
 can vanish everywhere. The situation is illustrated in Fig. 1
(left panel).

As anticipated the positivity of  is just a first condition
for the viability of the model.  A second request is the monotonicity
of  as a function of radius, and this reduces to the
determination of the minimum value  so that d /dr ≤
0. The explicit discussion of this additional restriction is given in
Appendix A, for the whole family of γγ models. In
particular we found that for 1≤γ<3 (the range containing JJ
models, or the analogous two-component Hernquist models), the
positivity and monotonicity conditions for  coincide,
i.e. (ξ) =(ξ).

Equation (18) allows to discuss the relative trend of DM and stars in
JJ models, both at large radii and near the center, as a function of
 and ξ.  For r→∞ and ξ >1 it is easy to show
that ∼ (ξ -1), and so in the outskirts DM and
stars are proportional. When ξ <1 instead the situation is more
complicated: while in non minimum halo models ∼ (/ξ
-1) and so DM and stars distributions are again proportional, in
the minimum halo case ∼ 2(1-ξ)/s∝ r^-5, so
that the galaxy is baryon-dominated in the external regions. The
situation inverts for r→ 0. In fact, in this case for ξ <1 we
have ∼ (/ξ -1) so that DM and stars mass are
locally proportional, but for ξ >1, while in non minimum halo
models ∼ (/ξ-1), in the minimum-halo models
∼ 2(1-1/ξ) s∝ r^-1, so that these models
are centrally baryon-dominated.

It can be of interest for applications to evaluate the relative amount
of dark and visible mass within a prescribed (spatial or projected)
radius. The minimum value for this quantity is easily calculated from
eqs. (2) and (13),

    (r)(r)≥(ξ)(1+s)/ξ+s-1,

where (r)=(r)-(r). In Fig. 1 (middle panel) the mass ratios
corresponding to three representative values of r are shown as a
function of ξ. For example in the case of a sphere of radius equal
to a half mass radius of the stellar distribution (i.e. r=),
the minimum value / is less than unity for ξ >1:
this is a significant improvement of JJ models with respect to the models
of CMZ09, where this ratio can not be less than unity (see Fig.3
therein).

A similar behavior is obtained for the ratio of projected
DM-to-visible mass within some prescribed aperture R, and from the
eq. (9) it is easy to show that

    M_pDM(R) M_ p*(R)≥(ξ) g(η/ξ) g(η)-1.

In Fig. 1 (right panel) we plot this quantity as a function of ξ
for three representative values of the aperture radius, i.e. /2,, and 2. Again the qualitative trend is the same as in
the other panels, with minimum value well below unity for ξ
>1. Note that for R= and considering the limit of eq. (20) for
ξ→∞, we obtain for the mass ratio the value ≃
1.43, in perfect agreement with the analogous result for CMZ09
models.

It is interesting to compare  the DM halo
  profile of JJ models in eq. (17) with the NFW profile (Navarro
  et al. 1997), that we rewrite for r < (the so-called
  truncation radius) as

    (r)=( -1) f(c)s( +s)^2,   f(c)=ln(1+c) -c 1+c,

where ≡ r_ NFW/ is the NFW scale-lenght in units of
, and c≡/r_ NFW: note that in equation above we impose that
the total halo mass  is the same as in eq. (17). From the
asymptotic expansion of  we already know that  and
 at small and large radii cannot in general be similar. Hovever,
in the case of minimum halo with ξ≥ 1, near the center
 increases as 1/r, so that  and  
can be made indentical for r→ 0 with the additional choice

    =√(ξ 2f(c)).

Therefore once a specific JJ minimum halo model is considered and a
radial range fixed, eqs. (21)-(22) allow to determine the best-fit NFW
profile with same total mass and central density profile of 
by tuning the value of c. For example, after a simple
“trial-and-error” exploration, we found that over a range extending
out to ≃ 4-8, a “best-fit” NFW profile can be made to
agree with a minimum halo  with ξ in the range ≃
2-5, with deviations < 10-20% (at large radii), and < 5% inside
≃ 4R_ e, adopting c in the range ≈ 10-20, and
resulting r_ NFW in the range ≈ 0.9 - 1.5.




§ THE PHASE-SPACE  DISTRIBUTION FUNCTION



Having established the structural limitations of the models,
before solving the Jeans equations, it is useful to discuss some basic
property of the phase-space distribution function (hereafter DF) of JJ
models, in order to exclude dynamically inconsistent
combinations of parameters (i.e., choices that would correspond to a
somewhere negative DF). Fortunately, as discussed extensively in CP92
(see also Ciotti 1996, 1999), it is possible to obtain lower bounds
for the OM anisotropy radius as a function of the density slope and
the total mass profile, without actually recovering the DF, which is
in general impossible in terms of elementary functions.  More
specifically, in CP92 a simple theorem was proved regarding the
necessary and sufficient limitations on  in multi-component OM
models. We also recall that the CP92 result has been shown to be just
a very special case of a class of important and more general
inequalities connecting the local density slope and the anisotropy
profile in consistent spherical models (the so-called Global Density
Slope - Anisotropy Inequality, GDSAI, e.g., see de Bruijne et
al. 1996, An & Evans 2006, Ciotti & Morganti 2009, 2010ab, van Hese
et al. 2011).

Thus, following the standard nomenclature (e.g., Binney & Tremaine 2008), we assume
for the stellar component a DF with the OM parameterization

    f=f(Q),   Q≡ - J^2 2^2,

where =-v^2/2 and J are the relative energy and
angular momentum modulus of each star (per unit mass), respectively,
and Ψ =-Φ is the relative potential; moreover the DF is
truncated as f(Q)=0 for
Q<0.  As a central BH of mass  is added at the center of the
galaxy, the total (relative) gravitational potential is
=+G/r, and from eq. (13)


    (r)≡ψ(s) = 
    μ s+ξlnξ+s s,
      μ=.

As well known the radial () and tangential () components of the
velocity dispersion tensor in OM models are related as

    β(r)≡1-^2(r)2^2(r)=r^2r^2+^2,

so that the fully isotropic case is obtained for →∞, while
for =0 the galaxy is supported by pure radial orbits.  For
finite values of , the velocity dispersion tensor becomes
isotropic for r→ 0 (in practice for r<), and fully
radially anisotropic for r→∞ (in practice for r>).
Introducing the augmented density

    ϱ(r)≡(r)(1+r^2^2),

the phase-space DF of the stellar component can be recovered from the
inversion integral

    f(Q)   =   1√(8)π^2d dQ∫_0^Qdϱ dd√(Q-)   =   1√(8)π^2∫_0^Qd^2ϱ d^2d√(Q-);

an analogous expression holds for the DF of the isotropic DM halo,
obtained by using ϱ=, and 
=∞ in eq. (23). 

In the integral above it is intended that ϱ is expressed in 
terms of , and the second identity follows from integration by 
parts when considering spatially untruncated profiles such those of JJ 
models. Note that the OM inversion for the CMZ09 model is somewhat
different (see eqs. [19]-[28] therein, and relative discussion),
because for these latter models Q is not defined in terms of the
relative potential (the potential in eq. [16] is purely logaritmic and
so diverges both r→ 0 and r→∞, making the introduction
of the relative potential useless), and f(Q) is not truncated as a
function of Q.

In Sects. 3.1 and 3.2, after a general discussion about
the limitations on the  imposed by the request of phase-space
consistency, i.e., f(Q)≥ 0 over the accessible phase-space, we
will see how far we can proceed analytically in the recovery of the DF
of JJ models with central BH.




 §.§ Necessary and sufficient conditions for consistency



Following CP92 a necessary condition for the positivity of the
DF of each of the mass components of JJ models (stars or DM) in the
total (galaxy plus central BH) potential is that

    dϱ (r) dr≤ 0  [NC]:

this condition is independent of the behavior of
the other density components of the system. A weak sufficient
  condition for consistency is obtained by requiring that the 
derivative inside the last integral in eq. (27) be positive. Also this
condition can be expressed as a function of radius as

    d dr[dϱ(r) drr^2 (r)]≥ 0,
       [WSC],

where the total mass profile is given by

    (r)=(r)+,

and (r) is given in eq. (13).  Therefore, a model failing
eq. (28) is certainly inconsistent, while a model obeying eq. (29) is
certainly consistent.  It follows that the true boundary in the
parameter space separating consistent and inconsistent models - that
in general can be only determined by direct inspection of the DF - is
“bracketed” by the NC and WSC limits.




Before embarking in the analysis of JJ models, some preliminary
consideration is in order.  First, about the effect of the central BH
on consistency. From eqs. (29) and (30) it follows quite easily that
if 1) the component is consistent for =0, and 2) d (r^2
dϱ/dr)/dr≥ 0, then the model with central BH is
certainly consistent. Note that point 2) is nothing else than the WSC
for the considered density profile interpreted as a tracer in the
gravitational field of the central BH itself; we will use this result
in the following discussion. A second consideration is about the
effect of anisotropy. When dealing with OM anisotropic systems, the
investigation of the NC and WSC, and the study of the DF positivity,
lead to cosider inequalities of the kind

    F+G^2≥ 0,   ≡,

that must hold over the domain  spanned by the functions'
argument.  In practice, the functions F and G are functions of r
(in the case of the NC and WSC) or functions of Q (in the case of
the DF). From inequality (31) it follows that all OM models can be divided in two
families.  When F is nowhere negative over  (e.g., in the
case of a consistent isotropic DF), consistency in the anisotropic
case is obtained for

    ≥≡√( max[0, sup_(-G F)]).

If G is also positive over , then =0 and the system can
be supported by radial orbits only.  In the second case F is
positive only over some proper subset  of , and
negative (or zero) over the complementary subset .  If G <0
somewhere[ In Ciotti (1999) and Ciotti (2000) it is
  erroneously stated that the model is inconsistent if G<0
  everywhere on .  All the results presented therein are
  however correct.]  on , then the condition (31) cannot be
satisfied and the model is inconsistent.  If G ≥ 0 on  one
must consider not only the lower limit  in eq. (32) evaluated
over , but also the condition

    ≤=√( inf_(-G F)),

and consistency is possible only if <.  Summarizing, if
F≥ 0 then ≥ for consistency. If F≤ 0 over
 and G≥ 0 there, then the inequality
≤≤ must be verified. Finally, if over  the
function G <0 somewhere, or <, then inequality (31)
cannot be satisfied and, in case of a DF analysis, the model must be
rejected as inconsistent.

The first application of eqs. (28)-(29) to JJ models concerns the
consistency of the DM halo.  Following the similar analysis in CMZ09,
for simplicity we restrict to the isotropic case, and then eq. (28)
shows the equivalence of the request of monotonicity of 
(Sect. 2.2) with the NC for a consistent DM halo.  Of course, the
restriction to isotropic case is quite arbitrary, as the virialized
end-states of N-body collapses are invariably characterized by some
amount of radial anisotropy (e.g., van Albada 1982; Nipoti, Londrillo
& Ciotti 2006), but for the present illustrative purposes this
assumption is fully justified.  The WSC for a fully isotropic DM halo
is worked out analytically in Appendix A.  In particular, when
restricting to the case of no central BH (μ = 0) we found, quite
surprisingly, that the condition imposed by the WSC to the halo is nothing else than the limit (18) imposed by positivity and
  monotonicity.  It remains to discuss the effect of a central BH.
Following the argument after eq. (30), it is not difficult to show
(Appendix A) that the addition of the central BH in case of isotropy
reinforces consistency, i.e., a DM halo that is consistent in
absence of central BH, it is certainly consistent when a BH is added.
Taken together, the two results above and those in Sect. 2.2 shows
that the isotropic DM halo of JJ models with central BH, once positivity only of  is assured, automatically satisfies the
NC and WSC conditions, and so it is supported by a nowhere negative
phase-space DF.

We now move to the more interesting case of the NC and WSC for the
stellar component of OM anisotropic JJ models. First, we recall that
NC of Jaffe models just reduces to have ≥ 0 while, from the
solution of a cubic equation the WSC for the one-component Jaffe model
gives ≥≃ 0.1068 (Ciotti 1999), marked by the red
solid red circle in Fig. 2.  Second, in Appendix A we show that the
WSC of the stellar component of JJ models is always in the case
described by eq. (32), i.e. only  exists. However, the function
at the r.h.s. of eq. (32) in the general case is sufficiently
complicated that only a numerical study is feasible. In any case, as
in the next Section we will determine the exact limit on 
obtained from the DF, here we just restrict to the case μ=0. The
resulting eq. (A8) is much simpler than the general one, and in
particular =(ξ), i.e. when μ =0 the limit on
anisotropy is independent of  (dotted red line in Fig. 2).  The
red triangle at ≃ 0.0487 marks the position of the WSC
limit for the CMZ09 model obtained by solving a cubic equation, and
that as expected is in accordance with the value of the red line for
ξ→∞. At the opposite limit we have the BH dominated case
(see Appendix A), with ≃ 0.31 marked by the red square,
coincident with the value of the red line for ξ→ 0, when the
total potential becomes that of a central point mass.  In practice,
from the arguments after eq. (30), we have now proved that the stellar
component of JJ models with central BH and OM anisotropy is certainly
consistent for > 0.31, independently of the mass of the central
BH and of the DM halo total mass and scale-lenght.  We conclude this
introductory analysis by noticing the fact that for JJ models, the
presence of a diffuse halo appears to increase the model ability to
sustain radial anisotropy, while for concentrated halos the
consistency of the stellar distribution requires a more isotropic
velocity dispersion tensor, as already found in other two-component OM
models (Ciotti 1996, 1999, CMZ09).




 §.§ Explicit phase-space DF



With the introduction of the dimensionless potential
ψ=/ and augmented density ϱ̃=ϱ/
from eqs. (24) and (26), respectively, eq. (27) writes

    f(q)   =   √(8)π^2^3/2∫_0^qd^2ϱ̃ dψ^2d ψ√(ψ -q)=      √(8)π^2^3/2[U(q)+V(q)^2],

where q≡ Q/. From eqs. (23)-(24) it follows that 0≤
q≤∞.

In Appendix B we show that it is possible to invert eq. (24) and
express analytically the radius as a function of the relative total
potential  by using the Lambert-Euler W function, obtaining

    s(ψ)= ξ Wμ-1,   
    W=W(0, μ e^ξψ +μ),

where W(0,z) is one of the two branches of the real determination of
the complex function W. In absence of the central BH we have
= and it can be shown that eq. (35) reduces to the elementary
function

    s(ψ)=ξ e^ψξ/-1,

in agreement with the solution of eq. (24) with μ =0.  With the
substitution (35) in eq. (26) we finally obtain the expression for
ϱ̃(ψ) to be used in eq. (27). The derivatives inside
the integral are evaluated from the exact relation in eq.  (B3). We
note that the field of application of the W function to physical
problems is rapidly expanding (e.g. see Valluri et al. 2000, Cranmer
et al. 2004, Ciotti & Bertin 2005 for an application to
self-consistent toroidal structures, Veberic 2012, Waters & Proga
2012, Herbst 2015, Ciotti & Pellegrini 2017 for the solution of
isothermal accretion on BHs at the center of galaxies).

In CMZ09 it is shown that for the stellar Jaffe model embedded in a
total singular isothermal density profile, and in absence of the
central BH, the functions U and V can be expressed as simple
linear combinations of exponentials and polylogarithms.  Here, not
surprisingly, the functions U and V cannot be expressed in terms
of known functions, even in absence of the central BH. However, it is
interesting to notice that in case of a dominant central BH (in
practice, sufficiently near to the center), the function f(q) can be
expressed by using simple functions (Appendix C).

We now determine numerically the lower limit for consistency of 
by inspection of the functions U and V.  Note that in absence of
the central BH (μ=0), the variable q in eq. (34) can be further
scaled as =q/, while a factor of ^-3/2 appears in the
functions U and V, as shown in eq. (C4). In particular, for these
models without BH the position of the maximum in eq. (32) depends on
 (and so in terms of q scales linearly with ), but the
value of  is independent of . The same situation occurs in
the CMZ09 models, and in the extreme case of a BH dominated JJ model,
where =q/μ (Appendix C, eqs. [C4]-[C5]), and the scaling
arguments above apply to the DF with  replaced by μ.  It is
numerically found that U≥ 0, so that eq. (32) applies and only
 exists: the solid line in Fig. 2 shows (ξ) determined
by the DF in absence of the central BH, for comparison with the other
curves presented.  Notice how the shape of the critical consistency
curve parallels the WSC condition (red dotted line), and how there are
consistent models failing the WSC. The black circle at ξ =1 marks
the value of the minimum value ≃ 0.02205 for the OM
one-component Jaffe model (Appendix C).  From the figure it i apparent
how the effect of a concentrated DM halo reduces the ability of the
stellar component to sustain radial orbits, while the opposite happens
for models with ξ >1. As an independent test of the derived DF,
the black triangle indicates the limit value of ≃ 0.0141
obtained in CMZ09 by numerical inspection of the DF (coincident, as
expected, with the limit value of the curve for ξ→∞), and
the black square the value of the BH dominated case ≃ 0.082,
coincident with the limit of the curve for ξ→ 0.





In Fig. 3 the DF of the stellar component of a selection of
representative JJ models is presented, in the isotropic (top panel)
and anisotropic (bottom panel, =0.1) cases. In both cases the DFs
are shown with and without the effect of the central BH (with μ =
10^-3), and for illustration, also the BH dominated DF (green
dashed line) is shown. It is clear how at high (relative) energies the
DF of JJ models with central BH is perfectly described by the BH
dominated DF. Also, it is apparent how at high energies the isotropic
and anisotropic DFs for models with the same structure are almost
identical, a property of OM anisotropy parameterization leading to
almost isotropic models in the central regions. It is also important
to note how the DFs of models without the central BH are higher
at high energies than in the analogous models with the central
BH. Also, notice how models with heavier and more extended halos and
so with higher velocity dispersions at large radii (Fig. 4, top panel) at low
relative energies have a lower DF. The physical reason of this
behavior is due to the fact that, qualitatively, the phase-space DF is
inversely proportional to the cube of velocity dispersion (because the
integral over the velocity space, at fixed position, must reproduce
the same value of the local density), so that, empirically, high
velocity dispersions corresponds to low values of the DF. This is
particularly apparent in the BH dominated case, with a low DF at high
energies and a high DF at low energies. This is also confirmed by the
low-energy tail of the DF, which is higher in the anisotropic
cases. In fact, from eqs. (25) and (35) it follows that for
r→∞, the total velocity dispersion profile is proportional to
( A+^2 I)/r^2, i.e., it is lower for smaller values
of . Finally, notice how orbital anisotropy produces a drop of
the DF at intermediate energies, with a depression that would be of
increasing depth for decreasing values of , finally leading to an
inconsistent DF. The curves relative to the anisotropic cases are very
similar to the analogous curves in Ciotti & Lanzoni (1997, Fig. 2),
and C99 (Figs. 2 and 3) and CMZ09 (Fig. 3), revealing the common
qualitative behavior of OM anisotropic DFs near the consistency limit,
i.e. the fact that the inconsistency manifests itself in general at
intermediate energies (see also Ciotti & Morganti 2008 for a
discussion).




§ JEANS EQUATIONS WITH OM ANISOTROPY



The Jeans equations for spherical systems with general (radial or
tangential) anisotropy has been discussed in Binney & Mamon (1982),
and in the OM case the solution can be written as

    (r)^2(r)   =   G r^2 +^2∫_r^∞(r) (r) (1+^2 r^2)dr   =   A(s)+ ^2 I(s) s^2+^2,

where (r) is given in eq. (30),  and the two radial functions

    I =  (s) +μ  (s),   A = 
     (s)+μ  (s),

are the isotropic and purely radial anisotropic components of the
velocity dispersion tensor, respectively.  In the formula above the
dimensionless mass factors,  and μ, have been explicitely
factorized. For →∞ we obtain the solution of the fully
isotropic case, while for =0 we obtain the purely radial case.




 §.§ The velocity dispersion profile



The integration of eq. (37) is elementary. In fact it is formally
equivalent to an integration already performed, for the different
class of two component Jaffe models in CLR96. where the OM Jeans
equation is integrated for a stellar Jaffe distribution, superimposed
to Jaffe DM halo of total mass , and length scale
r_DM. Therefore in CLR96 the combined contribution of the
stars and of the DM potential to the stellar velocity dispersion
profile is given by the sum of two different expressions. Here,
instead, only one integration is required because the total potential is assigned and, in practice, with a suitable
renaming of parameters, the formula in CLR96 for the DM halo
contribution could be used.  However, as we now consider also the
effect of the central BH, not included in the models in CLR96 and
CMZ09, we give the full set of formulae in homogeneous notation.

For the isotropic component

    =ln(ξ+s)ξ^3(ξ-1)^2+(3ξ-4)ln(1+s)
                 (ξ-1)^2-(3ξ^2+2ξ+1)ln sξ^3-2s^2(3ξ^2-ξ-1)+s(3ξ+2)(ξ-1)-ξ(ξ-1) 2ξ^2(ξ-1)s^2(1+s),
             -(6s^2+6s-1)(2s+1) 2s^2(1+s)-6 lns 1+s,

where the first expression holds for ξ 1, and the second for
ξ=1. As expected the two expressions agree with eqs. (A11) and
(A5) in CLR96, respectively[Due to a typo, the sign of the
  terms inside the square brackets of eq. (A11) in CLR96 should be,
  from left to right, plus, plus, minus, minus.].  The contribution of
the BH to the stellar isotropic velocity dispersion profile is given
by

    = 12s^3+6s^2-2s+1 3s^3(1+s)+4 lns 1+s.

Note that this expression could be formally obtained also by
considering the limit for ξ→ 0 of the function , because
from eq. (2)  for fixed r and → 0, becomes the
potential of a point mass.

For the anisotropic part we have

    =ln(ξ+s)ξ
                    (ξ-1)^2+(ξ-2)ln(1+s) (ξ-1)^2-ln
                    sξ-1 (1+s)(ξ-1),
             -2s+3 2(1+s)^2-lns 1+s,

where the first expression holds for ξ 1 and the second for
ξ =1, and they agree with eqs. (A10) and (A4) in CLR96,
respectively. The contribution of the central BH to the anisotropic
stellar velocity dispersion profile is

    =1+2s s(1+s)+2lns 1+s,

and again it is simple to prove that At= for ξ→ 0.
Following eq. (15), we also verified eqs. (39) and (41) considering
the limit for ξ→∞ of the functions ξ and ξ, and recovering eqs. (C2)-(C3) in CMZ09 evaluated for γ=2.

An insight of the behavior of  can be obtained by considering
the expansion for r→∞ and r→ 0 of the obtained formulae.
We begin with the outer galaxy regions.  A simple expansion of the
functions A and I shows that for r→∞ (in
practice, for r >>) the leading order term is the same for the
galaxy as for the BH, with

    ∼∼1 5s^5+O(s^-6),


    ∼∼1 3s^3+O(s^-4).

The coincidence of the leading term is just due to the fact that for
r→∞ the cumulative mass profile in eq. (13) converges to the
total galaxy mass, and for the Newton's theorem this leads to the same
contribution to the velocity dispersion as that of a central mass
.  Following the same approach adopted in CMZ09 we combine
eqs. (43)-(44), and the leading term of  in eq. (37) for r→∞ is obtained, for arbitrary value of , by retaining the
leading order term of the expansion of the much simpler expression

    (r) ^2(r) ∼ (+μ) 5s^2+3^2 15s^5 (s^2+^2).

In the case of finite  we have ^2 ∝ 1/(3s), while
in the fully isotropic case ^2 ∝ 1/(5s): as expected, the
isotropic  is lower than in case of finite , when the
outer regions become populated by radial orbits only. As expected
eq. (45) agrees with the analogous expression obtained for the
two-component model briefly discussed in Sect. 4.4 of CMZ09 (eq. [40]
therein, for γ=2 and for a dominant DM halo). This is at
variance with the behavior of the genuine CMZ09 model, where for
r→∞

    (r)^2(r)∼^22s^2+^2 4s^4(s^2+^2).

Therefore, although the full velocity profile of the CMZ09 model is
recovered from the the limit procedure in eq. (15) applied to
eqs. (39)-(41), the limit procedure applied to eq. (45) does not
converge to the asymptotic expansion of the velocity dispersion
profile in CMZ09. This is due to the fact that for ξ→∞
and r→∞ the integral (37) is not uniform in the
variables ξ and r, so that the two limits cannot be in general exchanged.

The other important region for observational and theoretical
works, is the galaxy center: here the velocity dispersion profile is
dominated by the BH contribution. In fact, at the leading order

    ∼1 2ξ s^2+O(s^-1),  ∼ -ln sξ+O(1),

so that for r→ 0 the galaxy contribution to the stellar velocity
dispersion profile is given by

    (r)^2(r) ∼ξ s^21 2,    >0 ;
             -ln s,
                                =0.

In particular, if =0, the central velocity dispersion diverges as
^2 ∝ -ln s, while for all values >0 the central
velocity dispersion converges to a finite value, coincident with that
of the isotropic case

    ^2 (0)= 2ξ.

This is relevant from the modelistic point of view, as it is well
known that self-gravitating isotropic γ models present a
depression of their velocity dispersion near the center with
(0)=0, except for the γ=0 and γ=2 models (e.g.,
see Bertin et al. 2002 for a general discussion of this phenomenon;
see also Binney & Ossipkov 2001). Notice that the value of the
central velocity dispersion, in the minimum halo model with ξ≥
1 is, according to eq. (18), independent of ξ, and coincident with that of the purely stellar Jaffe model. This shows
the danger of a “blind” use of (0) as a robust indicator of
the actual depth and shape of the galaxy potential well.

For the BH we obtain

    ∼1 3s^3+O(s^-2),  ∼1 s+O(ln s),

and the formula analogous to eq. (48) is

    (r) ^2(r) ∼μ s^31 3,
                                                            >0 ;
             1,
                                =0.

As expected, ^2 diverges as μ /r for r→ 0, and with a
factor of 3 of difference between the fully radially anisotropic case,
and all the other cases with >0 in agreement with the general
property of  in the central regions of γ models with a
BH (e.g., see C96, Baes & Dejonghe 2004, Baes et al. 2005). We
conclude by noticing that eqs. (48), (49) and (51) are also in
accordance with the analogous quantities for the CMZ09 model
(eqs. [C5]-[C6]) and the two-component models in Sect. 4.4 there
(eq. [44]), and with the results in the spherical (isotropic) limit of
one and two-component oblate power-law models with central BH in
Ciotti & Bertin (2005, eq. [C3]) and in Riciputi et al. (2005,
eq. [A4]), evaluated for γ=2



All the relevant properties of  described in this Section are
illustrated in Fig. 4 (top panel) by a selection of representative JJ models. In
particular, the effects of the central BH, of the DM halo, and of
orbital anisotropy, can be clealry seen near the center and at large
radii.




 §.§ Projected velocity dispersion



The projected velocity dispersion profile associated with a general
anisotropy function β(r) is given by

    Σ_*(R)^2(R)= 2∫_R^∞ [1-β(r)R^2r^2 ](r)^2(r) r√(r^2-R^2)dr,

(e.g., Binney & Tremaine 2008), and in the OM case β (r) is given
in eq. (25).

Unsurprisingly the projection integral cannot be evaluated
analytically for JJ models in terms of elementary functions.  However,
as for the spatial velocity dispersion profile interesting
informations can be obtained outside the core radius and near the
center.  In practice, in the external regions the stellar and total
density profiles can be approximated as a pure power-law of slope
-4. In this region the projection integral can be evaluated for
generic values of  and in analogy with eq. (45) the asymptotic
trend with radius of the projected profile can be obtained by
retaining the leading order term of the expansion of

    ^2(R)   ∼   8( +μ) 15πη[1+η^4 2^2(^2+η^2).
       -   .η^4(2^2+η^2)
     archsinh(/η) 2^3(^2+η^2)^3/2],

where η≡ R/. The expression in square parentheses
converges to 1 in the isotropic case, and to 1/3 for all finite
values of .  The analogous formula for the
CMZ09 limit models is

    ^2(R)∼^2(^2+η^2)^5/2-η^3(2^2+η^2)
                                 4^2(^2+η^2)^3/2,

and the same considerations made after eq. (46)  hold.

The case of the central regions is more complicated.  In fact, both
the integral (52) and the projected surface density Σ_* (see
eq. [7]) are asymptotically dominated by their integrands for r→
0, so that ^2 can be properly defined only as the limit for
R→ 0 of the ratio of two diverging quantities. For what concerns
the galaxy contribution, a simple calculation shows that for >0

    (0)=(0),

where (0) is given by eq. (49), again in agreement with
eq. (33) in CMZ09 for γ=2. For =0 instead the central
projected velocity dispersion diverges.Therefore, for the stellar
component of JJ models and >0, the projected central
  velocity dispersion coincides with the central radial component of
  the isotropic velocity dispersion. In presence of the central BH,
 is dominated by the BH contribution, and so it is the
projected velocity dispersion. With some care, from eqs. (51)-(52) it
can be shown that, from eqs. (51) and (52) and independently of
  the value of ≥ 0,

    ^2(R)∼2μ 3πη.

All the relevant properties of  expressed by the formulae in
this Section can be noticed in Fig. 4 (bottom panel), where we show
the projected velocity dispersion profiles for the same JJ models in
the top panel.  In particular Fig. 4 shows a well known consequence of
the OM parameterization, i.e., the fact that the isotropic 
profiles (black lines) in the outer regions are below those in the
corresponding radially anisotropic cases (red lines), while the
opposite holds for the sigp profiles, due to projection effects on
radial orbits in the outer regions, where the l.o.s. direction is
almost perpendicular to the stellar orbits.

We conclude this Section by noticing that CMZ09 (eq. [39]) briefly
commented on the spatial and projected velocity dispersion of a
two-component galaxy model made by the superposition of a stellar
distribution described by a γ model, and a DM halo described by
a Jaffe model.  Of course, when γ=2 this family reduces to JJ
models in CLR96: in turns it is easy to check the perfect
correspondance of eqs. (49) and (55) with eq. (42) in CMZ09 by
assuming there →∞ and β=ξ, when the model becomes
formally identical (in the limiting case of a DM halo “infinitely
massive”) to JJ models (without central BH).  The formulae (55) and
(56) also agree, as expected, with the projection formulae in the
spherical limit of the ellipsoidal models with γ=2 in Ciotti &
Bertin (2005, eqs. [C1] and [C7] therein).




§ VIRIAL, POTENTIAL, AND  KINETIC ENERGIES



Among the several global quantities that are associated with a stellar
system, those entering the Virial Theorem (hereafter VT) are certainly
the most interesting for many observational and theoretical studies
(e.g., Ciotti 2000, Binney & Tremaine 2008). For the stellar
component of JJ models we have

    2≡ - =--,

where =2π∫_0^∞ (^2+^2)r^2dr is the total kinetic energy of the stars,

    = -∫ <𝐱 ,∇>  d^3𝐱
    = - 4π G∫_0^∞ r (r)(r)dr,

is the interaction energy of the stars with the gravitational field of
the galaxy (stars plus DM), and finally

    =-4π G ∫_0^∞ r  (r) dr,

is the interaction energy of the stars with the central BH.
For a Jaffe galaxy  diverges, because the stellar
density profile diverges near the origin as r^-2; instead, this
quantity converges for γ models with 0≤γ<2. Therefore,
the VT implies that also the volume integral of
σ_*BH^2 diverges near the origin for a Jaffe
galaxy, as can be seen by direct integration of eq. (51)

The contribution of the total galaxy potential to 
= + (where  is due to the self-interaction of the
stellar distribution, and  to the effect of the DM halo)
is finite, with the remarkably simple result

    =-ξ-1-lnξ (ξ-1)^2,
                                                           ξ 1;
             1 2,
                                         ξ=1,

and taking the limit as in eq. (15), =-G^2, in accordance
with eq. (33) in CMZ09. More generally it can be shown that  is a
finite quantity for the stellar component of γγ models,
provided 0≤γ < 5/2 (e.g., for two component Hernquist model,
obtained for γ=1).  It follows that for this class of models it
is possible to define the (3-dimensional) galactic virial velocity
dispersion as σ^2_ V=-/: moreover, from eqs. (60),
(55) and (49) the value of σ_ V^2 is proportional to the
value of the central projected velocity dispersion ^2(0), and
the proportionality constant is a function of ξ only: for
ξ=1, σ_ V=(0). We also notice the interesting
behavior of  as a function of ξ in the minimum halo
case. While for increasing ξ≥ 1 it follows that ==ξ
increases correspondingly to arbitrarily large values, the
dimensionless coefficient in eq. (60) just increases from 1/2 for
ξ=1 to 1 for ξ→∞, due to the fact that in minimum-halo
case, more massive halos are necessarily more and more extended, with
a compensating effect on the depth of the total potential.

As well known, in multi component systems the virial energy W of a
given component is not the gravitational energy of the component
itself in the total potential. For this reason we now calculate
explicitly the different contributions to the potential energy 
of the stellar component of JJ models, and we also show how to obtain
the expression of  and  in a simple way.  As for the
interaction energy , also for the potential energy  holds the
decomposition

    =+,

where 

    =+=1 2∫ d^3𝐱+∫ d^3𝐱,

and

    =∫ d^3𝐱=-4π G ∫ r (r) dr=.

Therefore  diverges as . From a well known result, the
self-gravitational energy and the virial self energy of each density
component of a multi-component system coincide, and in our case 
from eq. (60) with =1 and ξ=1,

    =  = - 2,

so that we can compute = - without performing additional
integrations.  The evaluation of  is slightly more complicated,
because in principle it would require to substitute
=- in the second integral in eq. (60), and therefore
compute two integrals. But we adopt a different strategy, and we
compute the integral

    ≡∫ d^3𝐱=-lnξξ-1,
                                                           ξ 1;
             1,
                                     ξ=1.

so that from eq. (60)

    = -2,

and finally  is obtained by adding .

Note that  is not the gravitational energy  of the
stars in the galaxy total potential, . Yet,  is not just an
useful mathematical quantity, but it has an important physical
interpretation, and together  plays a fundamental role in the
theory of galactic winds and X-ray emission of early-type galaxies. In
fact, the energy per unit time to be provided to the ISM of early-type
galaxies (for example by supernova explosion, thermalization of
stellar winds, and AGN feedback) required to steadily extract the mass
losses of stars, injected over the galaxy body at the rate ρ̇_inj=α(t) is given by
L_grav=α(t)|| (e.g, see Pellegrini 2011, 2012,
Posacki et al. 2013). A nice feature of JJ models is that  is
finite and given by a remarkably simple expression, at variance with
the situation of CMZ09 models, where this quantity would diverge, or
other two-component models, where  is given by quite cumbersome
formulae. Therefore JJ models provide a very simple framework to
estimate the energetic of galactic gas flows hosted by X-ray emitting
early-type galaxies.




 §.§ Stability



Another particularly relevant application of the VT is in the field of model
stability, i.e. the determination of the conditions required to
prevent the onset of the so-called Radial Orbit Instability
(hereafter, ROI). In fact, it is well known that stellar systems
supported by a large amount of radial orbits are in general unstable
(e.g., Fridman & Polyachenko 1984, and references therein).  A stability analysis is
obviously well beyond the task of this work, but we can obtain some
quantitative information by investigating the value, as a function
of the model parameters, of the stability indicator

    Ξ≡2  = -4 2+/,

where  and = - are the total kinentic energes
of the stellar component of JJ models, associated with the radial and
tangential components of the velocity dispersion tensor, respectively,
and the last expression is obtained by evaluating  from the VT.
Of course, we exclude the effect of the central BH, due to the formal
divergence of the kinetic energy K_ *BH discussed in previous
Section.  From its definition Ξ→1 for →∞ (globally
isotropic models), while Ξ→∞ for → 0 (fully radially
anisotropic models).

Numerous investigations of one-component systems have confirmed that
the onset of ROI is in general prevented by the empirical requirement
that Ξ < 1.7± 0.25; the exact value of the limit is model
dependent (see, e.g., Merritt & Aguilar 1985; Bertin & Stiavelli
1989; Saha 1991, 1992; Bertin et al. 1994; Meza & Zamorano 1997;
Nipoti, Londrillo & Ciotti 2002).  Here we are considering
two-component systems, however N-body simulations have shown that the
presence of a DM halo does not change very much the situation with
respect to the one-component systems (e.g., see Stiavelli & Sparke
1991, Nipoti et al. 2002). In our case, we assume as a fiducial maximum value for
stability 1.7. 

Note that from eq. (60) and volume integration of eq. (37) with μ
=0, eq. (67) shows that Ξ is independent of .  Unfortunately
 cannot be expressed by using elementary functions, so that we
explore numerically the fiducial stability condition Ξ
(,ξ)=1.7. In Fig. 2 with green dashed curve we plot the
resulting lower bound for stability (ξ).  It is apparent that
the critical value of  increases with ξ, and this is due to
the fact that a spatially extended DM halo increases the contribution
to the kinetic energy of the velocity dispersion in the outer parts
that, in the OM case, are radially anisotropic. Therefore, in order to
guarantee stability in presence of an extended DM halo, the permitted
amount of radial orbits must correspondingly reduced, and larger
values of  are needed. As a limit case the green triangle marks
the position of the stability indicator for the limit models in CMZ09,
with a limit value of ≃ 1.78.  The opposite situation occurs
when the DM halo is more concentrated than the stellar component,
because in this case the velocity dispersion is increased
preferentially in the central regions, that in the OM case are in
practice isotropic, and so a larger amount of radial orbits can be
supported. All these trends nicely agree with those found for
different famillies of one and two-component γ models (Ciotti
1996, 1999, see also Carollo et al. 1995).  We finally notice how the
stability criterion requires minimum anisotropy radii appreciably
larger than those obtained from the consistency analysis (see
Sect. 3.2), and so it is likely that the maximally radially
anisotropic models with positive DF, would be prone to develop ROI.




§ CONCLUSIONS



The family of spherical, two-component galaxy models with the stellar
density distribution described by the Jaffe profile, embedded in a DM
halo such that the total density distribution is also a Jaffe
profile, is presented.  The DM halo is defined as the difference
between the total and the stellar density distributions.  A BH is
added at the center of to the system, and the dynamics of the stellar
component is described by the Osipkov-Merritt anisotropy profile. The
models are fully determined once the total stellar mass () and
scale length () are assigned, together with the total-to-stellar
mass ratio (), the total-to-stellar scale length ratio (ξ),
the BH-to-stellar mass ratio (μ), and finally the anisotropy
radius () of the stellar distribution.  These models represent a
generalization of the CMZ09 models, where the total density profile
was fixed at r^-2 at all radii. In fact JJ models, while retaining
interesting properties such as a realistic stellar density profile and
a total density profile that can be described an arbitrarily large
radial range by a r^-2 profile, have a finite total mass, and a
central BH.  At the same time, they still allow for an almost complete
analytical treatment, and several quantities of interest in
observational and theoretical works have remarkably simple explicit
expressions.  The main results can be summarized as follows.




  * After providing a summary of the structural quantities of
  observational interest for JJ models, for the more general family of
  two-component γγ models, we derive analitically the
  constraints on  and ξ needed to assure positivity and
  monotonicity of the DM halo density distribution.  For a given value
  of ξ, the model corresponding to the minimum value allowed for
   is called minimum halo model.  In JJ models (in which
  the positivity and monotonicity limits coincide), ≥max(ξ, 1/ξ).  Near the origin the density
  profile of the DM halo diverges as ∝ r^-2, but in
  the minimum halo model with ξ >1 the models are centrally
  “baryon dominated”, with ∝ r^-1.


  * It is shown that the models presented in CMZ09 are limit cases
  JJ models (in absence of the central BH), and we provide the
  framework to derive all the structural and dynamical quantities of
  the CMZ09 models from those of JJ models.


  * The minimum value of anisotropy radius , corresponding to a
  dynamically consistent stellar component (i.e., characterized by a
  nowhere negative DF), is first estimated by using the necessary and
  sufficient conditions given in CP92. It is shown that in absence of
  the central BH the minimum value of  so determined is a
  function of ξ only. The critical  decreases for increasing
  ξ, i.e., as already found in other two-component models, a DM
  halo more extended than the stellar distribution increases the
  ability of the stellar component to sustain radial anisotropy. On
  the contrary, more concentrated DM halos (and in particular a
  central BH), require a more isotropic orbital distribution. The
  preliminary consistency analysis is also performed for the DM halo,
  and it is proved that for isotropic DM halos in JJ models with
  μ=0 the conditions of positivity, monotonicity, and phase-space
  consistency coincide; the addition of a central BH reinforces
  consistency.


  * We then moved to study the phase-space DF for the stellar
  component as given by OM inversion.  We found that for JJ models it
  is possible to express analytically the dependence of radius on the
  total potential in terms of the Lambert-Euler W function, allowing
  for a fast and accurate recovery of the DF. In case of no BH
  (μ=0) the resulting expression reduces to elementary functions,
  and in the limit case of a dominat BH (or a very concentrated DM
  halo corresponding to ξ→ 0), the DF itself can be obtained in
  terms of elementary functions. After presenting a few representative
  cases of DFs, corresponding to different choices of , ξ,
  and , we determined numerically the (minimum) critical value of
   as a function of the model parameters, and we found that the
  obtained curve nicely parallels the bound given by the sufficient
  condition in CP92. We showed that in absence of the central BH, and
  in the case of a dominant BH, the critical  depends only on
  ξ, and it is independent of  and μ, respectively. In
  general, we confirmed that DM halos more extended than the stellar
  component increase the amount of radial anisotropy that can be
  supported by a positive DF, while the opposite happens in case of
  concentrated halos (or in presence of a central BH), again in
  accordance with previous findings relative to different
  two-component OM models. Quite unexpectedly, from the inspection of
  the analytical DF, and by independent numerical verification, we
  found that the single component Jaffe model cannot support
  purely radial orbits in the OM formulation, as detailed in Appendix
  C.


  * The Jeans equations for the stellar component are solved
  explicitely for generic values of the model parameters in terms of
  elementary functions.  The asymptotic expansions of  and
   for r→ 0 and r→∞ are obtained, and in particular
  it is shown that when μ=0 and for all values of >0
  (isotropic case included) ^2(0)=/(2ξ). In this
  case, by asymptotic expansion of the projection integral with
  >0, it is also shown that independently of the value of the
  anisotropy radius, (0)=(0). In presence of the BH, in
  the central regions ^2∝ r^-1 with a coefficient which
  is different for =0 or >0. In projection, due to a
  compensating effect, ^2(R)∼ 2μ/(3π R) for
  ≥ 0.


  * Finally, the analytical expressions of relevant quantities
  entering the Virial Theorem, such as the stellar kinetic energy, the
  virial energy interactions, the potential energies, are derived as a
  functions of the model parameters. With the aid of the obtained
  formulae we determined the minimum value of  corresponding to a
  value of ≃ 1.7 of the Friedmann-Poliachenko-Shuckman
  instability indicator, so that more anisotropic models are prone to
  the onset of Radial Orbit Instability. Again, in line with previous
  results, the minimum  for stability increases for increasing
  ξ, and in absence of the central BH its value depends only on
  ξ, being independent of .



We conclude by noting that JJ models, albeit highly idealized, suggest
a few interesting remarks of observational and theoretical character.
For example, after having fixed the properties of the models by using
available observational constraints (e.g., see Negri et al. 2014),
one could use JJ models to investigate how the so called sphere
  of influence of the BH depends on the galaxy properties and how its
definition is affected by orbital anisotropy. Following a preliminary
study (Ziaee Lorzad 2016), it is natural to define the radius of the
sphere of influence as the distance from the galaxy center where the
quantity

    Δσ^2 ≡σ_ *g^2+σ_ *BH^2-σ_ *g^2σ_ *g^2=
    μ (r) (r),

reaches some prescribed value (for example 20%, 50%, 100%) as a
function structural and dynamical properties of the galaxy itself.  JJ
models could also be used to obtain some preliminary estimate of
structural/dynamical properties of high-redshift galaxies (e.g., see
Sect. 4.4.1 in Vanzella et al. 2017), thanks to the very simple
expressions of their virial quantities.

Another interesting application of JJ models is in the field of BH
accretion because, as shown in Ciotti & Pellegrini (2017), it is
possible to solve analytically the generalized isothermal Bondi
accretion problem in Jaffe (or Hernquist) potentials with a central
BH. As the total density profile of JJ models is a Jaffe law, it
follows that for these models we can solve both the accretion
problem for the gas and the Jeans equations for the stellar component.
Moreover, JJ models allow for the computation of the stellar kinetic
energy, a quantity strictly related to the average temperature of the
ISM in early-type galaxies.  As the gas temperature determines the
location of the Bondi radius, JJ models represent a fully analytical
family of self consistent stellar dynamical-hydrodynamical models,
that will allow to compare the relative position of the sonic radius
and the radius of the sphere of influence as a function of the galaxy
properties.




§ ACKNOWLEDGMENTS


We thank the anonymous Referee for useful comments, and John Magorrian
and Silvia Pellegrini for interesting discussions on the models.
L.C. thanks G. Bertin, J. Binney, T. de Zeeuw, W. Evans,
D. Lynden-Bell, D. Merritt and S. Tremaine for useful comments on
Appendix C.  A.Z. acknowledges the Department of Physics and Astronomy
of Padua University, where a preliminary study of JJ models has been
the subject of her Master Thesis.







§ POSITIVITY AND MONOTONICITY OF THE DARK MATTER HALO IN ΓΓ MODELS



The condition for the positivity of the DM halo density profile
 in γγ-models with 0≤γ <3 is
established from eq. (17) as

    ≥(ξ+s)^4-γξ (1+s)^4-γ,      s≥ 0.

Therefore,  must be greater than or equal to the maximum
(ξ,γ) of the radial function at r.h.s.: note that  is
the minimum value of  in order to have a nowhere negative
DM halo. Simple algebra shows that the maximum is attained at infinity
for ξ <1, and at the origin for ξ >1, while for ξ=1 the
radial function is identically equal to 1.  From eq. (A1) it follows
that

    ≥(ξ,γ)=max(1ξ, ξ^3-γ),

and for γ=2 we obtain eq. (18).  

The monotonicity condition for  is obtained by requiring
that d/dr≤ 0, i.e.

    ≥(ξ+s)^5-γ(γ+4s) (1+s)^5-γξ(ξγ+4s),   s≥0.

Again we must determine the maximum (ξ,γ) of the r.h.s. of the equation above.
It is easy to show that for γ=0

    (ξ,0) =max(1ξ,ξ^4),

while for 0<γ<1

    (ξ,γ) =max(1ξ,f_+),

where f_+ is the value of the r.h.s. of eq. (A3) at the critical point

    s_+(ξ,γ)= √(γ[γξ^2+ξ(5-3γ)+γ]) -
        γ (1+ξ) 10:

f_+ →ξ^4 and f_+ →ξ^2 for γ→ 0^+ and γ→ 1^-, respectively. Finally
for 1≤γ <3 (and so in particular for JJ models, or for
two-component Hernquist models that could be constructed by using the
same approach of JJ models) it can be shown, quite surprisingly, that
the monotonicity condition coincides with the positivity condition,
and so  is given by eq. (A2).

The application of the WSC to the isotropic DM halo is obtained from
eq. (26) with →∞, i.e. ϱ=. The condition in absence
of the central BH (μ =0) reduces to  

    ≥2(s+ξ)^3[6s^3 + 4(1+2ξ)s^2  +(1+7ξ)s +2ξ] 
    4 ξ (1+s)^4 (3s^2 +3ξ s +ξ^2),    s≥ 0.

For ξ =1 the r.h.s. equals 1 independently of s. For ξ≠ 1
the determination of the maximum leads to study a fifth degree
equation. Fortunately, it can be proved by inspection that the
resulting expression with s≥ 0 is negative for ξ >1 (and thus
the maximum of eq. [A7] is reached at s=0), and positive for 0<ξ
<1 (and so the maximum is reached for s→∞). In the two
limits eq. (A7) evaluates to ξ and 1/ξ, respectively, and so
we conclude that the isotropic DM halo of JJ models (in absence of
central BH) is certainly consistent when  satisfies the
positivity and monotonicity condition in eq. (18). We are now in
position to consider the effect of the central BH. A direct analysis
would lead to a cumbersome expression, to be explored
numerically. However, by using the considerations after eq. (30), it
is simple to show that the additional term due to the BH is positive,
and so it reinforces the WSC when the positivity condition on 
is verified.

The application of the WSC to the OM anisotropic stellar component of
JJ models leads to the study of a seventh degree equation, and shows
that we are in the conditions pertinent to eq. (32). In absence of the
central BH the dependence on  disappears,

    ^2≥ -
    s^3 [s^2 +2 (ξ -1)s -ξ] 6s^3 +4 (1+2ξ)s^2 + (1+7ξ)s
      +2ξ,    s≥ 0,

and we should solve a fifth degree equation when searching for the
maximum of the r.h.s. In Sect. 3.1 we present the results obtained by
numerical inspection of equation above. Restricting further to the
case ξ =1 (i.e., reducing to the one-component Jaffe model), the
equation to be solved becomes cubic, with ≃ 0.1068 (Ciotti
1999). Finally, the case obtained for ξ→ 0 is formally
coincident with the case of a dominant central BH (i.e., only 
is retained in eqs. [29]-[30]), and for this limiting case the WSC
reduces again to a cubic equation, with solution ≃ 0.31.



§ THE LAMBERT-EULER W FUNCTION


As discussed in Sect.3, for JJ models with central BH it is possible
to invert eq. (24) and express the radius r in terms of the relative
total potential by using the Lambert-Euler W function.  The
integrand in the inversion integral (34) is then obtainedin explicit
and easily tractable form, without resorting to complicate numerical
procedures, because the W function is now fully implemented in the
most used computer algebra systems. The function W(z) (see, e.g.,
Corless et al. 1996) is a multivalued complex function defined
implicitly by the identity

    W e^W=z,

and the two real branches W(0,z) and W(-1,z) for real values of
z are shown in Fig. B1.



With the transformation of variables W=(1+ξ/s)μ/,
eq. (24) can be rewritten as 

    W+ln W = ξψ +μ +lnμ,

so that eq. (35) is obtained by exponentiation of eq. (B2) and
comparison with eq. (B1). It is immediate to conclude that for the
present problem the relevant branch is given by W(0,z), restricting
to points beyond point A.  In fact, when ψ→∞, the
argument in eq. (B2) is infinite, W(0,z)→∞, and from eq. (35)
s→ 0. Instead, when ψ→ 0, the argument tends to μ/
e^μ/, W →μ /, and so from eq. (35) s→∞ from
eq. (35), as it should [By definition of W, it follows that
  W(ye^y)=y.] .  Finally, note that the derivatives inside the
integral (34) can be expressed analytically in terms of W itself,
because from eq. (B1) it follows that

    dW dz=W z(1+W).





§ DF OF OM JAFFE MODEL WITH DOMINANT CENTRAL BH



We report the explicit phase-space DF of a Jaffe model with OM
anisotropy and with dominant central BH (i.e., the gravitational field
is produced by the BH only, and the stellar distribution is only a
tracer). The resulting expression can be interpreted as the asymptotic
limit of the DF at high relative energies, i.e. for galactic regions
sufficiently near the central BH. By using the nomenclature in
eq. (34), it is easy to show that the functions U and V can be
written as

    U(q)=U()μ^3/2,  
    V(q)=V()μ^3/2,  ≡qμ,

where

    U()=(16^3+40^2+18+9)√() 4 (1+q)^3
    -3(3+8) arcsenh√() 4 (1+)^7/2,

and

    V()=(13-2)√() 4 (1+)^3+
    3(1-4) arcsenh√() 4 (1+)^7/2.

The function U is nowehere negative in the range 0≤ <∞,
so the BH dominated Jaffe models are always in the first case
discussed in Sect. 3.1, and only  exists: a numerical evaluation
shows that ≃ 0.082, in agreement with the trend of the
solid line in Fig. 2 for ξ→ 0, when the DM halo “collapses” to
a central point mass.

For completeness we also report the explicit DF for the stellar
component of JJ models with ξ =1 and in absence of the central
BH, when the resulting expression reduces to the one-component DF in the OM
case. From eq. (34) we now have

    U(q)=U()^3/2,  
    V(q)=V()^3/2,  ≡q,

with 0≤ <∞ and

    U() 4√(2)=
     (√(2))+ (√(2)) -
    √(2)[ (√()) +  (√())],


    V() 4√(2)= (√(2)) -  (√())√(2),

where (x)=e^-x^2∫_0^xe^t^2dt is the Dawson's function,
and (x)=e^x^2∫_0^xe^-t^2dt=√(π)e^x^2
  Erf(x)/2. The functions above, when combined according to
eq. (34), are in perfect agreement with those given by Merritt (1985b,
eq. [6]) and Binney & Tremaine (2008).  The function U in eq. (C5)
is positive for alla values of , as shown by the WSC, but the
function V in eq. (C6) becomes negative for admissible values of
, so that  cannot be arbitrarily small.  Numerical
evaluation of eq. (32) shows that for consistency ≥≃
0.02205, in perfect agreement with the solid line in Fig. 2 (obtained
from the general DF) for ξ =1. 

From this result one could conclude that the purely radial model does
not exist. However the situation is not so simple. In fact, the DF of
a purely radial model can be written in all generality as
f=δ(J^2)h(), so that for a finite mass, spatially
untruncated model

    ρ() =2π r^2∫_0^h() d√(
        -),

(e.g., Ciotti 2000), and the inversion formula can be immediately
found (e.g., see Richstone & Tremaine 1984, Oldham & Evans 2016)

    h()   =   1√(2)π^2d d∫_0^ϱ  d√(-)   =   1√(2)π^2∫_0^dϱ dd√(Q-),

where ϱ =r^2ρ is expressed in terms of , and the
second identity follows from integration by parts when considering
spatially untruncated profiles such those of JJ models.  As shown by
Merritt (1985b, eq. 8) and Evans et al. (2015, eq. 31), for the purely
radial one-component Jaffe model 

    h()= 2^2π^2√()[√(2) (√()) -  (√(2))],

where ≡/. The function is positive at all energies, thus showing that the purely radial Jaffe model is consistent. 

These two seamingly contradictory results indicate that the purely
radial case, at least for the Jaffe model, is a singular limit for the
OM parameterization. In practice, we have shown that the non-existence
of the OM (or others) highly radial models cannot by itself exclude
the phase-space consistency of the purely radial configuration.  In
fact, the following argument, built by using the CP92 approach to the
purely radial case, reinforces this conclusion. From the second of
eq. (C8) it follows immediately that a sufficient condition for
consistency of the purely radial model is that the derivative inside
the integral be non-negative, i.e. in terms of radius

    d ϱ (r) dr≤ 0,  ϱ(r)=r^2ρ (r).

Therefore, in the purely radial model a density profile declining as
r^-2 or faster at all radii is a sufficient condition for
consistency (in agreement with the result obtained for the Jaffe
model), while the OM condition (28) (the analogous of eq. (C10) in the
limit of vanishing anisotropy radius) is only necessary for
phase-space consistency. The mathematical reason of the different
behavior is due to the fact that in eq. (C7) , at variance with the
corresponding expression in the OM case, the preparatory derivative of
the augmented density is not required to perform the Abel inversion.





99

 An, J.H., & Evans, N.W. 2006, ApJ, 642, 752

 Baes, M., Dejonghe, H. 2004, MNRAS, 351, 18

 Baes, M., Dejonghe, H., Buyle, P. 2005, A&A, 432, 411

 Bertin, G. 2000, Dynamics of Galaxies,
                   Cambridge University Press

 Bertin, G. &  Stiavelli, M. 1989, ApJ, 338, 723

 Bertin, G., Pegoraro, F., Rubini, F. & 
                   Vesperini, E. 1994, ApJ, 434, 94

 Bertin, G., Bertola, F., Buson, L.M., Danzinger, I.J., Dejonghe, H.,
                 Sadler, E.M., Saglia, R.P., de Zeeuw, P.T., 
                 Zeilinger, W.W. 1994, A&A, 292, 381

 Bertin, G., Ciotti, L., Del Principe, M. 2002, A&A, 386, 149

 Binney, J., & Mamon, G. 1982, MNRAS, 200, 361

 Binney, J., & Ossipkov, L.P. 2001, in 
           Stellar Dynamics: from Classic to Modern, 
           Proceedings of the International Conference, p. 317, 
           eds. Ossipkov, L.P. & Nikiforov, I.I. Sobolev Astronomical 
           Institute, Saint Petersburg)

 Binney, J., & Tremaine, S. 2008,
                 Galactic Dynamics, 2nd Ed.  (Princeton University Press, Princeton)

 de Bruijne, Jos H.J., van der Marel, R.P., de Zeeuw, P.T. 1996,
                 MNRAS, 282, 909

 Buyle, P., Hunter, C., Dejonghe, H.,  2007, MNRAS, 375, 773

 Cappellari, M., Emsellem, E., Bacon, R., Bureau, M.,
           Davies, R.L., de Zeeuw, P.T., Falcon-Barroso, J., Krajnovic, H., 
           Kuntschner, H., McDermid, R.M., Peletier, R.F., Sarzi, M., 
           van den Bosch, R.C.E., van den Ven, G., 2007, MNRAS, 379, 418

 Carollo, C.M., de Zeeuw, P.T., van der Marel, R.P. 1995,
                          MNRAS, 276, 1131

 Ciotti, L. 1996,  ApJ, 471, 68 

 Ciotti, L. 1999,  ApJ, 520, 574 

 Ciotti, L. 2000, Lecture Notes on Stellar Dynamics (Pisa: Scuola Normale Superiore Ed.)

 Ciotti, L., Bertin G. 2005, A&A, 437, 419

 Ciotti, L., Lanzoni, B., & Renzini, A. 1996, MNRAS, 282, 1 (CLR96)

 Ciotti, L., & Lanzoni, B. 1997,  A&A, 321, 724

 Ciotti, L., & Morganti, L. 2009, MNRAS, 393, 179

 Ciotti, L., & Morganti, L. 2010a, MNRAS, 401, 1091

 Ciotti, L., & Morganti, L. 2010b, MNRAS, 408, 1070

 Ciotti, L., & Pellegrini, S. 1992,  MNRAS, 255, 561 (CP92) 

 Ciotti, L., & Pellegrini, S. 2017,
         arXiv:1707.00741

 Ciotti, L., & Renzini, A., 1993, ApJ, 416, L49

 Ciotti, L., & Ostriker, J.P. 2012, AGN Feedback in
 Elliptical Galaxies: Numerical Simulations, in Hot interstellar
   Matter in Elliptical Galaxies, ASSL, vol. 378, p.83,
 eds. D.-W. Kim, S. Pellegrini

 Ciotti, L., Morganti, L., de Zeeuw, P.T.
                 2009, MNRAS, 393, 491 (CMZ09)

 Corless, R.M., Gonnet, G.H., Hare, D.E.G. et al. Adv Comput Math (1996) 5: 329. 

 Cranmer, S. R. 2004, Am. J. Phys, 72(11)

 Czoske, O., Barnabe, M., Koopmans, L.E.V., Treu, T., 
                 Bolton, A.S. 2008, ApJ, 384, 987

 Dehnen, W. 1993, MNRAS, 265, 250

 De Vaucouleurs, G., 1948, Ann.Astr., 11, 247

 Dye, S., Evans, N.W., Belokurow, V., Warren, S.J., Hewett, P. 2008, 
                  MNRAS, 388, 384

 Evans, N.W., An, J., Bowden, A., Williams, A.A. 2015, 
                  MNRAS, 450, 846

 Fridman, A.M., & Polyachenko, V.L. 1984, 
                  Physics of Gravitating Systems (Springer, New York)

 Gavazzi, R., Treu, T., Rhodes, J.D.. Koopmans, L.V.E., 
                  Bolton, A.S., Burles, S., Massey, R.J., Moustakas, L.A. 2007,
                  ApJ, 667, 176

 Gerhard, O., Kronawitter, A., Saglia, R.P., Bender, R. 2001, AJ, 121, 1936

 Herbst, R. S. 2015, PhD Thesis, Univ. of  Witwatersrand

 Hernquist, L. 1990, ApJ, 536, 359
           
 Hiotelis, N. 1994, A&A, 291, 725

 Jaffe, W. 1983, MNRAS, 202, 995

 Kochanek, C.S. 1994, ApJ, 436, 56

 Koopmans, L.V.E., Treu, T., Bolton, A.S., Burles, S., 
                        Moustakas, L.A. 2006, ApJ, 649, 599

 Kormendy, J., & Ho, L.C. 2013, ARA&A, 51, 511

 Magorrian, J., Tremaine, S., Richstone, D., et al. 1998,
       AJ, 115, 2285

 Merritt, D. 1985a, AJ, 90, 1027

 Merritt, D. 1985b, MNRAS, 214, Short Comm. 25P

 Merritt, D. &  Aguilar, L. A. 1985, MNRAS, 217, 787

 Meza, A. & Zamorano, N. 1997, AJ, 490, 136

 Naab, T. & Ostriker, J.P. 2007, MNRAS, 366, 899

 Navarro, J.F.,  Frenk, C.S., White, S.D.M.  1997, ApJ, 490, 493

 Negri, A.,  Posacki, S., Pellegrini, S., & Ciotti, L.,
  2014, MNRAS, 445, 1351

 Nipoti, C., Londrillo, P. & Ciotti, L. 2002, MNRAS, 332, 901

 Nipoti, C., Londrillo, P. & Ciotti, L. 2006, MNRAS, 370, 681

 Nipoti, C., Treu, T. & Bolton, A.S. 2008, MNRAS, 390, 349

 Oldham, L.J., & Evans, N.W. 2016, MNRAS, 462, 298

 Osipkov, L.P. 1979, Pis'ma Astron. Zh., 5, 77

 Poci, A., Cappellari, M., McDermid, R.M., arXiv:1612.05805

 Pellegrini, S. 2011, ApJ, 738, 57

 Pellegrini, S. 2012, Hot Gas Flows on Global and Nuclear Galactic Scales, in Hot interstellar
   Matter in Elliptical Galaxies, ASSL, vol. 378, p.21,
 eds. D.-W. Kim, S. Pellegrini

 Posacki, S., Pellegrini, S., & Ciotti, L., 2013, MNRAS,
 433, 2259

 Richstone, D.O., & Tremaine, S. 1984,  ApJ,   286, 27

 Riciputi, A., Lanzoni, B., Bonoli, S., Ciotti, L. 2005,  A&A,   443, 133

 Rix, H.W., de Zeeuw, P.T., Cretton, N., van der Marel, R.P., 
                  Carollo, C.M. 1997, ApJ, 488, 702
           
 Rusin, D. & Kochanek, C.S. 2005, ApJ, 623, 666

 Rusin, D., Kochanek, C.S., Keeton, C.R.  2003, ApJ, 595, 29

 Saha P., 1991, MNRAS, 148, 494

 Saha P., 1992, MNRAS, 254, 132

 Shankar, F., et al., 2017, ApJ, 840, 34

 Sersic, J.L., 1963, BAAA, 6, 99

 Stiavelli, M. & Sparke, L.S. 1991, ApJ, 382, 466

 Tremaine, S., et al. 1994, AJ, 107, 634

 Treu, T. & Koopmans, L.V.E. 2002, ApJ, 575, 87

 Treu, T. & Koopmans, L.V.E. 2004, ApJ, 611, 739

  Valluri, S. R. &  Jeffrey, D. J. &   Corless, R. M. 2000, Can. J. Physics, 78(9), 823

 van Albada, T.S. 1982, MNRAS, 201, 939

 van den Bosch, R.C.E., van de Ven, G., Verolme, E.K., 
                 Cappellari, M., de Zeeuw, P.T.  2008,   MNRAS, 385,
                 647

 van Hese, E., Baees, M., Dejonghe, H., 2011 ApJ, 726, 80

 Vanzella, E., et al. 2017, MNRAS, 467, 4304

 Veberic, D. 2012, Computer Physics Communications, 183,
        2622

 Waters, T.R., Proga, D., 2012, MNRAS, 426, 2239

 Ziaee Lorzad, A., 2016, Master Thesis, unpublished, Padua University





