© ESO, 2012

1. Introduction

The Λ cold dark matter (CDM) model, with (Ω_m, Ω_Λ) = ⁽0.3, 0.7⁾, has become the standard theory of cosmological structure formation. ΛCDM seems to agree with the observations on cluster-sized scales (Primack 2003); however, on galaxy/sub-galaxy scales there appears to be a discrepancy between observations and numerical simulations. High-resolution observations of rotation curves, in particular of low surface brightness (LSB) and dark matter dominated dwarf galaxies (de Blok et al. 2001; van den Bosch & Swaters 2001; Swaters et al. 2003; Weldrake et al. 2003; Donato et al. 2004; Gentile et al. 2005; Simon et al. 2005; Gentile et al. 2007; Banerjee et al. 2010) favour density profiles with a flat central core (e.g. Burkert 1995; Salucci & Burkert 2000; Gentile et al. 2004; Li & Chen 2009). In contrast, N-body (dark matter only) CDM simulations predict galactic density profiles that are too high in the centre (e.g. Navarro et al. 1996, 1997, NFW; Moore et al. 1999). This discrepancy is called the cusp-core problem; a complete review can be found in de Blok (2010).

Gravitational lensing is one of the most powerful tools in observational cosmology for probing the matter distribution of galaxies and clusters in the strong regime (Kochanek et al. 1989; Wambsganss & Paczynski 1994; Bartelmann 1996; Chae et al. 1998; Kochanek et al. 2000; Keeton & Madau 2001; Sand et al. 2002; Keeton 2002, 2003; Keeton & Zabludoff 2004; Limousin et al. 2008; Anguita et al. 2009; Zitrin et al. 2011a,b) and the weak regime (Kaiser & Squires 1993; Mellier 1999; Bartelmann & Schneider 2001; Hoekstra et al. 2004; Clowe et al. 2006; Mahdavi et al. 2007; Jee et al. 2009; Huang et al. 2011). Comparing these observations to theoretical models provides key information to help resolve the cusp-core problem. Evidently, one must use the most accurate density profile to obtain the best fit to observational data from strong- and weak-lensing studies.

Recently, N-body CDM simulations (Navarro et al. 2004; Merritt et al. 2006; Gao et al. 2008; Hayashi & White 2008; Stadel et al. 2009; Navarro et al. 2010) have found that certain three-parameter profiles provide an excellent fit to a wide range of dark matter haloes. One of these is the Einasto (1965) profile, a three-dimensional version of the two-dimensional Sérsic (1968) profile used to describe the surface brightness of early-type galaxies and the 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; Graham et al. 2006; Gadotti 2009). The Sérsic profile can be written as: ${\mathrm{\Sigma }}_{\mathrm{S}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{=}\mathrm{{\rm Y}}{\mathit{I}}_{\mathrm{e}}\exp \begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{-}{\mathit{b}}_{\mathrm{m}}\left[{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{e}}}\right)}^{\mathrm{1}\mathit{/}\mathit{m}}\mathrm{-}\mathrm{1}\right]\hspace{0.17em}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}$(1)where R is the distance in the sky plane, m the Sérsic index, Υ is the mass-to-light ratio, I_e is the luminosity density at the effective radius R_e, and b_m is a dimensionless function of m that can be determined from the condition that the luminosity inside R_e equals half of the total luminosity. Numerical solutions for b_m are given by Ciotti (1991), Moriondo et al. (1998), Prugniel & Simien (1997) and an asymptotic expansion ${\mathit{b}}_{\mathrm{m}}\mathrm{=}\mathrm{2}\mathit{m}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}\mathrm{+}\mathrm{4}\mathit{/}\mathrm{405}\mathit{m}\mathrm{+}\mathrm{𝒪}{}^{\mathrm{(}}\mathit{m}^{\mathrm{2}}{}^{\mathrm{)}}$ using analytical methods was obtained by Ciotti & Bertin (1999).

The Einasto profile (model) is characterised by a power-law logarithmic slope, $\mathit{\gamma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\equiv }\mathrm{-}\frac{\mathrm{d}\ln \mathit{\rho }}{\mathrm{d}\ln \mathit{r}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\propto }{\mathit{r}}^{\mathrm{1}\mathit{/}\mathit{n}}\mathit{,}$(2)with n, which we call the Einasto index, a positive number defining the steepness of the power law. Integrating leads to the general density profile $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{s}}\exp \begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{-}{\mathit{d}}_{\mathit{n}}\left[{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)}^{\mathrm{1}\mathit{/}\mathit{n}}\mathrm{-}\mathrm{1}\right]\hspace{0.17em}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}$(3)where r_s represents the radius of the sphere that contains half of the total mass, ρ_s is the mass density at r = r_s, and d_n is a numerical constant that ensures that r_s is indeed the half-mass radius. In the context of dark matter haloes, the density can also be expressed as $\mathit{\rho }\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}{\mathit{\rho }}_{-2}\hspace{0.17em}\exp \begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{-}\mathrm{2}\mathit{n}\left[{\left(\frac{\mathit{r}}{{\mathit{r}}_{-2}}\right)}^{\mathrm{1}\mathit{/}\mathit{n}}\mathrm{-}\mathrm{1}\right]\hspace{0.17em}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}$(4)where ρ_-2 and r_-2 are the density and radius at which $\mathit{\rho }{}^{\mathrm{\right(}}\mathit{r}^{\mathrm{\left)}}\mathrm{\propto }{\mathit{r}}^{-2}$. In the remainder of this paper, we will use yet another, equivalent version, $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\exp \left[\mathrm{-}{\left(\frac{\mathit{r}}{\mathit{h}}\right)}^{\mathrm{1}\mathit{/}\mathit{n}}\right]\mathit{,}$(5)where we introduced the central density ${\mathit{\rho }}_{\mathrm{0}}\mathrm{=}{\mathit{\rho }}_{\mathit{s}}\hspace{0.17em}{\mathrm{e}}^{{\mathit{d}}_{\mathit{n}}}\mathrm{=}{\mathit{\rho }}_{-2}\hspace{0.17em}{\mathrm{e}}^{\mathrm{2}\mathit{n}}\mathit{,}$(6)and the scale length $\mathit{h}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{s}}}{{\mathit{d}}_{\mathit{n}}^{\mathit{n}}}\mathrm{=}\frac{{\mathit{r}}_{-2}}{\mathrm{(}\mathrm{2}\mathit{n}{\mathrm{)}}^{\mathit{n}}}\mathrm{·}$(7)If a model is to describe real galactic systems, several conditions must be imposed on the model description functions, as discussed by Einasto (1969a). When a model is constructed, an initial descriptive function is chosen; the most practical choice is the density profile, because the main descriptive functions (cumulative mass profile, gravitational potential, surface mass density) are integrals of the density profile. Furthermore, a physical model has to satisfy several conditions: i) the density profile must be non-negative and finite; ii) the density must be a smoothly decreasing function that approaches zero at large radii; iii) some moments of the mass function must be finite, in particular moments that define the central gravitational potential, the total mass, and the effective radius of the system; and iv) the descriptive functions must not have jump discontinuities. Einasto (1969a) presented several families of valid descriptive functions, among which the so-called Einasto profile is a special case that agrees best with observations.

The Einasto profile was used by Einasto (1969b) to obtain a model of M31. Later, this model was applied by Einasto (1974) to several nearby galaxies, including M32, M87, Fornax and Sculptor dwarfs, and also M31 and the Milky Way. These models were multi-component ones, and each component has its own parameter set {ρ₀, h, n} representing certain physically homogeneous stellar populations.

Navarro et al. (2004) found that for haloes with masses from dwarfs to clusters, 4.54 ≲ n ≲ 8.33 with an average value of n = 5.88. Hayashi & White (2008) and Gao et al. (2008) found that n tends to decrease with mass and redshift, with n ~ 5.88 for galaxy-sized and n ~ 4.35 for cluster-sized haloes in the Millennium Run (MR) (Springel et al. 2005). Navarro et al. (2010) obtained similar results for galaxy-sized haloes in the Aquarius simulation (Springel et al. 2008). Also, Gao et al. (2008) showed that n ~ 3.33 for the most massive haloes of MR. Chemin et al. (2011) modelled the rotation curves of a spiral galaxies subsample from THINGS (The HI Nearby Galaxy Survey, de Blok et al. 2008), using the Einasto profile, and found that n tends to have lower values by a factor of two or more compared to those predicted by N-body simulations. Dhar & Williams (2012) fitted the surface brightness profiles of a large sample of elliptical galaxies in the Virgo cluster, using a multi-component Einasto profile consisting of two or three components for each galaxy with different n for each individual component. For the central component, they found values of n ≲ 1 for the most massive and shallow-cusp galaxies and n < 2 for steep-cusp galaxies, while 5 ≲ n ≲ 8 for the outer component of massive galaxies. These values for the outer component are consistent with the results from N-body simulations for Einasto dark matter haloes, and the authors argued that Einasto components with an n in this range could be dark matter dominated.

In the light of its increasing popularity to describe the density of simulated dark matter haloes, a detailed investigation of the properties of the Einasto model is of paramount importance. Some aspects of the Einasto model have been presented by several authors (Mamon & Łokas 2005; Cardone et al. 2005; Merritt et al. 2006; Dhar & Williams 2010). The most complete study of the properties of the Einasto model is the work by Cardone et al. (2005), who provided a set of analytical expressions for quantities such as the mass profile and gravitational potential and discussed the dynamical structure for both isotropic and anisotropic cases. Nevertheless, the Einasto model has not been studied analytically as extensively as the Sérsic models, and several properties still have to be further investigated in more detail. One area where progress is still to be made is on the value of the dimensionless parameter d_n.

The most important lacuna, however, concerns the surface density on the plane of the sky, an important quantity that defines the lensing properties of a dark matter halo. The surface mass density has an importance in theoretical predictions and observations. On the observational side, based mostly on the mass decomposition of rotation curves using cored haloes, various studies (Kormendy & Freeman 2004; Spano et al. 2008; Donato et al. 2009; Gentile et al. 2009) have shown that the product of the central density ρ₀ and the core radius r₀ is consistent with a universal value, independent of galaxy mass. The product ρ₀r₀ is directly proportional to the average surface density within r₀, and to the gravitational acceleration of dark matter at r₀. In view of a detailed comparison of these results with the outcome of the most recent dark matter simulations, it is therefore important to study the analytical properties of the surface density distribution of Einasto haloes. Recent works studying theoretical and observational aspects of dark matter surface densities include Boyarsky et al. (2010), Walker et al. (2010), Napolitano et al. (2010), and Cardone et al. (2011): efforts are being made to confirm or call into question the universality of the dark matter surface density within one dark halo scale length.

Cardone et al. (2005) showed that the surface density of the Einasto model cannot be expressed in terms of elementary functions and discussed the general properties using numerical integration. Dhar & Williams (2010) presented an analytical approximation for the Einasto surface brightness profile and demonstrated that this surface density profile is not Sérsic-like. A general analytical formula, which would enable an arbitrary-precision calculation and an analytical study of the asymptotic behaviour, is still lacking to date.

The most recent studies in strong and weak lensing (e.g. Donnarumma et al. 2011; Sereno & Umetsu 2011; Morandi et al. 2011; Umetsu et al. 2011; Oguri et al. 2012) use the NFW profile or its generalization (Zhao 1996; Jing & Suto 2000) to model the dark matter halo instead of the Einasto profile. One of the factors that limits its adoption is the absence of analytical formulas for its lensing properties. Therefore, its application in lensing studies is not as wide as for other profiles. A complete general set of analytical formulas for the lensing properties of the Einasto profile would help to increase its use in lensing studies.

In this paper, we extend the analytical study of the Einasto model using some of the techniques also employed in the extensive literature on the analytical properties of the Sérsic model (e.g. Ciotti 1991; Ciotti & Lanzoni 1997; Ciotti & Bertin 1999; Trujillo et al. 2001; Mazure & Capelato 2002; Cardone 2004; Graham & Driver 2005; Elíasdóttir & Möller 2007; Baes & Gentile 2011; Baes & van Hese 2011). In Sect. 2, we discuss an analytical expansion for the dimensionless parameter d_n. In Sect. 3, we derive an analytical expression for the Einasto surface mass density in terms of the Fox H function, using the Mellin transform-method. We then use the result for the projected surface mass density to calculate the cumulative surface mass, lens equation, deflection angle, deflection potential, magnification, shear and the critical curves for a spherically symmetric lens described by the Einasto profile in terms of this function. We also calculate explicit series expansions for the surface mass density and all the lensing properties. In Sect. 4, we derive some special cases for the lensing properties for integer and half-integer values of n in terms of the Meijer G function. Furthermore, we compare the Einasto and Sérsic surface mass densities for the same values of their respective indices. Finally, in Sect. 5 we discuss the asymptotic behaviour of the surface mass density and cumulative surface mass at small and large radii. We present our conclusions in Sect. 6.

2. Spatial properties

The total mass of an Einasto model with a mass density profile given by (5), is $\mathit{M}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}\mathit{n}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}\mathit{.}$(8)If we use this formula to replace the central density ρ₀ by the total mass M as a parameter in the definition of the Einasto models, we obtain $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{M}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}\mathit{n}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{s}}^{\mathrm{1}\mathit{/}\mathit{n}}}\mathit{,}$(9)where we have introduced the reduced radius $\mathit{s}\mathrm{=}\frac{{\mathrm{\left(}{\mathit{d}}_{\mathit{n}}\mathrm{\right)}}^{\mathit{n}}\hspace{0.17em}\mathit{r}}{{\mathit{r}}_{\mathit{s}}}\mathrm{·}$(10)At small radii, the density profile behaves as $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{M}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}\mathit{n}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{s}}^{\mathrm{1}\mathit{/}\mathit{n}}\mathrm{+}{\mathit{...}}^{\mathrm{)}}\mathit{.}$(11)The cumulative mass profile M(r) of a spherical mass distribution is found through the equation $\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathit{.}$(12)For the Einasto model, we find $\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{M}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n,}{\mathit{s}}^{\mathrm{1}\mathit{/}\mathit{n}}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\right]\mathit{,}$(13)where Γ(α,x) is the incomplete gamma function, $\mathrm{\Gamma }\mathrm{\left(}\mathit{\alpha ,x}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{x}}^{\mathrm{\infty }}{\mathit{t}}^{\mathit{\alpha }\mathrm{-}\mathrm{1}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{t}}\hspace{0.17em}\mathrm{d}\mathit{t}\mathit{.}$(14)Given that this is the radius of the sphere that encloses half of the total mass, we find that d_n is the solution of the equation $\mathrm{2}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n,}{\mathit{d}}_{\mathit{n}}\mathrm{\right)}\mathrm{=}\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}\mathit{.}$(15)This equation cannot be solved in a closed form; one option is to solve it numerically or use interpolation formulae. For example, Merritt et al. (2006), quoting G. A. Mamon (priv. comm.), propose ${\mathit{d}}_{\mathit{n}}\mathrm{\approx }\mathrm{3}\mathit{n}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{+}\mathrm{0.0079}\mathit{/}\mathit{n}$, but do not mention the origin of this approximation or quote its accuracy. An elegant way to determine an approximation for d_n builds on the work by Ciotti & Bertin (1999), who used an asymptotic expansion method to solve the general relation Γ(α,x) = Γ(α)/2. If we apply the resulting expansion to our problem, we obtain for the parameter d_n from the Einasto model $\begin{array}{ccc}{\mathit{d}}_{\mathit{n}}\mathrm{\approx }\mathrm{3}\mathit{n}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{+}\frac{\mathrm{8}}{\mathrm{1215}\hspace{0.17em}\mathit{n}}\mathrm{+}\frac{\mathrm{184}}{\mathrm{229635}\hspace{0.17em}{\mathit{n}}^{\mathrm{2}}}& & \\ \mathrm{+}\frac{\mathrm{1048}}{\mathrm{31000725}\hspace{0.17em}{\mathit{n}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{17557576}}{\mathrm{1242974068875}\hspace{0.17em}{\mathit{n}}^{\mathrm{4}}}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{n}}^{\mathrm{5}}}\right)\mathrm{·}& & \end{array}$(16)The gravitational potential of a spherically symmetric mass distribution ρ(r) can be found through the formula $\mathrm{\Psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi G}\left[\frac{\mathrm{1}}{\mathit{r}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}{\mathrm{\int }}_{\mathit{r}}^{\mathrm{\infty }}\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{r}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\right]\mathit{,}$(17)or equivalently $\mathrm{\Psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{r}}^{\mathrm{\infty }}\frac{\mathit{M}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}}{{\mathit{r}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{·}$(18)For the Einasto model, we find $\mathrm{\Psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{GM}}{\mathit{h}}\hspace{0.17em}{\mathit{s}}^{-1}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n,}{\mathit{s}}^{\mathrm{1}\mathit{/}\mathit{n}}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\mathrm{+}\frac{\mathit{s}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathrm{2}\mathit{n,}{\mathit{s}}^{\mathrm{1}\mathit{/}\mathit{n}}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\right]\mathrm{·}$(19)The Einasto model has a finite potential well, given by ${\mathrm{\Psi }}_{\mathrm{0}}\mathrm{=}\frac{\mathit{GM}}{\mathit{h}}\hspace{0.17em}\frac{\mathrm{\Gamma }\mathrm{\left(}\mathrm{2}\mathit{n}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\mathrm{·}$(20)At small radii, the potential decreases parabolically $\mathrm{\Psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{~}\frac{\mathit{GM}}{\mathit{h}}\hspace{0.17em}\left[\frac{\mathrm{\Gamma }\mathrm{\left(}\mathrm{2}\mathit{n}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\mathrm{-}\frac{{\mathit{s}}^{\mathrm{2}}}{\mathrm{6}\mathit{n}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathrm{3}\mathit{n}\mathrm{\right)}}\mathrm{+}\mathit{...}\right]\mathit{,}$(21)which is not surprising, given the finite density core of the Einasto density profile. At large radii, we obtain a Keplerian fall-off, $\mathrm{\Psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{~}\frac{\mathit{GM}}{\mathit{r}}\mathrm{+}\mathit{...}\mathit{,}$(22)which is also expected for a model with a finite total mass.

3. Lensing properties

3.1. Surface mass density and cumulative surface mass density

The surface mass density of a spherically symmetric lens is given by integrating along the line of sight of the 3D density profile $\mathrm{\Sigma }\left(\mathit{\xi }\right)\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\mathit{\rho }\left(\mathit{\xi ,r}\right)\hspace{0.17em}\mathrm{d}\mathit{z,}$(23)where ξ is the radius measured from the centre of the lens and $\mathit{r}\mathrm{=}\sqrt{{\mathit{\xi }}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$. This expression can also be written as an Abel transform (Binney & Tremaine 1987) $\mathrm{\Sigma }\left(\mathit{\xi }\right)\mathrm{=}\mathrm{2}\hspace{0.17em}{\mathrm{\int }}_{\mathit{\xi }}^{\mathrm{\infty }}\frac{\mathit{\rho }\mathrm{(}{\mathit{r}}^{\mathrm{)}}\hspace{0.17em}\mathit{r}\hspace{0.17em}\mathrm{d}\mathit{r}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{\xi }}^{\mathrm{2}}}}\mathrm{·}$(24)Inserting Eq. (5) into the above expression, we obtain $\mathrm{\Sigma }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\mathrm{2}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}{\mathrm{\int }}_{\mathit{x}}^{\mathrm{\infty }}\frac{{\mathrm{e}}^{\mathrm{-}{\mathit{s}}^{\mathrm{1}\mathit{/}\mathit{n}}}\hspace{0.17em}\mathit{s}\hspace{0.17em}\mathrm{d}\mathit{s}}{\sqrt{{\mathit{s}}^{\mathrm{2}}\mathrm{-}{\mathit{x}}^{\mathrm{2}}}}\mathit{,}$(25)where we have introduced the quantities x = ξ/h and s = r/h.

As discussed by Cardone et al. (2005) and Dhar & Williams (2010), the integral (25) cannot be expressed in terms of elementary or even the most regular functions for all the values of n. Only the central surface mass density can be evaluated analytically as $\mathrm{\Sigma }\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{n}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}\mathit{.}$(26)This situation is very similar to the deprojection of the Sérsic surface brightness profile (1). Indeed, the deprojection formula leads to an integral that is an inverse Abel transform, and also in this case it was long thought that no further analytical progress was possible. However, Mazure & Capelato (2002) used the computer package Mathematica to demonstrate that it is possible to write the Sérsic luminosity density in terms of the Meijer G function for all integer Sérsic indices m. Baes & Gentile (2011) and Baes & van Hese (2011) took this analysis one step further and showed that the deprojection of the Sérsic surface brightness profile for general values of m can be solved elegantly using Mellin integral transforms and gives rise to a Mellin-Barnes integral. The result is that the Sérsic luminosity density can be written compactly in terms of a Fox H function, which reduces to a Meijer G function for all rational values of m.

The obvious similarity between these two cases invites us to apply the same Mellin transform technique (Marichev 1982; Adamchick 1996; Fikioris 2007) to the integral (25). The basic idea behind this technique is that any definite integral $\mathit{f}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{g}\mathrm{\left(}\mathit{t,z}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{t,}$(27)can be written as $\mathit{f}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}{\mathit{f}}_{\mathrm{2}}\left(\frac{\mathit{z}}{\mathit{t}}\right)\hspace{0.17em}\frac{\mathrm{d}\mathit{t}}{\mathit{t}}\mathrm{·}$(28)This expression is exactly a Mellin convolution of two functions f₁ and f₂. Now we can apply the Mellin convolution theorem, which states that the Mellin transform of a Mellin convolution is equal to the products of the Mellin transforms of the original functions. As a result, the definite integral (27) can be written as the inverse Mellin transform of the product of the Mellin transforms of f₁ and f₂. The Mellin transform and its inverse are defined as $\begin{array}{ccc}{M}_{\mathit{f}}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\mathit{\phi }\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{f}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}{\mathit{z}}^{\mathit{u}\mathrm{-}\mathrm{1}}\hspace{0.17em}\mathrm{d}\mathit{z,}& & \\ M_{\mathit{\phi }}^{-1}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi i}}{\mathrm{\int }}_{\mathrm{ℒ}}\mathit{\phi }\mathrm{\left(}\mathit{u}\mathrm{\right)}\hspace{0.17em}{\mathit{z}}^{\mathrm{-}\mathit{u}}\hspace{0.17em}\mathrm{d}\mathit{u,}& & \end{array}$with ℒ a line integral over a vertical line in the complex plane. This implies that $\mathit{f}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi i}}{\mathrm{\int }}_{\mathrm{ℒ}}{M}_{{\mathit{f}}_{\mathrm{1}}}\mathrm{\left(}\mathit{u}\mathrm{\right)}\hspace{0.17em}{M}_{{\mathit{f}}_{\mathrm{2}}}\mathrm{\left(}\mathit{u}\mathrm{\right)}\hspace{0.17em}{\mathit{z}}^{\mathrm{-}\mathit{u}}\hspace{0.17em}\mathrm{d}\mathit{u}\mathit{.}$(31)It now turns out that the integral (31) is a Mellin-Barnes integral for large classes of the functions f₁ and f₂, and that this integral can be evaluated as a Fox H function or a Meijer G function in many cases.

We can immediately apply this formalism to the integral (25), with z = 1 and ${\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{2}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{t}}^{\mathrm{1}\mathit{/}\mathit{n}}}{\mathit{t}}^{\mathrm{2}}\mathit{,}$(32)and ${\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathrm{0}\mathrm{\le }\mathit{t}\mathrm{\le }{\mathit{x}}^{-1}\mathit{,}\\ & \mathrm{else}\mathit{.}\end{array}$(33)The Mellin transforms of these functions are readily calculated $\begin{array}{ccc}{M}_{{\mathit{f}}_{\mathrm{1}}}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\mathrm{2}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}\mathit{n}\hspace{0.17em}\mathrm{\Gamma }\hspace{0.17em}\mathrm{\left[}\mathit{n}\mathrm{\right(}\mathrm{2}\mathrm{+}\mathit{u}\mathrm{\left)}\mathrm{\right]}\mathit{,}& & \\ M_{{\mathit{f}}_{\mathrm{2}}}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\frac{\sqrt{\mathit{\pi }}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\frac{\mathrm{1}\mathrm{+}\mathit{u}}{\mathrm{2}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}{\frac{\mathit{u}}{\mathrm{2}}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{u}\hspace{0.17em}{\mathit{x}}^{\mathrm{1}\mathrm{+}\mathit{u}}}\mathrm{·}& & \end{array}$Combining these results, we obtain $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{x}\mathrm{\right)}& & \mathrm{=}\mathrm{2}\mathit{n}\hspace{0.17em}\sqrt{\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}\mathit{x}\hspace{0.17em}\\ & & \mathrm{\times }\hspace{0.17em}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathrm{i}}{\mathrm{\int }}_{\mathrm{ℒ}}\frac{\mathrm{\Gamma }\left(\mathrm{2}\mathit{ny}\right)\mathrm{\Gamma }\mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{)}}}{\mathrm{\Gamma }\left(\mathit{y}\right)}{\mathrm{[}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\mathrm{]}}^{\mathrm{-}\mathit{y}}\mathrm{d}\mathit{y}\mathit{.}\end{array}$(36)This Mellin-Barnes integral may be recognized as a Fox H function, which is generally defined as the inverse Mellin transform of a product of gamma-functions, $\begin{array}{ccc}{\mathit{H}}_{\mathit{p,q}}^{\mathit{m,n}}\left[\begin{array}{c}\\ \\ \end{array}\hspace{0.17em}|\hspace{0.17em}\mathit{z}\right]& & \mathrm{=}\\ \mathit{label}\mathit{eq}\mathrm{:}\mathit{defH}& & \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathrm{i}}{\mathrm{\int }}_{\mathrm{ℒ}}\frac{{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{m}}\mathrm{\Gamma }\mathrm{(}{\mathit{b}}_{\mathit{j}}\mathrm{+}{\mathit{B}}_{\mathit{j}}\mathit{s}\mathrm{)}{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathit{j}}\mathrm{-}{\mathit{A}}_{\mathit{j}}\mathit{s}\mathrm{)}}{{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathit{m}\mathrm{+}\mathrm{1}}^{\mathit{q}}\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{b}}_{\mathit{j}}\mathrm{-}{\mathit{B}}_{\mathit{j}}\mathit{s}\mathrm{)}{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathit{n}\mathrm{+}\mathrm{1}}^{\mathit{p}}\mathrm{\Gamma }\mathrm{(}{\mathit{a}}_{\mathit{j}}\mathrm{+}{\mathit{A}}_{\mathit{j}}\mathit{s}\mathrm{)}}\hspace{0.17em}{\mathit{z}}^{\mathrm{-}\mathit{s}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathit{.}\end{array}$(37)Using this definition, we can write the surface density of the Einasto model in the following compact form, $\mathrm{\Sigma }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\mathrm{2}\mathit{n}\hspace{0.17em}\sqrt{\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{H}}_{\mathrm{1}\mathit{,}\mathrm{2}}^{\mathrm{2}\mathit{,}\mathrm{0}}\left[\begin{array}{c}\\ \mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{\right)}\\ \mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{2}\mathit{n}\mathrm{\right)}\mathit{,}\mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\end{array}|\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\right]\mathit{.}$(38)The Fox H function is a general analytical function that is becoming more and more used both in mathematics and applied sciences; in fact, it is scheduled for inclusion in the Mathematica numerical library. While not the most mainstream special function, the Fox H is an extremely flexible function that contains very broad classes of elementary and special functions as particular cases. It is indeed a very powerful tool for analytical work. It has many general properties that allow one to manipulate expressions to equivalent forms, reduce the order for certain values of the parameters, etc. For details on its many useful properties, we refer the interested reader to Mathai (1978), Kilbas & Saigo (2004), or Mathai et al. (2009) and the references therein.

As an illustration of the power of the Fox H function as an analytical tool, and as a sanity check on the formula, we can calculate the total mass of the Einasto model by integrating the surface mass density over the plane of the sky, $\begin{array}{ccc}\mathit{M}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{\Sigma }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\hspace{0.17em}\mathit{\xi }\hspace{0.17em}\mathrm{d}\mathit{\xi }& \mathrm{=}& \mathrm{2}\mathit{n}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{H}}_{\mathrm{1}\mathit{,}\mathrm{2}}^{\mathrm{2}\mathit{,}\mathrm{0}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}\mathit{t}\right]\hspace{0.17em}{\mathit{t}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{d}\mathit{t}\mathit{.} \end{array}$(39)This integral can be calculated by setting $\mathit{s}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}$ and a = 1 in Eq. (2.8) in Mathai et al. (2009). We immediately find $\mathit{M}\mathrm{=}\mathrm{2}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}\mathit{n}\hspace{0.17em}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{3}{\mathit{n}}^{\mathrm{)}}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{\right)}}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}\mathit{n}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}\mathrm{3}{\mathit{n}}^{\mathrm{)}}\mathit{,}$(40)in agreement with Eq. (8).

An important quantity for gravitational lensing studies is the cumulative surface mass density, i.e. the total mass contained in a infinite cylinder with radius ξ, $\mathit{M}\left(\mathit{\xi }\right)\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\xi }}\mathrm{\Sigma }\mathrm{\left(}{\mathit{\xi }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{\xi }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}{\mathit{\xi }}^{\mathrm{\prime }}\mathit{.}$(41)We find $\begin{array}{ccc}\mathit{M}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \mathrm{2}\mathit{n}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}{\mathit{x}}^{\mathrm{3}}\hspace{0.17em}\\ & & \end{array}$(42)As another demonstration of the usefulness of the Fox H function, we can apply the residue theorem to the contour integral (37), and obtain explicit series expansions. The general procedure can be found in Kilbas & Saigo (1999), and a specific application to the deprojected Sérsic model is given in Baes & van Hese (2011). The analysis for the projected Einasto profile is completely analogous: again, the form of the series expansion depends on the multiplicity of the poles of the gamma functions Γ(b_j + B_js). For both Σ(x) and M(x), the poles of these gamma functions are −k₁/2n and 1/2 − k₂, with k₁ and k₂ any natural number. We encounter two cases.

Case 1: if n is either non-rational or a rational number p/q with an even denominator (and p,q coprime), all poles are simple and the expansion is a power series, $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{n}\sqrt{\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}[\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathit{k}\mathit{/}\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{2}\mathit{n}}& & \\ \mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathit{nk}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}\mathit{k}}]\mathit{,}& & \end{array}$(43)and $\begin{array}{ccc}\mathit{M}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{n}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}[\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathit{k}\mathit{/}\mathit{n}\mathrm{+}\mathrm{3}}}{\mathrm{2}\mathit{n}}& & \\ \mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathit{nk}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}\mathit{k}\mathrm{+}\mathrm{2}}]\mathit{.}& & \end{array}$(44)Case 2: if n is integer or a rational number p/q with an odd denominator, some poles are of second order, and the expansion is a logarithmic-power series. If we define ${\mathit{k}}_{\mathrm{0}}\mathrm{=}\frac{\mathit{q}\mathrm{-}\mathrm{1}}{\mathrm{2}}$, then we obtain after some algebra $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \mathrm{2}\mathit{n}\sqrt{\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}[\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{1}\mathit{k}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{p}\ne \mathrm{0}\end{array}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathit{k}\mathit{/}\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{2}\mathit{n}}\\ & & \mathrm{+}\frac{\mathrm{\Gamma }\mathrm{\left(}\mathit{n}\mathrm{\right)}}{\sqrt{\mathit{\pi }}}\mathrm{+}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{1}\mathrm{(}\mathit{k}\mathrm{+}{\mathit{k}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{q}\ne \mathrm{0}\end{array}\frac{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathit{nk}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}\mathit{k}}]\\ & & \mathrm{+}\mathrm{2}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{0}\mathrm{(}\mathit{k}\mathrm{+}{\mathit{k}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{q}\mathrm{=}\mathrm{0}\end{array}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{p}}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{)}\mathrm{!}}{\mathrm{(}\mathrm{2}\mathit{nk}\mathrm{-}\mathit{n}\mathrm{)}\mathrm{!}\hspace{0.17em}\mathit{k}\mathrm{!}\hspace{0.17em}\mathit{k}\mathrm{!}}\hspace{0.17em}{\left(\frac{\mathit{x}}{\mathrm{2}}\right)}^{\mathrm{2}\mathit{k}}\mathrm{\times }[\mathrm{-}\ln \left(\frac{\mathit{x}}{\mathrm{2}}\right)\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{k}}\mathrm{+}\mathit{\psi }\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}\mathit{n}\hspace{0.17em}\mathit{\psi }\mathrm{(}\mathrm{2}\mathit{nk}\mathrm{-}\mathit{n}\mathrm{)}\mathrm{-}\mathit{\psi }\mathrm{(}\mathrm{2}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}]\mathit{,}\end{array}$(45)and $\begin{array}{ccc}& & \mathit{M}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{n}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}[\mathrm{-}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{1}\mathit{k}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{p}\ne \mathrm{0}\end{array}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathit{k}\mathit{/}\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{2}\mathit{n}}\\ & & \hspace{0.17em}\mathrm{+}\frac{\mathrm{\Gamma }\mathrm{\left(}\mathit{n}\mathrm{\right)}}{\sqrt{\mathit{\pi }}}\mathrm{+}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{1}\mathrm{(}\mathit{k}\mathrm{+}{\mathit{k}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{q}\ne \mathrm{0}\end{array}\frac{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathit{nk}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}\mathit{k}}]\\ & & \hspace{0.17em}\mathrm{+}\mathrm{2}\mathit{\pi }\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{0}\mathrm{(}\mathit{k}\mathrm{+}{\mathit{k}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{q}\mathrm{=}\mathrm{0}\end{array}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{p}}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{)}\mathrm{!}}{\mathrm{(}\mathrm{2}\mathit{nk}\mathrm{-}\mathit{n}\mathrm{)}\mathrm{!}\hspace{0.17em}\mathit{k}\mathrm{!}\hspace{0.17em}\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}{\left(\frac{\mathit{x}}{\mathrm{2}}\right)}^{\mathrm{2}\mathit{k}}\mathrm{\times }[\mathrm{-}\ln \left(\frac{\mathit{x}}{\mathrm{2}}\right)\\ & & \hspace{0.17em}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{k}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{k}\mathrm{+}\mathrm{2}}\mathrm{+}\mathit{\psi }\mathrm{(}\mathit{k}\mathrm{+}\mathrm{2}\mathrm{)}\mathrm{+}\mathit{n}\hspace{0.17em}\mathit{\psi }\mathrm{(}\mathrm{2}\mathit{nk}\mathrm{-}\mathit{n}\mathrm{)}\mathrm{-}\mathit{\psi }\mathrm{(}\mathrm{2}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}]\mathit{,}\end{array}$(46) with ψ(k) the digamma function.

3.2. Lens equation and deflection angle

In the thin lens approximation, the lens equation for axially symmetric lens is $\mathit{\eta }\mathrm{=}\frac{{\mathit{D}}_{\mathrm{S}}}{{\mathit{D}}_{\mathrm{L}}}\hspace{0.17em}\mathit{\xi }\mathrm{-}{\mathit{D}}_{\mathrm{LS}}\hspace{0.17em}\mathit{\alpha ̂}\mathit{,}$(47)where the quantities η and ξ are the physical positions of the source in the source plane and an image in the image plane, respectively, is the deflection angle, and D_L, D_S and D_LS are the angular distances from observer to lens, from observer to source, and from lens to source, respectively.

With the dimensionless positions y = D_L η/D_S h and x = ξ/h, and dimensionless , the lens equation reduces to $\mathit{y}\mathrm{=}\mathit{x}\mathrm{-}\mathit{\alpha }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}$(48)The deflection angle for a spherical symmetric lens is (Schneider et al. 1992) $\mathit{\alpha }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\mathit{x}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{x}}{\mathit{x}}^{\mathrm{\prime }}\hspace{0.17em}\frac{\mathrm{\Sigma }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}}{{\mathrm{\Sigma }}_{\mathrm{crit}}}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathrm{2}}{\mathit{x}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{x}}{\mathit{x}}^{\mathrm{\prime }}\hspace{0.17em}\mathit{\kappa }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}$(49)where $\mathit{\kappa }\mathrm{=}\frac{\mathrm{\Sigma }\mathrm{\left(}\mathit{x}\mathrm{\right)}}{{\mathrm{\Sigma }}_{\mathrm{crit}}}\mathit{,}$(50)is the convergence and Σ_crit is the critical surface mass density defined by ${\mathrm{\Sigma }}_{\mathrm{crit}}\mathrm{\equiv }\frac{{\mathit{c}}^{\mathrm{2}}\hspace{0.17em}{\mathit{D}}_{\mathrm{S}}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}\mathit{G}\hspace{0.17em}{\mathit{D}}_{\mathrm{L}}\hspace{0.17em}{\mathit{D}}_{\mathrm{LS}}}\mathit{,}$(51)where c is the speed of light, and G is the gravitational constant. Evidently, the deflection angle is related to the integrated mass as $\mathit{\alpha }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathit{M}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{\pi }\hspace{0.17em}{\mathit{h}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{crit}}\mathit{x}}\mathrm{·}$(52)Introducing the central convergence, κ_c, a parameter that determines the lensing properties of the Einasto profile, ${\mathit{\kappa }}_{\mathrm{c}}\mathrm{\equiv }\frac{\mathrm{\Sigma }\mathrm{(}{\mathrm{0}}^{\mathrm{)}}}{{\mathrm{\Sigma }}_{\mathrm{crit}}}\mathrm{=}\frac{\mathrm{2}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}\mathit{n}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}{{\mathrm{\Sigma }}_{\mathrm{crit}}}\mathit{,}$(53)we can write $\mathit{\alpha }{}^{\mathrm{\right(}}\mathit{x}^{\mathrm{\left)}}$ in the form $\mathit{\alpha }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}\hspace{0.17em}\sqrt{\mathit{\pi }}}{\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}{\mathit{H}}_{\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{2}\mathit{,}\mathrm{1}}\left[\begin{array}{c}\\ \mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\mathit{,}\mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{\right)}\\ \mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{2}\mathit{n}\mathrm{\right)}\mathit{,}\mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\mathit{,}\mathrm{(}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\end{array}|\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\right]\mathit{,}$(54)with a completely analogous series expansion as for M(x). For a spherically symmetric lens that is capable of forming multiple images of the source, one sufficient condition is κ_c > 1 (Schneider et al. 1992). In the case κ_c ≤ 1, only one image of the source is formed. In addition to the condition κ_c > 1, multiple images are produced only if |y| ≤ y_crit (Li & Ostriker 2002), where y_crit is the maximum value of y when x < 0 or its minimum when x > 0. For singular profiles such as an NFW profile, the central convergence always is divergent, hence the condition κ_c > 1 is always met, this implies that an NFW profile is capable of forming multiple images for any mass. Nonsingular profiles, such as the Einasto profile, are not capable of forming multiple images for any mass. Instead, the condition κ_c > 1 sets a minimum value for lens mass required to form multiple images.

3.3. Deflection potential

The deflection potential $\mathit{\psi }{}^{\mathrm{\right(}}\mathit{x}^{\mathrm{\left)}}$ for a spherically symmetric lens is given by (Schneider et al. 1992) $\mathit{\psi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{x}}{\mathit{x}}^{\mathrm{\prime }}\hspace{0.17em}\mathit{\kappa }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\ln \left(\frac{\mathit{x}}{{\mathit{x}}^{\mathrm{\prime }}}\right)\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathrm{·}$(55)Inserting Eq. (38) into (55), we obtain again a result that can be re-expressed as a Fox H function $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}\mathit{x}\mathrm{\right)}& & \mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}\hspace{0.17em}\sqrt{\mathit{\pi }}}{\mathrm{2}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}{\mathit{x}}^{\mathrm{3}}\hspace{0.17em}\\ & & \mathrm{\times }\hspace{0.17em}{\mathit{H}}_{\mathrm{3}\mathit{,}\mathrm{4}}^{\mathrm{2}\mathit{,}\mathrm{2}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\right]\mathit{.}\end{array}$(56)For ψ(x) we have the same cases of expansions, a power series for simple poles in case 1, $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}\hspace{0.17em}\sqrt{\mathit{\pi }}}{\mathrm{2}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}[\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathit{k}\mathit{/}\mathit{n}\mathrm{+}\mathrm{3}}}{\mathrm{3}\mathit{n}\mathrm{+}\mathit{k}}& & \\ \mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathit{nk}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}\mathit{k}\mathrm{+}\mathrm{2}}}{\mathit{k}\mathrm{+}\mathrm{1}}]\mathit{,}& & \end{array}$(57)and a logarithmic-power series for multiple poles in case 2, $\begin{array}{ccc}& & \mathit{\psi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}\hspace{0.17em}\sqrt{\mathit{\pi }}}{\mathrm{2}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}[\mathrm{-}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{1}\mathit{k}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{p}\ne \mathrm{0}\end{array}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathrm{-}{\frac{\mathit{k}}{\mathrm{2}\mathit{n}}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathit{k}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathit{k}\mathit{/}\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{3}\mathit{n}\mathrm{+}\mathit{k}}\\ & & \mathrm{+}\frac{\mathrm{\Gamma }\mathrm{\left(}\mathit{n}\mathrm{\right)}}{\sqrt{\mathit{\pi }}}\mathrm{+}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{1}\mathrm{(}\mathit{k}\mathrm{+}{\mathit{k}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{q}\ne \mathrm{0}\end{array}\frac{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathit{nk}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{)}}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{k}}}{\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}\mathit{k}}}{\mathit{k}\mathrm{+}\mathrm{1}}]\\ & & \mathrm{+}\frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{2}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}\begin{array}{}\mathrm{\infty }\\ \sum \\ \mathit{k}\mathrm{=}\mathrm{0}\mathrm{(}\mathit{k}\mathrm{+}{\mathit{k}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{mod}\hspace{0.17em}\mathit{q}\mathrm{=}\mathrm{0}\end{array}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{p}}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{)}\mathrm{!}}{\mathrm{(}\mathrm{2}\mathit{nk}\mathrm{-}\mathit{n}\mathrm{)}\mathrm{!}\hspace{0.17em}\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}\hspace{0.17em}\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}{\left(\frac{\mathit{x}}{\mathrm{2}}\right)}^{\mathrm{2}\mathit{k}}\\ & & \mathrm{\times }[\mathrm{-}\ln \left(\frac{\mathit{x}}{\mathrm{2}}\right)\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{k}}\mathrm{+}\mathit{\psi }\mathrm{(}\mathit{k}\mathrm{+}\mathrm{2}\mathrm{)}\mathrm{+}\mathit{n}\hspace{0.17em}\mathit{\psi }\mathrm{(}\mathrm{2}\mathit{nk}\mathrm{-}\mathit{n}\mathrm{)}\mathrm{-}\mathit{\psi }\mathrm{(}\mathrm{2}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}]\mathit{.}\end{array}$(58)

3.4. Magnification, shear and the critical curves

Gravitational lensing effect preserves the surface brightness but it causes variations in shape and solid angle of the source. Therefore, the source luminosity is amplified by (Schneider et al. 1992) $\mathit{\mu }\mathrm{=}\frac{\mathrm{1}}{{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\kappa }}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}^{\mathrm{2}}}\mathit{,}$(59)where κ is the convergence and $\mathit{\gamma }\mathrm{=}\mathit{\gamma }{}^{\mathrm{\right(}}\mathit{x}^{\mathrm{\left)}}$ is the shear. The amplification μ has two contributions: one from the convergence, which describes an isotropic focusing of light rays in the lens plane, and the other which describes an anisotropic focusing caused by the tidal gravitational forces acting on the light rays, described by the shear. For a spherical symmetric lens, the shear is given by (Miralda-Escude 1991) $\mathit{\gamma }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{\mathrm{\Sigma ̅}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{-}\mathrm{\Sigma }\mathrm{(}{\mathit{x}}^{\mathrm{)}}}{{\mathrm{\Sigma }}_{\mathrm{crit}}}\mathrm{=}\mathit{\kappa ̅}\mathrm{-}\mathit{\kappa ,}$(60)where the average surface mass density within x is $\mathrm{\Sigma ̅}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{{\mathit{x}}^{\mathrm{2}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{x}}{\mathit{x}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{\Sigma }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{.}$(61)The magnification of the Einasto profile can be found combining Eqs. (38), (59)–(61). In the calculation of $\mathrm{\Sigma ̅}{}^{\mathrm{\right(}}\mathit{x}^{\mathrm{\left)}}$, we use again the Mellin-Barnes integral representation (36), and integrate it to obtain an expression in terms of Fox H function. Thus we get $\mathit{\mu }\mathrm{=}{\left[\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\kappa ̅}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\kappa ̅}\mathrm{-}\mathrm{2}{\mathit{\kappa }}^{\mathrm{)}}\right]}^{-1}\mathit{,}$(62)where $\mathit{\kappa }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}\hspace{0.17em}\sqrt{\mathit{\pi }}}{\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{H}}_{\mathrm{1}\mathit{,}\mathrm{2}}^{\mathrm{2}\mathit{,}\mathrm{0}}\left[\begin{array}{c}\\ \mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{\right)}\\ \mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{2}\mathit{n}\mathrm{\right)}\mathit{,}\mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\end{array}|\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\right]\mathit{,}$(63)and $\mathit{\kappa ̅}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}\hspace{0.17em}\sqrt{\mathit{\pi }}}{\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{H}}_{\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{2}\mathit{,}\mathrm{1}}\left[\begin{array}{c}\\ \mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\mathit{,}\mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{\right)}\\ \mathrm{\left(}\mathrm{0}\mathit{,}\hspace{0.17em}\mathrm{2}\mathit{n}\mathrm{\right)}\mathit{,}\mathrm{(}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\mathit{,}\mathrm{(}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{1}\mathrm{)}\end{array}|\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\right]\mathit{.}$(64)The magnification may be divergent for some image positions. The loci of the diverging magnification in the image plane are called the critical curves. We see from Eq. (62) that the Einasto profile has one pair of critical curves. The first curve, , is the tangential critical curve, which corresponds to an Einstein Ring with a radius, called the Einstein radius. The second curve, , is the radial critical curve, which also defines a ring and its corresponding radius. In both cases the equations must be solved numerically.

4. Integer and half-integer values of n

The expressions for the Einasto surface mass density (38) and its lensing properties (42), (54), (56), (62)–(64) in terms of the Fox H function for case n is a rational number, can be reduced to a Meijer G function, $\begin{array}{ccc}{\mathit{G}}_{\mathit{p,q}}^{\mathit{m,n}}\left[\begin{array}{c}\\ \\ \end{array}\hspace{0.17em}|\hspace{0.17em}\mathit{z}\right]& \mathrm{=}& \\ & & \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathrm{i}}{\mathrm{\int }}_{\mathrm{ℒ}}\frac{{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{m}}\mathrm{\Gamma }\mathrm{(}{\mathit{b}}_{\mathit{j}}\mathrm{+}\mathit{s}\mathrm{)}{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathit{j}}\mathrm{-}\mathit{s}\mathrm{)}}{{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathit{m}\mathrm{+}\mathrm{1}}^{\mathit{q}}\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{b}}_{\mathit{j}}\mathrm{-}\mathit{s}\mathrm{)}{\mathrm{􏽑}}_{\mathit{j}\mathrm{=}\mathit{n}\mathrm{+}\mathrm{1}}^{\mathit{p}}\mathrm{\Gamma }\mathrm{(}{\mathit{a}}_{\mathit{j}}\mathrm{+}\mathit{s}\mathrm{)}}\hspace{0.17em}{\mathit{z}}^{\mathrm{-}\mathit{s}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathit{.}\end{array}$(65)Meijer G function numerical routines had been implemented in computer algebraic systems, such as the commercial Maple, Mathematica and the free open-source Sage and mpmath library; in contrast, there is no Fox H function implementation available.

For case n integer or half-integer, we can write $\mathrm{\Sigma }{}^{\mathrm{\right(}}\mathit{x}^{\mathrm{\left)}}$ as a Meijer G function. Inserting the Gauss multiplication formula (Abramowitz & Stegun 1970) $\mathrm{\Gamma }\mathrm{\left(}\mathrm{2}\mathit{ny}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{2}\mathit{n}{\mathrm{)}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{ny}}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\mathit{n}}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathit{y}\mathrm{\right)}\hspace{0.17em}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{2}\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{\Gamma }\left(\frac{\mathit{j}}{\mathrm{2}\mathit{n}}\mathrm{+}\mathit{y}\right)\mathit{,}$(66)into Eq. (36) and comparing with the definition (65), we obtain an expression for the surface mass density of the Einasto profile in terms of the Meijer G function $\mathrm{\Sigma }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\sqrt{\mathit{n}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{G}}_{\mathrm{0}\mathit{,}\mathrm{2}\mathit{n}}^{\mathrm{2}\mathit{n,}\mathrm{0}}\left[\begin{array}{c}\\ \mathrm{-}\\ {b}\end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathrm{(}\mathrm{2}{\mathit{n}}^{\mathrm{)}}}^{\mathrm{2}\mathit{n}}}\right]\mathit{,}$(67a)where b is a vector of size 2n given by ${b}\mathrm{=}\{\frac{\mathrm{1}}{\mathrm{2}\mathit{n}}\mathit{,}\frac{\mathrm{2}}{\mathrm{2}\mathit{n}}\mathit{,}\mathit{...}\mathit{,}\frac{\mathrm{2}\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{2}\mathit{n}}\mathit{,}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\}\mathit{.}$(67b)Integrating the above expression for surface mass density according to (41) and using the integral properties of the Meijer G function (Eq. 07.34.21.0003.01 on the Wolfram Functions Site ¹), we obtain $\mathit{M}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{\sqrt{\mathit{n}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}}{\mathrm{2}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{2}}}\hspace{0.17em}{\mathit{x}}^{\mathrm{3}}\hspace{0.17em}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}^{\mathrm{2}\mathit{n,}\mathrm{1}}\left[\begin{array}{c}\\ \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\\ {b}\mathit{,}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathrm{(}\mathrm{2}{\mathit{n}}^{\mathrm{)}}}^{\mathrm{2}\mathit{n}}}\right]\mathrm{·}$(68)Inserting Eq. (65) into (49) and performing the integral of the Meijer G function, we find $\mathit{\alpha }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{2}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}\hspace{0.17em}\sqrt{\mathit{n}}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}^{\mathrm{2}\mathit{n,}\mathrm{1}}\left[\begin{array}{c}\\ \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\\ {b}\mathit{,}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathrm{(}\mathrm{2}{\mathit{n}}^{\mathrm{)}}}^{\mathrm{2}\mathit{n}}}\right]\mathrm{·}$(69)We obtain the deflection potential in terms of the Meijer G function combining Eqs. (55), (65) and again we integrate, applying the Meijer G function properties. We obtain $\mathit{\psi }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{4}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}\hspace{0.17em}\sqrt{\mathit{n}}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}{\mathit{x}}^{\mathrm{3}}\hspace{0.17em}{\mathit{G}}_{\mathrm{2}\mathit{,}\mathrm{2}\mathit{n}\mathrm{+}\mathrm{2}}^{\mathrm{2}\mathit{n,}\mathrm{2}}\left[\begin{array}{c}\\ \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\\ {b}\mathit{,}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{,}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathrm{(}\mathrm{2}{\mathit{n}}^{\mathrm{)}}}^{\mathrm{2}\mathit{n}}}\right]\mathrm{·}$(70)The convergence can be simply found dividing Eq. (65) by Σ_crit$\mathit{\kappa }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{2}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}\hspace{0.17em}\sqrt{\mathit{n}}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{G}}_{\mathrm{0}\mathit{,}\mathrm{2}\mathit{n}}^{\mathrm{2}\mathit{n,}\mathrm{0}}\left[\begin{array}{c}\\ \mathrm{-}\\ {b}\end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathrm{(}\mathrm{2}{\mathit{n}}^{\mathrm{)}}}^{\mathrm{2}\mathit{n}}}\right]\mathit{,}$(71)

the average convergence can be calculated by inserting Eq. (65) into Eq. (61), and then integrating, we find $\mathit{\kappa ̅}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{2}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}\hspace{0.17em}\sqrt{\mathit{n}}\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}{\mathit{n}}^{\mathrm{)}}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}^{\mathrm{2}\mathit{n,}\mathrm{1}}\left[\begin{array}{c}\\ \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\\ {b}\mathit{,}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathrm{(}\mathrm{2}{\mathit{n}}^{\mathrm{)}}}^{\mathrm{2}\mathit{n}}}\right]\mathrm{·}$(72)

4.1. Simple cases: n = 1 and n = $\frac{1}{2}$

The Einasto profile corresponding to n = 1 is a simple exponential model, characterised by the density profile $\mathit{\rho }\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\exp \left\{\mathrm{-}\left(\frac{\mathit{r}}{\mathit{h}}\right)\hspace{0.17em}\right\}\mathit{,}$(73)the surface density profile and lensing properties are readily found by setting n = 1 in Eqs. (65), (68), (69), (70), (71) and (72) $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{h}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{G}}_{\mathrm{0}\mathit{,}\mathrm{2}}^{\mathrm{2}\mathit{,}\mathrm{0}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{4}}\right]\mathit{,}\\ \mathit{M}\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \mathit{\pi }\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}{\mathit{x}}^{\mathrm{3}}\hspace{0.17em}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{3}}^{\mathrm{2}\mathit{,}\mathrm{1}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{4}}\right]\mathit{,}\\ \mathit{\alpha }\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{2}}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{3}}^{\mathrm{2}\mathit{,}\mathrm{1}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{4}}\right]\mathit{,}\\ \mathit{\psi }\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{4}}\hspace{0.17em}{\mathit{x}}^{\mathrm{3}}\hspace{0.17em}{\mathit{G}}_{\mathrm{2}\mathit{,}\mathrm{4}}^{\mathrm{2}\mathit{,}\mathrm{2}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{4}}\right]\mathit{,}\\ \mathit{\kappa }\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{2}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{G}}_{\mathrm{0}\mathit{,}\mathrm{2}}^{\mathrm{2}\mathit{,}\mathrm{0}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{4}}\right]\mathit{,}\\ \mathit{\kappa ̅}\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \frac{{\mathit{\kappa }}_{\mathrm{c}}}{\mathrm{2}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{3}}^{\mathrm{2}\mathit{,}\mathrm{1}}\left[\begin{array}{c}\\ \\ \end{array}|\hspace{0.17em}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{4}}\right]\mathrm{·}\end{array}$The specific Meijer G functions in these expressions can be reduced to more standard special functions. We obtain the equivalent expressions $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \\ \mathit{M}\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \mathrm{8}\mathit{\pi }\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{h}}^{\mathrm{3}}\hspace{0.17em}\left[\mathrm{1}\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}\hspace{0.17em}{\mathit{K}}_{\mathrm{2}}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\right]\mathit{,}\\ \mathit{\alpha }\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \frac{\mathrm{4}{\mathit{\kappa }}_{\mathrm{c}}}{\mathit{x}}\hspace{0.17em}\left[\mathrm{1}\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}\hspace{0.17em}{\mathit{K}}_{\mathrm{2}}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\right]\mathit{,}\\ \mathit{\psi }\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \mathrm{4}{\mathit{\kappa }}_{\mathrm{c}}\left[\ln \left(\frac{\mathit{x}}{\mathrm{2}}\right)\mathrm{+}\frac{\mathit{x}}{\mathrm{2}}\hspace{0.17em}{\mathit{K}}_{\mathrm{1}}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{+}{\mathit{K}}_{\mathrm{0}}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{+}\mathit{\gamma }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\right]\mathit{,}\\ \mathit{\kappa }\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& {\mathit{\kappa }}_{\mathrm{c}}\hspace{0.17em}\mathit{x}\hspace{0.17em}{\mathit{K}}_{\mathrm{1}}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathit{,}\\ \mathit{\kappa ̅}\mathrm{(}{\mathit{x}}^{\mathrm{)}}& \mathrm{=}& \frac{\mathrm{4}{\mathit{\kappa }}_{\mathrm{c}}}{{\mathit{x}}^{\mathrm{2}}}\hspace{0.17em}\left[\mathrm{1}\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}\hspace{0.17em}{\mathit{K}}_{\mathrm{2}}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\right]\mathit{,}\end{array}$where K_ν(x) is the modified Bessel function of the second kind of order ν, and γ ≈ 0.57721566 the Euler-Mascheroni constant. The expressions (80)–(85) can also be readily calculated by inserting the density profile (73) into the formula (24) and later carrying out all corresponding calculations for the lensing properties.