Issue 
A&A
Volume 634, February 2020



Article Number  A109  
Number of page(s)  12  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201937209  
Published online  19 February 2020 
The Nuker model for galactic nuclei
Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, 9000 Gent, Belgium
email: maarten.baes@ugent.be
Received:
29
November
2019
Accepted:
17
December
2019
The Nuker profile, characterised by an inner and outer powerlaw profile smoothly merged around a break radius, is a very popular model to describe the surface brightness profile of galactic nuclei. A disadvantage of this model for dynamical studies is that the spatial density distribution that corresponds to this surface brightness profile cannot be written in terms of elementary or regular special functions. We derive a compact and elegant analytical expression for the density of the Nuker model, based the Mellin integral transform method. We use this expression to discuss the general behaviour and asymptotic expansion of the density. We also discuss the special subclass of Nuker models with an infinitely sharp break and demonstrate that these models are always characterised by nonmonotonous and hence unphysical density profile. We extend our study to the dynamical structure of spherical isotropic galactic nuclei with a Nuker surface brightness profile. Based on this analysis, we extend and refine the classification of spherical isotropic galactic nuclei introduced by Tremaine et al. (1994, AJ, 107, 634). We demonstrate that both the inner density slope and the sharpness of the break between the inner and outer profiles critically determine the consistency and stability of the Nuker models.
Key words: methods: analytical / galaxies: nuclei / galaxies: structure / galaxies: kinematics and dynamics
© ESO 2020
1. Introduction
The first highresolution observations of galactic nuclei with the Hubble Space Telescope showed that their surface brightness profiles are generally characterised by powerlaw 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,
In this formula, R is the circular radius on the plane of the sky, R_{b} is the break radius that indicates the transition between the inner and the outer profile, I_{b} = I(R_{b}) 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 powerlaw 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 deprojection 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 powerlaw model, which is very similar to what is presented in Eq. (1). This general family of models, also called the generalised NavarroFrenkWhite (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 followup 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 leastsquares sense.
RindlerDaller 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 phasespace 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 threeparameter 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 = R_{b} = 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,
If α > 0, β > 2 and γ < 2, we find a finite luminosity
Using this expression and the convention of normalised units, we can rewrite the surface brightness profile of the Nuker model as
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 masstolightratio, the density profile ρ(r) can be found from the surface brightness profile I(R) through the standard deprojection formula (e.g. Binney & Tremaine 2008),
This gives
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 wellknown 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 wellknown. The method is surprisingly flexible and powerful, however, and it often yields closedform expressions that are very difficult to come up with using other methods. We start by rewriting expression (6) in the form
with
and
Expression (7) is a Mellin convolution of the two functions f_{1} and f_{2}. Similarly to the wellknown 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 f_{1} and f_{2}. Explicitly, we find
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 f_{1} and f_{2} can be calculated exactly,
This yields the following expression,
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 MellinBarnes integrals (Barnes 1910; Slater 1966). This particular MellinBarnes integral can be written in a compact and elegant form as
with the Fox H function, generally defined as
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 MittagLeffler 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 logarithmicpower series if there are multiple poles (Mathai et al. 2009). Explicit power and powerlogarithmic 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 doublyexponential 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 GPUenabled 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; RetanaMontenegro 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
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 nonnegative 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.
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 α. 
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 powerlawlike 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
This r^{−β − 1} powerlaw behaviour is expected, given the R^{−β} powerlaw 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
where E ≈ 0.57721566 is Euler’s constant^{1}, ψ(x) is the digamma function and
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 powerlaw falloff 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 bottomleft 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 powerlaw profile,
Apart from the “regular” Nuker models, the top panels of Fig. 1 also show the surface brightness profile corresponding to this broken powerlaw 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 MellinBarnes expression (13). After some calculation, the resulting expression is a simpler MellinBarnes integral that can be written as a Meijer G function,
Alternatively, we can also substitute expression (20) directly in the deprojection formula (5). After some algebra, we find
with the function 𝒮_{λ} defined as (see Appendix A)
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 powerlaw profile, as indicated in Eq. (22). At small radii,
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} powerlaw falloff that characterises the entire density profile for these models. For r ≲ 1, however, we obtain, thanks to Eq. (A.7),
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 nonmonotonic 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
where M(r) is the cumulative mass distribution
Starting from the MellinBarnes 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,
and
Analysing the latter expression, it is easy to demonstrate that the Nuker models are characterised by a finite potential well if γ < 1,
This result can also be obtained by plugging Eq. (4) into
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
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 nonmonotonic density profile, have a smooth potential.
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. 
3.2. Isotropic dynamical models
For each spherical potentialdensity 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,
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 Hfunction manipulation, we find
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.
Fig. 3. Intrinsic velocity dispersion profiles (top) and lineofsight 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. 
Apart from the intrinsic velocity dispersion, we can also investigate the lineofsight 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,
On the bottom row of Fig. 3 we show the lineofsight 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 lineofsight 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
A similar difference between intrinsic and lineofsight 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 lineofsight 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 lineofsight 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.
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. 
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 potentialdensity 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 selfconsistently 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. 1–2 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 nonradial 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 nonradial perturbations anymore. For γ = 0.6, the critical value of α that separates models with monotonic and nonmonotonic 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 foursteps recipe. Firstly, we reintroduce the scaling parameters R_{b} and L. Secondly, we set , with m and ϵ positive numbers that satisfy . Thirdly, we introduce a new length scale R_{e} = (ϵb)^{m}R_{b}, 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
To evaluate this expression, we use Wendel’s asymptotic relation (Wendel 1948; Qi & Luo 2012),
This results in
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 earlytype 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 foursteps 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),
Using the explicit MellinBarnes form (15) of the Fox H function, this limit can be evaluated as
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,
for m < 1 and a powerlaw density cusp,
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 foursteps recipe to Eqs. (29)–(32).
5. Discussion and summary
In this paper we construct a threeparameter family of spherical models characterised by a simple surface brightness profile and we discuss the corresponding potentialdensity 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 coreSérsic model proposed by Graham et al. (2003) have a broad applicability to model the entire surface brightness profile of earlytype 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 lineofsight 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 lineofsight 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 lineofsight 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^{−γ} powerlaw behaviour at small radii.
One limitation of the study of Tremaine et al. (1994) is that the oneparameter 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 threeparameter 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^{−γ} powerlaw 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.
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. 
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 lineofsight 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 nonzero central intrinsic and lineofsight 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 oneparameter family of (α, β, γ) = (2, β, 0) models, which we call the Plummerlike models. The potentialdensity pair of this subfamily 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 wellbehaved 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 nonradial 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.
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
 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]
 Antonov, V. A. 1962, Solution of the Problem of Stability of Stellar System (Leningrad: Vestnik Leningradskogo Universiteta) [Google Scholar]
 Baes, M. 2009, MNRAS, 392, 1503 [NASA ADS] [CrossRef] [Google Scholar]
 Baes, M., & Ciotti, L. 2019a, A&A, 626, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baes, M., & Ciotti, L. 2019b, A&A, 630, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baes, M., & Gentile, G. 2011, A&A, 525, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baes, M., & Van Hese, E. 2011, A&A, 534, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barnes, E. W. 1910, Q. J. Math., 41, 136 [Google Scholar]
 Bertin, G., Ciotti, L., & Del Principe, M. 2002, A&A, 386, 149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press) [Google Scholar]
 Byun, Y. I., Grillmair, C. J., Faber, S. M., et al. 1996, AJ, 111, 1889 [NASA ADS] [CrossRef] [Google Scholar]
 Casasent, D., & Psaltis, D. 1976, Opt. Eng., 15, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Chergui, H., Benjillali, M., & Alouini, M.S. 2019, IEEE Wirel. Commun. Lett., 8, 428 [CrossRef] [Google Scholar]
 Ciotti, L. 1991, A&A, 249, 99 [NASA ADS] [Google Scholar]
 Ciotti, L., & Lanzoni, B. 1997, A&A, 321, 724 [NASA ADS] [Google Scholar]
 Ciotti, L., & Pellegrini, S. 1992, MNRAS, 255, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Coffey, M. W., & Lettington, M. C. 2015, J. Number Theory, 148, 507 [CrossRef] [Google Scholar]
 Crane, P., Stiavelli, M., King, I. R., et al. 1993, AJ, 106, 1371 [NASA ADS] [CrossRef] [Google Scholar]
 de Ruiter, H. R., Parma, P., Capetti, A., et al. 2005, A&A, 439, 487 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Zeeuw, T. 1985, MNRAS, 216, 273 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, P. T., & Carollo, C. M. 1996, MNRAS, 281, 1333 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, T., & Franx, M. 1991, ARA&A, 29, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W. 1993, MNRAS, 265, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Dejonghe, H. 1986, Phys. Rep., 133, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Dejonghe, H. 1987, MNRAS, 224, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Dejonghe, H. 1989, ApJ, 343, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Doremus, J.P., Feix, M. R., & Baumann, G. 1971, Phys. Rev. Lett., 26, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Dullo, B. T., & Graham, A. W. 2014, MNRAS, 444, 2700 [NASA ADS] [CrossRef] [Google Scholar]
 Durret, F., Tarricq, Y., Márquez, I., Ashkar, H., & Adami, C. 2019, A&A, 622, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Epstein, B. 1948, Ann. Math. Statist., 19, 370 [CrossRef] [Google Scholar]
 Evans, N. W., & An, J. 2005, MNRAS, 360, 492 [NASA ADS] [CrossRef] [Google Scholar]
 Faber, S. M., Tremaine, S., Ajhar, E. A., et al. 1997, AJ, 114, 1771 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Ferrarese, L., Côté, P., Jordán, A., et al. 2006, ApJS, 164, 334 [NASA ADS] [CrossRef] [Google Scholar]
 Fikioris, G. J. 2007, Mellin Transform Method for Integral Evaluation: Introduction and Applications to Electromagnetics (Morgan & Claypool Publishers) [Google Scholar]
 Fox, C. 1957, Math. Proc. Camb. Philos. Soc., 53, 620 [NASA ADS] [CrossRef] [Google Scholar]
 Fox, C. 1961, Trans. Am. Math. Soc., 98, 395 [Google Scholar]
 Galassi, M., Gough, B., Rossi, F., et al. 2001, GNU Scientific Library: Reference Manual (Network Theory Limited) [Google Scholar]
 Graham, A. W., & Driver, S. P. 2005, PASA, 22, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Graham, A. W., Erwin, P., Trujillo, I., & Asensio Ramos, A. 2003, AJ, 125, 2951 [NASA ADS] [CrossRef] [Google Scholar]
 Haubold, H. J., Mathai, A. M., & Saxena, R. K. 2007, Bull. Astron. Soc. India, 35, 681 [NASA ADS] [Google Scholar]
 Hernquist, L., & Quinn, P. J. 1989, ApJ, 342, 1 [Google Scholar]
 Kilbas, A. A., & Saigo, M. 1999, J. Appl. Math. Stoch. Anal., 12, 191 [CrossRef] [Google Scholar]
 Laine, S., van der Marel, R. P., Lauer, T. R., et al. 2003, AJ, 125, 478 [Google Scholar]
 Lauer, T. R., Faber, S. M., Holtzman, J. A., et al. 1991, ApJ, 369, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Lauer, T. R., Faber, S. M., Lynds, R. C., et al. 1992, AJ, 103, 703 [NASA ADS] [CrossRef] [Google Scholar]
 Lauer, T. R., Ajhar, E. A., Byun, Y. I., et al. 1995, AJ, 110, 2622 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lauer, T. R., Faber, S. M., Gebhardt, K., et al. 2005, AJ, 129, 2138 [NASA ADS] [CrossRef] [Google Scholar]
 Lauer, T. R., Gebhardt, K., Faber, S. M., et al. 2007, ApJ, 664, 226 [NASA ADS] [CrossRef] [Google Scholar]
 Marichev, O. I. 1983, Handbook of Integral Transforms of Higher Transcendental Functions: Theory and Algorithmic Tables (Ellis Horwood Publisher) [Google Scholar]
 Mathai, A. M., Saxena, R. K., & Haubold, H. J. 2009, The HFunction: Theory and Applications (New York: Springer) [Google Scholar]
 Meijer, C. S. 1946, Proc. Nederl. Akad. Wetensch., 98, 227 [Google Scholar]
 Oliphant, T. E. 2007, Comput. Sci. Eng., 9, 10 [CrossRef] [PubMed] [Google Scholar]
 Pearson, J. W., Olver, S., & Porter, M. A. 2017, Numer. Algorithm, 74, 821 [CrossRef] [Google Scholar]
 Phillips, J. 1986, The Nag Library: A Beginner’s Guide (New York: Clarendon Press) [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [CrossRef] [Google Scholar]
 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]
 Qi, F., & Luo, Q.M. 2012, Banach J. Math. Anal., 6, 132 [CrossRef] [Google Scholar]
 Quillen, A. C., Bower, G. A., & Stritzinger, M. 2000, ApJS, 128, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Rest, A., van den Bosch, F. C., Jaffe, W., et al. 2001, AJ, 121, 2431 [NASA ADS] [CrossRef] [Google Scholar]
 RetanaMontenegro, E., van Hese, E., Gentile, G., Baes, M., & FrutosAlfaro, F. 2012a, A&A, 540, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RetanaMontenegro, E., FrutosAlfaro, F., & Baes, M. 2012b, A&A, 546, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RindlerDaller, T., Dejonghe, H., & Zeilinger, W. W. 2005, MNRAS, 356, 1403 [NASA ADS] [CrossRef] [Google Scholar]
 Saxena, R. K., Mathai, A. M., & Haubold, H. J. 2006, Ap&SS, 305, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Schäling, B. 2014, The Boost C++ Libraries (Fort Collins: XML Press) [Google Scholar]
 Schwarzschild, M. 1979, ApJ, 232, 236 [NASA ADS] [CrossRef] [Google Scholar]
 Schwarzschild, M. 1993, ApJ, 409, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Sérsic, J. L. 1968, Atlas de Galaxias Australes (Cordoba, Argentina: Observatorio Astronomico) [Google Scholar]
 Shafique Ansari, I., Yilmaz, F., Alouini, M.S., & Kucur, O. 2014, Trans. Emerg. Telecommun. Technol., 28, 1 [Google Scholar]
 Slater, L. J. 1966, Generalized Hypergeometric Functions (Cambridge: Cambridge University Press) [Google Scholar]
 Tremaine, S., Richstone, D. O., Byun, Y.I., et al. 1994, AJ, 107, 634 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo, I., Erwin, P., Asensio Ramos, A., & Graham, A. W. 2004, AJ, 127, 1917 [NASA ADS] [CrossRef] [Google Scholar]
 Van Hese, E., Baes, M., & Dejonghe, H. 2009, ApJ, 690, 1280 [NASA ADS] [CrossRef] [Google Scholar]
 Wendel, J. G. 1948, Am. Math. Mon., 55, 563 [CrossRef] [MathSciNet] [Google Scholar]
 Wu, C., Somervell, A., Haskell, T., & Barnes, T. 2003, Opt. Commun., 227, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Zaninetti, L. 2012, Rev. Mex. Astron. Astrofis., 48, 209 [NASA ADS] [Google Scholar]
 Zhao, H. 1996, MNRAS, 278, 488 [Google Scholar]
 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
For integer values of λ, this integral is readily evaluated in terms of elementary functions, for example
For noninteger values of λ this integration is more cumbersome. A general expression can be given in terms of the hypergeometric function,
The asymptotic expansion for x → 0 is
whereas for x → 1 we find
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
The density can also be written in terms of the hypergeometric function,
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 oneparameter family of Plummerlike models. In this case, we find a simple connection between the surface brightness profile and the density profile,
and the gravitational potential is
This oneparameter family of models is not only a subfamily of the general Nuker family, but also of the class of Zhao or generalised NFW models. These models are characterised by a double powerlaw 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 oneparameter family of models contains a number of wellknown simple models. Setting β = 4, one can recognise the Plummer model (Plummer 1911; Dejonghe 1987),
For β = 3 we recover the perfect sphere (de Zeeuw 1985),
All Figures
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 α. 

In the text 
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. 

In the text 
Fig. 3. Intrinsic velocity dispersion profiles (top) and lineofsight 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. 

In the text 
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. 

In the text 
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. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.