A Robust Mass Estimator for Dark Matter Subhalo Perturbations in Strong Gravitational Lenses

The perturbation in a lensed image produced by a subhalo is only evident if that subhalo lies near the tangential critical curve of the lens, distorting the critical curve. Near the warped critical curve, the surface brightness is altered, and additional images may even appear. The size of this critical curve distortion is comparable to the scale at which the lensed surface brightness is significantly perturbed by the subhalo.

We define the subhalo perturbation radius to be the distance from the subhalo to the critical curve along the direction where the magnification is perturbed the most by the subhalo. This radius, which we denote as $\delta {\theta }_{c}$, can be loosely thought of as the "radius of maximum warping" of the critical curve, i.e., the farthest distance by which the subhalo "pushes" the critical curve outward compared to the position of the critical curve without the subhalo present. As we show in Appendix A, $\delta {\theta }_{c}$ can be found by solving the following equation for r:

$$\begin{eqnarray}&&1-{\kappa }_{0}({{\boldsymbol{\theta }}}_{s}+{\boldsymbol{r}})-{{\rm{\Gamma }}}_{\mathrm{tot}}({{\boldsymbol{\theta }}}_{s}+{\boldsymbol{r}})={\bar{\kappa }}_{s}(r).\end{eqnarray} \tag{ 1 }$$

Here, ${\kappa }_{0}({\boldsymbol{\theta }})$ refers to the projected density profile of the host galaxy, ${{\rm{\Gamma }}}_{\mathrm{tot}}$ is the magnitude of the total unperturbed shear, i.e., the complex sum of the shear from the host galaxy ${{\rm{\Gamma }}}_{0}({\boldsymbol{\theta }})$ plus the external shear ${{\rm{\Gamma }}}_{\mathrm{ext}}$ created by neighboring galaxies, and ${\bar{\kappa }}_{s}(r)$ is the average projected density of the subhalo a distance r from the center of the subhalo. The direction of ${\boldsymbol{r}}$ corresponds to the direction along which the subhalo shear aligns with that of the total unperturbed shear ${{\boldsymbol{\Gamma }}}_{\mathrm{tot}}$. As long as the external shear is not too large (${{\rm{\Gamma }}}_{\mathrm{ext}}\lesssim 0.2$), then to good approximation, ${\boldsymbol{r}}$ simply lies along the radial direction with respect to the host galaxy (this is discussed in detail in Appendix B.1).

Since the projected subhalo mass enclosed within a given radius is simply ${m}_{\mathrm{sub}}(r)=\bar{\kappa }(r)\pi {r}^{2}{{\rm{\Sigma }}}_{\mathrm{cr}}$, it is evident that the solution $\delta {\theta }_{c}$ to Equation (1) depends on the subhalo mass enclosed within $r=\delta {\theta }_{c}$, regardless of how the subhalo mass is distributed within that radius. Thus, it follows that if a subhalo is detected in a gravitational lens and if the best-fit model reproduces both the host halo parameters and the subhalo perturbation radius $\delta {\theta }_{c}$ well, then the inferred subhalo mass within $\delta {\theta }_{c}$ will be accurate, regardless of the actual density profile of the subhalo.

In reality, unless the subhalo density profile perfectly matches the model, the host halo parameters will also be perturbed in order to better fit the lensed images in the vicinity of the subhalo. To get an idea of how this affects the inferred subhalo mass, let us first assume an idealized scenario where the host galaxy is axisymmetric and no external shear is present. In this case, Equation (1) can be written as

$$\begin{eqnarray}&&1-{\bar{\kappa }}_{0}({\theta }_{s}+r)={\bar{\kappa }}_{s}(r),\end{eqnarray} \tag{ 2 }$$

where we have used the fact that for an axisymmetric lens, ${\kappa }_{0}+{{\rm{\Gamma }}}_{0}={\bar{\kappa }}_{0}$. Note the direction of maximum warping is simply the radial direction in this case. If we assume that the projected density of the host galaxy follows a power-law profile, then we can write

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

where b is the Einstein radius of the host galaxy. We then have

$$\begin{eqnarray}&&{\bar{\kappa }}_{0}(\theta )={\left(\displaystyle \frac{b}{\theta }\right)}^{\alpha }.\end{eqnarray} \tag{ 4 }$$

In the vicinity of the subhalo, the host galaxy deflection profile ${\bar{\kappa }}_{0}\theta$ can be Taylor expanded to first order, which gives (after dividing by $\theta ={\theta }_{s}+r$)

$$\begin{eqnarray}&&{\bar{\kappa }}_{0}({\theta }_{s}+r)\approx {\left(\displaystyle \frac{b}{{\theta }_{s}}\right)}^{\alpha }\left[1-\alpha \displaystyle \frac{r}{{\theta }_{s}+r}\right].\end{eqnarray} \tag{ 5 }$$

Note that in the isothermal case ($\alpha =1$), this formula becomes exact; this is because the deflection is constant for an isothermal profile, which means that higher-order terms vanish anyway in that case (indeed, this is our motivation for expanding the deflection instead of expanding ${\bar{\kappa }}_{0}$ directly).

Let us assume that the subhalo is close to the critical curve, but not exactly on it; more precisely, if the distance between the subhalo and the unperturbed critical curve is ${\rm{\Delta }}\theta \equiv b-{\theta }_{s}$, then we assume that ${\rm{\Delta }}\theta /b\ll 1$. (If the subhalo is not close enough to the critical curve to satisfy this condition, its effects are unlikely to be observable in any case!) If we substitute $b={\theta }_{s}+{\rm{\Delta }}\theta$ into Equation (5), expand to first order in ${\rm{\Delta }}\theta /b$, and substitute the result into Equation (2), we find that

$$\begin{eqnarray}&&\alpha \left(-\displaystyle \frac{{\rm{\Delta }}\theta }{b}+\displaystyle \frac{r}{{\theta }_{s}+r}\right)\approx {\bar{\kappa }}_{s}(r),\end{eqnarray} \tag{ 6 }$$

and hence

$$\begin{eqnarray}&&1-{\bar{\kappa }}_{0,\mathrm{iso}}({\theta }_{s}+r)\approx \displaystyle \frac{{\bar{\kappa }}_{s}(r)}{\alpha },\end{eqnarray} \tag{ 7 }$$

where ${\kappa }_{0,\mathrm{iso}}$ refers to the isothermal ($\alpha =1$) profile. This equation tells us that if the host galaxy deviates from isothermal, then we will get the same subhalo perturbation radius if we have an isothermal host and the same subhalo mass scaled by $1/\alpha$. For the more general and realistic case where the host galaxy has some ellipticity and external shear is present, we show in Appendix B.2 that the inferred subhalo mass ${m}_{\mathrm{sub}}(\delta {\theta }_{c})$ still scales as $1/\alpha$, with only a small amount of error due to the external shear (see Equation (39) and following discussion). Thus, when fitting to lens data, as long as our model reproduces the position of the subhalo well, we can say that

$$\begin{eqnarray}&&\displaystyle \frac{{m}_{\mathrm{sub},\mathrm{fit}}(\delta {\theta }_{c})}{{\alpha }_{\mathrm{fit}}}\approx \displaystyle \frac{{m}_{\mathrm{sub},\mathrm{true}}(\delta {\theta }_{c})}{{\alpha }_{\mathrm{true}}}.\end{eqnarray} \tag{ 8 }$$

We therefore define the effective subhalo lensing mass ${\tilde{m}}_{\mathrm{sub}}\equiv {m}_{\mathrm{sub}}(\delta {\theta }_{c})/\alpha$. We demonstrate by fitting to simulated data in Section 5 that unlike the total subhalo mass, ${\tilde{m}}_{\mathrm{sub}}$ can indeed be inferred robustly even if the subhalo density profile is modeled inaccurately.

It is useful to get an idea of the scale of $\delta {\theta }_{c}$. If both the host galaxy and the subhalo are modeled as a singular isothermal spheres, then Equation (2) becomes

$$\begin{eqnarray}&&1-\displaystyle \frac{b}{{\theta }_{s}+r}=\displaystyle \frac{{b}_{s}}{r},\end{eqnarray} \tag{ 9 }$$

where b is the Einstein radius of the host galaxy and b_s is the Einstein radius of the subhalo. This leads to a quadratic equation in r if we assume that the subhalo is right on the critical curve (i.e., ${\theta }_{s}=b$) and that ${b}_{s}\ll b$, the solution is approximately

$$\begin{eqnarray}&&\delta {\theta }_{c}\approx \sqrt{{{bb}}_{s}}+\displaystyle \frac{1}{2}{b}_{s},\end{eqnarray} \tag{ 10 }$$

where b is the Einstein radius of the host galaxy and b_s is the Einstein radius of the subhalo. Although the subhalo projected mass may extend well beyond this radius, it is typically only in the neighborhood of $\delta {\theta }_{c}$ that the subhalo perturbation will be apparent and distinguishable from a smooth lens model.

Equation (10) is only accurate in the very idealized case we have described. In general, the subhalo will not sit exactly on the unperturbed critical curve, and the host galaxy will have some ellipticity as well as external shear from neighboring galaxies, all of which will tend to reduce the magnification and hence make $\delta {\theta }_{c}$ somewhat smaller. A general analytic formula is derived in Appendix B that takes all of these effects into account (see Equation (47)). For those who are not in the lens modeling business, however, it is useful to have a simpler approximate formula, which is accurate to within 30% in realistic lens scenarios. We recommend the following formula (which is derived in Appendix B):

$$\begin{eqnarray}&&\delta {\theta }_{c}\approx \sqrt{\displaystyle \frac{\epsilon {\theta }_{s}{b}_{s}}{\alpha }}\left(1-\displaystyle \frac{x}{2}+\displaystyle \frac{{x}^{2}}{8}\right)+\displaystyle \frac{\epsilon {b}_{s}}{2\alpha }(1-x),\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&x=\xi \sqrt{\displaystyle \frac{{\theta }_{s}}{b\epsilon \alpha }}.\end{eqnarray} \tag{ 12 }$$

In this formula, ${\theta }_{s}$ is the projected angular separation between the subhalo and the center of the host galaxy, and α is the log-slope of the host galaxy density profile. The constants and ξ are determined by the tidal radius of the subhalo. If the subhalo distance to the host galaxy center r₀ (including the component along the line of sight) is much larger than the Einstein radius b, then we have $\epsilon \approx 1$, $\xi \approx 0$. On the other hand, if we make the assumption that ${r}_{0}\approx b$, as is often assumed in lens modeling studies, then we have $\epsilon =0.879$, $\xi =0.293$. Otherwise, for a given assumed value for r₀ (or equivalently r_t), and ξ can be calculated using Equations (43)–(45).

We find that Equation (11) is accurate to within 1% if the subhalo sits exactly on the (unperturbed) critical curve, provided the external shear ${{\rm{\Gamma }}}_{\mathrm{ext}}\lesssim 0.2$. The farther away the subhalo is from the critical curve, the worse this approximation becomes—generally it will give a value for $\delta {\theta }_{c}$ that is too high. In the simulated lenses we examine in Section 5, Equation (11) is greater than the true $\delta {\theta }_{c}$ by about 17%–25% depending on the simulation; by comparison, the more precise Equation (47) is accurate to within 2% (see Table 4 in Appendix B for a comparison between the different estimates of $\delta {\theta }_{c}$). However, in the simulated lenses we examine in Section 5, we find that the optimal radius where the effective subhalo lensing mass can be robustly estimated is slightly larger than the actual value for $\delta {\theta }_{c}$ in each case, and this may tend to be generally the case. One can see from Figure 7 that Equation (11) works better than using the actual $\delta {\theta }_{c}$. Therefore, we generally recommend using this approximation and calculating the effective lensing mass within this radius. Lens modelers who desire a more precise value can solve Equation (1) numerically or use the analytic formula of Equation (47).

For the subhalo, we require a projected density profile with a variable log-slope and a smooth truncation at a specified tidal radius. To accomplish this, we choose the "cuspy halo model" of Muñoz et al. (2001), whose density profile is

$$\begin{eqnarray}&&\rho =\displaystyle \frac{{\rho }_{0}{r}_{t}^{n}}{{r}^{\gamma }{({r}^{2}+{r}_{t}^{2})}^{(n-\gamma )/2}}.\end{eqnarray} \tag{ 15 }$$

This model has an inner slope γ and an outer slope n. The pseudo-Jaffe profile corresponds to the special case $\gamma =2$, n = 4, so that the inner slope corresponds to the profile of an isothermal sphere, while the outer slope achieves a smooth truncation around the tidal radius r_t. For general $\gamma ,n$ values, the projected (spherically symmetric) density profile $\kappa (r)$ and deflection $\alpha (r)$ have analytic solutions in terms of the Gaussian hypergeometric function as given in Muñoz et al. (2001), which we use here. As with the pseudo-Jaffe profile, we fix n = 4 for the outer slope, but allow the inner slope to vary.

We also consider the effect of varying the subhalo distance from the host galaxy center (while keeping the projected position of the subhalo fixed). The resulting tidal radius of the subhalo can be estimated by considering the subhalo before truncation, whose density profile can be written as $\rho ={\rho }_{0}{({r}_{t}/r)}^{\gamma }$. Using the usual formula for the Jacobi radius, we find

$$\begin{eqnarray}&&{r}_{t}=\displaystyle \frac{{r}_{0}^{2}{\kappa }_{s,0}}{(3-\gamma )b},\end{eqnarray} \tag{ 16 }$$

where r₀ is the distance from the subhalo to the center of the host galaxy, b is the Einstein radius of the host, and ${\kappa }_{s,0}=2\pi {\rho }_{0}{r}_{t}/{{\rm{\Sigma }}}_{\mathrm{cr}}$ (our definition differs from that of Muñoz et al. 2001 by a factor of $2\pi$).

For model fitting, we follow in the footsteps of previous work and model the subhalo explicitly with a pseudo-Jaffe profile, which corresponds to $\gamma =2$, n = 4 and whose projected density profile has a simple analytic form,

$$\begin{eqnarray}&&{\kappa }_{s}(r)=\displaystyle \frac{{b}_{s}}{2}[{r}^{-1}-{({r}^{2}+{r}_{t}^{2})}^{-\tfrac{1}{2}}].\end{eqnarray} \tag{ 17 }$$

The parameter b_s is related to ${\kappa }_{s,0}$ by ${\kappa }_{s,0}={b}_{s}/{r}_{t};$ plugging this into Equation (16) gives ${r}_{t}={r}_{0}\sqrt{{b}_{s}/b}$. Rather than vary r_t explicitly, we set ${r}_{0}=b$, which then fixes ${r}_{t}=\sqrt{{{bb}}_{s}}$, recovering the usual formula for the tidal radius under this assumption. Thus, our model subhalo parameters to be varied are b_s along with the center coordinates ${x}_{c,s}$ and ${y}_{c,s}$.

To investigate the effects of modeling a subhalo with an incorrect density profile, we consider four different gravitational lens simulations, which are listed and described in Table 1. In each scenario, the "actual" subhalo density profile is given by Equation (15) and has the same projected position in the lens plane. For each simulated lens, we follow the footsteps of previous work and fit the subhalo with a pseudo-Jaffe model (Equation (17)), whose tidal radius is fixed by the condition that its distance to the host galaxy center r₀ is equal to the Einstein radius $b\approx 1\buildrel{\prime\prime}\over{.} 6\approx 7.4$ kpc. Compared to the model subhalo, the "actual" subhalo density profile is identical in Simulation 1, has a shallower inner slope γ in simulations 2–3, and has a much larger distance from the host galaxy r₀ = 100 kpc in Simulation 4. Note that a larger r₀ also implies a significantly larger tidal radius r_t according to Equation (16).

The simulated image and source planes for Simulation 1 are shown in Figure 2. In the first simulation (which is our "control" scenario since the actual subhalo profile matches the model), we choose the normalization of the density profile ${\kappa }_{s,0}$ so that the total subhalo mass is approximately ${10}^{9}\,{M}_{\odot }$. In the three remaining simulations, the normalization is chosen so that the mass within the subhalo perturbation radius $\delta {\theta }_{c}$ is roughly the same in each case, thus ensuring that the scale of the perturbation is not radically different in each case. Figure 3 plots the surface brightness and critical curves in each simulation. We see that the scale of the critical curve perturbation is indeed approximately the same, although the shape of the perturbation differs: the steeper the subhalo density profile, the more the critical curve is "pinched"; likewise, the critical curve is pinched slightly more for a larger tidal radius.

Zoom In Zoom Out Reset image size

For each scenario, we model the smooth lens component as a power-law ellipsoid (Barkana 1998; Tessore & Metcalf 2015), whose density profile can be expressed as

$$\begin{eqnarray}&&{\kappa }_{0}(x,y)=\displaystyle \frac{2-\alpha }{2}{b}^{\alpha }{({{qx}}^{2}+{y}^{2}/q)}^{-\tfrac{\alpha }{2}},\end{eqnarray} \tag{ 18 }$$

where q is the axis ratio. In this form, b can be thought of as the average Einstein radius, since the semimajor axis of the critical curve is equal to $b/\sqrt{q}$, while the semiminor axis is equal to $b\sqrt{q}$. Equation (18) is then rotated by an angle θ, where we define the θ to be the (counterclockwise) angle between the galaxy semimajor axis and the y-axis of the observer's coordinate system. Finally, in addition to the host galaxy, we add an external shear term that we parameterize by the shear components ${{\rm{\Gamma }}}_{1}$, ${{\rm{\Gamma }}}_{2}$.

To model our simulated lenses, we use the method of pixel-based source reconstruction, where the source is pixellated and the pixel values are inferred by linear inversion. Since this method is well described in many papers (see e.g., Warren & Dye 2003; Suyu et al. 2006; Tagore & Keeton 2014), we do not describe it in great detail in here, but rather point out distinct features of our algorithm.

For a given set of lens parameters, we use linear three-point interpolation for the ray-tracing to construct the lensing matrix L, which is well described in Tagore & Keeton (2014). We use curvature regularization for smoothing of the source, and perform a matrix inversion to obtain the inferred source surface brightness pixel map. The reconstructed source is pixellated using an adaptive Cartesian grid (Dye & Warren 2005), where the first-level pixels split into four smaller pixels if the magnification of a source pixel is higher than $4{\mu }_{0}$, where ${\mu }_{0}$ is called the splitting threshold. This process is then repeated recursively for the higher-level pixels. For simplicity, the number of first-level source pixels is fixed to 0.3 times the number of image pixels being fit. (We note, however, that more sophisticated adaptive source grids are possible; see e.g., Vegetti & Koopmans 2009 and Nightingale & Dye 2015.)

In addition to these points, we enumerate here the features that are unique in our implementation:

• (1)  
  In principle, the splitting threshold can be varied and the optimal value inferred from maximizing the Bayesian evidence. In practice, we find that ${\mu }_{0}\approx 5\mbox{--}6$ is typically preferred, so to save computation time, we fix the splitting threshold to ${\mu }_{0}=5$.
• (2)  
  Although our adaptive source grid is Cartesian, its borders are not fixed, but instead are redrawn each time the lens parameters are varied so that it is just large enough to include the inferred source signal. In this way, we exclude as many pixels as possible that only contain noise, which could degrade the fit. In practice, we accomplish this by drawing the box so that it just barely contains all pixels for whom the surface brightness is greater than 1.3 times the pixel noise. To be conservative, we then extend the length and width of the box by 15% to ensure that all the signal is being included. A side benefit of this is that computation time is saved, since typically fewer pixels are being evaluated compared to a fixed grid.
• (3)  
  As is common practice, we mask regions of the image that only contain noisy pixels, so that the majority of the pixels being fit to contain the actual signal. However, there is a risk that certain solutions produce images in the masked region (i.e., outside the region being fit) that should not exist, and this may be missed. To avoid this, after inverting to find the source surface brightness distribution, we then reconstruct the lensing matrix using all of the pixels in the image, including the masked regions, and find the resulting surface brightness in each pixel. We then impose a penalty function prior, so that if the surface brightness in the masked region rises above 30% of the brightest pixel in the image, a steep penalty is imposed and this solution is discouraged.
• (4)  
  In subhalo modeling studies that have used pixellated source reconstruction to date, an nonparametric approach (Vegetti & Koopmans 2009; Hezaveh et al. 2016) is typically first used to attempt to locate where the lensing potential is being perturbed by a subhalo, before proceeding to the full parameterized subhalo model. One very important reason for doing so, rather than a parametric fit alone, is that if the mass distribution of the smooth lens component is more complicated than the chosen parametric model, spurious detections can result in an attempt to improve the fit (Vegetti et al. 2014). In this paper, since we are working with simulated data, we skip this step entirely; from a practical standpoint, this is possible thanks to an adaptive Markov chain Monte Carlo (MCMC) algorithm called T-Walk (Christen & Fox 2010), which is uniquely suited to efficient exploration of multimodal distributions. The T-Walk algorithm uses a Metropolis–Hastings step, but at each step, certain chains perform a "jump" over the current point in another, randomly chosen chain in an attempt to find new modes in the posterior distribution. In addition, linear coordinate transformations are performed in order to better handle parameter degeneracies. While a substantial number of chains may be required to confidently sample the parameter space, we find that the T-Walk algorithm nevertheless converges much faster than a nested sampling algorithm even if there are many separate modes in the parameter space. One important benefit is that if there are multiple subhalo perturbations in the lensed image, the MCMC can find both modes even if only one subhalo is included in the model. One disadvantage of using T-Walk, however, is that it is not necessarily clear ahead of time how many chains will be needed to sample all the modes in the parameter space—for this reason, along with the reasons stated above, a nonparametric analysis to identify lens perturbations in actual data is still recommended before proceeding to the MCMC.

To investigate the bias in the inferred subhalo mass in each scenario, in Table 3 we show the "true" total subhalo mass and the best-fit subhalo mass for each scenario. As expected, the total subhalo mass is well recovered in Simulation 1, since the subhalo density profile matches the profile of the model. In Simulations 2 and 3, where the slope of the density profile differs but the subhalo distance matches the model (${r}_{0}=b$), the error is still fairly small, roughly 11% in each case. This is a fortunate consequence of the fact that the scale of the subhalo lensing perturbation is naturally similar to the tidal radius under this assumption, as discussed in Section 2. However, the situation is very different if the tidal radius differs significantly from the model, as is evident in the Simulation 4 results: here, the actual subhalo mass is nearly five times higher than the inferred mass, with an error of nearly 80%. The reason is clear: in this case, the subhalo tidal radius is approximately 2.6 arcsec, compared to 0.27 arcsec in the best-fit model. Thus, the subhalo density profile extends much farther out compared to the model profile, and beyond the critical curve perturbation radius $\delta {\theta }_{c}$, this extra mass can be "hidden" by adjusting the smooth lens component. This high systematic error in the total subhalo mass has also been pointed out by Vegetti et al. (2012) and Vegetti & Vogelsberger (2014) (indeed, they estimate this error by considering different projected distances to the lens plane, under the assumption that the subhalo number density traces the dark matter profile of the host galaxy).

Table 3.  Comparison of Subhalo Mass Estimators for the Actual vs. Best-fit Values in each Simulated Gravitational Lens (see Table 1 for a Descriptive Comparison of each Simulation)

	${m}_{T}/{10}^{9}\,{M}_{\odot }$	${m}_{T,\mathrm{fit}}/{10}^{9}\,{M}_{\odot }$	${\rm{\Delta }}{m}_{T}/{m}_{T}$	$\delta {\theta }_{c}$ (kpc)	${\tilde{m}}_{\mathrm{sub}}/{10}^{9}\,{M}_{\odot }$	${\tilde{m}}_{\mathrm{sub},\mathrm{fit}}/{10}^{9}\,{M}_{\odot }$	${\rm{\Delta }}{\tilde{m}}_{\mathrm{sub}}/{\tilde{m}}_{\mathrm{sub}}$
Simulation 1	1.027	1.058	3.02%	0.679	0.473	0.476	0.60%
Simulation 2	1.504	1.674	11.3%	0.716	0.494	0.518	4.87%
Simulation 3	2.341	2.100	10.3%	0.809	0.592	0.662	11.7%
Simulation 4	8.698	1.794	79.4%	0.744	0.519	0.556	7.12%

Note. m_T refers to the total subhalo mass, while the effective subhalo mass ${\tilde{m}}_{\mathrm{sub}}\equiv {m}_{\mathrm{sub}}(\delta {\theta }_{c})/\alpha$ is the projected subhalo mass enclosed within the subhalo perturbation radius $\delta {\theta }_{c}$ (which we list in parsecs in the middle column). Note that for Simulation 4, where the subhalo has a substantially larger tidal radius, the inferred total subhalo mass m_T is off by ∼80% while the effective subhalo lensing mass ${\tilde{m}}_{\mathrm{sub}}$ is recovered very well by comparison.

Download table as:  ASCIITypeset image

It has been suggested in Vegetti & Vogelsberger (2014) that the mass enclosed within the subhalo Einstein radius is robust to changes in the density profile, so we also test this idea on our simulated lenses. With the exception of Simulation 1 (where the subhalo density profile and tidal radius are correctly modeled), we find that at the subhalo Einstein radius b_s in the best-fit model (which is in the range 004–005; see Table 2), the enclosed subhalo mass is significantly overestimated. In the most extreme case, Simulation 3, the best-fit subhalo mass enclosed within its Einstein radius is ∼3.5 times higher than the actual mass, an error of 250%. It is also worth noting that the Einstein radius of the actual subhalo can also be dramatically different from the model if the slope of the density profile differs enough; in Simulation 3, the true Einstein radius of the subhalo is $7.3\times {10}^{-4}$ arcsec, much smaller than the Einstein radius of the best-fit model subhalo. Thus, we conclude that the Einstein radius of the subhalo itself plays no special role, and the inferred mass within this radius is not robust to changes in the density profile.

Having established that the bias in the inferred subhalo mass can be significant, we now investigate whether the effective subhalo lensing mass ${\tilde{m}}_{\mathrm{sub}}\equiv {m}_{\mathrm{sub}}(\delta {\theta }_{c})/\alpha$ (defined in Section 2) can indeed be considered a robust mass estimator. In Figure 7 we plot the subhalo mass profile of the best-fit model in each case (dashed line) divided by the log-slope α of the host galaxy, compared to the actual subhalo mass profile (solid line; we recall that $\alpha =1$ for the actual model in each case). The shaded region denotes the 95% confidence interval, which is obtained by calculating a derived posterior probability in the subhalo mass enclosed within each radius. The simulation recovers the mass profile very well, as expected; for the remaining cases, the mass enclosed within $\delta {\theta }_{c}$ (denoted by the dotted vertical line) divided by α is indeed quite close to the actual subhalo mass even though the mass elsewhere is different. In Table 3 the values for ${\tilde{m}}_{\mathrm{sub}}$ are tabulated along with the error. In Simulation 4, the error in this quantity is roughly 8%, a vast improvement compared to the discrepancy in the total subhalo mass. With the exception of the shallowest density profile $\gamma =1$ (Simulation 3), the error in ${\tilde{m}}_{\mathrm{sub}}$ is significantly smaller than the error in the total subhalo mass.

Zoom In Zoom Out Reset image size

Further evidence for the invariance of the effective subhalo mass ${\tilde{m}}_{\mathrm{sub}}$ can be seen in the posteriors themselves. In the parameter region where the data are well fit, if ${m}_{\mathrm{sub}}(\delta {\theta }_{c})/\alpha$ is invariant, then it follows that ${b}_{s}/\alpha$ must also be approximately invariant. Indeed, this can be seen in the approximate formula Equation (11): to highest order, $\delta {\theta }_{c}$ will remain the same if the ratio ${b}_{s}/\alpha$ is fixed. Thus we expect a linear degeneracy between the parameters b_s and α, which is shown beautifully in the joint posterior in these parameters in Figure 6. For the best-fit solution, we have ${b}_{s}\approx 0.05$, $\alpha \approx 1.4$, and thus the ratio ${b}_{s}/\alpha \approx 0.036$. If we pick a different point, for example where $\alpha =1.2$, then to maintain this ratio, we must therefore have ${b}_{s}\approx 0.036\times 1.2\approx 0.043$, and indeed this value falls within the allowed parameter region in the figure. The relationship is not satisfied quite as well for $\alpha =1$, where the posterior lies slightly above the expected value ${b}_{s}\approx 0.036;$ however, in this region, the subhalo position has drifted farther away from the critical curve, which can be seen in the joint α-${x}_{c,s}$ posterior at $\alpha =1$. Since the subhalo position has shifted, the effective subhalo lensing mass will change, and in order to produce the same critical curve perturbation, b_s needs to be slightly larger, as the posterior shows. Nevertheless, in the neighborhood of the best-fit region, it is clear that the α-b_s degeneracy is fully consistent with the invariance of the effective subhalo lensing mass ${\tilde{m}}_{\mathrm{sub}}$.

Although Figure 7 shows that the effective subhalo lensing mass ${\tilde{m}}_{\mathrm{sub}}$ is fairly robust, it is also biased slightly high in each case. Why is this? There are two reasons. First, although $\delta {\theta }_{c}$ in the best-fit model is close to the actual value in each case, it is not exactly the same—it differs either because the critical curve perturbation itself is larger or because the subhalo is slightly farther from the unperturbed critical curve compared to the actual subhalo. This is most evident in Simulation 3, where the best-fit subhalo position is offset by 0.03 arcsec to the right compared to the true position; this partly accounts for the larger $\delta {\theta }_{c}$ in the best-fit model for Simulation 3, since the subhalo is farther from the critical curve. This shifting of the subhalo away from the actual position of the critical curve is a compromise to avoid an overly "sharp" pinching of the critical curve (and a corresponding high magnification at the position of the subhalo) that results from the steeper density profile in the model subhalo. Likewise, the size of the critical curve perturbation itself may differ slightly as a compromise in order to better fit the perturbed surface brightness at different radii; given that the surface brightness variations will differ from lens to lens, the bias due to this is very difficult to quantify in general.

The second reason for the bias is due to differences in the smooth lens component. As is evident in Table 2, having a different log-slope α is not the only difference in the host galaxy parameters compared to their actual values. There is also a slightly different external shear ${{\rm{\Gamma }}}_{\mathrm{ext}}$ (specifically, the x-component ${{\rm{\Gamma }}}_{1}$) and axis ratio q, and additionally, the host galaxy is shifted significantly to the left. All of these change the left-hand side of Equation (1), and hence the subhalo mass enclosed within $\delta {\theta }_{c}$ is not strictly invariant. Part of this can be specifically accounted for: in Appendix B it is shown that assuming the subhalo position is well reproduced, the invariant quantity is actually ${m}_{\mathrm{sub}}(\delta {\theta }_{c})/\alpha \eta \,={\tilde{m}}_{\mathrm{sub}}/\eta$ (Equation (39)), where the parameter η is determined by the external shear ${{\rm{\Gamma }}}_{\mathrm{ext}}$ (Equation (34)). Thus, we should have ${\tilde{m}}_{\mathrm{sub},\mathrm{fit}}\approx {m}_{\mathrm{sub},\mathrm{true}}\tfrac{{\eta }_{\mathrm{fit}}}{{\eta }_{\mathrm{true}}}$. In this case, the direction of the subhalo ϕ and the direction of the external shear perturber ${\phi }_{p}$ differ by 90°, so we have $\eta \approx 1+{{\rm{\Gamma }}}_{\mathrm{ext}}\approx 1+{{\rm{\Gamma }}}_{1}$. Thus, using the best-fit values from Table 2, we have ${\eta }_{\mathrm{true}}\approx 1.036$, while ${\eta }_{\mathrm{fit}}\approx 1.07$ for Simulations 2–4. Thus, ${\tilde{m}}_{\mathrm{sub}}$ is inflated by a factor of $1.07/1.036\approx 1.03$, or about 3% in each case. While this does not account for all of the error in ${\tilde{m}}_{\mathrm{sub}}$, we can see in the final column of Table 3 that it accounts for roughly half the error in Simulations 2 and 4 (although it accounts for only a quarter of the error in Simulation 3; much of the remaining error is likely due to the offset position of the subhalo). Unfortunately, the factor η does not allow us to improve our mass estimator, because there is no way of knowing the true external shear (and hence ${\eta }_{\mathrm{true}}$), but it at least allows us to quantify the error due to this effect in the simulated cases.

Even if it were possible to have prior knowledge of the true parameters of the smooth lens component, fixing these to their correct values does not improve the fit; indeed, it makes it worse. The more freedom the smooth lens component has to partially mimic the effect of the subhalo beyond $\delta {\theta }_{c}$, the better the subhalo model will be able to reproduce the subhalo perturbation radius and hence the (scaled) mass within $\delta {\theta }_{c}$. Thus, having sufficient freedom to vary the smooth lens component is essential to achieving a good fit, even if a small amount of bias inevitably results in these parameters as a result of the mismatched subhalo density profile. In fact, it may be beneficial to include additional components to the smooth lens model beyond what we have done here—in particular, higher-order multipole moments might allow for a better fit in the region of the subhalo, but we have not explored this option in this work.

While we have demonstrated in Section 5 that our subhalo mass estimator ${\tilde{m}}_{\mathrm{sub}}$ can be determined robustly, we have seen in Section 6 that this condition depends on achieving a good fit, such that the subhalo model is favored by the Bayesian evidence. In particular, the fit must be good enough so that the inferred position of the subhalo is reasonably accurate (such that the error in the subhalo position is small compared to the subhalo distance to the critical curve, $\delta {\theta }_{c}$). One can see from Table 2 that this is true for each of the simulations analyzed in this paper. However, in Simulation 3 the subhalo position has the greatest error, with the best-fit position shifted by nearly 0.03 arcsec farther from the critical curve. This is a consequence of the true subhalo density profile having a much shallower slope ($\rho \sim {r}^{-1}$) compared to the model subhalo ($\rho \sim {r}^{-2}$). If the model subhalo is too close to the critical curve, it "pinches" the critical curve too much (we recall the difference in perturbation shape from Figure 3), resulting in a perturbation that cannot fit the data well, so a compromise is reached where the best-fit subhalo is sightly farther out. Indeed, by comparing the Bayesian evidence ${ \mathcal E }$ in the last column of Table 2, one sees that of the four simulations we modeled, Simulation 3 gave the poorest fit, and also gave the greatest error in the effective lensing mass (Table 3).

Is it possible to have a realistic scenario where the density profile and tidal radius of the subhalo differ so dramatically from the model that a good fit cannot be achieved at all? While we cannot conclude about this yet, it is conceivable that a typical NFW subhalo may fall in this category. In this paper we varied the suhalo log-slope and tidal radius separately to isolate each effect. However, if the true subhalo has an NFW profile, it can be expected to exhibit a combination of these two effects—that is, a shallow inner log-slope of −1, plus a larger tidal radius since we expect ${r}_{0}\gt b$. In this case, the discrepancy in the inferred total subhalo mass will be even more profound than it was for Simulation 4, since the density profile would not fall off as quickly beyond $\delta {\theta }_{c}$ and more mass would lie outside the subhalo perturbation radius. More worringly, the perturbation of the neighboring images differs so much compared to the pseudo-Jaffe model that it is questionable whether such a scenario could even be well fit by a pseudo-Jaffe subhalo at all. If so, the implication would be that the substructures detected thus far in gravitational lenses cannot correspond to NFW subhalos with realistic tidal radii. However, as we have not attempted to fit simulated data with an NFW subhalo in this paper, this remains an open question at present.

The effect we have seen on the log-slope of the host galaxy α is rather striking, and it is tempting to interpret an $\alpha \gt 1$ result from an actual data set as evidence that the subhalo density profile is indeed shallower than isothermal or that the tidal radius is underestimated. In Vegetti et al. (2010), the inferred log-slope $\alpha \approx 1.2$, which could be consistent with the scenarios discussed here. On the other hand, Vegetti et al. (2012) infer $\alpha \approx 1.05$, while for SDP.81, Hezaveh et al. (2016) find $\alpha \approx 1.06;$ in both cases, the slope of the host galaxy profile is just barely consistent with isothermal. If taken at face value, this might suggest that the subhalo density profile cannot be much shallower than isothermal. However, we must be cautious when interpreting these results. It may be the case that the host galaxy density profile may have a log-slope somewhat less than 1, in which case these results could still be consistent with a shallower subhalo profile. Indeed, Schneider & Sluse (2013) show that the log-slope of the host galaxy in the vicinity of the images can be expected to be biased due to the mass-sheet degeneracy if the galaxy profile differs from a strict power law at other radii. Furthermore, even if the true α is very close to 1 in these cases, the parameter-space solution may not necessarily be unique. As is clear in Figure 6, $\alpha =1$ may still be a viable region of parameter space; if α is fixed to 1, the subhalo may still be found in that parameter space even if the fit is poorer compared to having a steeper α. Thus, it is possible that in these cases, a mode with higher α might be found by a more complete multimodal sampling of the parameter space, e.g., using the T-Walk algorithm, and would produce an even better fit.

Ultimately, to minimize the bias in both the inferred subhalo properties and the host galaxy parameters, the ideal approach is to vary both the subhalo tidal radius r_t and inner log-slope γ as free parameters. However, while it remains to be seen how well a subhalo density profile can be constrained, e.g., in high-resolution ALMA images of submillimeter gravitational lenses, it seems likely that the tidal radius will be prior dominated given that the subhalo pertubation beyond $\delta {\theta }_{c}$ can be mimicked by perturbing the smooth lens component, and hence the exact tidal radius will be difficult to constrain. Given the difficulty in constraining these parameters (together with the added computational expense of varying them), we therefore argue that the effective subhalo lensing mass ${m}_{\mathrm{sub}}(\delta {\theta }_{c})/\alpha$ is a very useful quantity given that it is largely robust to changes in the subhalo density profile and tidal radius. Given several subhalo detections, one could use these mass measurements to test whether they are consistent with the expected mass function from ΛCDM, or whether an alternate dark matter scenario (e.g., WDM) would provide a better fit.

In this section, we provide a summary of procedures for calculating the subhalo perturbation radius $\delta {\theta }_{c}$ for a given lens model at different levels of approximation.

If the model subhalo is isothermal/pseudo-Jaffe and one is content with approximating $\delta {\theta }_{c}$ to 20%–30% accuracy, we recommend using Equation (11) because it is relatively easy to use. Despite the approximation, the subhalo mass within this radius is still quite robust (and in fact works slightly better than the exact values in our Simulations 2–4). For lens modelers who wish to calculate $\delta {\theta }_{c}$ to better accuracy, the more rigorous Equation (47) can be used, or Equation (1) can be solved numerically.

To solve for $\delta {\theta }_{c}$ numerically using Equation (1) (which works regardless of the model being used for the subhalo profile), one could adopt the following method. Calculate the direction ${\phi }_{{\rm{\Gamma }}}$ of the unperturbed shear ${{\boldsymbol{\Gamma }}}_{\mathrm{tot}}$ at the position of the subhalo. The solution to Equation (1) will lie somewhere along the line perpendicular to this direction, i.e., along ${\phi }_{r}={\phi }_{{\rm{\Gamma }}}+90^\circ$, because along this line the subhalo shear aligns with the shear of the unperturbed lens model. In this way, Equation (1) is reduced to an equation in r, which can be found with a root-finding algorithm.

To get a sense of the relative accuracy of these different methods, in Table 4 in Appendix B we list a comparison of the $\delta {\theta }_{c}$ value calculated numerically using this procedure, and the values calculated using both Equation (11) and the more exact Equation (47) for the four simulations considered in this paper.

Table 4.  Comparison between Different Solutions for the Subhalo Perturbation Radius $\delta {\theta }_{c}$ of the Best-fit Model in Each Simulated Lens Scenario (See Table 1 for a Description of Each Simulation).

	$\delta {\theta }_{c}$ (exact)	$\delta {\theta }_{c}$ (Equation (47))	$\delta {\theta }_{c}$ (approx, Equation (11))
Sim. 1	0149	0147	0185
Sim. 2	0155	0151	0188
Sim. 3	0175	0171	0202
Sim. 4	0161	0157	0193

Note. We list the numerical solution (first column), analytic solution Equation (47) (second column), and approximate analytic solution Equation (11) (third column).

Download table as:  ASCIITypeset image

Here we obtain a more rigorous result for the subhalo perturbation radius for the case where the subhalo is not exactly on a critical curve and both ellipticity and external shear are present. Suppose the subhalo is located at a radius ${\theta }_{s}$ (not necessarily on the critical curve) and angle ϕ with respect to the x-axis, and the host galaxy has an axis ratio q and major axis oriented at an angle ${\phi }_{0}$ with respect to the x-axis. Suppose also that there is an external shear ${{\rm{\Gamma }}}_{\mathrm{ext}}$ created by a distant perturbing mass whose direction is denoted by angle ${\phi }_{p}$.

The projected density profile of the host galaxy is modeled as a power-law ellipsoid. In a coordinate system ($x^{\prime} ,y^{\prime}$) where the major axis of the lens coincides with the $x^{\prime}$-axis, we write the profile as follows, using the notation ${\boldsymbol{\theta }}^{\prime} =(x^{\prime} ,y^{\prime} )$:

$$\begin{eqnarray}&&{\kappa }_{0}({\boldsymbol{\theta }}^{\prime} )=\displaystyle \frac{2-\alpha }{2}{b}^{\alpha }{({qx}{{\prime} }^{2}+y{{\prime} }^{2}/q)}^{-\tfrac{\alpha }{2}}.\end{eqnarray} \tag{ 25 }$$

In this formulation, the normalization parameter b can be thought of as the average Einstein radius, while the semimajor axis of the critical curve is $b/\sqrt{q}$. Another common choice is to make the semimajor axis of the critical curve the normalization parameter (call it $b^{\prime}$), in which case one should substitute $b^{\prime} =b/\sqrt{q}$ in the following formulas. In any event, in the actual coordinate system of the lens plane ($x,y$), we calculate ${\kappa }_{0}({\boldsymbol{\theta }})$ by rotating these coordinates by the angle ${\phi }_{0};$ in other words, by first calculating ${\theta }_{i}^{\prime} ={{ \mathcal R }}_{{ij}}{\theta }_{j}$, where ${{ \mathcal R }}_{{ij}}$ is the corresponding rotation matrix, and inserting into Equation (25).

We can write this in terms of polar coordinates in the lens frame $(\theta ,\phi )$ and separate out the angular dependence by rewriting it as follows:

$$\begin{eqnarray}&&{\kappa }_{0}({\boldsymbol{\theta }})=\displaystyle \frac{2-\alpha }{2}{\left(\displaystyle \frac{b\sigma }{\theta }\right)}^{\alpha },\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray}&&\sigma \equiv {\left(q{\cos }^{2}(\phi -{\phi }_{0})+\displaystyle \frac{1}{q}{\sin }^{2}(\phi -{\phi }_{0})\right)}^{-\tfrac{1}{2}}.\end{eqnarray} \tag{ 27 }$$

To solve Equation (24), we need an expression for the shear from the host galaxy ${{\rm{\Gamma }}}_{0}$. For an isothermal ellipsoid ($\alpha =1$), the magnitude of the shear is equal to ${\kappa }_{0}$ at that location, so that ${\kappa }_{0}+{{\rm{\Gamma }}}_{0}=2{\kappa }_{0}$. For other values of α, this is no longer true, but it can be shown that the ratio $g\equiv {{\rm{\Gamma }}}_{0}/{\kappa }_{0}$ depends only on the angle $\phi -{\phi }_{0};$ for a fixed angle, g is the same at all radii. This ratio can be found numerically using gravitational lensing software, but it can also be derived analytically using the results of Tessore & Metcalf (2015), where a closed-form solution for the deflection of a power-law ellipsoid lens is derived. The result is

$$\begin{eqnarray}g & = & \left|1-\displaystyle \frac{4(1-\alpha )\sqrt{q}}{(1+q)(2-\alpha )\sigma }{e}_{2}^{i({\rm{\Delta }}\phi ^{\prime} -{\rm{\Delta }}\phi )}{F}_{1}\right.\\ & & \left.\times \left(1;\displaystyle \frac{\alpha }{2};2-\displaystyle \frac{\alpha }{2};\displaystyle \frac{q-1}{q+1}{e}^{i2{\rm{\Delta }}\phi }\right)\right|,\end{eqnarray} \tag{ 28 }$$

where ${}_{2}{F}_{1}(...)$ is the complex Gaussian hypergeometric function, and ${\rm{\Delta }}\phi ^{\prime}$ is defined by $\tan {\rm{\Delta }}\phi ^{\prime} =\tfrac{1}{q}\tan {\rm{\Delta }}\phi$, where it is understood that ${\rm{\Delta }}\phi ^{\prime}$ must be in the same quadrant as ${\rm{\Delta }}\phi$. Note that in the isothermal case $\alpha =1$, Equation (28) reduces to g = 1, as expected.

If the axis ratio q = 1, Equation (28) reduces to the much simpler expression $g=\alpha /(2-\alpha )$, which is consistent with the fact that ${{\rm{\Gamma }}}_{0}={\bar{\kappa }}_{0}-{\kappa }_{0}$ when the host galaxy is axisymmetric. This expression may be used as a rough approximation for g in the regime $q\gt 0.8$, $0.8\lt \alpha \lt 1.2$ or so with only a few percent loss of accuracy in $\delta {\theta }_{c};$ for q and α outside this range, however, the error can become severe, so in that case, we recommend using Equation (28) or else calculating g using gravitational lensing software.

The direction of the shear is also important, since this determines the direction of maximum warping of the critical curve. In the $\alpha =1$ case, it can be shown algebraically that the shear from an isothermal ellipsoid is still along the tangential direction, even for low qvalues (and even though the critical curve is not along the tangential direction!). Thus, in the absence of external shear, the direction of maximum warping of the critical curve by a subhalo is still along the radial direction for an isothermal galaxy. For $\alpha \ne 1$, this is no longer strictly true; however, we find that as long as $\alpha \gt 0.7$, the difference compared to the radial direction is no more than 5°. The presence of external shear can also affect the shear direction, however. At a point given by polar coordinates $(\theta ,\phi )$, the total shear magnitude and direction are given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{tot}}\approx \sqrt{{g}^{2}{\kappa }_{0}^{2}+{{\rm{\Gamma }}}_{\mathrm{ext}}^{2}+2g{\kappa }_{0}{{\rm{\Gamma }}}_{\mathrm{ext}}\cos 2(\phi -{\phi }_{p})},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\phi }_{{\rm{\Gamma }}}\approx \displaystyle \frac{\pi }{2}+\displaystyle \frac{1}{2}{\sin }^{-1}\left(\displaystyle \frac{g{\kappa }_{0}\sin 2\phi +{{\rm{\Gamma }}}_{\mathrm{ext}}\sin 2{\phi }_{p}}{{{\rm{\Gamma }}}_{\mathrm{tot}}}\right),\end{eqnarray} \tag{ 30 }$$

where we have made the approximation that the shear of the host galaxy is along the radial direction (which is exactly true if $\alpha =1$). It should be emphasized that ${\phi }_{p}$ is the direction of the (distant) perturbing galaxy producing the external shear, not the direction of the external shear itself (which differs from ${\phi }_{p}$ by 90°). Note that for ${{\rm{\Gamma }}}_{\mathrm{ext}}=0$, the shear angle ${\phi }_{{\rm{\Gamma }},0}=\phi +\tfrac{\pi }{2}$, i.e., the shear is tangential in that case. However, the difference in the shear angle due to the external shear, ${\rm{\Delta }}{\phi }_{{\rm{\Gamma }}}={\phi }_{{\rm{\Gamma }}}-{\phi }_{{\rm{\Gamma }},0}$, is fairly small provided that ${{\rm{\Gamma }}}_{\mathrm{ext}}\lesssim 0.2;$ in this case, we find that typically ${\rm{\Delta }}{\phi }_{{\rm{\Gamma }}}\lesssim 10^\circ$, and the corresponding change in $\delta {\theta }_{c}$ due to this direction change was within a few percent. Therefore, in the following derivation we approximate the direction of $\delta {\theta }_{c}$ to still be along the radial direction, given that the effect of the perturbed shear direction is rather small.

To derive the formula for $\delta {\theta }_{c}$, we return to Equation (24) and rewrite it as follows:

$$\begin{eqnarray}&&\eta -(1+g){\kappa }_{0}({{\boldsymbol{\theta }}}_{s}+{\boldsymbol{r}})={\bar{\kappa }}_{s}(r),\end{eqnarray} \tag{ 31 }$$

where we have defined

$$\begin{eqnarray}\eta & \equiv & 1+{{\rm{\Gamma }}}_{0}-{{\rm{\Gamma }}}_{\mathrm{tot}}=1+g{\kappa }_{0}-g{\kappa }_{0}\\ & & \times \sqrt{1+{\left(\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{ext}}}{g{\kappa }_{0}}\right)}^{2}+2\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{ext}}}{g{\kappa }_{0}}\cos 2(\phi -{\phi }_{p})}.\end{eqnarray} \tag{ 32 }$$

As long as the external shear is not too small (${{\rm{\Gamma }}}_{\mathrm{ext}}\lesssim 0.3$ or so), η varies only a miniscule amount as r is varied (in fact, it is easy to show that if we choose the direction $\phi ={\phi }_{p}$ or $\phi ={\phi }_{p}+\tfrac{\pi }{2}$, it does not vary at all with r; only for directions unaligned with the external shear is there a slight variation). Because it varies so little, a very good approximation can be found by fixing η to its value at the (unperturbed) critical curve along the direction of the subhalo. Using the fact that on the unperturbed critical curve we have $1-{\kappa }_{0}-{{\rm{\Gamma }}}_{\mathrm{tot}}=0$, after some tedious algebra, we derive

$$\begin{eqnarray}\eta & \approx & \displaystyle \frac{1+g{{\rm{\Gamma }}}_{\mathrm{ext}}\cos 2(\phi -{\phi }_{p})}{g-1}\\ & & \times \left\{-1+\sqrt{1+\displaystyle \frac{({g}^{2}-1)(1-{{\rm{\Gamma }}}_{\mathrm{ext}}^{2})}{{(1+g{{\rm{\Gamma }}}_{\mathrm{ext}}\cos 2(\phi -{\phi }_{p}))}^{2}}}\right\}.\end{eqnarray} \tag{ 33 }$$

In general, this expression increases extremely slowly with g, so a very good approximation is achieved by taking the limiting case $g\to 1$ (which corresponds to the isothermal case):

$$\begin{eqnarray}&&\eta \approx \displaystyle \frac{1-{{\rm{\Gamma }}}_{\mathrm{ext}}^{2}}{1+{{\rm{\Gamma }}}_{\mathrm{ext}}\cos 2(\phi -{\phi }_{p})}.\end{eqnarray} \tag{ 34 }$$

In order for Equation (31) to be analytically tractable, we need to simplify the expression for ${\kappa }_{0}$ (Equation (26)). For an isothermal ellipsoid, the magnitude of the deflection (which is proportional to ${\kappa }_{0}\theta$) is independent of radius; while this is not strictly true for $\alpha \ne 1$, ${\kappa }_{0}\theta$ will nevertheless vary slowly with radius provided the slope α is not too extreme. Therefore, we expand ${\kappa }_{0}\theta$ to first order in θ around the position of the subhalo ${\theta }_{s}$, which gives (after dividing by θ)

$$\begin{eqnarray}&&{\kappa }_{0}({{\boldsymbol{\theta }}}_{s}+{\boldsymbol{r}})\approx {\kappa }_{0}({{\boldsymbol{\theta }}}_{s})\left[1-\alpha \displaystyle \frac{r}{{\theta }_{s}+r}\right].\end{eqnarray} \tag{ 35 }$$

Plugging this into Equation (31) results in the following equation, which can now be solved analytically for an assumed subhalo density profile:

$$\begin{eqnarray}&&\eta -(1+g){\kappa }_{0}({{\boldsymbol{\theta }}}_{s})+\displaystyle \frac{\alpha (1+g){\kappa }_{0}({{\boldsymbol{\theta }}}_{s})r}{{\theta }_{s}+r}={\bar{\kappa }}_{s}(r).\end{eqnarray} \tag{ 36 }$$

Before assuming a density profile for the subhalo, we note that in the absence of the subhalo (${\bar{\kappa }}_{s}=0$), Equation (31) tells us that the position of the (unperturbed) critical curve ${\theta }_{c}$ along the direction ϕ is given by $\eta =(1+g){\kappa }_{0}({\theta }_{c})$. Since the subhalo will be close to the unperturbed critical curve, we can expand $(1+g){\kappa }_{0}({{\boldsymbol{\theta }}}_{s})$ around ${\theta }_{c}$, which to second order is

$$\begin{eqnarray}&&(1+g){\kappa }_{0}({{\boldsymbol{\theta }}}_{s})\approx \eta \left[1-\alpha \displaystyle \frac{{\rm{\Delta }}{\theta }_{s}}{{\theta }_{c}}+\displaystyle \frac{\alpha (\alpha +1)}{2}{\left(\displaystyle \frac{{\rm{\Delta }}{\theta }_{s}}{{\theta }_{c}}\right)}^{2}\right],\end{eqnarray} \tag{ 37 }$$

where ${\rm{\Delta }}{\theta }_{s}={\theta }_{s}-{\theta }_{c}$. For the moment, we keep only first-order terms in ${\rm{\Delta }}{\theta }_{s}/{\theta }_{c}$ and $r/{\theta }_{s}$. Combining Equations (37) and (36) and keeping first-order terms, we find

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{\theta }_{s}}{{\theta }_{c}}+\displaystyle \frac{r}{{\theta }_{s}}\approx \displaystyle \frac{{\bar{\kappa }}_{s}(r)}{\alpha \eta }.\end{eqnarray} \tag{ 38 }$$

This result is very important: it implies that if a lens model reproduces the correct subhalo perturbation radius $r=\delta {\theta }_{c}$ (which is the solution to the above equation), and if the best-fit position of the subhalo is accurate, the quantity ${\bar{\kappa }}_{s}(r)/\alpha \eta$ will be well-reproduced regardless of the mass distribution of the subhalo, and even if the best-fit α and external shear ${{\rm{\Gamma }}}_{\mathrm{ext}}$ differ from their actual values. Since the average kappa is proportional to the projected subhalo mass m_sub(r) enclosed within a given radius, we can say that

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{m}_{\mathrm{sub}}(\delta {\theta }_{c})}{\alpha \eta }\right)}_{\mathrm{fit}}\approx {\left(\displaystyle \frac{{m}_{\mathrm{sub}}(\delta {\theta }_{c})}{\alpha \eta }\right)}_{\mathrm{true}}.\end{eqnarray} \tag{ 39 }$$

In general, after fitting a lens, it is impossible to know to what extent the true value for the external shear differs from the best-fit value. However, since $\eta \approx 1$ regardless, in practice, this affects Equation (39) only slightly, typically by not more than a few percent. The effect of α can be more dramatic, but since most gravitational lens galaxies are well fit by an isothermal profile, in practice, one can surmise that ${\alpha }_{\mathrm{true}}\approx 1$ and find the mass enclosed within the subhalo under this assumption. Thus, we define the effective subhalo lensing mass ${\tilde{m}}_{\mathrm{sub}}$ to be the projected mass enclosed within the subhalo perturbation radius, divided by α; in other words,

$$\begin{eqnarray}&&{\tilde{m}}_{\mathrm{sub}}\equiv {m}_{\mathrm{sub}}(\delta {\theta }_{c})/\alpha .\end{eqnarray} \tag{ 40 }$$

In Section 5 we show that this mass estimator is robust even if the assumed density profile and/or tidal radius of the subhalo are inaccurate.

Now that we have established the usefulness of the subhalo perturbation radius $\delta {\theta }_{c}$, we assume the subhalo has a pseudo-Jaffe density profile with a tidal radius r_t, for which we have

$$\begin{eqnarray}&&{\bar{\kappa }}_{s}(r)=\displaystyle \frac{{b}_{s}}{{r}^{2}}(r+{r}_{t}-\sqrt{{r}_{t}^{2}+{r}^{2}}).\end{eqnarray} \tag{ 41 }$$

In order to render the equation analytically tractable, we Taylor expand the deflection $r\bar{\kappa }$ to first order. Since a rough approximation for $\delta {\theta }_{c}$ is given by $r\approx \delta {\theta }_{c}\approx \sqrt{{{bb}}_{s}}$, we expand around this point. (The approximate solution we expand around can be improved upon, e.g., by taking the highest-order term in Equation (11), which typically leads to an extra 1% percent accuracy in $\delta {\theta }_{c};$ however, in the interest of simplicity, we choose $\sqrt{{{bb}}_{s}}$ instead.) The result is

$$\begin{eqnarray}&&{\bar{\kappa }}_{s}\,\approx \,\epsilon \displaystyle \frac{{b}_{s}}{r}-\xi \sqrt{{b}_{s}/b},\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray}&&\epsilon =1+2\beta -2{(1+{\beta }^{2})}^{\tfrac{1}{2}}+{(1+{\beta }^{2})}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\xi =\beta -{(1+{\beta }^{2})}^{\tfrac{1}{2}}+{(1+{\beta }^{2})}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 44 }$$

and β is defined as the ratio of the subhalo distance r₀ to the center of the host galaxy over the Einstein radius,

$$\begin{eqnarray}&&\beta \equiv \displaystyle \frac{{r}_{0}}{b}=\displaystyle \frac{{r}_{t}}{\sqrt{{{bb}}_{s}}}.\end{eqnarray} \tag{ 45 }$$

The latter equality is given by the formula for the Jacobi radius of an isothermal subhalo. For many subhalos, we can expect that $\beta \gg 1$, and so a reasonable approximation in many cases may be obtained by taking $\beta \to \infty$, which gives $\epsilon \approx 1$, $\xi \approx 0$.

As we discussed in Section 3, it has been customary to assume that ${r}_{0}\approx b$ (so that $\beta \approx 1$), based on the fact that the projected distance of the subhalo to the host galaxy center is comparable to the Einstein radius. We emphasize that this assumption is unlikely to hold because the Einstein radius is typically in the range of 5–10 kpc and a subhalo would be unlikely to survive the severe tidal stripping that would result from being so close to the host galaxy. However, for the sake of comparison to previous work, we consider the results using this assumption, in which case ${r}_{t}\approx \sqrt{{{bb}}_{s}}$, and we have

$$\begin{eqnarray}&&\epsilon =3\left(1-\displaystyle \frac{1}{\sqrt{2}}\right),\,\xi =\displaystyle \frac{\epsilon }{3}\,(\mathrm{for}\,{r}_{0}=\mathrm{Einstein}\ \mathrm{radius}).\end{eqnarray} \tag{ 46 }$$

Substituting Equations (42) and (37) into Equation (36) and solving the resulting quadratic equation in r, we find the following solution:

$$\begin{eqnarray}&&\delta {\theta }_{c}\approx \displaystyle \frac{1}{\mu }(\sqrt{{b}_{s}{\theta }_{s}\epsilon \mu +{(\zeta /2)}^{2}}+\zeta /2),\end{eqnarray} \tag{ 47 }$$

where

$$\begin{eqnarray}&&\mu =\alpha \eta \left[1+\displaystyle \frac{\xi }{\alpha \eta }\sqrt{\displaystyle \frac{{b}_{s}}{b}}+(1-\alpha )\psi \right]\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&\zeta ={b}_{s}\epsilon -\xi {\theta }_{s}\sqrt{\displaystyle \frac{{b}_{s}}{b}}-\alpha \eta {\theta }_{s}\psi ,\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&\psi =\displaystyle \frac{{\rm{\Delta }}{\theta }_{s}}{{\theta }_{c}}\left(1-\displaystyle \frac{1}{2}(\alpha +1)\displaystyle \frac{{\rm{\Delta }}{\theta }_{s}}{{\theta }_{c}}\right).\end{eqnarray} \tag{ 50 }$$

To use these formulae, we need to calculate the position of the unperturbed critical curve ${\theta }_{c}$ along the direction of the subhalo. This can be found numerically using gravitational lensing software, in which case, there is no need to calculate g at all since it does not appear explicitly in Equations (47)–(50). Alternatively, ${\theta }_{c}$ can be found analytically by solving the equation $\eta =(1+g){\kappa }_{0}({{\boldsymbol{\theta }}}_{s})$. By substituting Equation (26), one can show that

$$\begin{eqnarray}&&{\theta }_{c}=b\sigma {\left[\displaystyle \frac{1}{2\eta }(1+g)(2-\alpha )\right]}^{1/\alpha }.\end{eqnarray} \tag{ 51 }$$

We have tested Equation (47) in many different scenarios and find that it is accurate to within 2%, provided the external shear is not too large (${{\rm{\Gamma }}}_{\mathrm{ext}}\lesssim 0.2$). In order to test our formula on the simulated lenses we modeled in Section 5, in Table 4 we list the numerical solution (first column) for $\delta {\theta }_{c}$ along with the values given by Equation (47) (second column) for the best-fit model in each simulation.

Since Equation (47) is still cumbersome to use for non-lens modelers, it is useful to have a simpler approximate expression. If we approximate ${\rm{\Delta }}{\theta }_{s}\approx 0$ (and hence $\psi \approx 0$) in the above equations and expand to first order in $\sqrt{{b}_{s}/b}$, we find

$$\begin{eqnarray}&&\delta {\theta }_{c}\approx \sqrt{\displaystyle \frac{\epsilon {\theta }_{s}{b}_{s}}{\alpha \eta }}\left(1-\displaystyle \frac{x}{2}+\displaystyle \frac{{x}^{2}}{8}\right)+\displaystyle \frac{\epsilon {b}_{s}}{2\alpha \eta }(1-x),\end{eqnarray} \tag{ 52 }$$

where

$$\begin{eqnarray}&&x=\xi \sqrt{\displaystyle \frac{{\theta }_{s}}{b\epsilon \alpha \eta }}.\end{eqnarray} \tag{ 53 }$$

One may set $\eta \approx 1$ in Equations (52) and (53) with at most a few percent loss of accuracy, so we reproduce this formula in Section 2 with η set to one (Equation (11)). If ${\rm{\Delta }}{\theta }_{s}$ is indeed zero, we find that this formula (with η set to 1) is accurate to within ∼1% provided the external shear is not too large (${{\rm{\Gamma }}}_{\mathrm{ext}}\lesssim 0.1$). The farther away the subhalo is from the critical curve, the worse this approximation becomes—generally, it will give a value for $\delta {\theta }_{c}$ that is too high. However, in the simulated lenses we examine in Section 5, the optimal radius where the effective subhalo lensing mass is invariant is slightly larger than the actual value for $\delta {\theta }_{c}$, and this may be true in general. One can see from Figure 7 that Equation (52) works better than using the actual $\delta {\theta }_{c}$, even though it is greater by about 17%–25%, depending on the simulation. The values given by this approximation are listed in the third column of Table 4 for the best-fit model in each simulation considered in Section 5.

Incidentally, one may verify from Equation (52) that if the log-slope of the host galaxy α and external shear η are varied, then to obtain the same $\delta {\theta }_{c}$ (to highest order), the Einstein radius of the subhalo b_s must be scaled by a factor $1/\alpha \eta ;$ in other words, to keep the same perturbation scale, ${b}_{s}/\alpha \eta$ must be invariant. This is another manifestation of the fact that ${\bar{\kappa }}_{s}(\delta {\theta }_{c})/\alpha \eta$ is invariant, as follows from Equation (39).