Free Access
Issue
A&A
Volume 634, February 2020
Article Number A109
Number of page(s) 12
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201937209
Published online 19 February 2020

© ESO 2020

1. Introduction

The first high-resolution observations of galactic nuclei with the Hubble Space Telescope showed that their surface brightness profiles are generally characterised by power-law behaviour at small radii (Lauer et al. 1991, 1992; Crane et al. 1993; Ferrarese et al. 1994). Lauer et al. (1995) proposed a relatively simple model for the surface brightness profile of galactic nuclei,

(1)

In this formula, R is the circular radius on the plane of the sky, Rb is the break radius that indicates the transition between the inner and the outer profile, Ib = I(Rb) is the surface brightness at the break radius, β and γ correspond to the negative logarithmic surface brightness slopes at large and small radii respectively, and α is a parameter that sets the width of the transition between the inner and outer profiles. This model has become generally known in the extragalactic community as the Nuker model.

The availability of a simple model that accurately describes the surface brightness profile of galactic nuclei is important for many reasons. Numerous studies have used this model to parameterise the surface brightness profiles of the central regions of galaxies (e.g. Byun et al. 1996; Quillen et al. 2000; Rest et al. 2001; Laine et al. 2003; de Ruiter et al. 2005; Lauer et al. 2005, 2007; Durret et al. 2019). Based on this model, a bimodal distribution in the central structure of luminous elliptical galaxies has been discovered, with most galaxies either power-law systems with steep cusps, or core systems with shallow cusps interior to a steeper envelope brightness distribution (Lauer et al. 1995, 2007; Faber et al. 1997).

One disadvantage of the Nuker model, particularly for theoretical studies of the structure and dynamics of galactic nuclei, is that the spatial density distribution that corresponds to the surface brightness profile (1) is not easy to obtain, even in the simplifying assumption of spherical symmetry. Indeed, inserting the surface brightness profile into the standard de-projection formula yields an integral that cannot be evaluated in terms in elementary functions. This is unfortunate, as the spatial density is probably the most fundamental property for dynamical studies (Dejonghe 1986). Most popular families of dynamical models have a simple analytical expression for the spatial mass density as a starting point (e.g. Dehnen 1993; Tremaine et al. 1994; Zhao 1996; Evans & An 2005).

One way to deal with this problem is to use a model that approximates the Nuker model, but in which the density has a simple analytical expression. In particular, Zhao (1996) explored a family of models in which the spatial density profile, rather than the surface brightness profile, is characterised by a double power-law model, which is very similar to what is presented in Eq. (1). This general family of models, also called the generalised Navarro-Frenk-White (NFW) models, is sufficiently simple that the potential and many other interesting dynamical properties can be expressed in terms of elementary or standard special functions. In a follow-up paper, Zhao (1997) provided a recipe to approach the Nuker model with one of the generalised NFW models by matching the surface brightness profiles in a least-squares sense.

Rindler-Daller et al. (2005) took this approach one step further, using a linear combination of simple models to approach the Nuker model. More specifically, they used a quadratic programming technique (Dejonghe 1989) to minimise the difference of the Nuker surface brightness profile and the one corresponding to a linear combination of generalised NFW models. They showed that only a few components are necessary to obtain a good fit. Still, this approach is not guaranteed to have the correct behaviour at all radii, in particular around the transition region between the inner and outer profile.

In this paper, we use a different and direct approach to study the structure and dynamics of the spherical systems characterised by the Nuker surface brightness profile. We make use of an advanced analytical technique based on the Mellin transform method (Marichev 1983; Fikioris 2007). This approach has the advantage that it yields exact analytical expressions for the spatial density and related properties, which allows for a detailed study of the properties and asymptotic behaviour of the model. A similar approach has been adopted to explore the structure and dynamical properties of spherical models with a Sérsic surface brightness profile (Baes & Gentile 2011; Baes & Van Hese 2011; Baes & Ciotti 2019a,b).

This paper is organised as follows. In Sect. 2, we derive an analytical expression for the spatial density profile of the Nuker model using the Mellin integral transform method, we discuss the general behaviour and asymptotic expansion, and we look in detail at the special subclass of models with an infinitely sharp break. In Sect. 3, we discuss some important dynamical properties of the Nuker model, including the potential, velocity dispersion, and phase-space distribution function. In Sect. 4, we discuss a special subset of Nuker models, namely the family of Sérsic models. In Sect. 5, we discuss and summarise our results, where we mainly use the family of Nuker models to refine the classification of galactic nuclei proposed by Tremaine et al. (1994).

2. Spatial density of the Nuker model

2.1. Setting the scene

The Nuker model is characterised by the surface brightness profile (1). As the break radius and the surface brightness at that radius are just scaling parameters, we can consider this family of models a three-parameter family, characterised by the triplet (α, β, γ). In fact, in the remainder of this paper, in order to simplify the expressions, we use normalised units, with G = M = L = Rb = 1, where G is the gravitational constant, M the total mass, and L the total luminosity. The total luminosity of the Nuker model is found by integrating the surface brightness profile over the plane of the sky,

(2)

If α >  0, β >  2 and γ <  2, we find a finite luminosity

(3)

Using this expression and the convention of normalised units, we can rewrite the surface brightness profile of the Nuker model as

(4)

We limit the parameters to the ranges indicated above, with the exception that we only consider positive values of γ to avoid unrealistically increasing surface density profiles at small radii. Hence, in the remainder of this paper, we assume α >  0, β >  2 and 0 ≤ γ <  2.

2.2. Analytical expression for the density

Assuming a constant mass-to-light-ratio, the density profile ρ(r) can be found from the surface brightness profile I(R) through the standard deprojection formula (e.g. Binney & Tremaine 2008),

(5)

This gives

(6)

The integral in this expression can not be readily evaluated using standard methods. To find an analytical expression, we follow a method known as the Mellin transform method (Marichev 1983; Fikioris 2007). While the Mellin transform is a well-known integral transform used in a variety of fields as statistics (Epstein 1948; Fox 1957), analytic number theory (Coffey & Lettington 2015), signal and image analysis (Casasent & Psaltis 1976; Wu et al. 2003), or stellar dynamics (Dejonghe 1986), its application to analytically evaluate integrals is less well-known. The method is surprisingly flexible and powerful, however, and it often yields closed-form expressions that are very difficult to come up with using other methods. We start by rewriting expression (6) in the form

(7)

with

(8)

and

(9)

Expression (7) is a Mellin convolution of the two functions f1 and f2. Similarly to the well-known Fourier convolution theorem, the Mellin transform of a Mellin convolution is equal to the product of the Mellin transforms of the two original functions. Specifically, it implies that we can evaluate the expression (7) as the inverse Mellin transform of the product of the Mellin transforms of f1 and f2. Explicitly, we find

(10)

where 𝔐f(u) denotes the Mellin transform of a function f and the contour integration is along a vertical path ℒ in the complex plane (for more details on Mellin transforms and their inverse, see Fikioris 2007). The Mellin transforms of the functions f1 and f2 can be calculated exactly,

(11)

(12)

This yields the following expression,

(13)

This integral in this expression is a contour integral involving a product of a power function and a number of gamma functions, and such integrals are known as Mellin-Barnes integrals (Barnes 1910; Slater 1966). This particular Mellin-Barnes integral can be written in a compact and elegant form as

(14)

with the Fox H function, generally defined as

(15)

The contour ℒ in this contour integral is chosen such that it separates the poles of the two factors in the numerator. This function was introduced by Fox (1961) as a generalisation of the Meijer G function (Meijer 1946) and shares many of its properties. It is a universal, analytical function that contains many special functions, including generalised Bessel functions, elliptic integrals, generalised hypergeometric functions and the Mittag-Leffler function, as special cases. Different monographs are dedicated to the identities, asymptotic properties, expansion formulae, and integral transforms of the Fox H function (e.g. Mathai et al. 2009).

The numerical evaluation of the Fox H function is challenging and general implementations of the Fox H function are not (yet) available in popular numerical packages as NAG (Phillips 1986), GSL (Galassi et al. 2001), BOOST (Schäling 2014), or SciPy (Oliphant 2007). The Fox H function is a generalisation of the Meijer G function, itself a generalisation of the hypergeometric function. A large variety of techniques exist for the evaluation of hypergeometric functions, including Taylor series expansions, continued fractions, quadrature methods, and more. It turns out, however, that different methods are required for different parameter and argument regimes (for a review, see Pearson et al. 2017). It is, therefore, not surprising that a single implementation for the evaluation of the Fox H function for all values of the parameters and argument is not readily available.

The most direct approach to evaluating the Fox H function is to truncate the general series expansion after a number of terms, based on some stopping criterion. This approach is generally fast, but has the drawback that the exact series expansion depends on the multiplicity of the poles of the integrand in the definition (15) of the Fox H function. If all poles of the defining function are single poles, the Fox H function can be expressed as a power series, whereas it becomes a logarithmic-power series if there are multiple poles (Mathai et al. 2009). Explicit power and power-logarithmic series expansions are given by Kilbas & Saigo (1999) and summarised by Baes & Van Hese (2011). An additional drawback of this approach is that even the simple task of summing different terms in the series expansion can be accomplished in different ways and this can have a significant impact on stability and efficiency.

An alternative method to numerically evaluating the Fox H function is based on a direct numerical integration of the contour integral in Eq. (15). This has been proven to be an interesting approach for the evaluation of general hypergeometric functions (e.g., Press et al. 2002, Sect. 5.14). While it is generally less efficient than more standard approaches as truncated series evaluation, it has the advantage that it often is applicable to a large suite of input parameters and arguments. When adaptive Gaussian or doubly-exponential quadrature formulae are applied to the different segments of this contour, most Fox H functions can be evaluated to relatively high precision. Several experimental numerical implementations for the Fox H function based on this approach are available for various programming languages. Among these are the relatively simple and straightforward implementations in Mathematica and Python by Shafique Ansari et al. (2014) and Alhennawi et al. (2016), respectively, and a GPU-enabled Matlab implementation by Chergui et al. (2019).

The real strength of the Fox H function is its use for analytical studies, however, and as such it is gradually becoming more adopted in applied sciences, including astrophysics (e.g. Saxena et al. 2006; Haubold et al. 2007; Van Hese et al. 2009; Baes & Van Hese 2011; Retana-Montenegro et al. 2012a,b; Zaninetti 2012; Baes & Ciotti 2019a,b.

A direct check on expression (14) can be derived by calculating the total mass (or luminosity). Integrating the density profile over the entire 3D space, we obtain

(16)

Setting s  =  a  =  1 in Eq. (2.8) of Mathai et al. (2009) we can immediately evaluate this integral to recover just M  =  1, as required.

2.3. General behaviour and asymptotic expansion

For all values of the parameters in the ranges we consider (α >  0, β >  2 and 0 ≤ γ <  2), the Nuker model has a non-negative density. The top and middle rows of Fig. 1 show the surface brightness profile and the density profile corresponding to Nuker models with different values of α and γ, and with a fixed value β = 4. A full illustration of the Nuker model profiles for all possible values of the three parameters would be too much information and we have chosen to fix β because its value affects the shape of the surface brightness and density profiles the least. The particular value β = 4 is chosen as this is the value corresponding to the Plummer (1911) model, which is the Nuker model corresponding to (α, β, γ) = (2, 4, 0). For the Plummer model, all of the properties discussed in this paper can be calculated using elementary functions (Dejonghe 1987), which gives us an interesting comparison test.

thumbnail Fig. 1.

Surface brightness profile (top row) and density profile (second row) for the family of Nuker models. The three columns correspond to different values of γ and different lines within each panel correspond to different values of α. Panels on the bottom row: zoom in on the density profile around the break radius r = 1 and contain additional Nuker models with larger values of α.

Open with DEXTER

Comparing the corresponding panels on the first and second rows of Fig. 1, it appears that, to first order, the density profiles are similar to the surface brightness profiles, with a power-law-like behaviour at both small and large radii. We can use the explicit expression for ρ(r) to investigate this in detail. Indeed, under certain conditions, always satisfied for the Fox H functions considered in this paper, the asymptotic expansion of the Fox H function can be calculated using the residue theorem.

For the expansion at large radii, the dominant term is characterised by a single pole for all values of the parameters α, β and γ. Following the recipes outlined in Kilbas & Saigo (1999), we find

(17)

This rβ − 1 power-law behaviour is expected, given the Rβ power-law behaviour for the surface brightness at large projected radii.

For the behaviour at small radii, we can also directly apply the recipes from Kilbas & Saigo (1999). It takes some bookkeeping to determine the order and the multiplicity of the poles, which turn out to depend on the values of both γ and α. Ultimately we find after an extensive calculation

(18)

where E ≈ 0.57721566 is Euler’s constant1, ψ(x) is the digamma function and

(19)

There are several interesting aspects to this asymptotic expansion. A first observation is that all Nuker models with a cuspy surface brightness profile (γ >  0) also have a cusp in their density profile. The opposite is not necessarily true, however: models with γ = 0 and α <  1 do have a weak density cusp with ρ ∝ r−1 + α, but a finite central surface brightness. Nuker models with γ = 0 and α = 1 have a logarithmical density cusp and a finite central surface brightness.

Secondly, the asymptotic behaviour of the density for the subset of Nuker models with a density core, that is, with γ = 0 and α >  1, is still very diverse. It is particularly remarkable to look at the second term in the asymptotic expansion for these models. This term is negative for 1 <  α ≤ 2, indicating that the density decreases as a function of radius. For the models with α >  2, however, this second term has a positive sign, implying that the density increases with increasing radius in the central regions. The higher the value of α, or equivalently, the sharper the transition between the flat inner part of the surface brightness profile and the power-law fall-off beyond the break radius, the stronger this increase: it increases from slightly stronger than linear for α ≳ 2 to quadratic for α >  3. This curious behaviour can easily be spotted in the left panel on the middle row of Fig. 1 and even better in the bottom-left panel, where we zoom in on the region around the break radius r = 1.

2.4. Models with an infinitely sharp break

In the previous subsection, we show that Nuker models with γ = 0 and α >  2 are characterised by a density profile that is not monotonically decreasing as a function of radius. In fact, it turns out that this is not only the case for Nuker models with γ = 0. On the contrary, it is a general feature for all Nuker models, irrespective of the value of β and γ, as long as the value of α is large enough.

The optimal way to illustrate this special feature is to look at the limiting subclass of Nuker models characterised by α → ∞. As indicated in the Introduction section, the parameter α characterises the sharpness of the transition between the inner and outer parts of the surface brightness profile. For small values of α, this transition is smooth and extended, whereas the sharpness of the transition increases with increasing α (see top row of Fig. 1). In the limit α → ∞, the transition is infinitely sharp and the surface brightness profile reduces to a simple broken power-law profile,

(20)

Apart from the “regular” Nuker models, the top panels of Fig. 1 also show the surface brightness profile corresponding to this broken power-law model.

The density profile corresponding to the Nuker models with an infinitely sharp break can in principle be determined by taking the limit α → ∞ in expression (14). It is easier, however, to determine this limit for the explicit Mellin-Barnes expression (13). After some calculation, the resulting expression is a simpler Mellin-Barnes integral that can be written as a Meijer G function,

(21)

Alternatively, we can also substitute expression (20) directly in the deprojection formula (5). After some algebra, we find

(22)

with the function 𝒮λ defined as (see Appendix A)

(23)

The asymptotic expansion at both small and large radii is readily obtained from the asymptotic expansion of the Meijer G function in Eq. (21), or using the expression (22) and the expansion (A.6) of the 𝒮λ function. At large radii, we find a pure power-law profile, as indicated in Eq. (22). At small radii,

(24)

This last expression can also be obtained by directly taking the limit α → ∞ in the expressions (18) and (19).

The advantage of expression (22) is that it is easy to investigate the behaviour of the density profile around the break radius r = 1. For r ≳ 1 we obviously have the pure rβ − 1 power-law fall-off that characterises the entire density profile for these models. For r ≲ 1, however, we obtain, thanks to Eq. (A.7),

(25)

The coefficient of the second term in this expansion is always negative, indicating that, for r ≲ 1, the density profile of all Nuker models with an infinitely sharp break increases with increasing radius.

A consequence of this behaviour is that, for all values of β and γ, there is a maximum value of α such that the density profile of Nuker models with α smaller than this value are monotonically decreasing functions of r, whereas models with a larger value of α have a non-monotonic density profile. This is clearly illustrated in the bottom row of Fig. 1, where we show the density for an extended set of Nuker models that also includes models with high values of α. For models with a surface brightness core (γ = 0), we already demonstrated that all models with α >  2 have an increasing density profile at small radii. The middle and right panels clearly show that also for γ >  0, there is an upturn of the density profile close to the break radius if α sufficiently high. For (β, γ) = (4, 0.6) this happens for α >  13.9, for (β, γ) = (4, 1.2) the critical value is 31.4. The stronger the central cusp of a Nuker model, the sharper the break in the surface brightness profile needs to be to generate an upturn in the density profile.

3. Dynamical properties of the Nuker model

3.1. Gravitational potential

Starting from the density profile (14), we can calculate the corresponding gravitational potential using the formula

(26)

where M(r) is the cumulative mass distribution

(27)

Starting from the Mellin-Barnes form (13) for the density, we find after some algebraic calculation that both the cumulative mass profile and the gravitational potential can also be written in terms of the Fox H function,

(28)

and

(29)

Analysing the latter expression, it is easy to demonstrate that the Nuker models are characterised by a finite potential well if γ <  1,

(30)

This result can also be obtained by plugging Eq. (4) into

(31)

as indicated by Ciotti (1991). Models with γ ≥ 1 have a density cusp stronger than r−2, as can be seen from Eq. (18), and generate an infinitely deep potential well.

Using the formulae for the asymptotic expansion of the Fox H function from Kilbas & Saigo (1999), we can retrieve the explicit behaviour of Ψ(r) at small and large radii. After some calculation, we find at small radii

(32)

At large radii, we find a Keplerian decline, Ψ(r)∼1/r, as expected for a finite mass system.

Figure 2 shows the potential for Nuker models with β = 4 and with different values of α and γ. In all cases, the potential is a smoothly declining function of r. Even the models with a very sharp transition between inner and outer surface brightness profiles, which are characterised by non-monotonic density profile, have a smooth potential.

thumbnail Fig. 2.

Gravitational potential for the family of Nuker models. The selection of models and the meaning of the different lines are as in Fig. 1.

Open with DEXTER

3.2. Isotropic dynamical models

For each spherical potential-density pair, we can construct a unique isotropic dynamical model that generates this density profile using standard equations of galaxy dynamics (e.g. Binney & Tremaine 2008). The velocity dispersion profile σ(r) of the isotropic dynamical model corresponding to a density profile ρ(r) can be found using the solution of the Jeans equation,

(33)

For the general family of Nuker models, this integral cannot be evaluated analytically. We can, however, predict how the velocity dispersion profiles will behave, based on the analysis in Appendix C of Bertin et al. (2002). Nuker models with γ >  1, characterised by a r−(γ + 1) density cusp and an infinitely deep potential well, have a velocity dispersion profile diverging as r−(γ − 1)/2. On the other hand, Nuker models with a weak density cusp, that is, with 0 <  γ <  1 or γ = 0 and α ≤ 1, have a central hole in their velocity dispersion profile. Finally, the models with a finite central density, that is, with γ = 0 and α >  1 have a finite central velocity dispersion. For these models, the central dispersion σ0 can be calculated exactly. After some calculation and Fox H-function manipulation, we find

(34)

On the top row of Fig. 3 we show the velocity dispersion profiles of different Nuker models. The behaviour at small radii is indeed as indicated above: the left panel (γ = 0) shows both dispersion profiles converging to a finite value (α >  1) and profiles with a central hole (α ≤ 1), whereas the models in the central panel have velocity dispersion profiles with a central hole, and the right panels contains models with diverging dispersion profiles. Interesting is also the way the dispersion profiles change as the break become gradually sharper, for fixed values of β and γ. For models with a smooth transition, the dispersion profiles are correspondingly smooth and featureless. When the break gets sharper, the dispersion profile shows a distinctive dive just before the break radius, which corresponds to the increase in the density. The Nuker models with an infinitely sharp break even have a sharp downward peak at r = 1 in their dispersion profile.

thumbnail Fig. 3.

Intrinsic velocity dispersion profiles (top) and line-of-sight velocity dispersion profiles (bottom) for the family of isotropic Nuker models. The selection of models and the meaning of the different lines are as in Fig. 1.

Open with DEXTER

Apart from the intrinsic velocity dispersion, we can also investigate the line-of-sight velocity dispersion σp(R), which is the quantity that is observed in actual galaxies. It can be obtained by projecting the velocity dispersion along the line of sight,

(35)

On the bottom row of Fig. 3 we show the line-of-sight velocity dispersion profiles for the Nuker models. These profiles are similar to the intrinsic dispersion profiles, but there are some subtle differences. The main difference is the asymptotic behaviour of the velocity dispersion profile for subfamily of Nuker models with γ = 0. As argued above, the behaviour of the intrinsic velocity dispersion profile of these models depends on the steepness of the break: models with α >  1 have a dispersion profile that converges to a finite value, whereas models with smaller values of α have a central hole in their intrinsic dispersion profile. On the contrary, the line-of-sight velocity dispersion profiles of all Nuker models with γ = 0 converge to a finite value, even those models with a weak density cusp. This value can be calculated exactly and turns out to be

(36)

A similar difference between intrinsic and line-of-sight velocity dispersion has been noted for other families of models with weak density cusps (Dehnen 1993; Tremaine et al. 1994; Zhao 1996).

Apart from this subtle difference, the behaviour of the intrinsic and line-of-sight velocity dispersion profiles is qualitatively similar. In particular, the change in the shape of the profile for progressively sharper breaks is comparable. For models with a sharp break, the distinctive dip that is clearly visible in the intrinsic dispersion profiles is still present, but it is slightly smoothed out by the line-of-sight integration.

An important caveat on the discussion of the velocity dispersion profiles above is that they are formally derived from the solution of the Jeans equation, but it is not guaranteed a priori that they correspond to physically viable dynamical models. For a dynamical model to be physical or consistent, the distribution function f(r, v) must be positive over the entire phase space. It is well known that, for spherical every potential–density pair, the unique isotropic distribution function is a function of binding energy ℰ only and that it can be calculated using Eddington’s formula (Binney & Tremaine 2008). We have numerically calculated the isotropic distribution function f(ℰ) for all the Nuker models considered before and show the results on the top row of Fig. 4. The corresponding differential energy distributions 𝒩(ℰ), which describe the distribution of orbits with a given binding energy, are plotted on the bottom row.

thumbnail Fig. 4.

Distribution function (top) and differential energy distribution (bottom) for the family of isotropic Nuker models. The selection of models and the meaning of the different lines are as in Fig. 1. Dotted lines correspond to unphysical models, that is, models that have a negative distribution function in some region of phase space.

Open with DEXTER

We have already shown that, for any values for inner and outer slopes of the surface brightness profile, there is a maximum value of α that guarantees that the density is a monotonically decreasing function of radius. Any models with a sharper break have an upturn in their density profile at radii r ≲ 1 and this weird feature translates to a downward peak or a wiggle at similar radii in the dispersion profiles. It is not surprising that these Nuker potential-density pairs do not generate physically viable isotropic dynamical models: a radially decreasing density profile is a necessary condition for the positivity of the distribution function for spherical isotropic models (Ciotti & Pellegrini 1992). In particular, the Nuker models with an infinitely sharp break cannot be generated self-consistently by a physical isotropic distribution function. However, a monotonically decreasing density profile is not a sufficient condition for a positive distribution function and our family of Nuker models shows that explicitly. An example of this is the (α, β, γ) = (8, 4, 0.6) model, characterised by the yellow line in the central columns of Figs. 12 show that this particular model has a smoothly declining surface brightness profile, density profile and gravitational potential, without obvious signatures that would suggest a strange nature. The velocity dispersion profile of this model (top row of Fig. 3) does show a cumbersome curvature around the break radius. Figure 4 shows that the distribution function and differential energy distribution of this model are negative for binding energies between ℰ ∼ 0.81 and ℰ ∼ 0.88 and hence that this model is not a physically viable isotropic model.

In general, it turns out that for every β and γ, the energy structure of the isotropic models varies systematically as α increases. All models with a soft break have positive and monotonically increasing f(ℰ), which implies that the corresponding isotropic dynamical models are not only physical, but also stable to both radial and non-radial perturbations (Antonov 1962; Doremus et al. 1971; Binney & Tremaine 2008). The differential energy distributions of these Nuker models show a broad distribution of orbits over the different binding energies.

As α increases, the distribution 𝒩(ℰ) becomes more peaked and the peak gradually shifts to higher binding energies. When α increases even more, the distribution function starts to show a particular dip at binding energies slightly beyond the peak value. Examples of such models are the (α, β, γ) = (4, 4, 0.6) model, corresponding to the orange line in the central panels of Fig. 4, and the (α, β, γ) = (8, 4, 1.2) model, represented by the yellow line in the right panels. With a distribution function that is no longer a monotonically increasing function of binding energy, these models are not guaranteed to be stable against radial and non-radial perturbations anymore. For γ = 0.6, the critical value of α that separates models with monotonic and non-monotonic distribution functions is 3.8. For γ = 1.2 this value is 4.8.

Increasing α further strengthens the “dip” in the distribution function, until the point where it becomes negative. For γ = 0.6, this point is reached for α = 6.8, for γ = 1.2 at α = 11.4. All models with sharper break in their surface brightness distribution cannot support an isotropic distribution function. These models are indicated as dotted lines in Fig. 4.

4. The Sérsic model as a special case of the Nuker model

The Nuker model is a very flexible model that can fit a wide range of surface brightness profiles. A special subset of Nuker models is obtained using the following four-steps recipe. Firstly, we re-introduce the scaling parameters Rb and L. Secondly, we set , with m and ϵ positive numbers that satisfy . Thirdly, we introduce a new length scale Re = (ϵb)mRb, with b a dimensionless number. Finally, we take the limit ϵ → 0 in the resulting expressions.

Applying this recipe to the Nuker surface brightness profile (4) gives

(37)

To evaluate this expression, we use Wendel’s asymptotic relation (Wendel 1948; Qi & Luo 2012),

(38)

This results in

(39)

This is the surface brightness profile of the Sérsic model (Sérsic 1968), probably the most the popular model to describe the surface brightness distribution of early-type galaxies and bulges of spiral galaxies. As a result of its popularity, the properties of Sérsic model have been examined in great detail (e.g. Ciotti 1991; Ciotti & Lanzoni 1997; Graham & Driver 2005. In particular, we have used similar analytical methods as used in this paper to study the density and other intrinsic properties of the Sérsic model in a series of papers (Baes & Gentile 2011; Baes & Van Hese 2011; Baes & Ciotti 2019a,b).

The analysis above shows that the Sérsic models form a special subfamily of the general Nuker family. We can, therefore, study the properties of the Sérsic model using the general expressions derived in the two previous chapters. In principle we only need to follow the four-steps recipe outlined above. If we apply this recipe to expression (14) for the density, we find after some manipulation and application of Wendel’s asymptotic relation (38),

(40)

Using the explicit Mellin-Barnes form (15) of the Fox H function, this limit can be evaluated as

(41)

in agreement with expression (17) from Baes & Ciotti (2019a). Other expressions can be derived in a similar way and it usually involves an application of Wendel’s asymptotic relation. For example, from Eqs. (18) and (19) we can derive that the Sérsic model has a finite density,

(42)

for m <  1 and a power-law density cusp,

(43)

for m >  1, in agreement with Eqs. (19) and (20) from Baes & Ciotti (2019a). In a similar way, we can derive expressions for the potential of the Sérsic model by applying the four-steps recipe to Eqs. (29)–(32).

5. Discussion and summary

In this paper we construct a three-parameter family of spherical models characterised by a simple surface brightness profile and we discuss the corresponding potential-density pair and the consistency and stability of the corresponding isotropic dynamical models.

We are well aware of the obvious limitations to this set of models. One obvious limitation is that the models are spherically symmetric and hence do not offer the rich orbital structure that more realistic triaxial models possess (e.g. de Zeeuw & Franx 1991; Binney & Tremaine 2008). There are, however, several ways to turn spherically symmetric models into flattened or triaxial models (Schwarzschild 1979, 1993; Hernquist & Quinn 1989; de Zeeuw & Carollo 1996; Baes 2009). A second limitation is that the Nuker models only provide an accurate model for the central regions of galaxies and not for entire galaxies. Other models, such as the core-Sérsic model proposed by Graham et al. (2003) have a broad applicability to model the entire surface brightness profile of early-type galaxies (e.g. Trujillo et al. 2004; Ferrarese et al. 2006; Dullo & Graham 2014). The Nuker model, however, remains an attractive model by virtue of its simple and still flexible form. Finally, even though the expressions (14) and (29) for density and potential are formally written as closed analytical formulae, they still involve a rather obscure special function for which a general numerical implementation is not readily available in the most popular numerical libraries. This makes a practical application of the Nuker models in direct modelling efforts less obvious. The analytical expressions do, however, allow for a detailed analytical investigation of the properties of the Nuker models and we therefore believe that these models can be very useful for theoretical studies concerning galaxies and galactic nuclei.

In particular, they can be used to extend and refine the classification of galactic nuclei proposed by Tremaine et al. (1994). Based on the set of η-models they presented, the authors conclude that three types of central structure are possible. The three types of central structures are denoted as flat core structure (Type I), weak cusp structure (Type II), and strong cusp structure (Type III), respectively. Type I structures, characterised by γ = 0, have a finite central surface brightness, a density cusp as steep at r−1, a finite depth of the central potential well, and an asymptotically constant line-of-sight dispersion. Type II structures, corresponding to 0 <  γ <  1, have a density cusp ranging from r−1 to r−2, a finite potential well, and a central line-of-sight velocity that approaches zero near the centre. Finally, type III structures with 1 <  γ <  2 have a density cusp stronger than r−2, an infinitely deep potential well, and a line-of-sight velocity dispersion profile that diverges near the centre. Tremaine et al. (1994) conclude that this classification does not only apply to their family of η-models, but for any spherical galaxy with an isotropic distribution function and a surface brightness profile that follows a Rγ power-law behaviour at small radii.

One limitation of the study of Tremaine et al. (1994) is that the one-parameter family of models they considered is not fully representative for all central structures, even under the assumption of spherical symmetry and an isotropic dynamical structure. Our three-parameter set of models is more complete and contains a larger variance and we believe that our models fully cover the parameter space of spherical models with an isotropic distribution function and a surface brightness profile that follows a Rγ power-law behaviour at small radii. The analysis in this paper shows that γ is the prime parameter that governs the behaviour of the Nuker models, which supports the classification of central structure based on γ by Tremaine et al. (1994). However, we also demonstrate that γ alone is not sufficient to fully describe the structure of the Nuker models and that α as a secondary parameter has an important role. This is particularly so for the special case of models with a flat core (γ = 0), but also more generally, α turns out to be an important parameter for the consistency and stability of the isotropic dynamical models.

Based on the analysis presented in this paper, we propose to extend the classification of the central structure of spherical models as pictured in Fig. 5. The three main classes introduced by Tremaine et al. (1994), based on the value of γ, still remain. Based on the value of α we add some subdivision.

thumbnail Fig. 5.

Classification of the central structure of spherical galactic nuclei, based on the work of Tremaine et al. (1994). See text for discussion. The current figure corresponds to β = 4, but it has the classification has the same qualitative behaviour for all values of β >  2.

Open with DEXTER

For the flat core structures, we subdivide the range of models in three subclasses. Type Ia structures correspond to γ = 0 and 0 <  α <  1. They have a cored surface brightness profile, but a weak r−1 + α cusp in their density profile. Their isotropic velocity dispersion profiles tend to zero at small radii, but their line-of-sight dispersion profile converge to a finite value. The distribution function of the isotropic dynamical models is positive and monotonic, indicating a stable isotropic dynamical structure. Type Ib structures correspond to γ = 0 and 1 <  α ≤ 2. These systems are characterised by a core in both the surface brightness profile and the density. They have a finite potential density well and a non-zero central intrinsic and line-of-sight velocity dispersion. The isotropic dynamical models are consistent and stable. Finally, all models with γ = 0 and α >  2 unphysical: their surface brightness distribution decreases so slowly that it can only be generated by a density profile that increases as a function of radius. This leads to negative isotropic distribution functions. The border case is formed by the one-parameter family of (α, β, γ) = (2, β, 0) models, which we call the Plummer-like models. The potential-density pair of this sub-family of Nuker models simplifies substantially (Appendix B).

For the Type II and Type III galaxies, a further subdivision can be made based on the monotonicity of the density profile and the consistency and stability of the isotropic distribution function. As discussed in Sect. 2.4, there is, for every couple (β, γ), a limiting value of α above which the density profile shows an increase just short of the break radius. In other words, there is a limit to the sharpness of the break in the surface brightness profile if we want a physical, that is, monotonically decreasing density profile. Furthermore, the in Sect. 3.2 we showed that, for every β and γ, the energy structure of the isotropic models varies systematically as α increases. Going from small to large values of α, we first pass a region of well-behaved models that are consistent and stable. Subsequently we encounter models in which the distribution function f(ℰ) is still positive but no longer monotonically increasing, which turns them possible unstable for radial and non-radial perturbations. Increasing α even more we enter into a regime where the isotropic distribution function becomes negative in some part of phase space, which implies that the isotropic models are unphysical (but anisotropic models might still exist with this density profile). Finally, we enter the previously mentioned region where the density is no longer monotonically decreasing.


1

Euler’s constant is usually denoted by the symbol γ, but for obvious reasons we prefer to use a different symbol.

Acknowledgments

The anonymous referee is acknowledged for his/her very useful and constructive comments that improved the content and presentation of this paper. MB thanks the Fund for Scientific Research Flanders (FWO Vlaanderen) and the Belgian Science Policy Office (BELSPO) for financial support.

References

  1. Alhennawi, H. R., El Ayadi, M. M. H., Alouini, M.-S., Ismail, M. H., & Mourad, H.-A. M. 2016, IEEE Trans. Veh. Technol., 65, 1957 [CrossRef] [Google Scholar]
  2. Antonov, V. A. 1962, Solution of the Problem of Stability of Stellar System (Leningrad: Vestnik Leningradskogo Universiteta) [Google Scholar]
  3. Baes, M. 2009, MNRAS, 392, 1503 [NASA ADS] [CrossRef] [Google Scholar]
  4. Baes, M., & Ciotti, L. 2019a, A&A, 626, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Baes, M., & Ciotti, L. 2019b, A&A, 630, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Baes, M., & Gentile, G. 2011, A&A, 525, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Baes, M., & Van Hese, E. 2011, A&A, 534, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Barnes, E. W. 1910, Q. J. Math., 41, 136 [Google Scholar]
  9. Bertin, G., Ciotti, L., & Del Principe, M. 2002, A&A, 386, 149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press) [Google Scholar]
  11. Byun, Y. I., Grillmair, C. J., Faber, S. M., et al. 1996, AJ, 111, 1889 [NASA ADS] [CrossRef] [Google Scholar]
  12. Casasent, D., & Psaltis, D. 1976, Opt. Eng., 15, 258 [NASA ADS] [CrossRef] [Google Scholar]
  13. Chergui, H., Benjillali, M., & Alouini, M.-S. 2019, IEEE Wirel. Commun. Lett., 8, 428 [CrossRef] [Google Scholar]
  14. Ciotti, L. 1991, A&A, 249, 99 [NASA ADS] [Google Scholar]
  15. Ciotti, L., & Lanzoni, B. 1997, A&A, 321, 724 [NASA ADS] [Google Scholar]
  16. Ciotti, L., & Pellegrini, S. 1992, MNRAS, 255, 561 [NASA ADS] [CrossRef] [Google Scholar]
  17. Coffey, M. W., & Lettington, M. C. 2015, J. Number Theory, 148, 507 [CrossRef] [Google Scholar]
  18. Crane, P., Stiavelli, M., King, I. R., et al. 1993, AJ, 106, 1371 [NASA ADS] [CrossRef] [Google Scholar]
  19. de Ruiter, H. R., Parma, P., Capetti, A., et al. 2005, A&A, 439, 487 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. de Zeeuw, T. 1985, MNRAS, 216, 273 [NASA ADS] [CrossRef] [Google Scholar]
  21. de Zeeuw, P. T., & Carollo, C. M. 1996, MNRAS, 281, 1333 [NASA ADS] [CrossRef] [Google Scholar]
  22. de Zeeuw, T., & Franx, M. 1991, ARA&A, 29, 239 [NASA ADS] [CrossRef] [Google Scholar]
  23. Dehnen, W. 1993, MNRAS, 265, 250 [NASA ADS] [Google Scholar]
  24. Dejonghe, H. 1986, Phys. Rep., 133, 217 [NASA ADS] [CrossRef] [Google Scholar]
  25. Dejonghe, H. 1987, MNRAS, 224, 13 [NASA ADS] [CrossRef] [Google Scholar]
  26. Dejonghe, H. 1989, ApJ, 343, 113 [NASA ADS] [CrossRef] [Google Scholar]
  27. Doremus, J.-P., Feix, M. R., & Baumann, G. 1971, Phys. Rev. Lett., 26, 725 [NASA ADS] [CrossRef] [Google Scholar]
  28. Dullo, B. T., & Graham, A. W. 2014, MNRAS, 444, 2700 [NASA ADS] [CrossRef] [Google Scholar]
  29. Durret, F., Tarricq, Y., Márquez, I., Ashkar, H., & Adami, C. 2019, A&A, 622, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Epstein, B. 1948, Ann. Math. Statist., 19, 370 [CrossRef] [Google Scholar]
  31. Evans, N. W., & An, J. 2005, MNRAS, 360, 492 [NASA ADS] [CrossRef] [Google Scholar]
  32. Faber, S. M., Tremaine, S., Ajhar, E. A., et al. 1997, AJ, 114, 1771 [NASA ADS] [CrossRef] [Google Scholar]
  33. Ferrarese, L., van den Bosch, F. C., Ford, H. C., Jaffe, W., & O’Connell, R. W. 1994, AJ, 108, 1598 [NASA ADS] [CrossRef] [Google Scholar]
  34. Ferrarese, L., Côté, P., Jordán, A., et al. 2006, ApJS, 164, 334 [NASA ADS] [CrossRef] [Google Scholar]
  35. Fikioris, G. J. 2007, Mellin Transform Method for Integral Evaluation: Introduction and Applications to Electromagnetics (Morgan & Claypool Publishers) [Google Scholar]
  36. Fox, C. 1957, Math. Proc. Camb. Philos. Soc., 53, 620 [NASA ADS] [CrossRef] [Google Scholar]
  37. Fox, C. 1961, Trans. Am. Math. Soc., 98, 395 [Google Scholar]
  38. Galassi, M., Gough, B., Rossi, F., et al. 2001, GNU Scientific Library: Reference Manual (Network Theory Limited) [Google Scholar]
  39. Graham, A. W., & Driver, S. P. 2005, PASA, 22, 118 [NASA ADS] [CrossRef] [Google Scholar]
  40. Graham, A. W., Erwin, P., Trujillo, I., & Asensio Ramos, A. 2003, AJ, 125, 2951 [NASA ADS] [CrossRef] [Google Scholar]
  41. Haubold, H. J., Mathai, A. M., & Saxena, R. K. 2007, Bull. Astron. Soc. India, 35, 681 [NASA ADS] [Google Scholar]
  42. Hernquist, L., & Quinn, P. J. 1989, ApJ, 342, 1 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kilbas, A. A., & Saigo, M. 1999, J. Appl. Math. Stoch. Anal., 12, 191 [CrossRef] [Google Scholar]
  44. Laine, S., van der Marel, R. P., Lauer, T. R., et al. 2003, AJ, 125, 478 [NASA ADS] [CrossRef] [Google Scholar]
  45. Lauer, T. R., Faber, S. M., Holtzman, J. A., et al. 1991, ApJ, 369, L41 [NASA ADS] [CrossRef] [Google Scholar]
  46. Lauer, T. R., Faber, S. M., Lynds, R. C., et al. 1992, AJ, 103, 703 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lauer, T. R., Ajhar, E. A., Byun, Y. I., et al. 1995, AJ, 110, 2622 [NASA ADS] [CrossRef] [Google Scholar]
  48. Lauer, T. R., Faber, S. M., Gebhardt, K., et al. 2005, AJ, 129, 2138 [NASA ADS] [CrossRef] [Google Scholar]
  49. Lauer, T. R., Gebhardt, K., Faber, S. M., et al. 2007, ApJ, 664, 226 [NASA ADS] [CrossRef] [Google Scholar]
  50. Marichev, O. I. 1983, Handbook of Integral Transforms of Higher Transcendental Functions: Theory and Algorithmic Tables (Ellis Horwood Publisher) [Google Scholar]
  51. Mathai, A. M., Saxena, R. K., & Haubold, H. J. 2009, The H-Function: Theory and Applications (New York: Springer) [Google Scholar]
  52. Meijer, C. S. 1946, Proc. Nederl. Akad. Wetensch., 98, 227 [Google Scholar]
  53. Oliphant, T. E. 2007, Comput. Sci. Eng., 9, 10 [CrossRef] [PubMed] [Google Scholar]
  54. Pearson, J. W., Olver, S., & Porter, M. A. 2017, Numer. Algorithm, 74, 821 [CrossRef] [Google Scholar]
  55. Phillips, J. 1986, The Nag Library: A Beginner’s Guide (New York: Clarendon Press) [Google Scholar]
  56. Plummer, H. C. 1911, MNRAS, 71, 460 [NASA ADS] [CrossRef] [Google Scholar]
  57. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical Recipes in C++: The Art of Scientific Computing (Cambridge: Cambridge University Press) [Google Scholar]
  58. Qi, F., & Luo, Q.-M. 2012, Banach J. Math. Anal., 6, 132 [CrossRef] [Google Scholar]
  59. Quillen, A. C., Bower, G. A., & Stritzinger, M. 2000, ApJS, 128, 85 [NASA ADS] [CrossRef] [Google Scholar]
  60. Rest, A., van den Bosch, F. C., Jaffe, W., et al. 2001, AJ, 121, 2431 [NASA ADS] [CrossRef] [Google Scholar]
  61. Retana-Montenegro, E., van Hese, E., Gentile, G., Baes, M., & Frutos-Alfaro, F. 2012a, A&A, 540, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Retana-Montenegro, E., Frutos-Alfaro, F., & Baes, M. 2012b, A&A, 546, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Rindler-Daller, T., Dejonghe, H., & Zeilinger, W. W. 2005, MNRAS, 356, 1403 [NASA ADS] [CrossRef] [Google Scholar]
  64. Saxena, R. K., Mathai, A. M., & Haubold, H. J. 2006, Ap&SS, 305, 289 [NASA ADS] [CrossRef] [Google Scholar]
  65. Schäling, B. 2014, The Boost C++ Libraries (Fort Collins: XML Press) [Google Scholar]
  66. Schwarzschild, M. 1979, ApJ, 232, 236 [NASA ADS] [CrossRef] [Google Scholar]
  67. Schwarzschild, M. 1993, ApJ, 409, 563 [NASA ADS] [CrossRef] [Google Scholar]
  68. Sérsic, J. L. 1968, Atlas de Galaxias Australes (Cordoba, Argentina: Observatorio Astronomico) [Google Scholar]
  69. Shafique Ansari, I., Yilmaz, F., Alouini, M.-S., & Kucur, O. 2014, Trans. Emerg. Telecommun. Technol., 28, 1 [Google Scholar]
  70. Slater, L. J. 1966, Generalized Hypergeometric Functions (Cambridge: Cambridge University Press) [Google Scholar]
  71. Tremaine, S., Richstone, D. O., Byun, Y.-I., et al. 1994, AJ, 107, 634 [NASA ADS] [CrossRef] [Google Scholar]
  72. Trujillo, I., Erwin, P., Asensio Ramos, A., & Graham, A. W. 2004, AJ, 127, 1917 [NASA ADS] [CrossRef] [Google Scholar]
  73. Van Hese, E., Baes, M., & Dejonghe, H. 2009, ApJ, 690, 1280 [NASA ADS] [CrossRef] [Google Scholar]
  74. Wendel, J. G. 1948, Am. Math. Mon., 55, 563 [CrossRef] [MathSciNet] [Google Scholar]
  75. Wu, C., Somervell, A., Haskell, T., & Barnes, T. 2003, Opt. Commun., 227, 75 [NASA ADS] [CrossRef] [Google Scholar]
  76. Zaninetti, L. 2012, Rev. Mex. Astron. Astrofis., 48, 209 [NASA ADS] [Google Scholar]
  77. Zhao, H. 1996, MNRAS, 278, 488 [Google Scholar]
  78. Zhao, H. 1997, MNRAS, 287, 525 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: The function 𝒮λ

We define the function 𝒮λ(x), with λ ≥ 0 and 0 ≤ x ≤ 1 as

(A.1)

For integer values of λ, this integral is readily evaluated in terms of elementary functions, for example

(A.2)

(A.3)

(A.4)

For non-integer values of λ this integration is more cumbersome. A general expression can be given in terms of the hypergeometric function,

(A.5)

The asymptotic expansion for x → 0 is

(A.6)

whereas for x → 1 we find

(A.7)

Appendix B: Nuker models with α = 2

For general values of the parameters α, β and γ, the density and potential of the Nuker models can be written in terms of the Fox H function. For the subset of models characterised by α = 2, these formulae are reduced in complexity. Indeed, when α = 2, all components of the vectors A and B are equal to one and the Fox H functions reduce to Meijer G functions. Instead of the formulae (14) and (29), we immediately obtain

(B.1)

(B.2)

The density can also be written in terms of the hypergeometric function,

(B.3)

when γ = 1, this expression for the density can be further reduced in case β is also an integer number. For odd values of β, the density can be written purely in terms of elementary functions, for even values of β it involves complete elliptic integrals.

Particularly simple cases are those models with α = 2 and γ = 0, a one-parameter family of Plummer-like models. In this case, we find a simple connection between the surface brightness profile and the density profile,

(B.4)

(B.5)

and the gravitational potential is

(B.6)

This one-parameter family of models is not only a sub-family of the general Nuker family, but also of the class of Zhao or generalised NFW models. These models are characterised by a double power-law density profile similar to the Nuker profile, but for the spatial density rather than for the surface brightness on the plane of the sky (Zhao 1996). This one-parameter family of models contains a number of well-known simple models. Setting β = 4, one can recognise the Plummer model (Plummer 1911; Dejonghe 1987),

(B.7)

(B.8)

(B.9)

For β = 3 we recover the perfect sphere (de Zeeuw 1985),

(B.10)

(B.11)

(B.12)

All Figures

thumbnail Fig. 1.

Surface brightness profile (top row) and density profile (second row) for the family of Nuker models. The three columns correspond to different values of γ and different lines within each panel correspond to different values of α. Panels on the bottom row: zoom in on the density profile around the break radius r = 1 and contain additional Nuker models with larger values of α.

Open with DEXTER
In the text
thumbnail Fig. 2.

Gravitational potential for the family of Nuker models. The selection of models and the meaning of the different lines are as in Fig. 1.

Open with DEXTER
In the text
thumbnail Fig. 3.

Intrinsic velocity dispersion profiles (top) and line-of-sight velocity dispersion profiles (bottom) for the family of isotropic Nuker models. The selection of models and the meaning of the different lines are as in Fig. 1.

Open with DEXTER
In the text
thumbnail Fig. 4.

Distribution function (top) and differential energy distribution (bottom) for the family of isotropic Nuker models. The selection of models and the meaning of the different lines are as in Fig. 1. Dotted lines correspond to unphysical models, that is, models that have a negative distribution function in some region of phase space.

Open with DEXTER
In the text
thumbnail Fig. 5.

Classification of the central structure of spherical galactic nuclei, based on the work of Tremaine et al. (1994). See text for discussion. The current figure corresponds to β = 4, but it has the classification has the same qualitative behaviour for all values of β >  2.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.