© The Authors 2023

1. Introduction

Gravitational lensing allows for a direct measure of the mass of galaxies, whether or not that mass is visible, and as such it is a valuable tool to study galaxies and dark matter (DM). Lensing can also be a tool to study cosmology, because time delays between multiple images can be compared to those predicted from a lens model to constrain a time-delay distance D_Δt ∝ 1/H₀. These tools are only as accurate as the galaxy mass distribution models on which they depend (see reviews by e.g., Birrer et al. 2022; Shajib et al. 2022).

Galaxy lenses are generally massive elliptical galaxies that are often described using multiple mass components, with one tracing the light representing the baryon mass distribution and the other representing the DM mass distribution. Dark matter mass distributions are typically described using an Navarro Frenk White (NFW) profile (Navarro et al. 1996), which takes the following form, parameterized in 3D in terms of a characteristic density ρ_s and scale radius r_s:

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_{\mathrm{s}}}{(r/r_{\mathrm{s}}){(1+r/r_{\mathrm{s}})}^2}·\end{array}$$(1)

While the spherical 3D profile has a clean analytical form, the 2D-projected elliptical NFW profile, which is necessary for lensing applications, does not have an analytical form for the deflection angle or lensing potential. To circumvent this problem, Golse & Kneib (2002) showed that a general elliptical lensing mass density profile can be expressed analytically if the ellipticity is introduced in the lensing potential rather than the mass distribution itself. This allows for analytical calculations of lensing properties, and so this parameterization has been the established method to model NFW profiles for various applications, including time-delay cosmography (e.g., Wong et al. 2020; Rusu et al. 2020; Shajib et al. 2022). Alternative methods to introduce ellipticity in the surface mass density while keeping the computation feasible include precalculating the numerical integrals on a grid to be interpolated (Schramm 1990; Schneider & Weiss 1991; Keeton 2001), or expanding the mass profile into a series of analytically tractable components, such as multi-Gaussian expansion (MGE; van de Ven et al. 2010; Shajib 2019) or the sum of cored steep ellipsoids (CSE; Oguri 2021). Because these methods introduce ellipticity directly in the mass distribution, they are thought to be more physically realistic. In this work we compare the Golse & Kneib (2002) representation of lens galaxies that captures the ellipticity in the lens potential, to the CSE-based method of Oguri (2021), which has the ellipticity in the mass-density. We refer to these descriptions as the NFWp and NFWm parameterizations, respectively.

This work aims at quantifying the impact of these different parametrizations on the predicted lensing properties, with a specific interest in galaxy-scale time delay cosmography. The main feature of interest is that while both parameterizations have the same radial profile (defined as the 1D profile of convergence within circular annuli), the azimuthal structure is different because of the different ways in which ellipticity is introduced. Kochanek (2021) shows the azimuthal freedom of a lens model must be considered with great care, which encourages us to be wary of the exact prescription of the azimuthal structure within lens models. Furthermore, Van de Vyvere et al. (2022a,b) show in individual lenses how unmodeled azimuthal structures can bias lens models, even though as a population these effects appear to average out. Our concern is that the choice of NFW parameterization could represent a systematic bias in the treatment of the azimuthal structure. As DM is a significant mass component in galaxies and as it is in many cases represented with a NFW profile, we focus here on comparing two different parametrizations of this profile and quantifying their impact on lens models. To quantify this comparison, we used the H₀ value recovered from the lensing model, as it directly reflects the capacity of a model to capture differences in the Fermat potential at the location of the images.

While the azimuthal structure of interest comes from the NFW prescription, we tested the effect of this prescription in the case where the profile has two components where only the DM component is represented by an NFW profile, such that the test is more applicable to the analysis of observed systems. A sensible method to evaluate the impact of the choice of parametrization of the NFW on H₀ consists in emulating and modeling strongly lensed systems that resemble known ones. Gomer et al. (2022, hereafter TDCVIII) created a population of mock lens images from analytical profiles designed to match the observed population of Time-Delay COSMOgraphy (TDCOSMO) lenses. The profiles used to create these lenses were two-component profiles, with a Chameleon profile (Dutton & Treu 2014) representing the light and an NFWp profile representing the DM. TDCVIII then fit these systems with a power-law elliptical mass distribution (PEMD; Barkana 1998) model with external shear. For this work, we created mocks that are identical to those used by TDCVIII, except that we changed the NFW parameterization to create a systematic comparison between the NFWp and NFWm implementations. We fit these mocks in the same manner as TDCVIII and compared the results of the parameter inferences.

We quantified the impact of these parameterizations in the context of the mass sheet transformation (MST), detailed in Sect. 2.2, where lenses with different radial profile shapes can give the same observables (i.e. image positions and fluxes). We made use of the MST-independent quantity ξ₂ (defined in Sect. 2.2) to diagnose the degree to which the degeneracies in this paper can be attributed to the MST. The primary goal of this work is to evaluate the degree to which the choice NFW parameterization can play a role in the determination of H₀ and identify the source of that role in the context of lensing degeneracies.

The paper is structured as follows: Section 2 reviews the NFW profile parameterizations as well as ξ, Sect. 3 compares the mock populations and the results of the PEMD fits, Sect. 4 discusses these results, Sect. 5 expands this discussion to its implications for other works, and Sect. 6 summarizes and concludes this work. Appendix A gives an analytical description of the axis ratio of the mass of the NFWp parameterization and Appendix B details several subtleties regarding the calculation or ξ in this work. This work adopts a fiducial flat cold dark matter (ΛCDM) cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_M = 0.3.

2. Formalism

In this section, we review the key differences between the NFWp and NFWm parameterizations of the density profile (Sect. 2.1), and discuss the quantity ξ which we propose to use as a diagnostic of degeneracies between models beyond the MST (Sect. 2.2).

2.1. NFW parameterization

Mass distributions of gravitational lenses are described in terms of surface mass density profiles normalized by the critical lensing density (i.e. convergence) κ(r) = Σ(r)/Σ_crit. If one wishes to use the NFW profile for lensing applications, one must be able to accommodate ellipticity and project into the 2D plane of the sky. The convergence for a circular NFW profile can be expressed as (Bartelmann 1996)

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{NFW}}=\frac{2{\kappa }_0}{u^2-1}[1-F(u)],\end{array}$$(2)

with scaled radius u = r/r_s and normalization κ₀ = ρ_sr_s/Σ_crit, and where

$$\begin{array}{c}\hfill F(u)=\{\begin{array}{cc}\frac{1}{\sqrt{1-u^2}}\mathrm{arctanh}\sqrt{1-u^2}\hfill & (u<1)\hfill \\ \frac{1}{\sqrt{u^2-1}}\arctan \sqrt{u^2-1}\hfill & (u>1).\hfill \end{array}\end{array}$$

The corresponding lensing potential can then be expressed as (Meneghetti et al. 2003)

$$\begin{array}{c}\hfill {\psi }_{\mathrm{NFW}}=2{\kappa }_0{\theta }_{\mathrm{s}}^{2}h(u),\end{array}$$(3)

where θ_s is the scale radius expressed in arcseconds and

$$\begin{array}{c}\hfill h(u)=\{\begin{array}{cc}{\ln }^2\frac{u}{2}-{\mathrm{arccosh}}^2\frac{1}{u}\hfill & (u<1)\hfill \\ {\ln }^2\frac{u}{2}+{\arccos }^2\frac{1}{u}\hfill & (u>1).\hfill \end{array}\end{array}$$

Ellipticity can be added to a given 1D profile by replacing the r argument with an elliptical radius such as $r_{\mathrm{ell}}=\sqrt{q{x}^2+{\mathit{y}}^2/q}$, where q is the axis ratio of the ellipse. For some profiles such as the Chameleon profile, this ellipticity is added directly to the convergence. However, such a substitution is not guaranteed to result in an analytical form for the lensing potential: the NFW profile has no such analytical form for the lensing potential arising from an elliptical 2D-projected mass distribution, instead requiring expensive numerical integrals to approximate.

The NFWp parameterization solves this problem by instead adding ellipticity analytically into the lensing potential, a strategy which can be applied generically for any analytical lensing potential profile. Golse & Kneib (2002) show that one can add ellipticity to the potential by replacing r with $r_{ϵ}=r_{\mathrm{ell}}\sqrt{2q_{\psi }/(1+q_{\psi }^2)}$ in Eq. (3), giving analytical expressions for the convergence and shear. To minimize confusion between ellipticity conventions, we have introduced q_ψ as the axis ratio of the potential. Similarly, we use q_κ for the convergence, and reserve q to a general context. Because of the analytical form of the lensing potential, the lensing calculation is fast to compute. For low ellipticity values, the mass distribution contours have an approximately elliptical shape. However, for high ellipticity values this leads to the mass distribution taking on nonphysical dumbbell-shaped contours, and so modelers often restrict use of the NFWp profile to low values of ellipticity, although the exact definition of “low” may differ depending on the specific lensing application and is ultimately a choice on the part of the modeler. This parameterization is the conventional approach to model NFW profiles in most lensing contexts.

Alternatively, the NFWm parameterization, in which the ellipticity is directly implemented into the mass distribution of the NFW profile, is composed of a combination of cored steep ellipsoid (CSE) profiles with a joint centroid position, each of which has a core radius S​_i and amplitude A_i. Each CSE has the following convergence profile:

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{CSE}}=\frac{A_i{\kappa }_0}{2{(S_i^2+r_{\mathrm{ell}}^2)}^{3/2}},\end{array}$$(4)

where r_ell is used as the radial element. The CSE profile has ellipticity introduced in the mass distribution, and has analytical forms for the lensing potential and its derivatives (Keeton & Kochanek 1998). Oguri (2021) provided an approximation of the NFW profile using 44 CSE components such that the sum emulates the NFW profile to a precision of ∼0.01% in terms of the lensing potential, deflection, and convergence. Calculation of lensing quantities takes approximately ten times as long as the NFWp, but with the advantage that the elliptical mass contour shape is retained for all values of ellipticity.

The input ellipticity for the NFWp profile describes the potential and as such does not equate to the ellipticity of the mass distribution. Therefore, in order to make an apt comparison between the two parameterizations, we modify the input ellipticity to the NFWp profile such that the two parameterizations have the same mass ellipticity at small radii. In Appendix A, we describe analytically the relationship between the axis ratio of the potential and that of the mass for a general profile, and detail our process to match these ellipticities in the specific case of the NFWp profile.

After matching the mass ellipticities, we show in Fig. 1 the shape of the mass distribution for an NFW profile using both the NFWp (red) and NFWm (gray) parameterizations of ellipticity. We show three different values of q_κ (and corresponding q_ψ in the NFWp case) as an input. In the circular case, the two profiles agree. In the elliptical case, the two profiles are not identically shaped but both appear physical. As the ellipticity increases, the NFWp contours become more oval-shaped, and if we were to further increase the ellipticity, they would become a nonphysical dumbbell shape (Kassiola & Kovner 1993). This effect was also discussed in detail by Golse & Kneib (2002), and as such the NFWp profile is typically restricted in use to low values of ellipticity.

Fig. 1. NFW 2D mass isodensity contours of κNFW (top), ellipticity as a function of semimajor axis (middle), and the b4 multipole component as a function of semimajor axis (bottom), in units of NFW scale radius. The NFWp parameterization introduces ellipticity through the lensing potential and is shown in red while the CSE-based NFWm parameterization introduces ellipticity directly in the mass and is shown in gray. The panels from left to right indicate the comparison for three input axis ratios. To match the κNFW ellipticities, we use the procedure in Appendix A, which results in the NFWp potential having an input qψ as indicated above the top panels.

To quantify the azimuthal structure as a function of semimajor axis, we use the photutils software, which uses the method of Jedrzejewski (1987) to fit the elliptical structure of isophotes. The method iteratively fits the position angle and ellipticity of an ellipse. Once this fit is found, departures from ellipticity in terms of Fourier multipole components are calculated as

$$\begin{array}{c}\hfill I=I_0+a_{n}sin(nE)+b_{n}cos(nE),\end{array}$$(5)

where n = 3 or 4, I is the intensity of a given contour, I₀ is the intensity of the fit ellipse contour, and E is the eccentric anomaly. After a fit is achieved, the semimajor axis of the ellipse in question is increased and the process is repeated, such that the result is a description of the ellipticity, position angle, and Fourier amplitudes as a function of semimajor axis. We apply this fitting procedure to the isodensity contours and measure the ellipticity of the mass (expressed as 1 − q) as a function of semimajor axis, which we plot in the middle panels. In terms of the Fourier multipole components, the only component which returns nonzero values is b₄, a parameter of particular interest because nonzero b₄ is often measured in isophotes of real early type galaxies and is known to have an effect on H₀ recovery in lens systems (Van de Vyvere et al. 2022a). A positive value of b₄ indicates a “disky” shape and a negative value indicates a “boxy” shape. In the bottom panels, we plot the fit values of b₄, which indicate that the contours systematically become boxier as a function of semimajor axis.

Setting aside the known q mismatch effect (which we have accounted for) and the “dumbbell-shaped” effect (which does not apply in this ellipticity regime), in this work we identify a third difference in behavior between the two approaches; namely, there is a significant increase in ellipticity with radius. We refer to this effect as an “ellipticity gradient” in this work. Additionally, the mass distributions resulting from NFWp profiles present a nonzero boxiness which increases with radius, even in the case where the mass contours have physically realistic convex shapes. The presence of this type of azimuthal structure even in the low-ellipticity regime (i.e., before the dumbbell shapes arise) is somewhat overlooked in the literature, as in most use cases of the NFWp profile the intention is to mock an elliptical shape, and so an approximately constant elliptical shape is implicitly assumed in the low-ellipticity regime.

2.2. Mass sheet transformation and H₀

It is possible for two different mass distributions to reproduce the same imaging information. The most well studied of such degeneracies is the MST (Falco et al. 1985), where the convergence is rescaled by a constant factor of λ and a uniform sheet of mass¹ is correspondingly added:

$$\begin{array}{c}\hfill {\kappa }_{\lambda }(r)=\lambda \kappa (r)+(1-\lambda ).\end{array}$$(6)

In addition, the unobservable source position is rescaled by the same factor,

$$\begin{array}{c}\hfill {\beta }_{\lambda }=\lambda \beta .\end{array}$$(7)

Under the MST, image positions and relative fluxes are unchanged. Meanwhile, time delays are affected by a factor of λ.

The MST is critical for lensing in a cosmological context, because the true mass distribution of the lens is unknown and lensing cannot distinguish between two mass profiles which are within an MST of one another. Because a given lens could be fit with more than one possible model (e.g., PEMD or composite), it is useful to describe lens profiles in a MST-independent way, meaning to relate only the quantities which lensing directly constrains rather than those which come from model choice. One such MST-independent quantity is the Einstein radius, R_E, for which in this work we adopt the definition of the circular aperture within which the mean integrated surface mass density is equal to the critical density for lensing. The MST-independent nature of this quantity makes it one of the few quantities which lensing directly measures rather than infers from a model, and so it would be exceedingly useful to find more analogous quantities. The dimensionless quantity ξ, expressible using a combination of derivatives of the lensing potential ψ, has been developed for this purpose by Sonnenfeld (2018), Kochanek (2020), Birrer (2021). Together with the Einstein radius, R_E, ξ is an MST-invariant quantity, and so this work uses it as a metric to evaluate the degree to which two mass distributions are within an MST mapping from one to another.

This work primarily uses the Kochanek (2020)ξ₂ definition, which is defined by Taylor expanding deflection for a circular lens for an image near R_E:

$$\begin{array}{c}\hfill {\xi }_2\equiv R_{\mathrm{E}}\frac{\alpha {″}_{\mathrm{E}}}{1-{\kappa }_{\mathrm{E}}},\end{array}$$(8)

where κ_E is the mean convergence and ${\alpha }_{\text{E}}^{″}$ is the second derivative of the deflection, where the “E” subscript refers to the evaluation at r = R_E. The “2” subscript on ξ₂ refers to being derived from the second order term in the expansion, which we share in Appendix B.1. The analogous quantity derived by Sonnenfeld (2018) and shown in terms of the radial stretch factor by Birrer (2021), ξ_rad, is conceptually equivalent with a difference only of a factor of 2,

$$\begin{array}{c}\hfill {\xi }_{\mathrm{rad}}\equiv R_{\mathrm{E}}\frac{{\psi }_{\mathrm{E}}^{″\prime }}{1-\psi {″}_{\mathrm{E}}}=R_{\mathrm{E}}\frac{\alpha {″}_{\mathrm{E}}}{1-{\alpha }_{\mathrm{E}}^{\prime }}=\frac{1}{2}{\xi }_2,\end{array}$$(9)

noting that α = ψ′ and using for a circular lens (Bartelmann 2010),

$$\begin{array}{c}\hfill \kappa =\frac{1}{2}{\mathrm{\nabla }}^2\psi =\frac{1}{2}\left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \psi }{\partial r}\right)\right]=\frac{1}{2}(\psi ″+\frac{\alpha }{r})·\end{array}$$(10)

Evaluating at the Einstein radius and recognizing that α(R_E) = R_E, one can express κ_E in terms of the local first derivative of deflection,

$$\begin{array}{c}\hfill \kappa =\frac{1}{2}({\alpha }_{\mathrm{E}}^{\prime }+1),\end{array}$$(11)

reconciling the two definitions of ξ by inserting κ_E into Eq. (8).

If one mass model is within an MST of another, they will match the same R_E and ξ₂. If the fit model is a power law, the ξ₂ constraint dictates the radial slope of the profile: γ = 2 + ξ₂/2. The convergence of the power law can easily be derived as

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{E},\mathrm{PL}}=\frac{2-{\xi }_2}{4}·\end{array}$$(12)

The accuracy of the Fermat potential,

$$\begin{array}{c}\hfill \tau =\frac{{(\theta -\beta )}^2}{2}-\psi (\theta ),\end{array}$$(13)

can be shown to be proportional to (1 − κ_E) (Kochanek 2002). While Fermat potential is somewhat abstract, the observable quantity is the relative time delay between images Δt, which is given by the Fermat potential, scaled by cosmological distances such that $\mathrm{\Delta }t=D_{\mathrm{\Delta }t}{c}^{-1}\mathrm{\Delta }\tau \propto H_0^{-1}\mathrm{\Delta }\tau$. As such, errors in the Fermat potential are compensated by errors in the recovered value of H₀, such that the observable Δt is correctly matched. Therefore, one can use H₀ as a proxy to represent the accuracy of the Fermat potential, which we elect to use since it is a more intuitive quantity with one specific value for a given model (rather than for each pair of images, which is more cumbersome to express). Through this relation, the expected fractional error on H₀ can be expressed using a comparison between the true convergence (which is generally unknown except in simulated lenses) and that of the power-law model. This error is equivalent to the MST parameter λ if the difference between the profiles corresponds to an MST:

$$\begin{array}{c}\hfill \frac{1-{\kappa }_{\mathrm{E},\mathrm{PL}}}{1-{\kappa }_{\mathrm{E},\mathrm{true}}}=\frac{H_{0,\mathrm{PL}}}{H_{0,\mathrm{true}}}=\lambda .\end{array}$$(14)

Finally, one can estimate the width of the recovered H₀ posterior, ΔH_0, PL, through a propagation of errors on ξ₂:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{H}_{0,\mathrm{PL}}}{H_{0,\mathrm{true}}}=\lambda \frac{\mathrm{\Delta }{H}_{0,\mathrm{PL}}}{H_{0,\mathrm{PL}}}=\lambda \frac{\mathrm{\Delta }(1-\frac{2-{\xi }_2}{4})}{1-\frac{2-{\xi }_2}{4}}=\lambda \frac{\mathrm{\Delta }{\xi }_2}{2+{\xi }_2}\simeq \lambda \frac{\mathrm{\Delta }{\xi }_2}{2},\end{array}$$(15)

where we have used ξ₂ ∼ 0 ≪ 2, noting that galaxy-scale lens systems have mass profile slopes near isothermal, that is, γ = 2 (e.g., van de Ven et al. 2010; Auger et al. 2010; Shajib et al. 2021). For example, to reach a 1% accuracy on H₀, the target precision for a population of systems, ξ₂ must be known to within an absolute error ∼0.02 (assuming λ ∼ 1). Since ξ₂ ∼ 0, this description using absolute error in ξ₂ is more robust than a description using relative error (such as that of Kochanek 2021).

In practice, the utility of ξ₂ has been considered for two main purposes. First, because it is considered to be model-independent, a measurement of ξ₂ from one model could be used to help guide another model in the optimization process. One example of this idea was demonstrated by Shajib et al. (2021), who included lensing information in a kinematics analysis by folding the posterior distributions for R_E and ξ₂ into the Bayesian framework, rather than jointly modeling the lensing and kinematic parameters together, thereby under the implicit assumption that the stars and DM in the composite lens share the single ellipticity of the lens model. Secondly, ξ₂ has been considered for its application on the simulation side for systematics testing of time-delay cosmography. In many cases one wishes to consider the case where the true lens is a more complex mass distribution than the model, and evaluate the errors introduced by an overly simplistic lens model (e.g., Cao et al. 2022; Van de Vyvere et al. 2022a,b). Because the process of creating and fitting mock images is slow, it would be expedient to be able to estimate what the recovered H₀ would be from a given model directly from the mass distribution (see Xu et al. 2016; Tagore et al. 2018; Gomer & Williams 2020). We have considered using ξ₂ for this purpose (predicting the would-be PEMD fit and using Eq. (14)), which would enable one to create much larger samples of lenses for systematics checks. Our attempt to confirm this method using the NFWp mocks in TDCVIII did not match our expectation to the desired precision (approximately 3% discrepancy on H₀, see Sect. 3). Several explanations of this mismatch are possible, such as a possible bias in the estimate of κ_E, or higher order terms being required in the Taylor expansion used to define ξ₂ (see Appendix B). In the next section, we show that the main driver of the mismatch unveiled in TDCVIII is directly attributable to the implicit azimuthal structure introduced through the NFWp profile.

3. Comparing fits to both NFW parameterizations

The experiment in this work is designed to probe the effect of the choice of NFW ellipticity parameterization in the context of galaxy-scale strong lensing. The NFW component represents the DM of a galaxy, but since no galaxy is purely DM, the impact of the NFW parameterization should be studied in the framework of a composite model with both baryon and DM mass components. TDCVIII created a population of mock systems analogous to the TDCOSMO lens population using a Chameleon profile (Dutton & Treu 2014) for the baryon component and a NFWp profile for the DM component, constructed to have stellar and DM mass distributions typical of real lenses. The Chameleon profile has a convergence given as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\kappa }_{\mathrm{Cham}}=& \frac{A_0}{(1+q_{\kappa })}[\frac{1}{\sqrt{r_{\mathrm{ell}}^2+4w_{\mathrm{c}}^2/{(1+q_{\kappa })}^2}}\hfill \\ \hfill & -\frac{1}{\sqrt{r_{\mathrm{ell}}^2+4w_{\mathrm{t}}^2/{(1+q_{\kappa })}^2}}],\hfill \end{array}\end{array}$$(16)

with parameters w_c and w_t with w_t > w_c where A₀ sets the mass scale at zero radius. and is tailored to closely mimic a Sérsic profile with a given R_Sersic and n_Sersic, but with the benefit that its lensing potential is analytically expressible. We create 20 mocks following the same strategy as TDCVIII, but we replace the NFW component first with an NFWp component and later with an NFWm component. We note that the ellipticities of the convergence of the Chameleon and NFW components were not matched in TDCVIII, so we recreate the NFWp mocks using the ellipticity matching procedure discussed in Appendix A. This population of mock lenses and corresponding fits provides an excellent laboratory to explore the role of the implicit azimuthal structure embedded within the NFW parameterization.

With fits to both populations, we can compare the resulting values of H₀ and ξ₂ to determine if the parameterization plays a relevant role for time delay cosmography. In both cases, the composite profile will not recover the fiducial value of H₀ due to the MST, but if the MST is the only effect at play, the two should recover the same values of H₀ and ξ₂, biased from the fiducial H₀ by the amount predicted by Eq. (14). We show that the implicit azimuthal structure in the potential-based parameterization causes a deviation from this prediction, inconsistent with an MST.

An alternative design for a similar experiment is to use a Chameleon+NFW profile as the model to fit a lens, and to test if the NFW parameterization affects the recovery of H₀ and other parameters. Under an ideal setup, such an experiment would be more akin to how time-delay cosmography modeling is implemented and able to more directly attack the question of whether or not the NFW parameterization introduces a bias in these conditions. However, such a setup introduces complexities in that the input azimuthal structure is more difficult to describe and control. While it is possible to create an input population with a reasonable approximation of the azimuthal structure of real systems (see e.g., Van de Vyvere et al. 2022a,b), we instead elect for this more controlled experiment where lensing degeneracies can be more directly probed because the azimuthal structure can be traced exclusively to the NFW parameterization, with the acknowledgement that this setup is not directly equivalent to the common practice of using an NFW profile to fit a mock.

3.1. Experiment specifics

We use the same input systems as TDCVIII, which were constructed to match the observed lens population. A large number of two-component profiles were synthesized, and then a subset was selected which matched several observable quantities such as the ellipticities, Einstein radii, and effective half-light radii of real systems. By probing this population of parameters, we ensure that the results of this experiment hold across a range of realistic lenses. We share the parameters for our mock lenses in Table 1. For more details, see Gomer et al. (2022). The result is a set of input two-component profiles with a realistic distribution of parameters, including ellipticities ranging roughly from input q = 0.6 to q = 0.8. TDCVIII simply used this input q for both components, resulting in a mismatch between the q_κ of the Chameleon stellar component and the q_ψ of the NFWp DM component. We instead rescale the q_ψ of the NFWp component according to Eq. (A.7) to match the q_κ between the two components. We synthesize 20 mock lens images using both the NFWp and NFWm parameterizations.

We show three example composite profiles drawn from our population in Fig. 2 for both the NFWp and NFWm implementations. We note that the field of view differs from Fig. 1 to better view the structure at the Einstein radius, which is significantly interior to the NFW scale radius. Inside the Einstein radius, the Chameleon profile dominates, resulting in a nearly perfect elliptical shape. Only in the outer regions does the ellipticity gradient and nonzero boxiness resulting from the NFWp profile become apparent. Radially, the two NFW implementations have essentially identical structure.

Fig. 2. Three composite mass profiles drawn from our sample, with the Einstein radius (green, dotted) and NFW scale radius (blue, dotted) indicated. The implementation of these composite profiles is shown for NFWp (red dashed) and NFWm (gray solid) DM profiles. Top row shows isodensity contours; second row shows ellipticity as a function of semimajor axis; third row shows b4 as a function of semimajor axis. The bottom panels show the radial profile over a log scale, defined as the value of κ in a circular annulus at a given radius, rather than with respect to semimajor axis, expressed relative to the NFWm radial profile. The composite profiles are plotted as well as both components individually. Units of distance are now in arcseconds.

Similar to TDCVIII, we use lenstronomy² (Birrer & Amara 2018; Birrer et al. 2021) to create and fit our mock lenses. Our NFWp mocks are the same as the set used in TDCVIII (specifically the TDCOSMO-like set in that work), except that we rescale the ellipticity so that the mass components match. Our NFWm mocks simply substitute the NFW components for the NFWm parameterization. One last small change we make between our two parameterizations compared is to set the source position relative to the caustic in order to maintain the same image configuration; since the caustic changes slightly between the two parameterizations, this slightly moves the source position. This effect is further quantified in Sect. 5.1.

We show a comparison between the two mocks for an example image in Fig. 3. The images are quite similar, with the main differences coming from slightly different point source magnifications with no discernible differences in the arcs. The input Einstein radii and ξ₂ values closely match between both parameterizations: the Einstein radii match to within our numerical uncertainty while the input ξ₂ values match to better than 0.01 in all cases, with a mean difference of 0.001.

Fig. 3. Example image comparison between NFWp and NFWm ellipticity parameterizations (top), with resulting residuals from the PEMD+shear fits (bottom). Caustics for both parameterizations are shown in the bottom right with the source position rescaled as described in the text and indicated as a cross. This example image comes from a system with moderately large ellipticity with an input axis ratio of qκ = 0.64.

3.2. Fit results

Like TDCVIII, we fit our mock population with a PEMD+shear model. All systems are fit well, with typical residuals comparable to Fig. 3 (bottom). We find that the azimuthal structure attributable to the NFWp parameterization can be absorbed by the lens model, similar to Van de Vyvere et al. (2022a,b), who found that nonzero multipole components and ellipticity gradients can often be absorbed by the lens model so long as the deviations from a constant ellipse shape are not too extreme.

If the MST is the only effect at play, the recovered values of H₀ should match the predictions according to Eq. (14) based on the ξ₂ of the input mass distribution. We plot the fit values of H₀ for both the NFWm population and the NFWp population in Fig. 4 alongside the predicted H₀. Error bars are estimated via a Markov chain Monte Carlo (MCMC) estimation using emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013). The fits to the NFWp systems are systematically biased relative to the expectation by approximately 2.5%, while the fits to the NFWm systems lie on the expectation line. This result indicates that the Fermat potential recovered by the model is inaccurate by approximately 2.5%. This discrepancy appears to be directly caused by the systematic azimuthal structure introduced by the NFWp profile. The effect was originally reported in TDCVIII, although since that work did not match the q_κ values between the baryon and DM components and therefore used a more elliptical NFWp component, the magnitude of the effect was slightly larger than seen in this work, quoting a discrepancy of 3%.

Fig. 4. Expected H0 in km s−1 Mpc−1 calculated using the ξ2 calculated from the input mass profiles for 20 mock lenses compared to the value of H0 recovered from a PEMD fit. The dotted line indicates a 1:1 correspondence.

3.3. Numerical checks

For completeness, we consider the hypothetical possibility that κ_E could be numerically biased when comparing one parameterization to another, resulting in the observed H₀ discrepancy when Eq. (14) is applied. As such, we also plot directly the ξ₂ values of the PEMD fits compared to those of the input mass distributions in Fig. 5. Again the fits using NFWm parameterization fall tightly on the expected 1:1 line, while those of the NFWp parameterization are systematically biased by approximately 0.05, consistent with Eq. (15). This result confirms that the discrepancy arises because of implicit azimuthal structure in the NFWp parameterization, rather than numerical effects due to κ_E estimation.

Fig. 5. Input ξ2 for the 20 mock lens mass profiles compared to the recovered ξ2 from the PEMD fits, with a 1:1 relation indicated as the dotted line.

We have tested several other possibilities to explain this mismatch. Some sources of error we quantified are detailed in Appendix B. Namely, we quantify errors due to the truncation of the Taylor expansion used to define ξ₂, the effect of numerical errors in ξ₂ calculation, and the effect of azimuthal averaging. These effects cannot introduce systematic errors above the percent level, and so we ultimately conclude that the mismatch has been caused by the ellipticity parameterization within the NFW profile.

4. Discussion

We have found that the NFWp parameterization of the NFW profile, when fit with a PEMD+shear model, results in a mismatch between the input ξ₂ and that of the fit, ultimately resulting in an H₀ which is systematically underpredicted. Theory predicts that the fit value of ξ₂ will match the input, since it is the only quantity (along with R_E) which lensing is able to constrain in an MST-independent manner. This leads one to wonder what exactly this mismatch means.

Our interpretation of this result is that the azimuthal structure in the input was not adequately accounted for by the fit model, analogous to the example in Kochanek (2021). The PEMD model, even with external shear, has no ellipticity gradients or nonzero b₄ which we now know to be ubiquitous in the NFWp parameterization, shown in Fig. 1. Lacking the capacity to include this azimuthal structure, the PEMD recovered a biased value of ξ₂ compared to the value which would have been recovered if the fit had sufficient azimuthal freedom. When the mocks were replaced with those with constant ellipticity, the expected value of ξ₂ was recovered. The takeaway is that the expected ξ₂ of a fit model and the actual ξ₂ of a general mass distribution will not be equivalent unless the azimuthal structure of the model is able to match the truth.

In other words, this mismatch occurs because the mapping from the input to the PEMD fit cannot be represented by an MST alone, as ξ₂ is invariant under the MST by construction. Nonetheless, the resulting image fits are good with no residuals, meaning the input mass model with its additional azimuthal structure is degenerate with a PEMD+shear model to within the image noise. We state this explicitly to highlight that the degeneracies at play in this experiment go beyond the MST.

The Source Position Transformation (SPT; Schneider & Sluse 2014; Wagner 2018), which generalizes the MST to more complex transformations of the radial distribution of the lens, but also yields azimuthal change of the mass (Unruh et al. 2017), may appear to be the degeneracy at work in the present experiment. One can show from Eqs. (10) and (11) that ξ₂ can be expressed as in terms of ratios of derivatives of κ:

$$\begin{array}{c}\hfill {\xi }_2=R_{\mathrm{E}}\frac{{\kappa }_{\mathrm{E}}^{\prime }}{{\kappa }_{\mathrm{E}}}+1.\end{array}$$(17)

As Unruh et al. (2017) show, such ratios of κ derivatives are invariant under the MST, but can change under the SPT. Through ξ₂, we show that this ratio is unchanged for the fits to the NFWm parameterization, but changed for those of the NFWp parameterization, meaning that the NFWp profiles are not an MST away from the fit, but the transformation is consistent with an SPT. The SPT is an approximate global degeneracy, and so we check the values of the relative time delays between the input and the PEMD fit. Under an MST, all three time delays are rescaled by the same constant λ, but under an SPT, the three relative delays are not scaled by the same constant value (Wertz et al. 2018). We find exactly this result, that each individual time delay of the fit is off by the same constant factor for the NFWm mocks but each are off by different factors for each delay for the NFWp mocks, supporting the conclusion that the NFWm mocks are an MST away from a PEMD, while the NFWp mocks are an SPT away from a PEMD. However, the SPT is not a complete description of this degeneracy, because the SPT also transforms the source shape, while all mocks in this work are created and fit using a circular source. Therefore we suspect that there is a related degeneracy at work, likely in the form of a “shape degeneracy” as discussed by Saha & Williams (2006), who found similar results using pixelated mass profiles with ellipticity gradients, albeit limited to point sources. Like the SPT, this degeneracy does not uniformly affect time delays, consistent with our findings. Although it is difficult to be more quantitative about their exact form and contribution, it is clear that higher-order lensing degeneracies beyond the MST are certainly at play in this work.

It may look surprising that while most of the azimuthal changes in the input profile appear outside the Einstein radius, (i.e., where the NFW profile starts to dominate), the impact on ξ₂ and H₀ remains noticeable. This result indirectly shows that morphological assumptions on the density profile beyond the Einstein radius can have a substantial impact on the lensed images.

We find that the discrepancy between the recovered value of ξ₂ for the Chameleon+NFWp and Chameleon+NFWm profiles moderately correlates with the input axis ratio, with Pearson correlation R = −0.54. The discrepancy also correlates with the recovered value of external shear, with R = 0.70. We interpret this to mean that the deviation from an MST worsens with ellipticity, and that external shear can help to absorb this more complex degeneracy. These results add to the body of evidence that external shear can sometimes reflect absorbed degeneracies rather than a physical quantity (Etherington et al. 2023).

One open question for lensing theory is how one should describe the effects of these shape degeneracies in order to include them in lens models. It may be possible to construct a quantity which is independent of these degeneracies analogous to the way ξ₂ is constructed to be independent of the MST, but we were unable to derive such a quantity in the confines of this work. The problem is that outside the circular limit, Eq. (10) requires an azimuthal ∂²ψ/∂ϕ² term, entangling azimuthal dependence in any global expression of ξ. One should instead use a local quantity such as the stretch differentials discussed by Birrer (2021). The ultimate goal would be to have a description of how this local quantity, integrated over the imaging information, changes with azimuthal structure. General azimuthal structure is difficult to parameterize, but to test the concept, we created an experiment where we introduced a change in the input q_κ, but left it unmodeled in the fit, then calculated the expected change in the integrated stretch factor, and found that it corresponds to the amount by which the PEMD ξ recovery is biased. This is still an open field of research, but we believe a complete description of lensing degeneracies requires an exploration along these lines.

This result puts some limitations on the practical applications of ξ₂. Firstly, the practice of using the recovered value of ξ₂ from a simple model to place a constraint on a more sophisticated model is not precise to more than a few percent in the general case. If the true mass distribution has azimuthal structure, the simple fit will recover a biased value of ξ₂ which then places an inaccurate constraint on the more sophisticated model. This practice can only work if the true mass distribution lacks azimuthal structure, in which case the simple fit will recover an unbiased ξ₂. Secondly, the use of ξ₂ on the simulation side to estimate the systematic effects of simplistic lens models comes with limitations as well, because the value of ξ₂ from an input mass distribution describes what would be recovered by a model that shares the same azimuthal complexity as the data. To perform systematic tests for cases where the azimuthal structure of the mock and of the model differs, the need to create and fit mock lens systems cannot be circumvented.

5. Connections with other works

In this section we discuss the implications of these results on several related fields of study: namely, measurements of ellipticity and lensing cross section (Sect. 5.1), flux ratio anomalies (Sect. 5.2), and H₀ determination (Sect. 5.3).

5.1. Ellipticity and lensing cross section

Two-component lens models are often used to provide observational constraints on the DM components of lens systems. In a science case where one wishes to know the ellipticities of DM mass distributions, modelers are already wary not to conflate the ellipticity ascribed to the lensing potential with that of the mass (Kassiola & Kovner 1993; Barkana 1998; Golse & Kneib 2002). However, when two-component mass models are used to fit lens systems, the NFWp axis ratio q_ψ, is sometimes quoted alongside the baryon component q_κ with little distinction made between the quantities, opening the door for confusion if one were to take the quoted q at face value. Perhaps deemed irrelevant to a given science case, the distinction between the NFWm and NFWp profile has been largely neglected, and hence the azimuthal structure introduced in the mass distribution even in the physical low-ellipticity regime has been somewhat overlooked.

When the NFWp parameterization is used, this slightly changes the caustic size, which can play a role in selection effects affecting lensing studies (Baldwin & Schechter 2021). Studies concerning lensing cross section should be wary about the NFWp parameterization artificially increasing the cross section. For the example case in Fig. 3, we find the caustic size changes by 3% (6% by area) compared to the NFWm case. We note that this example system has an input q_κ of 0.64, which is quite elliptical for our sample. As such this change in caustic size is likely somewhere between the average case and the extreme case. Finally, we also note that quantities defined using 1D profiles, such as the DM fraction, are unchanged by the choice of NFW parametrization in the potential or in the mass.

Through the rescaling of ellipticity implemented in this work in Appendix A, it is now relatively simple to convert from a quoted NFWp axis ratio result to what the corresponding mass axis ratio is in the center of the mass distribution. Rather than repeating previous work, this conversion may suffice depending on one’s intended scientific application.

5.2. Flux ratio anomalies

A second important consideration is the consequences of this work on lensing studies involving calculations of image flux ratios. For cusp-configuration systems, in which the three coalescing images of the cusp each have a signed magnification μ, one can define (e.g. Keeton et al. 2003)

$$\begin{array}{c}\hfill R_{\mathrm{cusp}}=\frac{{\mu }_1+{\mu }_2+{\mu }_3}{|{\mu }_1|+|{\mu }_2|+|{\mu }_3|},\end{array}$$(18)

which approaches zero as the source approaches the caustic. Similarly for the two coalescing images in fold systems,

$$\begin{array}{c}\hfill R_{\mathrm{fold}}=\frac{{\mu }_{\min }+{\mu }_{\mathrm{saddle}}}{|{\mu }_{\min }|+|{\mu }_{\mathrm{saddle}}|},\end{array}$$(19)

which also approaches zero as the source approaches the caustic. Deviations from these theoretical cusp and fold relations can be used to diagnose substructure within a lens mass, which have been interpreted as DM subhalos (e.g., McKean et al. 2007; MacLeod et al. 2013; Nierenberg et al. 2014) or as evidence of galaxy group effects or more complex macro-scale mass distributions (Xu et al. 2015).

Considerate of this, we evaluate if the distribution of R_cusp and R_fold would be changed by the NFW parameterization, and so we check this distribution for one of our systems by generating 500 sources for the caustics in Fig. 3. Lensing these sources, we designate systems as folds or cusps according to the criteria of Keeton et al. (2005), based on how many images lie within ∼1 R_E of one another. We plot the resulting R_cusp and R_fold distributions in Fig. 6. We find that the distribution of these flux ratios does not significantly change between the two NFW parameterizations. Comparing the distributions using a Kolmogorov–Smirnov test, we find p-values of 0.16 for R_cusp and 0.72 for R_fold, far from the threshold for similarity typically set at p < 0.05. As such, we conclude that the NFW parameterization does not directly impact such studies.

Fig. 6. Distributions for Rcusp (left) and Rfold (right) for many realized source positions of an example mock.

5.3. H₀ determination

Because we used NFW profiles in the input mass, rather than fitting a system using an NFW profile, the experiment in this work cannot prove whether or not this effect biases H₀ in real systems, but only draws attention to the fact that the NFWp and NFWm make different assumptions. Nonetheless, we believe that the lessons learned in this work regarding azimuthal structures may have some implications for this common use of NFW profiles.

When real systems are modeled, the true mass distribution is unknown. Because the NFWp and NFWm parameterizations give different azimuthal prescriptions, a modeler must decide which description of ellipticity they wish to assume. After all, ellipticity gradients and nonzero b₄ exist within real galaxies, and so a modeler may want to select a model which includes such azimuthal structure. However, we note that NFWp azimuthal structure is not the same as what is found in real systems, as it systematically increases ellipticity and boxiness with radius as a nonphysical consequence of the lensing potential, rather than representing the physical tendency of galaxies to have ellipticity which can increase or decrease with radius and structures which can be disky as well as boxy. Van de Vyvere et al. (2022a,b) studied the effect of such azimuthal structures in lens models, broadly concluding that individual lenses can result in biased parameter recoveries while the population averages out to an unbiased determination of H₀. Modelers should therefore be wary about the azimuthal structure introduced by the NFWp model; since this effect always applies in the same direction, it will not average out over the population of lenses. If the intention is that a mass profile has a constant elliptical shape, using the NFWp parameterization is not consistent with this assumption and would introduce a systematic effect in the modeling.

However, this is not the whole story when it comes to modeling the azimuthal structure of real systems. When systems are modeled with a composite profile, they are done so using more complex models than we have used here. Centroid positions, position angles, and ellipticities may be allowed to be offset between the two components (e.g., Rusu et al. 2020; Shajib et al. 2022). It is entirely possible for choices in this parameter space to compensate for additional azimuthal structure in the NFWp component. As a simple example, using a circular NFW component would result in an ellipticity gradient going from elliptical in the center to circular in the outer regions, in the opposite direction as the NFWp ellipticity gradient, perhaps negating or reversing its effect.

In addition, nearby perturber galaxies are included in the model (Wong 2018; Birrer et al. 2019; Shajib et al. 2020) and pixelated corrections to the lensing potential can be implemented (Suyu et al. 2010). These numerous considerations add considerable azimuthal freedom to the model. This freedom may be sufficient to describe the true mass model, in which case H₀ would not be biased beyond the MST. This possibility has not been directly tested, but the fact that TDCOSMO recovers consistent H₀ values between their power law models and composite models supports this hypothesis (see Fig. 6 of Millon et al. 2020). We do however note that in such a case the particular values of the individual components may not correspond to the true mass distributions, instead seeking a compromise which compensates for the NFWp gradients.

Furthermore, the present work does not include stellar kinematic constraints, which play a vital role in breaking lensing degeneracies (Birrer et al. 2020; Yıldırım et al. 2023; Shajib et al. 2023). Finally, we note any implications of the systematic effects discussed in this work would only apply to NFW profiles and therefore only to the composite fits of TDCOSMO, having no bearing on the power-law fits also adopted by TDCOSMO.

6. Conclusion

The NFW profile is a key ingredient of realistic models of galaxies, but it cannot be described analytically in the elliptical case. For lensing applications, the ellipticity can either be added in the potential (NFWp parameterization) or in the mass via an approximated profile (NFWm parameterization). When ellipticity is introduced in the potential, it introduces an azimuthal structure in the form of ellipticity gradients and nonzero boxiness in the mass distribution. We created two populations of composite mocks using each parameterization and fit them with a PEMD+shear model. The mocks created using the potential-based parameterization result in fits that recover biased values of H₀ relative to those from the mocks without this introduced azimuthal structure. This result has several consequences:

• The use of a potential-based parameterization of ellipticity introduces a ubiquitous azimuthal structure in κ in the form of ellipticity gradients and nonzero b₄, even for low values of ellipticity when the distribution is not yet dumbbell-shaped. We note that the presence of artificial variations of ellipticities can be mistakenly absorbed by shear (Van de Vyvere et al. 2022b; Etherington et al. 2023). We advise lens modelers who wish to assume an azimuthal shape with constant ellipticity to use the NFWm parameterization (which may be based on CSEs as we have used in this work or another formulation that keeps the ellipticity constant with radius) in order to be consistent with this assumption. Azimuthal structures can still be implemented through multiple mass components with differing ellipticities or position angles, but it would be done so explicitly rather than unintentionally.

• Our fits to mocks with both NFW parameterizations resulted in values of H₀, which were discrepant with one another by 2.5% (systematic), indicating an inadequacy of the model to capture the true Fermat potential at this level. However, we cannot claim that the practice of fitting mocks with NFW models necessarily introduces a bias on H₀ at the same level due to the additional azimuthal freedom of TDCOSMO-like models, which may compensate for this effect.

• The MST-independent quantity ξ₂ is an accurate predictor of the recovered mass model in the case where the input and the fit have the same azimuthal prescription. However, when the azimuthal structure in the input mock is not captured by the PEMD model, the value of ξ₂ is biased, indicating that the mapping is not a simple MST, and may in fact be a more general SPT or even a shape degeneracy. Various tests and subtleties in the possible ways to calculate ξ₂ are described in Appendix B.

• The bias introduced by the NFWp parameterization is mostly caused by a deviation from elliptical isodensity contours that take place beyond the Einstein radius. This indirectly shows that morphological assumptions motivated solely by the shape of the lensing galaxy interior to the Einstein radius and/or a generic absence of ellipticity of the NFW component may introduce a bias as large as several percent in some lensing-inferred quantities such as H₀, q, or the external shear magnitude.

• As we have shown that ξ₂ differs between two models with a different azimuthal structure, it can only be used to convert one radial profile to another by keeping the same model assumption on the azimuthal structure. The use of ξ₂ derived from PEMD modeling has been considered as a proxy for constraining other models, such as a composite model, without directly optimizing the model on the lensed images. As an early example of such an application, Shajib et al. (2021) used the ξ₂ from a PEMD lens model to constrain a composite model’s radial profile in dynamical modeling, although in this example the model was spherical and as such the treatment of the azimuthal structure is irrelevant. More generally, the accuracy of this procedure is limited to the amount by which the original power-law lens model is able to capture the azimuthal structure: in our case having an inaccuracy of approximately 0.05 on ξ₂. If a high accuracy is required, we caution against the use of ξ₂ as a constraining diagnostic when comparing models with a differing azimuthal structure.

• For a test case with significant ellipticity (input q = 0.64), the introduced azimuthal structure also changes the cross section for quad lenses by ∼6%. This is particularly relevant for understanding the selection function of lensed systems in existing and upcoming large surveys such as Euclid (e.g., Sonnenfeld et al. 2023).

• Flux ratios remain broadly unaffected by the choice of NFW parameterization of the macro model of the lens.

While the true mass distributions of lenses are not exactly known, it is important to quantify the effects of implicit assumptions inherent in the choice of lens model. The exact prescription of the NFW ellipticity is one such assumption that we show can have an effect if a high accuracy is required.

Acknowledgments

In addition to those mentioned in the text, this work uses the following Python packages: Python (Oliphant 2007; Millman & Aivazis 2011), Astropy (Astropy Collaboration 2013, 2018), Numpy (van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Pandas (McKinney 2010; The pandas development team 2020), and Seaborn (Waskom 2021). We thank the referee, whose insights regarding the ellipticity of the NFWp mass distribution are reflected in Appendix A, improving the experiment design of this work. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 787886). Support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51492 awarded to AJS by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

References

Appendix A: Ellipticity matching to an elliptical potential

Given that a lensing potential results in elliptical contours, this gives rise to a non-elliptical shape of the convergence profile, including ellipticity gradients. Here we analytically show the origin of these features, as well as analytically calculate the mass axis ratio for an NFW potential in the limit of small radius. Starting with a general potential, consider any elliptical potential of the form

$$\begin{array}{c}\hfill \psi (x,y)=f(a),\end{array}$$(A.1)

where

$$\begin{array}{c}\hfill a=\sqrt{x^2+\frac{y^2}{q_{\psi }^2}}\end{array}$$(A.2)

is the semimajor axis which captures all the dependence on the x and y coordinates. We note that the elliptical radii discussed in this work, $r_{\mathrm{ell}}=a\sqrt{q_{\psi }}$ and $r_{ϵ}=r_{\mathrm{ell}}\sqrt{2q_{\psi }/(1+q_{\psi }^2)}$, both take on the form g(q_ψ)*a and as such capture all spatial dependence in terms of a, resulting in an elliptical shape. The convergence then takes the form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \kappa (x,y)& =\frac{1}{2}{\mathrm{\nabla }}^2\psi (x,y)\hfill \\ \hfill & =\frac{1}{2a}(1+\frac{1}{q_{\psi }^2})f^{\prime }(a)+\frac{1}{2a^2}(x^2+\frac{y^2}{q_{\psi }^4})(f″(a)-\frac{1}{a}{f}^{\prime }(a))\hfill \end{array}\end{array}$$(A.3)

A function of this general form

$$\begin{array}{c}\hfill \kappa (x,y)=A(a)+(x^2+\frac{y^2}{q_{\psi }^4})B(a)\end{array}$$(A.4)

is no longer a pure function of a, now containing additional x and y dependence, ergo no longer having purely elliptical contours. To calculate the axis ratio of this convergence, one can set x or y equal to zero to calculate the value of κ along the major axis or minor axis:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \kappa (x,y=0)& =\frac{1}{2}[\frac{f^{\prime }(a|a=x)}{x{q}_{\psi }^2}+f″(a|a=x)]\hfill \\ \hfill \kappa (x=0,y)& =\frac{1}{2}[\frac{q_{\psi }}{y}{f}^{\prime }\left(a|a=\frac{y}{q_{\psi }}\right)+\frac{1}{q_{\psi }^2}f″\left(a|a=\frac{y}{q_{\psi }}\right)].\hfill \end{array}\end{array}$$(A.5)

By setting κ(x_c, y = 0) = C and κ(x = 0, y_c) = C, one can solve for the x_c and y_c values corresponding to a particular isocontour with value C. This solution cannot be expressed for a general potential, but one can show that the NFW potential (Eq. (3)) can be expressed in the limit of small a as

$$\begin{array}{c}\hfill {\psi }_{\mathrm{NFW}}\simeq -a^{2}log\frac{a}{2},\end{array}$$(A.6)

resulting in invertible expressions for Eqs. A.5. Solving for x_c and y_c and taking the ratio of y_c/x_c gives the convergence axis ratio of a given contour, which in this limit is independent of the isocontour C,

$$\begin{array}{c}\hfill q_{\kappa ,a\to 0}=\frac{y_{\mathrm{c}}}{x_{\mathrm{c}}}=q_{\psi }exp\left(\frac{q_{\psi }^2-1}{q_{\psi }^2+1}\right).\end{array}$$(A.7)

We use this expression to calculate the input q_ψ for our NFWp profiles, guaranteeing that they have the same q_κ as the Chameleon profiles at innermost radii.

To more completely illustrate the relationship between q_ψ and q_κ, for an NFW profile, we recreate Fig. 2 of Golse & Kneib (2002), which gives the mass ellipticity as a function of the input potential ellipticity, where ϵ = (1 − q²)/(1 + q²). The mass axis ratio is calculated numerically by determining the point along the y axis which has the same convergence as a point on the x axis and taking q_κ = y_c/x_c. We plot this in the left panel of Fig. A.1. This figure shows that the relation between the two ellipticities is never equality, and the difference between them changes as a function of r: an ellipticity gradient. We also plot the same relation in terms of axis ratios q_ψ and q_κ (right panel). Plotted this way, one can see clearly that q_κ goes as Eq. (A.7) in the limit of small r, which is well approximated as $q_{\psi }^2$ for q_ψ > 0.6.

Fig. A.1. Relationship between the mass ellipticity and potential ellipticity for the NFWp profile. Left: ellipticity for different r values. Right: same relation in terms of axis ratios.

Appendix B: Systematics checks

B.1. Taylor expansion

We discuss here the Taylor expansion of Kochanek (2020) used to define ξ₂, and its accuracy as a function of the limiting order of the expansion. Expanding deflection for a circular lens for an image near R_E:

$$\begin{array}{c}\hfill \alpha (r)=R_{\mathrm{E}}+{\alpha }_{\mathrm{E}}^{\prime }(r-R_{\mathrm{E}})+\frac{1}{2}\alpha {″}_{\mathrm{E}}{(r-R_{\mathrm{E}})}^2+\frac{1}{6}{\alpha }_{\mathrm{E}}^{″\prime }{(r-R_{\mathrm{E}})}^3+\dots \end{array}$$(B.1)

From the lens equation (β = θ − α), where θ = r for a given image, β is given to second order in (r − R_E) as

$$\begin{array}{c}\hfill \beta (r)\simeq -2({\kappa }_{\mathrm{E}}-1)(r-R_{\mathrm{E}})-\frac{1}{2}\alpha {″}_{\mathrm{E}}{(r-R_{\mathrm{E}})}^2.\end{array}$$(B.2)

where we have used Eq. 11. Dividing by (1 − κ_E),

$$\begin{array}{c}\hfill \widehat{\beta }(r)\equiv \frac{\beta (r)}{1-{\kappa }_{\mathrm{E}}}\simeq 2(r-R_{\mathrm{E}})-\frac{1}{2}\frac{\alpha {″}_{\mathrm{E}}}{1-{\kappa }_{\mathrm{E}}}{(r-R_{\mathrm{E}})}^2.\end{array}$$(B.3)

One can see by substituting Eqs. 6 and 7 that $\widehat{\beta }(r)$ is invariant under the MST. As such, the right-hand side of the equation is also MST-invariant. Furthermore, since R_E is MST-invariant, the second-order term is MST-invariant as well. From here, ξ₂ is defined according to Eq. 8 by making the second term unitless via a factor of R_E.

Let us consider the error associated with the truncation of the Taylor expansion. Using any finite order for a Taylor expansion will introduce some discrepancy from the truth $\mathrm{\Delta }\widehat{\beta }$:

$$\begin{array}{c}\hfill \widehat{\beta }={\widehat{\beta }}_{\mathrm{Taylor}}+\mathrm{\Delta }\widehat{\beta },\end{array}$$(B.4)

where, for two terms, ${\widehat{\beta }}_{\mathrm{Taylor}}$ is defined by Eq. B.3. Since R_E is measured very precisely, consider that the error due to the truncation of the Taylor expansion will be interpreted as error in ξ₂. Rearranging this equation for a second order ${\widehat{\beta }}_{\mathrm{Taylor}}$,

$$\begin{array}{c}\hfill \widehat{\beta }=2(r-R_{\mathrm{E}})-\frac{1}{2R_{\mathrm{E}}}{(r-R_{\mathrm{E}})}^2[{\xi }_2-\frac{2R_{\mathrm{E}}}{{(r-R_{\mathrm{E}})}^2}\mathrm{\Delta }\widehat{\beta }].\end{array}$$(B.5)

We define the error in ξ₂ as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\xi }_2=-\frac{2R_{\mathrm{E}}}{{(r-R_{\mathrm{E}})}^2}\mathrm{\Delta }\widehat{\beta }.\end{array}$$(B.6)

We construct this quantity this way because we are curious how much the error associated with the truncation of the Taylor expansion can result directly in an error in ξ₂. This implicitly assumes that there is no error in R_E such that all of the truncation error is applied to ξ₂.

Though this quantity is defined using two terms in the Taylor expansion, it can be useful to calculate it using n terms in the expansion for ${\widehat{\beta }}_{\mathrm{trunc}}$ to estimate the increased accuracy of a higher order expansion. With three terms, for example, Δξ₂ under this construction assumes all error associated with the truncation still applies to ξ₂, which is strictly speaking inaccurate because ξ₃ should have error associated with it. Therefore, this representation serves as a conservative estimation of the maximum possible error in ξ₂, signifying by how much ξ₂ would need to change to compensate for the error in the expansion.

We plot this quantity in Fig. B.1 using a power law for which the true $\widehat{\beta }$ can be calculated analytically, using several different slopes and several values of n. We evaluate this error both in the case where r = 0.7R_E (interior to R_E) and when r = 1.3R_E (exterior to R_E), such that |r − R_E|/R_E = 0.3 in both cases. We find that for slopes between γ = 1.7 and γ = 2.3, the error is centered on zero with lessening scatter as the number of terms increases. The error is approximately 5 times larger for images interior to R_E than those exterior to R_E, but for n = 2 the error is always less than 0.04. This makes it unable to explain the systematic difference of ≃0.06 discussed in this paper. Furthermore, since the error is centered on zero, a population with an average slope of 2 will not have any systematic bias, although a population with a mean slope significantly different than 2 could result in a systematic bias of order 0.02, depending on the slope. Using additional Taylor terms decreases this error as expected, but to implement higher order terms in practice would require more careful accounting of errors on ξ₃ or ξ₄: which is a complexity that we have neglected. This result is in agreement with that of Birrer (2021), who quoted better than 1% accuracy on deflection using the second order expansion under similar conditions.

Fig. B.1. Error resulting from Taylor expansions of deflection of different order, interpreted as error in ξ2, using |r − RE|/RE = 0.3. In the top panel, |r|< RE, while in the bottom panel, |r|> RE. The ξ2 formalism is equivalent to using the second order expansion.

With a maximum error on H₀ ≃ 1% which is not systematic, the Taylor truncation cannot explain the 3% systematic discrepancy in this work. Furthermore, we find that the width of the region probed by the images does not correlate with the value of the H₀ mismatch (as one would expect errors to grow with |r − R_E|). We therefore conclude that the Taylor truncation error is not the cause of the discrepancy shown in Fig. 5.

B.2. Numerical robustness of ξ₂

As this work considers a comparison between similar values of ξ₂, it is important to quantify the effects of numerical precision. We calculate the ξ₂ of the input mass distribution in lenstronomy, which in its current implementation does so by sampling a ring of points at the Einstein radius and evaluating the radial derivatives of the lensing potential at each point, then taking an average over the set for a single value of ξ₂. We quantify the robustness with respect to the choice of the number of points used to sample the ring.

We first compare the calculation to the analytical power law case over the range of slope values in this work, which ranges approximately from ξ₂ = −0.3 to 0.3, based on Fig. 5. This corresponds to γ ∈ [1.85, 2.15]. In the circular case, lenstronomy calculates ξ₂ exactly, regardless of the number of points. In the elliptical case (where we set q = 0.6), we find that the numerical calculation is within 0.002 of the true value of ξ₂, with accuracy which improves as we increase the number of points up to 1000 points, where it remains constant at approximately 0.001. Interestingly, we find that for profiles with steeper slopes, the accuracy on ξ₂ does not continue to improve with an increased number of points, and has a systematic bias at the < 0.001 level.

We also test this effect on the composite profiles, although we cannot analytically calculate the truth value of ξ₂. Testing with the Chameleon+NFWp profiles, we find that the value is not robust (changes by approx. 0.02 or more) with fewer than 1000 points, but the robustness improves with more points. In particular, the change from 3000 points to 10000 points changes the evaluation of ξ₂ by 0.002 (median) with 0.004 standard deviation. We consequently chose to use 3000 points for this work. The Chameleon+NFWm profiles are even more robust, with a change of only 4 × 10⁻⁵ (median) with 0.0002 standard deviation when increasing from 3000 to 10000 points.

From these tests, we conclude that the numerical calculation of ξ₂ is robust to within at worst 0.004 statistical scatter (0.2% for H₀) with a systematic bias likely less than 0.002 (0.1% for H₀). This effect therefore cannot explain the discrepancy between the fitted ξ₂ values originating from the two profile parameterizations in this work.

B.3. Other explored effects

The derivation of ξ₂ comes from the circular limit, and so we consider some of the complexities that arise in the elliptical case. The first main effect we consider is that of propagating the uncertainty on the Einstein radius. The Einstein radius is a circularly averaged quantity and as such we must confirm that this averaging is robust in our noncircular mocks. Like the Einstein radius, ξ₂ is circularly averaged, although it may be that the regions where it is most accurately probed are the image locations rather than a uniform circle. As such, the second main effect we consider is the difference between a circularly averaged ξ₂ and that of an averaging based on the image positions.

An error in the calculation of the Einstein radius would result in evaluating κ_E and ξ₂ at a different location, which could in principle bias the result. With elliptical distributions, it is important to evaluate the effective Einstein radius, which differs from the normalization value that describes the Einstein radius for circular distributions. We evaluate the effective Einstein radius using a grid and by taking radial steps outward, calculating the circle within which the average density is equal to the critical density. We perform some robustness checks by evaluating the effective Einstein radius for several samples within the MCMC chain of the fit and find that the fit Einstein radius always lies within one radial step of the input, in our case approximately 0.013″, which we adopt as our uncertainty in R_E. We then evaluate the change in κ_E by evaluating the local convergence at this new Einstein radius and find that it changes by approximately 1% for both the input profile and the fit profile. The worst case scenario would be that the Einstein radius is recovered on the lower end for one and on the higher end for the other, resulting in the ratio in Eq. 14 being off by approximately 1%, although this appears equally likely for the NFWp parameterization compared to the NFWm parameterization and as such it is difficult to see how this could create a systematic bias. Similarly, we also evaluate the change in ξ₂ of the input profile associated with this change in evaluation radius. For the Chameleon+NFWm, the median change is approximately 0.002, with little spread. For the Chameleon+NFWp, the median change is about half as much, but with significantly more scatter (approximately 0.007). Neither of these effects can explain the observed discrepancy.

Another effect we consider is that ξ₂ itself is a circularly averaged property, which makes sense because the profile used to fit the lens to a particular ξ₂ is a global description. However, we were curious if there could be an effect due to the local measurements where the lensing information most directly probes, that is, at the image positions. As such, we evaluated ξ₂ using the image radial position instead of the Einstein radius for each of the four images, and took the mean of the four evaluations. We also tried taking a weighted average based on the brightness of the images. In either case, the value of ξ₂ changes from the traditional circularly averaged calculation by a median of less than 0.005, with spread of less than 0.025, insufficient to explain the discrepancy.

All Tables

All Figures

Fig. 1. NFW 2D mass isodensity contours of κNFW (top), ellipticity as a function of semimajor axis (middle), and the b4 multipole component as a function of semimajor axis (bottom), in units of NFW scale radius. The NFWp parameterization introduces ellipticity through the lensing potential and is shown in red while the CSE-based NFWm parameterization introduces ellipticity directly in the mass and is shown in gray. The panels from left to right indicate the comparison for three input axis ratios. To match the κNFW ellipticities, we use the procedure in Appendix A, which results in the NFWp potential having an input qψ as indicated above the top panels.

In the text

Fig. 2. Three composite mass profiles drawn from our sample, with the Einstein radius (green, dotted) and NFW scale radius (blue, dotted) indicated. The implementation of these composite profiles is shown for NFWp (red dashed) and NFWm (gray solid) DM profiles. Top row shows isodensity contours; second row shows ellipticity as a function of semimajor axis; third row shows b4 as a function of semimajor axis. The bottom panels show the radial profile over a log scale, defined as the value of κ in a circular annulus at a given radius, rather than with respect to semimajor axis, expressed relative to the NFWm radial profile. The composite profiles are plotted as well as both components individually. Units of distance are now in arcseconds.

In the text

Fig. 3. Example image comparison between NFWp and NFWm ellipticity parameterizations (top), with resulting residuals from the PEMD+shear fits (bottom). Caustics for both parameterizations are shown in the bottom right with the source position rescaled as described in the text and indicated as a cross. This example image comes from a system with moderately large ellipticity with an input axis ratio of qκ = 0.64.

In the text

	Fig. 4. Expected H0 in km s−1 Mpc−1 calculated using the ξ2 calculated from the input mass profiles for 20 mock lenses compared to the value of H0 recovered from a PEMD fit. The dotted line indicates a 1:1 correspondence.
In the text

	Fig. 5. Input ξ2 for the 20 mock lens mass profiles compared to the recovered ξ2 from the PEMD fits, with a 1:1 relation indicated as the dotted line.
In the text

	Fig. 6. Distributions for Rcusp (left) and Rfold (right) for many realized source positions of an example mock.
In the text

	Fig. A.1. Relationship between the mass ellipticity and potential ellipticity for the NFWp profile. Left: ellipticity for different r values. Right: same relation in terms of axis ratios.
In the text

	Fig. B.1. Error resulting from Taylor expansions of deflection of different order, interpreted as error in ξ2, using |r − RE|/RE = 0.3. In the top panel, |r|< RE, while in the bottom panel, |r|> RE. The ξ2 formalism is equivalent to using the second order expansion.
In the text