EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 525 (January 2011) A&A, 525 (2011) A136 Full HTML
Free Access
Issue
A&A
Volume 525, January 2011
Article Number A136
Number of page(s) 7
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201015716
Published online 08 December 2010

© ESO, 2010

1. Introduction

The Sérsic model (Sérsic 1968) is a three-parameter model for the surface brightness profile of galaxies that has been introduced as a generalization of the de Vaucouleurs R1/4 model (de Vaucouleurs 1948). It is defined as I(R)=I0exp[−b(RRe)1/m]\begin{equation} I(R) = I_0 \exp \left[-b\left(\frac{R}{R_{{\rm e}}}\right)^{1/m}\right] \label{SersicI} \end{equation}(1)where I(R) is the surface brightness at radius (on the plane of the sky) R, I0 the central surface brightness, Re the effective radius, and m the so-called Sérsic index that describes the index of the logarithmic slope power law. The parameter b is a dimensionless parameter that depends on the Sérsic index m, and its value can be derived from the requirement that the isophote corresponding to Re encloses half of the total flux.

Over the past two decades, the Sérsic model has become the standard for describing the surface brightness profiles of early-type galaxies and bulges of spiral galaxies (e.g. Davies et al. 1988; Caon et al. 1993; D’Onofrio et al. 1994; Cellone et al. 1994; Andredakis et al. 1995; Prugniel & Simien 1997; Möllenhoff & Heidt 2001; Graham & Guzmán 2003; Allen et al. 2006; Méndez-Abreu et al. 2008; Gadotti 2009). In the past few years, Sérsic-like models have also gained popularity as a model to describe the spherically averaged profiles for dark matter haloes. While models with a power-law behaviour at small and large radii were preferred in earlier simulations (e.g. Navarro et al. 1997; Moore et al. 1999), other models seem to fit the mass density distribution of more recent (and higher resolution) simulations better (Navarro et al. 2004, 2010; Merritt et al. 2005, 2006; Graham et al. 2006; Aceves et al. 2006; Gao et al. 2008; Duffy et al. 2008; Stadel et al. 2009). Among these models, the Sérsic model (i.e., a model where the projected surface density is described by a Sérsic law) has also been proposed as a universal description for simulated dark matter haloes.

Mainly as a result of its popularity for describing the surface brightness profiles of early-type galaxies, the properties of the Sérsic model have been examined in great detail. Ciotti (1991) and Ciotti & Lanzoni (1997) give a detailed description of the properties of the Sérsic model, including spatial and dynamical properties. Ciotti & Bertin (1999) provide a full asymptotic expansion of the dimensionless scale factor b of the Sérsic model. Graham & Driver (2005) present a compendium of mathematical formulae on photometric parameters such as Kron magnitues and Petrosian indices, and Cardone (2004) and Elíasdóttir & Möller (2007) investigate the lensing properties. A major problem with the Sérsic models is that the deprojection of the surface brightness profile to a luminosity density is in general non-analytical. In practice, one often uses approximations for the Sérsic models when the luminosity density (or the mass density when the Sérsic model is used to describe the distribution of dark matter) is necessary (e.g. Prugniel & Simien 1997; Lima Neto et al. 1999; Trujillo et al. 2002). An unexpected analytical progress was the work by Mazure & Capelato (2002), who demonstrate that it is possible to elegantly write the spatial luminosity density of the Sérsic model in terms of the Meijer G function. Unfortunately, their result only holds for integer values of the Sérsic index, which is a significant limitation for practical applications. Moreover, since their result fell as a deus ex machina out of the MathematicaR◯\hbox{$\tiny^{\textregistered}$} computer algebra package, it is hard to see whether it can be generalized to all Sérsic indices.

In this paper, we tackle the deprojection of the Sérsic surface brightness profile using analytical means. We apply an integration method based on Mellin integral transforms and derive an analytical expression for the luminosity density in terms of the Fox H function for all values of the Sérsic index m. For rational values of m, the luminosity density can be written in terms of the Meijer G function. As the Meijer G function is nowadays available both in symbolic computer algebra packages and as high-performance computing code, this opens up the possibility of calculating the luminosity density of the Sérsic models to arbitrary precision. The wide range of analytical properties of the Meijer G function also allow easy study of the asymptotic behaviour at small and large radii and computation of derivative quantities, such as the cumulative luminosity and the gravitational potential.

This paper is organized as follows. In Sect. 2 we discuss the general deprojection of the Sérsic surface brightness profile using the Mellin transform method and present a general expression in terms of the Fox H function. In Sect. 3 we present simpler expressions for integer, half-integer, and rational values of m in terms of the Meijer G function and discuss two special, interesting cases. In Sect. 4 we use the expressions for the luminosity density to calculate the total luminosity of the Sérsic models, which serves as a consistency check on the derived formulae. In Sect. 5 we investigate the asymptotic behaviour of the luminosity density, and in Sect. 6 we derive analytical expressions for the cumulative luminosity and the gravitational potential. Finally, in Appendix A we introduce the Meijer G and Fox H functions and discuss some of their most useful properties.

2. The luminosity density of the Sérsic model as a Fox H function

In spherical symmetry, the deprojected luminosity density ν(r) at the spatial radius r can be found from the surface brightness profile I(R) through the standard deprojection formula, ν(r)=−1π∫r∞dIdR(R) dRR2−r2·\begin{equation} \nu(r) = -\frac{1}{\pi} \int_r^\infty \frac{{{\rm d}}I}{{{\rm d}}R}(R)\, \frac{{{\rm d}}R}{\sqrt{R^2-r^2}}\cdot \label{deprojection} \end{equation}(2)If we introduce the reduced radial coordinate s=rRe,\begin{equation} s=\frac{r}{R_{\rm e}}, \end{equation}(3)we find ν(r)=I0Re bm π ∫s∞e−bt1/mt1m−1 dtt2−s2·\begin{equation} \nu(r) = \frac{I_0}{R_{\rm e}}\, \frac{b}{m\,\pi}\, \int_s^\infty \frac{{\rm e}^{-bt^{1/m}} t^{\frac1m-1}\,{{\rm d}}t} {\sqrt{t^2-s^2}}\cdot \label{depint} \end{equation}(4)The integral (4) cannot be evaluated in terms of elementary functions or even the standard special functions for general values of m. To evaluate it, we use a general method that builds on Mellin integral transforms and has become known as the Mellin transform method (Marichev 19839; Adamchik 1996; Fikioris 2007). The Mellin transform Mf(u) of a function f(z) is defined as Mf(u)=φ(u)=∫0∞f(z) zu−1 dz.\begin{equation} {\mathfrak{M}}_f(u) = \phi(u) = \int_0^\infty f(z)\,z^{u-1}\,{{\rm d}}z. \end{equation}(5)The inverse Mellin transform is found as Mφ-1(z)=f(z)=12πi∫ℒφ(u) z−u du,\begin{equation} {\mathfrak{M}}^{-1}_\phi(z) = f(z) = \frac{1}{2\pi {\rm i}} \int_{\mathcal{L}} \phi(u)\,z^{-u}\,{{\rm d}}u, \end{equation}(6)where the ℒ is a line integral over a vertical line in the complex plane. The driving force behind the Mellin transform method is the Mellin convolution theorem. The Mellin convolution of two functions f1(z) and f2(z) is defined as (f1⋆f2)(z)=∫0∞f1(t) f2(zt) dtt·\begin{equation} (f_1\star f_2)(z) = \int_0^\infty f_1(t)\,f_2\left(\frac{z}{t}\right)\, \frac{{{\rm d}}t}{t}\cdot \end{equation}(7)Similar to the well-known Fourier transform analogue, the Mellin convolution theorem states that the Mellin transform of a Mellin convolution is equal to the products of the Mellin transforms of the original functions, Mf1⋆f2(u)=Mf1(u) Mf2(u).\begin{equation} {\mathfrak{M}}_{f_1\star f_2}(u) = {\mathfrak{M}}_{f_1}(u)\,{\mathfrak{M}}_{f_2}(u). \end{equation}(8)Now it can be shown that any definite integral f(z)=∫0∞g(t,z) dt\begin{equation} f(z) = \int_0^\infty g(t,z)\,{{\rm d}}t \label{Marichev} \end{equation}(9)can be written as the Mellin convolution of two functions f1 and f2. As a result, the definite integral (9) can be transformed into an inverse Mellin transform, f(z)=12πi∫ℒMf1(u) Mf2(u) z−u du.\begin{equation} f(z) = \frac{1}{2\pi {\rm i}} \int_{\mathcal{L}} {\mathfrak{M}}_{f_1}(u)\,{\mathfrak{M}}_{f_2}(u)\, z^{-u}\,{{\rm d}}u. \label{MellinBarnes} \end{equation}(10)The power of this approach is that, if the functions f1 and f2 are of the hypergeometric type, which is true for many elementary functions and the majority of special functions, the integral (10) turns out to be a Mellin-Barnes integral. Depending on the involved coefficients, this integral can be evaluated as a Fox H function or, in simpler cases, as a Meijer G function (see Appendix A).

The form of Eq. (4) allows the Mellin tranform method to be applied, with z = 1 and f1(t)=I0Rebm π e−bt1/mt1m\begin{equation} f_1(t) = \frac{I_0}{R_{\rm e}} \frac{b}{m\,\pi}\, {\rm e}^{-bt^{1/m}} t^{\frac{1}{m}} \end{equation}(11)and

f2(t)={\begin{eqnarray} f_2(t) = \left\{ \begin{array}{lc} \frac{t}{\sqrt{1-s^{2}t^{2}}} &\qquad {\rm if }\,\, 0\leq t \leq s^{-1} \\ 0 & {\rm else}.\end{array} \right. \end{eqnarray}(12)The Mellin transforms of these functions are readily calculated Mf1(u)=I0Re 1π Γ(1+m u)bmuMf2(u)=π Γ(1+u2)Γ(u2) 1u s1+u·\begin{eqnarray} {\mathfrak{M}}_{f_1}(u) &=& \frac{I_0}{R_{\rm e}}\, \frac{1}{\pi}\, \frac{\Gamma(1+m\,u)}{b^{mu}} \\ {\mathfrak{M}}_{f_2}(u) &=& \frac{\sqrt{\pi}\,\Gamma\left(\frac{1+u}{2}\right)} {\Gamma\left(\frac{u}{2}\right)}\,\frac{1}{u\,s^{1+u}}\cdot \end{eqnarray}Substituting these values in the integral (10) and setting u = 2x, we obtain ν(r)=I0Re 1π s-112πi∫ℒΓ(1+2m x) Γ(12+x)Γ(1+x) (b2ms2)−x dx.\begin{equation} \nu(r) \! =\! \frac{I_0}{R_{\rm e}}\, \frac{1}{\sqrt{\pi}}\,s^{-1}\! \frac{1}{2\pi {\rm i}} \int_{\mathcal{L}} \frac{\Gamma(1+2m\,x)\,\Gamma\left(\frac12+x\right)} {\Gamma(1+x)}\, (b^{2m}s^2)^{-x}\, {{\rm d}}x. \label{nuHint} \end{equation}(15)If we compare this expression with the definition (A.11) of the Fox H function, we see that we can write the luminosity density of the Sérsic models in compact form as ν(r)=I0Re 1π s-1 H1,22,0[(1,1)(1,2m),(12,1) | b2ms2].\begin{equation} \nu(r) = \frac{I_0}{R_{\rm e}}\, \frac{1}{\sqrt{\pi}}\,s^{-1}\, H^{2,0}_{1,2}\!\left[\left. \begin{matrix} (1,1)\\ (1,2m), (\frac{1}{2},1) \end{matrix} \,\right|\, b^{2m}s^2 \right]. \label{nuH} \end{equation}(16)

3. Integer and rational values of the Sérsic index

Expression (16) represents a closed analytical expression for the luminosity density of the Sérsic model in terms of the Fox H function. While this function is gradually receiving more attention both in mathematics and applied sciences, ranging from astrophysics and earth sciences to statistics, its practical usefulness is still limited. In particular, no general numerical implementations of the Fox H function are available (to our knowledge). Fortunately, the Fox H function can be reduced to a Meijer G function in many cases. In this section we derive an analytical expression for the luminosity density of the Sérsic model in terms of the Meijer G function for all rational values of m. The Meijer G function is much better documented and can be an extremely useful tool for analytical work. The Meijer G function has many general properties that can manipulate expressions to equivalent forms and reduce the order for certain values of the parameters, among other advantages. Some of the most useful properties are listed in Appendix A, but there are many more. A particularly rich online source of information is the Wolfram Functions Site1. Moreover, commercial computer algebra systems such as MapleR◯\hbox{$\tiny^{\textregistered}$} and MathematicaR◯\hbox{$\tiny^{\textregistered}$} contain an implementation of the Meijer G function. It is also implemented in the open-source computer algebra package Sage and a freely available Python implementation of the Meijer G function to arbitrary precision is available from the Mpmath library2.

3.1. Integer and half-integer values of m

If m is an integer or half-integer value, we can simplify expression (16) by using the property 􏽙j=0N−1Γ(j+zN)=N12−z (2π)N−12 Γ(z).\begin{equation} \prod_{j=0}^{N-1} \Gamma\!\left(\frac{j+z}{N}\right) = N^{\frac12-z}\,(2\pi)^{\frac{N-1}{2}}\,\Gamma(z). \label{prodGamma} \end{equation}(17)Applying this recipe with N = 2m and z = 2m   x gives Γ(1+2m x)=(2m)12+2mx (2π)12−m ×  Γ(1+x) 􏽙j=12m−1Γ(j2m+x).\begin{eqnarray} \Gamma(1+2m\,x) &=& (2m)^{\frac12+2mx}\,(2\pi)^{\frac12-m}\, \nonumber\\ &&\times\,\, \Gamma(1+x)\,\prod_{j=1}^{2m-1}\Gamma\!\left(\frac{j}{2m}+x\right). \end{eqnarray}(18)By inserting this in expression (15), we find ν(r)=2 I0Re m(2π)m s-1 12πi×  ∫ℒΓ(12+x)􏽙j=12m−1Γ(j2m+x)[(b2m)2ms2]−x dx.\begin{eqnarray} \nu(r) &=& \frac{2\,I_0}{R_{\rm e}}\, \frac{\sqrt{m}}{(2\pi)^m}\,s^{-1}\, \frac{1}{2\pi {\rm i}} \nonumber\\ &&\times\,\, \int_{\mathcal{L}} \Gamma\!\left(\frac12+x\right) \prod_{j=1}^{2m-1}\Gamma\!\left(\frac{j}{2m}+x\right) \left[ \left(\frac{b}{2m}\right)^{2m}s^2 \right]^{-x}\, {{\rm d}}x. \end{eqnarray}(19)When comparing this expression with the definition (A.1) of the Meijer G function, we obtain the following compact expression for the luminosity density of the Sérsic model ν(r)=2 I0Re m(2π)m s-1 G0,2m2m,0[−b| (b2m)2ms2]% subequation 1164 0 \begin{equation} \nu(r) = \frac{2\,I_0}{R_{{\rm e}}}\, \frac{\sqrt{m}}{(2\pi)^m}\, s^{-1}\, G_{0,2m}^{2m,0} \left[ \begin{array}{l} - \\ {\vec{b}} \end{array} \left|\, \left(\frac{b}{2m}\right)^{2m} s^2 \right. \right] \end{equation}(20a)with b a vector with 2m elements given by b={12m,22m...,2m−12m,12}·% subequation 1164 1 \begin{equation} {\vec{b}} = \left\{\frac{1}{2m},\frac{2}{2m}\ldots, \frac{2m-1}{2m}, \frac12\right\}\cdot \end{equation}(20b)This expression is equivalent to (and actually even slightly simpler than) the expression obtained by Mazure & Capelato (2002) using the computer algebra package MathematicaR◯\hbox{$\tiny^{\textregistered}$}. Mazure & Capelato (2002) obtained their formula only for integer values of m, whereas our analysis shows that exactly the same expression also holds for half-integer values of m.

3.2. Rational values of m

Interestingly, these results can be generalized for all rational values of m. Setting m = p/q with p and q integer numbers, we can write expression (15) as ν(r)=I0Re 1π s-1 ×  12πi∫ℒq Γ(1+2p y) Γ(12+q y)Γ(1+q y) (b2ps2q)−y dy.\begin{eqnarray} \nu(r) &=& \frac{I_0}{R_{\rm e}}\, \frac{1}{\sqrt{\pi}}\,s^{-1}\, \nonumber\\ &&\times\,\, \frac{1}{2\pi {\rm i}} \int_{\mathcal{L}} \frac{q\,\Gamma(1+2p\,y)\,\Gamma\left(\frac12+q\,y\right)} {\Gamma(1+q\,y)}\, (b^{2p}s^{2q})^{-y}\, {{\rm d}}y. \end{eqnarray}(21)By multiple application of the identity (17), we can rewrite this expression in a format that leads to a Meijer G function. The result is ν(r)=2 I0Re p q(2π)p s-1 Gq−1,2p+q−12p+q−1,0[ab| (b2p)2ps2q]% subequation 1283 0 \begin{equation} \nu(r) = \frac{2\,I_0}{R_{\rm e}}\, \frac{\sqrt{p\,q}}{(2\pi)^p}\, s^{-1}\, G_{q-1,2p+q-1}^{2p+q-1,0} \left[ \begin{array}{l} {{\vec a}} \\ {{\vec b}} \end{array} \left|\, \left(\frac{b}{2p}\right)^{2p} s^{2q} \right. \right] \end{equation}(22a)with a a vector of dimension q − 1 with elements a={1q,2q,...,q−1q},% subequation 1283 1 \begin{equation} {{\vec a}} = \left\{ \frac{1}{q},\frac2q,\ldots,\frac{q-1}{q} \right\}, \end{equation}(22b)and b a vector with 2p + q − 1 elements given by b={12p,22p,...,2p−12p,12q,32q,...,2q−12q}·% subequation 1283 2 \begin{equation} {{\vec b}} = \left\{ \frac{1}{2p},\frac{2}{2p},\ldots,\frac{2p-1}{2p}, \frac{1}{2q},\frac{3}{2q},\ldots,\frac{2q-1}{2q} \right\}\cdot \end{equation}(22c)It is straightforward to check that the general expression (22) reduces to expression (20) for integer values of m. Using the order-reduction formulae (A.5) and (A.6) of the Meijer G function, one can also demonstrate that expression (22) reduces to (20) for half-integer values of m.

3.3. Special cases: m = 1 and m=12\hbox{$\textit{m} =\mathsfsl{ \frac12}$}

Among the family of Sérsic models, there are two well-known specific cases for which the luminosity density can be calculated explicitly in terms of elementary or special functions. The first of these two models is the exponential model, corresponding to m = 1. Exponential models are often used to describe the surface brightness profiles of dwarf elliptical galaxies (e.g. Faber & Lin 1983; Binggeli et al. 1984). If we introduce the notation h = Re/b, we can write the surface brightness profile as I(R)=I0 e−R/h.\begin{equation} I(R) = I_0\,{\rm e}^{-R/h}. \end{equation}(23)If we deproject this surface brightness profile using the deprojection formula (2), we recover the well-known result that the luminosity density can be written in terms of a modified Bessel function of the second kind, ν(r)=I0π h K0(rh)·\begin{equation} \nu(r) = \frac{I_0}{\pi\,h}\,K_0\left(\frac{r}{h}\right)\cdot \label{nu1} \end{equation}(24)If we set m = 1 in the expression (20), we obtain ν(r)=I0π Re s-1 G0,22,0[−12,12| b2s24]·\begin{equation} \nu(r) = \frac{I_0}{\pi\,R_{\rm e}}\, s^{-1}\, G_{0,2}^{2,0} \left[ \begin{array}{l} - \\ \frac12,\frac12 \end{array} \left|\, \frac{b^2s^2}{4} \right. \right]\cdot \end{equation}(25)Using formula (A.3), this expression reduces to expression (24).

Another interesting special case is m=12\hbox{$m=\frac12$}, which corresponds to a Gaussian surface brightness profile. Such profiles do not correspond to the observed surface brightness profiles of galaxies, but they are very useful as components in multi-Gaussian expansions: even with a relatively modest set of Gaussian components, realistic geometries can accurately be reproduced (e.g., Emsellem et al. 1994a,b; De Bruyne et al. 2001; Cappellari 2002). If we introduce σ=Re/2b\hbox{$\sigma=R_{\rm e}/\sqrt{2b}$} and we use the total luminosity instead of the effective intensity as a parameter, we can write the surface brightness profile as I(R)=L2πσ2exp(−R22σ2)·\begin{equation} I(R) = \frac{L}{2\pi\sigma^2} \exp\left(-\frac{R^2}{2\sigma^2}\right)\cdot \label{Igauss} \end{equation}(26)One of the key advantages of a multi-Gaussian expansion of an observed surface brightness profile is that the corresponding luminosity density has a simple analytical form. Indeed, substituting (26) into the deprojection formula (2), one can easily check that the deprojection of a Gaussian distribution on the sky is also a Gaussian distribution, ν(r)=L(2πσ2)3/2exp(−r22σ2)·\begin{equation} \nu(r) = \frac{L}{(2\pi\sigma^2)^{3/2}} \exp\left(-\frac{r^2}{2\sigma^2}\right)\cdot \label{nu12} \end{equation}(27)This result can also be found by setting m=12\hbox{$m=\frac12$} in Eq. (20), ν(r)=I0π Re s-1 G0,11,0[−12| bs2].\begin{equation} \nu(r) = \frac{I_0}{\sqrt{\pi}\,R_{\rm e}}\, s^{-1}\, G_{0,1}^{1,0} \left[ \begin{matrix} - \\ \frac12 \end{matrix} \left|\, bs^2 \right. \right]. \end{equation}(28)If we use Eq. (A.2), we easily recover expression (27).

4. The total luminosity

The total luminosity of the Sérsic model can be calculated by integrating the intensity on the plane of the sky, L=2π∫0∞I(R) R dR.\begin{equation} L = 2\pi\int_0^\infty I(R)\,R\,{{\rm d}}R. \end{equation}(29)By inserting Eq. (1) one readily finds L=π I0 Re2 1b2m Γ(2m+1).\begin{equation} L = \pi\,I_0\,R_{\rm e}^2\, \frac{1}{b^{2m}}\, \Gamma(2m+1). \label{L} \end{equation}(30)As a consistency check on the formula (22) and to illustrate of the power of the Meijer G function as an analytical tool, we can also calculate the luminosity by integrating the luminosity density ν(r) over the entire space, L=4π∫0∞ν(r) r2 dr.\begin{equation} L = 4\pi\int_0^\infty \nu(r)\,r^2\,{{\rm d}}r. \end{equation}(31)Inserting Eq. (22), we obtain L=2 I0 Re2(2π)p−1 pq∫0∞Gq−1,2p+q−12p+q−1,0[ab| (b2p)2pt]t1q−1 dt.\begin{equation} L = \frac{2\,I_0\,R_{\rm e}^2}{(2\pi)^{p-1}}\, \sqrt{\frac{p}{q}} \int_0^\infty G_{q-1,2p+q-1}^{2p+q-1,0} \left[ \begin{array}{l} {{\vec a}} \\ {{\vec b}} \end{array} \left|\, \left(\frac{b}{2p}\right)^{2p} t \right. \right] t^{\frac1q-1}\,{{\rm d}}t. \end{equation}(32)If we use the general property (A.8) of the Meijer G function, we can evaluate this integral as L=2 I0 Re2(2π)p−1 pq(2pb)2pq􏽑j=12p+q−1Γ(1q+bj)􏽑j=1q−1Γ(1q+aj)·\begin{equation} L = \frac{2\,I_0\,R_{\rm e}^2}{(2\pi)^{p-1}}\, \sqrt{\frac{p}{q}} \left(\frac{2p}{b}\right)^{\frac{2p}{q}} \frac{ \prod_{j=1}^{2p+q-1} \Gamma\!\left(\frac{1}{q}+b_j\right) }{ \prod_{j=1}^{q-1} \Gamma\!\left(\frac{1}{q}+a_j\right) }\cdot \label{Lpro} \end{equation}(33)The product in the denominator can be simplified to 􏽙j=1q−1Γ(1q+aj)=􏽑j=0q−1Γ(j+1q)Γ(1q)=q−12(2π)q−12Γ(1q),\begin{eqnarray} \prod_{j=1}^{q-1} \Gamma\!\left(\frac{1}{q}+a_j\right) = \frac{ \prod_{j=0}^{q-1}\Gamma\!\left(\frac{j+1}{q}\right) }{ \Gamma\!\left(\frac{1}{q}\right) } = \frac{q^{-\frac12} (2\pi)^{\frac{q-1}{2}}}{\Gamma\!\left(\frac{1}{q}\right)}, \label{onder} \end{eqnarray}(34)where the last transition follows from identity (17). Similarly, the product in the numerator of Eq. (33) can be simplified to 􏽙j=12p+q−1Γ(1q+bj)=􏽑j=02p−1Γ(j+2p/q2p)􏽑j=0q−1Γ(j+3/2q)Γ(1q)=π (2π)p+q2−1(2p)12−2pq2 q Γ(2pq)Γ(1q)·\begin{eqnarray} \prod_{j=1}^{2p+q-1} \Gamma\!\left(\frac{1}{q}+b_j\right) &=& \frac{ \prod_{j=0}^{2p-1}\Gamma\!\left(\frac{j+2p/q}{2p}\right) \prod_{j=0}^{q-1}\Gamma\!\left(\frac{j+3/2}{q}\right) } {\Gamma\!\left(\frac1q\right)} \nonumber \\ \label{boven} &=& \frac{\sqrt{\pi}\,(2\pi)^{p+\frac{q}{2}-1}(2p)^{\frac12-\frac{2p}{q}}} {2\,q}\, \frac{\Gamma\!\left(\frac{2p}{q}\right)}{\Gamma\!\left(\frac1q\right)}\cdot \end{eqnarray}(35)If we substitute (34) and (35) into expression (33) and use m = p/q, we recover the expression (30) for the total luminosity of the Sérsic model, as required.

5. Asymptotic behaviour

One of the most useful properties of the general expression (22) is that it can elegantly determine the asymptotic behaviour of the luminosity density of the Sérsic model at small and large radii. It is well-known (e.g. Ciotti 1991) that the Sérsic models have a cusp for m > 1 and a finite luminosity density core at m < 1. This can be seen immediately by evaluating the integral (4) for r = 0 (which converges only for m < 1), ν0=I0Re bm π ∫0∞e−bt1/mt1m−2 dt=I0πRe bm Γ(1−m)·\begin{equation} \nu_0 = \frac{I_0}{R_{\rm e}}\, \frac{b}{m\,\pi}\, \int_0^\infty {\rm e}^{-bt^{1/m}} t^{\frac1m-2}\,{{\rm d}}t = \frac{I_0}{\pi R_{\rm e}}\,b^m\,\Gamma(1-m)\cdot \label{nu0} \end{equation}(36)For a more detailed discussion on the behaviour of the luminosity density at small radii, we can use the asymptotic expansion (A.9) of the Meijer G function. In particular, this equation shows that the lowest order term in the expansion is determined by the smallest components bk in the vector b. This depends on the value of the Sérsic index m.

For m<13\hbox{$m<\frac13$}, the smallest component of the vector b is b2p=12q\hbox{$b_{2p}=\frac{1}{2q}$} and the second-smallest is b2p+1=32q\hbox{$b_{2p+1}=\frac{3}{2q}$}. After some algebra, which involves similar techniques as applied in Sect. 4, we find the asymptotic expansion ν(r)~I0 bπRe[Γ(1−m)b1−m+12 Γ(1−3m)b1−3m  s2].% subequation 1805 0 \begin{equation} \nu(r) \sim \frac{I_0\,b}{\pi R_{\rm e}} \left[ \frac{\Gamma(1-m)}{b^{1-m}} + \frac12\,\frac{\Gamma(1-3m)}{b^{1-3m}}\,\,s^2 \right]. \label{nuasymp0-13} \end{equation}(37a)For m=13\hbox{$m=\frac13$}, the smallest component is still b2p≡b2=16\hbox{$b_{2p}\equiv b_2=\frac{1}{6}$}, but now the two components b1 and b2p + 1 ≡ b3 are both equal to 12\hbox{$\frac{1}{2}$}. In this case we cannot use the expansion formula (A.9), since this formula is only valid if all components of the vector b are different. For m=13\hbox{$m=\frac{1}{3}$} we find the expansion ν(r)~I0 bπRe[Γ(23)b2/3+32ln(1s)s2−14 (3+2γ+2lnb−6ln2) s2]% subequation 1805 1 \begin{eqnarray} \nu(r) &\sim& \frac{I_0\,b}{\pi R_{\rm e}} \Bigg[ \frac{\Gamma(\frac23)}{b^{2/3}} + \frac{3}{2}\ln\!\left(\frac{1}{s}\right)s^2 \nonumber\\ \label{nuasymp13}&&- \frac14\,(3+2\gamma+2\ln b-6\ln 2)\,s^2 \Bigg] \end{eqnarray}(37b)where γ ≈ 0.57721566 is the Euler-Mascheroni constant. If 13<m<1\hbox{$\frac{1}{3}<m<1$}, b2p=12q\hbox{$b_{2p}=\frac{1}{2q}$} remains the smallest component of the vector b, but the second-smallest is now b1=12p\hbox{$b_1=\frac{1}{2p}$}. One finds ν(r)~I0 bπRe[Γ(1−m)b1−m+π2m Γ(12−12m)Γ(1−12m) s1m−1].% subequation 1805 2 \begin{equation} \nu(r) \sim \frac{I_0\,b}{\pi R_{\rm e}} \left[ \frac{\Gamma(1-m)}{b^{1-m}} + \frac{\sqrt{\pi}}{2m}\,\frac{\Gamma\left(\frac12-\frac1{2m}\right)} {\Gamma\left(1-\frac{1}{2m}\right)}\, s^{\frac1m-1} \right]. \label{nuasymp13-1} \end{equation}(37c)For m = 1, we again have two equal components in the vector b, and we cannot readily apply formula (A.9). The asymptotic expansion for small r now reads as ν(r)~I0 bπRe[ln(1s)−(γ+lnb−ln2)].% subequation 1805 3 \begin{equation} \nu(r) \sim \frac{I_0\,b}{\pi R_{\rm e}} \left[ \ln\!\left(\frac{1}{s}\right) - (\gamma+\ln b-\ln2) \right]. \end{equation}(37d)Finally, if m > 1, the smallest component is b1=12p\hbox{$b_{1}=\frac{1}{2p}$}, which leads to ν ( r ) ~ I 0   b R e   1 2 m π   Γ ( 1 2 − 1 2 m ) Γ ( 1 − 1 2 m )   s − ( 1 − 1 m ) . % subequation 1805 4 \begin{equation} \nu(r) \sim \frac{I_0\,b}{R_{\rm e}}\, \frac{1}{2m\sqrt{\pi}}\,\frac{\Gamma\left(\frac12-\frac1{2m}\right)} {\Gamma\left(1-\frac{1}{2m}\right)}\, s^{-(1-\frac1m)}. \end{equation}(37e)The five different asymptotic expansions (37) demonstrate the different behaviours of the luminosity density at small radii, depending on the value of the Sérsic index m. For m < 1 the Sérsic model has a finite density core with the central luminosity density given by Eq. (36). At m = 1 the model has a logarithmic cusp, and at m > 1 we have a power-law cusp with logarithmic slope 1−1m\hbox{$1-\frac{1}{m}$}. In particular, the de Vaucouleurs model has a luminosity density profile that behaves as s−3/4 at small radii (Young 1976; Mellier & Mathez 1987). Surprisingly, the Sérsic models with m<12\hbox{$m<\frac12$} do not have a monotonically decreasing luminosity density profile with increasing radius. In the expansions (37a) and (37b), the first non-constant term has a positive coefficient, and hence the luminosity density increases with increasing radius in the nuclear region. The same accounts for 13<m<12\hbox{$\frac13<m<\frac12$}, since the coefficient of the second term in the expansion (37c) is positive for 13<m<12\hbox{$\frac13<m<\frac12$} and negative for 12<m<1\hbox{$\frac12<m<1$}.

At large radii, a single formula for the asymptotic expansion holds for all Sérsic indices, ν(r)~I0Reb2πm e−bs1/m (1s)1−12m,\begin{equation} \nu(r) \sim \frac{I_0}{R_{\rm e}} \sqrt{\frac{b}{2\pi m}}\, {\rm e}^{-bs^{1/m}}\, \left(\frac{1}{s}\right)^{1-\frac{1}{2m}}, \end{equation}(38)in agreement with the result derived by Ciotti (1991).

6. Some other properties of the Sérsic model

The analytical expression (22) for the luminosity density of the Sérsic models allows other properties of this family to be expressed analytically in terms of the Meijer G function. The most important ones are the cumulative luminosity and the gravitational potential.

For a spherically symmetric system, they cumulative luminosity L(r) can be calculated as L(r)=4π∫0rρ(r′) r′2 dr′.\begin{equation} L(r) = 4\pi \int_0^r \rho(r')\,r'^2\,{{\rm d}}r'. \label{defLr} \end{equation}(39)After substitution of expression (22) in Eq. (A.7), we find L(r)=2 I0 Re2(2π)p−1 pq s2 Gq,2p+q2p+q−1,1[1−1q,ab,−1q |(b2p)2ps2q],\begin{equation} L(r) = \frac{2\,I_0\,R_{\rm e}^2}{(2\pi)^{p-1}}\, \sqrt{\frac{p}{q}}\, s^2\, G^{2p+q-1,1}_{q,2p+q} \!\left[\left. \begin{array}{l} 1-\frac{1}{q}, {{\vec a}} \\ {{\vec b}}, -\frac{1}{q} \end{array} \,\right|\left(\frac{b}{2p}\right)^{2p} s^{2q} \right], \end{equation}(40)which reduces to L(r)=2 I0 Re2 m(2π)m−1 s2 G1,2m+12m,1[0b,−1 |(b2m)2ms2]\begin{equation} L(r) = \frac{2\,I_0\,R_{\rm e}^2\,\sqrt{m}}{(2\pi)^{m-1}}\, s^2\, G^{2m,1}_{1,2m+1} \!\left[\left. \begin{array}{l} 0 \\ {{\vec b}}, -1 \end{array} \,\right|\left(\frac{b}{2m}\right)^{2m} s^2 \right] \end{equation}(41)for integer or half-integer values of the Sérsic index m. Again, this expression is equivalent to the expression found by Mazure & Capelato (2002). The asymptotic expansion of the cumulative luminosity for small r can be found in the same way as for the luminosity density in Sect. 5. One finds after some calculation for m < 1 L(r)~43 I0 Re2 bm Γ(1−m) s3.% subequation 2130 0 \begin{equation} L(r) \sim \frac43\,I_0\,R_{\rm e}^2\,b^m\,\Gamma(1-m)\,s^3. \end{equation}(42a)For m = 1 we obtain L(r)~4 I0 Re2 b3 [ln(1s)+(13−γ−lnb+ln2)]s3,% subequation 2130 1 \begin{equation} L(r) \sim \frac{4\,I_0\,R_{\rm e}^2\,b}{3}\, \left[ \ln\left(\frac{1}{s}\right) + \left(\frac13-\gamma-\ln b + \ln 2\right) \right] s^3, \end{equation}(42b)and for m > 1 L(r)~2π I0 Re2 b2m+1 Γ(12−12m)Γ(1−12m) s2+1m.% subequation 2130 2 \begin{equation} L(r) \sim \frac{2\sqrt{\pi}\,I_0\,R_{\rm e}^2\,b}{2m+1}\, \frac{\Gamma\left(\frac12-\frac1{2m}\right)} {\Gamma\left(1-\frac{1}{2m}\right)}\, s^{2+\frac1m}. \end{equation}(42c)These asymptotic expressions can also be found by directly inserting Eqs. (37) into formula (39).

If we assume that mass follows light (or in case the Sérsic model is used to describe the mass density), we can also calculate the (positive) gravitational potential Ψ(r). The most convenient way in the present case is to use the formula Ψ(r)=G Υ∫r∞L(r′) dr′r′2\begin{equation} \Psi(r) = G\,\Upsilon \int_r^\infty \frac{L(r')\,{{\rm d}}r'}{r'^2} \end{equation}(43)where the Υ is the mass-to-light ratio. The result reads as Ψ(r)=G Υ I0 Re(2π)p−1 pq3/2 s×  Gq+1,2p+q+12p+q,1[ |(b2p)2ps2q],\begin{eqnarray} \Psi(r) &=& \frac{G\,\Upsilon\,I_0\,R_{\rm e}}{(2\pi)^{p-1}}\, \frac{\sqrt{p}}{q^{3/2}}\, s \nonumber\\ &&\times\,\, G^{2p+q,1}_{q+1,2p+q+1} \!\left[\left. \begin{array}{l} 1-\frac{1}{q}, 1-\frac{1}{2q}, {{\vec a}} \\ {{\vec b}}, -\frac{1}{2q}, -\frac{1}{q} \end{array} \,\right|\left(\frac{b}{2p}\right)^{2p} s^{2q} \right], \end{eqnarray}(44)or if we introduce the total mass M = Υ   L using expression (30) Ψ(r)=GMRe b2pq(2π)ppq  Γ(2pq) s×  Gq+1,2p+q+12p+q,1[ |(b2p)2ps2q].\begin{eqnarray} \Psi(r) &=& \frac{GM}{R_{\rm e}}\, \frac{b^{\frac{2p}{q}}}{(2\pi)^p\sqrt{pq}\,\, \Gamma\left(\frac{2p}{q}\right)}\, s \nonumber\\ &&\times\,\, G^{2p+q,1}_{q+1,2p+q+1} \Bigg[ \begin{array}{l} 1-\frac{1}{q}, 1-\frac{1}{2q}, {{\vec a}} \\ \label{Psigen} {{\vec b}}, -\frac{1}{2q}, -\frac{1}{q} \end{array} \,\Bigg|\left(\frac{b}{2p}\right)^{2p} s^{2q} \Bigg]. \end{eqnarray}(45)For integer and half-integer values of the Sérsic index m, this expression simplifies to Ψ(r)=GMRe b2m(2π)mm  Γ(2m) s×  G2,2m+22m+1,1[ |(b2m)2ms2].\begin{eqnarray} \Psi(r) & =& \frac{GM}{R_{\rm e}}\, \frac{b^{2m}}{(2\pi)^m\sqrt{m}\,\,\Gamma(2m)}\, s \nonumber\\ &&\times\,\, G^{2m+1,1}_{2,2m+2} \!\left[\left. \begin{array}{l} 0, \frac12 \\ {{\vec b}}, -\frac12, -1 \end{array} \,\right|\left(\frac{b}{2m}\right)^{2m} s^2 \right]. \end{eqnarray}(46)This expression can be reduced slightly further since the coefficient 12\hbox{$\frac{1}{2}$} appears in both the a and b coefficient vectors. Applying Eq. (A.5), the final result reads as Ψ(r)=GMRe b2m(2π)mm  Γ(2m) s G1,2m+12m,1[0b′ |(b2m)2ms2]% subequation 2289 0 \begin{equation} \Psi(r) = \frac{GM}{R_{\rm e}}\, \frac{b^{2m}}{(2\pi)^m\sqrt{m}\,\,\Gamma(2m)}\, s \, G^{2m,1}_{1,2m+1} \!\left[\left. \begin{array}{l} 0 \\ {{\vec b}}' \end{array} \,\right|\left(\frac{b}{2m}\right)^{2m} s^2 \right] \end{equation}(47a)with b′ a vector with 2m + 1 elements given by b′={12m,22m...,2m−12m,−12,−1}.% subequation 2289 1 \begin{equation} {{\vec b}}' = \left\{\frac{1}{2m},\frac{2}{2m}\ldots, \frac{2m-1}{2m}, -\frac12,-1\right\}. \end{equation}(47b)This expression is somewhat simpler than, but equivalent to, expression (28) in Mazure & Capelato (2002).

Since the luminosity density of the Sérsic models never falls more steeply than r-1 at small radii, it is no surprise that all Sérsic models have a finite potential well for all values of m. Ciotti (1991) has already derived an expression for the depth of the potential well using the general expression Ψ0=−4G Υ∫0∞dIdR(R) R dR.\begin{equation} \Psi_0 = -4G\,\Upsilon\int_0^\infty \frac{{{\rm d}}I}{{{\rm d}}R}(R)\,R\,{{\rm d}}R. \end{equation}(48)Applied to the Sérsic model surface brightness profile, the result reads as (his Eq. (12)) Ψ0=4G Υ I0 Re Γ(m+1)bm·\begin{equation} \Psi_0 = 4G\,\Upsilon\,I_0\,R_{\rm e}\,\frac{\Gamma(m+1)}{b^m}\cdot \label{Psi0} \end{equation}(49)Taking the limit r → 0 for the expression (45), we find Ψ0=GMRe 2bmπ Γ(m)Γ(2m),\begin{equation} \Psi_0 = \frac{GM}{R_{\rm e}}\, \frac{2b^{m}}{\pi}\, \frac{\Gamma(m)}{\Gamma(2m)}, \end{equation}(50)equivalent to (49). Using the expansion formulae for the Meijer G function, we can calculate the asymptotic expansion for the potential at small radii. Not surprisingly, we again obtain different expansions for m smaller than, equal to, and larger than 1. After a lengthly calculation, one finds for m < 1 a quadratically decreasing potential, Ψ(r)~GMRe[2bmπ Γ(m)Γ(2m)−b3m3πm Γ(1−m)Γ(2m) s2]·\begin{equation} \Psi(r) \sim \frac{GM}{R_{\rm e}} \left[ \frac{2b^{m}}{\pi}\, \frac{\Gamma(m)}{\Gamma(2m)} - \frac{b^{3m}}{3\pi m}\,\frac{\Gamma(1-m)}{\Gamma(2m)}\,s^2 \right]\cdot \end{equation}(51)For the exponential model m = 1, one obtains Ψ(r)~GMRe{2bπ−b33π[log(1s)+  (12−γ+lnb−ln2−12ψ3/2+12ψ5/2)]s2}\begin{eqnarray} \Psi(r) &\sim& \frac{GM}{R_{\rm e}} \left\{ \frac{2b}{\pi} - \frac{b^3}{3\pi} \left[ \log\left(\frac{1}{s}\right) \right. \right. \nonumber\\ &&+\,\, \left. \left. \left(\frac12 - \gamma + \ln b - \ln 2 -\frac12\psi_{3/2} +\frac12\psi_{5/2} \right) \right] s^2 \right\} \end{eqnarray}(52)where ψz is the digamma function. For m > 1 the potential decreases more softly than quadratically, Ψ(r)~GMRe[2bmπ Γ(m)Γ(2m)−b2m+1π (m+1) (2m+1) Γ(12−12m)Γ(1−12m)Γ(2m) s1+1m].\begin{eqnarray} \Psi(r) &\sim& \frac{GM}{R_{\rm e}} \bigg[ \frac{2b^{m}}{\pi}\, \frac{\Gamma(m)}{\Gamma(2m)} \nonumber\\ &-& \frac{b^{2m+1}}{\sqrt{\pi}\,(m+1)\,(2m+1)}\, \frac{\Gamma\left(\frac12-\frac{1}{2m}\right)} {\Gamma\left(1-\frac{1}{2m}\right)\Gamma(2m)}\, s^{1+\frac{1}{m}} \bigg]. \end{eqnarray}(53)Finally, at large radii, the potential of all Sérsic models falls off as Ψ(r)~GMr,\begin{equation} \Psi(r) \sim \frac{GM}{r}, \end{equation}(54)as required for a system with a finite mass.

7. Conclusions

We have used the Mellin transform technique to derive a closed, analytical expression for the spatial luminosity density ν(r) of the Sérsic model. For general values of the Sérsic parameter m, this expression is a Fox H function. We derived simplified expressions for ν(r) in terms of the Meijer G function for all rational values of m, and for integer values of m our results are equivalent to the expressions found by Mazure & Capelato (2002). Our analytical calculations complement other theoretical studies of the Sérsic model (Ciotti 1991; Ciotti & Lanzoni 1997; Ciotti & Bertin 1999; Trujillo et al. 2001; Cardone 2004; Graham & Driver 2005; Elíasdóttir & Möller 2007) and, given the extended literature on the analytical properties of the Meijer G function, can be used to further examine the properties of this model analytically. We have investigated the asymptotic behaviour of the luminosity density at small and large radii, and find a rich variety in behviour depending on the value of m. We also derived analytical expression for derived quantities, in particular the cumulative luminosity and the gravitational potential. Our results can also be used in practical calculations: as implementations of the Meijer G function are nowadays available both in symbolic computer algebra packages and as high-performance computing code, our results open up the possibility of calculating the luminosity density of the Sérsic models to arbitrary precision.

1

The Wolfram Functions Site (http://functions.wolfram.com/) is a comprehensive online compendium that provides a huge collection of formulas and graphics about mathematical functions. It is created with MathematicaR◯\hbox{$\tiny^{\textregistered}$} and is developed and maintained by Wolfram Research with partial support from the National Science Foundation. A compendium of formulae on the Meijer G function can be found at http://functions.wolfram.com/PDF/MeijerG.pdf.

2

Mpmath (http://code.google.com/p/mpmath/) is a free pure-Python library for multiprecision floating-point arithmetic. It provides an extensive set of transcendental functions, unlimited exponent sizes, complex numbers, interval arithmetic, numerical integration and differentiation, root-finding, linear algebra, and much more.

References

  1. Aceves, H., Velázquez, H., & Cruz, F. 2006, MNRAS, 373, 632 [NASA ADS] [CrossRef] [Google Scholar]
  2. Adamchik, V. 1996, Math. Educ. Res., 5, 16 [Google Scholar]
  3. Allen, P. D., Driver, S. P., Graham, A. W., et al. 2006, MNRAS, 371, 2 [NASA ADS] [CrossRef] [Google Scholar]
  4. Andredakis, Y. C., Peletier, R. F., & Balcells, M. 1995, MNRAS, 275, 874 [NASA ADS] [CrossRef] [Google Scholar]
  5. Binggeli, B., Sandage, A., & Tarenghi, M. 1984, AJ, 89, 64 [NASA ADS] [CrossRef] [Google Scholar]
  6. Caon, N., Capaccioli, M., & D’Onofrio, M. 1993, MNRAS, 265, 1013 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cappellari, M. 2002, MNRAS, 333, 400 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cardone, V. F. 2004, A&A, 415, 839 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Cellone, S. A., Forte, J. C., & Geisler, D. 1994, ApJS, 93, 397 [NASA ADS] [CrossRef] [Google Scholar]
  10. Ciotti, L. 1991, A&A, 249, 99 [NASA ADS] [Google Scholar]
  11. Ciotti, L., & Bertin, G. 1999, A&A, 352, 447 [NASA ADS] [Google Scholar]
  12. Ciotti, L., & Lanzoni, B. 1997, A&A, 321, 724 [NASA ADS] [Google Scholar]
  13. Davies, J. I., Phillipps, S., Cawson, M. G. M., & Disney, M. J., & Kibblewhite, E. J. 1988, MNRAS, 232, 239 [NASA ADS] [CrossRef] [Google Scholar]
  14. De Bruyne, V., Dejonghe, H., Pizzella, A., Bernardi, M., & Zeilinger, W. W. 2001, ApJ, 546, 903 [NASA ADS] [CrossRef] [Google Scholar]
  15. de Vaucouleurs, G. 1948, Ann. Astrophys., 11, 247 [Google Scholar]
  16. D’Onofrio, M., Capaccioli, M., & Caon, N. 1994, MNRAS, 271, 523 [NASA ADS] [CrossRef] [Google Scholar]
  17. Duffy, A. R., Schaye, J., Kay, S. T., & Dal la Vecchia, C. 2008, MNRAS, 390, L64 [NASA ADS] [CrossRef] [Google Scholar]
  18. Elíasdóttir, Á., Möller, O. 2007, J. Cosmol. Astro-Part. Phys., 7, 6 [Google Scholar]
  19. Emsellem, E., Monnet, G., & Bacon, R. 1994a, A&A, 285, 723 [NASA ADS] [Google Scholar]
  20. Emsellem, E., Monnet, G., Bacon, R., & Nieto, J.-L. 1994b, A&A, 285, 739 [NASA ADS] [Google Scholar]
  21. Faber, S. M., & Lin, D. N. C. 1983, ApJ, 266, L17 [NASA ADS] [CrossRef] [Google Scholar]
  22. Fikioris G., 2007, Mellin Transform Method for Integral Evaluation, Introduction and Application to Electromagnetics (Morgan & Claypool) [Google Scholar]
  23. Fox, C. 1961, Transactions of the American Mathematical Society, 98, 395 [Google Scholar]
  24. Gadotti, D. A. 2009, MNRAS, 393, 1531 [NASA ADS] [CrossRef] [Google Scholar]
  25. Gao, L., Navarro, J. F., Cole, S. et al. 2008, MNRAS, 387, 536 [NASA ADS] [CrossRef] [Google Scholar]
  26. Graham, A. W., & Driver, S. P. 2005, Publ. Astron. Soc. Austral., 22, 118 [Google Scholar]
  27. Graham, A. W., & Guzmán, R. 2003, AJ, 125, 2936 [NASA ADS] [CrossRef] [Google Scholar]
  28. Graham, A. W., Merritt, D., Moore, B., Diemand, J., & Terzić, B. 2006, AJ, 132, 2701 [NASA ADS] [CrossRef] [Google Scholar]
  29. Lima Neto, G. B., Gerbal, D., & Márquez., I. 1999, MNRAS, 309, 481 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kilbas, A. A., & Saigo, M. 2004, H-Transforms: Theory and Applications (CRC Press) [Google Scholar]
  31. Marichev O., Handbook of Integral Transforms of Higher Transcendental Functions, Theory and Algorithmic Tables (Horwood, Chichester) [Google Scholar]
  32. Mathai, A. M., & Saxena R. M. 1978, The H Function with Applications in Statistics and Other Disciplines (Wiley) [Google Scholar]
  33. Mathai, A. M., Saxena R. M., & Haubold H. J. 2009, The H-Function: Theory and Applications (Springer) [Google Scholar]
  34. Mazure, A., & Capelato, H. V. 2002, A&A, 383, 384 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Meijer, C. S. 1936, Über Whittakersche bezw. Besselsche Funktionen und deren Produkte, Nieuw Archief voor Wiskunde, 18, 10 [Google Scholar]
  36. Mellier, Y., & Mathez, G. 1987, A&A, 175, 1 [NASA ADS] [Google Scholar]
  37. Méndez-Abreu, J., Aguerri, J. A. L., Corsini, E. M., & Simonneau, E. 2008, A&A, 478, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Merritt, D., Navarro, J. F., Ludlow, A., & Jenkins, A. 2005, ApJ, 624, L85 [NASA ADS] [CrossRef] [Google Scholar]
  39. Merritt, D., Graham, A. W., Moore, B., Diemand, J., & Terzić, B. 2006, AJ, 132, 2685 [NASA ADS] [CrossRef] [Google Scholar]
  40. Möllenhoff, C., & Heidt, J. 2001, A&A, 368, 16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Moore, B., Quinn, T., Governato, F., Stadel, J., & Lake, G. 1999, MNRAS, 310, 1147 [NASA ADS] [CrossRef] [Google Scholar]
  42. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
  43. Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039 [NASA ADS] [CrossRef] [Google Scholar]
  44. Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21 [NASA ADS] [CrossRef] [Google Scholar]
  45. Prugniel, P., & Simien, F. 1997, A&A, 321, 111 [NASA ADS] [Google Scholar]
  46. Sérsic, J. L. 1968, Cordoba, Argentina: Observatorio Astronomico [Google Scholar]
  47. Stadel, J., Potter, D., Moore, B., et al. 2009, MNRAS, 398, L21 [NASA ADS] [CrossRef] [Google Scholar]
  48. Trujillo, I., Graham, A. W., & Caon, N. 2001, MNRAS, 326, 869 [NASA ADS] [CrossRef] [Google Scholar]
  49. Trujillo, I., Asensio Ramos, A., Rubiño-Martín, J. A., et al. 2002, MNRAS, 333, 510 [Google Scholar]
  50. Young, P. J. 1976, AJ, 81, 807 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: The Meijer G and Fox H functions

The Meijer G function (Meijer 1936) is a universal, analytical function that was introduced as a generalization of the hypergeometric series. Nowadays it is more commonly defined as an inverse Mellin transform, i.e. a path integral in the complex plane, Gp,qm,n[ | z]≡Gp,qm,n[ | z]=12πi\appendix \setcounter{section}{1} \begin{eqnarray} && G^{m,n}_{p,q}\!\left[\left. \begin{array}{l} {\vec{a}}\\{\vec{b}} \end{array} \,\right|\, z \right] \equiv G^{m,n}_{p,q}\!\left[\left. \begin{array}{l} a_1,\ldots,a_p\\b_1,\ldots,b_q \end{array} \,\right|\, z \right] = \frac{1}{2\pi i}\nonumber\\ &&\,\,\,\times \int_{\cal{L}} \frac{\prod_{j=1}^m \Gamma(b_j+s)\,\prod_{j=1}^n \Gamma(1-a_j-s)} {\prod_{j=m+1}^q \Gamma(1-b_j-s)\,\prod_{j=n+1}^p \Gamma(a_j+s)} \, z^{-s}\,{{\rm d}}s \label{MeijerG} \end{eqnarray}(A.1)with ℒ a path in the complex plane, and a and b are two vectors of dimension p and q respectively.

The general Meijer G function can reproduce many commonly used special functions, including Bessel functions, elliptic integrals, and hypergeometric functions. Some of the special cases are \appendix \setcounter{section}{1} \begin{eqnarray} &&G^{1,0}_{0,1}\!\left[\left. \begin{array}{l}-\\b\end{array} \,\right|\,z\right] = {\rm e}^{-z}\,z^b \label{G1001} \\ && G^{2,0}_{0,2}\!\left[\left. \begin{array}{l}-\\b_1,b_2\end{array}\,\right|\,z\right] = 2\,z^{\frac12(b_1+b_2)}\,K_{b_1-b_2}\left(2\sqrt{z}\right) \label{G2002} \end{eqnarray}where Kν(z) is the modified Bessel function of the second kind.

The Meijer G function has numerous useful properties that allow transforming expressions to equivalent expressions. For example, one can shift all parameters by a given number, Gp,qm,n[a1,...,apb1,...,bq | z]=z−α Gp,qm,n[a1+α,...,ap+αb1+α,...,bq+α | z].\appendix \setcounter{section}{1} \begin{equation} G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} a_1,\ldots,a_p\\ b_1,\ldots,b_q \end{array} \,\right|\,z\right] = z^{-\alpha}\, G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} a_1+\alpha,\ldots,a_p+\alpha\\ b_1+\alpha,\ldots,b_q+\alpha \end{array} \,\right|\,z\right]. \label{shift} \end{equation}(A.4)An important property is that one can reduce the order of the Meijer G function if one parameter appears in both the upper and lower parameter vectors (depending on the position). For example, if one of the ak with n < k ≤ p equals one of the bj with 1 ≤ j ≤ m, then G p,q m,n [   |   z ] = \appendix \setcounter{section}{1} \begin{eqnarray} &&G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} a_1,\ldots,a_p\\ b_1,\ldots,b_q \end{array} \,\right|\,z\right] =\nonumber\\ &&G^{m-1,n}_{p-1,q-1} \!\left[\left. \begin{array}{l} a_1,\ldots,a_{k-1},a_{k+1},\ldots,a_p \\ b_1,\ldots,b_{j-1},b_{j+1},\ldots,b_q \end{array} \,\right|\,z\right]. \label{schrap} \end{eqnarray}(A.5)Another powerful property that allows the order of the Meijer G function to be reduced in certain cases is G p,q m,n [ a b   |   z ] = k 1 + ν + ( p − q ) / 2 ( 2 π ) ( k − 1 ) δ   G kp,kq km,kn [ a ′ b ′   |   z k k k ( q − p ) ] , \appendix \setcounter{section}{1} % subequation 2754 0 \begin{equation} G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} {\vec {a}} \\ {\vec{b}} \end{array} \,\right|\,z\right] = \frac{k^{1+\nu+(p-q)/2}}{(2\pi)^{(k-1)\delta}}\, G^{km,kn}_{kp,kq} \!\left[\left. \begin{array}{l} {\vec{a}}' \\ {\vec{b}}' \end{array} \,\right|\,\frac{z^k}{k^{k(q-p)}}\right], \label{mult} \end{equation}(A.6a)where k is a positive integer number, δ and ν are defined as δ = m + n + p + q 2 ,   ν = ∑ j = 1 q b j − ∑ j = 1 p a j , \appendix \setcounter{section}{1} % subequation 2754 1 \begin{eqnarray} \delta \!=\! m+n+\frac{p+q}{2}, \qquad \nu \!=\! \sum_{j=1}^{q}b_j - \sum_{j=1}^p a_j, \end{eqnarray}(A.6b)and the vectors a′ and b′ are defined asa′={a1k,a1+1k,...,a1+k−1k,...apk,...,ap+k−1k}b′={b1k,b1+1k,...,b1+k−1k,...bqk,...,bq+k−1k}·\appendix \setcounter{section}{1} % subequation 2754 2 \begin{eqnarray} &&{\vec {a}}' \!=\! \left\{ \frac{a_1}{k}, \frac{a_1+1}{k},\ldots,\frac{a_1+k-1}{k}, \ldots \frac{a_p}{k}, \ldots,\frac{a_p+k-1}{k} \right\} \\ && {\vec {b}}' \!=\! \left\{ \frac{b_1}{k}, \frac{b_1\!+\!1}{k},\ldots,\frac{b_1\!+\!k-1}{k}, \ldots \frac{b_q}{k},\ldots,\frac{b_q\!+\!k-1}{k} \right\}\cdot \end{eqnarray}One of the most powerful properties of the Meijer G function as an analytical tool is that several integrals involving Meijer functions can be evaluated in terms of higher order Meijer functions. For example, ∫Gp,qm,n[ | zw]zα−1 dz=\appendix \setcounter{section}{1} \begin{eqnarray} &&\int G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} a_1,\ldots,a_p \\ b_1,\ldots,b_q \end{array} \,\right|\,zw\right] z^{\alpha-1}\, {{\rm d}}z =\nonumber\\ &&z^{\alpha} G^{m,n+1}_{p+1,q+1} \!\left[\left. \begin{array}{l} 1-\alpha,a_1,\ldots,a_p \\ b_1,\ldots,b_q,-\alpha \end{array} \,\right|\,zw\right] \label{intMeijer} \end{eqnarray}(A.7)

The corresponding definite integral can be evaluated as ∫∞0Gp,qm,n[ | zw]zα−1 dz=\appendix \setcounter{section}{1} \begin{eqnarray} &&\int_0^\infty G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} a_1,\ldots,a_p \\ b_1,\ldots,b_q \end{array} \,\right|\,zw\right] z^{\alpha-1}\, {{\rm d}}z =\nonumber\\[2mm] &&\frac{ \prod_{j=1}^m \Gamma(b_j+\alpha)\, \prod_{j=1}^n \Gamma(1-a_j-\alpha) } { \prod_{j=n+1}^p \Gamma(a_j+\alpha)\, \prod_{j=m+1}^q \Gamma(1-b_j-\alpha) } w^{-\alpha}. \label{defintMeijer} \end{eqnarray}(A.8)

Another useful property is the asymptotic expansion of the Meijer function for small z, in the case of p ≤ q and simple poles, Gp,qm,n[ | z]=∑k=1m􏽑j=1,j≠kmΓ(bj−bk) 􏽑j=1nΓ(1−aj+bk)􏽑j=n+1pΓ(aj−bk) 􏽑j=m+1qΓ(1−bj+bk)\appendix \setcounter{section}{1} \begin{eqnarray} && G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} {\vec {a}} \\ {\vec {b}} \end{array} \,\right|\,z\right] = \sum_{k=1}^m \frac{\prod_{j=1,j\ne k}^m \Gamma(b_j-b_k)\, \prod_{j=1}^n \Gamma(1-a_j+b_k)} {\prod_{j=n+1}^p \Gamma(a_j-b_k)\, \prod_{j=m+1}^q \Gamma(1-b_j+b_k)} \nonumber\\[2mm] &&\times z^{b_k} \left[ 1 + (-1)^{-m-n+p} \frac{\prod_{j=1}^p (1-a_j+b_k)}{\prod_{j=1}^q (1-b_j+b_k)}\, z + {\cal{O}}(z^2) \right]. \label{asymp} \end{eqnarray}(A.9)For the asymptotic expansion at large z, one can use the identity Gp,qm,n[a1,...,apb1,...,bq | z]=Gq,pn,m[1−b1,...,1−bq1−a1,...,1−ap | 1z]·\appendix \setcounter{section}{1} \begin{equation} G^{m,n}_{p,q} \!\left[\left. \begin{array}{l} a_1,\ldots,a_p \\ b_1,\ldots,b_q \end{array} \,\right|\,z\right] = G^{n,m}_{q,p} \!\left[\left. \begin{array}{l} 1-b_1,\ldots,1-b_q \\ 1-a_1,\ldots,1-a_p \end{array} \,\right|\,\frac{1}{z}\right]\cdot \end{equation}(A.10)The Fox H function (Fox 1961) is a generalization of Meijer G function and is also defined as an inverse Mellin transform, Hp,qm,n[ | z]≡Hp,qm,n[ | z]=12πi\appendix \setcounter{section}{1} \begin{eqnarray} &&H^{m,n}_{p,q}\!\left[\left. \begin{array}{l} ({\vec {a}},{\vec{A}}) \\ ({\vec{b}},{\vec{B}}) \end{array} \,\right|\, z \right] \equiv H^{m,n}_{p,q}\!\left[\left. \begin{array}{l} (a_1,A_1),\ldots,(a_p,A_p)\\(b_1,B_1),\ldots,(b_q,B_q) \end{array} \,\right|\, z \right] \nonumber\\ &&= \frac{1}{2\pi {\rm i}} \nonumber\\ &&\times\int_{\cal{L}}\!\! \frac{\prod_{j=1}^m \Gamma(b_j+B_js)\,\prod_{j=1}^n \Gamma(1-a_j-A_js)} {\prod_{j=m+1}^q \Gamma(1-b_j-B_js)\prod_{j=n+1}^p \Gamma(a_j+A_js)} z^{-s}\,{{\rm d}}s. \label{defH} \end{eqnarray}(A.11)Not surprisingly, the Fox H function shares many of the properties of the Meijer G functions, and complete volumes have been written about its identities, asymptotic properties, expansion formulae, and integral transforms (e.g. Mathai & Saxena 1978; Kilbas & Saigo 2004; Mathai et al. 2009).

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.