# Optimal Dose-Limited Phase Estimation without Entanglement

Stewart A. Koppell<EMAIL_ADDRESS>Mark A. Kasevich  Physics Department,
Stanford University

###### Abstract

Phase estimation is one of the most important facets of quantum metrology,
with applications in sensing, microscopy, and quantum computation. When
estimating a phase shift in a lossy medium, there is an upper bound on the
attainable information per particle sent through the phase shift. Previously,
only entanglement-enhanced measurements have been shown to saturate this
bound. We introduce a measurement scheme which can saturate the bound without
relying on entanglement.

††preprint: APS/123-QED

##  I Introduction

A central theme of quantum metrology is that entanglement in the probe state
can enhance measurement precision. The NOON state interferometer, for example,
achieves Heisenberg-limited precision using a maximally entangled $n$-particle
probe[1]. While Heisenberg-limited measurements are often thought of as a
uniquely quantum phenomenon, the same precision scaling is possible using a
classical probe passed $m$ times through the interferometer. In some
treatments, $n$ and $m$ are defined as equivalent resources and multi-pass
measurements are said to be Heisenberg-limited[2, 3, 4]. To clarify when a
quantum advantage exists, it is important to carefully define the physical
measurement resources and limitations.

When measuring a dynamical quantity such as in magnetometry or gravitational
field detection, measurement bandwidth limits the maximum time the probe can
coherently accrue its phase shift. Once this limit is saturated, the only way
to increase sensitivity is with entanglement in the probe. In these
applications, $m$ and $n$ may be optimized independently and the Heisenberg
limit is exclusive to quantum measurements.

Phase estimation is also used in quantum computing as a subroutine for many
important algorithms (e.g. Shor’s algorithm). Kitaev’s phase estimation
algorithm [5] and the iterative phase estimation algorithm (IPEA) [6, 7] can
achieve a precision scaling like $1/2^{nm}$ using $n$ ancillary qubits and
$m\log{nm}$ iterations [7]. Thus spatial resources (qubits) can be traded for
temporal resources (runtime) while maintaining an exponential advantage.

The discrimination of unitary operations is another noteworthy measurement
task where entanglement and multi-passing can be compared. In general, any
single-pass measurement with a separable probe state repeated a finite number
of times has a non-zero probability of failure. Perfect discrimination was
first described as an application of entanglement [8, 9]. Later, in a letter
which posed the question “What kind of tasks can be achieved without
entanglement?”, it was shown that a multi-pass (MP) measurement can also
suffice [10].

The primary interests of this Letter are applications like microscopy where
the measurement is performed on a material sample which may be destroyed or
altered by inelastic interactions with the probe. Assuming the sample can be
modeled as a static phase shift, the goal is to maximize the information per
dose $d$ (probe intensity incident on the sample). For single-pass
interferometer measurements, the dose metric simply counts the number of
particles which enter the sample channel. A similar resource accounting is
performed when the measurement has access to a strong local oscillator (i.e.
none of the probe intensity needs to be allocated to a reference channel). For
MP measurements, the dose metric takes into account the losses on each pass,
which will reduce the dose inflicted on subsequent passes (in general,
$d<n\times m$).

In optical microscopy, both entanglement-enhanced measurements[11, 12, 13] and
MP measurements[14] have been shown to increase imaging sensitivity per dose.
Similar schemes for improving the dose-limited resolution of transmission
electron microscopes have been proposed[15, 16, 17]. However until now it has
been unclear whether entanglement is necessary to achieve the best possible
measurement sensitivity per dose.

We will employ the Quantum Fisher Information (QFI) formalism to address this
question. Via the Quantum Cramer Rao Bound (QCRB), the QFI provides a lower
limit on the error $\Delta_{\theta}$ in an estimate of an unknown parameter
$\theta$[18, 19, 20]. If a measurement attaining QFI $\mathcal{J}$ is repeated
$N$ times, the QCRB is

$\Delta_{\theta}^{2}\geq\frac{1}{N\mathcal{J}}\ .$ (1)

The QFI for a lossless NOON-state interferometer is proportional to $n^{2}$.
This quadratic scaling, characteristic of the Heisenberg limit, is not the
only hallmark of quantum advantage. In fact, in the presence of any amount of
loss, the QFI cannot asymptotically scale better than $n$ (the standard
quantum limit)[21, 22]. Instead, quantum measurements may confer a constant
advantage over the classical alternatives.

Using dose as a measurement resource, the figure of merit to optimize is
$\xi=\mathcal{J}/d$, which we will refer to as the dose efficiency. For a
sample with transmissivity $\eta$, the quantum limit is [21]

$\xi_{\text{QL}}=\frac{4\eta}{(1-\eta)}\ .$ (2)

Here we will describe a measurement which can saturate this limit without an
entanglement-enhanced probe.

##  II Dose as a Measurement Resource

It will be instructive to begin by calculating $\xi$ for several well-known
measurements. A standard approach to calculating the QFI in the presence of
loss involves Kraus operators, which describe the dynamics of open quantum
systems[23, 24, 25]. We can take simpler approach when the probe is a pure
single-particle state or a NOON state.

As a result of possible absorption, the exit state will generally be mixed:
$\rho=\sum p_{i}\left|\psi_{i}\right\rangle\left\langle\psi_{i}\right|$. If
$\mathcal{J}[\psi_{i}]$ is the QFI for $\psi_{i}$, then
$\mathcal{J}[\rho]\leq\sum p_{i}\mathcal{J}[\psi_{i}]$ (i.e. the QFI is
convex). The equality holds only when each $\left|\psi_{i}\right\rangle$ is
measured optimally in the same basis. When the absorption of a single particle
projects the probe into a state which yields no information, then
$\mathcal{J}=p\tilde{\mathcal{J}}$ where $\tilde{\mathcal{J}}$ is the QFI
conditional on no absorption, which happens with probability $p$.

We begin with a lossy Mach-Zehnder interferometer (MZI) using a single-
particle probe. If $\alpha$ is the amplitude sent to the sample arm, the
normalized exit wavefunction is

$\psi(\theta)=\frac{1}{\sqrt{p}}\begin{bmatrix}\sqrt{1-\alpha^{2}}\\\
\alpha\sqrt{\eta}e^{i\theta}\end{bmatrix}\ .$ (3)

with $p=1-(1-\eta)\alpha^{2}$. The QFI for a pure state is [18]

$\mathcal{J}=4\langle\dot{\psi}|\dot{\psi}\rangle-4|\langle\dot{\psi}|\psi\rangle|^{2}\
,$ (4)

where the dot indicates a derivative with respect to $\theta$. So the the
lossy MZI has

$\tilde{\mathcal{J}}=\frac{4}{p}|\alpha|^{2}\eta\left(1-\frac{1}{p}|\alpha|^{2}\eta\right)\
.$ (5)

The expected QFI per dose $d=|\alpha^{2}|$ is maximized in the limit
$\alpha\rightarrow 0$, yielding

$\xi=\frac{p\tilde{\mathcal{J}}}{d}=4\eta\ .$ (6)

The principle that the dose efficiency is maximized when the intensity in the
sample channel is infinitesimal will be true for each of the measurements we
discuss below (when $\alpha$ is strictly zero, the measurements provide no
information). Compared to the quantum limit, the MZI yields less information
per dose by a factor of $1/(1-\eta)$, which diverges as $\eta\rightarrow 1$.

In the presence of loss, NOON state interferometers are sub-optimal because
they fail when a single particle is absorbed. Nevertheless, an unbalanced NOON
state (with probability $|\alpha|^{2}\ll 1$ of finding $n$ particles in the
sample channel and probability $1-|\alpha|^{2}$ of finding $n$ particles in
the reference channel) is a significant improvement over a single particle
probe for a MZI. For $\eta\sim 1$, $\xi_{\text{NOON}}/\xi_{\text{QL}}\gtrsim
1/e$ for optimal $n$. Explicitly, $\mathcal{J}=4n^{2}\eta^{n}|\alpha|^{2}$ and
$d=n|\alpha|^{2}$, so

$\xi_{\text{NOON}}=4n\eta^{n}\ .$ (7)

which is maximized by choosing $n\sim-1/\ln{\eta}$.

MP and NOON state interferometers provide the same $\mathcal{J}$ when $n=m$
and are both spoiled by a single absorption. However, the dose from a MP
measurement is limited: there can be no more than a one absorption. As a
result, MP measurements can be more dose efficient. For optimal $m$,
$\xi_{\text{MP}}/\xi_{\text{QL}}\gtrsim 0.65$, and in general
$J=4m^{2}\eta^{m}$ and $d=\sum_{k=0}^{m-1}\eta^{k}$ so

$\xi_{\text{MP}}=4m^{2}\eta^{m}\frac{1-\eta}{(1-\eta^{m})}\ .$ (8)

While a MP interferometer does not fully achieve the quantum limit, the
simplicity of the scheme may make it the practical choice in many
applications.

The expected number of particles in a pure single mode Gaussian state is
$\left\langle n\right\rangle=\alpha^{2}+\sinh^{2}(r)$ where $\alpha$ is the
displacement parameter and $r$ is the squeezing parameter (measured in dB, the
squeezing is $20\ln{r}/10$). The number of particles contributing to the
squeezing is $n_{\text{sq}}=\sinh^{2}(r)$. If $\left\langle n\right\rangle\gg
n_{\text{sq}}$, then [26]

$\xi_{\text{SQZ}}(n_{\text{sq}})=\frac{4\eta}{2\eta
n_{\text{sq}}-2\eta\sqrt{(n_{\text{sq}}(n_{\text{sq}}+1))}+1}\ .$ (9)

When $n_{\text{sq}}$ is large, the denominator approaches $1-\eta$ and this
measurement saturates the quantum limit for dose efficiency.

Fig. 1 Shows $\xi/\xi_{\text{QL}}$ for each of these measurements vs (A) the
absorption probability $1-\eta$ and (B) the relevant measurement parameter
($m$, $n$, or $n_{\text{sq}}$) with $\eta=0.9$.

Figure 1:  QFI per dose ($\xi$) relative to the quantum limit
($\xi_{\text{QL}}$) for measurements with: a single particle probe and single-
pass (SP), a maximally entangled probe of size $n$ (NOON), a single-particle
passed $m$ times through the sample (MP), an $m-$stage chain interferometer
with constant beamsplitter angles (CIC), an $m$-stage chain interferometer
with optimized beamsplitter angles (CIO), and a squeezed Gaussian probe state
(SQZ) with $n_{\text{sq}}$ particles contributing to the squeezing. (A)
$\xi/\xi_{\text{QL}}$ vs probability of absorption by the sample (for each
particle on each pass). For the NOON, MP, and CPC measurements, $m$ and $n$
are optimal. The CIO measurement performance is calculated for $m$=4, 32, and
128 stages. The striped background shows the amount of squeezing (in dB)
needed for a Gaussian probe to have equivalent performance. (B)
$\xi/\xi_{\text{QL}}$ vs the relevant measurement parameter for $\eta=0.9$
(marked by a vertical black line in (A)).

##  III Phase Estimation with a Chain Interferometer

Figure 2:  Chain interferometer scheme for estimation of an unknown phase
$\theta$. In each of $m$ stages, the reference channel (black) and sample
channel (red) are weakly coupled with a beamsplitter (BS) before the probe
interacts with the unknown phase shift in the sample arm and a programmed
phase shift in the reference arm. The result is output wavefunction
$\psi(\theta)$. The QFI of $\psi(\theta)$ is related to the variance in the
estimate of $\theta$ using an optimal measurement $M$.

The most general possible $m$-stage entanglement-free phase estimation scheme
could use up to $m$ auxiliary channels in addition to the sample channel, with
unitary transformations applied to couple all the channels between each stage.
This would likely be infeasible to implement as a microscope architecture. We
will show that an optimal measurement is possible using only two channels in a
configuration we call a chain interferometer (CI). Figure 2 shows a CI with
$m$ stages. At each stage, a beamsplitter couples the reference channel to a
sample channel. In practice, the measurement may be applied to a single,
unique sample by using fast-switching optics to fold the optical path onto
itself, in which case only one beamsplitter is required. This architecture has
been proposed before in the context of interaction-free measurement (IFM),
where the presence of an absorber in the sample channel can be detected with
an arbitrarily small chance of damage (loss)[27].

The CI can be optimized by adjusting $m$, the beamsplitter transmissivity $T$,
and reference phase shift $\gamma$. When $\gamma=0$ and $T=\sin(\pi/2m)^{2}$,
the CI implements IFM. However this configuration does not saturate the
quantum limit for dose-limited phase estimation. This architecture has been
numerically optimized in Ref. [26], which found that when
$T\approx(0.01/2m)^{2}$ and $\gamma=0$ (assuming small $\theta$), the CI could
attain $\xi/\xi_{\text{QL}}>0.9$.

We will take a more analytical approach and allow the beamsplitter
transmissivity to change for each stage. We could also allow the reference
phase shift to change for each stage, but this turns out to be unnecessary.
Given the discussion in Section II, we will allocate the majority of the probe
amplitude to the reference arm with the ansatz $T_{k}<\epsilon^{2}$ for some
$\epsilon\ll 1$. We will also set $\gamma\sim\theta$ (this choice requires
some prior knowledge of the $\theta$, which is always an implicit assumption
of the QFI formalism[18]).

Unlike a simple MP interferometer, the exit wavefunction of a CI contains
terms corresponding to different numbers of passes through the sample.
Interference between these terms provides extra information. The details of
calculating $\mathcal{J}$ and $d$ are given in the supplementary material. To
order $\epsilon^{2}$,

$\mathcal{J}=4\left|\sum_{k=0}^{m-1}(m-k)\eta^{(m-k)/2}\tau_{k}\right|^{2}$
(10)

where $\tau_{k}=\sqrt{T_{k}}$ and

$d=\sum_{k=0}^{m-1}\left|\sum_{k^{\prime}=0}^{k}\eta^{(k-k^{\prime})/2}\tau_{k^{\prime}}\tau_{k}\right|^{2}\
.$ (11)

We found the set of $\tau_{k}$ which maximize $\xi$ by analytically solving
the $m=2$ and $m=3$ case and extrapolating. The result is

$\tau_{k}=\begin{cases}\epsilon\eta^{m/2}&k=0\\\
\epsilon(1-\eta)\eta^{(m-k)/2}&k>0\end{cases}\ .$

Using this prescription, we find

$\lim_{m\rightarrow\infty}\mathcal{J}=\epsilon^{2}\frac{4\eta^{2}}{(1-\eta)^{2}}\quad\text{and}\quad\lim_{m\rightarrow\infty}d=\epsilon^{2}\frac{\eta}{1-\eta}\
.$ (12)

so that $\xi\rightarrow\xi_{QL}$ for large $m$.

In Fig. 1 (A) and (B), the lines labeled CIC are $\xi/\xi_{\text{QL}}$ for a
CI with a constant, infinitesimal beamsplitter transmissivity using the
optimal number of stages. The lines labeled CIO show $\xi/\xi_{\text{QL}}$ for
a CI with optimized beamsplitter transmissivities. For $\eta=0.9$, a
10(32)-stage CI (using either strategy) is more dose efficient than a standard
interferometer using a Gaussian probe squeezed by 10(20) dB.

##  IV Optical Losses

So far we have neglected losses other than those caused by the sample. We will
now consider the effects of loss during probe preparation (probability
$1-\eta_{\text{P}}$), loss between each stage of a MP or chain interferometer
(probability $1-\eta_{\text{rt}}$), and loss during detection (probability
$1-\eta_{\text{D}}$).

Loss during probe preparation decreases intensity without affecting the dose
efficiency for single-particle measurements. However this type of loss can
still reduce the effectiveness of entangled probes. For squeezed states, this
loss injects vacuum noise, putting a ceiling on the maximum effective
squeezing. For the same reason, entanglement-enhanced measurements also suffer
disproportionately from loss at the detection stage.

The dose efficiencies for each of the measurements discussed so far are
derived in the supplemental document. For CI and squeezed-state measurements,
the dose efficiency has upper bound

$\xi\leq\frac{4\eta\eta_{\text{D}}}{(1-\eta\eta^{*})}\ ,$ (13)

which is saturated in the limit $m\rightarrow\infty$ or
$n_{\text{sq}}\rightarrow\infty$, respectively. For a CI,
$\eta^{*}=\eta_{\text{rt}}$, and for a squeezed state interferometer,
$\eta^{*}=\eta_{P}\eta_{\text{D}}$. When
$\eta_{\text{rt}}<\eta_{\text{P}}\eta_{\text{D}}$, the CI has superior
performance. The optimal beamsplitter transmissivities $\tau_{k}$ are adjusted
in the presence of loss as described in the supplementary material.

It may be possible to build a MP interferometer with lower loss than is
possible with a CI. If, for example, $\eta_{\text{rt}}=0.99$ for MP, then it
becomes more effective then CIO with $\eta_{\text{rt}}=0.95$ for $\eta>0.94$.
Low MP loss also makes it feasible to combine multi-passing and squeezing. A
multi-pass squeezed state interferometer (with $\eta_{\text{rt}}=0.99$)
surpasses CIO (with $\eta_{\text{rt}}=0.95$) for $\eta>0.84$.

In Fig 3 we show the dose efficiency of measurements with
$\eta_{\text{rt}}=0.95$ and $\eta_{\text{P}}=\eta_{\text{D}}=0.9$ relative to
the quantum limit for $\eta_{\text{rt}}=\eta_{\text{P}}=\eta_{\text{D}}=1$.

Figure 3:  QFI per dose with optical losses outside of the sample. We assume
10% loss during both probe preparation and detection and 5% loss at each stage
of the MP and CI measurements. The measurements are as described in Fig. 1.
The curve labeled SQZ represents the upper bound on squeezing enhancement
($n_{\text{sq}}\rightarrow\infty$), and the dotted line labeled MPSQZ is the
bound measurements with both squeezing and multiple passes.

##  V Conclusion

We have identified a new type of phase estimation protocol which attains the
quantum limit for information per dose without using entanglement. This makes
it possible, in principle, to employ a quantum-optimal measurement in
microscopy using only basic optical elements. When taking into account losses
in the optics for forming the probe and performing the measurement of the exit
wavefunction, we found single-particle multi-stage measurements can
potentially outperform the best known entanglement-enhanced measurements.

##  Appendix A Dose efficiency for a CI

Let $\psi^{(k)}_{j}$ be the wavefunction just before encountering the sample
on pass $k$ (numbered $k=0$ to $k=m-1$), with index $j=0$ and $j=1$ referring
to the reference and sample channels, respectively. We will assume the probe
is a single particle in initial state

$\psi^{0}=\begin{bmatrix}1\\\ 0\end{bmatrix}\ .$ (14)

The output wavefunction will be written $\psi_{j}$. We will assume that the
coupling between the reference and sample channels is infinitesimal. As a
result, the QFI can be written

$\mathcal{J}=4p\left\langle\left.\dot{\psi}_{j}\right|\dot{\psi}_{j}\right\rangle$
(15)

where $1-p$ is the probability that the sample absorbs the probe. The dose is
defined as

$d=\sum_{k=0}^{m-1}\left|\psi^{(k)}_{1}\right|^{2}\ .$ (16)

Let $B_{k}$ and $S$ be operators representing actions of the beamsplitters and
sample, respectively:

$B_{k}=\begin{bmatrix}\sqrt{1-\tau_{k}^{2}}&\tau_{k}\\\
\tau_{k}&\sqrt{1-\tau_{k}^{2}}\end{bmatrix}$ (17)

and

$S_{k}=\begin{bmatrix}1&0\\\ 0&\sqrt{\eta}e^{i\theta}\end{bmatrix}\ .$ (18)

With the assumption that the beamsplitter transmissivity is small
($\tau_{k}<\epsilon\ll 1$), we can ignore terms of order $\epsilon^{3}$ and
higher. To this order, the output wavefunction is

$\begin{bmatrix}\psi_{0}\\\
\psi_{1}\end{bmatrix}=\frac{1}{\sqrt{p}}\begin{bmatrix}1-\frac{1}{2}\underset{k}{\sum}\tau_{k}^{2}-\underset{k>k^{\prime}}{\sum}\tau_{k}\tau_{k^{\prime}}\eta^{(k^{\prime}-k)/2}e^{i(k-k^{\prime})\theta}\\\
-\sum_{k}\tau_{k}\eta^{(m-k)/2}e^{i(m-k)\theta}\end{bmatrix}\ .$ (19)

Using equations 15 and 16, we arrive at the expressions for $\mathcal{J}$ and
$d$ found in the main text. Notice that

$|\langle\dot{\psi}_{j}|\psi_{j}\rangle|^{2}\sim\mathcal{O}(\epsilon^{4})\ ,$

justifying the approximation in 15.

##  Appendix B Formulas for $\xi$ with optical loss

If $1-\eta_{P}$ and $1-\eta_{D}$ are the losses during probe preparation and
detection, respectively, then an MZI with intensity $|\alpha|^{2}\ll 1$ in the
sample arm has

$\mathcal{J}=4\eta_{\text{P}}\eta\eta_{\text{D}}|\alpha|^{2},\quad\text{and}\quad
d=\eta_{P}|\alpha|^{2}$

so

$\xi=4\eta\eta_{\text{D}}\ .$ (20)

If the loss between passes in a MP interferometer is $1-\eta_{\text{rt}}$,
then

$\displaystyle J$
$\displaystyle=4m^{2}\eta_{\text{P}}\eta^{m}\eta_{\text{rt}}^{m-1}\eta_{\text{D}}\quad\text{and}$
$\displaystyle d$
$\displaystyle=\eta_{\text{P}}\sum_{k=0}^{m-1}(\eta_{\text{rt}}\eta)^{k}\ ,$

so

$\xi_{\text{MP}}=4m^{2}\eta^{m}\eta_{\text{rt}}^{m-1}\eta_{\text{D}}\frac{1-\eta_{\text{rt}}\eta}{1-(\eta_{\text{rt}}\eta)^{m}}$
(21)

Notice $\xi$ is independent of $\eta_{P}$ for MP and SP measurements: loss
during probe preparation prolongs the measurement but does not affect dose
efficiency. This will also be the case for a CI, but not for measurements with
an entangled probe.

For single-pass entanglement-enhanced measurements, the location of the
absorption will not affect $\mathcal{J}$ (though it may effect $d$). For a
NOON state interferometer with $n$-particle amplitude $\alpha$ in the sample
arm,

$\mathcal{J}=4n^{2}(\eta_{\text{P}}\eta\eta_{\text{D}})^{n}\quad\text{and}\quad
d=\eta_{P}n|\alpha|^{2}$

so

$\xi_{\text{NOON}}=4n\eta_{\text{P}}^{n-1}(\eta\eta_{\text{D}})^{n}\ .$ (22)

For a squeezed-state interferometer with total loss $\eta_{\text{tot}}$ (with
large $n_{\text{sq}}$),

$\mathcal{J}=\frac{4|\alpha|^{2}\eta_{\text{tot}}}{1-\eta_{\text{tot}}}$

Using $\eta_{\text{tot}}=\eta_{\text{P}}\eta\eta_{\text{D}}$ and
$d=|\alpha|^{2}\eta_{\text{P}}$,

$\xi_{\text{SQZ}}=\frac{4\eta\eta_{\text{D}}}{1-\eta_{\text{P}}\eta\eta_{\text{D}}}$
(23)

When there is no loss outside of the sample, then single-pass squeezing is an
optimal strategy. However, depending on the size of $\eta_{\text{rt}}$
relative to $\eta_{\text{P}}$ and $\eta_{\text{D}}$, a MP squeezed-state
measurement can outperform squeezing or MP alone. The total loss is
$\eta_{\text{tot}}=\eta_{\text{P}}\eta^{m}\eta_{\text{rt}}^{m-1}\eta_{\text{D}}$
and

$\displaystyle\mathcal{J}$
$\displaystyle=\frac{4|\alpha|^{2}m^{2}\eta_{\text{tot}}}{1-\eta_{\text{tot}}}\
,$ $\displaystyle d$
$\displaystyle=|\alpha|^{2}\eta_{P}\frac{1-(\eta_{\text{rt}}\eta)^{m}}{1-\eta_{\text{rt}}\eta}\
.$

Finally, a lossy CI with has

$\displaystyle\mathcal{J}$
$\displaystyle=4\eta_{\text{P}}\eta_{\text{rt}}^{m-1}\eta_{\text{D}}\left|\sum_{k=0}^{m-1}(m-k)\eta^{(m-k)/2}\tau_{k}\right|^{2},$
$\displaystyle d$
$\displaystyle=\eta_{\text{P}}\sum_{k=0}^{m-1}\eta_{\text{rt}}^{k}\left|\sum_{k^{\prime}=0}^{k}\eta^{(k-k^{\prime})/2}\tau_{k^{\prime}}\right|^{2}$

and the set of $\tau_{k}$ which maximize $\xi$ are:

$\tau_{0}=\epsilon(\eta_{\text{rt}}\sqrt{\eta})^{m}\quad\tau_{k>0}=\epsilon(1-\eta_{\text{rt}}\eta)(\eta_{\text{rt}}\sqrt{\eta})^{m-k}$

giving

$\xi_{\text{CI}}=\frac{4\eta\eta_{\text{D}}}{1-\eta\eta_{\text{rt}}}$ (24)

in the limit of large $m$.

## References

  * Holland and Burnett [1993] M. J. Holland and K. Burnett, Interferometric detection of optical phase shifts at the heisenberg limit, Phys. Rev. Lett. 71, 1355 (1993).
  * Higgins  _et al._ [2007] B. Higgins, D. W. Berry, S. D. Bartlett, H. M. Wiseman, and G. J. Pryde, Entanglement-free heisenberg-limited phase estimation, Nature 450, 393 (2007).
  * de Burgh and Bartlett [2005] M. de Burgh and S. D. Bartlett, Quantum methods for clock synchronization: Beating the standard quantum limit without entanglement, Phys. Rev. A 72, 042301 (2005).
  * Luis [2002] A. Luis, Phase-shift amplification for precision measurements without nonclassical states, Phys. Rev. A 65, 025802 (2002).
  * Kitaev [1996] A. Y. Kitaev, Quantum measurements and the abelian stabilizer problem, Electron. Colloquium Comput. Complex. 3 (1996).
  * Griffiths and Niu [1996] R. B. Griffiths and C.-S. Niu, Semiclassical fourier transform for quantum computation, Phys. Rev. Lett. 76, 3228 (1996).
  * Dobšíček  _et al._ [2007] M. Dobšíček, G. Johansson, V. Shumeiko, and G. Wendin, Arbitrary accuracy iterative quantum phase estimation algorithm using a single ancillary qubit: A two-qubit benchmark, Phys. Rev. A 76, 030306 (2007).
  * Acín [2001] A. Acín, Statistical distinguishability between unitary operations, Phys. Rev. Lett. 87, 177901 (2001).
  * D’Ariano  _et al._ [2001] G. M. D’Ariano, P. Lo Presti, and M. G. A. Paris, Using entanglement improves the precision of quantum measurements, Phys. Rev. Lett. 87, 270404 (2001).
  * Duan  _et al._ [2007] R. Duan, Y. Feng, and M. Ying, Entanglement is not necessary for perfect discrimination between unitary operations., Physical review letters 98 10, 100503 (2007).
  * Humphreys  _et al._ [2013] P. C. Humphreys, M. Barbieri, A. Datta, and I. A. Walmsley, Quantum enhanced multiple phase estimation, Physical review letters 111, 070403 (2013).
  * Israel  _et al._ [2014] Y. Israel, S. Rosen, and Y. Silberberg, Supersensitive polarization microscopy using noon states of light, Physical review letters 112, 103604 (2014).
  * Casacio  _et al._ [2021] C. A. Casacio, L. Madsen, A. Terrasson, M. Waleed, K. Barnscheidt, B. Hage, M. Taylor, and W. Bowen, Quantum-enhanced nonlinear microscopy., Nature 594 7862, 201 (2021).
  * Juffmann  _et al._ [2016] T. Juffmann, B. B. Klopfer, T. L. Frankort, P. Haslinger, and M. A. Kasevich, Multi-pass microscopy, Nature Communications , 1 (2016).
  * Okamoto and Nagatani [2013] H. Okamoto and Y. Nagatani, Entanglement-assisted electron microscopy based on a flux qubit, Applied Physics Letters 104 (2013).
  * Putnam and Yanik [2009] W. P. Putnam and M. F. Yanik, Noninvasive electron microscopy with interaction-free quantum measurements, Physical Review A 80, 040902 (2009).
  * Juffmann  _et al._ [2017] T. Juffmann, S. A. Koppell, B. B. Klopfer, C. Ophus, R. M. Glaeser, and M. A. Kasevich, Multi-pass transmission electron microscopy, Scientific Reports 7, 1 (2017).
  * Holevo [1982] A. S. Holevo,  _Probabilistic and Statistical Aspects of Quantum Theory_ (Amsterdam: North- Holland, 1982).
  * Helstrom [1967] C. W. Helstrom, Minimum mean-squared error of estimates in quantum statistics, Physics Letters A 25, 101 (1967).
  * Braunstein and Caves [1994] S. L. Braunstein and C. M. Caves, Statistical distance and the geometry of quantum states, Phys. Rev. Lett. 72, 3439 (1994).
  * Escher  _et al._ [2012] B. Escher, R. Filho, and L. Davidovich, General framework for estimating the ultimate precision limit in noisy quantum-enhanced metrology, Nature Physics 7 (2012).
  * Demkowicz-Dobrzański  _et al._ [2012] R. Demkowicz-Dobrzański, J. J. Kolodynski, and M. Guta, The elusive heisenberg limit in quantum-enhanced metrology, Nature Communications 3, 1063 (2012).
  * Dorner  _et al._ [2009] U. Dorner, R. Demkowicz-Dobrzanski, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Optimal quantum phase estimation, Phys. Rev. Lett. 102, 040403 (2009).
  * Demkowicz-Dobrzanski  _et al._ [2009] R. Demkowicz-Dobrzanski, U. Dorner, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Quantum phase estimation with lossy interferometers, Phys. Rev. A 80, 013825 (2009).
  * Demkowicz-Dobrzański  _et al._ [2015] R. Demkowicz-Dobrzański, M. Jarzyna, and J. Kołodyński, Chapter four - quantum limits in optical interferometry (Elsevier, 2015) pp. 345–435.
  * Birchall  _et al._ [2017] P. M. Birchall, J. L. O’Brien, J. C. F. Matthews, and H. Cable, Quantum-classical boundary for precision optical phase estimation, Phys. Rev. A 96, 062109 (2017).
  * Kwiat  _et al._ [1995] P. Kwiat, H. Weinfurter, T. Herzog, A. Zeilinger, and M. A. Kasevich, Interaction-free measurement, Physical Review Letters 74, 4763 (1995).


# Inferring properties of dark galactic halos using strongly lensed
gravitational waves

Eungwang Seo  SUPA, School of Physics and Astronomy, University of Glasgow,
Glasgow G12 8QQ, United Kingdom  Department of Physics, The Chinese University
of Hong Kong, Shatin, NT, Hong Kong  Tjonnie G. F. Li  Department of Physics,
The Chinese University of Hong Kong, Shatin, NT, Hong Kong  Institute for
Theoretical Physics, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven, Belgium
Department of Electrical Engineering (ESAT), KU Leuven, Kasteelpark Arenberg
10, B-3001 Leuven, Belgium  Martin A. Hendry  SUPA, School of Physics and
Astronomy, University of Glasgow, Glasgow G12 8QQ, United Kingdom

###### Abstract

Gravitational waves (GWs) can be distorted, similarly to electromagnetic (EM)
waves, by massive objects through the phenomenon of gravitational lensing. The
importance of gravitational lensing for GW astronomy is becoming increasingly
apparent in the GW detection era, in which nearly 100 events have already been
detected. As current ground-based interferometers reach their design
sensitivities, it is anticipated that these detectors may observe a few GW
signals that are strongly lensed by the dark halos of intervening galaxies or
galaxy clusters. Analyzing strong lensing effects on GW signals is, thus,
becoming important to understand the lens’ properties and correctly infer the
intrinsic GW source parameters. However, one cannot accurately infer lens
parameters for complex lens models with only GW observations because there are
strong degeneracies between the parameters of lensed waveforms. In this paper,
we discuss how to conduct parameter estimation of strongly lensed GW signals
and infer the lens parameters with the help of EM observations, such as the
shape and the redshift of the lens. We find that for simple spherically
symmetric lens models, the lens parameters can be well recovered using only GW
information. On the other hand, recovering the lens parameters requires
systems in which four or more GW images are detected with additional EM
observations for non-axially symmetric lens models. Combinations of GW and EM
observations can further improve the inference of the lens parameters.

patch/longtable=false

††thanks<EMAIL_ADDRESS>
##  1 Introduction

Measuring the masses of the dark halos of galaxies or galaxy clusters is one
of the most important challenges in astrophysics and cosmology. The halo
masses are a key factor for understanding important phenomena relevant to the
physics of galaxies, from galaxy formation and evolution (Behroozi et al.,
2010; Fanidakis et al., 2012; Wetzel et al., 2012) to galaxy quenching (Woo et
al., 2013; Zu & Mandelbaum, 2016; Wang et al., 2018; Behroozi et al., 2019).
Most methods developed to date for measuring masses of galaxies depend only on
electromagnetic (EM) wave observations. Application of photometry to
gravitational lensing events and spectroscopy for velocity dispersion
measurements of galaxies and clusters have been widely used – and have
estimated the mass of galaxies with mild error ranges (Broadhurst et al.,
1994; Hoekstra, 2003; Mandelbaum et al., 2006; Hansen et al., 2009; Elahi et
al., 2018; Kafle et al., 2018). However, EM waves can be easily absorbed and
re-emitted by baryonic matter. Also, accurate measurements of galaxy
components’ dynamics are necessary, which are challenging.

In this regard, one promising approach to measuring dark halo masses is to use
gravitational waves (GWs), which are not subject to some of the same
limitations as EM waves. GWs can also be gravitationally lensed by massive
objects like EM waves – resulting in larger angular areas on the lens plane
(Wang et al., 1996; Nakamura, 1998; Takahashi & Nakamura, 2003; Sathyaprakash
& Schutz, 2009). Unlike EM waves, the universe is nearly transparent to GWs,
and lens mass measurements using GWs are, in principle, more straightforward
by analyzing their strain amplitudes.

Among the different kinds of GW lensing that may be considered, strong lensing
occurs when there are massive objects along the line of sight, such as
galaxies and galaxy clusters, which have a larger Schwarzschild radius than
the wavelengths of the propagating GW signals. This phenomenon splits the
original unlensed GW signal into multiple lensed signals, each of which has a
different amplitude, arrival time, and phase. Until now, confident evidence of
strong lensing has not been found in the O1, O2, and O3 events (Hannuksela et
al., 2019; Abbott et al., 2021, 2023). In the near future, however, it is
expected that up to $\mathcal{O}(1\sim 2)$ strongly lensed compact binary
coalescences per year may be observed as ground-based GW detectors reach their
design sensitivities (Biesiada et al., 2014; Ding et al., 2015; Ng et al.,
2018; Li et al., 2018; Mukherjee et al., 2021; Wierda et al., 2021).

To investigate the properties of the lens from analysis of lensed GW signals,
a lens model has to be specified, which characterizes the properties of dark
galactic halos. In previous works, axially symmetric lens models were mainly
considered for simplicity (Cao et al., 2014; Hannuksela et al., 2019; Hou et
al., 2020; Robertson et al., 2020; Abbott et al., 2021; Cheung et al., 2021;
Seo et al., 2022). However, non-axially symmetric lens models should be
considered for more realistic cases.

In this paper, we perform Bayesian parameter estimation and employ rejection
sampling techniques to infer the properties of the dark halos of intervening
galaxies using various lens models: 1) simple lens models that have previously
been used, such as a point mass and singular isothermal sphere (SIS) and 2)
more complicated lens models that can be non-isothermal and have ellipticity.
Furthermore, we also calculate the lensing cross-section of each density
profile to investigate which lens model is more effective to generate GW
strong lensing.

Throughout this paper, we assume the so-called “Planck 15” cosmology,
$H_{0}=67.8$ km $s^{-1}Mpc^{-1}$, $\Omega_{m}=0.308$, $\Omega_{K}=0$ and
$\Omega_{\Lambda}=0.692$ (Ade et al., 2016).

##  2 Strong lensing configuration

The thin lens approximation holds in astrophysical gravitational lensing
systems since distances between the observer and the source and lens are very
large. Thus, one can describe a gravitational lensing system using a two-
dimensional plane. Figure. 1 shows the geometry of the system in detail, where
$D_{L}$, $D_{S}$, and $D_{LS}$ are the angular diameter distances between
observer and lens, observer and source, and source and lens, respectively. It
is convenient to express the image and source positions in dimensionless
values.

$x=\frac{|\vec{\xi}|}{\xi_{0}},\>\>y=\frac{D_{L}|\vec{\eta}|}{D_{S}\xi_{0}},$
(1)

where $\vec{\xi}$ and $\vec{\eta}$ are the displacements from the lens to an
image and from the line of sight to the source, respectively and $\xi_{0}$ is
a normalization constant which makes $x$ and $y$ dimensionless. Since we
assume a dark galactic halo as our lens scale, which has a larger
Schwarzschild radius than the wavelength of the propagating GWs from typical
binary black hole systems detected by current ground-based interferometers, we
can apply the geometrical optics approximation to calculate the lensed GW
signal (Wang et al., 1996; Takahashi & Nakamura, 2003; Dai & Venumadhav,
2017). Within the geometrical optics regime, the lensed GW signal $h_{l}$ and
unlensed GW signal $h_{u}$, both defined in the frequency domain $f$, have the
following relationship,

$h_{l}(f,\boldsymbol{\theta_{s}},\boldsymbol{\theta_{l}},M^{z}_{l})=F_{\rm
geo}(f,\boldsymbol{\theta_{l}},M^{z}_{l})\times
h_{u}(f,\boldsymbol{\theta_{s}}),$ (2)

where $\boldsymbol{\theta_{s}}$ and $\boldsymbol{\theta_{l}}$ are parameters
of the GW source and model-dependent lens parameters respectively, and
$M^{z}_{l}=M_{l}(1+z_{l})$ is the redshifted lens mass of the intervening
halo. $F_{\rm geo}$ is the amplification factor, which is given by (see
Takahashi & Nakamura (2003))

$F_{\rm
geo}(f,\boldsymbol{\theta_{l}},M^{z}_{l})=\sum_{j}|\mu_{j}|^{0.5}\exp{i\pi\left(\frac{8fGM^{z}_{l}}{c^{3}}T_{j}-n_{j}\right)},$
(3)

where $\mu_{j}$ is the magnification factor of the $j$-th image at source
position $x_{j}$, and $T_{j}$ and $n_{j}$ are the Fermat potential and Morse
index of the $j$-th image, respectively.

Figure 1: Schematic of the gravitational lensing system. GW signals from a
binary black hole on the source plane go through different paths to produce
multiple images by a dark galactic halo on the lens plane. $\xi$ is the
distance between lens and image positions, the so-called impact parameter, and
$\eta$ is the source position from the line of sight. $D_{L}$ and $D_{S}$ are
the angular diameter distances between the observer and source and the
observer and lens plane, respectively, and $D_{LS}$ is the angular diameter
distance between the source and lens plane.

Magnification factors and Fermat potentials can have different values for the
same $\boldsymbol{\theta_{l}}$ depending on what lens model is assumed. A lens
model describes how much an image is distorted at an image position, with its
deflection potential. Thus, it is useful to show how deflection potentials
determine the two parameters above. The magnification factor and the Fermat
potential for the $j$-th image can be written as

$\displaystyle\mu_{j}=\frac{1}{\operatorname{det}\left(\delta_{ab}-\frac{\partial^{2}\psi_{j}}{\partial
x_{a}\partial x_{b}}\right)},$ $\displaystyle
T_{j}=\frac{1}{2}\left(|\vec{x}_{j}-\vec{y}|\right)^{2}-\psi_{j}+\phi_{m},$
(4)

where $\psi_{j}$ is the deflection potential at the $j$-th image and $x_{a}$
indicates the $a$-component of $x_{j}$. The $\phi_{m}$ term is the arrival
time of the unlensed signal unchanged by lensing effects.

Using Eq. 4, one can also calculate the differences between ($\mu_{j},T_{j}$)
and ($\mu_{k},T_{k}$), for the $j$-th image and $k$-th image, which we
identify as the relative magnification factor $\mu_{\rm rel}$ and
dimensionless time delay $\Delta T$.

$\mu_{\rm
rel}\equiv\frac{\mu_{j}}{\mu_{k}}=\frac{\operatorname{det}\left(\delta_{cd}-\frac{\partial^{2}\psi_{k}}{\partial
x_{c}\partial
x_{d}}\right)}{\operatorname{det}\left(\delta_{ab}-\frac{\partial^{2}\psi_{j}}{\partial
x_{a}\partial x_{b}}\right)}.$ (5) $\displaystyle\Delta T$
$\displaystyle\equiv\absolutevalue{T_{k}-T_{j}}$
$\displaystyle=\absolutevalue{\frac{1}{2}\left(|x_{k}-y|^{2}-|x_{j}-y|^{2}\right)+\psi_{k}-\psi_{j}}.$
(6)

These two quantities, which we call _lensing observables_ , can be obtained
numerically or analytically depending on the assumed lens model. Furthermore,
an actual time delay between the two images can be derived from Eq. 2 with a
given redshifted lens mass;

$\Delta t=\frac{4GM^{z}_{l}}{c^{3}}\Delta T.$ (7)

Note that we do not deal with the Morse phase because it has just three
discrete values $n_{j}=0/0.5/1$ when a Type I / Type II / Type III image
occurs, respectively (Takahashi & Nakamura, 2003; Dai & Venumadhav, 2017).

##  3 Lens model

We explore both relatively simple (spherically symmetric) lens models to
complex (non-axially symmetric) lens models to simulate various GW lensing
systems on a galactic scale. For some simple lens models, the lensing
observables are analytically obtained from $\boldsymbol{\theta_{l}}$. On the
other hand, it is non-trivial to calculate the lensing observables from
$\boldsymbol{\theta_{l}}$ for more realistic lens models. One needs to solve
the lens equations $y=x_{j}-\nabla\psi_{j}$ to obtain the image positions and
calculate the deflection potentials of each image in order to obtain the
lensing observables.

###  3.1 Axially symmetric lens

For axially symmetric lens models, we introduce three mass distributions as
described in Takahashi & Nakamura (2003).

Firstly, we consider the point mass (PM) model, which is one of the simplest
density profiles for small objects whose Einstein radii are much larger than
their sizes, such as compact objects and dwarf star clusters. The lens
parameter of an isolated point mass is the source position
($\boldsymbol{\theta_{l}}=y$). Within geometrical optics, this model always
generates double images. The lensing observables of the PM lens are given by

$\displaystyle\mu_{\mathrm{rel}}=\frac{2\left(y^{2}+2\right)}{y\sqrt{y^{2}+4}+y^{2}+2}-1,$
$\displaystyle\Delta
T=\frac{y\sqrt{y^{2}+4}}{2}+\ln\left(\frac{\sqrt{y^{2}+4}+y}{\sqrt{y^{2}+4}-y}\right).$
(8)

Note that the PM model is unsuitable for describing galactic halos, but we
consider this model to show the robustness of the inference technique
presented in this paper.

Secondly, the singular isothermal sphere (SIS) is the simplest model to
describe the distribution of matter in a gravitational system, such as a
galaxy or a galaxy cluster, which assumes the velocity dispersion of particles
in the system is constant at all radii. Like the PM model, the SIS has one
lens parameter ($\boldsymbol{\theta_{l}}=y$) and analytic lensing observables.
The relative magnification factor and the dimensionless time delay between two
images are written as

$\displaystyle\mu_{\rm rel}=\frac{1-y}{1+y}$ $\displaystyle\Delta T=2y.$ (9)

On the other hand, the SIS generates double images only if $y<1$.

Lastly, the Navarro-Frank-White (NFW) model is a density profile identified in
numerical N-body simulations, assuming $\Lambda$CDM cosmology (Navarro et al.,
1997). The density profile of the model is described by

$\rho(r)=\frac{\rho_{s}}{\frac{r}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{2}}~{},$
(10)

where $\rho_{s}$ and $r_{s}$ are the characteristic density and the scale
radius of the lens. Due to its dependence on cosmological parameters, the
lensing effects of the NFW model are determined by various parameters, such as
lens redshift ($z_{l}$), concentration ($c$), and a critical source position
($y_{cr}$) (Wright & Brainerd, 2000; Takahashi & Nakamura, 2003; Oguri, 2019).
Nevertheless, the dimensionless surface density
$\kappa_{s}=\rho_{s}r_{s}/\Sigma_{cr}$ includes information about the above
parameters. Thus, we define the lens parameters of the NFW model as
$\boldsymbol{\theta_{l}}=(y,\kappa_{s})$.

When $y<y_{cr}$, where $y=y_{cr}$ indicates the radial caustic of the lens,
three images are generated for the NFW model. The NFW model does not have an
analytic form for the lensing observables, unlike the above simple lenses.
Therefore, one needs to numerically derive the lensing observables using Eqs.
5 and 2, which requires solving the aforementioned lens equations. The
deflection angle $\alpha=\nabla\psi$ of the NFW model is given by

$\alpha_{\rm
NFW}(x)=\frac{4\kappa_{s}}{x}\times\left\\{\begin{array}[]{rl}\ln\frac{x}{2}+\frac{2}{\sqrt{x^{2}-1}}\arctan\sqrt{\frac{x-1}{x+1}}&(x>1)\\\
\ln\frac{x}{2}+\frac{2}{\sqrt{1-x^{2}}}\operatorname{arctanh}\sqrt{\frac{1-x}{1+x}}&(x<1)\\\
\ln\frac{1}{2}+1&(x=1)\end{array},\right.$ (11)

where the image position $x=|\vec{\xi}|/r_{s}$.

###  3.2 Non-axially symmetric lens

For more realistic scenarios, it is reasonable to use non-axially symmetric
lens models with more lens parameters, such as the minor to major axis ratio,
$q$. Normally, unlike simple lenses, lensing observables of elliptical lens
models cannot be analytically calculated from $\boldsymbol{\theta_{l}}$.
Similarly to the NFW model, the lensing observables should be numerically
calculated from image positions and deflection angles. However, the deflection
angle depends on the position angle with regard to the center of the halo as
well as the displacement to the image position. Thus, it is more demanding to
solve the lens equations for non-axially symmetric lens models.

We introduce three elliptical mass distribution models generally used in EM
strong lensing analyses (e.g., Treu (2010); Tessore & Metcalf (2015); Oguri
(2019); Nightingale et al. (2019); Shajib et al. (2021)). First, the singular
isothermal ellipsoid (SIE) is the simplest elliptical lens model (Kormann et
al., 1994). The SIE model assumes that the SIS model still approximately
describes the mass distribution of the galaxy or cluster of galaxies, but with
some ellipticity. Thus, the lensing effects of the SIE model are described by
two parameters, $\boldsymbol{\theta_{l}}=(q,y)$. Double or quadruple images
are generated depending on the combination of the two lens parameters. To
obtain lensing observables for the SIE model, one needs to solve the lens
equation using the deflection angle as follows,

$\begin{split}\alpha_{\rm
SIE}(x_{1},x_{2})=&\sqrt{\frac{q}{1-q^{2}}}\left[\operatorname{arsinh}\left(\frac{\sqrt{1-q^{2}}}{q}\frac{x_{1}}{\sqrt{x^{2}_{1}+x^{2}_{2}}}\right)\right.\\\
&+\left.\arcsin\left(\sqrt{1-q^{2}}\frac{x_{2}}{\sqrt{x^{2}_{1}+x^{2}_{2}}}\right)\right]~{},\end{split}$
(12)

where $x_{1}$ and $x_{2}$ are the $x$\- and $y$\- components of the image
position normalized by the Einstein radius $\theta_{\rm E}$ of the lens.

A more general lens model can no longer be isothermal. The singular power-law
ellipsoidal mass distribution (SPEMD) is one of the generally used elliptical
and non-isothermal density profiles (Barkana, 1998), which has an additional
lens parameter, the slope ($t$) of halo density. Even though lensing effects
due to the SPEMD must be calculated using numerical methods, Tessore & Metcalf
(2015) derived a pseudo-analytic form of the deflection angle of SPEMD with a
complex formalism. Assuming that $\alpha=\alpha_{x_{1}}+i\alpha_{x_{2}}$, the
deflection angle at the image position $(x_{1},x_{2})$ of the SPEMD model is
given by

$\begin{split}\alpha_{\rm
SPEMD}(x_{1},x_{2})={}&\frac{2\sqrt{q}}{1+q}\left(\frac{\sqrt{q}}{\sqrt{q^{2}x^{2}_{1}+x^{2}_{2}}}\right)^{t-1}\\\
&\times\mathrm{e}^{\mathrm{i}\varphi}{}_{2}F_{1}\left(1,\frac{t}{2};2-\frac{t}{2};-\frac{1-q}{1+q}\mathrm{e}^{\mathrm{i}2\varphi}\right)~{},\end{split}$
(13)

where $\varphi=\arctan(x_{2}/qx_{1})$ is the elliptical angle, and
${}_{2}F_{1}$ is the Gaussian hypergeometric function. The SPEMD model can be
reduced to an SIE when the density slope $t=1$. When $t\geq 1$, double and
quadruple images can be generated like for the SIE. Triple and quintuple
images, on the other hand, can be generated when $t<1$. This is because for a
SPEMD lens with $t<1$, the critical lines and caustics are more complicated
than $t\geq 1$ case (See Fig.2 in Tessore & Metcalf (2015)). For both cases,
the Einstein radius of a halo normalizes the image position and deflection
angle.

Finally, we consider the elliptical Navarro-Frank-White (ENFW) model, which is
a modification of the standard spherical NFW profile described in Sec. 3.1.
One can compute the deflection angles of elliptical lens models by introducing
an elliptical coordinate system. We use
$x_{\epsilon}=\sqrt{x_{1\epsilon}^{2}+x_{2\epsilon}^{2}}=\sqrt{a_{1\epsilon}x_{1}^{2}+a_{2\epsilon}x_{2}^{2}}$
to express the deflection angle of ENFW, as in (Golse & Kneib, 2002),

$\begin{split}\alpha_{\rm ENFW}(x_{1},x_{2})={}&\alpha_{\rm
NFW}\left(x_{\epsilon}\right)\sqrt{a_{1\epsilon}}\cos\phi_{\epsilon}\\\
&+\alpha_{\rm
NFW}\left(x_{\epsilon}\right)\sqrt{a_{2\epsilon}}\sin\phi_{\epsilon}~{},\end{split}$
(14)

where $\phi_{\epsilon}=\arctan{(x_{2\epsilon}/x_{1\epsilon})}$ is the
elliptical angle, same as $\varphi$ in Eq. 13. For consistency with the other
elliptical lens models above, we define $a_{1\epsilon}=q$ and
$a_{2\epsilon}=1/q$. Depending on whether the source is located in caustics or
not, triple or quintuple images can be generated.

For all lens models in this work, we calculate $\mu_{\rm rel}$ and $\Delta T$
from physically possible $\boldsymbol{\theta_{l}}$ values for the six lens
models. We then save, in a _lensing observable bank_ , ordered triplets of
lens parameters and lensing observables
($\Theta=\\{\boldsymbol{\theta_{l}},\mu_{\rm rel},\Delta T\\}$) satisfying the
multiple-image condition (i.e., the necessary condition for strong lensing to
occur).

##  4 Strong lensing cross section

It is informative to investigate how efficiently various lens models create
multiple GW images when the same source redshift ($z_{s}$), lens redshift
($z_{l}$) and lens mass ($M_{l}$) are given. The strong lensing cross-section
($\sigma$) can be defined as the area where multiple images form, magnified
above a certain magnification threshold, $\mu_{\rm th}$ (Schneider et al.,
2006; Fedeli et al., 2006; Robertson et al., 2020). For a simple point-like
lens, a disk with a radius equal to the Einstein radius ($r_{E}$) of the lens
represents the strong lensing cross-section. However, the strong lensing
cross-sections of the complex lens models considered in Sec. 3 cannot be
analytically calculated.

Since we can regard a GW source as a point source, the strong lensing cross-
section is simply the area on the source plane enclosed by the caustics of the
lens, where multiple images are created. In this work, however, we impose an
additional condition on the definition of the lensing cross-section for
realistic detections, requiring that all images are magnified, i.e.
$\mu_{j}>1.0,\forall\,j$. When a source is placed within an area $A_{s}$
satisfying the above conditions, the lensing cross-section can be written as

$\sigma_{1.0}=\frac{\xi_{0}D_{S}}{D_{L}}\int_{A_{\mathrm{s}}}\mathrm{d}^{2}y~{},$
(15)

where $y$ is the dimensionless unit on the source plane (See Eq.1) (Fedeli et
al., 2006).

(a) $\sigma_{\rm 1.0}$ where at least two images are created and magnified

(b) $\sigma_{\rm 1.0}$ where at least four images are created and magnified

Figure 2: Left panel: Strong lensing cross section as a function of
$z_{l}/z_{s}$, where at least two GW images are magnified (i.e., $\mu>1.0$)
for (E)NFW (red), SPEMD with $t=0.5,~{}1.5,~{}1.99$ (orange, green, blue), SIE
(yellow), and PM (black) lens models. Coloured solid lines indicate the
$\sigma_{\rm 1.0}$ of spherical lenses ($q=1.0$), and dashed lines do
$\sigma_{\rm 1.0}$ of elliptical lenses with $q=0.1$. Generally, it is more
probable for denser lens models to create magnified images. Furthermore, a
SPEMD lens closer to the observer has larger $\sigma_{\rm 1.0}$, while the
value of $\sigma_{1.0}$ for an (E)NFW lens decreases when the lens is closer
to the observer because $\sigma_{1.0}$ depends on the critical density
$\Sigma_{\rm cr}$. Moreover, a lens with high ellipticity is less efficient in
generating magnified images when $t\geq 1$, while a lens with $t<1$ has the
contrary aspect. Right panel: Same as the left panel, but at least four lensed
GW signals are magnified. Coloured dotted lines indicate $\sigma_{1.0}$ of
elliptical lenses with $q=0.5$. We set the maximum value of the axis ratio as
$q=0.9$ because a lens with $q>0.9$ has a diminutive inner caustic area. The
function of $\sigma_{1.0}$ are similar to those of the double-image cases
shown in Fig. 2(a), but the value of $\sigma_{1.0}$ changes more steeply with
the ellipticity.

We adopt the ray-shooting technique to calculate the area $A_{s}$. Firstly, we
calculate the largest caustic (cut)111The cut is an area where multiple lensed
images can be created when a source is located within it, even though the area
is not enclosed by a caustic. A cut will convert into a caustic when the
singularity at the centre of a lens is removed. area $B_{s}$ on the
dimensionless source plane, given the lens models and their $\theta_{l}$
except for $y$. Note that $y$ determines not the lens properties but the
lensing effects on propagating GWs. By using the ray-shooting technique, the
corresponding image positions and magnification factors can be calculated for
all source positions within $B_{s}$. The ratio between the area filled with
the source positions satisfying $\mu(y)>1.0$ and the entire $B_{s}$ will be
the dimensionless cross-section of interest, $A_{s}$.

Next, we adopt the Einstein radius $r_{E}$ (for PM and SPEMD-like models) and
the scale radius $r_{s}$ (for E(NFW) models) as the normalisation constant
$\xi_{0}$ that converts dimensionless values into real values. The Einstein
radius of a PM lens is a function of Einstein mass, $M_{E}$, but the Einstein
radii of the SIS, SIE, and SPEMD lenses depend explicitly on the velocity
dispersion, $\sigma_{v}$, of a galaxy rather than the $M_{E}$ of the galaxy.
For this reason, we assume a plausible galaxy that has an Einstein mass
$M_{E}\sim 10^{12}M_{\odot}$ and velocity dispersion $\sigma_{v}\sim 200$
km/s, based on Zahid et al. (2016).

For E(NFW) models, we assume another, different galactic halo whose mass
within the scale radius is $M(r_{s})=10^{13}M_{\odot}$. This is because the
E(NFW) lenses with $M(r_{s})=10^{12}M_{\odot}$ create multiple images only
vanishingly rarely due to their inefficiency at image splitting.

The strong lensing cross-sections with at least two magnified images are shown
for various lens models in Fig. 2(a). We include the six lens models
considered. Among these lens models, elliptical models can have a double-image
system for any value of $q$, but we consider only $q\geq 0.1$ for
practicality. For the SPEMD, we input three different density slope parameters
($t=0.5,1.5,1.99$) to see how the density of a lens affects the lensing cross-
section.

The SPEMD models, including SIS and SIE at a lower $z_{l}$ tend to have a
higher $\sigma_{1.0}$ value. Furthermore, $\sigma_{1.0}$ increases with the
slope $t$ of the density, which means that a denser lens model is more
efficient at generating multiple signals. Note that the function of
$\sigma_{1.0}$ of the SPEMD with $t=1.99$ is very similar to that of PM
because when $t\rightarrow 2$, the SPEMD lens becomes a point mass (Tessore &
Metcalf, 2015). Moreover, a SPEMD lens with $t<1$ is more efficient in
splitting and magnifying GW signals when the lens has high ellipticity (low
$q$). A SPEMD lens with $t\geq 1$ and a higher ellipticity, on the other hand,
has less efficiency in generating multiple magnified signals.

Non-axially symmetric lenses have multiple caustics, in which a propagating GW
signal can be split into quadruple or more signals. In this regard, we also
calculate $\sigma_{1.0}$ for lens systems which can create quadruple or more
lensed GW signals.

Fig. 2(b) shows strong lensing cross-sections $\sigma_{1.0}$, generating and
magnifying at least four lensed GW signals. The shapes of the plots are the
same as the plots in Fig. 2(a) since the unit areas of each lens model are the
same as $\theta_{E}^{2}$. However, the rate of change of $\sigma_{1.0}$
becomes more pronounced as the ellipticity varies. Also, $\sigma_{1.0}$ has
the maximum value at $q=0.5$, attributed to the limited region capable of
producing magnified images when $q\rightarrow 0.1$, and the shrinking of the
caustic area as $q\rightarrow 1$. Note that we exclude the SPEMD with $t=1.99$
and ENFW models as their inner caustics are too small to create quadruple or
more images in a galaxy-scale gravitational lensing scenario.

##  5 Inferring the parameters of a galaxy lens

We simulate nine GW lensing systems based on the lens models described in the
previous section. For simplicity, we inject the same intrinsic unlensed signal
$h_{\rm{u}}$ (GW150914-like event with $D_{L}=3~{}\rm{Gpc}$) and redshifted
lens mass value $M^{z}_{l}=10^{12}\rm{M}_{\odot}$ for all lensing systems
except for NFW and ENFW lens cases. For the two lens models, the strong
lensing cross sections ($\sigma_{1.0}$) with the above setting are too narrow
and thus, we inject a GW150914-like event with $D_{L}=10~{}\rm{Gpc}$ and
redshifted lens mass value $M^{z}_{l}=10^{13}\rm{M}_{\odot}$.

We use IMRPhenomXPHM (Pratten et al., 2021) to generate $h_{\rm{u}}$ and set
possible source-lens systems, from double-image to quintuple-image systems
where the number of the image is determined by the model-dependent parameters.
Considering the detectability of the LIGO–Virgo detectors (Abbott et al.,
2020), we assume that all lensed GW signals ($h_{l}(f)$) are detected at
network signal-to-noise ratios (SNR) $\rho_{\rm net}>8$ except for NFW and
ENFW cases 222Type III images of NFW and ENFW mostly have low SNRs.
Nevertheless, we set the SNR of the two type III images $\rho>4$, which can be
found by sub-threshold search (Li et al., 2023).

ID | Model | Injected $\boldsymbol{\theta_{l}}$ | No. of images
---|---|---|---
1 | PM | $y=0.22$ | Double
2 | SIS | $y=0.22$ | Double
3 | NFW | $(y,~{}\kappa_{s})=(0.004,~{}0.167)$ | Triple
4 | SIE | $(q,~{}y)=(0.85,~{}0.22)$ | Double
5 | SIE | $(q,~{}y)=(0.63,~{}0.22)$ | Quadruple
6 | SPEMD | $(q,~{}y,~{}t)=(0.85,~{}0.22,~{}1.5)$ | Double
7 | SPEMD | $(q,~{}y,~{}t)=(0.58,~{}0.17,~{}0.75)$ | Qunituple
8 | ENFW | $(q,~{}y,~{}\kappa_{s})=(0.6,~{}0.011,~{}0.167)$ | Triple
9 | ENFW | $(q,~{}y,~{}\kappa_{s})=(0.82,~{}0.008,~{}0.167)$ | Quintuple
Table 1: Injected $\boldsymbol{\theta_{l}}$ values used to generate the nine
GW lensing systems considered in this work. The parameters are randomly chosen
from a pool of $\boldsymbol{\theta_{l}}$ that can generate the specified
multiple (double to quintuple) images . See Sec. 3 for the definition of lens
parameters for each lens model.

For each system, we choose the injected $\boldsymbol{\theta_{l}}$ shown in
Table. 1 to generate multiple lensed signals. The SNR ranges of lensed GW
signals are $\rho=(11,16)$ for Systems 1,2 and 4$\sim$7, and $\rho=(4,10)$ for
Systems 3, 8 and 9.

Next, we infer the relative magnification factors $\mu_{\rm{rel}}$ and time
delays $\Delta t$ among the lensed GW signals by performing joint parameter
estimation on pairs of the first lensed signal ($h_{l,1}$) and the $j$-th
lensed signal ($h_{l,j}$). If the system has $n$ lensed GW signals, the
relative magnification factors between $h_{l,1}$ and $h_{l,j}$ can be written
as

$\mu_{\rm{rel},1\it{j}}=\frac{\mu_{j}}{\mu_{1}}=\left(\frac{D_{L,1}}{D_{L,j}}\right)^{2}\text{for}\quad
j=2,...,n-1~{},$ (16)

where $D_{L,j}$ is the apparent luminosity distance inferred from the $j$-th
lensed signal. Also, the actual time delays can be written as

$\Delta t_{1j}=|t_{c,j}-t_{c,1}|\quad\text{for}\quad j=2,...,n-1~{},$ (17)

where $t_{c,j}$ is the apparent merger time inferred from the $j$-th lensed
signal.

Figure 3: The posteriors of $\mu_{\rm rel}$ and $\Delta t$ from the joint
parameter estimation for two lensed GW signals. An isolated point mass with
$(y,M^{z}_{l})=(0.22,10^{12}M_{\odot})$ is assumed. Gold lines indicate the
true values. Vertical dashed lines mark the upper and lower bound of 90%
credible intervals. The lensing observables are well constrained, and the
posteriors peak at the true values ($\mu_{\rm rel}\simeq 0.64$ and $\Delta
t\simeq 8645972.23\rm{s}$).

To conduct joint parameter estimation, we use the golum pipeline (Janquart et
al., 2021) which in turn makes use of the Bayesian inference tools bilby
(Ashton et al., 2019; Romero-Shaw et al., 2020) and dynesty (Speagle, 2020)
for nested sampling. Simulated lensed GW signals are injected into Gaussian
noise characterized by the representative LIGO-Virgo O4 power spectral density
(PSD) (Abbott et al., 2020). We assume uniform prior distributions for the
relative magnification factor $\mu_{\rm rel}$ and actual time delay $\Delta
t$. For the source parameters, we adopt the prior distributions for precessing
binary black holes in Romero-Shaw et al. (2020).

Fig. 3 shows the posteriors of the $\mu_{\rm rel}$ and $\Delta t$ from the
lensing system assuming the SIS model (see Table. 1). The posteriors peak at
the injected values with well-constrained credible intervals. Although the
peaks for the other lensing system cases cannot be precisely at the injected
values due to the effects of noise, we confirm that the peaks of the posterior
become nearer to the injected values when we use zero noise. We will use these
posteriors of $\mu_{\rm rel}$ and $\Delta t$ as constraints to infer the
$\boldsymbol{\theta_{l}}$.

Except for a few simple models like PM and SIS lenses (See Eqs. 3.1, 3.1), the
$\mu_{\rm rel}$ and $\Delta t$ do not have one-to-one relationships with
$\boldsymbol{\theta_{l}}$. Furthermore, the inferred $\mu_{\rm rel}$ takes the
form of a probability distribution rather than a point estimate with a
specific value, which may hinder us from accurately recovering the lens
parameters. In this sense, we can use rejection sampling technique (MacKay &
JC, 2003) to infer the $\boldsymbol{\theta_{l}}$ since we can specifically
determine the value of $\mu_{\rm rel}$ from a given $\boldsymbol{\theta_{l}}$
and have the probability distribution of $\mu_{\rm rel}$.

We use the lensing observable bank mentioned in Sec 3, which contains ordered
triplets ($\Theta=\\{\boldsymbol{\theta_{l}},\mu_{\rm rel},\Delta T\\}$), in
which $\boldsymbol{\theta_{l}}$ create at least two lensed GW signals. In this
study, we use lenstronomy (Birrer & Amara, 2018) to calculate $\mu_{\rm rel}$
and $\Delta T$ from a given $\boldsymbol{\theta_{l}}$ for the rejection
sampling. Thus, each sampled $\mu_{\rm rel}$ value is associated with a
specific set of lensing parameter values, $\boldsymbol{\theta_{l}}$. Fig. 4
shows an example case for the posterior distribution $P(\mu_{\rm rel})$, which
includes accepted (green) and rejected (red) points as a result of the
rejection sampling.

Figure 4: Illustration of rejection sampling. The solid blue line is the
proposal distribution (uniform) to sample points, and the solid black line is
the target distribution $P(\mu_{\rm rel})$ equal to the posterior of obtained
$\mu_{\rm rel}$ in Figure 3. The solid yellow line is a new proposal
distribution made by multiplying a constant $M$ by the initial uniform
distribution $Q(\mu_{\rm rel})$, which satisfies with $P(\mu_{\rm rel})\leq
MQ(\mu_{\rm rel})$. Based on the acceptance criterion, green and magenta
points indicate the accepted and rejected $\boldsymbol{\theta_{l}}$,
respectively.

Finally, the $\boldsymbol{\theta_{L}}$ points from the accepted triplets
(green dots) may be used to represent the posteriors of these parameters.
Also, the associated $\Delta T$ points may be converted to the corresponding
$M^{z}_{l}$ samples by using Eq. 7.

It is noteworthy that the computation time of the rejection sampling technique
is less than 1 CPU hour for a double-image system (first and second image
pair). For triple to quintuple-image systems, using rejection sampling would
require a few CPU hours as it involves conducting rejection sampling for the
first and $j^{\rm th}$ image pair only after obtaining the result for the
first and $i^{\rm th}$ pair ($i<j$).

##  6 Results

###  6.1 Axially symmetric lens cases

To check the validity of the rejection sampling technique, we first compare
the $\boldsymbol{\theta_{l}}$ posteriors inferred from rejection sampling to
analytically obtained posteriors using Eqs. 3.1 and 3.1 for PM and SIS
respectively. For these two lens models, the relative magnification factors
$\mu_{\rm rel}$ are univariate functions of the lens parameter
$\boldsymbol{\theta_{l}}=y$. Fig. 5 shows that the recovered posteriors of
$M^{z}_{l}$ and $y$ for System 1 (PM) and System 2 (SIS) using the above two
methods peak at the injected values with great accuracy, and they are almost
identical. These results demonstrate that rejection sampling can be used to
recover the lens parameters from the lensing observables. Note that this does
not mean rejection sampling can successfully recover all
$\boldsymbol{\theta_{l}}$ regardless of the degeneracies between them.

Figure 5: Posteriors of $\boldsymbol{\theta_{l}}$ for Systems 1 (PM double)
and 2 (SIS double) in Table 1 inferred from the relative magnification factor
$\mu_{\rm rel}$ using rejection sampling (blue and violet). The solid black
lines indicate that the injected values and the vertical dashed lines are 90%
symmetric credible intervals for each parameter. The $\boldsymbol{\theta_{l}}$
posteriors analytically obtained using Eq. 3.1 and Eq. 3.1 are shown with
dashed grey and black lines, respectively. In both cases, the posteriors of
$y$ and $M^{z}_{L}$ inferred using rejection sampling and analytically are
almost identical, which demonstrates the validity of the rejection sampling
technique.

Among the considered spherically symmetric lens models, the NFW model does not
have an analytic expression for the lensing observables. Thus, we use the
rejection sampling technique for System 3 (NFW triple). After identifying the
estimated $\mu_{\rm rel}$ posteriors with the detected triple images as a
target distribution, we draw $\mu_{\rm rel}$ samples from the lensing
observable bank of the NFW to conduct rejection sampling.

Fig. 6 shows the $\boldsymbol{\theta_{l}}$ and $M^{z}_{l}$ posteriors for the
NFW case. Even though injected values are contained within the 90% credible
intervals, all posteriors are poorly constrained. Even $\kappa_{s}$ is not
converged well and does not change from the prior. This is because many of the
$(y,\kappa_{s})$ pairs in the prior pool, which can create multiple images,
have similar $\mu_{\rm rel}$ values under galaxy-scale lens scenarios (i.e.,
only very small values of $y$ can generate GW strong lensing due to low
surface density). The obtained $\mu_{\rm rel}$ posterior is not well enough
constrained to break this degeneracy between $y$ and $\kappa_{s}$.

Nevertheless, we can better recover the lens parameters of NFW when $y$ and
$\kappa_{s}$ have large values (i.e., as will occur in galaxy cluster scale
lens scenarios).

Figure 6: Corner plot showing the posteriors of the lens parameters
$\boldsymbol{\theta_{l}}=(y,\kappa_{s})$ and redshifted lens mass $M^{z}_{l}$
for System 3. The solid black lines indicate the injected values, and vertical
dashed lines mark the upper and lower bound of 90% credible intervals. Each 1D
posterior of $y$ and $\kappa_{s}$ is poorly constrained because the degeneracy
between the two parameters is significant. Correspondingly, the $M^{z}_{l}$
does not peak at the true value.

###  6.2 Non-axially symmetric lens cases

The ellipticity of a lens makes the lens model more complicated and
degeneracies between $\boldsymbol{\theta_{l}}$ stronger. Thus, the
$\boldsymbol{\theta_{l}}$ and $M^{z}_{l}$ posteriors for non-axially symmetric
lens models are more poorly constrained compared to those of spherically
symmetric lens models.

(a) System 4 (SIE double)

(b) System 6 (SPEMD double)

(c) System 8 (ENFW triple)

Figure 7: Corner plots showing the posteriors of $\boldsymbol{\theta_{l}}$ and
$M^{z}_{l}$ for Systems 4, 6, and 8. Solid black lines indicate the injected
values and vertical dashed lines mark the upper and lower bound of 90%
credible intervals. Even though all injections are included within the 90%
credible intervals, the $\boldsymbol{\theta_{l}}$ posteriors do not converge
well. Meanwhile, the $M^{z}_{l}$ posteriors are not highly biased from the
injection because the mass strongly depends on the time delays between images,
which can be constrained much better than relative magnification factors.
Overall, there are limitations on accurately estimating dark halo properties
with only double or triple-image systems.

Fig. 7 shows the estimated posteriors of $\boldsymbol{\theta_{l}}$ and
$M^{z}_{l}$ for Systems 4, 6, and 8 in Table 1 where double images (SIE and
SPEMD) and triple images (ENFW) are detected. Even though, for all three
models, the true values of the parameters lie within the 2$\sigma$ credible
regions of the joint 2-D posteriors of the lens parameters, the marginal 1-D
posterior of each parameter does not peak at the injected value and is not
well constrained. This is because complex elliptical lenses have more lens
parameters than spherically symmetric lenses, which results in stronger
degeneracies between the parameters of $\boldsymbol{\theta_{l}}$. Therefore,
the properties of an elliptical dark galactic halo cannot be accurately
estimated when the GW lensing system, having the halo, generates only double
or triple images considering our lens model assumptions.

When more than double images are created and detected, however, the
degeneracies between $\boldsymbol{\theta_{l}}$ of elliptical lens models can
be further reduced thanks to there being more constraints. In Fig. 8, we
present the posteriors of $\boldsymbol{\theta_{l}}$ and $M^{z}_{l}$ for
Systems 5, 7, and 9. These systems produce quadruple or quintuple GW images,
but it is plausible that we may only observe a subset of the complete GW image
ensemble. In light of this, we present the outcomes of two scenarios: 1) the
detection of only pairs of GW images and (grey) 2) the detection of all GW
images (coloured).

Regarding the degree of convergence to the injected values, the posteriors
obtained from more than double images show better results. This is especially
the case for SPEMD, where all $\boldsymbol{\theta_{l}}$ posteriors converge
well to the injected values. For the ENFW case, on the other hand, the
degeneracy between $y$ and $\kappa_{s}$ is significant, similar to the NFW
lens case. Thus, for non-axially symmetric models, the posteriors of
$\boldsymbol{\theta_{l}}$ and $M^{z}_{l}$ are generally not well constrained
from GW information alone.

(a) System 5 (SIE quadruple)

(b) System 7 (SPEMD quadruple)

(c) System 9 (ENFW quintuple)

Figure 8: Same configurations as in Fig. 7, but showing the posteriors of
$\boldsymbol{\theta_{l}}$ and $M^{z}_{l}$ for Systems 5, 7, and 9. Histograms
filled with colour represent the posteriors obtained from all sets of
quadruple (or quintuple) GW images, whereas histograms depicted as solid grey
lines illustrate the posteriors obtained exclusively from double images, which
constitute a subset of the entire image collection. Since more constraints
(i.e., more $\mu_{\rm rel}$ and $\Delta t$) are provided, the inferred
posteriors for $\boldsymbol{\theta_{l}}$ are better constrained than those in
the case where a double image is detected. Specifically, for all models, the
widths of the posteriors for $q$ are narrowed down, which can be helpful to
localize the lens galaxy in the inferred sky area.

###  6.3 Non-axially symmetric lens + indicative EM observation cases

Figure 9: Same configurations as in Fig. 8, but with the assumption that
quadruple (quintuple) image systems are detected with the indicative
electromagnetic observations and thereby, axis ratios and image positions are
confined with $\pm$ 5% and $0.1^{"}$ error ranges. For easier comparison to
the results without constraints from indicative EM observations, we plot the
posteriors in Fig. 8 with solid grey lines. The inferred posteriors of
$\boldsymbol{\theta_{l}}$ and $M^{z}_{l}$ for all non-axially symmetric lens
models peak around the injected values (solid black lines) within 90% credible
intervals. Notably, the recovery of $M^{z}_{l}$ shows improvement, enabling
more precise calculations of $M_{l}$ with the estimated $z_{l}$ of the host
galaxy.

In this sub-section, we consider how extra constraints from EM observations
can reduce the degeneracies between parameters further, so that $M^{z}_{l}$
can be better constrained. For instance, the axis ratio of the lens galaxy and
EM image positions of the lensed host galaxy can be directly estimated if one
explicitly identifies the lens galaxy. However, it is challenging to search
for the lens galaxy in the GW localization area, which may contain numerous
galaxies.

In this context, Yu et al. (2020); Hannuksela et al. (2020); Wempe et al.
(2022) propose that one can localize the host galaxy-lens galaxy system using
LSST surveys (Abell et al., 2009) subject to several assumptions, such as the
detection of both lensed GW signals and EM signals from the host
galaxy333Generally, the host galaxy of a GW source can also be lensed by
intervening dark galactic halos.. Based on the magnification factors and time
delays obtained from the GW and EM observations, one can reduce the number of
lens galaxy candidates within the sky area determined by the GW localization.

In particular, Hannuksela et al. (2020) argues that detecting quadruple GW
images can effectively narrow down the sky area given the current LIGO-Virgo
detectors. Moreover, Fig. 8 shows that the axis ratio of a lens galaxy can be
more precisely determined when four or more gravitational wave (GW) images are
detected. This GW information can be utilized to identify the lens galaxy with
a well-constrained axis ratio.

Therefore, we assume here that the axis ratios of the lens galaxies and image
positions can be confined to within an error range of $\pm$5% and $0.1^{"}$
respectively, once each lens has been specified based on the observations of
quadruple or more GW images. Fig. 9 shows the posteriors of
$\boldsymbol{\theta_{l}}$ and $M^{z}_{l}$ of non-symmetric lens galaxies in
each GW lensing system (See Table. 1) derived under this assumption. Thanks to
the reduced degeneracies between $\boldsymbol{\theta_{l}}$, all their
posteriors are now better constrained and peak at the injected values. Above
all, the injected $M^{z}_{l}$ are successfully recovered. One can calculate
the lens mass $M_{l}$ with the inferred redshifted lens mass $M^{z}_{l}$ and
estimated redshift $z_{l}$ from the corresponding EM observation.

##  7 Discussion

Dark galactic halos are keys to understanding galaxy evolution, formation of
large-scale structures, and constraining cosmological parameters. Precisely
measuring their masses and structures is essential. Gravitational lensing of
GWs can significantly contribute to GW astronomy as increasing numbers of
lensed GWs are observed in the near future by the next generations of
detectors.

One can complement the 3D map of the mass distribution from lensed GW signals
by constraining the halo mass function. It is important to note that, for NFW
and ENFW models, we can constrain both the lens mass $M^{z}_{l}$ and the
dimensionless surface density $\kappa_{s}$, which are directly related to
cosmological parameters. The relationship between the two parameters can also
be investigated to constrain the concentration ($c$) and critical density
($\rho_{c}$).

In this work, we have proposed how to estimate lens and source parameters from
the lensed GW signals using various lens models. We assume that all lensed GW
signals are detected and the correct lens model is used to analyze the
signals. Our results demonstrate that we can retrieve lens parameters from the
information contained in the lensed GW signals, even where these parameters do
not have trivial relationships with each other.

Furthermore, we can make use of EM observations, which become viable when
quadruple or more GW image systems are detected, to confine the axis ratios of
the lenses, thereby enhancing the accuracy of their mass measurements.
Although only a BBH source is discussed in this work, we can also consider BNS
and NSBH signals, where we could, in principle, observe both GW and EM lensing
simultaneously, which facilitates localizing the host galaxy (Smith et al.,
2023; Ma et al., 2023).

However, we recognise that the situation could be more complicated in real GW
detection scenarios. Firstly, we may observe only a subset of all the lensed
GW signals. For example, only two lensed signals might be observed when an SIE
galaxy generates four lensed signals from a given GW source. Assuming that the
GW detectors keep operating and our lens model is accurate, there are two
scenarios: i) The two non-detected signals have not yet arrived because of
their longer time delays. ii) The two non-detected signals have relatively
lower SNR than the two detected signals. Therefore, one should perform
rejection sampling assuming both double and quadruple-image triplets to
investigate which scenario is more plausible and which number of images is
preferred.

Next, we cannot instantly determine the lens density profile when we detect a
set of lensed GW signals from a given GW source. However, by reversely
incorporating the rejection sampling technique into joint parameter estimation
steps, one can conduct hypothesis testing to determine which model is more
preferred to describe the lensed GW signals. This process will also be helpful
in the search for strongly lensed GW signals because Çalışkan et al. (2023a);
Janquart et al. (2023) argue that lens model-dependent joint parameter
estimations are necessary to overcome false alarms and improve the search.

Finally, we note that the substructure within a lens galaxy, such as star
clusters and dark matter subhalos, could induce single or multiple
microlensing signals to occur in a strongly-lensed GW event (Diego et al.,
2019; Mishra et al., 2021; Seo et al., 2022; Yeung et al., 2023). As a result,
additional lens parameters would be required to describe the lensed GW
signals, also necessitating the application of wave optics to fully understand
the entire lensing effects. 444When contemplating a space-based
interferometer, which may enable us to observe significantly higher SNR
events, the enhanced degeneracies among the lens parameters can in principle
be resolved within the wave-optics regime. (Tambalo et al., 2022;
Çal$\iota$şkan et al., 2023; Çalışkan et al., 2023b) Since the rejection
sampling technique is highly dependent on the magnification factors and time
delays, ignoring these additional microlensing effects could result in
incorrect estimation of the lens parameters. Thus, the development of
parameter estimation techniques that incorporate macrolens-microlens lensing
systems is necessary. In future work we will expand our study to address the
above problems.

## Acknowledgements

E. S. is supported by grants from the Research Grants Council of the Hong Kong
(Project No. CUHK 24304317), the Croucher Foundation of Hong Kong, the
Research Committee of the Chinese University of Hong Kong, and the College of
Science and Engineering of the University of Glasgow. T. G. F. is supported by
FWO (Research Foundation Flanders) Research Project nr. G086722N and FWO
International Research Infrastructure nr. I002123N. M. H. is supported by the
Science and Technology Facilities Council (Grant Ref. ST/V005634/1). The
authors are grateful for computational resources provided by the LIGO
Laboratory and supported by the National Science Foundation Grants PHY-0757058
and PHY-0823459

## References

  * Abbott et al. (2020) Abbott, B. P., Abbott, R., Abbott, T., et al. 2020, Living reviews in relativity, 23, 1
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, The Astrophysical Journal, 923, 14
  * Abbott et al. (2023) Abbott, R., Abe, H., Acernese, F., et al. 2023, arXiv preprint arXiv:2304.08393
  * Abell et al. (2009) Abell, P. A., Allison, J., Anderson, S. F., et al. 2009
  * Ade et al. (2016) Ade, P. A., Aghanim, N., Arnaud, M., et al. 2016, Astronomy & Astrophysics, 594, A13
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, The Astrophysical Journal Supplement Series, 241, 27
  * Barkana (1998) Barkana, R. 1998, The Astrophysical Journal, 502, 531
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, Monthly Notices of the Royal Astronomical Society, 488, 3143
  * Behroozi et al. (2010) Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, The Astrophysical Journal, 717, 379
  * Biesiada et al. (2014) Biesiada, M., Ding, X., Piórkowska, A., & Zhu, Z.-H. 2014, Journal of Cosmology and Astroparticle Physics, 2014, 080
  * Birrer & Amara (2018) Birrer, S., & Amara, A. 2018, Physics of the Dark Universe, 22, 189
  * Broadhurst et al. (1994) Broadhurst, T. J., Taylor, A., & Peacock, J. 1994, Mapping cluster mass distributions via gravitational lensing of background galaxies
  * Çal$\iota$şkan et al. (2023) Çal$\iota$şkan, M., Ji, L., Cotesta, R., et al. 2023, Physical Review D, 107, 043029
  * Çalışkan et al. (2023a) Çalışkan, M., Ezquiaga, J. M., Hannuksela, O. A., & Holz, D. E. 2023a, Physical Review D, 107, 063023
  * Çalışkan et al. (2023b) Çalışkan, M., Kumar, N. A., Ji, L., et al. 2023b, arXiv preprint arXiv:2307.06990
  * Cao et al. (2014) Cao, Z., Li, L.-F., & Wang, Y. 2014, Physical Review D, 90, 062003
  * 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
  * Dai & Venumadhav (2017) Dai, L., & Venumadhav, T. 2017, arXiv preprint arXiv:1702.04724
  * Diego et al. (2019) Diego, J. M., Hannuksela, O. A., Kelly, P. L., et al. 2019, Astronomy & Astrophysics, 627, A130
  * Ding et al. (2015) Ding, X., Biesiada, M., & Zhu, Z.-H. 2015, Journal of Cosmology and Astroparticle Physics, 2015, 006
  * Elahi et al. (2018) Elahi, P. J., Power, C., Lagos, C. d. P., Poulton, R., & Robotham, A. S. 2018, Using velocity dispersion to estimate halo mass: Is the Local Group in tension with $\Lambda$CDM?, Oxford University Press
  * Fanidakis et al. (2012) Fanidakis, N., Baugh, C., Benson, A., et al. 2012, The evolution of active galactic nuclei across cosmic time: what is downsizing?, Oxford University Press
  * Fedeli et al. (2006) Fedeli, C., Meneghetti, M., Bartelmann, M., Dolag, K., & Moscardini, L. 2006, Astronomy & Astrophysics, 447, 419
  * Golse & Kneib (2002) Golse, G., & Kneib, J.-P. 2002, Astronomy & Astrophysics, 390, 821
  * Hannuksela et al. (2019) Hannuksela, O., Haris, K., Ng, K., et al. 2019, The Astrophysical Journal Letters, 874, L2
  * Hannuksela et al. (2020) Hannuksela, O. A., Collett, T. E., Çalışkan, M., & Li, T. G. 2020, Monthly Notices of the Royal Astronomical Society, 498, 3395
  * Hansen et al. (2009) Hansen, S. M., Sheldon, E. S., Wechsler, R. H., & Koester, B. P. 2009, The Astrophysical Journal, 699, 1333
  * Hoekstra (2003) Hoekstra, H. 2003, How well can we determine cluster mass profiles from weak lensing?, The Royal Astronomical Society
  * Hou et al. (2020) Hou, S., Fan, X.-L., Liao, K., Zhu, Z.-H., et al. 2020, Physical Review D, 101, 064011
  * Janquart et al. (2021) Janquart, J., Hannuksela, O. A., Haris, K., & Van Den Broeck, C. 2021, Monthly Notices of the Royal Astronomical Society, 506, 5430
  * Janquart et al. (2023) Janquart, J., More, A., & Van Den Broeck, C. 2023, Monthly Notices of the Royal Astronomical Society, 519, 2046
  * Kafle et al. (2018) Kafle, P. R., Sharma, S., Lewis, G. F., Robotham, A. S., & Driver, S. P. 2018, Monthly Notices of the Royal Astronomical Society, 475, 4043
  * Kormann et al. (1994) Kormann, R., Schneider, P., & Bartelmann, M. 1994, Astronomy and Astrophysics, 284, 285
  * Li et al. (2023) Li, A. K., Lo, R. K., Sachdev, S., et al. 2023, Physical Review D, 107, 123014
  * Li et al. (2018) Li, S.-S., Mao, S., Zhao, Y., & Lu, Y. 2018, Monthly Notices of the Royal Astronomical Society, 476, 2220
  * Ma et al. (2023) Ma, H., Lu, Y., Guo, X., Zhang, S., & Chu, Q. 2023, Monthly Notices of the Royal Astronomical Society, 518, 6183
  * MacKay & JC (2003) MacKay, & JC, D. 2003, Information theory, inference and learning algorithms (Cambridge university press)
  * Mandelbaum et al. (2006) Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, Monthly Notices of the Royal Astronomical Society, 368, 715
  * Mishra et al. (2021) Mishra, A., Meena, A. K., More, A., Bose, S., & Bagla, J. S. 2021, Monthly Notices of the Royal Astronomical Society, 508, 4869
  * Mukherjee et al. (2021) Mukherjee, S., Broadhurst, T., Diego, J. M., Silk, J., & Smoot, G. F. 2021, Monthly Notices of the Royal Astronomical Society, 501, 2451
  * Nakamura (1998) Nakamura, T. T. 1998, Physical review letters, 80, 1138
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. 1997, The Astrophysical Journal, 490, 493
  * Ng et al. (2018) Ng, K. K., Wong, K. W., Broadhurst, T., & Li, T. G. 2018, Physical Review D, 97, 023012
  * Nightingale et al. (2019) Nightingale, J. W., Massey, R. J., Harvey, D. R., et al. 2019, Monthly Notices of the Royal Astronomical Society, 489, 2049
  * Oguri (2019) Oguri, M. 2019, Reports on Progress in Physics, 82, 126901
  * Pratten et al. (2021) Pratten, G., García-Quirós, C., Colleoni, M., et al. 2021, Physical Review D, 103, 104056
  * Robertson et al. (2020) Robertson, A., Smith, G. P., Massey, R., et al. 2020, Monthly Notices of the Royal Astronomical Society, 495, 3727
  * Romero-Shaw et al. (2020) Romero-Shaw, I. M., Talbot, C., Biscoveanu, S., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 3295
  * Sathyaprakash & Schutz (2009) Sathyaprakash, B. S., & Schutz, B. F. 2009, Living reviews in relativity, 12, 1
  * Schneider et al. (2006) Schneider, P., Kochanek, C., & Wambsganss, J. 2006, Gravitational lensing: strong, weak and micro: Saas-Fee advanced course 33, Vol. 33 (Springer Science & Business Media)
  * Seo et al. (2022) Seo, E., Hannuksela, O. A., & Li, T. G. 2022, The Astrophysical Journal, 932, 50
  * Shajib et al. (2021) Shajib, A. J., Treu, T., Birrer, S., & Sonnenfeld, A. 2021, Monthly Notices of the Royal Astronomical Society, 503, 2380
  * Smith et al. (2023) Smith, G. P., Robertson, A., Mahler, G., et al. 2023, Monthly Notices of the Royal Astronomical Society, 520, 702
  * Speagle (2020) Speagle, J. S. 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132
  * Takahashi & Nakamura (2003) Takahashi, R., & Nakamura, T. 2003, The Astrophysical Journal, 595, 1039
  * Tambalo et al. (2022) Tambalo, G., Zumalacárregui, M., Dai, L., & Cheung, M. H.-Y. 2022, arXiv preprint arXiv:2212.11960
  * Tessore & Metcalf (2015) Tessore, N., & Metcalf, R. B. 2015, Astronomy & Astrophysics, 580, A79
  * Treu (2010) Treu, T. 2010, Annual Review of Astronomy and Astrophysics, 48, 87
  * Wang et al. (2018) Wang, H., Mo, H., Chen, S., et al. 2018, The Astrophysical Journal, 852, 31
  * Wang et al. (1996) Wang, Y., Stebbins, A., & Turner, E. L. 1996, Physical review letters, 77, 2875
  * Wempe et al. (2022) Wempe, E., Koopmans, L. V., Wierda, A., Hannuksela, O. A., & Broeck, C. v. d. 2022, arXiv preprint arXiv:2204.08732
  * Wetzel et al. (2012) Wetzel, A. R., Tinker, J. L., & Conroy, C. 2012, Monthly Notices of the Royal Astronomical Society, 424, 232
  * Wierda et al. (2021) Wierda, A. R. A., Wempe, E., Hannuksela, O. A., Koopmans, L. V., & Van Den Broeck, C. 2021, The Astrophysical Journal, 921, 154
  * Woo et al. (2013) Woo, J., Dekel, A., Faber, S., et al. 2013, Monthly Notices of the Royal Astronomical Society, 428, 3306
  * Wright & Brainerd (2000) Wright, C. O., & Brainerd, T. G. 2000, The Astrophysical Journal, 534, 34
  * Yeung et al. (2023) Yeung, S. M., Cheung, M. H., Seo, E., et al. 2023, Monthly Notices of the Royal Astronomical Society, 526, 2230
  * Yu et al. (2020) Yu, H., Zhang, P., & Wang, F.-Y. 2020, Monthly Notices of the Royal Astronomical Society, 497, 204
  * Zahid et al. (2016) Zahid, H. J., Geller, M. J., Fabricant, D. G., & Hwang, H. S. 2016, The Astrophysical Journal, 832, 203
  * Zu & Mandelbaum (2016) Zu, Y., & Mandelbaum, R. 2016, Mapping stellar content to dark matter haloes–II. Halo mass is the main driver of galaxy quenching, The Royal Astronomical Society