A precise analytical approximation for the deprojection of the Sérsic profile

The Sérsic model shows a close ﬁt to the surface brightness (or surface density) proﬁles of elliptical galaxies and galaxy bulges, and possibly also those of dwarf spheroidal galaxies and globular clusters. The deprojected density and mass proﬁles are important for many astrophysical applications, in particular for mass-orbit modeling of these systems. However, the exact deprojection formula for the Sérsic model employs special functions that are not available in most computer languages. We show that all previous analytical approximations to the 3D density proﬁle are imprecise at low Sérsic index ( n (cid:46) 1 . 5). We derived a more precise analytical approximation to the deprojected Sérsic density and mass proﬁles by ﬁtting two-dimensional tenth-order polynomials to the residuals of the analytical approximations by Lima Neto et al. (1999, MNRAS, 309, 481; LGM) for these proﬁles, relative to the numerical estimates. Our LGM-based polynomial ﬁts have typical relative precision better than 0.2% for both density and mass proﬁles, for Sérsic indices 0 . 5 ≤ n ≤ 10 and radii 0 . 001 < r / R e < 1000. Our approximation is much more precise than previously published approximations (except, in some models, for a few discrete values of the index). An appendix compares the deprojected Sérsic proﬁles with those of other popular simple models.


Introduction
The Sérsic model (Sérsic 1963;Sersic 1968) is the generalization of the R 1/4 law (de Vaucouleurs 1948) to describe the surface brightness profiles of elliptical galaxies (Caon et al. 1993) and the bulges of spiral galaxies (Andredakis et al. 1995). It has also been used to describe the surface density profiles of nuclear star clusters (Graham & Spitler 2009), resolved dwarf spheroidal galaxies (Battaglia et al. 2006), and globular clusters (Barmby et al. 2007).
The surface (mass or number) density (or equivalently surface brightness) of the Sérsic model is where R is the projected distance to the source center, R e is the effective radius containing half of the projected luminosity, n is the Sérsic index, and Σ 0 is the central surface density. The term b n is a function of n, obtained by solving the equation: where γ(a, x) = x 0 t a−1 e −t dt is the lower incomplete gamma function.
Since the Sérsic model accurately represents astronomical objects viewed in projection, it is important to know its corresponding three-dimensional (3D) density and mass profiles. These serve as a reference for comparison with other possible observational tracers, as well as to dark matter. Moreover, the 3D density profile is required for modeling the kinematics of spherical structures because it appears in the Jeans equation of local dynamical equilibrium. Since the Jeans equation also contains the total mass profile, the 3D mass profiles of stellar components are required to estimate the dark matter mass profile of elliptical and dwarf spheroidal galaxies.
Many authors assume that simple three-dimensional models resemble Sérsic models for certain values of the Sérsic index: It is often assumed that massive ellipticals and spiral bulges are well represented by the Hernquist (1990) model (e.g., Widrow & Dubinski 2005). On the other hand, dwarf spheroidal galaxies are often described with the Plummer (1911) model (e.g., Muñoz et al. 2018, who also tried Sérsic and other models), while ultra diffuse galaxies have been described with the Einasto (Einasto 1965;Navarro et al. 2004) 1 model (Nusser 2019). Łokas & Mamon (2001) noted that the projected Navarro et al. (1996, hereafter NFW) model resembles an n = 3 Sérsic for reasonable concentrations. Finally, n = 4 Sérsic models are considered to resemble the Jaffe (1983) model (Ciotti et al. 2019). In Appendix A, we compare these models to the deprojected Sérsic.
Unfortunately, the deprojection of the Sérsic surface density profile to a 3D (mass or number) 2 density profile, through Abel (1826) inversion 1 Navarro et al. (2004) showed how the Einasto model accurately represents the density profiles of dark matter halos in dissipationless cosmological simulations, while Merritt et al. (2005) were the first to note its similar form to the Sérsic model, and Merritt et al. (2006) were first to realize that this model had been previously introduced by Einasto. 2 The number profile always has the same form as the mass profile, and is obtained by simply replacing M(r) with N(r) and M ∞ with N ∞ , e.g., in Eqs. (7), (8), and (17).
(e.g., Binney & Mamon 1982), as well as the corresponding 3D mass (or number) profile both involve the complicated Meijer G special function (Mazure & Capelato 2002 for integer values of n, and Baes & Gentile 2011 for general values of n) or the other, complicated Fox H function (Baes & van Hese 2011), neither of which are available in popular computer languages. Following the shape of the analytical approximation to the R 1/4 law by Mellier & Mathez (1987), Prugniel & Simien (1997, hereafter, PS) proposed an analytical approximation for the 3D density of the Sérsic profile: where p n is a function depending only on n: Equation (5) yields a simple analytical form for the 3D mass profile, Lima Neto et al. (1999, hereafter LGM) later perfected the approximation of Eq. (6) with According to LGM, Eq. (9) has 5% relative accuracy for 0.56 ≤ n ≤ 10 and −2 < log(r/R e ) < 3. However, the power-law approximation at small radii is unjustified for small n. Indeed, as shown by Baes & Gentile (2011), the central density profile converges to a finite value for n < 1 (and the inner density profile diverges only logarithmically for n = 1), as we illustrate in Sect. 3. Simonneau & Prada (1999 proposed the quasi-Gaussian expansion for the density profile where where x j and w j are ten fit parameters. The individual SP density profiles (the terms inside the sum of Eq. (10)) have a similar (but not the same) form to that of the PS/LGM one, hence the similar shape of the mass profile:  Emsellem & van de Ven (2008) for the deprojected Sérsic density profile (filled circles). The solid and dotted curves show the spline cubic and linear interpolations, respectively. At small n, the parameters vary abruptly and the interpolations (both linear and cubic) are therefore uncertain. Trujillo et al. (2002) proposed an ellipsoid formula, which in the limit of spherical symmetry becomes where K ν (x) is the modified Bessel function of the second kind 3 of index ν, while ν n , p n , a n,0 , a n,1 , and a n,2 are empirical functions of index n. Trujillo et al. (2002) only provided their results for integer and half-integer values of n for n ≤ 5 and only integer values of n beyond. Emsellem & van de Ven (2008, hereafter EV) repeated their analysis on a finer grid of n, with steps of 0.1 for 0.5 ≤ n ≤ 1.5 and with one more term, a n,3 , in Eq. (14) involving 168 parameters. However, as shown in Fig. 1 -Both PS and LGM (and Márquez et al. 2000, which is the same as LGM, but with a slightly different last term for p n , which was a typo) are inappropriate for low n (Baes & Gentile 2011) and less precise than claimed (Emsellem & van de Ven 2008). -SP is limited to n ≥ 1, and is generally less precise than EV. -Trujillo et al. (2002) is only given for half-integer values of n and their parameters vary wildly with n for n ≤ 1.5. These authors do not provide a formula for the mass profile. -EV also suffers from discrete values of n, even though the grid is finer (∆n = 0.1 for n ≤ 1.5). These authors also did not provide a formula for the mass profile. In this article, we provide polynomial fits to the log residuals of the LGM approximation, which allow high accuracy to be reached for both the 3D density and 3D mass profiles in a wide range of Sérsic indices. In Sect. 2, we present the mathematical formalism and briefly explain our numerical integration method. We then show in Sect. 3 how our polynomial plus LGM approximation is orders of magnitude more precise than the formulae of LGM, SP, and Trujillo et al. (2002), as well as that of EV for low n, and is only slightly less precise for n 3. We conclude and discuss our results in Sect. 4.

Equations using dimensionless profiles
We express the general surface density, 3D density, and 3D mass (or number) profiles in terms of dimensionless functions: Hereafter, we use x = r/R e and X = R/R e . For the Sérsic model (see Graham & Driver 2005 for a thorough review of the Sérsic profile), the dimensionless surface density profile is while for the PS model, one can write the dimensionless 3D density and mass profiles as It is easy to show that the deprojection Eq.
(3) becomes The dimensionless mass profile is where Eq. (24) is obtained by inversion of the order of integration in the second equality of Eq. (23).
We performed our analysis using either the highly accurate approximations for b n of Ciotti & Bertin (1999, hereafter, CB) or the exact (numerical) solutions of Eq. (2). We noticed that the difference between these two approaches was negligible (see Sect. 3).
We then fit two-dimensional polynomials to both log ρ LGM (x, n)/ ρ(x, n) and log M LGM (x, n)/ M(x, n) , for geometrically spaced x and n, writing with polynomial orders 2 ≤ k ≤ 12. For this, we used Python's package numpy.linalg.lstsq. We found the smallest residuals for tenth-order polynomials when using both the b n approximation of CB and b n by numerically solving Eq.
(2). The coefficients are provided in Tables B.1-B.4. In the rest of the paper, we present the results relative to the CB approximation, since it is a simpler and more used model, and also because our tenthorder polynomial fits the exact b n case remarkably well.
2.3. Numerical precision: tests for known simple analytical deprojections (n = 0.5 and 1) For Sérsic indices n = 0.5 and n = 1, there are analytical solutions for the 3D density profile:  28), with green-colored titles) as a function of both Sérsic index (abscissae) and radius (ordinates). The color scale given in the vertical color bars is linear for log ratios between −0.001 and 0.001 and logarithmic beyond. The gray region and green curves in the upper left of the density panels are for regions where the numerical integration reached the underflow limit or density 10 −30 times ρ(R e ), respectively, because of the very rapid decline of density at large radii for low n, and also covers n < 1, which is not covered by the SP model. We note that the EV and Trujillo+02 models perform better at specific values of n that are often missed in our grid.
where K 0 (x) is the modified Bessel function of the second kind of index 0. We can therefore verify the numerical integration of Eq. (21) for these two Sérsic indices. For the interval −3 ≤ log(r/R e ) ≤ 3, we compared the densities from numerical integration with the analytical formulae of Eq. (27) using the CB approximation for b n . The match is very good, with root-mean-square (rms) values of log ( ρ ana / ρ num ) of 1.5 × 10 −7 and 2 × 10 −8 for n = 0.5 and n = 1, respectively. The same comparison using the exact b n yields 7 × 10 −5 and 2 × 10 −8 , respectively (with one particular value of r causing the higher rms for n = 0.5).

Results
As seen in Fig. 2, the 3D density profiles depart from the power laws proposed by LGM at low n, especially for low radii, as expected by the asymptotic expansions of Baes & Gentile (2011) for n < 1. Interestingly, the LGM formula is also inadequate at low radii for n = 1.25 and 2.25, although the asymptotic expansion of Baes & Gentile (2011) indicate power-law behavior at small radii. This poor accuracy of the LGM approximation at low radii is a serious concern when performing kinematic modeling of systems with possible central massive black holes. For example, Gaia second data release (DR2) positions and proper motions for stars in nearby globular clusters extend inwards to 0.7 arcsec from the center, which translates to 0.002 R e .
We now compare the accuracy of the different analytical approximations for the 3D density and 3D mass profiles. Figure 3 displays the ratio log f model / f , for f = ρ and f = M, for the main analytical approximations available in the literature, along with our new model where f is either the 3D density or 3D mass profile. We see that our model presents a more continuous behavior over the full range of Sérsic indices and radii. Our approximation displays the smallest residuals among all models for n 3 (except that SP outperforms our model for mass estimates at r > 3 R e for n > 1.3).
The variation of accuracy with Sérsic index can be seen in more detail in Fig. 4, which displays the rms of log f model / f , over the radial domain where ρ(r) > 10 −30 ρ(R e ), of the main analytical approximations using 1000 log-spaced Sérsic indices. Figure 4 indicates that the SP and EV approximations for density are less accurate than 2.3% (0.01 for log f model / f ) for n < 1.6 and n < 1.3, respectively. Our approximation (Eq. (28)) is more accurate than SP for n < 4.3 (density) and n < 3.1 (mass), and is more accurate than EV for n < 3.4 (density), except for their particular choices of n. Figure 4 shows that the EV approximation is much more accurate at specific values of n (we note that our grid does not contain all of these values precisely, and therefore the EV approximation is even more accurate at these specific values of n). However, these specific values of n represent a negligible measure compared to the full continuous range of 0.5 ≤ n ≤ 10. Therefore, the EV approximation at low n is not reliable for estimating the 3D density profile.
We analyzed the results shown in Fig. 4 using b n from CB or by numerically solving Eq. (2), and the results were very similar. In fact, the results are similar if we adopt one form of b n in the numerical integration and the other in the analytical approximations. This can be explained by the fact that log f LGM / f is practically the same for both estimates of b n , yielding a very similar fit of Eq. (26).
Finally, Table 1 provides the rms accuracies computed over the full range of radii −3 ≤ log(r/R e ) ≤ 3 and 0.5 ≤ n ≤ 10, except for the SP formula, which does not allow n ≤ 1, and also avoiding the domain where ρ(r) < 10 −30 ρ(R e ). We see that, averaging over all Sérsic indices, our approximation is much more accurate than all others (with over ten times lower rms).

Conclusions and discussion
The Sérsic model is usually considered to provide excellent fits to the surface density (or surface brightness) profiles of elliptical galaxies, spiral bulges, and even dwarf spheroidal galaxies and globular clusters. In the past, many authors have used simple analytical models to describe these systems, arguing that their models, once projected, resemble Sérsic models. It is more relevant to compare the physically meaningful three-dimensional density profiles of these simple models to the deprojected Sérsic model.
This comparison is made in Appendix A for the Plummer, Jaffe, Hernquist, Einasto, and NFW models. As seen in Fig. A.1, most of the simple models do not provide good fits to the deprojected Sérsic model, even for narrow ranges of the Sérsic index. The Plummer model requires a low index at small radii, but a much higher index at large radii, and the normalized density profile fits poorly at most radii. The Hernquist model resembles the n = 2.8 deprojected Sérsic model at low radii and the n = 5.7 Sérsic at large radii. The NFW models resemble the n = 2.8 deprojected Sérsic at low radii (consistent with the  Trujillo et al. (2002) 0.1496 - Emsellem & van de Ven (2008) 0.0382 -New 0.0005 0.0007 Hybrid-1 (optimized for n cut ) 0.0004 0.0005 Hybrid-2 (optimized for r cut ) 0.0004 0.0005 Notes. The rms accuracies are computed over the full range of radii −3 ≤ log(r/R e ) ≤ 3 (100 steps) and 0.5 ≤ n ≤ 10 (50 steps), except for the SP formula, which does not allow n ≤ 1, and also avoiding the domain where ρ(r) < 10 −30 ρ(R e ). Trujillo et al. (2002) and EV do not provide analytical mass profiles. The lower two rows display hybrid models, both with our new approximation ρ new for n < 3.4 and ρ EV for n ≥ 3.4. The first hybrid model has a mass profile M new for n < 3 and M SP for n ≥ 3, while in hybrid model 2, the mass profile is M new for r < R e and M LGM for r ≥ R e . The values in bold highlight the accuracies with our approximations.
similarity of the projected Sérsic with NFW discovered by Łokas & Mamon 2001), but have a shallower slope at large radii than even the shallowest (n = 8) deprojected Sérsic model. On the other hand, the Jaffe model resembles the n = 5.7 model at all radii. Moreover, as seen in Fig. A.2, the Einasto model provides a fair representation (rms difference of density profiles normalized to value at half-mass radius less than 0.1 dex) of the deprojected Sérsic model for n > 6.5. We reconsidered the different analytical deprojections of the Sérsic surface brightness (or surface density) profile. We found that the analytical approximations present in the literature do not show satisfying results when the Sérsic index is in the range 0.5 ≤ n 1.5 (apart from the specific values of n given by Emsellem & van de Ven 2008). In particular, the power-law times exponential density profile of Prugniel &Simien (1997) andLima Neto et al. (1999) fails to reproduce the inner density profiles for low n, even up to n = 2.25 despite the power-law behavior expected at small radii for n > 1 (Baes & Gentile 2011).
With a tenth-order two-dimensional polynomial fit, we propose a new analytical approximation (Eq. (28)) that is precise over the range log 0.5 ≤ log n ≤ 1, for −3 ≤ log (r/R e ) ≤ 3. Our approximation provides the highest precision when averaging over all values of Sérsic indices and radii (Table 1). While the approximations of Simonneau & Prada (1999 on one hand and of Emsellem & van de Ven (2008) on the other are more accurate than ours for n > 4.3 and 3.4, respectively, ours is more accurate at lower Sérsic indices. This is important for the study of astronomical sources with low Sérsic indexes, such as galaxy bulges, nuclear star clusters, dwarf spheroidal galaxies, and globular clusters. Moreover, our approximation of Eq. (28) to the density profile is sufficiently accurate for most scientific analyses for n > 3. Nevertheless, the user could use a hybrid approximation, combining either the Simonneau & Prada (2004) or Emsellem & van de Ven (2008) approximations for n ≥ 3.4 and ours for n < 3.4 (as shown in the last rows of Table 1). Finally, our analysis has the advantage of also providing a precise approximation for the mass profile, whereas no analytical expression can be derived from the density profile of Emsellem & van de Ven (2008). Our Python 3 codes are available at 4 along with coefficients of Tables B.1 and B.2.
These results will be useful in future mass-orbit modeling analyses of low-mass spherical systems, as we are preparing for globular clusters (Vitral & Mamon, in prep.).     In Eq. (A.2) for NFW, c = r max /a, where r max is the maximum allowed radius (because, contrary to all other models discussed here, the NFW model has logarithmically divergent mass). Also, for Einasto, P (−1) (a, y) is the inverse regularized lower incomplete gamma function, i.e., x = P (−1) (a, y) satisfies γ(a, x)/Γ(a) = y 5 . For Sérsic, the conversion was done by fitting a third-order polynomial and recovering the relation r h /R e = 3 j=0 a i log i n, where {a 0 , a 1 , a 2 , a 3 } = {1.32491, 0.0545396, −0.0286632, 0.0035086}. The Einasto model, which is the 3D analog of the Sérsic model, resembles the deprojected Sérsic model. Figure (26) and (28), for f = M and b n calculated from Ciotti & Bertin (1999) accurate approximation.    (26) and (28), for f = M and b n calculated from Eq. (2).