# Galaxy lens reconstruction based on strongly lensed gravitational waves:
similarity transformation degeneracy and mass-sheet degeneracy

Jason S.C. Poon1, Stefano Rinaldi2,3,4,5, Justin Janquart6,7, Harsh Narola6,7,
Otto A. Hannuksela1
1Department of Physics, The Chinese University of Hong Kong, Shatin, New
Territories, Hong Kong
2Dipartimento di Fisica “E. Fermi”, Università di Pisa, Largo Bruno Pontecorvo
3, 56127 Pisa, Italy
3INFN, Sezione di Pisa, Largo Bruno Pontecorvo 3, 56127 Pisa, Italy
4Institut für Theoretische Astrophysik, ZAH, Universität Heidelberg, Albert-
Ueberle-Str. 2, 69120 Heidelberg, Germany
5Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, Via
Marzolo 8, 35121 Padova, Italy
6Institute for Gravitational and Subatomic Physics (GRASP), Department of
Physics, Utrecht University, Princetonplein 1, 3584 CC Utrecht, The
Netherlands
7Nikhef - National Institute for Subatomic Physics, Science Park 105, 1098 XG
Amsterdam, The Netherlands E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Gravitational wave (GW) galaxy lens reconstruction is a crucial step for many
GW lensing science applications. However, dark siren GW lensing (i.e. lensed
GW without observed electromagnetic (EM) counterpart) suffers from similarity
transformation degeneracy and mass-sheet degeneracy. We review these two
degeneracies and discuss their implications on GW-based lens reconstruction
and two well-known GW lensing science cases: the Hubble constant measurement
and test for modified GW propagation. Building upon previous works, our
conclusions are: 1) GWs can only infer the scale-free lens mass model
parameters, the dimensionless source position, the GW luminosity distance and
the time delay scaling (a combination of Einstein radius, lens redshift, and
cosmology). 2) Lens reconstruction (of singular isothermal ellipsoid lens)
with only two GW signals is unlikely to yield a complete lens model, while
four (three) signals can measure all the above parameters accurately (with
large uncertainties). 3) The similarity transformation degeneracy causes the
lens redshift/Einstein radius/cosmology to be degenerate in dark siren
measurements. Breaking the degeneracy can be achieved by supplementing the GWs
with EM observation of lens redshifts/Einstein radius (source redshift is not
required). 4) The mass-sheet degeneracy causes the GW luminosity distance to
be entirely degenerate with a constant mass sheet. 5) Contrary to expectation,
the Hubble constant is degenerate with the mass-sheet even when supplemented
with lens reconstruction/redshift/Einstein radius and can only be lifted with
lens galaxy velocity dispersion measurement, while modified GW propagation
test discussed in prior literature is unaffected by the degeneracy. These
properties highlight the need for GW observations to be supplemented by EM
observations, which could become accessible through a lens archival search or
a rapid EM follow-up.

###### keywords:

gravitational lensing: strong – gravitational waves – cosmological parameters

††pubyear: 2024††pagerange: Galaxy lens reconstruction based on strongly
lensed gravitational waves: similarity transformation degeneracy and mass-
sheet degeneracy–21

##  1 Introduction

The gravitational lensing of light has been captured by electromagnetic (EM)
telescopes and studied for many decades, and it has become a standard yet
powerful tool in modern astronomy (Wambsganss, 1998; Bartelmann, 2010). Its
mature development and methodology allow scientists to probe the universe,
with some of its science applications including cosmological studies (such as
Hubble constant measurement (Williams & Schechter, 1997)), dark matter studies
(for instance, testing if massive compact halo object (MACHO) is a possible
candidate (Basak et al., 2022)), lens galaxy density profile and substructures
studies, as well as galaxy cluster structure studies (Bartelmann, 2010).

Meanwhile, according to general relativity, gravitational waves (GWs) are
expected to be lensed like light when they propagate near massive objects (De
Paolis et al., 2002). GW science itself is undergoing rapid development since
its first direct detection by Laser Interferometer Gravitational-Wave
Observatory (LIGO) in 2015 (Abbott et al., 2016): it opened up new windows to
probe and study countless scientific problems, including neutron star and
black hole physics (such as their structure and formation channels), testing
general relativity and cosmological studies such as expansion of universe and
dark energy (Sathyaprakash & Schutz, 2009; Lattimer, 2012; Gair et al., 2013;
Yunes & Siemens, 2013; Will, 2014; Berti et al., 2015; Lasky, 2015; Abbott et
al., 2019a, b; Baiotti, 2019; Abbott et al., 2021a; Bailes et al., 2021;
Abbott et al., 2022). Thanks to technological advancement, current GW
detectors (i.e. the second generation detectors), such as aLIGO (Aasi et al.,
2015), aVIRGO (Acernese et al., 2014) and KAGRA (Akutsu et al., 2020; Aso et
al., 2013; Somiya, 2012), are undergoing substantial upgrades (Cahillane &
Mansell, 2022; Bagnasco, 2023), and more new and advanced detectors are being
planned, such as the ground-based detector Einstein Telescope (Punturo et al.,
2010) and Cosmic Explorer (Reitze et al., 2019), as well as the space-based
detector Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al.,
2017).

One avenue that has garnered recent interest is the gravitational lensing of
GWs. Several groups have forecasted that, at their design sensitivity, second-
generation detectors would detect strongly-lensed GW with a detection rate of
around one such event per year, assuming standard formation channels (Ng et
al., 2018; Li et al., 2018; Ng et al., 2018; Oguri, 2018; Wierda et al., 2021;
Xu et al., 2022; Smith et al., 2023; Gupta et al., 2023). These forecasts have
helped motivate a recent effort to search for GW lensing (see Hannuksela et
al. (2019); McIsaac et al. (2020); Dai et al. (2020); Liu et al. (2021); Diego
et al. (2021); Abbott et al. (2021b, 2023); Janquart et al. (2023b)) and
develop scientific applications for them.

Indeed, if detected, GW lensing might allow several applications in
fundamental physics, astrophysics, and cosmology. For example, GW lensing
might allow sub-arcsecond localization of GW sources (Hannuksela et al., 2020;
Wempe et al., 2022). Other applications include tests of modified GW
propagation (Finke et al., 2021; Narola et al., 2023), testing the speed of
gravity (Collett & Bacon, 2017; Fan et al., 2017), GW polarization test (Smith
et al., 2017; Goyal et al., 2021a; Hernandez, 2022), other GW properties (Yang
et al., 2019; Mukherjee et al., 2020a, b; Chung & Li, 2021; Goyal et al.,
2021a; Hernandez, 2022; Iacovelli et al., 2022; Colaço et al., 2023),
cosmography (Sereno et al., 2010; Liao et al., 2017; Li et al., 2019; Yang et
al., 2019; Hou et al., 2021; Cao et al., 2022; Wang et al., 2022; Wu et al.,
2023) including high-redshift cosmography (Jana et al., 2023), probing dark
matter, subhalos, and compact objects (Takahashi, 2006; Dai et al., 2018; Liao
et al., 2018; Diego, 2020; Oguri & Takahashi, 2020; Choi et al., 2021; Wang et
al., 2021b; Basak et al., 2022; Cao et al., 2022; Guo & Lu, 2022; Oguri &
Takahashi, 2022; Fairbairn et al., 2023; Liu, 2023; Liu et al., 2023; Seo et
al., 2023; Tambalo et al., 2023), associating fast radio burst source and GW
source (Singh et al., 2023), constraining binary system formation channel
(Mukherjee et al., 2021). Indeed, detections could enable a plethora of
applications.

For many science applications, such as those in cosmography and tests of
modified GW propagation, "reconstruction" of the gravitational lens is one of
the crucial steps in the analysis. In particular, lensed GWs encode
information about the mass profile of the objects that gravitationally lensed
them. The attempt to reconstruct this lens mass profile and other extrinsic
properties of the lens is collectively called "lens reconstruction." Such
reconstructions allow one to study the lens model itself, including its
profile and the possible substructures within.

However, strong GW lensing suffers from two major degeneracies which limit
these reconstructions. The first one is the similarity transformation
degeneracy(Gorenstein et al., 1988), which refers to the fact that lensed GWs
cannot disentangle the size of the gravitational lens from the angular
diameter distances. I.e. the lens size is degenerate with the distances
between the lens, source and observer. The second degeneracy is mass-sheet
degeneracy (Falco et al., 1985), which refers to the fact that adding a
"constant mass sheet" to the lens profile will leave the image positions and
flux ratios invariant. In GW observations, as we will see, such a constant
mass sheet is degenerate with the inferred luminosity distance and the overall
time delay. These degeneracies limit the scientific applications of GW lensing
and what can be done with GWs alone.111Note that wave optics lensing and GWs
supplemented by EM information are an exception (Hannuksela et al., 2020;
Cremonese et al., 2021) (Chen et al., in prep.).

The similarity transformation degeneracy was originally reviewed by Gorenstein
et al. (1988) and Saha (2000), but the effect has yet to be discussed
significantly in the context of lensed GW observations. In this work, we would
introduce this degeneracy, but more in the context of GWs, based on their
findings. Due to the similarity transformation degeneracy, lens systems with
different redshifts, sizes and cosmological parameters (such as the Hubble
constant) can produce identical lensing observables. In particular, the size
of the lens is degenerate with the distances in the source-lens-observer
system. The degeneracy can be broken when the "size" of the lens can be
obtained from complementary EM observation. For example, a complementary
telescope image of the lensing object can break the degeneracy (because the
lens size and a redshift of the source and the lens can be directly observed
in the telescope). Such a telescope image might be available if the lens
system is located by an archival strong lens catalog222Such as, with lenscat
https://github.com/lenscat/lenscat. search (Hannuksela et al., 2020; Wempe et
al., 2022; Shan et al., 2023), which could be performed without an EM
counterpart. Alternatively, a rapid EM followup, if an EM counterpart to the
event is available, would also allow us to retrieve such a telescope image
(Smith et al., 2023). If no supplementary EM information is available, the
degeneracy limits GW applications. Thus, it is crucial to understand the
effect.

The second degeneracy is the mass-sheet degeneracy, which leaves the image
positions, flux ratios, and time-delay ratios unchanged when a uniform and
infinitely-spanning mass sheet is added to the lens model (Falco et al.,
1985). Such mass sheet can approximate scenarios where a galaxy lens is
embedded in low-mass galaxy cluster, there are massive structures along the
line of sight, or the lens has a positive power-law slope that differs from an
isothermal profile (Falco et al., 1985; Schneider & Sluse, 2013; Birrer et
al., 2016; Oguri, 2019). It has been a known degeneracy affecting strong lens
modelling (of light) for many decades (Schneider & Sluse, 2013), and it
affects GW lens reconstruction the same way. In particular, the constant mass
sheet will be degenerate with the luminosity distance and time delay measured
by the GWs. However, as we will see, complementary information about the
source redshift can alleviate the degeneracy, but not in cosmological
applications. Even when the source redshift is known, the degeneracy causes
the Hubble constant to be degenerate with the mass-sheet density. Another
intriguing way out of the degeneracy is by measuring the wave optics effects
when the lens size is comparable to the wavelength of the GW (Cremonese et
al., 2021) (Chen et al., in prep.). Without wave optics effects or
complementary EM redshift information, the degeneracy limits many GW
applications, so it is important to understand the effect.

This work reviews the idea of GW-based lens reconstruction/modelling of
galaxy-type lenses. We also discuss the nature and effect of the two
degeneracies and how one should consider them in the GW lensing context. In
Section 2, we first review the basics of strong gravitational lensing for both
light and GWs, before focusing more on its effect on GWs. Section 3 explains
lens reconstruction for lensed GW signal using a Bayesian framework. Section 4
introduces the similarity transformation degeneracy and discusses its impact
on GW-based lens reconstruction. In Section 5, we discuss the mass-sheet
degeneracy problem in the context of GW lensing. In Section 6, we present the
impact of the degeneracies on the test for modified GW propagation using lens
reconstruction. We also discuss their implications for various other science
cases.

Throughout the paper, we assume a flat $\Lambda$CDM universe with cosmology
from Plank18 data: $\Omega_{m}=0.308,\Omega_{K}=0,\Omega_{\Lambda}=0.692$
(Aghanim et al., 2020), with $H_{0}$ being a free (and unknown) parameter
unless stated.

##  2 Basics of Gravitational Wave Strong Lensing

In section 2.1, we will first review the basics of strong gravitational
lensing in geometrical optics, which applies to both EM waves and GWs as long
as the wavelength is much shorter than the scale of the lens object. If the
wavelength is comparable or longer than the characteristic lens scale, wave
optics (which includes diffraction and interference effect) must be considered
(Nakamura & Deguchi, 1999; Takahashi & Nakamura, 2003). However, because the
typical size of strong galaxy lenses is much larger than the wavelength of GWs
measurable with ground-based detectors, we limit our discussion to the
geometrical optics effects. We will also neglect microlensing effects.333Note
that microlensing _could_ become relevant in actual galaxy lensing scenarios
and will likely require further investigation (Diego et al., 2019; Cheung et
al., 2021; Mishra et al., 2021; Meena et al., 2022; Yeung et al., 2023; Mishra
et al., 2024). Then, in section 2.2, we will review the effect of strong
gravitational lensing on GWs, especially on the observation end.

###  2.1 Strong gravitational lensing

Figure 1: Geometry of a lens system: $D_{l},D_{s},D_{ls}$ are the angular
diameter distances from observer to lens, observer to source and lens to
source, respectively. Aligning our principal axis with the centre of the lens,
$\vec{\beta},\vec{\theta}$ are the angular position of the source and image
respectively, with $\vec{\eta}=D_{s}\vec{\beta}$ and
$\vec{\xi}=D_{l}\vec{\theta}$ being the (Cartesian) source and image
position(s) (on source plane and lens plane) respectively.

Fig. 1 illustrates the typical setup of a lens system. We denote the angular
diameter distance from the observer to the lens, observer to source and lens
to source as $D_{l},D_{s},D_{ls}$ respectively. The redshift of lens and
source are denoted as $z_{l}$ and $z_{s}$, respectively. We also denote the
angular position of the source and image by $\vec{\beta}$ and $\vec{\theta}$,
respectively. The characteristic scale of the lens is the angular Einstein
radius $\theta_{E}$ (some literature may refer to this as the Einstein angle),
which depends on the lens mass and the angular diameter distances in all
generality.

In this work, we should work in dimensionless coordinates of the lens system:
we define the dimensionless source and image positions as
$\vec{y}=(y_{1},y_{2})=\vec{\beta}/\theta_{E}$ and
$\vec{x}=(x_{1},x_{2})=\vec{\theta}/\theta_{E}$ respectively. The
dimensionless deflection/lensing potential $\Psi$ (deflection potential scaled
by $\theta_{E}^{2}$) is then defined by

$\Psi(\vec{x})=\frac{1}{\pi}\int{\kappa(\vec{x}^{\prime})\text{ln}|\vec{x}-\vec{x}^{\prime}|d^{2}x^{\prime}},$
(1)

where the convergence $\kappa$ is defined by

$\kappa(\vec{x})=\frac{\Sigma(\vec{\xi}=D_{l}\theta_{E}\vec{x})}{\Sigma_{cr}},$
(2)

where $\Sigma(\vec{\xi})$ is the surface mass density of the lens (in
cartesian coordinate $\vec{\xi}$ on the lens plane) and
$\Sigma_{cr}={c^{2}D_{s}}/{4\pi GD_{l}D_{ls}}$ is the critical density
(Takahashi & Nakamura, 2003; Tambalo et al., 2023). This particular definition
highlights the fact that $\Psi(\vec{x})$ does not contain scaling information,
such as the total mass of the lens. As we will see, the notation is
advantageous because the dimensionless quantities are unaffected by the
similarity transformation degeneracy.

The lens equation in dimensionless form can then be written as

$\vec{y}=\vec{x}-\nabla_{\vec{x}}\Psi(\vec{x}),$ (3)

which solving for $\vec{x}$ will give us the (dimensionless) image position(s)
under a known source position and lens model (there can be multiple images in
general).

Meanwhile, the time delay $t_{d}$ of an image compared to an unlensed ray is
given by

$t_{d}(\vec{x})=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}\left[\frac{|\vec{x}-\vec{y}|^{2}}{2}-\Psi(\vec{x})\right],$
(4)

where ${|\vec{x}-\vec{y}|^{2}}/{2}$ accounts for the delay due to extra
distance travelled by the lensed ray, and $\Psi(\vec{x})$ accounts for the
delay due to gravitational time dilation effect (Takahashi & Nakamura, 2003;
Bulashenko & Ubach, 2022). Also, the magnification $\mu$ of the image is given
by

$\mu(\vec{x})=\left[\det\left(\frac{\partial\vec{y}}{\partial{\vec{x}}}\right)\right]^{-1},$
(5)

where $\frac{\partial\vec{y}}{\partial{\vec{x}}}$ is given by

$\frac{\partial\vec{y}}{\partial{\vec{x}}}=\begin{pmatrix}1-\frac{\partial^{2}\Psi}{\partial
x_{1}^{2}}&-\frac{\partial^{2}\Psi}{\partial x_{1}\partial x_{2}}\\\
-\frac{\partial^{2}\Psi}{\partial x_{2}\partial
x_{1}}&1-\frac{\partial^{2}\Psi}{\partial x_{2}^{2}}\end{pmatrix}.$ (6)

We note that the magnification of an image is scale-free, which means it only
depends on the dimensionless system.

###  2.2 Effect of strong lensing on gravitational waves

Denoting $\tilde{h}_{0}(f)$ the original unlensed GW waveform in the frequency
domain and $\tilde{h}_{L}(f)$ the observed lensed ones, the relation between
the two waveforms is

$\tilde{h}_{L}(f)=\tilde{h}_{0}(f)\sum_{j}|\mu(\vec{x}_{j})|^{1/2}\exp{\left[2\pi
ift_{d}(\vec{x}_{j})-i\pi n_{j}\text{sign}(f)/2\right]},$ (7)

where the index $j$ denotes the $j$-th lensed image/signal, and $n_{j}$ is its
Morse index, with $n_{j}=\\{0,1,2\\}$ (referred as type I, II and III image)
when the image $\vec{x}_{j}$ is at a local minimum, saddle point and local
maximum of the time delay surface $t_{d}(\vec{x}_{j})$ respectively (Takahashi
& Nakamura, 2003). We also note that the emergence of Morse indices results
from stationary phase approximation in geometrical optics (Nakamura & Deguchi,
1999).

Therefore, in time domain, if $n_{j}=0$ or $2$, the lensed signal $h_{L}(t)$
and unlensed signal $h_{0}(t)$ are related by

$h_{L}(t)=\sum_{j}|\mu(\vec{x}_{j})|^{1/2}h_{0}\left[t-t_{d}(\vec{x}_{j})\right](-1)^{n_{j}/2},$
(8)

and if $n_{j}=1$, the relation is

$h_{L}(t)=\sum_{j}|\mu(\vec{x}_{j})|^{1/2}\int^{+\infty}_{-\infty}{-\tilde{h}_{0}(f)\text{sign}(f)ie^{-i2\pi
f\left[t-t_{d}(\vec{x}_{j})\right]}},$ (9)

which different images arrive at different times and have different amplitudes
(due to magnification) (Dai & Venumadhav, 2017; Ezquiaga et al., 2021).444If
the GW system is symmetric and contains no higher-order modes, the Morse phase
is absorbed into the phase of coalescence (Dai & Venumadhav, 2017; Ezquiaga et
al., 2021). However, they will have the same frequency-domain evolution (up to
magnification and the phase shift (Morse phase) caused by the Morse index)
(Takahashi & Nakamura, 2003), which is what strong lensing GW searches rely
on, trying to find events with compatible characteristics (Haris et al., 2018;
Goyal et al., 2021b; Janquart et al., 2021a; Kim et al., 2021; Janquart et
al., 2022; More & More, 2022; Janquart et al., 2023a; Li et al., 2023a, b;
Chakraborty & Mukherjee, 2024). Such searches have already been carried out on
real data, but there is no universally accepted evidence for GW lensing
detections yet, although many authors have considered the situation
(Hannuksela et al., 2019; Liu et al., 2021; Abbott et al., 2021b; Bianconi et
al., 2023; Diego et al., 2021; Janquart et al., 2023b; Abbott et al., 2023).

For a lens and GW source system, assuming we successfully identified all of
its lensed images (i.e. we identify the lensed GW signals from the detector
data stream), we can then measure the lensing observables, namely the
magnifications, time delays, and Morse phases of the strong lensing images.
Firstly, since GW is a standard siren (Holz & Hughes, 2005), from the lensed
GW images (denote the i-th image by subscript hereafter), we can observe the
effective (or apparent) luminosity distance $d_{L,i}^{\text{eff}}$, which is
related to the image magnification $\mu_{i}$ by

$d_{L,i}^{\text{eff}}=\frac{d_{L}}{\sqrt{|\mu_{i}|}},$ (10)

where $d_{L}$ is the true luminosity distance without the lensing effect. Note
that, in many traditional strong lensing scenarios, the absolute flux of the
light source is not know.555Exceptions to this would be strong lensing of
standard candles such as supernovae. In such scenarios, one can only measure
the relative magnification ratios ${\mu_{i}}/{\mu_{j}}$ between two images
(from the flux ratios) (Congdon & Keeton, 2018). The absolute flux or the
effective luminosity distance are not accessible. However, in lensed GWs (and
standard candles), one can measure the effective luminosity distance. This
allows for an extra degree of information on $d_{L}$. Secondly, in GW lensing,
the relative time delay between two images, $\Delta t_{d,ij}=t_{d,i}-t_{d,j}$
can also be measured with submillisecond accuracy (Fan et al., 2017). Hence,
the precise time delay can further constrain the lens model. Thirdly, by
comparing the waveform of the lensed GW signals, we can also observe the
relative Morse index (or difference in Morse index) $\Delta
n_{ij}=n_{i}-n_{j}$ between different images (Janquart et al., 2021b). The
direct GW observables are then the effective luminosity distances, relative
time delays, and relative (or in some cases, absolute 666Often we can only
accurately resolve the _relative_ Morse phase due to a degeneracy between the
Morse phase and GW coalescence phase (Dai et al., 2020; Lo & Magaña Hernandez,
2023; Janquart et al., 2021b), which however can be lifted in the presence of
higher-order modes.) Morse phases of the strong lensing images.

Figure 2: Illustration of dark siren GW lensing: the left panel shows the
setup of the lens system, which the GWs will arrive the observer (obs.) at
different time, with different amplitude and different Morse index. The
observables in a dark siren GW lensing scenario are shown on the right panel:
from the data recorded by the GW detector, one can measure the relative time
delay $\Delta t_{d,ij}$ and the relative Morse index $\Delta n_{ij}$ between
images, and since GW is a standard siren, from the amplitude of the GW signal
(together with waveform analysis), one can measure the effective/apparent
luminosity distance $d_{L,i}^{\text{eff}}$ of the images.

In many traditional strong lensing scenarios, the telescope angular resolution
is sufficiently good to locate and resolve the multiple strong lensing images.
However, due to the broad sky localization area of a GW signal (order of ten
to hundreds of square degrees in the near future (Pankow et al., 2018)),
unlike in many traditional lensing observation, we would not be able to
directly localize the precise sky location of the GW signals/images
individually from the GW data alone (Hannuksela et al., 2020). Therefore, we
will henceforth name “dark siren GW lensing” as the scenario in which we only
have observation of lensed GW from the GW detectors but do not have any EM
information on the lens and source system due to lack of EM counterpart for
the compact binary merger or broad sky localization 777However, we note that
it is still possible to identify the dark siren GW lens system to the same
lens system that lensed the GW source’s host galaxy if it is observed by EM
telescopes. The identification would require statistical comparison of the
lensing observables measured in the two channels. It will be discussed in
Section 3.. All the GW data will be denoted $\vec{d}_{GW}$. The data includes
information regarding the images’ effective luminosity distance
$d_{L,i}^{\text{eff}}$, the relative time delay between image pairs $\Delta
t_{d,ij}$ and the relative Morse index $\Delta n_{ij}$ between different
images (Fig. 2).

##  3 Lens Reconstruction Using Bayesian Analysis

This section reviews the basics of lens reconstruction using a Bayesian
analysis framework for GW lensing. In particular, we will focus on galaxy
lenses. We note that similar methods have been adopted by various works
previously, for example, Hannuksela et al. (2020); Wempe et al. (2022); Wright
& Hendry (2022); Seo et al. (2023); Wright et al. (2023).

There are various ways to model lenses, with some approaches being parametric
and others non-parametric. Based on strong galaxy lensing quasar measurements,
it was realized early on that even relatively simple smooth lens profiles
(such as the singular isothermal ellipsoid (SIE) model (Kormann et al., 1994))
fitted the time delays of quasars relatively well, but the fluxes were
difficult to be reproduced due to so-called flux-ratio anomalies (Mao &
Schneider, 1998; Metcalf & Madau, 2001). Since these early works, more
sophisticated parametric models have been developed to study strong galaxy
lenses, accounting for additional degrees of freedom such as parameters to
describe substructures, shearing, convergence, and non-singular density
profiles (Congdon & Keeton, 2018). However, to first order (and neglecting
shears), one could assume that the general morphology of many strong galaxy
lenses is modelled reasonably well using a singular isothermal ellipsoid (SIE)
when supplemented with conservative modelling errors to account for the
additional degrees of freedom not included in the SIE model (Keeton et al.,
1997).888Ultimately, corrections to the SIE model are necessary; for example,
the addition of shear will introduce a larger spread on the expected, allowed
magnification values, which would need to be marginalised over. However, to
understand the degeneracies and the general idea behind strong GW lens
modelling, an SIE model with generous modelling uncertainties is expected to
provide a fair qualitative description and starting point.

After choosing a lens model, the gravitational lens is characterized by the
parameters of the lens. We denote the lens parameters $\vec{\Phi}_{l}$. The
lens parameters encode all information strictly related to the lens, including
the lens redshift $z_{l}$ and the parameters describing the lens surface mass
density $\Sigma$ (i.e. the mass distribution of the lens). For instance, for
point mass lens, the main lens parameter would be its mass $M_{\text{PML}}$
(i.e. $\vec{\Phi}_{l}=(z_{l},M_{\text{PML}})$). For SIE, we have
$\vec{\Phi}_{l}=(z_{l},f,v)$ where $f$ is the lens’ minor to major axis ratio
and $v$ is the velocity dispersion. In this work, we restrict our discussion
to the SIE model as an approximate representation of a galaxy lens.999The
effect of making the lens mass profile more complex will be discussed in a
later section.

The second set of parameters encompasses the GW source information
$\vec{\Phi}_{s}$, including the source position $\vec{\beta}$ and the source
redshift $z_{s}$. Since we only focus on the effects induced by lensing but
not the nature of the GW waveform, we do not need the binary orbital
information (such as its chirp mass and mass ratio).101010The other GW
waveform parameters are nevertheless important, but they do not enter into
lens modelling unless there is a relationship between, for example, the total
mass of the binary and the brightness of the galaxy (see Wempe et al. (2022),
for more information). The last set of parameters contains the cosmological
information $\vec{\Phi}_{c}$. For simplicity, we only include the Hubble
constant $H_{0}$ as a free parameter here (and we also only consider a flat
$\Lambda$CDM universe).111111In general, one could also include radiation
energy density $\Omega_{r}$, matter density $\Omega_{m}$ and dark energy
density $\Omega_{\Lambda}$ as free parameters. However, we limit the analysis
to only $H_{0}$ for simplicity. Table 1 summarises the different parameter
sets.

Type: |  $\vec{\Phi}_{l}$ (lens) |  $\vec{\Phi}_{s}$ (source) |  $\vec{\Phi}_{c}$ (cosmology)
---|---|---|---
Free parameters: | $z_{l},\Sigma$ | $z_{s},\vec{\beta}$ | $H_{0}$
Table 1: Three types of parameters in a GW lens system: the lens parameters
$\vec{\Phi}_{l}$, the GW source parameter $\vec{\Phi}_{s}$ and the
cosmological parameters $\vec{\Phi}_{c}$. Within $\vec{\Phi}_{l}$, $z_{l}$ is
the lens redshift and $\Sigma$ is the mass distribution of the lens, for which
one should assume a mass model (in this work, we assume the SIE model with
conservative time-delay and magnification modelling uncertainties, which
characterize the mass distribution by two parameters: the minor to major axis
ratio $f$ and galaxy velocity dispersion $v$). For $\vec{\Phi}_{s}$, note that
we are only considering the effect induced by lensing (which are magnification
and delayed arrival time of the signal, but not the detail of the signal or
waveform itself), hence for the GW source, only its redshift $z_{s}$ and
source position $\vec{\beta}$ are relevant. For $\vec{\Phi}_{c}$, in this
work, we only consider the Hubble constant $H_{0}$ as a free parameter in our
flat $\Lambda$CDM universe.

The idea of lens reconstruction is to reconstruct the lens information
$\vec{\Phi}_{l}$ based on given observables: in the case of dark siren GW
lensing, the observables are the effective luminosity distances
$d_{L,i}^{\text{eff}}$, relative time delays $\Delta t_{d,ij}$ and relative
Morse index (Morse index difference) $\Delta n_{ij}$ (Fig. 2). If the lens
system is located through, e.g., archival search, then a telescope image of
the source galaxy would be available. The telescope image would provide the
flux ratios (and therefore the magnification ratios ${\mu_{i}}/{\mu_{j}}$) as
a function of the image positions $\vec{\theta}_{i}$, $\vec{\theta}_{j}$. In
the case that an EM counterpart is available, the GW image positions and
fluxes could also be retrieved. Other supplementary information would include
the photometric/spectroscopic redshifts of the lens galaxy $z_{l}$ and the
source galaxy $z_{s}$, and the measurement of the galaxy lens velocity
dispersion $\sigma_{v}$. Nevertheless, without supplementary EM information,
the lens reconstruction can only rely on the effective luminosity distances,
time delays, and relative Morse phases retrieved from the GWs.

Aside from the lens profile, inferring source and cosmological information
$\vec{\Phi}_{s}\text{ and }\vec{\Phi}_{c}$ is also interesting. Lens
reconstruction turns out to be a very important tool or middle-step in various
science applications of lensing, including sub-arcsecond localization of GW
source, measuring Hubble constant and testing modified gravity effect on the
propagation of GWs (Hannuksela et al., 2020; Finke et al., 2021; Narola et
al., 2023). However, as we will see, even in the presence of a supplementary,
high-resolution telescope image and photometric/spectroscopic redshifts, the
mass-sheet degeneracy will limit cosmological applications, while, contrary to
expectations, the degeneracy does not affect modified GW propagation tests.

We use standard Bayesian parameter estimation to reconstruct the lens. In
particular, our goal is to infer the posterior distribution of the lens and
source parameters, and also cosmological information given measured data:

$P(\vec{\Phi}_{l},\vec{\Phi}_{s},\vec{\Phi}_{c}|\vec{d}_{GW}),$

where $\vec{d}_{GW}=\\{d_{L,i}^{\text{eff}},\Delta t_{d,ij},\Delta n_{ij}\\}$
is the GW lensing observables that we measured. For simplicity, denote all
system parameters as
$\vec{\Phi}_{all}=\\{\vec{\Phi}_{l},\vec{\Phi}_{s},\vec{\Phi}_{c}\\}$. As
mentioned, in this work, we will adopt parametric models121212Primarily the
SIE lens model. in reconstructing the lens. Hence, the dimensionless
deflection potential $\Psi(\vec{x})$ and Einstein radius $\theta_{E}$ will be
parametric131313I.e. $\Psi(\vec{x})\rightarrow\Psi(\vec{x};\vec{\Phi}_{l})$
and $\theta_{E}\rightarrow\theta_{E}(\sigma_{v},D_{l},D_{ls},D_{s})$.. So, we
denote the parameters that characterise the selected dimensionless lens
profile as $\vec{\Phi}_{\Psi}$ (for example, for an SIE lens,
$\vec{\Phi}_{\Psi}=f$), and now the lens parameters are represented by
$\vec{\Phi}_{l}=(\theta_{E},\vec{\Phi}_{\Psi},z_{l})$. Then, the system
parameters are effectively
$\vec{\Phi}_{all}=(\vec{\beta},z_{s},\theta_{E},\vec{\Phi}_{\Psi},z_{l},H_{0})$.
As a word of caution, however, some of these parameters will turn out to be
degenerate with one another when only GW observations (i.e. dark siren) are
available. We will nevertheless first proceed with the Bayesian formulation
using the full set of parameters.

The Bayes theorem allows one to formulate the inference problem. In
particular, using the Bayes theorem, we can write

$P(\vec{\Phi}_{all}|\vec{d}_{GW})=\frac{P(\vec{d}_{GW}|\vec{\Phi}_{all})P(\vec{\Phi}_{all})}{P(\vec{d}_{GW})},$
(11)

where $P(\vec{d}_{GW}|\vec{\Phi}_{all})$ is the likelihood of obtaining
measured data given system parameters, $P(\vec{\Phi}_{all})$ is the prior
(probability distribution $\vec{\Phi}_{all}$ set by prior knowledge) and
$P(\vec{d}_{GW})$ is the evidence. One can obtain the posterior
$P(\vec{\Phi}_{all}|\vec{d}_{GW})$ by nested sampling.141414Nested sampling
can be be used to infer the posterior parameters as long as a likelihood
$P(\vec{d}_{GW}|\vec{\Phi}_{all})$ (which included the measured data
$\vec{d}_{GW}$), and a prior $P(\vec{\Phi}_{all})$ is defined. The nested
sampling procedure would then output both the posterior
$P(\vec{\Phi}_{all}|\vec{d}_{GW})$ and the evidence $P(\vec{d}_{GW})$. The
evidence is important for model comparison (say, comparing the point mass lens
model against the SIS lens model as demonstrated in Wright et al. (2023), or
comparing the SIE lens model versus PEMD (power law elliptical mass density
profile) lens model), which we do not do here. We focus on the lens
reconstruction (inference of the lensing parameters) instead of model
selection/comparison.

We note that the likelihood $P(\vec{d}_{GW}|\vec{\Phi}_{all})$ cannot be
written down directly, as the GW data only contains the measured image
properties $d_{L,i}^{\text{eff}},\Delta t_{d,ij}$ and $\Delta n_{ij}$, but no
information about $\vec{\Phi}_{all}$ directly. Hence, we would need a mapping
to relate $\vec{d}_{GW}$ and $\vec{\Phi}_{all}$. We denote the image
properties as $\vec{\nu}_{im}=(d_{L,i}^{\text{eff}},\Delta t_{d,ij},\Delta
n_{ij})$. Note that image properties are functions of the system parameters:
$\vec{\nu}_{im}=\vec{\nu}_{im}(\vec{\Phi}_{all})$. Assuming a suitable noise
model, we can obtain $P(\vec{d}_{GW}|\vec{\nu}_{im})$ easily (as demonstrated
in previous works such as Liu et al. (2021); Janquart et al. (2023b); Lo &
Magaña Hernandez (2023)). To obtain the likelihood
$P(\vec{d}_{GW}|\vec{\Phi}_{all})$, we note that
$P(\vec{d}_{GW}|\vec{\Phi}_{all})=\int{P(\vec{d}_{GW}|\vec{\nu}_{im})P(\vec{\nu}_{im}|\vec{\Phi}_{all})d\vec{\nu}_{im}}$.
Since the mapping from $\vec{\Phi}_{all}$ to $\vec{\nu}_{im}$ is a one to one
mapping, we have
$P(\vec{\nu}_{im}|\vec{\Phi}_{all})=\delta(\vec{\nu}_{im}-\vec{\nu}_{im}(\vec{\Phi}_{all}))$.
Therefore, to sample the likelihood, one can set all the system parameters
$\vec{\Phi}_{all}$ as free parameters, then solve the image positions from the
lens equation (Eq.(3)), then calculate the image properties
$\vec{\nu}_{im}(\vec{\Phi}_{all})$ using Eq.(4), Eq.(5) and Eq.(10), and
relate it to the measured data by the likelihood
$P(\vec{d}_{GW}|\vec{\nu}_{im})$.

As mentioned, in dark siren GW lensing observation, due to various
degeneracies, we would not be able to infer the posterior
$P(\vec{\Phi}_{all}|\vec{d}_{GW})$ of all the parameters $\vec{\Phi}_{all}$ in
the lens system. Instead, only posterior distributions of some of all free
parameters and a combination of some free parameters can be inferred.
Nevertheless, the idea of lens reconstruction in the Bayesian framework
discussed here remains unchanged.

After dark siren GW lens reconstruction, one could attempt to localize the
corresponding lens-source system in a follow-up archival search (Hannuksela et
al., 2020). Such a joint analysis could lead to a positive identification of
the lens system in EM channels if the lens galaxy is found and it is
sufficiently unique such that a joint analysis between the EM lens
reconstruction and the GW data could rule out a sufficient number of lenses.
The methodology was demonstrated in Hannuksela et al. (2020) and Wempe et al.
(2022). If the source and host galaxy are localized, one could obtain the
redshifts of the source and the lens through photometric/spectroscopic
measurement, and a supplementary lens reconstruction that could resolve some
of the quantities that cannot be resolved from GWs, such as the Einstein
radius, in addition to giving a supplementary constraints on the other lensing
parameters (such as the axis ratio). However, the localization starts with a
dark siren GW lens reconstruction analysis. Furthermore, understanding what
type of information GW-based strong lensing reconstructions can provide by
themselves and what type of caveats or degeneracies they are subject to
without supplementary information, will help in outlining complementary
avenues between the two. Thus, in this work, we perform lens reconstruction
based on (injected) dark siren lensed GW observations.

To infer the lens parameters, we use the GW Bayesian inference tool Bilby
(Ashton et al., 2019) and the nested sampler PyMultiNest (Buchner et al.,
2014). To model the gravitational lensing effects, we use the multi-purpose
gravitational lensing tool lenstronomy (Birrer & Amara, 2018).

##  4 Similarity Transformation Degeneracy

The similarity transformation degeneracy is a degeneracy between the (angular
diameter) distances in the lens system and the angular scale (i.e. the
Einstein radius $\theta_{E}$) of the lens, and the main reason for this
degeneracy is due to the scale-free nature of many lensing observables (such
as the magnification and time delay ratio) (Saha, 2000). It was first
discussed in Gorenstein et al. (1988) from the perspective of lens equation,
and Saha (2000) provided a comprehensive review of the degeneracy from the
point of view of time delay surfaces. Here, we revisit this degeneracy in the
context of GW lensing and formulate it using the dimensionless lens system.
Because of the similarity transformation degeneracy, all (dark siren) strongly
lensed GW observations will be unable to disentangle the size of the lens from
the distances in the source-lens-observer system and cosmology. Complementary
information (either from EM observations or otherwise) breaks the degeneracy.

We recall that the time delay of an image (relative to unlensed ray) is given
by

$t_{d}(\vec{x})=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}\left[\frac{|\vec{x}-\vec{y}|^{2}}{2}-\Psi(\vec{x})\right],$

which we emphasize again, $\vec{x}$ and $\vec{y}$ are dimensionless image and
source position, and $\Psi(\vec{x})$ is the dimensionless deflection
potential. Also, recall that the solutions of the dimensionless image position
can be solved entirely from the dimensionless lens equation (Eq. (3)). Hence,
one can write down the time delay in the following form:

$t_{d}(\vec{x})=T_{*}\tau(\vec{x}),$ (12)

where we define $\tau(\vec{x})$ as the dimensionless (scale free) time delay:

$\tau(\vec{x})\equiv\left[\frac{|\vec{x}-\vec{y}|^{2}}{2}-\Psi(\vec{x})\right],$
(13)

and $T_{*}$ as the time delay scaling:

$T_{*}\equiv\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}.$ (14)

Note that the time delay scaling $T_{*}$ is closely related to the time delay
distance $D_{\Delta t}$ commonly used by other works by

$T_{*}=\frac{D_{\Delta t}}{c}\theta_{E}^{2},$ (15)

as $D_{\Delta t}$ is given by

$D_{\Delta t}\equiv(1+z_{l})\frac{D_{l}D_{s}}{D_{ls}}.$ (16)

Furthermore, assuming a flat universe (i.e. $\Omega_{k}=0$) with fixed
cosmology but a variable Hubble constant, we can see
$\frac{D_{l}D_{s}}{D_{ls}}$ would be a function of source and lens redshift
and the Hubble constant, as the angular diameter distance $D$ (at redshift $z$
away) is given by

$D(z,H_{0})=\frac{c}{(1+z)H_{0}}\int_{0}^{z}\frac{dz^{\prime}}{\sqrt{\Omega_{r}(1+z^{\prime})^{4}+\Omega_{m}(1+z^{\prime})^{3}+\Omega_{\Lambda}}},$
(17)

where $\Omega_{r}$, $\Omega_{m}$ and $\Omega_{\Lambda}$ are the (normalized)
present-day radiation energy density, matter density, and dark energy density,
respectively. Hence, the time delay scaling is a function of the lens and
source redshift, the (angular) Einstein radius, as well as the Hubble
constant: $T_{*}=T_{*}(z_{l},z_{s},\theta_{E},H_{0})$. All the information
containing the scale of the system (i.e. the redshifts, the angular diameter
distances and the Einstein radius of the lens), as well as cosmological
information (such as $H_{0}$) are coupled together in the time delay scaling
$T_{*}$. Therefore, in GW lensing, when we measure the relative time delay
between two images $i$ and $j$, $\Delta t_{d,ij}=T_{*}\tau_{ij}$, with
$\tau_{ij}=(\tau_{i}-\tau_{j})$ and $\tau_{i}=\tau(\vec{x}_{i})$, we will have
a common factor $T_{*}$ for all the image pairs. Since the factor $T_{*}$
enters as an overall multiplication in the time delay, one cannot disentangle
the angular diameter distances, the Einstein radius $\theta_{E}$ and the
Hubble constant $H_{0}$ from each other, as they are coupled together inside
$T_{*}$. In other words, different sets of
$\\{z_{l},z_{s},\theta_{E},H_{0}\\}$ can give the same value of $T_{*}$, and
we would not be able to separate the two by just examining $\Delta t_{d,ij}$.

Meanwhile, in GW lensing, we also measure the effective luminosity distance of
each image $d_{L,i}^{\text{eff}}={d_{L}}/{\sqrt{|\mu_{i}|}}$.151515We stress
again that $\mu_{i}$ is scale-free. We then also note that the original
luminosity distance would be a function of source redshift and Hubble
constant: $d_{L}=d_{L}(z_{s},H_{0})$, as

$d_{L}=(1+z_{s})\frac{c}{H_{0}}\int_{0}^{z_{s}}\frac{dz^{\prime}}{\sqrt{\Omega_{r}(1+z^{\prime})^{4}+\Omega_{m}(1+z^{\prime})^{3}+\Omega_{\Lambda}}}.$
(18)

We note that in GW observation of compact binary coalescence, we cannot
measure $z_{s}$ directly due to the source mass-redshift degeneracy.161616It
may be possible for binary neutron star merger if the equation of state is
known and tidal effects are analyzed, but we neglect these indirect methods in
the following analysis (Messenger et al., 2014; Pang et al., 2020). The role
of $d_{L}$ is very similar to $T_{*}$: for all the images and their observed
effective luminosity distances, $d_{L}$ would appear as a common factor, and
both the redshift $z_{s}$ and the Hubble constant $H_{0}$ are coupled inside
$d_{L}$. Therefore, the redshifts and the cosmology are typically inferred
using the directly measurable GW parameters.

Type of parameters |  $\vec{\Phi}_{l}$ (SIE lens) |  $\vec{\Phi}_{s}$ (GW source) |  $\vec{\Phi}_{c}$ (cosmology)
---|---|---|---
Setup 1 | $z_{l}=0.4,v=6\times 10^{-4}c,f=0.7$ | $z_{s}=1.8,\vec{\beta}\approx(4.45\times 10^{-3},6.35\times 10^{-3})\text{ (in arcseconds)}$ | $H_{0}=67.8\text{ km}\text{ s}^{-1}\text{ Mpc}^{-1}$
Setup 2 | $z_{l}=0.8,v=5.83\times 10^{-4}c,f=0.7$ | $z_{s}=1.8,\vec{\beta}\approx(2.63\times 10^{-3},3.76\times 10^{-3})\text{ (in arcseconds)}$ | $H_{0}=67.8\text{ km}\text{ s}^{-1}\text{ Mpc}^{-1}$
Setup 3 | $z_{l}=1.1,v=6.18\times 10^{-4}c,f=0.7$ | $z_{s}=1.7,\vec{\beta}\approx(1.72\times 10^{-3},2.45\times 10^{-3})\text{ (in arcseconds)}$ | $H_{0}=63.2\text{ km}\text{ s}^{-1}\text{ Mpc}^{-1}$

Lensing observables |  $d_{L,1}^{\text{eff}},d_{L,2}^{\text{eff}},d_{L,3}^{\text{eff}},d_{L,4}^{\text{eff}}$ (in Mpc) |  $\Delta t_{d,12},\Delta t_{d,23},\Delta t_{d,34}$ (in seconds)
---|---|---
Common result for setup 1, 2, 3 | $6383,6216,7260,7624$ | $-33629,-228960,-42685$

Table 2: Three different lens system setups leading to the same lensing
observables under similarity transformation degeneracy: for the SIE lens
parameters, $z_{l},v,f$ are the redshift, velocity dispersion and SIE axis
ratio respectively and $c$ is the speed of light. For the GW source,
$z_{s},\vec{\beta}$ are the redshift and the (angular) source position in
arcseconds respectively. Lastly, $H_{0}$ is the Hubble constant. The first and
second setup share the same source redshift and Hubble constant, but the
second setup has a larger lens redshift and smaller velocity dispersion
compared to the first. Meanwhile, the first setup has relatively low lens
redshift and also a smaller lens galaxy velocity dispersion $v$ compared to
the third setup. For the angular Einstein radius, the first setup has a larger
angular Einstein radius of $0.635$ arcseconds versus the second setup’s
$0.376$ arcseconds and the third setup’s $0.245$ arcseconds. Since the three
systems’ dimensionless system are the same (i.e. they share the same
dimensionless source position $\vec{y}$ and SIE axis ratio $f$), and they
share the same value of time delay scaling $T_{*}$ and luminosity distance
$d_{L}$, so they share the same observed effective luminosity distances
$d_{L,i}^{\text{eff}}$ and Morse index for each of the four images (with the
first two images having $n=0$ and the last two images having $n=1$), as well
as the relative time delay between the image pairs $\Delta t_{d,ij}$ as shown
above. The physical picture of the degeneracy represented by the first two
setups is illustrated in Fig. 3. Figure 3: Illustrative example of similarity
transformation degeneracy in GW lensing: the upper and lower lens system on
the left panel corresponds to Setup 1 and 2 described in Table 2 respectively.
They share the same dimensionless system, source redshift and Hubble constant,
but they have different lens redshift $z_{l}$ ($0.4$ versus $0.8$), velocity
dispersion $v$ ($6\times 10^{-4}c$ versus $5.83\times 10^{-4}c$) and hence
Einstein radius $\theta_{E}$ ($0.635$ arcseconds versus $0.376$ arcseconds).
However, as shown on the right panel and in Table 2, they lead to the same
measurement of relative time delays and effective luminosity distances of the
images as their dimensionless system are the same and the values of $T_{*}$
and $d_{L}$ also match.

###  4.1 Degenerate transformations

####  4.1.1 $z_{l}$-$\theta_{E}$ degeneracy

Let us illustrate the degeneracy between the distance and the size of the
lens. Suppose, that we introduce a transformation such that the distance to
lens changes in such a way that the overall time-delay distance changes by
factor $a$. Suppose, furthermore, that the Einstein radius of the lens changes
by factor $1/a^{1/2}$ with all other quantities kept constant.171717The
dimensionless source/image positions are kept fixed. Because the lens equation
only depends on these dimensionless quantities, they are not affected by the
transformation.  That is, we introduce the following transformation:

$\begin{split}&z_{l}\rightarrow z_{l}^{\prime}=x_{0}z_{l}\,,\\\
&\theta_{E}\rightarrow\theta_{E}^{\prime}=\theta_{E}/a^{1/2}\,,\\\ &D_{\Delta
t}\rightarrow D_{\Delta t}^{\prime}=D_{\Delta t}a\,.\end{split}$

(Here $x_{0}$ is chosen such that the Einstein radius gains a factor
$a^{-1/2}$.) In this case, the overall time delay does not change:

$T_{*}\rightarrow T_{*}^{\prime}=\frac{D_{\Delta
t}^{\prime}}{c}(\theta_{E}^{\prime})^{2}=\frac{D_{\Delta
t}}{c}\theta_{E}^{2}=T_{*}$

Since the magnification depends only on the dimensionless quantities, they do
not change. For this reason, the effective luminosity distances also do not
change. Thus, under the transformation, all GW lensing observables are
unchanged:

$\begin{split}&d_{L,i}^{\text{eff}}\rightarrow
d_{L,i}^{\text{eff},\prime}=d_{L,i}^{\text{eff}}\\\ &\Delta
t_{d,ij}\rightarrow\Delta t_{d,ij}^{\prime}=\Delta t_{d,ij}\,,\\\
&n_{j}\rightarrow n_{j}^{\prime}=n_{j}\,.\end{split}$ (19)

That is, the lens system before/after produces identical GW observables.

####  4.1.2 $H_{0}$-$\theta_{E}$ degeneracy

Alternatively, one can introduce a transformation in the Hubble constant and
the angular size of the lens (all other quantities kept constant):

$\begin{split}&H_{0}\rightarrow H_{0}^{\prime}=H_{0}/a\,,\\\
&\theta_{E}\rightarrow\theta_{E}^{\prime}=\theta_{E}/a^{1/2}\,,\\\ &D_{\Delta
t}\rightarrow D_{\Delta t}^{\prime}=D_{\Delta t}a\,,\end{split}$

which is likewise degenerate since

$\begin{split}&d_{L,i}^{\text{eff}}\rightarrow
d_{L,i}^{\text{eff},\prime}=d_{L,i}^{\text{eff}}\\\ &\Delta
t_{d,ij}\rightarrow\Delta t_{d,ij}^{\prime}=\Delta t_{d,ij}\,,\\\
&n_{j}\rightarrow n_{j}^{\prime}=n_{j}\,.\end{split}$ (20)

####  4.1.3 Partial $z_{s}$-$z_{l}$-$\theta_{E}$-$H_{0}$ degeneracy

Finally, there exists a non-trivial transformation involving the source
redshift transformations, the Hubble constant, the scale of the lens, when
keeping the luminosity distance fixed:

$\begin{split}&z_{s}\rightarrow z_{s}^{\prime}=x_{1}z_{s}\,,\\\
&H_{0}\rightarrow H_{0}^{\prime}=H_{0}/a\,,\\\ &d_{L}\rightarrow
d_{L}^{\prime}=d_{L}(z_{s}^{\prime},H_{0}^{\prime})=d_{L}\\\
&\theta_{E}\rightarrow\theta_{E}^{\prime}=\theta_{E}/b^{1/2}\,,\\\ &D_{\Delta
t}\rightarrow D_{\Delta t}^{\prime}=D_{\Delta t}b\,,\end{split}$

where $x_{1}$ is chosen such that $d_{L}^{\prime}=d_{L}$. The above
transformation likewise leads to the same GW observables:

$\begin{split}&d_{L,i}^{\text{eff}}\rightarrow
d_{L,i}^{\text{eff},\prime}=d_{L,i}^{\text{eff}}\\\ &\Delta
t_{d,ij}\rightarrow\Delta t_{d,ij}^{\prime}=\Delta t_{d,ij}\,,\\\
&n_{j}\rightarrow n_{j}^{\prime}=n_{j}\,.\end{split}$ (21)

However, it is worth noting that this final degeneracy is a _partial_
degeneracy, because $d_{L}$ must be kept fixed in the transformation. For this
reason, the degeneracy can be broken by measuring $z_{l}$ and $\theta_{E}$
only, even without $z_{s}$ measurement. 181818More technically, if there are
$4$ degrees of freedom in allowing for $z_{s}$, $z_{l}$, $\theta_{E}$, and
$H_{0}$ to vary, then the condition $d_{L}^{\prime}=d_{L}$ produces one
”constraint”, such that there are really only $4-1=3$ degrees of freedom, so a
measurement of two out of the four degenerate quantities is sufficient to
break the degeneracy. Section 4.4 will provide a more detailed formalism for
this partial degeneracy.

####  4.1.4 What quantities are degenerate?

From the above considerations, one should conclude that it is possible that
for two systems with very different absolute scale and distances as well as
different value of the Hubble constant (i.e. different set of
$\\{z_{l},z_{s},\theta_{E},H_{0}\\}$). That is, the scale of the lens
$\theta_{E}$, the redshift to the lens and source $z_{l},z_{s}$, and the
Hubble constant $H_{0}$ are entirely degenerate with one another. To further
illustrate the effect, Table 2 shows three very different SIE lens systems
that would lead to the same lensing observables. The first two systems
considered are illustrated in Fig. 3. Saha (2000) referred to the
transformation that induced the degeneracy presented above as the similarity
transformation; therefore, we call this degeneracy the similarity
transformation degeneracy. To back up these arguments, we will give further
evidence of these degeneracies in our lens reconstruction results, where the
degeneracies between $\\{z_{l},z_{s},\theta_{E},H_{0}\\}$ are obvious.

###  4.2 Lens Reconstruction with Similarity Transformation Degeneracy

What can GWs tell us about the galaxies that lens them? For dark siren GW
lensing (i.e. observation without EM counterpart or any other supplementary
information), we can only measure the relative time delays $\Delta t_{d,ij}$,
effective luminosity distances $d_{L,i}^{\text{eff}}$ and relative Morse
indices $\Delta n_{ij}$ of the images 191919The absolute Morse phase of type-
II images may also be measurable when the GW experiences higher-order-modes
(Dai & Venumadhav, 2017; Ezquiaga et al., 2021; Janquart et al., 2021b; Wang
et al., 2021a).. In this case, only a limited amount of information can be
reconstructed. We focus on galaxy lenses here: as the image with Morse index
$n_{j}=2$ (i.e. the type III image) is a highly demagnified (central) image
that can rarely be observed in galaxy lensing (Dahle et al., 2013; Collett &
Bacon, 2016; Collett et al., 2017), we neglect it. Meanwhile, the images with
Morse index $n_{j}=\\{0,1\\}$ correspond to positive and negative value of
magnification $\mu$ respectively (Vijaykumar et al., 2023). In addition, for a
typical galaxy lens system with four images, the arrival order of images is
type I-I-II-II (Schneider et al., 1992; Wierda et al., 2021). Hence, assuming
all four images are observed, we can sort out the image types by arrival
order, which would be included in the lens reconstruction process. GWs can
thus yield the effective luminosity distances $d_{L,i}^{\text{eff}}$, relative
time delays $\Delta t_{d,ij}$, and absolute Morse phases $n_{i}$ as
observables. However, it is important to understand what lensing quantities
these observables can map to; it is clear that due to the similarity
transformation degeneracy, not all lensing quantities can be
measured.202020For example, can we retrieve the Einstein radius? How about the
axis ratio? What about the redshifts?

From a reverse engineering point of view, after we obtain the GW observables
($d_{L,i}^{\text{eff}},\Delta t_{d,ij},n_{i}$), we can obtain the absolute
value of the magnification ratio of images by

$\left|\frac{\mu_{i}}{\mu_{j}}\right|=\left(\frac{d_{L,j}^{\text{eff}}}{d_{L,i}^{\text{eff}}}\right)^{2},$
(22)

which is completely determined by the dimensionless lens system212121That is,
it only depends on the scaled image/source positions and the scaled lens
potential, but not on the angular diameter distances or the Einstein radius
(scale of the lens). after $d_{L}$ is cancelled out. We note that, for our
galaxy lensing case where we ignore type III images, one can also obtain the
magnification ratio (signs included) by taking the Morse index into account:

$\frac{\mu_{i}}{\mu_{j}}=\left(\frac{d_{L,j}^{\text{eff}}}{d_{L,i}^{\text{eff}}}\right)^{2}(-1)^{\Delta
n_{ij}}.$ (23)

Meanwhile, we can also make use of the ratio of relative time delay of
different image pairs and obtain the ratio of dimensionless time delay:

$\frac{\tau_{ij}}{\tau_{mn}}=\frac{\Delta t_{d,ij}}{\Delta t_{d,mn}},$ (24)

which is completely fixed by the dimensionless lens system, meaning it does
not depend on the time-delay scaling $T_{*}$. Suppose we have sufficient
number of observed images and a parametric lens model (for instance, four
images and an SIE lens mass profile). In that case, we can reconstruct the
dimensionless lens system from the effective luminosity distance ratios and
relative time delay ratios. For instance, if we choose SIE to be our lens
model, using the magnification ratios ${\mu_{i}}/{\mu_{j}}$ and dimensionless
time delay ratios ${\tau_{ij}}/{\tau_{mn}}$, one would be able to reconstruct
the axis ratio $f$ (which is part of the dimensionless system), as well as the
dimensionless source position $\vec{y}$ (hence also the dimensionless image
positions $\vec{x}$ by Eq. (3)). In Section 4.5, we use a SIS lens to
demonstrate this idea analytically.

After reconstructing the dimensionless lens system, using Eq. (13), we can
calculate the dimensionless time delay $\tau_{i}$ of the i-th image. Then we
can obtain the value of $T_{*}$ from the observed $\Delta t_{d,ij}$ using Eq.
(12). Similarly, we can calculate the i-th image’s magnification $\mu_{i}$
from Eq. (5) and use the relation $d_{L}=d_{L,i}^{\text{eff}}\sqrt{|\mu_{i}|}$
to calculate $d_{L}$.

What, therefore, can GWs alone infer? GWs can be used to recover the
dimensionless system: dimensionless source position $\vec{y}$ and
dimensionless deflection potential parameters $\vec{\Phi}_{\Psi}$ from the
ratios ${\mu_{i}}/{\mu_{j}}$ and ${\tau_{ij}}/{\tau_{mn}}$. However, in terms
of the scale of the system, we can only recover the time-delay factor/scaling
$T_{*}$ and the absolute luminosity distance $d_{L}$. Therefore, we will not
be able to infer the actual physical scales and sizes (such as $D_{s}$,
$D_{l}$ and $\theta_{E}$) and the actual mass distribution $\Sigma$. In other
words, we can only reconstruct the posterior
$P(\vec{y},\vec{\Phi}_{\Psi},T_{*},d_{L}|\vec{d}_{GW})$ (under a selected lens
model and form of $\Psi(\vec{x})$), but not the posterior of all free
parameters $P(\vec{\Phi}_{all}|\vec{d}_{GW})$ that would allow inferring the
cosmological parameters from the lens system. Therefore, in the process of
lens reconstruction from dark siren lensed GW in the Bayesian framework, we
should only set $\\{\vec{y},\vec{\Phi}_{\Psi},T_{*},d_{L}\\}$ as the
(irreducible) free parameters in our lens model, instead of
$\vec{\Phi}_{all}=(\vec{\beta},z_{s},\theta_{E},\vec{\Phi}_{\Psi},z_{l},H_{0})$.
That is, the inference of $\\{\vec{y},\vec{\Phi}_{\Psi},T_{*},d_{L}\\}$ is a
data-driven, while inference of $\vec{\Phi}_{all}$ is a prior-driven.222222In
particular, prior information about $z_{l}$ and $\theta_{E}$ are required to
perform meaningful reconstructions with $\vec{\Phi}_{all}$ (such a prior-
information-driven approach was demonstrated in Jana et al. (2023), which
supports the assertion).

Type of parameters |  $\vec{\Phi}_{l}$ (SIE lens) |  $\vec{\Phi}_{s}$ (GW source) |  $\vec{\Phi}_{c}$ (cosmology)
---|---|---|---
Injection values | $z_{l}=0.5,v=0.0006c,f=0.5$ | $z_{s}=2,\vec{\beta}\approx(0.00592,0.00887)\text{ (in arcseconds)}$ | $H_{0}=67.8\text{ km}\text{ s}^{-1}\text{ Mpc}^{-1}$
Table 3: Injection values of the system parameters for our lens reconstruction example: 1. for the lens parameters, $z_{l},v,f$ are the redshift, velocity dispersion and SIE axis ratio, respectively. Note that $c$ is the speed of light and we align the long axis of the SIE with the vertical axis (i.e. the $y_{2}$ or $x_{2}$ axis). 2. For the source, $z_{s},\vec{\beta}$ are the redshift and the (angular) source position in arcseconds, respectively. 3. For cosmology, $H_{0}$ is the Hubble constant. Observables | Effective luminosity distance: $d_{L,1}^{\text{eff}},d_{L,2}^{\text{eff}},d_{L,3}^{\text{eff}},d_{L,4}^{\text{eff}}$ (in Mpc) | Relative time delay between images: $\Delta t_{d,12},\Delta t_{d,23},\Delta t_{d,34}$ (in s)
---|---|---
Injected values | $9404,9237,12596,13169$ | $-55817,-491735,-66594$
Table 4: GW lensing observables of injected system setup shown in Table 3.
Note that we have four images, and the image types are I-I-II-II. Figure 4:
Result of SIE lens reconstruction of dark siren lensed GW using injected
parameters shown in Table 3 and Table 4, which an uncorrelated Gaussian
likelihood with $5\%$ error on observed relative time delays and effective
luminosity distances is assumed. After taking similarity transformation
degeneracy into account, the five free parameters in the inference are
$\\{\vec{y}=(y_{1},y_{2}),f,T_{*},d_{L}\\}$, which are the dimensionless
source position, the SIE axis ratio, the time delay scaling
$T_{*}\equiv\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}$ and the
(unlensed) luminosity distance respectively. In this figure, we show the
posterior distribution of the parameters that represent the dimensionless lens
system: $\vec{y}=(y_{1},y_{2})$ and $f$. We can see, the parameters are
recovered fairly well. We note that the dimensionless source position suffers
from multiple solutions only because SIE lens has two degrees of reflectional
symmetry. Other than that, there is no degeneracy among the parameters.
Figure 5: Result of SIE lens reconstruction of dark siren lensed GW using
injected parameters shown in Table 3 and Table 4, which an uncorrelated
Gaussian likelihood with $5\%$ error on observed relative time delays and
effective luminosity distances is assumed. After taking similarity
transformation degeneracy into account, the five free parameters are
$\\{\vec{y}=(y_{1},y_{2}),f,T_{*},d_{L}\\}$, which are the dimensionless
source position, the SIE axis ratio, the time delay scaling
$T_{*}\equiv\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}$ and the
(unlensed) luminosity distance respectively. In this figure, we show the
posterior distribution of the parameters that represent the scale of the lens
system: $T_{*}\text{ (in seconds)}$ and $d_{L}\text{ (in Mpc)}$. We can see,
the two parameters are recovered fairly well.

Below, we demonstrate parameter estimation results using our "irreducible"
parameters, which accounts for the similarity transformation degeneracy. We
inject an SIE lens with the parameters described in Table 3. After obtaining
the injected observables for the four images as described in Table 4, we
perform Bayesian parameter estimation of the lens system. For the likelihood
$P(\vec{d}_{GW}|\vec{\nu}_{im})$, we assume a simple and uncorrelated 5%
Gaussian error on the measured $d_{L,j}^{\text{eff}}$ and $\Delta
t_{d,ij}$232323This is to roughly account for both measurement and lens model
uncertainties. Note that a correlated probability distribution of the measured
$d_{L,j}^{\text{eff}}$ and $\Delta t_{d,ij}$ that emerge from direct GW data
analysis (such as Haris et al. (2018); Hannuksela et al. (2019); Goyal et al.
(2021b); Janquart et al. (2021a, 2022, 2023a, 2023b); Lo & Magaña Hernandez
(2023); Abbott et al. (2023)) should be used in a quantitative analysis; we
are primarily interested in understanding the various lensing degeneracies
qualitatively, where a rough representative error estimate should be
sufficient. Under the similarity transformation degeneracy for dark siren GW
lensing, we would only sample five free parameters in this case: the two-
dimensional dimensionless source position $\vec{y}=(y_{1},y_{2})$, the SIE
axis ratio $f$, the luminosity distance $d_{L}$ and the time delay scaling
$T_{*}$. We also assume a uniform prior for all parameters for simplicity.
Fig. 4 and Fig. 5 show the lens reconstruction result of the system’s
dimensionless parameters ($\vec{y},f$) and parameters that represent the scale
of the system ($T_{*},d_{L}$) respectively. Firstly, there are four peaks in
the posterior for the dimensionless source position $\vec{y}$. This is to be
expected due to the reflectional symmetry of the SIE lens. Then, we can see
for all five parameters, the injected values are included in the recovered
posteriors. The full parameter estimation result of all five parameters is
shown in Fig. 15 in the Appendix.

We then repeat the exercise with 7 free parameters:
$\\{\vec{y}=(y_{1},y_{2}),f,\theta_{E},z_{s},z_{l},H_{0}\\}$, which does not
account for similarity transformation degeneracy. The setup is identical to
the run presented in Fig. 4 and Fig. 5, with a uniform prior on all sampling
parameters. Fig. 6 shows the posterior distribution of the parameters that
describe the scale of the system: namely the Einstein radius $\theta_{E}$, the
source and lens redshifts $z_{s},z_{l}$ and the Hubble constant $H_{0}$. The
full result of the parameter estimation is shown in Fig. 16 in Appendix. As
expected, the four quantities that describe the actual scale of the system,
$\\{\theta_{E},z_{s},z_{l},H_{0}\\}$, are poorly recovered and constrained due
to the degeneracies that are present among these parameters: for instance, the
$z_{s}\text{-}H_{0}$ relation is spanned across the entire parameter space
diagonally, implying that the two are fully degenerate (due to their
combination within $d_{L}$). Similarly, the $z_{l}\text{-}\theta_{E}$ pair is
also highly degenerate. Meanwhile, we note that the dimensionless system
parameters are still recovered as shown in Fig. 16. Therefore, despite having
7 constraints from the dark siren GW observations (4 $d_{L,j}^{\text{eff}}$
and 3 $\Delta t_{d,ij}$), that, naively, could be matched to the 7 free
lensing parameters ($f$, $v$, $\vec{y}$, $z_{s}$, $z_{s}$, $H_{0}$), one
cannot constrain the full system due to the similarity transformation
degeneracy. Indeed, these results support our earlier discussion.

Figure 6: Result of SIE lens reconstruction of dark siren lensed GW using
injected parameters shown in Table 3 and Table 4, without accounting for
similarity transformation degeneracy. An uncorrelated Gaussian likelihood with
$5\%$ error on observed relative time delays and effective luminosity
distances is also assumed. The seven free parameters are
$\\{\vec{y}=(y_{1},y_{2}),f,\theta_{E}\text{ (in
arcseconds)},z_{s},z_{l},H_{0}\text{ (in km}\text{ s}^{-1}\text{
Mpc}^{-1}\text{)}\\}$, which are the dimensionless source position, the SIE
axis ratio, the Einstein radius, the source and lens redshift and the Hubble
constant respectively. Here we show the posterior distribution of the four
parameters that represent the scale of the system: $\theta_{E}\text{ (in
arcseconds)},z_{s},z_{l},H_{0}\text{ (in km}\text{ s}^{-1}\text{
Mpc}^{-1}\text{)}$. They are poorly recovered (which can be sensitive to prior
choice) and showed high degree of degeneracies, which is expected from
similarity transformation degeneracy.

A detected, strongly lensed GWs would give us a measurement of the relative
time delays $\Delta t_{d,ij}$ and effective luminosity distances
$d_{L,i}^{\text{eff}}$. With just $P(\Delta t_{d,ij})$ and
$P(d_{L,i}^{\text{eff}})$, as demonstrated above, we would suffer from
similarity transformation degeneracies and can only infer a reduced number of
parameters of the system. This may not be very useful for science application
that rely on lens reconstruction. However, if the lens system is localised
through an archival search or a rapid EM follow-up, complementary EM
information of the corresponding lens system, especially the lens and source
redshift and Einstein radius of the lens, can break the similarity
transformation degeneracy. In Section 4.4, we demonstrate that such
identification and EM information is crucial to the Hubble constant
measurement.

###  4.3 Lens Reconstruction with Two or Three Images Only

Here, we provide additional examples of dark siren GW SIE lens reconstruction
attempt with only two or three of all four GW images are used in the process
of parameter estimation. In practice, this may have to be done when the
remaining lensed images have not arrived yet, or they are too dim to be
detected. We also assume the order of the two or three images are identified
or assumed correctly within all images. 242424To elaborate, we first obtain
the injected observables from, for instance, the first three images (in
practice, since the first two images are type I image and the third one is
type II image, we would not be able to distinguish if we observed the first
three images or the first two and fourth images from the lens system). Then we
attempt lens reconstruction by (correctly) assuming that the three images we
obtained are the first three of the four images. The setup of the lens system
is presented in Table 3 and Table 4.

Figure 7: Result of SIE lens reconstruction of dark siren lensed GW using
injected parameters shown in Table 3 and the observables of the first two
images from Table 4, which an uncorrelated Gaussian likelihood with $0.001\%$
error on observed relative time delays and effective luminosity distances is
assumed. After taking similarity transformation degeneracy into account, the
five free parameters in the inference are
$\\{\vec{y}=(y_{1},y_{2}),f,T_{*},d_{L}\\}$, which are the dimensionless
source position, the SIE axis ratio, the time delay scaling
$T_{*}\equiv\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}$ and the
(unlensed) luminosity distance respectively. In this figure, we show the
posterior distribution of the parameters that represent the dimensionless lens
system: $\vec{y}=(y_{1},y_{2})$ and $f$. We can see, the parameters cannot be
recovered and the posterior is uninformative.

Fig. 7 presented the result of parameter estimation of the dimensionless
system parameters $\vec{y}=(y_{1},y_{2})$ and $f$ in the case of only using
the first two images. The full result is shown in Fig. 17 in the Appendix. We
note that an uncorrelated Gaussian likelihood with only $0.001\%$ error on the
observables is assumed in order to minimize uncertainty. As shown, the
injected parameters are not recovered and the posteriors appears
uninformative. This is expected, as there is a mismatch in degree of freedom:
we only have three observables ($d_{L,1}^{\text{eff}},d_{L,2}^{\text{eff}}$
and $\Delta t_{d,12}$), but the entire lens system has five irreducible free
parameters ($\\{\vec{y}=(y_{1},y_{2}),f,T_{*},d_{L}\\}$).

Figure 8: Result of SIE lens reconstruction of dark siren lensed GW using
injected parameters shown in Table 3 and the observables of the first three
images from Table 4, which an uncorrelated Gaussian likelihood with $5\%$
error on observed relative time delays and effective luminosity distances is
assumed. After taking similarity transformation degeneracy into account, the
five free parameters in the inference are
$\\{\vec{y}=(y_{1},y_{2}),f,T_{*},d_{L}\\}$, which are the dimensionless
source position, the SIE axis ratio, the time delay scaling
$T_{*}\equiv\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}$ and the
(unlensed) luminosity distance respectively. In this figure, we show the
posterior distribution of the parameters that represent the dimensionless lens
system: $\vec{y}=(y_{1},y_{2})$ and $f$. We can see, the parameters can be
recovered, but with very large uncertainty.  Figure 9: Result of SIE lens
reconstruction of dark siren lensed GW using injected parameters shown in
Table 3 and the observables of the first three images from Table 4, which an
uncorrelated Gaussian likelihood with $0.01\%$ error on observed relative time
delays and effective luminosity distances is assumed. After taking similarity
transformation degeneracy into account, the five free parameters in the
inference are $\\{\vec{y}=(y_{1},y_{2}),f,T_{*},d_{L}\\}$, which are the
dimensionless source position, the SIE axis ratio, the time delay scaling
$T_{*}\equiv\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}$ and the
(unlensed) luminosity distance respectively. In this figure, we show the
posterior distribution of the parameters that represent the dimensionless lens
system: $\vec{y}=(y_{1},y_{2})$ and $f$. We can see, the parameters are
recovered fairly well. We note that the dimensionless source position suffers
from multiple solutions only because SIE lens has two degrees of reflectional
symmetry. Other than that, there is no degeneracy among the parameters.

Fig. 8 and Fig. 9 show the result of the inference of the lens system’s
dimensionless parameters using the first three images, with the difference
being they assume a $5\%$ and $0.01\%$ error on the (uncorrelated) Gaussian
likelihood of the observables respectively. The full results are also shown in
Fig. 18 and Fig. 19 in the Appendix respectively. From Fig. 8, we can see the
injection overlaps with the posterior, although the uncertainties are quite
broad. From Fig. 9, which the error of the Gaussian likelihood is small and
the parameter estimation uncertainty is minimized, we recover the injection
value accurately. This result is also expected, as the five degrees of
information in the observables
($d_{L,1}^{\text{eff}},d_{L,2}^{\text{eff}},d_{L,3}^{\text{eff}},\Delta
t_{d,12}$ and $\Delta t_{d,23}$) match the five degrees of freedom in the
irreducible free parameters ($\\{\vec{y}=(y_{1},y_{2}),f,T_{*},d_{L}\\}$).
However, we note that the inference accuracy will not be as good as the case
of using all four images. Nonetheless, it may still allow us to predict the
image property of the fourth image and help for preparing or searching for the
fourth image in the GW data, given that we identify the order of the images
correctly, as discussed in Ng et al. (2024).

###  4.4 Science Application: Hubble Constant Measurement

Here, we review the procedures of measuring Hubble constant from lensed GW in
the context of the various degeneracies. Demonstrations of such Hubble
constant measurements have been done by various authors in the context of
supplementary EM data Liao et al. (2017); Hannuksela et al. (2020) and
statistical studies Jana et al. (2023).

Suppose one has access to the EM counterpart or observation of the lens galaxy
and the lensed host galaxy of the GW source, for example, through an archival
search or a rapid EM follow-up Hannuksela et al. (2020); Smith et al. (2023).
In this case, we have access to the source and lens redshifts ($z_{s}$ and
$z_{l}$). With this additional information, one can measure the Hubble
constant $H_{0}$ as follows. Since the (unlensed) luminosity distance of the
source $d_{L}$ and angular diameter distance to the source $D_{s}$ are related
by $D_{s}={d_{L}}/{(1+z_{s})^{2}}$, $T_{*}$ can be written as

$T_{*}=\frac{1+z_{l}}{c}\frac{D_{l}}{D_{ls}}\frac{d_{L}}{(1+z_{s})^{2}}\theta_{E}^{2}=d_{L}\left[\frac{D_{l}}{D_{ls}}\frac{1+z_{l}}{c(1+z_{s})^{2}}\theta_{E}^{2}\right].$
(25)

For our flat universe, we note that the Hubble constant information is encoded
in the luminosity distance $d_{L}$ as shown in Eq.(18). Then, we write

$\begin{split}\Delta
t_{d,ij}&=T_{*}\tau_{ij}=d_{L}\left[\frac{D_{l}}{D_{ls}}\frac{1+z_{l}}{c(1+z_{s})^{2}}\theta_{E}^{2}\right]\tau_{ij}\\\
&=d_{L}G\tau_{ij},\end{split}$ (26)

where $G$ is defined as

$G=G(z_{l},z_{s},\theta_{E})\equiv\frac{D_{l}}{D_{ls}}\frac{1+z_{l}}{c(1+z_{s})^{2}}\theta_{E}^{2}.$
(27)

We stress that $G$ is independent of the Hubble constant (as it cancels out in
the factor ${D_{l}}/{D_{ls}}$). In addition, some may argue that the angular
Einstein radius $\theta_{E}$ may be dependent on $H_{0}$, for instance, for a
point mass lens with mass $M$,
$\theta_{E}=\sqrt{{4GMD_{ls}}/{c^{2}D_{l}D_{s}}}$. But we note that for some
lenses, such as the SIE (and SIS) that we are considering here,
$\theta_{E}={4\pi v^{2}D_{ls}}/{c^{2}D_{s}}$ is independent of the Hubble
constant. Therefore, one should treat $\theta_{E}$ as just an angular scale of
the lens that one can observe with an EM telescope (González et al., 2020).

To summarize, from lensed GWs, we would obtain $\Delta
t_{d,ij}=d_{L}G\tau_{ij}$ and
$d_{L,i}^{\text{eff}}={d_{L}}/{\sqrt{|\mu_{i}|}}$. Previously, we mentioned we
would be inferring $T_{*}$ and $d_{L}$ in GW lens reconstruction, but here, we
note it being equivalent to inferring the value of $d_{L}$ and $G$ instead.
This would collectively describe the degeneracies outlined in Section 4.1,
especially the partial degeneracy as a result of
$\\{z_{l},z_{s},\theta_{E}\\}$ being coupled in the two separate quantities
$d_{L}(z_{s},H_{0})$ and $G(z_{l},z_{s},\theta_{E})$.

If the gravitational lens system is identified in the EM band, the EM
telescope image data, as well as source and lens redshift information from
photometric/spectroscopic redshift measurement, will yield a measurement for
$z_{s}$, $z_{l}$ and $\theta_{E}$ directly. In principle, it may be possible,
in certain rare scenarios, that the lens system is identified without
identification of the source galaxy. In such a case, $z_{l}$ and $\theta_{E}$
might be known and one would be able to solve for $z_{s}$ from the inferred
value of $G$. That is, $z_{s}$ can be inferred and is not required for
cosmological applications. Indeed, one can solve for $H_{0}$ from the inferred
luminosity distance $d_{L}(z_{s},H_{0})$ together with the time-delay
distance. That is, the source redshift is not strictly speaking required for
cosmological applications.

To illustrate this, we perform a post-lens reconstruction analysis using the
result from the last section (i.e. the posterior
$P(\vec{y},f,T_{*},d_{L}|\vec{d}_{GW})$ showed in Fig. 4, Fig. 5 and Fig. 15),
from which we infer $H_{0}$. We assume the correct lens system is identified
from EM observations, and we further assume the precise value of
$z_{l}=0.5,\theta_{E}\approx 0.592\text{ arcseconds}$ and $z_{s}=2$ are
obtained without uncertainty. Fig. 10 shows the result, in which we recovered
the injected value of the Hubble constant without source redshift. Therefore,
the similarity transformation degeneracy can be broken, and the Hubble
constant can be measured when one can observe the lens and the host galaxy
through another channel.

Figure 10: Result of Hubble constant $H_{0}$ inference using result from the
run shown in Fig. 4, Fig. 5 and Fig. 15, after assuming accurate information
of $z_{l}=0.5,\theta_{E}\approx 0.592\text{ arcseconds}$ and $z_{s}=2$
obtained from EM observation and lens system identification. The parameters
shown: $\\{T_{*}\text{ (in seconds)},d_{L}\text{ (in Mpc)},H_{0}\text{ (in
km}\text{ s}^{-1}\text{ Mpc}^{-1}\text{)}\\}$ are the time delay scaling, the
luminosity distance and Hubble constant respectively. We note that the peak of
the inferred $H_{0}$ overlaps with the injection value of $67.8\text{
km}\text{ s}^{-1}\text{ Mpc}^{-1}$.

In our discussion, we stated that, due to the similarity transformation
degeneracy, we could only perform precise cosmological studies in GW lensing
when complemented with EM information (source and lens redshift, as well as
Einstein radius). But we stress that, it is still possible to constrain
cosmological parameters (such as $H_{0}$) to some degree if a known
statistical distribution of the source and lens redshift ($z_{s}$ and
$z_{l}$), as well as the velocity dispersion $v$ of the lens are applied in
the analysis (such statistical distribution can be obtained by non-lensing
related survey). Such methodology is demonstrated in Jana et al. (2023) (in
which they also studied the constraint on the matter density $\Omega_{m}$). In
our work, we restrict ourselves to the case of precise measurements, where we
do not use the statistical distributions of the source and lens parameters
$z_{s}$, $z_{l}$ and $v$. This explains the apparent discrepancy between our
conclusions.

###  4.5 Analytical Example: SIS Lens

Here, we provide an analytical example to explicitly demonstrate the idea of
dark siren GW lens reconstruction under similarity transformation degeneracy.
Consider a singular isothermal sphere (SIS) lens, with its surface density
given by

$\Sigma(\vec{\xi})=\frac{v^{2}}{2\xi},$ (28)

where $v$ is the velocity dispersion of the lens galaxy and $\vec{\xi}$ is the
(Cartesian) impact parameter on the lens plane (i.e.
$\vec{\xi}=\vec{x}\theta_{E}D_{l}$ in our notation). Its (angular) Einstein
radius is given by $\theta_{E}=4\pi v^{2}{D_{ls}}/{D_{s}}$ and the
dimensionless deflection potential is $\Psi(\vec{x})=\lvert\vec{x}\rvert$
(Takahashi & Nakamura, 2003).

Therefore, the dimensionless lens equation is now

$\vec{y}=\vec{x}-\nabla_{\vec{x}}\lvert\vec{x}\rvert=\vec{x}\left(1-\frac{1}{\lvert\vec{x}\rvert}\right).$
(29)

Since this is a symmetric system, we expect the source and images to lie on
the same line (which goes through the centre of the lens). Hence it becomes a
one dimensional problem: $\vec{y}\rightarrow y$, $\vec{x}\rightarrow x$.
Moreover, we can set $y>0$ without loss of generality. Then, the dimensionless
lens equation becomes

$y=x\left(1-\frac{1}{\lvert x\rvert}\right),$ (30)

which has two solutions for the dimensionless image position when $0<y<1$:

$\begin{split}x_{+}&=y+1,\\\ x_{-}&=y-1.\end{split}$ (31)

We then have the time delay of the $x_{+}$ image as:

$\begin{split}t_{d}(x_{+})&=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}\left[\frac{|y+1-y|^{2}}{2}-(y+1)\right]\\\
&=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}\left(-y-\frac{1}{2}\right),\end{split}$
(32)

and that of the $x_{-}$ image as:

$t_{d}(x_{-})=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}\left(y-\frac{1}{2}\right).$
(33)

Hence, the relative time delay between the two images is

$\Delta
t_{d}=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}(2y)=T_{*}(2y).$
(34)

Meanwhile, the image magnifications are given by

$\mu_{\pm}=\pm 1+{1}/{y}.$ (35)

The corresponding effective luminosity distances are

$d_{L,\pm}^{\text{eff}}={d_{L}}/{\sqrt{\lvert\mu_{\pm}\rvert}}={d_{L}}/{\sqrt{\lvert\pm
1+{1}/{y}\rvert}}.$ (36)

Now, assume we only have dark siren lensed GW observation so that we only
observe $\Delta t_{d}$ and $d_{L,\pm}^{\text{eff}}$. Furthermore, assume the
observables we obtained are certain. One can then solve for the dimensionless
source position $y$ fully in terms of the observables: by realising
${\mu_{+}}/{\mu_{-}}=-\left({d_{L,-}^{\text{eff}}}/{d_{L,+}^{\text{eff}}}\right)^{2}$,
we have

$y=\frac{\mu_{+}/\mu_{-}-1}{\mu_{+}/\mu_{-}+1}=\frac{\left({d_{L,-}^{\text{eff}}}/{d_{L,+}^{\text{eff}}}\right)^{2}+1}{\left({d_{L,-}^{\text{eff}}}/{d_{L,+}^{\text{eff}}}\right)^{2}-1}.$
(37)

Then, from Eq.(34), $T_{*}$ can be calculated from the observables directly as

$T_{*}=\frac{\Delta t_{d}}{2y}=\frac{\Delta
t_{d}}{2}\frac{\left({d_{L,-}^{\text{eff}}}/{d_{L,+}^{\text{eff}}}\right)^{2}-1}{\left({d_{L,-}^{\text{eff}}}/{d_{L,+}^{\text{eff}}}\right)^{2}+1},$
(38)

and $d_{L}$ can be calculated by

$d_{L}=d_{L,\pm}^{\text{eff}}\sqrt{\left|\pm
1+\frac{\left({d_{L,-}^{\text{eff}}}/{d_{L,+}^{\text{eff}}}\right)^{2}-1}{\left({d_{L,-}^{\text{eff}}}/{d_{L,+}^{\text{eff}}}\right)^{2}+1}\right|}.$
(39)

At this stage, we have used up all the information we have, finding the value
of $T_{*}$ but not being able to find the values of its components, including
$\theta_{E}$ (which is then related to the velocity dispersion $v$ of the
galaxy). Therefore, in a more general dark siren SIS lens reconstruction
scenario, we would only be able to infer $P(y,T_{*},d_{L}|\vec{d}_{GW})$,
instead of the full joint posterior distribution
$P(v,z_{l},\beta,z_{s},H_{0}|\vec{d}_{GW})$, which would be the optimal
result.

##  5 Mass-sheet degeneracy

The mass-sheet degeneracy (or magnification degeneracy) was first pointed out
by Falco et al. (1985). It refers to a transformation of the lens system that
leaves observables such as image positions and flux/magnification ratio of
images invariant, except for the quantity $H_{0}t_{d}$. This degeneracy is a
well-known issue in lensing observation and analysis on the EM side, and it
can lead to a biased estimation of the Hubble constant $H_{0}$ from lensing if
unaccounted for (Schneider & Sluse, 2013).

The mass sheet transformation modifies the convergence of lens model as
follows (Schneider & Sluse, 2013):

$\bar{\kappa}(\vec{\theta})=(1-\kappa_{0})\kappa(\vec{\theta})+\kappa_{0},$
(40)

where $\bar{\kappa}$ denotes the transformed convergence of the new lens model
and $\kappa$ denotes the original convergence of the original lens model.
Physically, it can be interpreted as adding a uniform and infinitely-spanning
mass sheet with convergence $\kappa_{0}$ at the lens while reducing the mass
density of the original lens by the factor of $(1-\kappa_{0})$. The presence
of such mass sheet may be due to external convergence, such as the lens being
embedded in low-mass galaxy cluster, or due to structures along the line of
sight (Falco et al., 1985; Schneider & Sluse, 2013; Birrer et al., 2016).
Schneider & Sluse (2013) also demonstrated that the infinitely spanned uniform
mass sheet has very similar effect as a finite mass sheet that is located
within $2\theta_{E}$ from the centre of the lens, hence the mass sheet
transformation can also approximate a similar transformation of the lens model
just near the core region of the lens. Note that the Einstein radius
$\theta_{E}$ is not modified 252525In fact, for a circular lens, the mass-
sheet degeneracy also preserves the total mass within the Einstein radius
(Saha, 2000).. We also note that, the transformed dimensionless deflection
potential, $\bar{\Psi}(\vec{x})$, becomes:

$\bar{\Psi}(\vec{x})=(1-\kappa_{0})\Psi(\vec{x})+\frac{\kappa_{0}}{2}|\vec{x}|^{2}.$
(41)

The transformation also changes the old source position $\vec{\beta}$ to the
new source position $\vec{\bar{\beta}}$ by

$\vec{\bar{\beta}}=(1-\kappa_{0})\vec{\beta}.$ (42)

In dimensionless coordinates, the transformation is simply

$\vec{\bar{y}}=(1-\kappa_{0})\vec{y},$ (43)

where $\vec{\bar{y}}$ is the new dimensionless source position.

Therefore, the transformed lens equation is:

$\vec{\bar{y}}=\vec{x}-\nabla_{\vec{x}}\bar{\Psi}(\vec{x}).$ (44)

Directly substituting the transformation above, we have the following steps:

$\displaystyle(1-\kappa_{0})\vec{y}$
$\displaystyle=\vec{x}-\nabla_{\vec{x}}\left[(1-\kappa_{0})\Psi(\vec{x})+\frac{\kappa_{0}}{2}|\vec{x}|^{2}\right]$
$\displaystyle=\vec{x}-(1-\kappa_{0})\nabla_{\vec{x}}\Psi(\vec{x})-\frac{\kappa_{0}}{2}\nabla_{\vec{x}}|\vec{x}|^{2}$
$\displaystyle=\vec{x}-(1-\kappa_{0})\nabla_{\vec{x}}\Psi(\vec{x})-\kappa_{0}\vec{x}$
$\displaystyle=(1-\kappa_{0})\vec{x}-(1-\kappa_{0})\nabla_{\vec{x}}\Psi(\vec{x}),$

which, after cancelling $(1-\kappa_{0})$ from both sides, reduces to the
original lens equation of the system before mass sheet transformation:
$\vec{y}=\vec{x}-\nabla_{\vec{x}}\Psi(\vec{x})$. Hence, the solution of
dimensionless image position $\vec{x}$ remains the same, and since the
Einstein radius is not changed, the actual image positions of the lens system
are not altered by the mass sheet transformation. Therefore, it is impossible
to lift the degeneracy from image position alone, even in the case of EM
lensing observation (Congdon & Keeton, 2018).

Under the mass sheet transformation, the new image magnification is

$\bar{\mu}_{i}=\frac{\mu_{i}}{(1-\kappa_{0})^{2}},$ (45)

where $\mu$ is the original magnification before mass sheet transformation and
the subscript $i$ denotes the $i$-th image. We note that, since all image
magnifications are affected by the same factor of $(1-\kappa_{0})^{-2}$, the
magnification ratios are preserved under mass sheet transformation

$\frac{\bar{\mu}_{i}}{\bar{\mu}_{j}}=\frac{\mu_{i}}{\mu_{j}}.$ (46)

Hence, in the case of EM lensing observation, as long as the actual source
flux is unknown and we only have access to flux ratio between images, we
cannot lift the mass-sheet degeneracy.

For GW lensing, since it is a standard siren, in the case of no lens and
source redshift measurement, after mass sheet transformation, we will be
measuring effective luminosity distance of the images:

$\bar{d}^{\text{eff}}_{L,i}=\frac{d_{L}}{\sqrt{|\bar{\mu}_{i}|}}=\frac{d_{L}(1-\kappa_{0})}{\sqrt{|\mu_{i}|}}=\frac{d_{L}^{\kappa_{0}}}{\sqrt{|\mu_{i}|}}$
(47)

where we define $d_{L}^{\kappa_{0}}\equiv(1-\kappa_{0})d_{L}$ as the mass-
sheet-degenerated luminosity distance. We also note that the Morse index of
the images remain unchanged.

Moreover, the new relative time delay is

$\Delta\bar{t}_{d,ij}=(1-\kappa_{0})\Delta t_{d,ij},$ (48)

where $\Delta t_{d,ij}=T_{*}\tau_{ij}$ is the original relative time delay
between image $i$ and $j$, and again, we emphasize that $\tau_{ij}$ is
difference of the dimensionless time delay for image pair $i,j$ of the
original system. Hence, we can define the mass-sheet-degenerated time delay
scaling as

$T_{*}^{\kappa_{0}}\equiv(1-\kappa_{0})T_{*}$ (49)

We can then write the new relative time delay as

$\Delta\bar{t}_{d,ij}=T_{*}^{\kappa_{0}}\tau_{ij}.$ (50)

In addition, we note that, similar to the invariant magnification ratio of
different images, the relative time delay ratios of different image pairs:

$\frac{\Delta\bar{t}_{d,ij}}{\Delta\bar{t}_{d,kl}}=\frac{\tau_{ij}}{\tau_{kl}},$
(51)

are also invariant under mass sheet transformation due to the common factor of
$T_{*}^{\kappa_{0}}$ on all time delays. Fig. 11 summarizes the physical
picture and the invariant observables of a multi-messenger lens system under
mass sheet transformation.

Figure 11: Illustrative example of mass-sheet degeneracy in multi-messenger
lensing: the upper and lower lens system on the left panel corresponds to a
lens system before and after mass sheet transformation respectively. We note
that the source emits both EM waves and GW waves (for instance, a binary
neutron star merger, or a binary black hole merger and its host galaxy). The
middle panel illustrates the image position and magnification observed through
EM telescope, and the right panel illustrates the relative time delay between
images observed by GW detectors. Due to mass-sheet degeneracy, the two systems
lead to the same image positions, flux ratio and relative time delay ratio of
different image pairs.

###  5.1 Lens Reconstruction with Similarity Transformation and mass-sheet
degeneracy

Here, for dark siren GW lensing, we can only measure relative time delays
$\Delta\bar{t}_{d,ij}$ and effective luminosity distances
$\bar{d}^{\text{eff}}_{L,i}$ of the images. In lens reconstruction, since
$T_{*}^{\kappa_{0}}$ and $d_{L}^{\kappa_{0}}$ are common factor of all the
images’ relative time delay and effective luminosity distance respectively,
similar to the case of similarity transformation degeneracy, we can eliminate
them and obtain the magnification ratios as

$\frac{\mu_{i}}{\mu_{j}}=\left(\frac{\bar{d}_{L,j}^{\text{eff}}}{\bar{d}_{L,i}^{\text{eff}}}\right)^{2}(-1)^{\Delta
n_{ij}},$ (52)

and the dimensionless time delay ratios as

$\frac{\tau_{ij}}{\tau_{mn}}=\frac{\Delta\bar{t}_{d,ij}}{\Delta\bar{t}_{d,mn}}.$
(53)

From these ratios, we can only reconstruct the dimensionless system that does
not have mass sheet transformation taken into account, namely the original
(pre-mass-sheet-transformation) lens profile $\Psi(\vec{x})$ (for instance,
the axis ratio $f$ in SIE lens) and also the original dimensionless source
position $\vec{y}$. Then, we can only infer the values of the mass-sheet-
degenerated time delay scaling $T_{*}^{\kappa_{0}}$ and the mass-sheet-
degenerated luminosity distance $d_{L}^{\kappa_{0}}$, and we emphasize that we
cannot resolve any terms inside of them. In summary, we can only reconstruct
$P(\vec{y},\Psi(\vec{x}),T_{*}^{\kappa_{0}},d_{L}^{\kappa_{0}}|\vec{d}_{GW})$
with dark siren GW lensing.

Type of parameters |  $\vec{\Phi}_{l}$ (SIE lens with mass sheet) |  $\vec{\Phi}_{s}$ (GW source) |  $\vec{\Phi}_{c}$ (cosmology)
---|---|---|---
Injection values | $z_{l}=0.5,v=0.0006c,f=0.5,\kappa_{0}=0.05$ | $z_{s}=2,\vec{\beta}\approx(0.00592,0.00887)\text{ (in arcseconds)}$ | $H_{0}=67.8\text{ km}\text{ s}^{-1}\text{ Mpc}^{-1}$
Table 5: Injection values of the system parameters for our lens reconstruction example under mass sheet transformation: 1. for the lens parameters, $z_{l},v,f,\kappa_{0}$ are the redshift, velocity dispersion, SIE axis ratio and mass sheet convergence respectively. Note that $c$ is the speed of light and we align the long axis of the SIE with the vertical axis (i.e. the $y_{2}$ or $x_{2}$ axis). 2. For the source, $z_{s},\vec{\beta}$ are the redshift and the (angular) source position in arcseconds (before mass sheet transformation is applied) respectively. 3. For cosmology part, $H_{0}$ is the Hubble constant. Observables | Effective luminosity distance: $d_{L,1}^{\text{eff}},d_{L,2}^{\text{eff}},d_{L,3}^{\text{eff}},d_{L,4}^{\text{eff}}$ (in Mpc) | Relative time delay between images: $\Delta t_{d,12},\Delta t_{d,23},\Delta t_{d,34}$ (in s)
---|---|---
Injected values | $8934,8775,11966,12510$ | $-53026,-467148,-63264$
Table 6: GW lensing observables of injected system setup shown in Table 5.
Note that we have four images and the image types are I-I-II-II. Figure 12:
Result of SIE lens reconstruction of dark siren lensed GW using injected
parameters shown in Table 5 and Table 6. After taking similarity
transformation degeneracy and mass-sheet degeneracy into account, the five
free parameters are
$\\{\vec{y}=(y_{1},y_{2}),f,T_{*}^{\kappa_{0}},d_{L}^{\kappa_{0}}\\}$, which
are the dimensionless source position (before mass sheet transformation), the
SIE axis ratio, the mass-sheet-degenerated time delay scaling
$T_{*}^{\kappa_{0}}\equiv(1-\kappa_{0})\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}\theta_{E}^{2}$
and the mass-sheet-degenerated luminosity distance respectively. Here, we show
the posterior distribution for the parameters that describe the scale of the
system: $T_{*}^{\kappa_{0}}\text{ (in seconds)},d_{L}^{\kappa_{0}}\text{ (in
Mpc)}$. We can see, the two parameters are recovered fairly well without
degeneracy.

Here we show an example of dark siren GW lens reconstruction under both
similarity transformation degeneracy and mass-sheet degeneracy. Similar to
Section 4, we inject an SIE lens with the system parameters shown in Table 5,
which gives injected observables presented in Table 6. Same as the likelihood
$P(\vec{d}_{GW}|\vec{\nu}_{im})$ described in Section 4.2, we assume an
uncorrelated 5% Gaussian error on the measured $d_{L,j}^{\text{eff}}$ and
$\Delta t_{d,ij}$. Fig. 12 shows the result of lens reconstruction of the free
parameters $\\{T_{*}^{\kappa_{0}},d_{L}^{\kappa_{0}}\\}$ (the full result of
parameter estimation of the parameters
$\\{\vec{y},f,T_{*}^{\kappa_{0}},d_{L}^{\kappa_{0}}\\}$ is shown in Fig. 20 in
the Appendix), and we are able to recover them accurately. However, it is not
possible to distinguish the effect of mass sheet convergence $\kappa_{0}$ from
the absolute time-delay scaling $T_{*}$ or luminosity distance $d_{L}$. We
note that, in Section III A of Cremonese et al. (2021), they also demonstrated
the scaling relation of mass sheet density and the GW amplitude (which is
related to the effective luminosity distance defined in Eq. (47)), as well as
mass sheet degeneracy itself from the perspective of GW waveform in the
geometric optics.

###  5.2 Inability to Resolve mass-sheet degeneracy Even with a Telescope
Image and Spectroscopic Redshifts

Similar to Section 4.4, we assume we obtain the source and lens redshift
($z_{s}$ and $z_{l}$) through identifying the corresponding lens system in EM
observations. We furthermore assume that the Einstein radius $\theta_{E}$ can
be obtained from a high-resolution image of the lens. Since GW is a standard
siren, it may intuitively make sense that one should be able to break the
mass-sheet degeneracy and measure both $\kappa_{0}$ and $H_{0}$ at the same
time. Unfortunately, this is not the case. Under the mass sheet
transformation, the relative time delay between images are given by

$\begin{split}\Delta\bar{t}_{d,ij}&=(1-\kappa_{0})T_{*}\tau_{ij}=\left[(1-\kappa_{0})d_{L}\right]\left[\frac{D_{l}}{D_{ls}}\frac{1+z_{l}}{c(1+z_{s})^{2}}\theta_{E}^{2}\right]\tau_{ij}\\\
&=d_{L}^{\kappa_{0}}G\tau_{ij},\end{split}$ (54)

where
$d_{L}^{\kappa_{0}}=d_{L}^{\kappa_{0}}(z_{s},H_{0},\kappa_{0})\equiv(1-\kappa_{0})d_{L}(z_{s},H_{0})$
and $G(z_{l},z_{s},\theta_{E})$ is defined in Eq. (27). In a flat universe, we
have

$d_{L}^{\kappa_{0}}=\frac{c(1-\kappa_{0})}{H_{0}}\int_{0}^{z_{s}}\frac{dz^{\prime}}{\sqrt{\Omega_{r}(1+z^{\prime})^{4}+\Omega_{m}(1+z^{\prime})^{3}+\Omega_{\Lambda}}}.$
(55)

Thus, after mass sheet transformation,
$\bar{d}^{\text{eff}}_{L,i}=\frac{d_{L}^{\kappa_{0}}}{\sqrt{|\mu_{i}|}}$ and
$\Delta\bar{t}_{d,ij}=d_{L}^{\kappa_{0}}G\tau_{ij}$. Recall that for the
absolute scale of the lens system, we can only infer the values of the pair
$(T_{*}^{\kappa_{0}},d_{L}^{\kappa_{0}})$. Given the relationship between the
luminosity distance and the angular diameter distance to the source, we can
conclude that, effectively or equivalently, one can only infer the values of
the pair $d_{L}^{\kappa_{0}}(z_{s},H_{0},\kappa_{0})$ and
$G(z_{l},z_{s},\theta_{E})$.

As mentioned in Section 4.4, with EM information, we will have $z_{s}$,
$z_{l}$ and $\theta_{E}$ information. Under fixed cosmology, we can infer
$G(z_{l},z_{s},\theta_{E})$. However, after we infer the value of
$d_{L}^{\kappa_{0}}(z_{s},H_{0},\kappa_{0})$, even with knowledge of $z_{s}$,
we still have a coupled pair of unknown $(H_{0},\kappa_{0})$ hidden in
$d_{L}^{\kappa_{0}}$ in the form of ${(1-\kappa_{0})}/{H_{0}}$ (Eq. (55)),
making the two fully degenerate. The reason behind is the effect of mass sheet
transformation modifies the effective luminosity distance and time delay in
the same way, namely multiplying the factor $(1-\kappa_{0})$. Therefore,
lensed GWs with EM counterpart (or observation of the lens and the lensed host
galaxy) will not be able to break mass-sheet degeneracy with just lensing
observables.

Fig. 13 illustrates parameter estimation attempts trying to infer $\kappa_{0}$
and $H_{0}$ with a known (fixed) $z_{s}$, $z_{l}$ and $\theta_{E}$. Both
panels have the same setup as the one described in Tables 5 and 6. We assume
the same uniform prior for $H_{0}$. The difference between the two panels is
that the prior of $\kappa_{0}$ in the simulation shown in the upper panel is
set to be uniform from $0$ to $0.9$, while the simulation shown in the lower
panel is set to be uniform from $-0.9$ to $0.9$.262626This is unphysical since
we include the possibility of negative mass sheet, but mathematically, it is
still valid, and we intend to use this prior to check the validity of the
result shown in the upper panel. For both results, we can see there is a
$H_{0}$-$\kappa_{0}$ degeneracy spanning across the parameter space.

In the upper panel, it may look as if the peak of the posterior of $H_{0}$
still roughly overlaps with the injected values. However, it is just a
coincidence due to assumption of priors for $H_{0}$ and $\kappa_{0}$.When we
assume a different prior on $\kappa_{0}$, as shown in the lower panel, the
posterior peak of $H_{0}$ is far off from injection, and the posterior
probability of $\kappa_{0}$ drops as $\kappa_{0}$ increases, even in region
where $\kappa_{0}$ is negative. It becomes clear that the inferred posterior
distribution of $\kappa_{0}$ is not informative and affected by prior choice.
These results indicate that constraining both $\kappa_{0}$ and $H_{0}$
simultaneously from GWs strongly lensed by galaxy is not possible even when an
EM lens reconstruction is available, due to the mass-sheet degeneracy.

Figure 13: Result of the Hubble constant $H_{0}$ and mass sheet convergence
$\kappa_{0}$ inference for lensed GWs. Injection model and parameters shown in
Table 5 and Table 6 are used. In the run presented in the upper panel, the
prior of $\kappa_{0}$ is set to be uniform from $0$ to $0.9$, while that of
the lower panel is set to be uniform from $-0.9$ to $0.9$. We assume $z_{s}$,
$z_{l}$ and $\theta_{E}$ are known from identifying the corresponding lens
system in EM telescope observation. As we can see, the Hubble constant $H_{0}$
and mass sheet convergence $\kappa_{0}$ are poorly constrained, degenerate,
and sensitive to the choice of prior.

If we perform Hubble constant measurement ignoring the possible presence of a
mass-sheet (as demonstrated in Section 4.4), our measured Hubble constant
would be biased from the true Hubble constant by a factor of
$\frac{1}{1-\kappa_{0}}$. Such issues in Hubble constant measurement via
observing lensing of EM waves are further discussed in Schneider & Sluse
(2013).

In addition, the above analysis can be applied to similar situations in
lensing of light. For instance, for time varying light sources that may not be
standard candle, such as quasars, if one assume a fixed cosmology and known
$H_{0}$, mass-sheet degeneracy can be broken and the mass sheet density
$\kappa_{0}$ can be found as follows. Despite there will not be measurement of
the effective luminosity distance $\bar{d}_{L}^{\text{eff}}$ (as the total
flux is not known), from the image position and time delay measurement from
the EM telescope, one should be able to measure the Einstein radius
$\theta_{E}$, as well as reconstruct the dimensionless source position and
calculate the dimensionless time delay $\tau_{ij}$ between images. Together
with redshift measurement of the lens and source, $G(z_{s},H_{0},\theta_{E})$
can be calculated, hence the value of $d_{L}^{\kappa_{0}}$ can be calculated
with eq. (54). As cosmology is fixed and $z_{s}$ is known, $\kappa_{0}$ can be
found from eq. (55) easily. This would be helpful if one would like to study
the lens galaxy, such as it being possibly embedded in a galaxy cluster or a
possible finite mass sheet located near the centre of the lens that can mimic
the mass-sheet transformation (Falco et al., 1985; Schneider & Sluse, 2013;
Birrer et al., 2016), without standard candle sources (but an accurate value
of $H_{0}$ from other measurement is required).

##  6 Science applications and implications of the two degeneracies

###  6.1 Test of Modified Gravitational Wave Propagation

Some theories modifying general relativity on cosmological scales have as
typical signature a modified GW dispersion relation (propagation effect)
(Deffayet & Menou, 2007; Belgacem et al., 2018a, 2019; Lagos et al., 2019;
Finke et al., 2021). Moreover, modified GW propagation may be one of the most
promising window to observationally test such theories (Mirshekari et al.,
2012; Finke et al., 2021). The idea of testing modified GW propagation effect
on the luminosity distance of GW (which we will refer as GW luminosity
distance, $d_{L}^{\text{gw}}$, below) by lensing was previously discussed in
Finke et al. (2021); Iacovelli et al. (2022); Narola et al. (2023), on which
the results presented in this section build. In addition, here we ask the
question: "Given that the proposed modified GW propagation tests rely on
luminosity distance measurements, could the mass-sheet degeneracy alter their
results?" The answer appears to be no.

The relationship between GW and EM luminosity distance $d_{L}^{\text{em}}$
(which is just the original/standard luminosity distance as defined in Eq.
(18)) is written as

$d_{L}^{\text{gw}}=\Xi(z_{s})d_{L}^{\text{em}},$ (56)

where for a wide range of modified gravity theories, the GW to EM luminosity
distance ratio $\Xi(z_{s})$ can be parametrized as in Belgacem et al. (2018b):

$\Xi(z_{s})=\Xi_{0}+\frac{1-\Xi_{0}}{(1+z_{s})^{n}},$ (57)

where the parameters $\Xi_{0}$ and $n$ are determined by the specific theory
(Belgacem et al., 2019).

We now include mass-sheet degeneracy in our analysis, using same convention as
Section 5. The measured effective luminosity distance for the $i$-th GW lensed
image is given by

$\bar{d}_{L,i}^{\text{eff}}=\frac{d_{L}^{\text{gw}}(1-\kappa_{0})}{\sqrt{|\mu_{i}|}}=\frac{d_{L}^{\text{gw},\kappa_{0}}}{\sqrt{|\mu_{i}|}},$
(58)

where we defined the mass-sheet-degenerated GW luminosity distance as
$d_{L}^{\text{gw},\kappa_{0}}\equiv
d_{L}^{\text{gw}}(1-\kappa_{0})=\Xi(z_{s})d_{L}^{em}(1-\kappa_{0})$.

Meanwhile, the relative time delay between images is not affected, which means
we still have $\Delta\bar{t}_{d,ij}=T_{*}^{\kappa_{0}}\tau_{ij}$ as before.
Since the angular diameter distance to source and EM luminosity distance
relation is $D_{s}={d_{L}^{\text{em}}}/{(1+z_{s})^{2}}$, we have

$T_{*}^{\kappa_{0}}=\left[d_{L}^{\text{em}}(1-\kappa_{0})\right]G(z_{l},z_{s},\theta_{E}).$
(59)

Therefore, in dark siren GW lens reconstruction, besides from all the pre-mass
sheet transformation dimensionless quantities of the lens system (i.e. the
dimensionless source position $\vec{y}$ and the SIE axis ratio $f$ in our
example), we can only infer the values of $d_{L}^{\text{gw},\kappa_{0}}$ and
$T_{*}^{\kappa_{0}}$. Assume $z_{l}$, $z_{s}$ and $\theta_{E}$ are known from
EM observation after identifying the lens system, we can calculate the value
of $G(z_{l},z_{s},\theta_{E})$ 272727We note that unlike in Section 4.4 or 5.2
where GW and EM luminosity distance are the same, $G(z_{l},z_{s},\theta_{E})$
here cannot be inferred directly., and the GW to EM luminosity distance ratio
can then be calculated by

$\Xi(z_{s})=\frac{d_{L}^{\text{gw},\kappa_{0}}}{T_{*}^{\kappa_{0}}}G(z_{l},z_{s},\theta_{E}),$
(60)

which means that we do not even need to have knowledge on the Hubble constant
$H_{0}$ or the mass sheet convergence $\kappa_{0}$. We only require a good
lens reconstruction process that infer accurate values of
$d_{L}^{\text{gw},\kappa_{0}}$ and $T_{*}^{\kappa_{0}}$, as well as $z_{l}$,
$z_{s}$ and $\theta_{E}$ information of the system. In other words, contrary
to a naive expectation, this type of modified gravity test is immune to mass-
sheet degeneracy. It can be understood as a result of the modified GW
propagation only affecting the effective luminosity distance measurement (Eq.
(58)) while not affecting the angular diameter distances and hence the time
delay scaling (Eq. (59)) and time delay measurement.

Using the injection parameters given in Table 5, also injecting
$n=2,\Xi_{0}=1.8$, which leads to $\Xi(z_{s})=1.71$, we show the recovery of
$\Xi(z_{s})$ for an SIE model under unknown mass sheet transformation. The
full lens reconstruction result of the free parameters
$\\{\vec{y}=(y_{1},y_{2}),f,T_{*}^{\kappa_{0}},d_{L}^{\text{gw},\kappa_{0}}\\}$
is shown in Fig. 21 in the Appendix. Assuming $z_{s},z_{l}$ and $\theta_{E}$
are known from EM observation, we can infer the value of $\Xi(z_{s})$ using
Eq. (60). The recovered posterior distribution of $T_{*}^{\kappa_{0}}$ and
$d_{L}^{\text{gw},\kappa_{0}}$ is shown in Fig. 14.

Figure 14: Result of $\Xi(z_{s})$ inference from lensed GW under injected
parameters shown in Table 5 (with the addition of $n=2,\Xi_{0}=1.8$). After a
posterior distribution of $T_{*}^{\kappa_{0}}$ and
$d_{L}^{\text{gw},\kappa_{0}}$ is obtained in the lens reconstruction run
shown in Fig. 21 in the Appendix, assuming $z_{s},z_{l}$ and $\theta_{E}$ are
known from EM observation and lens system identification, Eq. (60) is used to
infer $\Xi(z_{s})$. The peak of the posterior overlaps with the injection
value of $\Xi(z_{s})=1.71$.

###  6.2 Implications on Selected Other Applications of Strongly Lensed
Gravitational Waves

As demonstrated, the modified GW propagation test is affected by the
similarity transformation degeneracy but not the mass-sheet degeneracy. In
short, the two degeneracies affect different strong lensing applications
differently. Studies that do not require lens reconstruction should not be
affected by both degeneracies. Two of such examples are testing speed of GW
against light (Collett & Bacon, 2017; Fan et al., 2017) and associating fast
radio burst sources with GW sources (Singh et al., 2023): they are model
independent and only require comparing the time delay between GWs and EM
counterpart. In addition, GW polarization tests are also unaffected, as they
rely on comparing the polarization content of the lensed GW signals that are
each detected by a different detector configuration (Smith et al., 2017; Goyal
et al., 2021a; Hernandez, 2022) and do not require lens modelling.

Meanwhile, cosmography studies with strongly lensed GWs, such as Hubble
constant measurement as demonstrated above are affected by the similarity
transformation degeneracy, as they generally require lens reconstruction to
single out the effect induced by lensing and a further decoupling of the
astronomical quantities (such as redshifts and the angular diameter distances)
inside the time delay scaling $T_{*}$ and luminosity distance $d_{L}$.
Specific applications such as Hubble constant measurement would suffer from
mass-sheet degeneracy as well. Non-lensing related means from EM telescope
observation, for instance, independent estimates of the lens mass using
spectroscopic measurement on velocity dispersion (Schneider & Sluse, 2013),
are required to break such degeneracy. However, Cremonese et al. (2021); Chen
et al. (in prep.) showed that if GW lensing is in the wave optics regime, it
would be possible to break the mass-sheet degeneracy from GW observation
alone.

##  7 Summary and Outlook

| Can we measure …? | Can we break … degeneracy? | Can we test …?
---|---|---|---
Observations | $\vec{y}$ | $f$ | $T_{*}$ | $d_{L}$ | $\theta_{E}$ | Similarity transformation | Mass-sheet | Modified GW propagation
2 GW images | ✗ | ✗ | ✗ | ✗ | ✗ | ✗ | ✗ | ✗
3 GW images (correct ordering) | ✓* | ✓* | ✓* | ✓* | ✗ | ✗ | ✗ | ✗
4 GW images | ✓ | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | ✗
4 GW images & EM obs (w/o $z_{l},z_{s}$) | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✗ | ✗
4 GW images & EM obs (w/ $z_{l}$, w/o $z_{s}$) | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✗
4 GW images & EM obs (w/ $z_{l},z_{s}$) | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✓
4 GW images & EM obs (w/ $z_{l},z_{s}$ & $v$) | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
Table 7: Summary table of the measurability of SIE lens system parameters (the
dimensionless observables position $\vec{y}$, the axis ratio $f$, the time
delay scaling $T_{*}$, the luminosity distance $d_{L}$ and the angular
Einstein radius $\theta_{E}$), can the similarity transformation degeneracy
and the mass-sheet ($\kappa_{0}-H_{0}$) degeneracy be broken, and can the
modified GW propagation test be done under seven different observation
scenarios, which from top to bottom are: 1. when only two of the four lensed
GWs are observed; 2. when only three of the four lensed GWs are observed, with
the ordering of the three images identified or assumed correctly within the
four; 3. all four lensed GWs are observed without complementary information
from EM telescope (dark siren); 4. all four lensed GWs are observed and the
lens system is identified in EM observation, but redshifts of lens and source
are both unavailable; 5. all four lensed GWs are observed and the lens system
is identified in EM observation, only redshift of lens is available; 6. all
four lensed GWs are observed and the lens system is identified in EM
observation, with redshifts of both lens and source available; 7. all four
lensed GWs are observed and the lens system is identified in EM observation,
with redshifts of both lens and source available, and the lens galaxy’s
velocity dispersion can be measured.
✓: can measure/break the degeneracy/test.
*: can measure, but with large uncertainty.
✗: cannot measure/break the degeneracy/test.  | Can we break … degeneracy?
---|---
Observations (all four images observed) | Similarity transformation | Mass-sheet
Dark standard siren (GW) | ✗ | ✗
Bright standard siren (GW & EM obs w/ $z_{l},z_{s}$, w/o $v$) | ✓ | ✗
Bright standard siren (GW & EM obs w/ $z_{l},z_{s},v$) | ✓ | ✓
Standard Candle (w/ $z_{l},z_{s}$, w/o $v$) | ✓ | ✗
Standard Candle (w/ $z_{l},z_{s},v$) | ✓ | ✓
Time-varying, non-standard candle light source (w/ $z_{l},z_{s}$, w/o $v$) | ✓ | ✗
Time-varying, non-standard candle light source (w/ $z_{l},z_{s},v$) | ✓ | ✓
Non-time-varying, non-standard candle light source (w/ $z_{l},z_{s},v$) | ✓* | ✗
Table 8: Summary table of can the similarity transformation degeneracy and the
mass-sheet degeneracy be broken under eight different lensing of GW and/or
lensing of light observation scenarios (which we assumed a SIE lens system,
with all four images are observed for all cases). From top to bottom are: 1.
dark siren lensed GWs are observed; 2. lensed GWs observation with
complementary EM information (redshift of source and lens) (i.e. bright
siren); 3. lensed GWs are observed with complementary EM information (redshift
of source and lens), with the addition of lens galaxy velocity dispersion $v$
measurement; 4. lensed standard candle (such as supernovae) are observed, with
redshift measurement of source and lens; 5. lensed standard candle (such as
supernovae) are observed, with measurement of redshift of source and lens, as
well as velocity dispersion $v$ of lens galaxy; 6. lensing of a time-varying,
but non-standard candle light source is observed, with redshift measurement of
source and lens; 7. lensing of a time-varying, but non-standard candle light
source is observed, with measurement of redshift of source and lens, as well
as velocity dispersion $v$ of lens galaxy; 8. lensing of a non-time-varying
and non-standard candle light source is observed, with measurement of redshift
of source and lens, as well as velocity dispersion $v$ of lens galaxy.
✓: can break the degeneracy.
*: can break the degeneracy, but $H_{0}$ cannot be measured even neglecting mass sheet.
✗: cannot measure/break the degeneracy.

We have presented a framework for performing lens reconstruction based on dark
siren strongly lensed GWs using a Bayesian analysis framework and
"irreducible" lensing parameters. GW-only lensing will always suffer from two
major degeneracies, which we have detailed in this work. The first one is the
similarity transformation degeneracy, which is a degeneracy between the lens
and source redshift, as well as Einstein radius of the lens and cosmological
parameters such as the Hubble constant. The second one is the mass-sheet
degeneracy, which can alter the lens model and together with modification on
the Hubble constant, they can keep the lensing observables unchanged. Not
accounting for the mass-sheet degeneracy will lead to biased $H_{0}$
measurements.

The similarity transformation degeneracy can be broken once we have a EM
counterpart of the GW signal or we localize and identify the corresponding
lens system (when there is no EM counterpart from the GW event) as
demonstrated in Hannuksela et al. (2020). Unfortunately, the mass-sheet
degeneracy will not be broken in our strong lensing case even if we have the
redshift information and the Einstein radius from high-resolution telescope
images. Instead, other observation such as spectroscopic measurement on lens
galaxy’s velocity dispersion is needed to break the degeneracy (Schneider &
Sluse, 2013).

Furthermore, we have also reviewed two GW lensing science applications and how
they are affected by these degeneracies. Table 7 summarizes the quantities
that we can measure (in an SIE lens system), the degeneracies that we can
break, and if modified GW dispersion relations can be tested for different
observation scenarios, such as dark siren GW observation and GW observation
complemented by redshift and lens galaxy’s velocity dispersion measurement. In
Table 8, we also summarized whether similarity transformation and mass-sheet
degeneracy can be broken under eight different observational scenario of GW
and/or EM lensing. In addition, for Hubble constant measurements, one
advantage of GW strong lensing is that one can take advantage of the
millisecond accuracy of arrival time measurements afforded by the detectors
Fan et al. (2017), which can complement the more traditional time-delay
measurements (Tewes, M. et al., 2013). For this Hubble constant measurement,
one needs to have knowledge of the lens redshift and size/Einstein radius, but
the source redshift is, intriguingly, not required.

In addition, we also summarised the resolvability of the similarity
transformation and mass-sheet degeneracy under different scenarios of
observation of GW and/or light sources (such as standard candle, non-standard
candle but time-varying light sources) in Table LABEL:table:_summary_table_2..

Secondly, we also showed the modified GW propagation test is immune to mass-
sheet degeneracy (and does not require knowledge on Hubble constant), which
can be beneficial for testing theories modifying general relativity.

Besides, understanding these GW-based lens reconstructions and their
degeneracies will be important also for other science applications, for
instance, broader cosmography study. Understanding the galaxy modelling
capabilities will also play a role in understanding how GWs will help probe
dark matter and substructure within galaxy lenses. Indeed, we hope that our
work will help motivate further studies into complementary avenues in
multimessenger lensing, where GWs and EM can provide complementary information
that will help in understanding of strong lenses and forming a rich scientific
case.

##  8 Acknowledgements

We thank Hemantakumar Phurailatpam, Laura Uronen, Leo Ng for valuable
discussions, feedback and comments. Jason S.C. Poon is supported by the Hong
Kong PhD Fellowship Scheme (HKPFS) from the Hong Kong Research Grants Council
(RGC). Jason S.C. Poon and Otto A. Hannuksela acknowledge support by grants
from the ResearchGrants Council of Hong Kong (Project No. CUHK 14304622 and
14307923), the start-up grant from the Chinese University of Hong Kong, and
the Direct Grant for Research from the Research Committee of The Chinese
University of Hong Kong. Stefano Rinaldi acknowledges financial support from
the European Research Council for the ERC Consolidator grant DEMOBLACK, under
contract no. 770017 and from the EU Horizon 2020 Research and Innovation
Program under the Marie Sklodowska-Curie Grant Agreement No. 734303. Justin
Janquart and Harsh Narola are supported by the research programme of the
Netherlands Organisation for Scientific Research (NWO).

##  9 Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding authors.

## References

  * Aasi et al. (2015) Aasi J., et al., 2015, Classical and Quantum Gravity, 32, 074001
  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2019a) Abbott B., et al., 2019a, Physical Review D, 100
  * Abbott et al. (2019b) Abbott B., et al., 2019b, Physical Review Letters, 123
  * Abbott et al. (2021a) Abbott B. P., et al., 2021a, The Astrophysical Journal, 909, 218
  * Abbott et al. (2021b) Abbott R., et al., 2021b, ApJ, 923, 14
  * Abbott et al. (2022) Abbott R., et al., 2022, Physical Review D, 106
  * Abbott et al. (2023) Abbott R., Abe H., Acernese F., Ackley K., Adhicary S., Adhikari N., Adhikari R. X., et al. 2023, arXiv e-prints, p. arXiv:2304.08393
  * Acernese et al. (2014) Acernese F., et al., 2014, Classical and Quantum Gravity, 32, 024001
  * Aghanim et al. (2020) Aghanim N., et al., 2020, Astronomy & Astrophysics, 641, A6
  * Akutsu et al. (2020) Akutsu T., et al., 2020, Progress of Theoretical and Experimental Physics, 2021, 05A101
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, Laser Interferometer Space Antenna (arXiv:1702.00786)
  * Ashton et al. (2019) Ashton G., et al., 2019, The Astrophysical Journal Supplement Series, 241, 27
  * Aso et al. (2013) Aso Y., Michimura Y., Somiya K., Ando M., Miyakawa O., Sekiguchi T., Tatsumi D., Yamamoto H., 2013, Phys. Rev. D, 88, 043007
  * Bagnasco (2023) Bagnasco S., 2023, The Ligo-Virgo-KAGRA Computing Infrastructure for Gravitational-wave Research (arXiv:2311.12559)
  * Bailes et al. (2021) Bailes M., et al., 2021, Nature Reviews Physics, 3, 344
  * Baiotti (2019) Baiotti L., 2019, Progress in Particle and Nuclear Physics, 109, 103714
  * Bartelmann (2010) Bartelmann M., 2010, Classical and Quantum Gravity, 27, 233001
  * Basak et al. (2022) Basak S., Ganguly A., Haris K., Kapadia S., Mehta A. K., Ajith P., 2022, The Astrophysical Journal Letters, 926, L28
  * Belgacem et al. (2018a) Belgacem E., Dirian Y., Foffa S., Maggiore M., 2018a, Phys. Rev. D, 98, 023510
  * Belgacem et al. (2018b) Belgacem E., Dirian Y., Foffa S., Maggiore M., 2018b, Physical Review D, 98
  * Belgacem et al. (2019) Belgacem E., et al., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 024
  * Berti et al. (2015) Berti E., et al., 2015, Classical and Quantum Gravity, 32, 243001
  * Bianconi et al. (2023) Bianconi M., et al., 2023, Monthly Notices of the Royal Astronomical Society, 521, 3421–3430
  * Birrer & Amara (2018) Birrer S., Amara A., 2018, Physics of the Dark Universe, 22, 189–201
  * Birrer et al. (2016) Birrer S., Amara A., Refregier A., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 020
  * Buchner et al. (2014) Buchner J., et al., 2014, Astronomy & Astrophysics, 564, A125
  * Bulashenko & Ubach (2022) Bulashenko O., Ubach H., 2022, Journal of Cosmology and Astroparticle Physics, 2022, 022
  * Cahillane & Mansell (2022) Cahillane C., Mansell G., 2022, Galaxies, 10
  * Cao et al. (2022) Cao M.-D., Zheng J., Qi J.-Z., Zhang X., Zhu Z.-H., 2022, The Astrophysical Journal, 934, 108
  * Chakraborty & Mukherjee (2024) Chakraborty A., Mukherjee S., 2024, GLANCE – Gravitational Lensing Authenticator using Non-Modelled Cross-Correlation Exploration of Gravitational Wave Signals (arXiv:2403.03982)
  * Cheung et al. (2021) Cheung M. H., Gais J., Hannuksela O. A., Li T. G., 2021, Monthly Notices of the Royal Astronomical Society, 503, 3326
  * Choi et al. (2021) Choi H. G., Park C., Jung S., 2021, Physical Review D, 104
  * Chung & Li (2021) Chung A. K.-W., Li T. G., 2021, Physical Review D, 104
  * Colaço et al. (2023) Colaço L. R., Holanda R. F. L., Nunes R. C., Gonzalez J. E., 2023, Forecasts analysis on varying-$\alpha$ theories from gravitational wave standard sirens (arXiv:2201.04073)
  * Collett & Bacon (2016) Collett T. E., Bacon D. J., 2016, MNRAS, 456, 2210
  * Collett & Bacon (2017) Collett T. E., Bacon D., 2017, Physical Review Letters, 118
  * Collett et al. (2017) Collett T. E., et al., 2017, ApJ, 843, 148
  * Congdon & Keeton (2018) Congdon A. B., Keeton C. R., 2018, Principles of Gravitational Lensing : Light Deflection as a Probe of Astrophysics and Cosmology, 1st ed. 2018. edn. Astronomy and Planetary Sciences, Springer International Publishing, Cham
  * Cremonese et al. (2021) Cremonese P., Ezquiaga J., Salzano V., 2021, Physical Review D, 104
  * Dahle et al. (2013) Dahle H., et al., 2013, ApJ, 773, 146
  * Dai & Venumadhav (2017) Dai L., Venumadhav T., 2017, On the waveforms of gravitationally lensed gravitational waves (arXiv:1702.04724)
  * Dai et al. (2018) Dai L., Li S.-S., Zackay B., Mao S., Lu Y., 2018, Physical Review D, 98
  * Dai et al. (2020) Dai L., Zackay B., Venumadhav T., Roulet J., Zaldarriaga M., 2020, arXiv preprint arXiv:2007.12709
  * De Paolis et al. (2002) De Paolis F., Ingrosso G., Nucita A. A., Qadir A., 2002, A&A, 394, 749
  * Deffayet & Menou (2007) Deffayet C., Menou K., 2007, The Astrophysical Journal, 668, L143–L146
  * Diego (2020) Diego J. M., 2020, Physical Review D, 101
  * Diego et al. (2019) Diego J. M., Hannuksela O. A., Kelly P. L., Pagano G., Broadhurst T., Kim K., Li T. G. F., Smoot G. F., 2019, Astronomy & Astrophysics, 627, A130
  * Diego et al. (2021) Diego J. M., Broadhurst T., Smoot G. F., 2021, Phys. Rev. D, 104, 103529
  * Ezquiaga et al. (2021) Ezquiaga J. M., Holz D. E., Hu W., Lagos M., Wald R. M., 2021, Physical Review D, 103
  * Fairbairn et al. (2023) Fairbairn M., Urrutia J., Vaskonen V., 2023, Journal of Cosmology and Astroparticle Physics, 2023, 007
  * Falco et al. (1985) Falco E. E., Gorenstein M. V., Shapiro I. I., 1985, ApJ, 289, L1
  * Fan et al. (2017) Fan X.-L., Liao K., Biesiada M., Piórkowska-Kurpas A., Zhu Z.-H., 2017, Physical Review Letters, 118
  * Finke et al. (2021) Finke A., Foffa S., Iacovelli F., Maggiore M., Mancarella M., 2021, Physical Review D, 104
  * Gair et al. (2013) Gair J. R., Vallisneri M., Larson S. L., Baker J. G., 2013, Living Reviews in Relativity, 16
  * González et al. (2020) González J. D. R., et al., 2020, The Astrophysical Journal, 902, 44
  * Gorenstein et al. (1988) Gorenstein M. V., Falco E. E., Shapiro I. I., 1988, ApJ, 327, 693
  * Goyal et al. (2021a) Goyal S., Haris K., Mehta A. K., Ajith P., 2021a, Physical Review D, 103
  * Goyal et al. (2021b) Goyal S., D. H., Kapadia S. J., Ajith P., 2021b, Physical Review D, 104
  * Guo & Lu (2022) Guo X., Lu Y., 2022, Physical Review D, 106
  * Gupta et al. (2023) Gupta I., et al., 2023, arXiv preprint arXiv:2307.10421
  * Hannuksela et al. (2019) Hannuksela O. A., Haris K., Ng K. K. Y., Kumar S., Mehta A. K., Keitel D., Li T. G. F., Ajith P., 2019, The Astrophysical Journal, 874, L2
  * Hannuksela et al. (2020) Hannuksela O. A., Collett T. E., Ç alışkan M., Li T. G. F., 2020, Monthly Notices of the Royal Astronomical Society, 498, 3395
  * Haris et al. (2018) Haris K., Mehta A. K., Kumar S., Venumadhav T., Ajith P., 2018, Identifying strongly lensed gravitational wave signals from binary black hole mergers (arXiv:1807.07062)
  * Hernandez (2022) Hernandez I. M., 2022, Measuring the polarization content of gravitational waves with strongly lensed binary black hole mergers (arXiv:2211.01272)
  * Holz & Hughes (2005) Holz D. E., Hughes S. A., 2005, The Astrophysical Journal, 629, 15–22
  * Hou et al. (2021) Hou S., Fan X.-L., Zhu Z.-H., 2021, Monthly Notices of the Royal Astronomical Society, 507, 761–771
  * Iacovelli et al. (2022) Iacovelli F., Finke A., Foffa S., Maggiore M., Mancarella M., 2022, Modified gravitational wave propagation: information from strongly lensed binaries and the BNS mass function (arXiv:2203.09237)
  * Jana et al. (2023) Jana S., Kapadia S. J., Venumadhav T., Ajith P., 2023, Physical Review Letters, 130
  * Janquart et al. (2021a) Janquart J., Hannuksela O. A., Haris K., Van Den Broeck C., 2021a, Monthly Notices of the Royal Astronomical Society, 506, 5430–5438
  * Janquart et al. (2021b) Janquart J., Seo E., Hannuksela O. A., Li T. G. F., Van Den Broeck C., 2021b, The Astrophysical Journal Letters, 923, L1
  * Janquart et al. (2022) Janquart J., Hannuksela O. A., Haris K., Broeck C. V. D., 2022, GOLUM: A fast and precise methodology to search for, and analyze, strongly lensed gravitational-wave events (arXiv:2203.06444)
  * Janquart et al. (2023a) Janquart J., Haris K., Hannuksela O. A., Van Den Broeck C., 2023a, Monthly Notices of the Royal Astronomical Society, 526, 3088–3098
  * Janquart et al. (2023b) Janquart J., et al., 2023b, Monthly Notices of the Royal Astronomical Society, 526, 3832–3860
  * Keeton et al. (1997) Keeton C. R., Kochanek C. S., Seljak U., 1997, The Astrophysical Journal, 482, 604–620
  * Kim et al. (2021) Kim K., Lee J., Yuen R. S. H., Hannuksela O. A., Li T. G. F., 2021, The Astrophysical Journal, 915, 119
  * Kormann et al. (1994) Kormann R., Schneider P., Bartelmann M., 1994, A&A, 284, 285
  * Lagos et al. (2019) Lagos M., Fishbach M., Landry P., Holz D. E., 2019, Phys. Rev. D, 99, 083504
  * Lasky (2015) Lasky P. D., 2015, Publications of the Astronomical Society of Australia, 32, e034
  * Lattimer (2012) Lattimer J. M., 2012, Annual Review of Nuclear and Particle Science, 62, 485–515
  * Li et al. (2018) Li S.-S., Mao S., Zhao Y., Lu Y., 2018, Monthly Notices of the Royal Astronomical Society, 476, 2220–2229
  * Li et al. (2019) Li Y., Fan X., Gou L., 2019, The Astrophysical Journal, 873, 37
  * Li et al. (2023a) Li A. K. Y., Chan J. C. L., Fong H., Chong A. H. Y., Weinstein A. J., Ezquiaga J. M., 2023a, TESLA-X: An effective method to search for sub-threshold lensed gravitational waves with a targeted population model (arXiv:2311.06416)
  * Li et al. (2023b) Li A. K., Lo R. K., Sachdev S., Chan J., Lin E., Li T. G., Weinstein A. J., 2023b, Physical Review D, 107
  * Liao et al. (2017) Liao K., Fan X.-L., Ding X., Biesiada M., Zhu Z.-H., 2017, Nature Communications, 8
  * Liao et al. (2018) Liao K., Ding X., Biesiada M., Fan X.-L., Zhu Z.-H., 2018, The Astrophysical Journal, 867, 69
  * Liu (2023) Liu W., 2023, Gravitational Lens System as a Long Baseline Detector of Extremely Low Frequency Primordial Gravitational Wave (arXiv:2212.02221)
  * Liu et al. (2021) Liu X., Hernandez I. M., Creighton J., 2021, The Astrophysical Journal, 908, 97
  * Liu et al. (2023) Liu A., Wong I. C. F., Leong S. H. W., More A., Hannuksela O. A., Li T. G. F., 2023, Monthly Notices of the Royal Astronomical Society, 525, 4149–4160