# Host group degeneracy in gravitational lensing time delay determination of
$H_{0}$

Luca Teodori<EMAIL_ADDRESS>Weizmann Institute, Department of
Particle Physics and Astrophysics, Rehovot, Israel 7610001  Kfir Blum
<EMAIL_ADDRESS>Weizmann Institute, Department of Particle Physics
and Astrophysics, Rehovot, Israel 7610001

###### Abstract

Massive elliptical galaxies, that serve as lenses in gravitational lensing
time delay measurements of the Hubble parameter $H_{0}$, often reside in a
host group. We consider degeneracies in the modeling of the group halo. When
the group effect on imaging can be summarized by its flexion (the next order
term beyond shear in the tidal expansion), the posterior likelihood map can
develop disjoint local minima, associated with an approximate discrete
symmetry of a dominant flexion term. Monte-Carlo Markov Chain (MCMC)
algorithms that are not designed to explore a rich posterior landscape can
miss some of the minima, introducing systematic bias. We study mock data and
demonstrate that the bias in $H_{0}$ can exceed 10%, and pulls the inference
value of $H_{0}$ above its truth value, for a reason that can be traced to the
structure of a mismodeled flexion term. MCMC algorithms that are designed to
cope with a rich posterior landscape can uncover the structure. If the group
is X-ray bright enough, X-ray data may also help to resolve the degeneracy, by
pinpointing the group’s center of mass. Finally, we show that some
implementations in the literature used an inaccurate kinematical prior, mis-
modeling the group velocity dispersion by as much as $20\%$.

###### Contents

  1. I Introduction
  2. II Host group as a core-MSD
  3. III Illustration with mock data
    1. III.1 Mock setup description
    2. III.2 External priors: tracer galaxy kinematics and theoretical input from N-body simulations
    3. III.3 Results
  4. IV Origin of $H_{0}$ bias in a displaced group
  5. V Summary
  6. A Deflection angle expansion, NFW profile
    1. A.1 MSD
    2. A.2 $G$-term degeneracy
    3. A.3 Time delays
    4. A.4 Effect of ellipticity
  7. B Kinematics constraints for host group
    1. B.1 Power law
    2. B.2 NFW
    3. B.3 Kinematical and cosmological priors
    4. B.4 Sample variance
  8. C Estimating the probability of an MCMC to fall into a displaced minimum
  9. D How common is the flexion degeneracy?
  10. E Full corner plots.

##  I Introduction

Gravitationally lensed quasars allow a determination of the Hubble parameter
$H_{0}$ Refsdal (1964); Suyu  _et al._ (2013); Treu and Marshall (2016); Suyu
_et al._ (2017); Grillo  _et al._ (2020), and the results of such measurements
Rusu  _et al._ (2020); Birrer  _et al._ (2019); Shajib  _et al._ (2020); Chen
_et al._ (2019); Wong  _et al._ (2020); Millon  _et al._ (2020) were widely
considered as tests of the cosmological model Riess  _et al._ (2019); Verde
_et al._ (2019); Di Valentino  _et al._ (2021). However, systematic
degeneracies are a limiting factor in the interpretation of lensing data Falco
_et al._ (1985); Kochanek (1991, 2002); Liesenborgs and De Rijcke (2012);
Schneider, Peter and Sluse, Dominique (2013); Kochanek (2021); Fleury  _et
al._ (2021); Teodori  _et al._ (2022). Relaxing some of the modeling
assumptions made in Rusu  _et al._ (2020); Birrer  _et al._ (2019); Shajib
_et al._ (2020); Chen  _et al._ (2019); Wong  _et al._ (2020); Millon  _et
al._ (2020), a possible tension between the value of $H_{0}$ inferred via
lensing and via cosmic microwave background (CMB) and large-scale structure
(LSS) Aghanim  _et al._ (2020); D’Amico  _et al._ (2020); Ivanov  _et al._
(2020) analyses may be attributed to a core feature in the density profile in
or around the lenses Schneider, Peter and Sluse, Dominique (2013); Blum  _et
al._ (2020); Birrer  _et al._ (2020); Blum and Teodori (2021).

A core feature is an approximate mass-sheet degeneracy (MSD) Falco  _et al._
(1985). It could be an intrinsic characteristic of the lens galaxy itself, on
distances of dozens of kpc Blum  _et al._ (2020). However, the effect could
also come from larger scales. On intermediate scales, in between cosmology and
lens internal structure, it is noteworthy that massive galaxies like the
lenses of Rusu  _et al._ (2020); Birrer  _et al._ (2019); Shajib  _et al._
(2020); Chen  _et al._ (2019); Wong  _et al._ (2020); Millon  _et al._ (2020)
are often members of a group111Ref. Wilson  _et al._ (2016) finds that half of
their sample of 26 galaxy lenses can be associated with a group., and the dark
matter halo of a host group may act to some extent as a core-MSD. Ref. Wilson
_et al._ (2016, 2017) studied the impact of host and line of sight (LOS)
groups on lensing systems. Among the systems considered, PG1115+080,
RXJ1131-1231, HE0435-1223, and WFI2033-4723 featured in the $H_{0}$ campaign
of Millon  _et al._ (2020)222At least two additional systems from Millon  _et
al._ (2020) (SDSS1206 Birrer  _et al._ (2019) and DES J0408 Buckley-Geer  _et
al._ (2020)), do not feature in Wilson  _et al._ (2016, 2017), but are known
to also be associated with host groups. Thus at least 6 out of the 7 lensed
quasars from Millon  _et al._ (2020) involve a group.. We show these systems
in Fig. 1. Interestingly, PG1115+080 yields a high central value of $H_{0}$
($81^{+8}_{-7}$ km/s/Mpc) Millon  _et al._ (2020); at the same time, this lens
resides in a massive group, inducing convergence $\kappa_{\rm g}\approx 0.2$
Wilson  _et al._ (2017). The group convergence goes directly into the
inference of $H_{0}$, via $\delta H_{0}/H^{\rm
truth}_{0}\approx-\delta\kappa_{\rm g}$, where $\delta\kappa_{\rm
g}=\kappa_{\rm g}^{\rm truth}-\kappa_{\rm g}^{\rm model}$ is any error in the
model determination of $\kappa_{\rm g}$, and $\delta H_{0}=H_{0}^{\rm
truth}-H_{0}^{\rm model}$. We do not know if the lensing results for
PG1115+080 are indeed biased by its group modeling (and will not make such a
claim in this paper); but clearly, it is important to understand to what
accuracy can lensing analyses determine $\kappa_{\rm g}$.

Figure 1: Lensing-inferred $H_{0}$ (from Millon  _et al._ (2020), with
“$H_{0}$(CMB)” taken from Aghanim  _et al._ (2020)) vs. estimated LOS
convergence (from Wilson  _et al._ (2017)).

In this paper we explore the group effect. The outline and main results are as
follows. In Sec. II we consider simple analytic estimates. We show that when
the distance separating the group’s centroid from the primary lens is large in
comparison to the primary lens’s Einstein angle, attempts to model the group
halo effects based on imaging data suffer from an approximate MSD, where the
mass sheet comes from the group halo itself. This version of the MSD persists
even when second-order tidal effects (flexion terms Goldberg and Bacon (2005);
Bacon  _et al._ (2006)) are clearly detectable. Kinematics and astrometry of
tracer galaxies are needed to break the degeneracy.

In Sec. III we perform numerical mock data analysis, motivated by a realistic
example. We show that the posterior likelihood of the lensing reconstruction
problem exhibits three disjoint minima. These minima reflect an approximate
discrete degeneracy, related to the transformation properties of a dominant
flexion term under coordinate rotations around the primary lens. We show that
a naive implementation of a commonly used Monte-Carlo Markov Chain (MCMC)
algorithm tends to fall into one of the minima, missing the others.
Interestingly, and dangerously, the best fit $H_{0}$ obtained in any one of
the wrong minima is systematically biased high, for a reason that we explain
in Sec. IV. The size of the $H_{0}$ bias can reach $\sim 10\%$. An MCMC
implementation that is especially designed to probe a rich likelihood
landscape, uncovers the full degeneracy structure.

We summarize in Sec. V. Many details are kept to appendices, including
relevant calculations that exist in the literature (for completeness of some
of our arguments), and also a few results that we did not see elsewhere. In
App. A we spell out some details of the lensing potential of a Navarro-Frenk-
White (NFW) profile Navarro  _et al._ (1997). In App. B we review the relation
between group member velocity dispersion and halo model. We note that some
lensing studies implemented an incorrect kinematics prior for the NFW group
model. Depending on selection cuts and some other factors, the systematic
error in the interpretation of group velocity dispersion can reach $\sim
20\%$. In App. C we give a rough estimate of the probability that an MCMC will
actually fall into a wrong minimum; for the example of PG1115+080 Millon  _et
al._ (2020) this probability is not large, on the order of $\sim 10\%$. In
App. D we attempt a (very) crude estimate of the fraction of strong galaxy
lensing systems where flexion degeneracy may be expected to become a concern,
finding this fraction to be in the ballpark of $\sim 10\%-30\%$. In App. E we
collect expanded versions of MCMC triangle plots along with some sanity checks
of our analysis.

We should comment that throughout the analysis, we do not consider stellar
kinematic measurements of the primary lens itself (see Yıldırım  _et al._
(2021); Shajib  _et al._ (2023) for state-of-the-art). Primary lens kinematics
can constrain the group-induced MSD, if it can reach a sensitivity for over-
all scaling of the primary lens mass model at the level of
$1-\delta\kappa_{g}$.

##  II Host group as a core-MSD

In this section we provide simple analytic estimates that clarify the group
effect on the lensing problem. The discussion is useful in understanding
features of the numerical analysis of the next section.

In systems like those considered in Millon  _et al._ (2020), the Einstein
angle of the primary lens is of the order of $\theta_{\rm E}\sim
1^{\prime\prime}$, and most of the imaging information lies at angular
separation $\theta\sim\theta_{\rm E}$ around the lens. This angular scale
projects onto a physical separation of the order of $\lesssim 10$ kpc at the
redshifts $z_{\rm l}\sim 0.1-1$ of typical lenses. In comparison, a typical
separation of any galaxy (with the possible exception of the brightest group
galaxy (BGG)) from the host group’s center of mass is $\gtrsim 100$ kpc, that
is, angular distance $h\gtrsim 10^{\prime\prime}$. Therefore, for a rough
estimate of the impact of a group on imaging analyses, it is sensible to
expand the group’s lensing potential in powers of $\theta/h$.

We adopt complex notation for 2D angle vectors on the sky Schramm and Kayser
(1995); Schneider and Er (2008), defining, e.g.,

$\displaystyle\theta$ $\displaystyle=$
$\displaystyle\theta_{1}+\mathrm{i}\theta_{2},\;\;\;\;\beta\;=\;\beta_{1}+\mathrm{i}\beta_{2},$
(1)

etc. We set the origin of coordinates at the center of the primary lens. With
this formalism, the lensing equation can be written as

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle\theta-\alpha_{\rm
l}(\theta)-\Delta\alpha(\theta)-\kappa_{\rm ext}\theta-\gamma_{\rm
ext}\theta^{*}.$ (2)

Here, $\alpha_{\rm l}(\theta)$ is the deflection angle due to the primary
lens, $\Delta\alpha(\theta)$ is the deflection due to the host group or
cluster, and $\beta$ is the source position. The external convergence and
shear, $\kappa_{\rm ext}$ and $\gamma_{\rm ext}$, contain a combination of
different line of sight (LOS) contributions333For example, in terms of
observer-source, observer-lens, and lens-source LOS terms, we have
$\kappa_{\rm ext}=\kappa_{\rm s}+\kappa_{\rm l}-\kappa_{\rm ls}$. We also
absorb some external convergence and shear terms into the definition of the
“primary lens”, “group”, and “source position” terms Teodori  _et al._
(2022)..

In what follows, when we refer to a host group, we consider the group’s
central dark matter halo, rather than individual member galaxies. For a group
located at center of mass position $\mathrm{e}^{\mathrm{i}\phi}h$, with
$h\gg|\theta|$, $\Delta\alpha$ can be expanded as a power series in $\theta$
Schneider and Er (2008) (see App. A for more details):

$\displaystyle\Delta\alpha(\theta)$ $\displaystyle=$
$\displaystyle\Delta\beta+\Delta\kappa\,\theta+\Delta\gamma\,\theta^{*}$ (3)
$\displaystyle+$
$\displaystyle\frac{1}{4}F^{*}\theta^{2}+\frac{1}{2}F\theta\theta^{*}+\frac{1}{4}G\theta^{*2}+\ldots\,.$

In the case of an axisymmetric group halo profile, it is possible to decompose
the expansion coefficients as

$\Delta\beta=\Delta\beta_{0}\mathrm{e}^{\mathrm{i}\phi},\
\Delta\gamma=\Delta\gamma_{0}\mathrm{e}^{2\mathrm{i}\phi},\
F=F_{0}\mathrm{e}^{\mathrm{i}\phi},\ G=G_{0}\mathrm{e}^{3\mathrm{i}\phi},$ (4)

where $\Delta\beta_{0},\Delta\gamma_{0},F_{0},G_{0}$ are real numbers that
depend on $h$ but not on $\phi$. Note that $\Delta\kappa$ is independent on
$\phi$.

Axisymmetry is mildly broken in realistic elliptic profiles; we discuss this
point in App. A.4. However, as long as the ellipticity is small, the $\phi$
dependence of $\Delta\beta_{0},\Delta\gamma_{0},F_{0},G_{0}$ is weak, and does
not affect our main point.

Of the expansion coefficients in Eq. (3), $\Delta\beta$ is degenerate with
$\beta$, which is a free parameter in the modeling, and has no effect on time
delays. $\Delta\kappa$ is degenerate with $\kappa_{\rm ext}$, a free
parameter.444A cosmological prior can sometimes be imposed on $\kappa_{\rm
ext}$, by comparing the system’s field to N-body simulations. Varied along the
MSD, the combination $\Delta\kappa+\kappa_{\rm ext}$ is invisible to imaging,
but affects time delays. $\Delta\gamma$ is degenerate with $\gamma_{\rm ext}$,
a free parameter. Thus, as far as imaging data is concerned, the leading order
effect that constrains the group comes from the MSD-invariant combinations
$F/(1-\Delta\kappa-\kappa_{\rm ext})$ and $G/(1-\Delta\kappa-\kappa_{\rm
ext})$, related to what is known in the literature as reduced flexion Goldberg
and Bacon (2005); Bacon  _et al._ (2006).

An attempt to constrain the group model via imaging data, suffers from the
following three main difficulties. The first obvious point is that the flexion
terms are small; in the $1/h$ expansion, the flexion terms are parametrically
suppressed as $F,G\propto\Delta\kappa/h$. The second point is that direct
modeling of the group halo can still leave room for a residual MSD: a model to
capture correctly the imaging distortion produced by the flexion, while
mismodeling the convergence. The third point, which may be the most important
in practice, is that the posterior likelihood of the group halo model exhibits
a discrete approximate degeneracy, related to the $2\pi/3$ phase degeneracy of
the $G$ term (see Eq. (4)).

It is useful to study a simple example. Consider a group halo described by an
isotropic power-law (PL) density profile with 3D slope $\gamma_{\rm g}$ and
Einstein angle $\theta_{\rm Eg}$. In this case, the coefficients of Eqs. (3)
and (4) are given by

$\displaystyle\Delta\kappa$ $\displaystyle=$ $\displaystyle\frac{3-\gamma_{\rm
g}}{2}\left(\frac{\theta_{\rm Eg}}{h}\right)^{\gamma_{\rm
g}-1},\;\;\;\Delta\gamma_{0}=\frac{1-\gamma_{\rm g}}{3-\gamma_{\rm
g}}\Delta\kappa,$ (5) $\displaystyle F_{0}$ $\displaystyle=$
$\displaystyle\left(\gamma_{\rm
g}-1\right)\frac{\Delta\kappa}{h},\;\;\;G_{0}=-\frac{\gamma_{\rm
g}^{2}-1}{3-\gamma_{\rm g}}\frac{\Delta\kappa}{h},$ (6)
$\displaystyle\Delta\beta$ $\displaystyle=$
$\displaystyle-\frac{2\,h}{3-\gamma_{\rm g}}\Delta\kappa.$ (7)

For $\gamma_{\rm g}\approx 2$ (singular isothermal sphere (SIS)), the line of
sight velocity dispersion (LOSVD) in the center of the group is related to the
group’s Einstein radius via555See App. B for a review of the derivation. For
$\gamma_{\rm g}$ close, but not equal to 2, the RHS of Eq. (8) is rescaled by
an $\mathcal{O}(1)$ factor, equal to $\sim 0.5$ ($\sim 1.5$) for $\gamma_{\rm
g}=1.5$ ($\gamma_{\rm g}=2.2$), and a weak dependence arises on the scale
radius of the profile.

$\displaystyle\theta_{\rm Eg}$ $\displaystyle\approx$ $\displaystyle\frac{4\pi
d_{\rm ls}}{d_{\rm s}}\sigma^{2}_{\rm
los}\,\approx\,4.5^{\prime\prime}\,\frac{d_{\rm ls}}{d_{\rm
s}}\left(\frac{\sigma_{\rm los}}{400~{}\rm km/s}\right)^{2}.$ (8)

We can therefore estimate the convergence,

$\displaystyle\Delta\kappa_{\rm SIS}$ $\displaystyle\approx$ $\displaystyle
0.5\frac{d_{\rm ls}}{d_{\rm s}}\left(\frac{\sigma_{\rm los}}{400~{}\rm
km/s}\right)^{2}\left(\frac{20^{\prime\prime}}{h}\right),$ (9)

and to the flexion terms,

$F_{0,\rm
SIS}\,\approx\,\frac{0.025}{1^{\prime\prime}}\left(\frac{20^{\prime\prime}}{h}\right)\Delta\kappa_{\rm
SIS},\;\;\;G_{0,\rm SIS}\,\approx\,-3\,F_{0,\rm SIS}.$ (10)

For example, $\Delta\kappa\approx 0.1$ produced by a SIS group with $d_{\rm
ls}/d_{\rm s}=0.5$, $\sigma_{\rm los}=400$ km/s, and $h=50^{\prime\prime}$,
would cause (if not modeled) a $\sim 10\%$ bias in the inference of $H_{0}$,
while the imaging distortion produced by the group’s flexion field at
$\theta\sim 1^{\prime\prime}$ would only amount to $\delta\theta\approx
0.0075^{\prime\prime}$, typically dominated by the $G$ term. Next order terms
in the expansion (beyond the flexion) are further suppressed by $\sim
1^{\prime\prime}/h\sim\mathcal{O}(1/10)$, and can be neglected.

This discussion highlights the obvious hierarchy between flexion and
convergence, but as was mentioned earlier, there is also MSD. In Eqs. (5-6),
the leading effect of the MSD (at small $\Delta\kappa\ll 1$) is seen by
noticing that even if one can fix $F$ and $G$ exactly from imaging data, this
still allows $\Delta\kappa$ to vary freely, as long as $h$ is varied
simultaneously with $\Delta\kappa\propto h$. Thus, unless we have good
external prior data on the group’s center of mass (parameterised by $h$),
determining the flexion from imaging data alone does not fix the convergence.
(More precisely, the degeneracy is somewhat regulated by the fact that the
true MSD invariant quantities are $F/(1-\Delta\kappa-\kappa_{\rm ext})$ and
$G/(1-\Delta\kappa-\kappa_{\rm ext})$, rather than $F$ and $G$ themselves. The
exact MSD is then captured via the freedom to adjust $\Delta\kappa$ and $h$
while keeping $(\Delta\kappa/h)/\left(1-\kappa_{\rm ext}-\Delta\kappa\right)$
constant.)

Finally, another important point is the phase degeneracy of the $G$ term. This
degeneracy means that models of the halo in which the direction to the group’s
center is changed by $\phi\to\phi\pm 2\pi/3$ yield identical $G$ terms.
Although the rotated models do not reproduce the truth value of the $F$ term,
the $G$ degeneracy can produce isolated local minima in the posterior
likelihood, that can trap unweary MCMC chains. Interestingly, the offset $F$
term in the “wrong” minima causes the model to pull towards a biased estimate
of the group’s convergence, thereby biasing $H_{0}$. This issue is seen to be
an important point in the next section.

##  III Illustration with mock data

The host group is often modeled explicitly if the system is known to reside in
a group (see e.g. Millon  _et al._ (2020); Momcheva  _et al._ (2006); Wilson
_et al._ (2016, 2017)). We now study such modeling using mock data. Our
implementation is based on the package lenstronomy Birrer and Amara (2018);
Birrer  _et al._ (2015, 2021).

###  III.1 Mock setup description

We chose PG1115+080 as a reference object to guide our study. Ref. Millon  _et
al._ (2020) inferred $H_{0}=81^{+8}_{-7}$ km/s/Mpc from this system, quite
high compared to the CMB value Aghanim  _et al._ (2020). At the same time, the
system is known to reside in a group Wilson  _et al._ (2017) with $\sigma_{\rm
los}=390^{+50}_{-60}~{}{\rm km/s}$ (an earlier study found $\sigma_{\rm
los}=440^{+90}_{-80}~{}{\rm km/s}$ Momcheva  _et al._ (2006)), estimated to
induce $\kappa_{\rm g}\sim 0.2$. The group’s projected center of mass is not
far from the primary lens ($h$ of the order of $10^{\prime\prime}$). This
suggests that flexion terms are probably not negligible for this system,
distorting the image on scales $\Delta\theta\sim\kappa_{\rm g}\theta_{\rm
E}^{2}/h\sim 0.01^{\prime\prime}$. An image of the field is shown in Fig. 2.

We consider the following mock setup. For the primary lens, we consider an
elliptic power-law density profile with 3D slope $\gamma=2.17$, Einstein angle
$\theta_{\rm E}=1.08^{\prime\prime}$, and ellipticity parameters
$e_{1}=-0.2,\,e_{2}=0.05$, corresponding to $q=(1-e)/(1+e)\approx 0.66$, where
$e=\sqrt{e^{2}_{1}+e^{2}_{2}}$. For the group, we consider an elliptic NFW
profile Jing and Suto (2002); Vega-ferrero  _et al._ (2017) with
$e_{1}=-0.07$, $e_{2}=0.03$, compatible with the findings of Schrabback  _et
al._ (2021). Fig. 3 illustrates the setup, including the truth position of the
group halo center of mass. This setup could mimic PG1115+080 (Fig. 2) if the
BGG – or the X-ray blob found by Grant  _et al._ (2004) – happens to indicate
the group’s center of mass position.

Figure 2: The field of PG1115+080 (image from SIMBAD Wenger  _et al._ (2000)).
The primary lens, GL, is marked by a cross, surrounded by the quasar images.
White ellipse shows the 68% CL prior derived in Wilson  _et al._ (2016) and
used in Chen  _et al._ (2019) to constrain the group’s center of mass. Green
and blue ellipses give a rough illustration of the X-ray emission associated
to the group, found in Grant  _et al._ (2004) and Fassnacht  _et al._ (2008),
respectively, using different subtraction schemes for the quasar emission.
Figure 3: Illustration of the mock setup. Contours show $\log_{10}$ of the
convergence for the primary lens and the NFW halo in black and red,
respectively .

Following Chen  _et al._ (2019), we choose our inference pipeline to include
only a simple spherical group model, ignoring group halo ellipticity in the
modeling.

###  III.2 External priors: tracer galaxy kinematics and theoretical input
from N-body simulations

As in Chen  _et al._ (2019), we include external priors on the group halo,
dictated by cosmological N-body simulations and by kinematics data. We defer
most of the details to App. B, notably Secs. B.2 and B.3. However, we would
like to point out a potential inaccuracy in the kinematics analyses of some
previous works.

Galaxy groups are often assumed to follow the NFW density profile Navarro  _et
al._ (1997),

$\displaystyle\rho(r)=\frac{\rho_{0}R^{3}_{\rm s}}{r(r+R_{\rm s})^{2}},$ (11)

In terms of the parameters $\rho_{0}$ and $R_{\rm s}$, the LOSVD of tracer
galaxies (with number density assumed to follow the same profile as the dark
matter density) can be expressed as

$\displaystyle\sigma_{\rm los}^{2}(\theta)$ $\displaystyle=$ $\displaystyle
G\rho_{0}R_{\rm s}^{2}\,f\left(\frac{\theta}{\theta_{\rm s}}\right),$ (12)

where $f(a)$ is a dimensionless function derived in Sec. B.2.

While Eq. (12) (averaged as needed within some aperture cut of the
observations) gives the correct translation between LOSVD data and the NFW
model parameters, Ref. Chen  _et al._ (2019) (after Koranyi and Geller (2000);
Wong  _et al._ (2011)) considered a different expression as a proxy for the
LOSVD data:

$\displaystyle\bar{\sigma}^{2}$ $\displaystyle=$ $\displaystyle\frac{GM_{\rm
vir}}{3R_{\rm vir}},$ (13)

where $R_{\rm vir}$ and $M_{\rm vir}=M(R_{\rm vir})$ are the virial radius and
virial mass, and $c_{\rm vir}=R_{\rm vir}/R_{s}$ is the NFW concentration
parameter.In App. B we show that identifying the quantity $\bar{\sigma}^{2}$
with the observable $\sigma^{2}_{\rm los}$ introduces an error of up to $\sim
20\%$ (the precise error depends on the analysis aperture).

One reason for introducing the auxiliary quantities $R_{\rm vir}$ and $M_{\rm
vir}$, is that N-body simulations provide theoretically-motivated priors that
are often presented in terms of these quantities Eke  _et al._ (2001); Maccio’
_et al._ (2008). These cosmological priors are, however, quantitatively and
conceptually decoupled from the kinematics data interpretation. Instead, as we
review in App. B, the cosmological priors dictate a certain redshift-dependent
relation between the NFW parameters $\rho_{0}$ and $R_{\rm s}$. There is no
obstacle to implement this theoretical prior while still maintaining the
correct kinematical expression, Eq. (12).

In our main implementation of the MCMC, we set the standard deviation for
$\sigma_{\rm los}$ to $120$ km/s. This doubles the nominal uncertainty quoted
by Wilson  _et al._ (2017) for the group of PG1115+080, but we believe that
such cautionary procedure is reasonable. To be clear, we are not suggesting to
doubt the observational LOSVD from Wilson  _et al._ (2017). Rather, the
uncertainties we worry about concern the theoretical interpretation within
simplified halo models. The systematic error due to using $\bar{\sigma}$
instead of $\sigma_{\rm los}$, as done in Chen  _et al._ (2019), is of the
order of $\delta\sigma_{\rm los}^{2}/\sigma_{\rm los}^{2}\sim 20\%$, so
$\delta\sigma_{\rm los}/\sigma_{\rm los}\sim 10\%$, or $\delta\sigma_{\rm
los}\sim 40$ km/s, quite comparable to the “bare” observational uncertainty.
There are additional plausible errors: the group is likely to be aspherical
Schrabback  _et al._ (2021), the velocity distribution need not be isotropic,
and the group may not be fully virialised. Each of these could cause
systematic shifts of tens of percent in the kinematics interpretation.

For completeness, in App. E we check how a LOSVD uncertainty of $60$ km/s
changes the results. We find the difference to be quantitatively insignificant
for our main results.

###  III.3 Results

First, to obtain a global view of the posterior “landscape”, we run the zeus
MCMC algorithm Karamanis and Beutler (2021); Karamanis  _et al._ (2021), which
is designed to cope with multiple likelihood minima.

Figure 4: Mock analysis. Truth set-up: elliptic group located at $c_{x}^{\rm
nfw}=18^{\prime\prime}$, $c_{y}^{\rm nfw}=-10^{\prime\prime}$. Top left panel:
zeus run, exposing the global likelihood landscape. Other panels: emcee runs,
falling into different local minima. The emcee runs are initiated with
different priors for the position of the group center, keeping the same
Gaussian standard deviation of $16^{\prime\prime}$. See also Fig. 5 .

We show the main result in the top-left panel of Fig. 4. Notice the 3 disjoint
minima in the posterior likelihood in the $c_{x}^{\rm nfw}-c_{y}^{\rm nfw}$
panel, where $c_{x,y}^{\rm nfw}$ are the angular coordinates of the halo
center of mass. The origin of this threefold degeneracy is the $\phi\to\phi\pm
2\pi/3$ degeneracy of the flexion $G$ term. This run has no prior on
$c_{x}^{\rm nfw}$ and $c_{y}^{\rm nfw}$. More details and comprehensive MCMC
corner plots are provided in App. E.

The local minima could trap an MCMC, if the scanning algorithm is not suited
to probe multimodal posteriors. To see this, we repeat the analysis, this time
using the emcee algorithm Foreman-Mackey  _et al._ (2013). The results are
shown in the top-right and bottom panels of Fig. 4. Indeed, emcee tends to
discover only one of the local minima, missing the others. In the top-right
and bottom-right panels, emcee converges on a biased group halo position. The
biased local minima yield an $H_{0}$ posterior that is on the high side,
pulling $H_{0}$ above its truth value by $\sim 3\sigma$ and $\sim 1.8\sigma$,
respectively. The convergence pattern associated to the three minima is shown
in Fig. 5.

Performing more MCMC runs, we find that emcee consistently converges into just
one of the minima, missing the others. The choice of the minimum found by
emcee mostly depends on the initial position of the MCMC walkers in the
$c_{x}^{\rm nfw}-c_{y}^{\rm nfw}$ parameter space. In our runs, the initial
walker allocation is guided by a Gaussian prior with a similar 68%CL radius as
that used in Chen  _et al._ (2019) (see Fig. 2), and with centers shown by
star symbols in Fig. 5.

After performing a series of numerical trials we emphasize that, at least for
the given, rather broad (but realistic) width of the prior, the starting point
of the MCMC walkers appears to be a more significant factor in determining
which of the three minima traps the chain, than the prior center itself. These
considerations coincide if the MCMC walkers are initiated at the prior center.

Given the above discussion, we can make a rough estimate of the probability of
the MCMC to fall into a displaced minimum, by matching it with the probability
of the group center prior to be nearer a false minimum than the truth one. In
Sec. C we estimate this probability from mock realizations of samples of
tracer galaxies. The result depends on the analyses aperture, the number of
galaxies in a sample, and the underlying group profile. For 13 members in the
NFW profile, with an aperture of $\sim 10R_{s}$ (i.e., around 3 virial radii),
the false minimum probability we find is $\sim 10\%$.

Figure 5: Illustration of the three posterior likelihood minima probed by the
MCMC inference. Contours show $\log_{10}$ convergence isocontours, for the 3
group halo solutions (red, orange, and blue, respectively). Dots and stars
highlight the inferred source position and the center of the group prior,
respectively. Green squares denote the position of the images. The parameter
values used in the plot are read from the emcee run best fit, see Fig. 4 .

We comment that X-ray data may help to pinpoint the centroid of a massive
group, resolving the threefold degeneracy. For example, Ref. Grant  _et al._
(2004) found an X-ray blob that appears to be centered around the BGG of the
host group of PG1115+080 (green contour in Fig. 2). If the X-ray emission can
be associated with the group’s centroid, it may provide a narrower prior than
that derived from tracer galaxies. The X-ray analysis may be complicated by
blending with the lensed quasar: using a different method to mask the quasar,
Ref. Fassnacht  _et al._ (2008) found a shifted, more extended, and brighter
group emission. Even so, it seems plausible that X-ray data could help
narrowing the group’s centroid prior666In principle, X-ray luminosity
measurements can also constrain the group’s velocity dispersion and mass. In
practice, however, the limited accuracy of X-ray luminosity measurements
Fassnacht  _et al._ (2008), combined with the observed scatter in
luminosity–mass or luminosity–dispersion correlations Lovisari  _et al._
(2021), suggests that the corresponding constraints may not lead to
significant improvement over kinematics alone..

It is natural to ask whether the discrete modeling bias tends to come along
with an anomalously large external shear estimate. While this could indeed
happen, and might serve as a useful “alarm bell” if the $H_{0}$ bias is larger
than 10% or so, our numerical examples also show false minima solutions with
acceptable values of external shear. For example, Fig. 11 shows a false
minimum solution with $\sim$15% bias in $H_{0}$ but with external shear values
that are compatible with those found in Chen  _et al._ (2019) for PG1115+080
(see Figs. 7,8 there).

Finally, we emphasize again that PG1115+080 was selected to guide our mock
specifically because of its massive, near-by group association. This set-up,
while we believe it deserves study on its own right, may be an un-
representative outlier among lensed quasar systems, and it is therefore
natural to ask just how uncommon it is. A detailed analysis of the fraction of
systems that may exhibit flexion degeneracy is beyond the scope of this work,
but in App. D we attempt a crude estimate using the data from Wilson  _et al._
(2017). Our results suggest that detectable flexion distortion from LOS groups
(including, not limited to, the host of the primary lens) may affect $\sim
10\%-30\%$ of lensed quasar systems.

##  IV Origin of $H_{0}$ bias in a displaced group

As we explained, the three-fold approximate degeneracy manifest in Fig. 4 is
due to the behavior under rotations of the $G$ flexion. The $F$ flexion,
however, behaves like a vector. Hence, when the inference falls into a wrong
minimum, it attempts to minimize the difference between the truth $F$
deflection and the wrong-minimum inference value of $F$. This can be achieved
by keeping $G_{\rm eff}$ the same, but reducing $F_{\rm eff}$ as much as
possible, where

$F_{\rm eff}=\frac{F_{0}}{1-\Delta\kappa}\ ,\ G_{\rm
eff}=\frac{G_{0}}{1-\Delta\kappa}$ (14)

are the reduced flexion terms (for clarity, here we omit the cosmological term
$\kappa_{\rm ext}$). To see this point, denote the reduced flexion of the
inference model by $F_{\rm eff}$, and denote the truth flexion by $F^{(\rm
truth)}_{\rm eff}$.777The truth and model values of the $G$ term are assumed
to approximately coincide, $G^{(\rm truth)}_{\rm eff}\approx G_{\rm eff}$. We
expect the MCMC to minimize

$\displaystyle\absolutevalue{F_{\rm eff}^{(\rm truth)}-F_{\rm
eff}\mathrm{e}^{\frac{2\mathrm{i}n\pi}{3}}}=$ $\displaystyle\left|F^{(\rm
truth)}_{\rm eff}\right|\sqrt{1-2\frac{F_{\rm eff}}{F^{(\rm truth)}_{\rm
eff}}\cos\left(\frac{2n\pi}{3}\right)+\left(\frac{F_{\rm eff}}{F^{(\rm
truth)}_{\rm eff}}\right)^{2}},\;\;$ (15)

with $n=\pm 1$ selecting the position of the false minimum. Since $\cos(\pm
2\pi/3)=-1/2$, the above expression is minimized for $F_{\rm eff}=0$.

Consider the PL model of Sec. II. In this model, we have $F_{\rm eff}/G_{\rm
eff}=(\gamma_{\rm g}-3)/(\gamma_{\rm g}+1)$. Therefore, for $1<\gamma_{\rm
g}<3$ (the range of interest), decreasing $|F_{\rm eff}|$ at fixed $G_{\rm
eff}$ entails increasing $\gamma_{\rm g}$, while adjusting the other
parameters of the model so as to keep $G_{\rm eff}$ constant. Those other
parameters were introduced in Eqs. (5-6) as $h$ and $\theta_{\rm Eg}$, but we
can equally well replace $\theta_{\rm Eg}$ by $\Delta\kappa$. Now, we have
$G_{\rm eff}=-\frac{\gamma^{2}_{\rm g}-1}{3-\gamma_{\rm
g}}\frac{\Delta\kappa}{h\left(1-\Delta\kappa\right)}$. The $\gamma_{\rm
g}$-dependent factor, $\frac{\gamma^{2}_{\rm g}-1}{3-\gamma_{\rm g}}$,
increases with increasing $\gamma_{\rm g}$; to compensate for this and keep
$G_{\rm eff}$ constant, the factor
$\frac{\Delta\kappa}{h\left(1-\Delta\kappa\right)}$ needs to decrease. For
small $\Delta\kappa\ll 1$, this means that near any one of the displaced
likelihood minima, the MCMC will attempt to decrease $\Delta\kappa/h$ in
comparison to its truth value. Part of this adjustment entails decreasing the
model value of $\Delta\kappa$, which therefore biases $H_{0}$ high, as $\Delta
H_{0}/H_{0}\approx-\left(\Delta\kappa^{(\rm truth)}-\Delta\kappa^{(\rm
model)}\right)$.

In App. A we show that a similar analysis holds also for the NFW model, used
in the MCMC implementation: also in that case, falling into a displaced group
minimum causes the fit to underestimate of $\Delta\kappa$, leading to an
overestimate of $H_{0}$. This analysis clarifies the trend seen in Fig. 4.

##  V Summary

Lens galaxies in quasar lensing time delay measurements are often members of
galaxy groups, that must be modeled for an accurate determination of $H_{0}$.
The group modeling exhibits approximate versions of the MSD (Sec. II).
Essentially, it is a displaced-center version of the core-MSD considered in
Schneider, Peter and Sluse, Dominique (2013); Blum  _et al._ (2020).

At leading order in the tidal approximation, the group halo enters imaging
through the flexion. We showed an approximate threefold discrete modeling
degeneracy, associated with rotating the assumed position of the group
centroid by an angle of $2\pi/3$ around the primary lens (Secs. II and III).
This produces a posterior likelihood with three disjoint minima. MCMC
algorithms that fail to expose this structure may fall into a displaced
minimum. The inferred value of $H_{0}$ found in a displaced minimum is
systematically biased high (Sec. IV). Using numerical mock data experiments
motivated by a realistic system, we demonstrated that the $H_{0}$ bias can
reach $\sim 10\%$.

The choice of the minimum detected by the MCMC can strongly depend on the
starting position of the walkers in the space of group centroid coordinates.
If the starting point is chosen as the centroid prior center, then the
probability for the MCMC to land in a displaced minimum may be rather small.
For a sample of 13 tracer galaxies (relevant for PG1115+080) with an aperture
of about 3 virial radii in an NFW halo, this probability is $\sim 10\%$.

Our analysis suggests the following recommendations.

  1. 1.

Bayesian cosmography analyses should explore the full posterior likelihood
landscape. Awareness of this possible three disjoint minima structure, if not
already there, is needed.

  2. 2.

X-ray data may help to pinpoint the centroid of a massive group, resolving the
degeneracy.

  3. 3.

As an aside, we note that some cosmography analyses (e.g. Chen  _et al._
(2019)) used an incorrect kinematics prior to constrain the group model. The
error in the interpretation of tracer galaxy velocity dispersion depends on
the analysis aperture, and can reach $\sim 20\%$. A correct version of the
kinematics prior is reviewed in App. B.

The number of strong lensing time delay systems is expected to increase by
more than an order of magnitude in the near future Oguri and Marshall (2010);
Liao  _et al._ (2015); Dobke  _et al._ (2009), an important step towards
possibly reaching a few percent lensing determination of $H_{0}$ Birrer and
Treu (2021). This program could be further assisted by many resources Treu
_et al._ (2018); Bacon  _et al._ (2020); Blanchard  _et al._ (2020); Abell
_et al._ (2009); Gardner  _et al._ (2006). Our study highlights some pitfalls
(and suggests solutions) that need to be taken into account if the precision
goal for $H_{0}$ should also be accurate.

###### Acknowledgements.

We thank Simon Birrer, Marko Simonović, and Raphael Flauger for useful
discussions. This work made use of the following public software packages:
lenstronomy Birrer and Amara (2018); Birrer  _et al._ (2015, 2021), emcee
Foreman-Mackey  _et al._ (2013), zeus Karamanis and Beutler (2021); Karamanis
_et al._ (2021), corner Foreman-Mackey (2016), astropy Robitaille  _et al._
(2013); Price-Whelan  _et al._ (2018). This research has made use of the
SIMBAD database, operated at CDS, Strasbourg, France Wenger  _et al._ (2000).
This work was supported by the Israel Science Foundation grant 1784/20, and by
MINERVA grant 714123. LT wishes to acknowledge association with the
International Helmholtz-Weizmann Research School for Multimessenger Astronomy.

##  Appendix A Deflection angle expansion, NFW profile

We start with some general preliminaries; for the discussion of the NFW
profile, the reader can skip to Eq. (11).

For our purpose, which is to analyze the lensing equation in the vicinity of a
particular galaxy member of the group, it is convenient to use a coordinate
system that is centered on the primary lens galaxy, and displaced from the
group center of mass by a separation angle $\vec{h}$. In these coordinates the
group lensing potential reads

$\displaystyle\Psi_{\rm halo}(\theta_{1},\theta_{2}):=\Psi_{\rm
halo}(|\vec{\theta}-\vec{h}|).$ (16)

We can expand with respect to the small parameters $\theta_{1}/h,\
\theta_{2}/h$,

$\Psi_{\rm halo}(\theta_{1},\theta_{2})\simeq\Psi_{\rm
halo}(0,0)+\gradient_{\vec{\theta}}\Psi_{\rm
halo}(0,0)\cdot\vec{\theta}+\ldots\ .$ (17)

Recall that the lensing potential and the deflection angle are related as
$\gradient\Psi=\vec{\alpha}$; in particular, we can write the lens equation as

$\beta_{i}=\theta_{i}-\partial_{i}\Psi_{\rm
lens}(\vec{\theta})-\partial_{i}\Psi_{\rm
halo}(0,0)-\partial_{i}\partial_{j}\Psi_{\rm halo}(0,0)\theta_{j}+\ldots\ .$
(18)

As discussed on the main text, in complex notation $\vec{h}\to
h\mathrm{e}^{\mathrm{i}\phi}$. The lensing potential $\Psi$ in this formalism
is real. One can obtain the complex deflection angle by means of the
derivative operator

$\gradient_{\rm
c}:=\partial_{1}+\mathrm{i}\partial_{2}\implies\alpha=\gradient_{\rm c}\Psi\
.$ (19)

In the expansion of the lensing potential, it is easy to see that

$\Delta\beta=\gradient_{\rm c}\Psi_{\rm halo},\
\Delta\kappa=\frac{1}{2}\gradient_{\rm c}\gradient^{*}_{\rm c}\Psi_{\rm
halo},\ \Delta\gamma=\frac{1}{2}\gradient_{\rm c}\gradient_{\rm c}\Psi_{\rm
halo},$ (20)

where the derivatives are computed at the origin. Analogously, we can express
third order derivatives of the lensing potential as derivatives of the shear,

$G=\gradient_{\rm c}\Delta\gamma,\ F=\gradient^{*}_{\rm c}\Delta\gamma,$ (21)

where $G$, $F$ are the flexion terms.

In the axisymmetric case, a group halo center at $\phi=0$ and fixed $h$ would
yield the same expansion coefficients of a halo located at a generic $\phi$,
granted that we rotate our coordinate system accordingly. We can pass to
cylindrical coordinates888With this definition, when $\theta=0$, $\phi$
corresponds to the angle between the origin and the position of the halo
center.

$\displaystyle\vec{h}-\vec{\theta}=r(\cos\phi,\sin\phi)^{\top},$ (22)

and write

$\partialderivative{\theta_{i}}\Psi_{\rm
halo}=\partialderivative{r}{\theta_{i}}\partialderivative{r}\Psi_{\rm
halo}(r)=\frac{\theta_{i}-h_{i}}{r}\Psi^{\prime}_{\rm halo}(r).$ (23)

In complex notation, we thus have

$\gradient_{\rm c}\Psi_{\rm
halo}=-\mathrm{e}^{\mathrm{i}\phi}\Psi^{\prime}_{\rm halo}(r).$ (24)

Notice the minus sign, which is there due to our choice of coordinates; it is
easy to see that

$\partialderivative{\phi}{\theta_{1}}=\frac{\sin\phi}{r}\ ,\
\partialderivative{\phi}{\theta_{2}}=-\frac{\cos\phi}{r},$ (25)

which has opposite sign with respect to the usual choice of cylindrical
coordinate. We can thus write

$\gradient_{\rm
c}=-\mathrm{e}^{\mathrm{i}\phi}\quantity(\partialderivative{r}+\frac{\mathrm{i}}{r}\partialderivative{\phi}),$
(26)

so that

$\displaystyle\Delta\kappa=\frac{1}{2}\quantity(\Psi^{\prime\prime}_{\rm
halo}+\frac{\Psi^{\prime}_{\rm halo}}{r}),$ (27)
$\displaystyle\Delta\gamma=\frac{\mathrm{e}^{2\mathrm{i}\phi}}{2}\quantity(\Psi^{\prime\prime}_{\rm
halo}-\frac{\Psi^{\prime}_{\rm halo}}{r}),$ (28) $\displaystyle
G=-\frac{\mathrm{e}^{3\mathrm{i}\phi}}{2}\quantity(\Psi^{\prime\prime\prime}_{\rm
halo}-3\frac{\Psi^{\prime\prime}_{\rm halo}}{r}+3\frac{\Psi^{\prime}_{\rm
halo}}{r^{2}}),$ (29) $\displaystyle
F=-\frac{\mathrm{e}^{\mathrm{i}\phi}}{2}\quantity(\Psi^{\prime\prime\prime}_{\rm
halo}+\frac{\Psi^{\prime\prime}_{\rm halo}}{r}-\frac{\Psi^{\prime}_{\rm
halo}}{r^{2}})=-\mathrm{e}^{\mathrm{i}\phi}\derivative{\Delta\kappa}{r}.$ (30)

We now specify to the NFW model. Here we discuss the spherical model without
ellipticity, commenting on ellipticity App. A.4. Using a coordinate system
centered on the group, the NFW lensing potential is Bartelmann (1996)

$\displaystyle\Psi_{\rm NFW}(\theta)=2\tilde{\kappa}\,\theta^{2}_{\rm
s}\quantity(\log^{2}\quantity(\frac{\theta}{2\theta_{\rm
s}})-\rm{arccosh}^{2}\quantity(\frac{\theta_{\rm s}}{\theta})),$ (31)
$\displaystyle\theta_{\rm s}:=\frac{R_{\rm s}}{d_{\rm l}},\
\tilde{\kappa}:=\frac{\rho_{0}R_{\rm s}}{\Sigma_{\rm c}},\ \Sigma_{\rm
c}=\frac{d_{\rm s}}{4\pi Gd_{\rm l}d_{\rm ls}}.$ (32)

Defining $x=\theta_{\rm s}/h$, we have:

$\displaystyle\Delta\beta_{0}=4h\,\tilde{\kappa}\,x^{2}\quantity(\log
2x-\frac{\mathrm{arccosh}\,x}{\sqrt{1-1/x^{2}}})\ ,$ (33)
$\displaystyle\Delta\kappa=\tilde{\kappa}\frac{2x^{2}}{1-x^{2}}\quantity(1-\frac{\mathrm{arccosh}\,x}{\sqrt{1-1/x^{2}}})\
,$ (34)
$\displaystyle\Delta\gamma_{0}=\tilde{\kappa}\frac{2x^{2}}{1-x^{2}}\Bigg{(}1+\frac{(2x^{2}-3)\mathrm{arccosh}\,x}{\sqrt{1-1/x^{2}}}$
$\displaystyle+2\quantity(1-x^{2})\log(2x)\Bigg{)},$ (35) $\displaystyle
F_{0}=\frac{\tilde{\kappa}}{h}\frac{2x^{2}}{(1-x^{2})^{2}}\quantity(2+x^{2}-\frac{3\mathrm{arccosh}\,x}{\sqrt{1-1/x^{2}}}),$
(36) $\displaystyle
G_{0}=\frac{\tilde{\kappa}}{h}\frac{2x^{2}}{(1-x^{2})^{2}}\Bigg{(}-\frac{15-20x^{2}+8x^{4}}{\sqrt{1-1/x^{2}}}\mathrm{arccosh}\,x$
$\displaystyle+6-3x^{2}+8\left(1-x^{2}\right)^{2}\log(2x)\Bigg{)}.$ (37)

Note that applying these expressions for $x<1$ requires using the logarithmic
definition of $\rm{arccosh}(z)=\ln(z+\sqrt{z^{2}-1})$ and allowing complex
$z$.

It is useful to define $f_{\kappa}(x),\,f_{G}(x),\,f_{F}(x)$ such that
$\Delta\kappa=\tilde{\kappa}\,f_{\kappa}(x)$,
$F_{0}=\frac{\tilde{\kappa}}{h}\,f_{F}(x)$,
$G_{0}=\frac{\tilde{\kappa}}{h}\,f_{G}(x)$. These functions are shown in Fig.
6.

We can use these expressions to analyze the modeling constraints obtained from
imaging data. As we already remarked in Sec. II, of the expansion terms,
$\Delta\beta_{0}$ is exactly degenerate with the unknown source position;
$\Delta\kappa$ is degenerate with external convergence, and can be absorbed by
the MSD (changing the inference of $H_{0}$); $\Delta\gamma_{0}$ is degenerate
with external shear; and so, only $F_{0}$ and $G_{0}$ can produce useful
modeling constraints. However, since the axisymmetric NFW model contains three
free parameters $\tilde{\kappa},\,h,\,x=\theta_{\rm s}/h$, even a perfect
measurement of the $F_{0}$ and $G_{0}$ terms still leaves a degeneracy.

###  A.1 MSD

Consider the usual MSD transformation, induced by a parameter $\lambda$:

$\displaystyle\Delta\kappa$ $\displaystyle\to$
$\displaystyle\Delta\kappa_{\lambda}=\lambda\Delta\kappa+1-\lambda,$ (38)
$\displaystyle F_{0}$ $\displaystyle\to$ $\displaystyle F_{0,\lambda}=\lambda
F_{0},\;\;\;G_{0}\,\to\,G_{0,\lambda}=\lambda G_{0}$ (39)

(with matching transformation on the other terms, that are not relevant here).
The combinations $F_{\rm eff}$ and $G_{\rm eff}$ of Eq. (14) are MSD-
invariant, and are the quantities that constrain the model parameters. Suppose
then that we are given precise determination of $F_{\rm eff}$ and $G_{\rm
eff}$ from the imaging. In this case, the ratio $F_{\rm eff}/G_{\rm eff}$
determines $x$ . Having fixed $x$, there remains a degeneracy in
$\Delta\kappa$, that we can express as

$\displaystyle\Delta\kappa$ $\displaystyle=$ $\displaystyle\frac{h\,G_{\rm
eff}\,f_{\kappa}(x)}{f_{G}(x)+h\,G_{\rm eff}\,f_{\kappa}(x)}.$ (40)

In Eq. (40), we think of $G_{\rm eff}$ and $x$ as fixed by the imaging data,
while the model parameter $h$ is free to vary (up to possible external priors,
discussed in the main text).

For $|h\,G_{\rm eff}|\ll 1$, the parameter regime of most interest for us, we
have $\Delta\kappa\approx h\,G_{\rm eff}\,f_{\kappa}(x)/f_{G}(x)$. Namely, an
imaging determination of the flexion terms $G_{\rm eff}$ and $F_{\rm eff}$
cannot determine $\Delta\kappa$, which is almost directly degenerate with a
change in the model parameter $h$ while holding $x=\theta_{\rm s}/h$ and
$\tilde{\kappa}/h$ fixed.

Figure 6: $x$-dependence of lensing expansion terms.

In terms of the original NFW model parameters, the $\Delta\kappa$ degeneracy
maps to adjusting $R_{\rm s}\propto h$ while holding $\rho_{0}$ fixed.
External priors are needed to break this degeneracy. A prior on $h$, the
cluster center of mass position, obviously ameliorates it. A prior on the
group velocity dispersion also ameliorates it, since $\sigma_{\rm
los}^{2}\propto\rho_{0}R_{\rm s}^{2}$. The important point is that the
external priors are crucial: without them, the MSD associated with an NFW
group is not broken by the imaging data, even with explicit modeling of the
group.

###  A.2 $G$-term degeneracy

The threefold $G$-term degeneracy can also be clarified using the analytic
expansion. To this end, it is useful to replace the model parameter
$\tilde{\kappa}$ by $\Delta\kappa$. As we have seen in Sec. IV, an MCMC
trapped near a displaced likelihood minimum will attempt to reduce $|F_{\rm
eff}|$ while keeping $G_{\rm eff}$ fixed. This amounts to flowing towards $x$
that minimizes $\left|F_{\rm eff}/G_{\rm eff}\right|=|f_{F}(x)/f_{G}(x)|$,
while adjusting $\Delta\kappa$ and $h$ so as to keep $G_{\rm
eff}=(f_{G}/f_{\kappa})(\Delta\kappa/h)/(1-\Delta\kappa)$=Const. As can be
seen from Fig. 6, minimizing $|f_{F}(x)/f_{G}(x)|$ pulls the fit towards
smaller $x$. In turn, this pulls $|f_{G}(x)/f_{\kappa}(x)|$ to a larger value,
meaning that in order to compensate and keep $G_{\rm eff}$ constant, the
combination $\left(\Delta\kappa/h\right)/\left(1-\Delta\kappa\right)$ is
pulled to a smaller value. For $\Delta\kappa\ll 1$, this tends to make the fit
pull towards a model in which $\Delta\kappa$ is smaller than its truth value,
causing a positive upwards bias in $H_{0}$.

###  A.3 Time delays

Time delays between image A and B can be written as (see e.g. Schneider  _et
al._ (1992); Bartelmann and Schneider (2001))

$\displaystyle\begin{aligned} \Delta t_{\rm AB}&=D_{\rm
dt}\quantity(\frac{\theta^{2}_{\rm A}}{2}-\vec{\beta}\cdot\vec{\theta}_{\rm
A}-\Psi(\vec{\theta}_{\rm A})-(\mathrm{A}\leftrightarrow\mathrm{B}))\\\
&=:D_{\rm dt}(\tau_{\rm A}-\tau_{\rm B})\ ,\ D_{\rm dt}=(1+z_{\rm
l})\frac{d_{\rm l}d_{\rm s}}{d_{\rm ls}}.\end{aligned}$ (41)

One can find the time delay from the complex lens equation, by noticing

$\gradient_{\rm
c}=2\partialderivative{\theta^{*}}\implies\tau=\frac{1}{2}\int\differential{\theta^{*}}\alpha+f(\theta)\
,$ (42)

where the integral is a definite integral. Notice that the integral misses a
possible function of $\theta$ alone; this function can be recovered by
imposing $\tau=\tau^{*}$.

Focusing on $\tau_{\rm A}$, in complex notation we can write

$\displaystyle\begin{aligned} \tau_{\rm A}&=\frac{\theta_{\rm
A}\theta^{*}_{\rm A}}{2}-\frac{1}{2}\quantity[(\beta+\gradient_{\rm
c}\Psi_{\rm nfw})\theta^{*}_{\rm A}+\mathrm{c.c.}]-\Psi_{\rm
lens}(\vec{\theta}_{\rm A})\\\ &-\Psi_{\rm
halo}(0,0)-\frac{1}{2}\kappa\theta_{\rm A}\theta^{*}_{\rm
A}-\frac{1}{4}\quantity[(\gamma+\gamma^{\rm ext}){\theta^{*}_{\rm
A}}^{2}+\mathrm{c.c.}]\\\ &-\frac{1}{24}\quantity(3F{\theta^{*}_{\rm
A}}^{2}\theta_{\rm A}+G{\theta^{*}_{\rm A}}^{3}+\mathrm{c.c.}).\end{aligned}$
(43)

In real notation,

$\displaystyle\begin{aligned} \tau_{\rm A}&=\frac{\theta^{2}_{\rm
A}}{2}-(\beta_{i}+\partial_{i}\Psi_{\rm nfw})\theta_{i}^{\rm A}-\Psi_{\rm
lens}(\vec{\theta}_{\rm A})-\Psi_{\rm halo}(0,0)\\\
&-\frac{1}{2}\kappa\theta^{2}_{\rm A}-\frac{1}{2}(\gamma_{1}+\gamma^{\rm
ext}_{1})({\theta^{\rm A}}^{2}_{1}-{\theta^{\rm
A}}^{2}_{2})-(\gamma_{2}+\gamma^{\rm ext}_{2})\theta^{\rm A}_{1}\theta^{\rm
A}_{2}\\\ &-\frac{1}{12}\Big{(}(3F_{1}+G_{1}){\theta^{\rm
A}}^{3}_{1}+(3F_{2}-G_{2}){\theta^{\rm A}}^{3}_{2}\\\
&+3(F_{1}-G_{1})\theta^{\rm A}_{1}{\theta^{\rm
A}}^{2}_{2}+3(F_{2}+G_{2})\theta^{\rm A}_{2}{\theta^{\rm
A}}^{2}_{1}\Big{)}.\end{aligned}$ (44)

###  A.4 Effect of ellipticity

We can implement ellipticity in the NFW profile by using pseudo-elliptical NFW
lens models Golse and Kneib (2002); Dumet-Montoya  _et al._ (2012), which is a
reliable approximation as long as the ellipticity is not too large. This
approximation, implemented in lenstronomy, uses the spherical NFW lensing
potential, computed in elliptical coordinates,

$\Psi_{\epsilon}(\vec{\theta})=\Psi(\vec{\theta}_{\epsilon})\ ,\
\vec{\theta}_{\epsilon}=\begin{pmatrix}\cos\varphi\sqrt{1-\epsilon}&\sin\varphi\sqrt{1-\epsilon}\\\
-\sin\varphi\sqrt{1+\epsilon}&\cos\varphi\sqrt{1+\epsilon}\\\
\end{pmatrix}\vec{\theta},$ (45)

where

$\epsilon=\frac{2e}{1+e^{2}},\ e=\sqrt{e_{1}^{2}+e_{2}^{2}},\
\varphi=\frac{1}{2}\arctan\frac{e_{1}}{e_{2}}.$ (46)

Notice that under rotations of the coordinate system, $\epsilon$ is invariant
but $\varphi$ changes. The expansion of the potential
$\Psi_{\mathrm{halo},\epsilon}$ follows Eq. (17); however, since the potential
is not axisymmetric, it is not possible to factor out the $\phi$ dependence on
$\Delta\gamma_{\epsilon}$, $G_{\epsilon}$, $F_{\epsilon}$ as we did in Eq.
(4). Nevertheless, it is possible to express these quantities as functions of
$\Delta\gamma_{0}$, $F_{0}$ and $G_{0}$.

Notice that we can write (here, $\vec{\theta}=(x,y)^{\top}$)

$\gradient_{\rm
c}=\mathrm{e}^{\mathrm{i}\varphi}\quantity(\sqrt{1-\epsilon}\partialderivative{x_{\epsilon}}+\mathrm{i}\sqrt{1+\epsilon}\partialderivative{y_{\epsilon}}).$
(47)

As an example, consider

$\displaystyle\gamma_{\epsilon}$ $\displaystyle=\frac{1}{2}\gradient_{\rm
c}\gradient_{\rm
c}\Psi_{\epsilon}=\frac{\mathrm{e}^{2\mathrm{i}\varphi}}{2}\Big{(}(1-\epsilon)\partialderivative{{}^{2}}{x^{2}_{\epsilon}}-(1+\epsilon)\partialderivative{{}^{2}}{y^{2}_{\epsilon}}$
$\displaystyle+2\mathrm{i}\sqrt{1-\epsilon^{2}}\partialderivative{x_{\epsilon}}\partialderivative{y_{\epsilon}}\Big{)}\Psi_{\epsilon}$
$\displaystyle=\mathrm{e}^{2\mathrm{i}\varphi}(\gamma_{1}(\vec{\theta}_{\epsilon})+\mathrm{i}\sqrt{1-\epsilon^{2}}\gamma_{2}(\vec{\theta}_{\epsilon})-\epsilon\kappa(\vec{\theta}_{\epsilon})),$
(48)

where the quantities without ϵ subscript refer to the spherical functions
computed at the elliptic coordinate. By defining the displacement vector in
elliptical coordinates in the complex notation, $\vec{h}\leftrightarrow
h_{\epsilon}\mathrm{e}^{\mathrm{i}\phi_{e}}$, we can express

$\gamma=\gamma_{0}(h_{\epsilon}(h,\phi))\mathrm{e}^{\mathrm{i}\phi_{\epsilon}(\phi)}.$
(49)

For the NFW example, $\gamma_{0}$ is the expression in Eq. (35). For $\kappa$,
one has

$\kappa_{\epsilon}=\kappa(h_{\epsilon}(h,\phi))-\epsilon\gamma_{1}(h_{\epsilon}(h,\phi)).$
(50)

A similar reasoning follows for $F_{\epsilon}$, $G_{\epsilon}$. As one would
expect, ellipticity modifies the effective flexion and convergence by
correction terms of order $\epsilon$. Our numerical MCMC analysis includes the
full effect, and moderate or small ellipticity has only a minor effect on the
results.

##  Appendix B Kinematics constraints for host group

Here we review constraints on a host group, obtainable by measurements of a
sample of galaxy members. We only consider spherical systems with an isotropic
velocity distribution. We start with a PL mass and tracer galaxy distribution,
and go on to consider the NFW profile. For clarity, our discussion repeats a
number of statements from the main text, completing these with details and
derivations.

###  B.1 Power law

We start with the PL profile,
$\rho(r)=\rho_{0}\left(r/R_{0}\right)^{-\gamma}$. The surface density of this
profile is

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\pi}\Gamma\left(\frac{\gamma-1}{2}\right)\rho_{0}R_{0}}{\Gamma\left(\frac{\gamma}{2}\right)}\left(\frac{d_{\rm
l}\theta}{R_{0}}\right)^{1-\gamma},$ (51)

resulting with the deflection angle and convergence,

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\left(\frac{\theta}{\theta_{\rm
E}}\right)^{1-\gamma}\theta,\,\,\,\,\kappa\,=\,\frac{3-\gamma}{2}\left(\frac{\theta}{\theta_{\rm
E}}\right)^{1-\gamma},$ (52)

where

$\displaystyle\theta_{\rm E}$ $\displaystyle=$
$\displaystyle\left[\frac{2\sqrt{\pi}\Gamma\left(\frac{\gamma-1}{2}\right)}{(3-\gamma)\Gamma\left(\frac{\gamma}{2}\right)}\right]^{\frac{1}{\gamma-1}}\frac{R_{0}}{d_{\rm
l}}\left(\frac{\rho_{0}R_{0}}{\Sigma_{c}}\right)^{\frac{1}{\gamma-1}}.$ (53)

The circular velocity is

$\displaystyle v_{\rm circ}^{2}$ $\displaystyle=$
$\displaystyle\frac{GM(r)}{r}\,=\,\frac{4\pi
G\rho_{0}R_{0}^{2}}{3-\gamma}\left(\frac{d_{\rm
l}\theta}{R_{0}}\right)^{2-\gamma}.$ (54)

The circular velocity is not directly measurable. What is measurable is the
LOSVD Binney and Tremaine (2008),

$\displaystyle\sigma_{\rm los}^{2}$ $\displaystyle=$
$\displaystyle\frac{2G}{s_{*}(r)}\int_{1}^{\infty}\frac{dy\,y}{\sqrt{y^{2}-1}}\,\int_{y}^{\infty}\frac{dx}{x^{2}}n_{*}(xr)M(xr)$
$\displaystyle=$ $\displaystyle\frac{8\pi
G\rho_{0}R_{0}^{3}r^{3-\gamma}}{(3-\gamma)s_{*}(r)}\int_{1}^{\infty}\frac{dy\,y}{\sqrt{y^{2}-1}}\,\int_{y}^{\infty}dx\,n_{*}(xr)\,x^{1-\gamma},$

where $n_{*}$ and $s_{*}$ are the galaxy number density and surface density,
respectively. Eq. (B.1) holds for an arbitrary isotropic profile, while Eq.
(LABEL:eq:losvdPL1) applies to the PL. If the galaxy number density is
distributed similarly to the mass density in the group,
$n_{*}(r)\propto(r/R_{0})^{-\gamma}$, we find (valid for
$\frac{3}{2}<\gamma<\frac{5}{2}$):

$\displaystyle\sigma_{\rm los}^{2}$ $\displaystyle=$
$\displaystyle\frac{\Gamma^{2}\left(\frac{\gamma}{2}\right)\Gamma\left(\gamma-\frac{3}{2}\right)}{4\sqrt{\pi}\,\Gamma^{2}\left(\frac{\gamma-1}{2}\right)\Gamma\left(\gamma\right)}\frac{d_{\rm
s}}{d_{\rm ls}}\theta_{\rm E}\left(\frac{\theta}{\theta_{\rm
E}}\right)^{2-\gamma}.$ (57)

For the SIS ($\gamma=2$), $\theta_{\rm E}=8\pi^{2}\frac{d_{\rm ls}}{d_{\rm
s}}G\rho_{0}R_{0}^{2}$ and $\sigma_{\rm los}^{2}=\frac{d_{\rm s}}{d_{\rm
ls}}\frac{\theta_{\rm E}}{4\pi}=\frac{v_{\rm circ}^{2}}{2}$.

The relation $\sigma_{\rm los}^{2}=\frac{d_{\rm s}}{d_{\rm
ls}}\frac{\theta_{\rm E}}{4\pi}$ is sometimes adopted by lensing analyses.
However, the expressions for $\sigma_{\rm los}^{2}$ depend on a series of
simplifying assumptions, including: PL mass distribution; same PL galaxy
number distribution; virial state (following Jeans equation); spherical
symmetry; isotropic velocity distribution; and to satisfy $\sigma_{\rm
los}^{2}=\frac{d_{\rm s}}{d_{\rm ls}}\frac{\theta_{\rm E}}{4\pi}$, the
specific case of the SIS. We do not expect those assumptions to hold precisely
for individual galaxy groups. Lensing analyses would be prudent to assign a
larger systematic uncertainty than the “bare” observational uncertainty on
$\sigma^{2}_{\rm los}$.

To explore one of those effects, consider varying $\gamma$ in the relation
between $\sigma_{\rm los}^{2}$ and $\theta_{\rm E}$. For non-SIS systems,
$\sigma_{\rm los}^{2}(\theta)$ depends on $\theta$, and a relevant observable
is the brightness-weighted average of $\sigma_{\rm los}^{2}$ in some aperture
$\theta_{A}$. We can compare such averaged LOSVD to the SIS relation:

$\displaystyle\frac{\langle\sigma_{\rm los}^{2}\rangle_{A}}{\frac{d_{\rm
s}}{d_{\rm ls}}\frac{\theta_{\rm E}}{4\pi}}$ $\displaystyle=$
$\displaystyle\frac{1}{\frac{d_{\rm s}}{d_{\rm ls}}\frac{\theta_{\rm
E}}{4\pi}}\frac{\int_{0}^{\theta_{A}}d\theta\theta
s_{*}(\theta)\sigma^{2}_{\rm los}(\theta)}{\int_{0}^{\theta_{A}}d\theta\theta
s_{*}(\theta)}$ (58) $\displaystyle=$
$\displaystyle\frac{3-\gamma}{5-2\gamma}\frac{\sqrt{\pi}\Gamma^{2}\left(\frac{\gamma}{2}\right)\Gamma\left(\gamma-\frac{3}{2}\right)}{\Gamma^{2}\left(\frac{\gamma-1}{2}\right)\Gamma\left(\gamma\right)}\left(\frac{\theta_{\rm
A}}{\theta_{\rm E}}\right)^{2-\gamma}.\;\;\;\;\;$ (59)

We illustrate Eq. (59) in Fig. 7.

Figure 7: $\gamma$\- and $\theta_{A}$\- dependence of aperture-averaged LOSVD.

###  B.2 NFW

For the NFW profile, assuming again that the tracer galaxy density follows the
mass density, Eq. (B.1) can be written as:

$\displaystyle\sigma_{\rm los}^{2}(\theta)$ $\displaystyle=$ $\displaystyle
G\rho_{0}R_{\rm s}^{2}\,f\left(\frac{\theta}{\theta_{\rm s}}\right),$ (60)
$\displaystyle f(a)$ $\displaystyle=$ $\displaystyle\frac{8\pi
a\int_{1}^{\infty}\frac{dyy}{\sqrt{y^{2}-1}}\int_{ya}^{\infty}\frac{dx}{x^{3}\left(1+x\right)^{2}}\left(\ln\left(1+x\right)-\frac{x}{1+x}\right)}{2\int_{0}^{\infty}\frac{dz}{\sqrt{a^{2}+z^{2}}(1+\sqrt{a^{2}+z^{2}})^{2}}}.$

Some of these integrals can be done in closed form, but that is not
particularly illuminating. The function $f(\theta/\theta_{\rm s})$ is shown in
Fig. 8.

Figure 8: The LOSVD $\theta$-dependence for the NFW profile.

Assuming the NFW model, a kinematics prior should use Eq. (60) (aperture-
averaged as needed) to translate a LOSVD measurement into a constraint on the
combination of $\rho_{0}$ and $R_{\rm s}$ appearing in the equation. However,
analyses in the literature took a somewhat different route. Ref. Chen  _et
al._ (2019) (after Koranyi and Geller (2000); Wong  _et al._ (2011))
considered the following expression as a proxy for the LOSVD of the NFW
profile:

$\displaystyle\bar{\sigma}^{2}$ $\displaystyle=$ $\displaystyle\frac{GM_{\rm
vir}}{3R_{\rm vir}},$ (62)

where $R_{\rm vir}=c_{\rm vir}\,R_{s}$, $M_{\rm vir}=M(R_{\rm vir})$, and
$c_{\rm vir}$ is the NFW concentration parameter, for which it is possible to
extract a theoretical prediction from N-body simulations Maccio’  _et al._
(2008). Namely, instead of using the physical relation, Eq. (60), to convert
the measured LOSVD into a constraint on the group halo model, Ref. Chen  _et
al._ (2019) made the identification $\sigma_{\rm los}^{2}\to\bar{\sigma}^{2}$,
and then used Eq. (62) to constrain $M_{\rm vir}$ and $R_{\rm vir}$. The
dependence of the RHS of Eq. (62) on $c_{\rm vir}$ can be clarified by noting
that

$\displaystyle\bar{\sigma}^{2}$ $\displaystyle=$ $\displaystyle
G\rho_{0}R_{s}^{2}\frac{4\pi}{3c_{\rm vir}}\left(\ln(1+c_{\rm
vir})-\frac{c_{\rm vir}}{1+c_{\rm vir}}\right).$ (63)

The dependence is not strong: the factor $\frac{4\pi}{3c_{\rm
vir}}\left(\ln(1+c_{\rm vir})-\frac{c_{\rm vir}}{1+c_{\rm vir}}\right)$ is
equal to $\left\\{0.89,0.85,0.80,0.62\right\\}$ when varying $c_{\rm
vir}=\left\\{3,4,5,10\right\\}$, respectively. However, it is noteworthy that
this dependence is not physical, but rather introduced artificially. The red
solid line in Fig. 9 compares the physical aperture-weighted LOSVD to the
quantity $\bar{\sigma}^{2}$, using $c_{\rm vir}=3.5$, relevant for low-
redshift massive galaxies.

Figure 9: Aperture-averaged LOSVD for NFW, compared to the quantity
$\bar{\sigma}^{2}=GM_{\rm vir}/(3R_{\rm vir})$. Red line: using Eq. (60) in
Eq. (58), replacing $\frac{d_{\rm s}}{d_{\rm ls}}\frac{\theta_{\rm E}}{4\pi}$
by $\bar{\sigma}^{2}$, with $c_{\rm vir}=3.5$. Red shaded band: varying
$c_{\rm vir}\to 3.5\times 10^{\pm 0.14}$. Dots: direct calculation of
$\sigma_{\rm los}^{2}$ from mock groups, assuming 13 tracer galaxies per
realization (see Sec. B.4). Blue line (somewhat hidden under the red) and blue
shaded band: mean and $1\sigma$ region for the mocks.

The meaning of Fig. 9 is that using Eq. (62) and identifying $\sigma^{2}_{\rm
los}\to\bar{\sigma}^{2}$ can introduce a systematic error. The magnitude of
the error depends on the aperture and on the value adopted for $c_{\rm vir}$,
but can amount to $\sim 20\%$.

###  B.3 Kinematical and cosmological priors

Here we describe our attempt to define a kinematics+cosmology prior, following
a similar procedure as used in Chen  _et al._ (2019) for PG1115+080.

The first step taken by Chen  _et al._ (2019) is to identify $\sigma_{\rm
los}\to\bar{\sigma}$.

The next step is to invoke theoretical expectations based on N-body
simulations. Following App. A of Wong  _et al._ (2011), the virial mass and
the virial radius of the halo are related to a characteristic redshift-
dependent overdensity:

$\displaystyle M_{\rm vir}$ $\displaystyle=$
$\displaystyle\frac{4\pi}{3}R^{3}_{\rm vir}\Delta_{\rm c}(z)\rho_{\rm c}(z),$
(64)

where $\Delta_{\rm c}(z)=178\quantity(\rho_{\rm m}(z)/\rho_{\rm c}(z))^{0.45}$
is the expected halo overdensity at the virial radius and $\rho_{\rm c}(z)$
and $\rho_{\rm m}(z)$ are the cosmological critical density and matter density
at redshift $z$ Eke  _et al._ (2001). Combining Eqs. (64) and (62) turns
$\bar{\sigma}$ into separate priors on $R_{\rm vir}$ and $M_{\rm vir}$:

$\displaystyle R_{\rm vir}$ $\displaystyle=$
$\displaystyle\frac{3\bar{\sigma}}{2\sqrt{\pi G\Delta_{\rm c}\rho_{\rm c}}},$
(65) $\displaystyle M_{\rm vir}$ $\displaystyle=$
$\displaystyle\frac{9\bar{\sigma}^{3}}{2\sqrt{\pi G^{3}\Delta_{\rm c}\rho_{\rm
c}}}.$ (66)

Next, Ref. Chen  _et al._ (2019) adds a prior on $c_{\rm vir}$, using the
$M_{\rm vir}-c_{\rm vir}$ relation from Ref. Maccio’  _et al._ (2008):

$\displaystyle\log_{10}c_{\rm vir}$ $\displaystyle\approx$ $\displaystyle
0.97-0.09\log_{10}\quantity(\frac{M_{\rm vir}\,h}{10^{12}M_{\odot}}).$ (67)

Treating $c_{\rm vir}$ as a function of $\bar{\sigma}$ via Eqs. (67) and (66),
one has now produced separate priors on the parameters $\rho_{0}$ and $R_{\rm
s}$:

$\displaystyle\rho_{0}$ $\displaystyle=$ $\displaystyle\frac{\Delta_{\rm
c}\rho_{\rm c}c^{3}_{\rm vir}}{3}\quantity(\ln(1+c_{\rm vir})-\frac{c_{\rm
vir}}{1+c_{\rm vir}})^{-1},$ (68) $\displaystyle R_{\rm s}$ $\displaystyle=$
$\displaystyle\frac{R_{\rm vir}}{c_{\rm vir}}.$ (69)

The uncertainty can be estimated by combining the observational uncertainty on
$\sigma_{\rm los}$ and the theoretical scatter on the $M_{\rm vir}-c_{\rm
vir}$ relation; the latter is taken in Ref. Chen  _et al._ (2019) as $\pm
0.14$ on the $\log_{10}c_{\rm vir}$ expression in Eq. (67).

In practice, our implementation of these priors in the MCMC mock analysis is
as follows:

  1. 1.

The parameters that are directly sampled by the chain, defining the posterior
likelihood space, are angular variables that are linearly proportional to
$\rho_{0}$ and $R_{\rm s}$. In each step of the calculation, we translate the
sampled value of $\rho_{0}$ into an expected value of $c_{\rm vir}$, using Eq.
(68).

  2. 2.

Having obtained $c_{\rm vir}$, we convert Eq. (67) into an equation for
$R_{\rm vir}$, replacing $M_{\rm vir}$ by $R_{\rm vir}$ using Eq. (64). Now,
we implement the estimated theoretical scatter of $\pm 0.14$ dex on $c_{\rm
vir}$ in the $M_{\rm vir}-c_{\rm vir}$ relation, to define a range of
acceptable values for $R_{\rm s}$ via $R_{\rm s}\subset R_{\rm
vir}/(10^{+0.14},10^{-0.14})$. In the chain, we discard sampled values of
$R_{\rm s}$ that fall outside of this range. This step therefore enforces a
redshift-dependent correlation between $\rho_{0}$ and $R_{\rm s}$. Up to this
point, we made no connection to kinematics.

  3. 3.

Finally we come to the kinematics. Using Eq. (63), we convert the model point
represented by $\rho_{0}$ and $R_{\rm s}$, along with the central value of
$c_{\rm vir}$, into a model prediction for $\bar{\sigma}$. This prediction is
then compared to the measured value of $\sigma_{\rm los}$ using the nominal
observational uncertainty to obtain a likelihood factor.

###  B.4 Sample variance

Sample variance is the dominant nominal source of uncertainty in
$\sigma^{2}_{\rm los}$, given that typically only a handful of galaxies are
measured as tracers of the group. Ref. Wilson  _et al._ (2017) used bootstrap
to estimate sample variance directly from the measured sample of galaxies;
here, we complement this route by generating random sets of $N_{\rm g}$
galaxies tracing an NFW halo.

The realizations are drawn from an equilibrium phase space distribution
function $f(r,v)$. For a spherical halo with a statistically static
distribution function, we have $f(r,v)=f(\varepsilon(r,v))$, where
$\varepsilon=\psi(r)-\frac{v^{2}}{2}$, $\psi(r)=-\phi(r)$, and $\phi(r)$ is
the halo Newtonian potential. $f(\varepsilon)$ is given by999See Binney and
Tremaine (2008), Ch.4.3.1.

$\displaystyle f(\varepsilon)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{8}\pi^{2}}\int_{0}^{\varepsilon}\frac{d\psi}{\sqrt{\varepsilon-\psi}}\frac{d^{2}\rho}{d\psi^{2}}.$
(70)

We calculate $f(\varepsilon)$ numerically, and use it to draw samples of
$N_{\rm g}$ tracer galaxies that fall within projected aperture $\theta_{\rm
A}$. The results of this exercise for $N_{\rm g}=13$ (as in PG1115+080) and
for different apertures are shown as blue dots in Fig. 9. For each value of
$\theta_{\rm A}/\theta_{\rm s}$ we generate 100 mock samples. For each sample
we calculate $\sigma_{\rm los}^{2}$ directly as the variance of LOS velocity
across the $N_{\rm g}$ galaxies. The mean and standard deviation of
$\sigma_{\rm los}^{2}$ are shown by the thick blue line and shaded region.

The variance we find for $\sigma^{2}_{\rm los}$ in Fig. 9 is roughly
consistent with sample variance of a normal distribution:
$\frac{\delta\sigma^{2}_{\rm los}}{\sigma^{2}_{\rm
los}}\approx\sqrt{\frac{2}{N_{g}-1}}\approx 0.4$ for $N_{g}=13$. For
comparison with Ref. Wilson  _et al._ (2016), we also calculate uncertainty
estimates for $\sigma_{\rm los}^{2}$ using the bootstrap method. As a rule,
the bootstrap method provides a slightly lower uncertainty estimate than the
variance found with mock realizations, but the difference is small: whereas
direct sample variance predicts $\Delta\sigma_{\rm los}^{2}/\sigma_{\rm
los}^{2}\approx 0.4$, bootstrap predicts $\Delta\sigma_{\rm
los}^{2}/\sigma_{\rm los}^{2}\approx 0.35$. Thus, we reproduce the sample
variance-dominated uncertainty estimate of Ref. Wilson  _et al._ (2017) for
PG1115+080.

##  Appendix C Estimating the probability of an MCMC to fall into a displaced
minimum

As noted in Sec. III.3, at least for a wide group centroid prior (as derived
in Ref. Wilson  _et al._ (2017) for PG1115+080), trial and error with emcee
Foreman-Mackey  _et al._ (2013) suggests that the initial placement of the
walkers is a key factor in deciding which likelihood minimum will attract the
fit. We can (roughly) estimate the probability of falling into a wrong minimum
by the probability for the walker placement to start off closer to a false
minimum than to the truth one. Let us assume that the initial placement of the
walkers is chosen to coincide with the prior’s center (this seems like a
natural choice). In this case, the probability to fall into a wrong minimum is
roughly given by the probability of the group centroid prior to be nearer a
false minimum than the truth one.

We can estimate this probability using mock samples of tracer galaxies as in
Sec. B.4. Consider a sample of $N_{\rm g}$ galaxies, and choose a random
member to be the “primary lens”. From the same sample, derive a group center
prior as the center of mass of the members. In general, of course, the prior
center does not coincide with the center of the halo used to generate the
mock.101010The statistical distribution of the mismatch, which goes to
defining the prior width, can be estimated either by bootstrap, as done in
Wilson  _et al._ (2017) for the actual data, or by repeated mocks. We test
both, and find them to be compatible. Now, draw a line connecting the primary
lens with the true halo center, and another line connecting the primary lens
with the prior center: the prior center is closer to a false minimum if the
smaller angle between these lines is larger than $60^{\rm o}$.

The result of this calculation depends on the number of member galaxies
$N_{\rm g}$, the group halo profile, and the analysis aperture. For $N_{\rm
g}=13$ in the NFW model with an aperture of 10 $R_{\rm s}$ (around 3 virial
radii), we find the wrong minimum probability to be $\sim 10$%.

It should be clear that the simple estimate we described here ignores a range
of possible selection effects, both natural (such as mass segregation, and the
bias of more massive group members to become primary lenses) as well as
analysis-specific (such as luminosity-dependent contamination and
incompleteness). Such effects can probably modify our estimates at the
$\mathcal{O}(1)$ level, but should not change the order of magnitude of the
result.

##  Appendix D How common is the flexion degeneracy?

PG1115+080 is proof of concept that lensing set-ups like the one we analyzed
are observationally relevant. But since PG1115+080 was selected specifically
because of its massive, near-by group association, this set-up may be quite
uncommon. Most of the lensed quasars analyzed for time-delay cosmography may
either not be associated with LOS groups; or, if they are, may lie far from
the group’s projected center of mass, making flexion terms less important than
for PG1115+080. A detailed analysis of the fraction of systems that may
exhibit flexion degeneracy is beyond the scope of this work. Nevertheless we
can make a crude estimate using the data from Wilson  _et al._ (2017), as
follows.

Tab. 1 in Ref. Wilson  _et al._ (2017) lists reconstructed properties of LOS
groups (including primary lens host where exists) of their 26 lens systems.
For each LOS group, estimated values of $M_{\rm vir}$ and $R_{\rm vir}$ are
quoted, along with the measured LOSVD, denoted there by $\sigma_{\rm grp}$. We
insert $M_{\rm vir}$ into Eq. (67) to estimate $c_{\rm vir}$; then use
$\sigma_{\rm grp}$, $R_{\rm vir}$, and $c_{\rm vir}$ in Eq. (63) to estimate
$\rho_{0}$. Using the quoted group redshift to calculate the relevant value of
$\Sigma_{\rm c}$, Eqs. (32) and (34) allow us to obtain $\Delta\kappa$. For
each system, we estimate the rough magnitude of the flexion by the root-mean-
square $\sqrt{\sum\limits_{i}\left(\Delta\kappa_{i}/h_{i}\right)^{2}}$,
including all LOS groups $i$, and estimate the corresponding deflection angle
near the primary lens by
$\Delta\theta\sim\sqrt{\sum\limits_{i}\left(\Delta\kappa_{i}/h_{i}\right)^{2}}\theta_{\rm
E}^{2}$, taking $\theta_{\rm E}=1^{\prime\prime}$ for definiteness.

Of the 26 systems considered in Wilson  _et al._ (2017), PG1115+080 does
indeed come up with the highest flexion estimate of $\Delta\theta\sim
0.02^{\prime\prime}$. The next highest system is HE0435, with
$\Delta\theta\sim 0.005^{\prime\prime}$. Altogether, 3 systems (7 systems) out
of the 26 have $\Delta\theta\gtrsim 0.0025^{\prime\prime}$
($\Delta\theta\gtrsim 1$ mas).

Numerical experiments with mock analyses mimicking the main features of the
pipeline of Millon  _et al._ (2020) suggest that flexion terms are
quantitatively important in the fit if the flexion-induced deflection angle is
larger than $\sim 1$ mas. Very roughly, we can expect that flexion degeneracy
should be a concern in the same parametric regime. It is interesting to note
that all of the 4 systems that take part in both of the cosmography Millon
_et al._ (2020) and kinematics Wilson  _et al._ (2017) campaigns, turn out to
exhibit $\Delta\theta\gtrsim 1$ mas.

##  Appendix E Full corner plots.

In this Appendix we collect some detailed results from the MCMC analysis.

In Fig. 10 we show triangle plots in which the cosmological prior on $c_{\rm
vir}$ is not included. We do this exercise in order to investigate the impact
of this prior on the results. The main point to notice is that omitting the
$c_{\rm vir}$ prior, the bias on $H_{0}$ becomes somewhat more pronounced
(compare Fig. 5, that includes this prior). At the same time, without this
prior, the best fit result for $R_{\rm s}$ in displaced (false) posterior
likelihood minima is driven to small values. This point is shown by a
comprehensive triangle plot in Fig. 11.

Figure 10: Mock analysis: zeus run (upper-right) and emcee runs, all with a
kinematics prior but no $c_{\rm vir}$ prior. We note that in false minima
fits, the inferred $\theta_{\rm s}$ (not shown here) is driven to very low
values . Figure 11: emcee run, multi-parameter view, corresponding to the run
in the bottom-left panel of Fig. 10. The absence of a $c_{\rm vir}$ prior
causes the inference to pull towards low values of $R_{s}$ near a displaced
minimum. The parameter $\tilde{\rho}_{0}$ is defined as
$\tilde{\kappa}/\theta_{\rm s}$ (which is independent on $R_{\rm s}$) .

In Fig. 12 we show triangle plots in which the kinematics prior is enforced
with a standard deviation of $60$ km/s on $\sigma_{\rm los}$. These results
can be compared to Fig. 5 from the main text, where, as noted in Sec. III.2,
the standard deviation on $\sigma_{\rm los}$ was taken as $120$ km/s. We do
not find a significant difference.

Figure 12: Mock analysis: emcee runs using a standard deviation of $60$ km/s
on $\sigma_{\rm los}$. These results can be compared with Fig. 5 from the main
text, in which the uncertainty in $\sigma_{\rm los}$ was doubled .

Fig. 13 gives a more complete perspective on the degeneracies and the global
structure of the likelihood as exposed by a zeus run.

Figure 13: Mock analysis: zeus run, comprehensive view . A cautionary comment
is needed: the run has not converged. The main point of this run is to
qualitatively show the three minima structure. No quantitative results are
obtained from this run.

## References

  * Refsdal (1964) S. Refsdal, Monthly Notices of the Royal Astronomical Society 128, 307 (1964), https://academic.oup.com/mnras/article-pdf/128/4/307/8073517/mnras128-0307.pdf .
  * Suyu  _et al._ (2013) S. H. Suyu _et al._ , Astrophys. J. 766, 70 (2013), arXiv:1208.6010 [astro-ph.CO] .
  * Treu and Marshall (2016) T. Treu and P. J. Marshall, Astron. Astrophys. Rev. 24, 11 (2016), arXiv:1605.05333 [astro-ph.CO] .
  * Suyu  _et al._ (2017) S. H. Suyu _et al._ , Mon. Not. Roy. Astron. Soc. 468, 2590 (2017), arXiv:1607.00017 [astro-ph.CO] .
  * Grillo  _et al._ (2020) C. Grillo, P. Rosati, S. H. Suyu, G. B. Caminha, A. Mercurio, and A. Halkola, Astrophys. J. 898, 87 (2020), arXiv:2001.02232 [astro-ph.CO] .
  * Rusu  _et al._ (2020) C. E. Rusu _et al._ , Mon. Not. Roy. Astron. Soc. 498, 1440 (2020), arXiv:1905.09338 [astro-ph.CO] .
  * Birrer  _et al._ (2019) S. Birrer _et al._ , Mon. Not. Roy. Astron. Soc. 484, 4726 (2019), arXiv:1809.01274 [astro-ph.CO] .
  * Shajib  _et al._ (2020) A. J. Shajib _et al._ (DES), Mon. Not. Roy. Astron. Soc. 494, 6072 (2020), arXiv:1910.06306 [astro-ph.CO] .
  * Chen  _et al._ (2019) G. C. F. Chen _et al._ , Mon. Not. Roy. Astron. Soc. 490, 1743 (2019), arXiv:1907.02533 [astro-ph.CO] .
  * Wong  _et al._ (2020) K. C. Wong _et al._ , Mon. Not. Roy. Astron. Soc. 498, 1420 (2020), arXiv:1907.04869 [astro-ph.CO] .
  * Millon  _et al._ (2020) M. Millon _et al._ , Astron. Astrophys. 639, A101 (2020), arXiv:1912.08027 [astro-ph.CO] .
  * Riess  _et al._ (2019) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Astrophys. J. 876, 85 (2019), arXiv:1903.07603 [astro-ph.CO] .
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A. G. Riess, Nature Astron. 3, 891 (2019), arXiv:1907.10625 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021), arXiv:2103.01183 [astro-ph.CO] .
  * Falco  _et al._ (1985) E. E. Falco, M. V. Gorenstein, and I. I. Shapiro, Astrophys. J. Lett. 289, L1 (1985).
  * Kochanek (1991) C. S. Kochanek, Astrophys. J. 373, 354 (1991).
  * Kochanek (2002) C. S. Kochanek, Astrophys. J. 578, 25 (2002), arXiv:astro-ph/0205319 .
  * Liesenborgs and De Rijcke (2012) J. Liesenborgs and S. De Rijcke, Mon. Not. Roy. Astron. Soc. 425, 1772 (2012), arXiv:1207.4692 [astro-ph.CO] .
  * Schneider, Peter and Sluse, Dominique (2013) Schneider, Peter and Sluse, Dominique, A&A 559, A37 (2013).
  * Kochanek (2021) C. S. Kochanek, Mon. Not. Roy. Astron. Soc. 501, 5021 (2021), arXiv:2003.08395 [astro-ph.CO] .
  * Fleury  _et al._ (2021) P. Fleury, J. Larena, and J.-P. Uzan, JCAP 08, 024 (2021), arXiv:2104.08883 [astro-ph.CO] .
  * Teodori  _et al._ (2022) L. Teodori, K. Blum, E. Castorina, M. Simonović, and Y. Soreq, JCAP 07, 027 (2022), arXiv:2201.05111 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A1 (2020), arXiv:1807.06205 [astro-ph.CO] .
  * D’Amico  _et al._ (2020) G. D’Amico, J. Gleyzes, N. Kokron, K. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, JCAP 05, 005 (2020), arXiv:1909.05271 [astro-ph.CO] .
  * Ivanov  _et al._ (2020) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, JCAP 05, 042 (2020), arXiv:1909.05277 [astro-ph.CO] .
  * Blum  _et al._ (2020) K. Blum, E. Castorina, and M. Simonović, Astrophys. J. Lett. 892, L27 (2020), arXiv:2001.07182 [astro-ph.CO] .
  * Birrer  _et al._ (2020) S. Birrer _et al._ , Astron. Astrophys. 643, A165 (2020), arXiv:2007.02941 [astro-ph.CO] .
  * Blum and Teodori (2021) K. Blum and L. Teodori, Phys. Rev. D 104, 123011 (2021), arXiv:2105.10873 [astro-ph.CO] .
  * Wilson  _et al._ (2016) M. L. Wilson, A. I. Zabludoff, S. M. Ammons, I. G. Momcheva, K. A. Williams, and C. R. Keeton, Astrophys. J. 833, 194 (2016), arXiv:1710.09908 [astro-ph.GA] .
  * Wilson  _et al._ (2017) M. L. Wilson, A. I. Zabludoff, C. R. Keeton, K. C. Wong, K. A. Williams, K. D. French, and I. G. Momcheva, Astrophys. J. 850, 94 (2017), arXiv:1710.09900 [astro-ph.GA] .
  * Buckley-Geer  _et al._ (2020) E. J. Buckley-Geer _et al._ (DES), Mon. Not. Roy. Astron. Soc. 498, 3241 (2020), arXiv:2003.12117 [astro-ph.GA] .
  * Goldberg and Bacon (2005) D. M. Goldberg and D. J. Bacon, Astrophys. J. 619, 741 (2005), arXiv:astro-ph/0406376 .
  * Bacon  _et al._ (2006) D. J. Bacon, D. M. Goldberg, B. T. P. Rowe, and A. N. Taylor, Mon. Not. Roy. Astron. Soc. 365, 414 (2006), arXiv:astro-ph/0504478 .
  * Navarro  _et al._ (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107 [astro-ph] .
  * Yıldırım  _et al._ (2021) A. Yıldırım, S. H. Suyu, G. C. F. Chen, and E. Komatsu, (2021), arXiv:2109.14615 [astro-ph.CO] .
  * Shajib  _et al._ (2023) A. J. Shajib _et al._ , (2023), arXiv:2301.02656 [astro-ph.CO] .
  * Schramm and Kayser (1995) T. Schramm and R. Kayser, Astron. Astrophys. 299, 1 (1995), arXiv:astro-ph/9408064 .
  * Schneider and Er (2008) P. Schneider and X. Er, Astron. Astrophys. 485, 363 (2008), arXiv:0709.1003 [astro-ph] .
  * Momcheva  _et al._ (2006) I. Momcheva, K. A. Williams, C. R. Keeton, and A. I. Zabludoff, Astrophys. J. 641, 169 (2006), arXiv:astro-ph/0511594 .
  * Birrer and Amara (2018) S. Birrer and A. Amara, Physics of the Dark Universe 22, 189 (2018).
  * Birrer  _et al._ (2015) S. Birrer, A. Amara, and A. Refregier, Astrophys. J. 813, 102 (2015), arXiv:1504.07629 [astro-ph.CO] .
  * Birrer  _et al._ (2021) S. Birrer _et al._ , J. Open Source Softw. 6, 3283 (2021), arXiv:2106.05976 [astro-ph.CO] .
  * Jing and Suto (2002) Y. P. Jing and Y. Suto, Astrophys. J. 574, 538 (2002), arXiv:astro-ph/0202064 .
  * Vega-ferrero  _et al._ (2017) J. Vega-ferrero, G. Yepes, and S. Gottlöber, Mon. Not. Roy. Astron. Soc. 467, 3226 (2017), arXiv:1603.02256 [astro-ph.CO] .
  * Schrabback  _et al._ (2021) T. Schrabback _et al._ , Astron. Astrophys. 646, A73 (2021), arXiv:2010.00311 [astro-ph.CO] .
  * Grant  _et al._ (2004) C. E. Grant, M. W. Bautz, G. Chartas, and G. P. Garmire, Astrophys. J. 610, 686 (2004), arXiv:astro-ph/0305137 .
  * Wenger  _et al._ (2000) M. Wenger _et al._ , Astron. Astrophys. Suppl. Ser. 143, 9 (2000), arXiv:astro-ph/0002110 .
  * Fassnacht  _et al._ (2008) C. D. Fassnacht, D. D. Kocevski, M. W. Auger, L. M. Lubin, J. L. Neureuther, T. E. Jeltema, J. S. Mulchaey, and J. P. McKean, Astrophys. J. 681, 1017 (2008), arXiv:0711.2066 [astro-ph] .
  * Koranyi and Geller (2000) D. M. Koranyi and M. J. Geller, Astron. J. 119, 44 (2000), arXiv:astro-ph/9909514 .
  * Wong  _et al._ (2011) K. C. Wong, C. R. Keeton, K. A. Williams, I. G. Momcheva, and A. I. Zabludoff, Astrophys. J. 726, 84 (2011), arXiv:1011.2504 [astro-ph.CO] .
  * Eke  _et al._ (2001) V. R. Eke, J. F. Navarro, and M. Steinmetz, Astrophys. J. 554, 114 (2001), arXiv:astro-ph/0012337 .
  * Maccio’  _et al._ (2008) A. V. Maccio’, A. A. Dutton, and F. C. v. d. Bosch, Mon. Not. Roy. Astron. Soc. 391, 1940 (2008), arXiv:0805.1926 [astro-ph] .
  * Karamanis and Beutler (2021) M. Karamanis and F. Beutler, Stat. Comput. 31, 61 (2021), arXiv:2002.06212 [stat.ML] .
  * Karamanis  _et al._ (2021) M. Karamanis, F. Beutler, and J. A. Peacock, Mon. Not. Roy. Astron. Soc. 508, 3589 (2021), arXiv:2105.03468 [astro-ph.IM] .
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013).
  * Lovisari  _et al._ (2021) L. Lovisari, S. Ettori, M. Gaspari, and P. A. Giles,  (2021), 10.3390/universe7050139, arXiv:2106.13256 [astro-ph.CO] .
  * Oguri and Marshall (2010) M. Oguri and P. J. Marshall, Mon. Not. Roy. Astron. Soc. 405, 2579 (2010), arXiv:1001.2037 [astro-ph.CO] .
  * Liao  _et al._ (2015) K. Liao _et al._ , Astrophys. J. 800, 11 (2015), arXiv:1409.1254 [astro-ph.IM] .
  * Dobke  _et al._ (2009) B. M. Dobke, L. J. King, C. D. Fassnacht, and M. W. Auger, Mon. Not. Roy. Astron. Soc. 397, 311 (2009), arXiv:0904.1437 [astro-ph.CO] .
  * Birrer and Treu (2021) S. Birrer and T. Treu, Astron. Astrophys. 649, A61 (2021), arXiv:2008.06157 [astro-ph.CO] .
  * Treu  _et al._ (2018) T. Treu _et al._ (DES), Mon. Not. Roy. Astron. Soc. 481, 1041 (2018), arXiv:1808.04838 [astro-ph.CO] .
  * Bacon  _et al._ (2020) D. J. Bacon _et al._ (SKA), Publ. Astron. Soc. Austral. 37, e007 (2020), arXiv:1811.02743 [astro-ph.CO] .
  * Blanchard  _et al._ (2020) A. Blanchard _et al._ (Euclid), Astron. Astrophys. 642, A191 (2020), arXiv:1910.09273 [astro-ph.CO] .
  * Abell  _et al._ (2009) P. A. Abell _et al._ (LSST Science, LSST Project), (2009), arXiv:0912.0201 [astro-ph.IM] .
  * Gardner  _et al._ (2006) J. P. Gardner _et al._ , Space Sci. Rev. 123, 485 (2006), arXiv:astro-ph/0606175 .
  * Birrer  _et al._ (2015) S. Birrer, A. Amara, and A. Refregier, Astrophys. J. 813, 102 (2015), arXiv:1504.07629 [astro-ph.CO] .
  * Foreman-Mackey (2016) D. Foreman-Mackey, The Journal of Open Source Software 1 (2016), 10.21105/joss.00024.
  * Robitaille  _et al._ (2013) T. P. Robitaille _et al._ (Astropy), Astron. Astrophys. 558, A33 (2013), arXiv:1307.6212 [astro-ph.IM] .
  * Price-Whelan  _et al._ (2018) A. M. Price-Whelan _et al._ , Astron. J. 156, 123 (2018), arXiv:1801.02634 .
  * Bartelmann (1996) M. Bartelmann, Astron. Astrophys. 313, 697 (1996), arXiv:astro-ph/9602053 .
  * Schneider  _et al._ (1992) P. Schneider, J. Ehlers, and E. E. Falco,  _Gravitational Lenses_ (1992).
  * Bartelmann and Schneider (2001) M. Bartelmann and P. Schneider, Phys. Rept. 340, 291 (2001), arXiv:astro-ph/9912508 .
  * Golse and Kneib (2002) G. Golse and J.-P. Kneib, Astron. Astrophys. 390, 821 (2002), arXiv:astro-ph/0112138 .
  * Dumet-Montoya  _et al._ (2012) H. S. Dumet-Montoya, G. B. Caminha, and M. Makler, Astron. Astrophys. 544, A83 (2012), arXiv:1208.5682 [astro-ph.CO] .
  * Binney and Tremaine (2008) J. Binney and S. Tremaine,  _Galactic Dynamics: Second Edition_ (2008).