An Accurate Analytic Mass Model for Lensing Galaxies

By integrating the volume density profile ρ(r) along the line of sight, which is taken to be the Z direction here, we obtain the surface density profile,

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Sigma }}(R)=2{\displaystyle \int }_{0}^{\infty }\rho (r){dZ}\\ \ =\,\left\{\begin{array}{ll}{ \mathcal B }(\alpha ){\rho }_{c}{r}_{c}^{\alpha }{R}^{1-\alpha }-2{\rho }_{c}{r}_{c}\times & \\ \tilde{z}\left[F\left(\tfrac{\alpha }{2},1;\tfrac{3}{2};{\tilde{z}}^{2}\right)-F\left(\tfrac{{\alpha }_{c}}{2},1;\tfrac{3}{2};{\tilde{z}}^{2}\right)\right] & \mathrm{if}\ R\leqslant {r}_{c}\\ { \mathcal B }(\alpha ){\rho }_{c}{r}_{c}^{\alpha }{R}^{1-\alpha } & \mathrm{if}\ R\geqslant {r}_{c},\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

where ${r}^{2}={R}^{2}+{Z}^{2}$, F() is the Gauss hypergeometric function, $\tilde{z}=\sqrt{1-{R}^{2}/{r}_{c}^{2}}$, and

$$\begin{eqnarray}&&{ \mathcal B }(\alpha )=\mathrm{Beta}\left(\displaystyle \frac{1}{2},\displaystyle \frac{\alpha -1}{2}\right)=\sqrt{\pi }\displaystyle \frac{{\rm{\Gamma }}(\tfrac{\alpha -1}{2})}{{\rm{\Gamma }}(\tfrac{\alpha }{2})},\end{eqnarray} \tag{ 5 }$$

where Γ(x) and Beta(x, y) are the complete gamma function and the beta function, respectively. The total mass within the projected radius R is then

$$\begin{eqnarray}&&\begin{array}{l}{M}_{2{\rm{D}}}(R)=2\pi {\displaystyle \int }_{0}^{R}{\rm{\Sigma }}(R)R\ {dR}\\ \ =\,\left\{\begin{array}{ll}2\pi \tfrac{{ \mathcal B }(\alpha )}{3-\alpha }{\rho }_{c}{r}_{c}^{\alpha }{R}^{3-\alpha }+{m}_{0}+\tfrac{4\pi {\rho }_{c}}{3}{r}_{c}^{3}\times & \\ {\tilde{z}}^{3}\left[F\left(\tfrac{\alpha }{2},1;\tfrac{5}{2};{\tilde{z}}^{2}\right)-F\left(\tfrac{{\alpha }_{c}}{2},1;\tfrac{5}{2};{\tilde{z}}^{2}\right)\right] & \mathrm{if}\ R\leqslant {r}_{c}\\ 2\pi \tfrac{{ \mathcal B }(\alpha )}{3-\alpha }{\rho }_{c}{r}_{c}^{\alpha }{R}^{3-\alpha }+{m}_{0} & \mathrm{if}\ R\geqslant {r}_{c}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 6 }$$

In lensing analyses, what we care about is the convergence, which is the surface mass density scaled by the critical surface mass density,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{crit}}=\displaystyle \frac{{c}^{2}}{4\pi G}\displaystyle \frac{{D}_{s}}{{D}_{d}{D}_{{ds}}},\end{eqnarray} \tag{ 7 }$$

where D_s, D_d, and D_ds are the angular diameter distances from the observer to the background source and to the lens, and from the lens to the source, respectively. By further introducing a scale radius b, i.e., 

$$\begin{eqnarray}&&{b}^{\alpha -1}=\displaystyle \frac{{ \mathcal B }(\alpha )}{{{\rm{\Sigma }}}_{\mathrm{crit}}}\displaystyle \frac{2}{3-\alpha }{\rho }_{c}{r}_{c}^{\alpha },\end{eqnarray} \tag{ 8 }$$

we can then obtain the dimensionless convergence,

$$\begin{eqnarray}&&\begin{array}{l}\kappa (R)={\rm{\Sigma }}(R)/{{\rm{\Sigma }}}_{\mathrm{crit}}\\ \ =\,\left\{\begin{array}{ll}\tfrac{3-\alpha }{2}{\left(\tfrac{b}{R}\right)}^{\alpha -1}-\tfrac{3-\alpha }{{ \mathcal B }(\alpha )}{\left(\tfrac{b}{{r}_{c}}\right)}^{\alpha -1}\times & \\ \tilde{z}\left[F\left(\tfrac{\alpha }{2},1;\tfrac{3}{2};{\tilde{z}}^{2}\right)-F\left(\tfrac{{\alpha }_{c}}{2},1;\tfrac{3}{2};{\tilde{z}}^{2}\right)\right] & \mathrm{if}\ R\leqslant {r}_{c}\\ \tfrac{3-\alpha }{2}{\left(\tfrac{b}{R}\right)}^{\alpha -1} & \mathrm{if}\ R\geqslant {r}_{c}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 9 }$$

and the mean convergence within the radius R,

$$\begin{eqnarray}&&\begin{array}{l}\bar{\kappa }(R)=\displaystyle \frac{1}{\pi {R}^{2}}\displaystyle \frac{{M}_{2{\rm{D}}}(R)}{{{\rm{\Sigma }}}_{\mathrm{crit}}}\\ =\,\left\{\begin{array}{ll}{\left(\tfrac{b}{R}\right)}^{\alpha -1}+{\bar{\kappa }}_{0}+\tfrac{2}{3}\tfrac{3-\alpha }{{ \mathcal B }(\alpha )}{\left(\tfrac{b}{{r}_{c}}\right)}^{\alpha -1}{\left(\tfrac{{r}_{c}}{R}\right)}^{2}\times & \\ {\tilde{z}}^{3}\left[F\left(\tfrac{\alpha }{2},1;\tfrac{5}{2};{\tilde{z}}^{2}\right)-F\left(\tfrac{{\alpha }_{c}}{2},1;\tfrac{5}{2};{\tilde{z}}^{2}\right)\right] & \mathrm{if}\ R\leqslant {r}_{c}\\ {\left(\tfrac{b}{R}\right)}^{\alpha -1}+{\bar{\kappa }}_{0} & \mathrm{if}\ R\geqslant {r}_{c},\end{array}\right.\end{array}\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\bar{\kappa }}_{0}=\displaystyle \frac{1}{\pi {R}^{2}}\displaystyle \frac{{m}_{0}}{{{\rm{\Sigma }}}_{\mathrm{crit}}}=-\displaystyle \frac{2}{{ \mathcal B }(\alpha )}\displaystyle \frac{\alpha -{\alpha }_{c}}{3-{\alpha }_{c}}{\left(\displaystyle \frac{b}{{r}_{c}}\right)}^{\alpha -1}{\left(\displaystyle \frac{{r}_{c}}{R}\right)}^{2}.\end{eqnarray} \tag{ 11 }$$

From the above functions, we can see that the mass density profile considered here is actually a combination of a power-law mass distribution with a negative or positive mass distribution in the central region. Thus, the convergence can also be written as

$$\begin{eqnarray}&&\kappa (R)={\rm{\Sigma }}(R)/{{\rm{\Sigma }}}_{\mathrm{crit}}={\kappa }_{1}(R)+{\kappa }_{2}(R),\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{\kappa }_{1}(R)=\displaystyle \frac{3-\alpha }{2}{\left(\displaystyle \frac{b}{R}\right)}^{\alpha -1}\end{eqnarray} \tag{ 13 }$$

is the power-law part, and κ₂(R) denotes the second part in the central region,

$$\begin{eqnarray}&&\begin{array}{l}{\kappa }_{2}(R)=-\displaystyle \frac{3-\alpha }{{ \mathcal B }(\alpha )}{\left(\displaystyle \frac{b}{{r}_{c}}\right)}^{\alpha -1}\\ \times \,\left\{\begin{array}{ll}\tilde{z}\left[F\left(\tfrac{\alpha }{2},1;\tfrac{3}{2};{\tilde{z}}^{2}\right)-F\left(\tfrac{{\alpha }_{c}}{2},1;\tfrac{3}{2};{\tilde{z}}^{2}\right)\right] & \mathrm{if}\ R\leqslant {r}_{c}\\ 0 & \mathrm{if}\ R\geqslant {r}_{c}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 14 }$$

In Figure 1, we present several examples of the BPL convergence profiles with α = 2 and different values of α_c. As shown, the inner part of the density profile can take the shape of either a flat core or a steep cusp, depending on the value of α_c.

Zoom In Zoom Out Reset image size

More generally, in order to describe the elliptical surface mass distributions, we can generalize the circular radius R to an elliptical radius ${R}_{\mathrm{el}}=\sqrt{{{qx}}^{2}+{y}^{2}/q}$, where q is the axial ratio of the isodensity ellipses. This definition of elliptical radius conserves the area and the total mass within R_el.

The geometry of the deflection of a light ray can be described by the concise lens equation (Schneider et al. 2006; Tessore & Metcalf 2015),

$$\begin{eqnarray}&&{z}_{s}=z-\alpha (z),\end{eqnarray} \tag{ 15 }$$

which shows the mapping of a light ray from the position z = x + iy on the image plane to the position ${z}_{s}={x}_{s}+{\text{}}{{iy}}_{s}$ on the source plane, where $\alpha (z)={\alpha }_{x}+{\text{}}i{\alpha }_{y}$ is the scaled deflection angle and i denotes the imaginary unit. For an elliptical surface mass distribution, the corresponding complex conjugate of α(z) can be calculated using the formulation of BK75,

$$\begin{eqnarray}&&{\alpha }^{* }(z)=\displaystyle \frac{2}{z}{\int }_{0}^{{R}_{\mathrm{el}}}\displaystyle \frac{\kappa (R)R\ {dR}}{\sqrt{1-{\zeta }^{2}{R}^{2}}},\end{eqnarray} \tag{ 16 }$$

with ${\zeta }^{2}=(1/q-q)/{z}^{2}$. In the special case of q = 1, we find

$$\begin{eqnarray}&&{\alpha }^{* }(z)=\displaystyle \frac{{R}^{2}\bar{\kappa }(R)}{z}=\bar{\kappa }(R){z}^{* },\end{eqnarray} \tag{ 17 }$$

where z* is the complex conjugate of z.

Inserting Equation (12) into Equation (16), we then obtain the deflection field for the BPL model, i.e., 

$$\begin{eqnarray}&&{\alpha }^{* }(z)={\alpha }_{1}^{* }(z)+{\alpha }_{2}^{* }(z),\end{eqnarray} \tag{ 18 }$$

with

$$\begin{eqnarray}&&{\alpha }_{1}^{* }(z)=\displaystyle \frac{{R}_{\mathrm{el}}^{2}}{z}{\left(\displaystyle \frac{b}{{R}_{\mathrm{el}}}\right)}^{\alpha -1}F\left(\displaystyle \frac{1}{2},\displaystyle \frac{3-\alpha }{2};\displaystyle \frac{5-\alpha }{2};{\zeta }^{2}{R}_{\mathrm{el}}^{2}\right)\end{eqnarray} \tag{ 19 }$$

for the SPL mass distribution κ₁ and

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{2}^{* }(z) & = & \displaystyle \frac{{r}_{c}^{2}}{z}\displaystyle \frac{3-\alpha }{{ \mathcal B }(\alpha )}{\left(\displaystyle \frac{b}{{r}_{c}}\right)}^{\alpha -1}\left[\displaystyle \frac{2}{3-{\alpha }_{c}}{\mathscr{F}}\left(\displaystyle \frac{3-{\alpha }_{c}}{2},{ \mathcal C }\right)\right.\\ & & -\,\left.\displaystyle \frac{2}{3-\alpha }{\mathscr{F}}\left(\displaystyle \frac{3-\alpha }{2},{ \mathcal C }\right)-{{ \mathcal S }}_{0}\right]\end{array}\end{eqnarray} \tag{ 20 }$$

for the complementary part κ₂, where

$$\begin{eqnarray}&&\begin{array}{l}{\mathscr{F}}(a,z)={}_{3}{F}_{2}\left(a,\displaystyle \frac{1}{2},1;a+1,\displaystyle \frac{3}{2};z\right)\\ =\,\left\{\begin{array}{ll}\tfrac{1}{1-2a}\left[F(a,1;a+1;z)-2{aF}\left(\tfrac{1}{2},1;\tfrac{3}{2};z\right)\right] & \mathrm{if}\ a\ne \tfrac{1}{2}\\ \tfrac{{\mathrm{Li}}_{2}(\sqrt{z})-{\mathrm{Li}}_{2}(-\sqrt{z})}{2\sqrt{z}} & \mathrm{if}\ a=\tfrac{1}{2}\end{array}\right.\end{array}\end{eqnarray} \tag{ 21 }$$

is a special case of the generalized hypergeometric function ₃F₂(), where it is reduced to the Gauss hypergeometric function F() here, and Li₂() denotes Spence's function, i.e., ${\mathrm{Li}}_{2}(z)\,=-{\int }_{0}^{z}\tfrac{\mathrm{ln}(1-t){dt}}{t}$, and

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal S }}_{0}(\alpha ,{\alpha }_{c},{\tilde{z}}_{\mathrm{el}},{ \mathcal C })\\ =\,\left\{\begin{array}{ll}\tfrac{1}{\sqrt{1-{ \mathcal C }}}\sum _{n=0}^{\infty }\tfrac{{\left(\tfrac{{\alpha }_{c}}{2}\right)}^{(n)}-{\left(\tfrac{\alpha }{2}\right)}^{(n)}}{{\left(\tfrac{3}{2}\right)}^{(n)}}\tfrac{2{\tilde{z}}_{\mathrm{el}}^{2n+3}}{2n+3}\times & \\ F\left(\tfrac{1}{2},\tfrac{2n+3}{2},\tfrac{2n+5}{2},\tfrac{{ \mathcal C }{\tilde{z}}_{\mathrm{el}}^{2}}{{ \mathcal C }-1}\right) & \mathrm{if}\ {R}_{\mathrm{el}}\leqslant {r}_{c}\\ 0 & \mathrm{if}\ {R}_{\mathrm{el}}\geqslant {r}_{c}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 22 }$$

where ${ \mathcal C }={r}_{c}^{2}{\zeta }^{2}$, ${\tilde{z}}_{\mathrm{el}}=\sqrt{1-{R}_{\mathrm{el}}^{2}/{r}_{c}^{2}}$, and x⁽ⁿ⁾ denotes the rising factorial of x.

The series ${{ \mathcal S }}_{0}$ converges rapidly for pixels close to the break radius, due to the fact that ${\tilde{z}}_{\mathrm{el}}\to 0$ as R_el → r_c, but the convergence becomes slower in the very central region where ${\tilde{z}}_{\mathrm{el}}\to 1$. This is, however, not an issue, because there are only a few pixels in the central region of a lens. On the other hand, for real observations, the very central region usually suffers from larger noise than other regions because of the light contamination of the foreground lens. Therefore, in actual data analysis, the very central region of a lens could be masked to speed up the evaluation of deflection angles if there is no clearly visible central lensed image.

Given the analytical deflection field, it is straightforward to derive the lensing shear,

$$\begin{eqnarray}&&{\gamma }^{* }(z)=\displaystyle \frac{\partial {\alpha }^{* }}{\partial z},\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray*}&&\displaystyle \frac{\partial }{\partial z}=\displaystyle \frac{1}{2}\left[\displaystyle \frac{\partial }{\partial x}-i\displaystyle \frac{\partial }{\partial y}\right]\end{eqnarray*}$$

is the Wirtinger derivative (Wirtinger 1926; Kormann et al. 1994).

For the ${\kappa }_{1}$ part of the BPL model, as shown in the paper of Tessore & Metcalf (2015), the conjugate of the shear is

$$\begin{eqnarray}&&{\gamma }_{1}^{* }(z)=\displaystyle \frac{\partial {\alpha }_{1}^{* }}{\partial z}=(2-\alpha )\displaystyle \frac{{\alpha }_{1}^{* }}{z}-{\kappa }_{1}(z)\displaystyle \frac{{z}^{* }}{z}.\end{eqnarray} \tag{ 24 }$$

Similarly, the shear for the second part κ₂ is

$$\begin{eqnarray}&&\begin{array}{l}{\gamma }_{2}^{* }(z)=\displaystyle \frac{\partial {\alpha }_{2}^{* }}{\partial z}\\ \,=\,-{\kappa }_{2}(z)\displaystyle \frac{q| z{| }^{2}-(1+{q}^{2}){r}_{c}^{2}}{{{qz}}^{2}-(1-{q}^{2}){r}_{c}^{2}}+2\displaystyle \frac{{r}_{c}^{2}}{{z}^{2}}\displaystyle \frac{3-\alpha }{{ \mathcal B }(\alpha )}{\left(\displaystyle \frac{b}{{r}_{c}}\right)}^{\alpha -1}\\ \times \,\left[\displaystyle \frac{2-{\alpha }_{c}}{3-{\alpha }_{c}}{\mathscr{F}}\left(\displaystyle \frac{3-{\alpha }_{c}}{2},{ \mathcal C }\right)-\displaystyle \frac{2-\alpha }{3-\alpha }{\mathscr{F}}\left(\displaystyle \frac{3-\alpha }{2},{ \mathcal C }\right)-{{ \mathcal S }}_{2}\right],\end{array}\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal S }}_{2}(\alpha ,{\alpha }_{c},{\tilde{z}}_{\mathrm{el}},{ \mathcal C })\\ =\,\left\{\begin{array}{ll}\tfrac{1}{{\left(1-{ \mathcal C }\right)}^{\tfrac{3}{2}}}\sum _{n=0}^{\infty }\tfrac{{\left(\tfrac{{\alpha }_{c}}{2}\right)}^{(n)}-{\left(\tfrac{\alpha }{2}\right)}^{(n)}}{{\left(\tfrac{3}{2}\right)}^{(n)}}\tfrac{2n+2}{2n+3}\times & \\ {\tilde{z}}_{\mathrm{el}}^{2n+3}F\left(\tfrac{1}{2},\tfrac{2n+3}{2},\tfrac{2n+5}{2},\tfrac{{ \mathcal C }{\tilde{z}}_{\mathrm{el}}^{2}}{{ \mathcal C }-1}\right) & \mathrm{if}\ {R}_{\mathrm{el}}\leqslant {r}_{c}\\ 0 & \mathrm{if}\ {R}_{\mathrm{el}}\geqslant {r}_{c}\end{array}\right..\end{array}\end{eqnarray} \tag{ 26 }$$

Note that the series ${{ \mathcal S }}_{2}$ here also suffers from the same convergence problem as the series ${{ \mathcal S }}_{0}$ in the very central region, where the calculations should be done with caution.

If there exists a central black hole, its effect on the shear field can be added. Based on Equations (17) and (23), the shear induced by a central black hole is

$$\begin{eqnarray}&&{\gamma }_{b}^{* }=-\displaystyle \frac{1}{\pi }\displaystyle \frac{{m}_{b}}{{{\rm{\Sigma }}}_{\mathrm{crit}}}\displaystyle \frac{1}{{z}^{2}},\end{eqnarray} \tag{ 27 }$$

where m_b is the mass of the central black hole.

The magnification μ for the BPL density profile with a central black hole can be calculated according to

$$\begin{eqnarray}&&{\mu }^{-1}={\left(1-{\kappa }_{1}-{\kappa }_{2}\right)}^{2}-| {\gamma }_{1}^{* }+{\gamma }_{2}^{* }+{\gamma }_{b}^{* }{| }^{2}.\end{eqnarray} \tag{ 28 }$$

By solving equation μ⁻¹ = 0, the critical curves can then be found. The corresponding caustics can be obtained by mapping the critical curves on the lens plane to the source plane using the lens equation.

In Figure 2, we present examples of the critical curves and caustics for several cases of BPL density profiles with b = 15, α = 2, and α_c = 0.6 but different ellipticity q, break radius r_c, or black hole mass m_b. The top and bottom panels are for the spherical (q = 1) and elliptical cases (q = 0.7), respectively. From left to right, the values of the core radius r_c and black hole mass m_b are shown for each case in the top-left corner. In order to demonstrate the effect of a black hole in units of ${M}_{\odot }$, all of the lenses are assumed to be at redshift z_d ≃ 0.178 with a source at z_s = 0.6.

Zoom In Zoom Out Reset image size

For the spherical cases, it is shown in the first panel that there is only one critical curve for the singular isothermal sphere density profile, i.e., the spherical power-law profile with α = 2. The second panel demonstrates that the radial critical curve will appear if there is a fairly flat core. Furthermore, the inclusion of a central massive black hole will produce the third critical curve in the region much closer to the center. However, if the central black hole is too massive, all of the inner critical curves will disappear and only the outmost tangential critical curve remains (see examples in Mao et al. 2001).

For the elliptical cases, the same number of critical curves are presented as the corresponding spherical cases. However, the critical curves become elliptical-like, and the caustics change accordingly. What is obvious is that the caustics for the outmost tangential critical curves are turned into astroid-like curves from a point or a circle.

The critical curves and caustics provide references to speculate on the image configurations. However, we know that the pattern of lensed images is not only sensitive to the lensing mass distributions but also the source properties. In Figures B.1–B.3 in B, we illustrate the complex dependence of the lensed images on the lens and source properties.

In addition to the lensing observations, stellar kinematics can provide complementary information to the mass distributions of galaxies, especially the 2D kinematics from the integral-field spectroscopy (IFS; Bundy et al. 2015; Cappellari 2016). However, detailed IFS observations are currently only available for galaxies in the local universe whereas most of the SL systems are at redshift higher than 0.1 (Bolton et al. 2008; Shu et al. 2015). For higher redshift galaxies, the most efficient way now is to measure the 1D velocity dispersions using the optical single-fiber spectroscopy. In this subsection, we investigate the velocity dispersions in detail based on the BPL model for the lensing mass.

We assume that the galaxies are spherical in the modeling of the dynamics of galaxies. According to the spherical Jeans equation and assuming a constant velocity anisotropy parameter β, the radial velocity dispersion for the stars is (Binney & Tremaine 2008)

$$\begin{eqnarray}&&{\sigma }_{r}^{2}(r)=\displaystyle \frac{G{\int }_{r}^{\infty }{{dr}}^{{\prime} }j({r}^{{\prime} })M({r}^{{\prime} }){\left({r}^{{\prime} }\right)}^{2\beta -2}}{{r}^{2\beta }j(r)},\end{eqnarray} \tag{ 29 }$$

where j(r) is the 3D luminosity density profile related to the stellar number density profile and M(r) corresponds to the total mass within the 3D radius r. The LOSVD at the projected radius R is then given by

$$\begin{eqnarray}&&{\sigma }_{\parallel }^{2}(R)=\displaystyle \frac{{\int }_{-\infty }^{\infty }{dZ}\,j(r)(1-\beta {R}^{2}{/r}^{2}){\sigma }_{r}^{2}(r)}{{\int }_{-\infty }^{\infty }{dZ}\,j(r)}.\end{eqnarray} \tag{ 30 }$$

For single-fiber spectroscopic observations based on ground-based telescopes, the effect of fiber aperture and atmospheric seeing should be taken into account. In this case, the observed velocity dispersion can be modeled as

$$\begin{eqnarray}&&\langle {\sigma }_{\parallel }^{2}\rangle =\displaystyle \frac{{\int }_{0}^{\infty }{dR}\,{Rw}(R)I(R){\sigma }_{\parallel }^{2}(R)}{{\int }_{0}^{\infty }{dR}\,{Rw}(R)I(R)},\end{eqnarray} \tag{ 31 }$$

where w(R) is the weighting function accounting for the fiber aperture size and the seeing effect, and

$$\begin{eqnarray}&&I(R)={\int }_{-\infty }^{\infty }{dZ}\,j(r)\end{eqnarray} \tag{ 32 }$$

is the surface brightness distribution. Rewriting Equation (36) by inserting Equations (30) and (32), one gets

$$\begin{eqnarray}&&\langle {\sigma }_{\parallel }^{2}\rangle =\displaystyle \frac{{\int }_{0}^{\infty }{dR}\,{Rw}(R){\int }_{-\infty }^{\infty }{dZ}\,j(r)(1-\beta {R}^{2}{/r}^{2}){\sigma }_{r}^{2}(r)}{{\int }_{0}^{\infty }{dR}\,{Rw}(R){\int }_{-\infty }^{\infty }{dZ}\,j(r)},\end{eqnarray} \tag{ 33 }$$

which is hereafter named as the aperture and luminosity (AL)-weighted LOSVD.

In practice, as discussed in Schwab et al. (2010), we can use a Gaussian smoothing function to approximate the weighting function w(R) to some extent by assuming

$$\begin{eqnarray}&&w(R)\approx \exp \left(-\displaystyle \frac{{R}^{2}}{2{\sigma }_{\mathrm{fib}}^{2}}\right),\end{eqnarray} \tag{ 34 }$$

with

$$\begin{eqnarray}&&{\sigma }_{\mathrm{fib}}\approx {\sigma }_{\mathrm{see}}\sqrt{1+{\chi }^{2}/4+{\chi }^{4}/40}\end{eqnarray} \tag{ 35 }$$

and $\chi ={R}_{\mathrm{fib}}/{\sigma }_{\mathrm{see}}$, where ${\sigma }_{\mathrm{see}}$ is the Gaussian standard deviation of the seeing function equivalent to the FWHM of the seeing divided by $2\sqrt{2\mathrm{ln}2}$.

If we adopt a Gaussian form of w(R), i.e., Equation (34), by changing the order of integrations, Equation (33) can be reduced to

$$\begin{eqnarray}&&\begin{array}{l}\langle {\sigma }_{\parallel }^{2}\rangle \\ =\,\displaystyle \frac{{\displaystyle \int }_{0}^{\infty }{dr}\,{r}^{2}j(r){\sigma }_{r}^{2}(r)\left[{\rm{\Phi }}\left(1,\tfrac{3}{2};-\tfrac{{r}^{2}}{2{\sigma }_{\mathrm{fib}}^{2}}\right)-\tfrac{2\beta }{3}{\rm{\Phi }}\left(2,\tfrac{5}{2};-\tfrac{{r}^{2}}{2{\sigma }_{\mathrm{fib}}^{2}}\right)\right]}{{\displaystyle \int }_{0}^{\infty }{dr}\,{r}^{2}j(r){\rm{\Phi }}\left(1,\tfrac{3}{2};-\tfrac{{r}^{2}}{2{\sigma }_{\mathrm{fib}}^{2}}\right)}\end{array}\end{eqnarray} \tag{ 36 }$$

with only 1D integrals, where the function ${\rm{\Phi }}({a}_{1},{a}_{2};-x)\,={{\rm{e}}}^{-x}{\rm{\Phi }}({a}_{2}-{a}_{1},{a}_{2};x)$ being Kummer's confluent hypergeometric function.

Looking into Equation (36), we notice that both j(r) and σ_fib can be treated as known quantities: j(r) can be inferred from fitting to the surface brightness distribution, and σ_fib can be estimated from Equation (35). The AL-weighted LOSVD $\langle {\sigma }_{\parallel }^{2}\rangle$ is a direct observable. The mass distribution that we are paying attention to is implicit in the radial velocity dispersion ${\sigma }_{r}^{2}(r)$.

To stay consistent with the BPL lensing mass model, we adopt the Sérsic profile with a power-law inner density profile to describe the luminosity density profile of galaxies. It is the power-law Sérsic (PL-Sérsic) profile developed by Terzić & Graham (2005). The 3D form of the PL-Sérsic profile is written as

$$\begin{eqnarray}j(r)=\left\{\begin{array}{ll}{j}_{c}{\left(r/{r}_{c}\right)}^{-{\alpha }_{c}} & \mathrm{if}\ r\leqslant {r}_{c}\\ {j}_{0}{\left(\displaystyle \frac{r}{s}\right)}^{-u}\exp \left[-{\left(\displaystyle \frac{r}{s}\right)}^{\nu }\right] & \mathrm{if}\ r\geqslant {r}_{c}\end{array}\right.,\end{eqnarray} \tag{ 37 }$$

where

$$\begin{eqnarray}&&{j}_{0}={j}_{c}{\left(\displaystyle \frac{{r}_{c}}{s}\right)}^{u}\exp \left[{\left(\displaystyle \frac{{r}_{c}}{s}\right)}^{\nu }\right],\end{eqnarray} \tag{ 38 }$$

j_c is the luminosity density at r_c, $s={R}_{\mathrm{eff}}/{k}^{n}$ is a scale radius defined by the 2D effective radius R_eff for the single Sérsic profile and the Sérsic index n (note that k here is a function of n and its expression can be found in Ciotti & Bertin 1999 and MacArthur et al. 2003), ν = 1/n, and u = 1 − 0.6097ν + 0.054635ν² (Lima Neto et al. 1999; Márquez et al. 2001).

The surface luminosity density profile corresponding to the 3D PL-Sérsic profile is thus

$$\begin{eqnarray}&&\begin{array}{l}I(R)=2{\displaystyle \int }_{R}^{\infty }\displaystyle \frac{j(r)r\ {dr}}{\sqrt{{r}^{2}-{R}^{2}}}\\ =\,\left\{\begin{array}{ll}2{j}_{c}{r}_{c}\tilde{z}F\left(\tfrac{{\alpha }_{c}}{2},1;\tfrac{3}{2};{\tilde{z}}^{2}\right)+2{\displaystyle \int }_{{r}_{c}}^{\infty }\tfrac{j(r){rdr}}{\sqrt{{r}^{2}-{R}^{2}}} & \mathrm{if}\ R\leqslant {r}_{c}\\ {I}_{0}\exp \left[-{\left(\tfrac{R}{s}\right)}^{\nu }\right] & \mathrm{if}\ R\geqslant {r}_{c}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}&&{I}_{0}=2{{sj}}_{0}\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{3-u}{\nu }\right)}{{\rm{\Gamma }}\left(\tfrac{2}{\nu }\right)}.\end{eqnarray} \tag{ 40 }$$

In the following analyses, we simply assume that the BPL mass profile and the PL-Sérsic light profile have the same break radius r_c and the inner density slope α_c. This assumption simplifies the velocity dispersion calculations.

Inserting Equations (2) and (37) into Equation (29), we then derive an analytical form of the radial velocity dispersion for r ≥ r_c, i.e., 

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{r}^{2}(r) & = & G\displaystyle \frac{{r}^{-2\beta }}{J(r)}\left\{\displaystyle \frac{4\pi {\rho }_{c}{r}_{c}^{\alpha }}{3-\alpha }\displaystyle \frac{{s}^{\eta }}{\nu }{\rm{\Gamma }}\left[\displaystyle \frac{\eta }{\nu },{\left(\displaystyle \frac{r}{s}\right)}^{\nu }\right]\right.\\ & & +\,\left.{m}_{0}\displaystyle \frac{{s}^{\lambda }}{\nu }{\rm{\Gamma }}\left[\displaystyle \frac{\lambda }{\nu },{\left(\displaystyle \frac{r}{s}\right)}^{\nu }\right]\right\}\\ & = & { \mathcal A }\displaystyle \frac{{r}^{-2\beta }}{J(r)}\left\{\displaystyle \frac{{s}^{\eta }}{\nu }{\rm{\Gamma }}\left[\displaystyle \frac{\eta }{\nu },{\left(\displaystyle \frac{r}{s}\right)}^{\nu }\right]\right.\\ & & -\,\left.{r}_{c}^{3-\alpha }\displaystyle \frac{\alpha -{\alpha }_{c}}{3-{\alpha }_{c}}\displaystyle \frac{{s}^{\lambda }}{\nu }{\rm{\Gamma }}\left[\displaystyle \frac{\lambda }{\nu },{\left(\displaystyle \frac{r}{s}\right)}^{\nu }\right]\right\},\end{array}\end{eqnarray} \tag{ 41 }$$

where $J(r)={r}^{-u}\exp \left[-{\left(\tfrac{r}{s}\right)}^{\nu }\right]$, $\eta =2-u-\alpha +2\beta$, $\lambda \,=-1-u+2\beta$, and

$$\begin{eqnarray}&&{ \mathcal A }=\displaystyle \frac{4\pi G{\rho }_{c}{r}_{c}^{\alpha }}{3-\alpha }=\displaystyle \frac{{c}^{2}}{2}\displaystyle \frac{{D}_{s}}{{D}_{d}{D}_{{ds}}}\displaystyle \frac{{b}^{\alpha -1}}{{ \mathcal B }(\alpha )}\end{eqnarray} \tag{ 42 }$$

showing the relation between the mass density profile and the lensing-related quantities. For r < r_c, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{r}^{2}(r) & = & \displaystyle \frac{4\pi G{\rho }_{c}{r}_{c}^{{\alpha }_{c}}}{3-{\alpha }_{c}}\displaystyle \frac{{r}_{c}^{\mu }-{r}^{\mu }}{\mu }{r}^{{\alpha }_{c}-2\beta }+{\sigma }_{r}^{2}({r}_{c}){\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{{\alpha }_{c}-2\beta }\\ & = & { \mathcal A }\displaystyle \frac{3-\alpha }{3-{\alpha }_{c}}\displaystyle \frac{{r}_{c}^{\mu }-{r}^{\mu }}{\mu }\displaystyle \frac{{r}^{{\alpha }_{c}-2\beta }}{{r}_{c}^{\alpha -{\alpha }_{c}}}+{\sigma }_{r}^{2}({r}_{c}){\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{{\alpha }_{c}-2\beta },\end{array}\end{eqnarray} \tag{ 43 }$$

where μ = 2 −2α_c + 2β. Note that ${\sigma }_{r}^{2}(r)$ is still finite even if μ → 0, due to the fact that ${\mathrm{lim}}_{\mu \to 0}\tfrac{{r}_{c}^{\mu }-{r}^{\mu }}{\mu }=\mathrm{ln}\tfrac{{r}_{c}}{r}$.

When r_c = 0, the mass distribution has only the SPL part, and the light distribution follows the pure Sérsic profile. In this case, the radial velocity dispersion becomes

$$\begin{eqnarray}&&{\sigma }_{r}^{2}(r)={ \mathcal A }\displaystyle \frac{{r}^{-2\beta }}{J(r)}\displaystyle \frac{{s}^{\eta }}{\nu }{\rm{\Gamma }}\left[\displaystyle \frac{\eta }{\nu },{\left(\displaystyle \frac{r}{s}\right)}^{\nu }\right].\end{eqnarray} \tag{ 44 }$$

Generally, a black hole can exist at the center of a massive galaxy. In this case, the black hole can be regarded as a point mass which may have detectable effects on the velocity dispersion. With the PL-Sérsic light profile considered here, the contribution of a black hole to the radial velocity dispersion is given by

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{{\rm{b}},r}^{2}(r)\\ =\,\left\{\begin{array}{ll}{{Gm}}_{{\rm{b}}}\tfrac{{r}_{c}^{{\lambda }^{{\prime} }}-{r}^{{\lambda }^{{\prime} }}}{{\lambda }^{{\prime} }}{r}^{{\alpha }_{c}-2\beta }+{\sigma }_{{\rm{b}},r}^{2}({r}_{c}){\left(\tfrac{r}{{r}_{c}}\right)}^{{\alpha }_{c}-2\beta } & \mathrm{if}\ r\leqslant {r}_{c}\\ {{Gm}}_{{\rm{b}}}\tfrac{{r}^{-2\beta }}{J(r)}\tfrac{{s}^{\lambda }}{\nu }{\rm{\Gamma }}\left[\tfrac{\lambda }{\nu },{\left(\tfrac{r}{s}\right)}^{\nu }\right] & \mathrm{if}\ r\geqslant {r}_{c}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 45 }$$

where m_b denotes the mass of the central black hole and ${\lambda }^{{\prime} }=-1-{\alpha }_{c}+2\beta$.

The Illustris project consists of a series of large-scale hydrodynamic simulations, which incorporates various kinds of baryonic physics including gas cooling; star formation and evolution; feedback from active galactic nuclei, supernovae, and supermassive black holes; and so on (Genel et al. 2014; Vogelsberger et al. 2014; Nelson et al. 2015). We adopt in this work the highest resolution run, named the Illustris-1 simulation, to generate our mock catalogs.

The Illustris-1 simulation follows the dynamics of 1820³ dark matter particles, with 1820³ hydrodynamical cells initially in a periodic box with 75${h}^{-1}\,\mathrm{Mpc}$ a side. The mass resolution for dark matter particles is $\sim 6.26\times {10}^{6}\,{M}_{\odot }$, and the baryonic matter has an initial mass resolution of $\sim 1.26\times {10}^{6}\,{M}_{\odot }$. Different gravitational softening lengths are applied to different types of particles. For dark matter particles, the softening length is fixed to be a comoving value of 1 ${h}^{-1}\,\mathrm{kpc}$. For stars and black holes, the softening length is limited to a maximum value of 0.5 ${h}^{-1}\,\mathrm{kpc}$ in physical scale. For the gas cells, an adaptive softening length is defined according to the fiducial cell size and a floor given by the collisionless baryonic particles.

Focusing on the galaxies with stellar mass larger than ${10}^{10}{h}^{-1}\,{M}_{\odot }$, we finally extract 5343 galaxies in Illustris-1 at redshift zero. In order to generate mock galaxies with redshift consistent with the current observed strong lenses, we artificially put the galaxies at redshift z_d ≃ 0.178, which is close to the median redshift of lenses identified by the Sloan Lens ACS (SLACS) Survey (Bolton et al. 2008; Shu et al. 2015). For the lens redshift z_d ≃ 0.178, 1'' corresponds to $3\,\mathrm{kpc}$ for the cosmology adopted by the Illustris project.⁷

We classify the resolved 5343 galaxies into two types based on their Sérsic indices and stellar dynamical properties. The adopted Sérsic index here is denoted as n₀, which is derived from the 3D fitting of the PL-Sérsic profile to the stellar mass distribution (see Section 4.1 for details).

For the dynamical properties, we refer to the fraction of kinetic energy invested in the ordered rotation (Sales et al. 2012; Penoyre et al. 2017), which is given by

$$\begin{eqnarray}&&{k}_{\mathrm{rot}}=\displaystyle \frac{{K}_{\mathrm{rot}}}{K}=\displaystyle \frac{1}{K}\sum \displaystyle \frac{1}{2}{m}_{i}{\left(\displaystyle \frac{{{\boldsymbol{j}}}_{i}}{\cdot }\hat{J}| {{\boldsymbol{r}}}_{i}\times \hat{J}| \right)}^{2},\end{eqnarray} \tag{ 46 }$$

where m_i is the mass of the ith stellar particle, ${{\boldsymbol{j}}}_{i}={{\boldsymbol{r}}}_{i}\times {{\boldsymbol{v}}}_{i}$ represents the specific angular momentum, $\hat{J}$ denotes the direction of the total angular momentum, and $K=\sum \tfrac{1}{2}{m}_{i}| {{\boldsymbol{v}}}_{i}{| }^{2}$ is the total kinetic energy of the galaxy. In the calculation of k_rot, in order to avoid possible bias due to satellites at the outskirts of the galaxy, we only consider the stellar particles within the spherical radius r₉₀, which is the radius enclosing 90% of the total stellar mass.

We define the galaxies with k_rot < 0.5 and n₀ > 1.5 as "elliptical" galaxies while the rest as "disk" galaxies. We thus obtain 1362 elliptical galaxies which make up about 25% of galaxies with stellar mass larger than ${10}^{10}{h}^{-1}\,{M}_{\odot }$. This fraction of elliptical galaxies is roughly consistent with that found in observations (see Vulcani et al. 2011; De Lucia et al. 2012; Wilman et al. 2013).

The surface mass distribution for a galaxy is obtained by projecting its 3D mass distribution along the x-direction. All of the matter components are taken into account for a galaxy, including the dark matter, stars, gas, and black holes.

Among the lensing-related quantities, what we are more interested in is the convergence map. We thus generate convergence maps by scaling the surface mass distributions directly using the critical surface mass density Σ_crit. For a lens system with the lens redshift z_d = 0.178 and source redshift z_s = 0.6, we have ${{\rm{\Sigma }}}_{\mathrm{crit}}\simeq 4.0\times {10}^{15}\,{M}_{\odot }\,{\mathrm{Mpc}}^{-2}$ for the Illustris-adopted cosmology.

By looking into the 5343 convergence maps, we find that only a small fraction of galaxies can produce observable SL images. For example, there are only 358 galaxies (most of them are elliptical galaxies, as also seen in observations) with Einstein radius larger than 05. The lensing probabilities for these Illustris galaxies are likely to be underestimated because of the smoothing effect of gravitational softening on the density profiles (see Section 3.1 for the smoothing lengths for the stellar and dark matter particles). In the following analyses, we mainly focus on the density profile fittings to the galaxies regardless of their lensing probabilities.

As for the light distributions, they are directly approximated using the stellar mass distributions with a constant stellar mass-to-light ratio. For further simplification, we use the stellar mass distribution to represent the light distribution, as we are more concerned about the general shape of the light distribution rather than its total luminosity or amplitude. So, hereafter, when we refer to the "light distribution" in this paper, it is actually the mass distribution of the stellar matter. The "mass distribution" thus represents the overall mass distribution including all matter components.

Consistent with the Hubble images, the resolution of 005 is taken to pixelize the mass and light distributions where the triangular-shaped cloud algorithm is applied (Hockney & Eastwood 1981).

We calculate the AL-weighted LOSVDs for the Illustris galaxies as follows:

$$\begin{eqnarray}&&{\sigma }_{\parallel ,\mathrm{sim}}^{2}=\displaystyle \frac{\sum {\omega }_{i}{v}_{\parallel ,i}^{2}}{\sum {\omega }_{i}}-{\left(\displaystyle \frac{\sum {\omega }_{i}{v}_{\parallel ,i}}{\sum {\omega }_{i}}\right)}^{2},\end{eqnarray} \tag{ 47 }$$

where ${v}_{\parallel ,i}$ is the velocity component parallel to the line of sight for the ith stellar particle and ${\omega }_{i}={m}_{i}\exp \left(-{R}_{i}^{2}/2{\sigma }_{{\rm{fib}}}^{2}\right)$ is the weighting function, with m_i and R_i denoting the mass and projected radius of the ith stellar particle, respectively. In this paper, we use σ_fib = 1.15, which is derived according to Equation (35) in view of fiber radius 15 and a typical seeing of 169 for the SLACS lenses.

In the theoretical modeling of the AL-weighted LOSVD, a constant velocity anisotropy parameter β is assumed for each galaxy. In simulations, β is evaluated by the global anisotropy (Cappellari et al. 2007),

$$\begin{eqnarray}&&\beta =1-\displaystyle \frac{{{\rm{\Pi }}}_{\theta \theta }+{{\rm{\Pi }}}_{\phi \phi }}{2{{\rm{\Pi }}}_{{rr}}},\end{eqnarray} \tag{ 48 }$$

with ${{\rm{\Pi }}}_{{kk}}={\sum }_{i=1}^{N}{M}_{i}{\sigma }_{k,i}^{2}$, the total energy from random motions along the k-direction (i.e., the r, θ, and ϕ directions in the spherical coordinate system), where the sum runs over all the radial bins within the radius r₉₀, and M_i and σ_k,i are the total mass and velocity dispersion along the k-direction in the ith radial bin, respectively.

Based on Equation (48), we find that the global anisotropy β is not sensitive to the radial binning, because of the mass weighting in each bin, although there exists a radial variation for β in general. For instance, the estimated values of β are almost the same for 30 bins linearly spaced in r, starting from ${r}_{\min }=0\ \mathrm{kpc}$, or logarithmically spaced starting from ${r}_{\min }=0.15\ \mathrm{kpc}$. We use the β value estimated from the 30 linear bins for each galaxy in the following analyses.

For the 3D fittings to the mass and light distributions, 30 radial bins are equally spaced logarithmically in the range from $r=0.15\ \mathrm{kpc}$ to r₉₀. The radius of the ith bin is denoted as r_i, which corresponds to the mean of $\mathrm{log}(r)$ in the ith bin. We then compute the spherically averaged mass density ρ_i and light density j_i at r_i in the ith bin, and the mean mass density ${\bar{\rho }}_{i}(\lt {r}_{i})$ within r_i. In view of the larger density uncertainties in the central region due to fewer particles and the possible center offsets between the different matter components, we only use the radial bins with radius larger than $0.3\ \mathrm{kpc}$ for the 3D fittings to avoid systematics.

The total χ² to be minimized is made of three pieces, i.e.,

$$\begin{eqnarray}&&{\chi }_{3{\rm{D}}}^{2}={\chi }_{\bar{\rho }}^{2}+{\chi }_{\rho }^{2}+{\chi }_{j}^{2},\end{eqnarray} \tag{ 49 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\bar{\rho }}^{2} & = & \sum _{i=1}^{1}{\left[\mathrm{ln}{\bar{\rho }}_{i}-\mathrm{ln}{\bar{\rho }}_{i,\mathrm{BPL}}({\rho }_{c},{r}_{c},{\alpha }_{c},\alpha ,{m}_{b})\right]}^{2}\\ {\chi }_{\rho }^{2} & = & \sum _{i}{\left[\mathrm{ln}{\rho }_{i}-\mathrm{ln}{\rho }_{i,\mathrm{BPL}}({\rho }_{c},{r}_{c},{\alpha }_{c},\alpha )\right]}^{2}\\ {\chi }_{j}^{2} & = & \sum _{i}{\left[\mathrm{ln}{j}_{i}-\mathrm{ln}{j}_{i,\mathrm{PL}-{\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}}({j}_{c},{r}_{c},{\alpha }_{c},{R}_{\mathrm{eff}},n)\right]}^{2},\end{array}\end{eqnarray} \tag{ 50 }$$

with χ² denoting the mean density profile $\bar{\rho }$ for only the innermost bin used to constrain the black hole mass, ρ the mass density profile, and j the light density profile.

Note that the mass and light density profile models share the same parameters r_c and α_c, because we assume that they have compatible inner density profiles. There are in total eight free parameters when the center is fixed, e.g., at the particle position with minimum gravitational potential energy. In order to estimate the Sérsic index n₀ to better clarify the galaxy types, the light density profiles are also fitted by just minimizing the ${\chi }_{j}^{2}$.

For 2D fittings, the maximum field of view (FoV) is chosen to be 14'' × 14'' with 281 × 281 pixels. This FoV is large enough for galaxy-scale SL observations because the Einstein radii are typically less than 3'' for almost all galaxy-scale lenses that have been discovered so far.

Similar to the 3D fittings, the 2D fittings are carried out by minimizing the following χ²,

$$\begin{eqnarray}&&{\chi }_{2{\rm{D}}}^{2}={\chi }_{\bar{\kappa }}^{2}+{\chi }_{\kappa }^{2}+{\chi }_{I}^{2},\end{eqnarray} \tag{ 51 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\bar{\kappa }}^{2} & = & \sum _{i}{\left[{\bar{\kappa }}_{i}-{\bar{\kappa }}_{i,\mathrm{BPL}}(b,{r}_{c},{\alpha }_{c},\alpha ,{m}_{b})\right]}^{2}\\ {\chi }_{\kappa }^{2} & = & \sum _{i}{\left[{\kappa }_{i}-{\kappa }_{i,\mathrm{BPL}}(b,{r}_{c},{\alpha }_{c},\alpha ,q,\phi )\right]}^{2}\\ {\chi }_{I}^{2} & = & \sum _{i}{\omega }_{I}^{2}{\left[{I}_{i}-{I}_{i,\mathrm{PL}-{\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}}({j}_{c},{r}_{c},{\alpha }_{c},{R}_{\mathrm{eff}},n,{q}_{I},{\phi }_{I})\right]}^{2},\end{array}\end{eqnarray} \tag{ 52 }$$

where the data value of ${\bar{\kappa }}_{i}(\lt {R}_{\mathrm{el},i})$ at the ith pixel is the mean of the pixels of convergence within the elliptical radius ${R}_{\mathrm{el},i}=\sqrt{{{qx}}_{i}^{2}+{y}_{i}^{2}/q}$, and κ_i and I_i are the convergence and surface brightness at the ith pixel, respectively. In ${\chi }_{I}^{2}$, ${\omega }_{I}=\tfrac{M}{L}\tfrac{1}{{{\rm{\Sigma }}}_{\mathrm{crit}}}$ is applied to make the amplitude of the surface brightness comparable to the convergence map. As previously mentioned, the light intensity I here is actually the surface mass distribution of stars, i.e., with M/L = 1. As for the model parameters, b is the scale radius defined in Equation (8), q and ϕ are the axis ratio and position angle of the surface mass distribution, respectively, while q_I and ϕ_I are for the surface brightness distribution. In addition to the center (x₀, y₀), there are in total 14 free parameters in the 2D elliptical fitting to a galaxy.

Note that for the ${\chi }_{2{\rm{D}}}^{2}$ defined above, the central nine pixels of the convergence maps are masked in order to reduce the possible influence of the black hole mass on the central mass distribution due to pixelization. For the ${\chi }_{\bar{\kappa }}^{2}$ which can constrain the black hole mass, only the 16 pixels surrounding the central 9 pixels are considered in the 2D fittings.

As we know, the actual mass density profiles are steeper at larger radius and likely to be truncated at a certain radius (e.g., Navarro et al. 1997; Springel & White 1999; Drakos et al. 2017). The BPL model proposed in this paper is a model for describing the mass distribution in the relatively central region of galaxies. The BPL fittings to the surface mass distributions may be sensitive to the fitting area. To examine the effect of the fitting area, we also pay attention to an FoV of 6'' × 6'' in the 2D elliptical fittings.

In order to quantify the possible bias more accurately, we further inspect the "2D radial" fittings to the azimuthally averaged surface density profiles, which are directly calculated based on the matter particles. For the 2D radial fittings, the ${\chi }_{2{\rm{D}}}^{2}$ expression defined in Equation (51) is adopted but with

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\bar{\kappa }}^{2} & = & \sum _{i=1}^{1}{\left[\mathrm{ln}{\bar{\kappa }}_{i}-\mathrm{ln}{\bar{\kappa }}_{i,\mathrm{BPL}}(b,{r}_{c},{\alpha }_{c},\alpha ,{m}_{b})\right]}^{2}\\ {\chi }_{\kappa }^{2} & = & \sum _{i}{\left[\mathrm{ln}{\kappa }_{i}-\mathrm{ln}{\kappa }_{i,\mathrm{BPL}}(b,{r}_{c},{\alpha }_{c},\alpha )\right]}^{2}\\ {\chi }_{I}^{2} & = & \sum _{i}{\left[\mathrm{ln}{I}_{i}-\mathrm{ln}{I}_{i,\mathrm{PL}-{\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}}({j}_{c},{r}_{c},{\alpha }_{c},{R}_{\mathrm{eff}},n)\right]}^{2},\end{array}\end{eqnarray} \tag{ 53 }$$

where i indicates the ith radial bin and the center is fixed at the position corresponding to the 3D center. A total of 30 bins are equally spaced in logarithmic scale in the range of $0.15\,\mathrm{kpc}$ to the projected 90% light radius R₉₀. Only the radial bins with a projected radius R larger than $0.3\,\mathrm{kpc}$ are used in the analysis. Similar to the 3D fittings, only the innermost bin is adopted to calculate ${\chi }_{\bar{\kappa }}^{2}$. For reference, Table 1 briefly summarizes all the fitting methods investigated above.

Table 1.  3D and 2D Density Profile Fittings

Method	χ2 Performance	Fitting Range or Area
3D (Radial)	${\chi }_{\bar{\rho }}^{2}+{\chi }_{\rho }^{2}+{\chi }_{j}^{2}$ using logarithmic radial bins	$0.3\,\mathrm{kpc}\,\to \,{r}_{90}$
2D (Elliptical)	${\chi }_{\bar{\kappa }}^{2}+{\chi }_{\kappa }^{2}+{\chi }_{I}^{2}$ on pixelated maps	14'' × 14'' or 6'' × 6'' (central 9 pixels are masked)
2D Radial	${\chi }_{\bar{\kappa }}^{2}+{\chi }_{\kappa }^{2}+{\chi }_{I}^{2}$ using logarithmic radial bins	$0.3\,\mathrm{kpc}\,\to \,{R}_{90}$

Note. For the 3D and 2D radial fittings, only the innermost bin is considered for the calculation of ${\chi }_{\bar{\rho }}^{2}$ and ${\chi }_{\bar{\kappa }}^{2}$. For the 2D elliptical fittings, only the 16 pixels around the central 9 pixels are used for the ${\chi }_{\bar{\kappa }}^{2}$.

Download table as:  ASCIITypeset image

Figure 3 displays the 3D fittings for a couple of typical galaxies. We notice that, for most of the galaxies in Figure 3, the BPL and PL-Sérsic models work equally well to describe the relevant density profiles within r₉₀. However, obvious deviations exist for some galaxies in the extremely central region or around the break radius. The large fluctuations in the central region are likely due to the limited resolution of the Illustris simulation and the possible center offset between different matter components. The deviations around the break radius are expected because the BPL model is a piecewise function which is continuous but not smooth at the break radius.

Zoom In Zoom Out Reset image size

In Figure 4, we show the statistical results of the 3D fittings. As shown in the left panel of Figure 4, the biases of the BPL fittings to $\bar{\rho }(\lt r)$ (corresponding to the total mass distribution) are typically less than 5% for elliptical galaxies and 10% for disk galaxies. If we look into the corresponding fittings to the volume density profiles ρ(<r) shown in the second panel, an obviously increasing trend is found for the biases at larger radius. The reason is that the true density profile is steeper at larger radius whereas the slope of the BPL model here is mainly determined by the mass distribution in the relatively inner region. Thus, the density profile at a larger radius tends to be overestimated by the BPL model. The biases can reach 30% and 10% at r₉₀ from a negative bias of about −10% for the disk and elliptical galaxies, respectively.

Zoom In Zoom Out Reset image size

The third panel shows the relative deviations of the PL-Sérsic profile fittings from the light distributions. It is shown that the deviation fluctuations are somewhat larger than the BPL fittings to the mass distributions, but still reasonable in consideration of the interplay between the mass and light density profiles in the ${\chi }_{3{\rm{D}}}^{2}$ fittings.

The slope distributions in the fourth panel illustrate the obvious differences between the inner and outer density profiles within r₉₀. The inner slope is about 0.5 for both the elliptical and disk galaxies. The outer slope is about 2 with a scatter of ∼0.16 for the elliptical galaxies. However, for disk galaxies, the fitted outer slope is much flatter than for elliptical galaxies and has a larger dispersion. In the last panel, we find that the distributions of the break radii, if normalized, are nearly the same for elliptical and disk galaxies. The modes of the r_c distributions are about $2\,\mathrm{kpc}$, which is larger than the softening length of dark matter particles.

To sum up, the 3D fittings inspected in this subsection demonstrate that the BPL model is feasible to describe the mass density profiles of Illustris galaxies within a certain radius and the PL-Sérsic model is also good enough to measure the light density profiles.

Figure 5 shows two examples of the 2D elliptical fittings to the convergence maps with an FoV of 14'' × 14''. As shown in Figure 5, the BPL model can fit very well the 2D mass distributions of the two galaxies inspected here. For both galaxies, the residuals are sufficiently small, and from the rightmost panels, we can see that the BPL model can provide excellent fits to the elliptical radial density profiles.

Zoom In Zoom Out Reset image size

By comparing the parameter values derived from 2D fittings (black numbers) and 3D fittings (magenta numbers), we find that they are not always consistent, especially for the disk galaxies. This can be attributed to the projection effect caused by the limitation of the BPL model and the nonspherical shape of galaxies. We present more detailed comparisons between the 2D and 3D BPL fittings in Section 5.3.

We now move to the statistical results of the 2D fittings, which are displayed in Figure 6. It is shown that, for all three methods of the 2D fitting, the BPL model can estimate the mean convergence $\bar{\kappa }$ maps of elliptical galaxies very well with a negligible bias within the radius R₉₀. However, for the disk galaxies, the mean convergence maps tend to be overestimated in the inner region, especially for the fittings with larger FoV.

Zoom In Zoom Out Reset image size

For the BPL fittings to the convergence maps, as expected, biases always exist at much larger radius both for elliptical and disk galaxies. A smaller FoV may cause a worse fitting of convergence in the outer region. However, as indicated by the left panels, the averaged convergence within R₉₀ is not very susceptible to the overestimation of the convergence at a relatively larger radius.

As for the PL-Sérsic fittings shown in the third column, the scatters are relatively larger for the 2D elliptical fittings. One reason is that the adopted FoV is not large enough for some massive galaxies. Therefore, the large deviations may emerge for massive galaxies when the radius is scaled by R₉₀.

In the fourth column of Figure 6, it is shown that the inner and outer density profiles can be still clearly separated by the slopes estimated from 2D fittings. The slope distributions are more consistent between the elliptical and disk galaxies. However, one may realize that both the inner and outer slopes are systematically higher than those from 3D fittings. For example, the modes of the inner slope distributions increase to about 0.8 from 0.5 for both the elliptical and disk galaxies. The mode of the outer slope distribution is biased to be about 2.3 (for 2D 14'' × 14'' fittings ) or 2.1 (for 2D 6'' × 6'' and radial fittings) for disk galaxies. However, for the elliptical galaxies, the outer slope distribution is only slightly biased.

In the last column, we display the break radius distributions, which are much wider than the distributions from 3D fittings, especially for the disk galaxies. The break radii for elliptical galaxies are systematically smaller than those for disk galaxies

Based on the results presented in this subsection, we realize that the performance of the 2D BPL fittings is sensitive to the binning or weighting methods. Even so, the inner and outer density profiles can be clearly differentiated whichever fitting method is used. More importantly, we find that the mean convergence maps can be recovered very well by the BPL fittings, especially for the elliptical galaxies.

As a model used to describe the relatively central region of galaxies, both the 2D and 3D BPL model fittings tend to overestimate the true mass density profiles outside a certain radius, because the outer slope of the BPL model is fixed while the true mass density profile usually decreases more and more rapidly with increasing radius. In addition to the nonspherical shape of galaxies, inconsistency must exist between the 2D and 3D BPL model fittings.

Figure 7 presents the one-to-one comparisons of parameter values of ρ_c, r_c, α_c, and α between the 2D and 3D fittings, where the ρ_c for a 2D fitting is directly derived from b, α, and r_c according to Equation (8). In these comparisons, the parameter values of 2D elliptical fittings are adopted directly to be compared by ignoring the nonspherical shape of galaxies.

Zoom In Zoom Out Reset image size

From Figure 7, we can realize that the 2D and 3D parameter values are not always strongly correlated, where large scatters and biases may exist, especially for the disk galaxies. We find, compared to the 3D fittings, that the values of ρ_c estimated by 2D fittings are slightly underestimated in general but with a larger break radius r_c, and steeper inner and outer slopes.

Figure 8 presents the comparisons between the 3D density profiles (i.e., ${\bar{\rho }}_{2{\rm{D}}}$, ρ_2D, and j_2D) predicted from 2D fittings and the true 3D density profiles. We notice that the scatter is large for the deprojected profiles of the 2D fittings. We find that the inner mass density profiles for disk galaxies are significantly upturned due to the steeper inner slope of 2D fittings. However, for elliptical galaxies, the fitting bias is not significant in the inner region.

Zoom In Zoom Out Reset image size

One may realize that a small negative bias exists around the 90% light radius for the ρ_2D predicted from 2D elliptical fittings with FoV 14'' × 14'' and 2D radial fittings. This is because the slopes around r₉₀ are relatively steeper for the 2D BPL fittings than for the corresponding 3D density profiles. However, this does not indicate that the mass density profiles are still underestimated at radius much larger than r₉₀. We examine the fitting bias at radius larger than two or three times r₉₀, e.g., and find that the true 3D mass density profiles are overestimated in general by the 2D fittings beyond a certain large radius.

In addition to the comparisons shown in Figure 8, we also compare the projections of 3D fittings (i.e., ${\bar{\kappa }}_{3{\rm{D}}}$, κ_3D, and I_3D) with the true surface density profiles in Figure 9. It is noticeable that, for disk galaxies, the surface mass density profiles are significantly biased high from inner to outer regions by the projections of 3D BPL fittings, and the positive biases are more serious at larger radius. While, for elliptical galaxies, the true surface mass density profiles can be described very well by ${\bar{\kappa }}_{3{\rm{D}}}$ and κ_3D within the projected half-light radius, although overestimation still appears at larger radius.

Zoom In Zoom Out Reset image size

The increasing trend of biases shown in Figure 9 can be understood easily by considering the flatter outer slopes of 3D BPL model fittings. Because the true density profiles decrease faster than those estimated by BPL model fittings at larger radius, the relative biases are expected to be more significant at larger radius. Compared to disk galaxies, the biases for elliptical galaxies are systematically smaller because their outer slopes in 3D fittings are about 2 and much steeper than the slopes (∼1.5) for disk galaxies. Note that, as shown in Figure 4, the BPL model does not always overestimate the volume density profiles and may underestimate them within the fitting range. Therefore, by projection, the bias in the relatively inner region can be reduced and may be not significant, e.g., for elliptical galaxies.

In contrast to BPL model fittings, the last panel of Figure 9 demonstrates the good performance of PL-Sérsic model fittings to light density profiles.

In short, this subsection further illustrates the effects of projection, fitting ranges, and methods on density profile fittings. We compare in detail the true volume density profiles with the deprojections of the 2D fittings and the true surface density profiles with the projections of the 3D fittings. We find that the BPL model can give a reasonable prediction of the 3D density profiles within r₉₀ from the deprojections of the 2D fittings (see Figure 8). However, the BPL model has limitations. The 3D BPL fittings cannot be directly projected to estimate the surface mass density profiles, especially for disk galaxies (see Figure 9). As a lens mass model, the BPL model is mainly proposed to describe the surface mass distributions and the relevant lensing quantities. So, we care more about the accuracy of the 2D fittings and their consistency with the true volume density profiles within a certain radius. The large biases presented by the projections of the 3D BPL fittings are just for theoretical investigation and not likely to happen in actual lensing analyses.

Velocity dispersions can help us improve the constraints on the mass distributions of galaxies. In this subsection, we inspect whether the AL-weighted LOSVDs ${\sigma }_{\parallel ,\mathrm{pred}}$ predicted from the BPL mass model fittings are consistent with the directly "observed" ones ${\sigma }_{\parallel ,\mathrm{sim}}$.

In Figure 10, we show the comparisons between the "predicted" ${\sigma }_{\parallel ,\mathrm{pred}}$ and the true "observed" ${\sigma }_{\parallel ,\mathrm{sim}}$ for the 3D and 2D fittings. It is found that the predicted and true AL-weighted LOSVDs are strongly correlated. For the 3D fittings, the bias is negligible, and the scatter is about 7%. However, for the 2D fittings, small positive biases exist as indicated by the positive shift of the scatter diagrams and the statistical distributions of ${\sigma }_{\parallel ,\mathrm{pred}}/{\sigma }_{\parallel ,\mathrm{sim}}$. The biases for the 2D fittings are larger than 3% but typically less than 6%.

Zoom In Zoom Out Reset image size

We find that the velocity dispersion bias can be corrected by accounting for the imperfect fittings of the BPL model and the projection effect. The BPL model can fit well the 2D mass distributions. However, the deprojection of the 2D fittings has significant scatter compared to the true 3D density profiles. The projection effect can complicate not only the reconstruction of mass distributions but also the prediction of AL-weighted LOSVDs.

For the 3D fittings, we conclude that there is no velocity dispersion bias (i.e., b_σ = 1) for the BPL model. For the 2D fittings, we find that the velocity dispersion bias can be roughly estimated by

$$\begin{eqnarray}&&{b}_{\sigma }\simeq 1.015{q}_{* }^{-0.07},\end{eqnarray} \tag{ 54 }$$

where q_* (<1) is the axial ratio of the stellar mass distribution. The quantity q_* can be estimated by the 2D elliptical PL-Sérsic profile fitting or the inertia tensor of stellar particles within a certain area, e.g., enclosed by the elliptical 90% light radius ${R}_{\mathrm{el},90}$. We find that the b_σ–q_* relation is not very sensitive to the measuring methods of q_* that we have investigated. The index −0.07 demonstrates the weak dependence of b_σ on the observed ellipticity of galaxies. More spherical galaxies are less subject to the projection effect but still suffer from the limitation of mass models.

By scaling the predicted ${\sigma }_{\parallel ,\mathrm{pred}}$ by b_σ, we can then get the corrected AL-weighted LOSVD ${\tilde{\sigma }}_{\parallel ,\mathrm{pred}}={\sigma }_{\parallel ,\mathrm{pred}}/{b}_{\sigma }$. In Figure 11, we present the comparisons between ${\tilde{\sigma }}_{\parallel ,\mathrm{pred}}$ and ${\sigma }_{\parallel ,\mathrm{sim}}$. Evidently, the bias of ${\sigma }_{\parallel ,\mathrm{pred}}$ from ${\sigma }_{\parallel ,\mathrm{sim}}$ is well corrected for by the bias factor b_σ for BPL model fittings. The velocity bias is finally reduced to be no more than 2% for the 2D fittings, and the intrinsic scatter is typically about 6%.

Zoom In Zoom Out Reset image size

It should be mentioned that the velocity anisotropy parameter β used above for each galaxy is assumed to be a known quantity in the calculation of the AL-weighted LOSVD. However, in observations, β cannot be easily determined even if the velocity dispersion profile is observed because of the degeneracy between the mass distribution and velocity anisotropy (Gerhard 1993; Kronawitter et al. 2000; Magorrian & Ballantyne 2001; Mamon & Boué 2010). By looking at the excellent consistency between ${\tilde{\sigma }}_{\parallel ,\mathrm{pred}}$ and ${\sigma }_{\parallel ,\mathrm{sim}}$ illustrated in Figure 11, we propose to measure an effective velocity anisotropy β_eff by solving the equation ${\sigma }_{\parallel ,\mathrm{pred}}(\beta )={b}_{\sigma }{\sigma }_{\parallel ,\mathrm{sim}}$.

In Figure 12, we show the comparisons between the predicted effective velocity anisotropy β_eff and the directly measured β for all of the fitting methods. According to the scatter diagrams, we are aware that the uncertainties for the inferred β_eff are significant, demonstrating the difficulty in inferring the velocity anisotropy for an individual galaxy by resorting to the SL mass measurement and AL-weighted LOSVD observation. However, as indicated by the black lines, the median of the β_eff values for a sample of galaxies with nearly the same β is basically unbiased for most of the fitting methods.

Zoom In Zoom Out Reset image size

There exist negative biases (especially for the disk galaxies) for the 2D elliptical fittings with FoV 14'' × 14''. These biases may be corrected for by introducing a more general velocity bias factor b_σ, which also depends on the fitting area. However, we know that such a large FoV of 14'' × 14'' may not be necessary for galaxy-scale lensed image reconstructions in real observations.

We thus argue that it is possible to measure the distribution of velocity anisotropy parameters for certain types of galaxies. For instance, the second and bottom rows of Figure 12 present the comparisons between the distributions of β_eff (thick histograms) and β (thin histograms) for elliptical (red) and disk (blue) galaxies, respectively. We can see that the distributions of β_eff and β are very well consistent with each other for most of the fittings. For the 2D fittings to the disk galaxies with FoV 14'' × 14'', an obvious negative bias exists for the predicted β_eff distributions because the velocity dispersion bias is not well corrected.