11institutetext: Department of Astronomy, Beijing Normal University, 100875,
Beijing, China; 22institutetext: Advanced Institute of Natural Sciences,
Beijing Normal University at Zhuhai 519087, China; 33institutetext: Department
of Physics, College of Sciences, Northeastern University, Shenyang 110004,
China; 44institutetext: National Centre for Nuclear Research, Pasteura 7,
02-093 Warsaw, Poland; 55institutetext: School of Science, Chongqing
University of Posts and Telecommunications, Chongqing 400065, China

# Direct measurement of the distribution of dark matter with strongly lensed
gravitational waves

Shuo Cao 1122 Jingzhao Qi 33 Zhoujian Cao 1122 Marek Biesiada 44 Wei Cheng 55
Zong-Hong Zhu 1122<EMAIL_ADDRESS>
(Received XXX/ Accepted XXX)

In this Letter, we present a new idea of probing the distribution of dark
matter exhibiting elastic and velocity-independent self-interactions. These
interactions might be revealed in multiple measurements of strongly lensed
gravitational waves, which can be observationally explored to determine the
strength of self-scatterings. Specifically, each individual galactic-scale
strong-lensing system whose source is a coalescing compact binary emitting
gravitational waves will provide a model-independent measurement of the shear
viscosity of dark matter along the line of sight. These individual
measurements could be a probe of large-scale distribution of dark matter and
its properties. Our results indicate that with 10-1000 strongly lensed
gravitational waves from ET and DECIGO, robust constraints on the large-scale
distribution of self-interacting dark matter might be produced. More stringent
limits on the dark matter scattering cross-section per unit mass
($\sigma_{\chi}/m_{\chi}$) relevant to galaxy and cluster scales are also
expected, compared with the conservative estimates obtained in the
electromagnetic domain. Finally, we discuss the effectiveness of our method in
the context of self-interacting dark matter particle physics.

###### Key Words.:

Gravitational waves — Gravitational lensing — Dark matter

##  1 Introduction

It is well known that dark matter (DM) – the dominant component of virialized
objects (galaxies and clusters) – is one of the largest open questions in
modern astrophysics. Especially the collisionless cold DM (CCDM), which is
widely recognized by the astrophysics community and strongly supported by the
current observations of large-scale structures (e.g. galaxy clusters), follows
the homogeneous and isotropic distribution in the Universe (Schaye et al.,
2015; Springel et al., 2018). However, the CCDM paradigm still suffers from
the well-known small-scale problem that concerns some observed features of
galaxies and dwarf galaxies (e.g., missing satellite, core-cusp ,or too-big-
to-fail problems) (Read et al., 2017). Interestingly, there are suggestions
that self-interacting (SI) DM (SIDM) might successfully explain these
discrepancies. This provides an interesting alternative clue for the current
cosmological model (Tulin & Yu, 2018). Many ideas have been proposed to
explore the possibilities that DM SI might generate cosmic accelerated
expansion (Zimdahl et al., 2001; Cao et al., 2011) and non-zero cosmological
shear viscosity (Hawking, 1966). More interestingly, SIDM with velocity-
dependent SI has gathered growing attention because it can provide a
consistent fit on the scattering cross sections per unit mass that is
applicable to different scales (Hayashi et al., 2020). Despite these
advantages, the nature of SIDM remains unknown. It is in particular still
debated whether SIDM traces CCDM in the Universe because no direct detection
of its distribution has been performed.

According to General Relativity, absorption and dispersion of gravitational
waves (GWs) could be neglected in a perfect-fluid Universe (Ehlers et al.,
1996). This theoretical point of view has been widely applied in some recent
works (Qi et al., 2019). However, the hypothesis of transparent GWs in the EM
domain remains untested experimentally (Abbott et al., 2016). In this Letter,
we propose a new idea of measuring the properties of DM with elastic and
velocity-independent SI through the multiple measurements of galactic-scale
strong gravitational lensing systems with neutron star (NS) mergers acting as
background sources. The advantage of our method is that I) viscously damped
gravitational waves from these standard sirens would reach the observer,
accompanied by the electromagnetic radiation in the form of short and intense
bursts of $\gamma$ rays (Abbott et al., 2017). II) Each individual GW - galaxy
strong-lensing system will provide a model-independent measurement of the
fluid shear viscosity of dark matter along the line of sight. Therefore, with
a sample of viscosity measurements at different positions on the sky, the
distribution of dark matter in the Universe might be tested directly. III) The
strongly lensed NS-NS systems are ideal laboratories for deriving robust
constraints on the the scattering cross sections per unit mass relevant for
different scales, from ultra-dwarf galaxies to relaxed galaxy clusters.

##  2 Method

The LIGO and Virgo Collaborations have so far released the measurements of
luminosity distances for 90 GW sources (Abbott et al., 2021). However, we
stress that the luminosity distances $D_{L,obs}$ were derived from the
observed strain $h(t)$ and frequency evolution $f(t)$ under the assumption
that DM in the Universe is treated as a perfect fluid. On the other hand, if
the DM can be better represented as a non-ideal fluid characterized by a shear
viscosity term $\eta$, then the relation $\beta=\frac{16\pi G}{c^{3}}\eta$
between the GW damping rate ($\beta$) and DM viscosity ($\eta$) is simple.
This relation, which was originally proposed in Hawking (1966) with the
$c=8\pi G=1$ convention, was also derived in the framework of geometric units
(Ehlers et al., 1996; Atreya et al., 2018). In this analysis, we reintroduce
the fundamental constants in the formulae. Considering the inverse relation
between the amplitude of GW waveform and the luminosity distance, when the GW
damping rate is taken into account, the viscosity-free luminosity distance
($D_{L}$) inferred from the standard siren GW signal is modified to

$\displaystyle D_{L,eff}(z,\beta)=D_{L}(z)e^{\beta D(z)/2},$ (1)

where $D_{L,eff}$ and $D$ represent the so-called effective luminosity
distance and comoving distance, respectively (Goswami et al., 2017; Cao et
al., 2021). More recently, a method of measuring the viscosity of DM in the
cosmological context has been suggested and implemented for the currently 11
GW events that were released by LIGO and Virgo Collaborations (Lu et al.,
2018). However, a number of concerns were also raised about the robustness of
the above test: precise measurements of distances that are unaffected by
viscosity in the light of problems with host galaxy identification and
redshift determination for coalescing binary black holes (BH) are not
available. This motivates the search for other methods to probe the DM
viscosity. In this Letter we propose strongly lensed GWs produced by
coalescing NSs for this purpose.

Figure 1: Individual measurements of the damping rate of gravitational waves
based on future observations of strongly lensed NS-NS systems from ET (upper
panel) and DECIGO (lower panel).

As one of the successful predictions of General Relativity, strong
gravitational lensing by galaxies has become one of the most important tools
in studying cosmology (Cao et al., 2015) and the galaxy structure and
evolution (Cao et al., 2016). With the dawn of GW astronomy, the robust
prediction suggested that a considerable number of GW signals from
inspiralling NSs would be gravitationally lensed, focusing on the third
generation of ground-based GW detectors (Ding, Biesiada & Zhu, 2015) and on
the second generation of space-based GW detectors (Piórkowska et al., 2021),
respectively. In this Letter, our idea relies on the derived multiple
distances in a specific strong-lensing system (with the background GW source
at redshift $z_{s}$ and the lensing galaxy at redshift $z_{l}$), with strongly
lensed GW signals observed together with their electromagnetic (EM)
counterparts. We outline below three approaches to measure the angular
diameter distances between the observer and the lens ($D_{A}(z_{l})$), between
the observer and the source ($D_{A}(z_{s})$), and between the lens and the
source ($D_{A}(z_{l},z_{s})$).

Firstly, with the measurements of the central velocity dispersion and the
location of multiple images (Einstein radius) of a strong-lensing system, the
ratio of two angular diameter distances $D_{A}(z_{l},z_{s})/D_{A}(z_{s})$ can
be precisely assessed. Specifically, by assuming a power-law model to describe
the lens mass distribution ($\rho\sim r^{-\gamma}$), the lens potential
relates to the moments of the stellar distribution function through the Jeans
equation (Koopmans et al., 2005). The combination of the lens mass and the
dynamical mass inside the Einstein radius ($\theta_{E}$) yields (Cao et al.,
2015)

$\frac{D_{A}(z_{l},z_{s})}{D_{A}(z_{s})}=\frac{\theta_{E}}{4\pi}\frac{c^{2}}{\sigma_{ap}^{2}}\left(\frac{\theta_{E}}{\theta_{ap}}\right)^{\gamma-2}f(\gamma)^{-1},$
(2)

where the luminosity-averaged velocity dispersion $\sigma_{ap}$ is measured
within the aperture $\theta_{ap}$ (see Cao et al. (2015) for the expression of
$f(\gamma)$ in the form of the radial mass profile slope). Therefore, the
observations of $\sigma_{ap}$, $\theta_{E}$, $\theta_{ap}$ , and $\gamma$
enables precise measurements of the distance ratio
$D_{A}(z_{l},z_{s})/D_{A}(z_{s})$. Secondly, the time-delay distance in a
specific GW - galaxy strong-lensing system can be measured. The time-delay
distance is defined as

$D_{\mathrm{\Delta
t}}(z_{l},z_{s})\equiv\frac{D_{A}(z_{l})D_{A}(z_{s})}{D_{A}(z_{l},z_{s})}=\frac{c}{1+z_{l}}\frac{\Delta
t_{i,j}}{\Delta\phi_{i,j}}.$ (3)

It could be measured precisely and accurately for transient sources such as GW
signals from coalescing binaries through the precise measurement of the time
delay and the well-reconstructed Fermat potential difference
($\Delta\phi_{i,j}$) between multiple images, that is, the multiple GW
signals. The EM counterpart would allow the host galaxy and lens galaxy
identification, and dedicated lens modelling techniques would enable a precise
reconstruction of the lens mass distribution. More importantly, the time
difference ($\Delta t_{i,j}$) in the arrival times of two signals (at angular
coordinates $\boldsymbol{\theta}_{i}$ and $\boldsymbol{\theta}_{j}$ on the
sky) can also be measured with unprecedented accuracy (Liao et al., 2017). We
stress that lensing time delays are not affected by GW damping effect due to
DM viscosity. This advantage of strongly lensed GW signals has been widely
discussed in fundamental physics (Collett & Bacon, 2017) and cosmology (Cao et
al., 2019). Now, the angular diameter distance to the lens can be calculated
as

$D_{A}(z_{l})=D_{\mathrm{\Delta t}}\frac{D_{A}(z_{l},z_{s})}{D_{A}(z_{s})}.$
(4)

We note that based on the Etherington reciprocity relation, which is of
fundamental importance in modern cosmology (Etherington, 1933), combined with
distance sum rule, the luminosity distance to the source can be expressed as

$\displaystyle D_{L}(z_{s})=$
$\displaystyle=(1+z_{s})(1+z_{l})D_{\mathrm{\Delta
t}}\frac{D_{A}(z_{l},z_{s})/D_{A}(z_{s})}{1-\left(D_{A}(z_{l},z_{s})/D_{A}(z_{s})\right)(1+z_{s})^{-2,}}$

which is unaffected by GW damping effect and DM SIs. Thirdly, another
advantage of strongly lensed GW events comes from their ability to provide
simultaneous measurements of the effective luminosity distance to the source,
that is, coalescing binary NSs can be calibrated as standard sirens. However,
due to the magnification effect of lensing, the luminosity distance measured
from the strongly lensed GW signals is not the true $D_{L,eff}$, but rather
$D_{L,eff}/A$, where $A$ is the lensing amplification factor. If the lensing
magnification $\mu$ could be independently measured through photometric EM
observations, we would be able to correct the amplitude of the GW signal for
the amplification factor $A=\sqrt{\mu}$. _Therefore, the true effective
luminosity distance from the observer to the source, which is affected by
viscous DM damping, could be derived as_

$D_{L,eff}(z_{s},\beta)=D_{L}(z_{s})e^{\beta
D_{L}(z_{s})/2(1+z_{s})^{2}}/\sqrt{\mu},$ (5)

_based on a detailed analysis of the images of the host galaxy and the
accompanying EM counterpart._

The function $\beta(z_{l},z_{s})$ can now be directly obtained through the
combined analysis with Eqs. (2)-(5) for individual lenses. This general
function is applicable to any form of SIDM along the line of sight. In the
following section, we apply this method to the simulated data of future
ground-based and space-based GW detectors.

##  3 Simulated data and constraints

A particular single strong-lensing system possesses its own sensitivity to DM
viscosity. The recent analysis of Cutler & Holz (2009) revealed that the
DECIGO, a proposed space-based GW detector, can discover up to 1000 lensed GW
events in the accumulated data of five to ten years. This is an improvement of
more than one order of magnitude over the three nested ET interferometers in
the redshift range $z\leq 5.00$ (Li et al., 2018). Following the approach
proposed by Collett (2015), we have simulated a realistic population of
strongly lensed GWs (10 and 1000 events in ET and in DECIGO, respectively),
the sampling distribution of which follows the intrinsic merger rate of double
compact objects (DCO) calibrated by strong-lensing effects (Ding, Biesiada &
Zhu, 2015). The lensed GW systems are randomly sampled on the sky, with the
masses of the NSs uniformly sampled in the interval of
$[1M_{\odot},2M_{\odot}]$. The number density of lenses, assumed to have a
Gaussian distribution of the central velocity dispersion similar to the SL2S
sample ($\sigma_{0}=210\pm 50$ km/s) (Sonnenfeld et al., 2013), is
characterized by the velocity dispersion function (VDF) of early-type galaxies
from SDSS DR5 (Choi et al., 2007). Based on the uncertainties of multiple
measurements in strongly lensed gravitational waves (see the appendix for
details), we numerically propagated the uncertainties of different observables
(time delays, mass density profile, Einstein radius, velocity dispersion,
line-of-sight contamination, etc.) to the uncertainty of different distances,
and finally to the uncertainties of the DM viscosity on which it depends.

Following the method outlined above, we show the measurements of
$\beta(z_{l},z_{s})$ expected to be detected by future ET and DECIGO in Fig.
1. Our simulation is based on the concordance $\Lambda$CDM model
(collisionless CDM) with the parameters determined by Planck 2018 data
(Aghanim et al., 2020).

_The question now is whether these measurements are sufficient to trace the
distribution of DM and its properties._ We find that both ET and DECIGO would
provide precise 3D measurements of the GW damping rate $\beta(z)$ and DM
viscosity $\eta(z)$ at different redshifts and different positions on the sky.
The significant discrepancy between these individual measurements could be a
probe of the large-scale distribution of DM and its properties. In particular,
the second-generation space-based detectors (e.g. DECIGO) would provide a
catalog of more precise $\beta(z)$ measurements that would be larger by one
order of magnitude. This would increase the possibility of determining
different $\beta(z_{l},z_{s})$ for different pairs $(z_{l},z_{s})$.
Summarizing the multiple $\beta$ measurements through the inverse-variance
weighting, a constraint on the $\beta$ parameter at the precision of
$\Delta\beta=10^{-4}\,{\rm Mpc}^{-1}$ might be expected, with ten multi-
messenger strongly lensed GW events from ET. For comparison, we also
considered the constraining power of the same technique in space-borne
detectors, using simulated data representative of the future observations of
lensed GW and EM signals. The result is shown in Fig. 3. Lensed systems seen
jointly in GW and EM signals by DECIGO give much more stringent constraints,
the uncertainty of the damping of GWs in a viscous Universe being
$\Delta\beta=10^{-6}\,{\rm Mpc}^{-1}$. This can be understood because the
number of strongly lensed standard sirens as well as of precise distance
measurements at much higher redshifts is higher. We note that statistically
significant differences between particular strong-lensing systems can be used
to distinguish between DM and modified gravity. Namely, as discussed for
instance in some classes of modified gravity theories (Belgacem et al., 2018a,
b; Mastrogiovanni, Steer & Barsuglia, 2020), GW propagation involves a
friction term. In consequence, GW amplitude is attenuated in a similar manner
as by DM shear viscosity. However, this effect should be isotropic, while
differences in DM distribution along the line of sight could manifest
themselves as discussed in this section.

Figure 2: Fits on the damping rate of GWs based on strongly lensed NS-NS
systems from ET and DECIGO. Figure 3: Limits on the DM scattering cross-
section per unit mass ($\sigma_{\chi}/m_{\chi}$) on the scales of galaxy and
galaxy clusters, based on the current and future observations in the EM and GW
domain.

Now the main question now is whether this accuracy is enough to distinguish
between possible DM shear viscosity and SI cross-sections necessary to solve
CCDM problems at galaxy and cluster scales. SIDM is parametrized by the
scattering cross-section per unit mass ($\sigma_{\chi}/m_{\chi}$), which is
generally a function of the relative velocity of the DM particles ($\chi$)
(Tulin & Yu, 2018; Atreya et al., 2018),

$\displaystyle\frac{\left\langle\sigma_{\chi}\right\rangle}{m_{\chi}}=\frac{6.3\pi
G\left\langle v\right\rangle}{c^{3}\beta}.$ (6)

In particular, because the collision rate becomes negligible in low-density
regions, SIDM particles behave like collisionless CDM at larger length scales.
These newly proposed velocity-dependent SIs provide a consistent fit to the
CDM paradigm at different scales, provided $\sigma_{\chi}/m_{\chi}\sim 1$
cm2/g on the scale of galaxies (Dooley et al., 2005) and
$\sigma_{\chi}/m_{\chi}\sim 0.1$ cm2/g on the scale of galaxy clusters
(Markevitch et al., 2004). However, the class of velocity-dependent SIDM
models remains largely unconstrained. Similarly, in this paper, two SIDM
scenarios are also considered with our method, with the phenomenological SIDM
halo model proposed in Kaplinghat et al. (2016). More specifically, focusing
on galactic-scale DM halos, ten lensed GWs detected by ET will provide a
stringent limit on its SI cross section as

$\Delta(\sigma_{\chi}/m_{\chi})\sim 10^{-4}cm^{2}/g.$ (7)

A much stronger limit on the DM SI cross section,

$\Delta(\sigma_{\chi}/m_{\chi})\sim 10^{-6}cm^{2}/g,$ (8)

can be obtained with 1000 lensed GW events detected by DECIGO. Previously, by
analyzing the rotation curves of five dwarf galaxies in the THINGS sample (Dh
et al., 2011) and seven LSB galaxies from Kuzio et al. (2008), Kaplinghat et
al. (2016) derived a limit on the cross section
$\Delta(\sigma_{\chi}/m_{\chi})\sim 5\times 10^{-1}$ cm2/g for all 12
galaxies. In the present paper, using 10-1000 strongly lensed GWs from the
third-generation ground-based GW detectors and space-based detectors, our
constraints are $10^{3}\sim 10^{5}$ times tighter than the results obtained
before. For the cluster-scale case, the limits we derive here are
$\Delta(\sigma_{\chi}/m_{\chi})\sim 10^{-5}$ cm2/g and
$\Delta(\sigma_{\chi}/m_{\chi})\sim 10^{-7}$ cm2/g in the framework of ET and
DECIGO, respectively. Therefore, we could obtain very optimistic limits on the
cluster-scale DM SI cross section, which represents an improvement of three to
five orders of magnitude over the results using six clusters of galaxies
($\Delta(\sigma_{\chi}/m_{\chi})\sim 2\times 10^{-2}$ cm2/g) (Kaplinghat et
al., 2016). The results are shown in Fig. 3. We now briefly describe the
effectiveness of our method in the context of SIDM particle physics. To
elaborate this idea, we take the simplest SIDM model, which could be a real
scalar field $s$ with quartic SIs ($\lambda_{s}s^{4}/4!$) (de Martino et al.,
2020). Although there are several ways to build this type of models, the self-
coupling strength in the dark sector is characterized by $\lambda_{s}$, which
should be in the range of ($0<\lambda_{s}<\sqrt{\pi/9}$) due to the stability
of the potential and perturbation limits (Cheng & Bian, 2019). Fig. 4 shows
the preferred region of the mass of the DM particle $m_{s}$ and the self-
coupling strength $\lambda_{s}$ on the scales of galaxies and galaxy clusters,
based on the simulated data regarding strongly lensed GW signals from the ET.
Our final results demonstrate that in the framework of SIDM capable of solving
the small-scale issues of the CDM paradigm, the DM particle mass (within some
specific SIDM models) can be assessed with an accuracy of $\sim
O(0.1)\rm{GeV}$. This conclusion, which could be extended to space-based
detectors (DECIGO) and other SIDM models (e.g. fermionic DM particles that
self-interact via a light mediator) (Tulin et al., 2013a, b), will strengthen
the constraining power of our method to inspire new observational programs and
theoretical work in the future.

Figure 4: Limits on the mass of DM particles and the coupling constant
reflecting the strength of SIs on the scales of galaxies and galaxy clusters,
based on lensed GW signals in ET. The dotted lines denotes the central value,
and the 68.3% C.L. is shown by the red (blue) shaded region.

##  4 Summary and outlook

The impact of the results presented in this work is expected to increase
significantly as more strongly lensed GWs are observed. However, based on our
results and arguments, the final question still is whether it is technically
possible to achieve uncertainties this small for the measurements of DM self-
scattering. Many technical problems would be involved from the observational
point of view. For instance, the current low detection rate of inspiraling
double compact object detection and lensing rate may be due to the lack of
either the necessary high temporal resolution or a wide enough field of view
in the current EM and GW detectors. Fortunately, the upcoming and planned GW
detectors (e.g. ET and DECIGO) and wide-area and deep surveys (e.g. the LSST),
with a much wider field of view and higher sensitivity, will be able to
discover and precisely localize a large number of strongly lensed GW and EM
events. However, it is very likely that only some fraction of the GWs with EM
counterparts will be measured. This means that we would observe some lensed GW
signals from NS-NS systems without an electromagnetic counterpart. Moreover,
accurate lens reconstruction could also be difficult when the Einstein radius
of the lens is too small or the background galaxy is too dim for a redshift
measurement. Therefore, we need to be cautious about the number of strongly
lensed GW signals that allow both the determination of EM counterparts and
precise lens reconstruction. Secondly, the systematic floor needs to be
further studied for future works in order to let us know how powerful this
lensing method might be in probing the distribution of DM exhibiting elastic
and velocity-independent SIs. For instance, we still need precise measurements
of the image magnification (or the amplification factor of GW signals), which
calls for more knowledge about the AGN accretion models and local image
environments from astrophysics inputs. Finally, if the estimates of the mean
collision velocity for DM particles have large uncertainties, the results
shown in Fig. 3-4 might be affected. Accounting for the uncertainty of the
mean collision velocity on the scales of galaxy and cluster, Eq. (6) might
make these limits less stringent. To summarize, with more abundant
observational information in the future, we will have a better understanding
of the SI nature of DM at both galaxy and cluster scales, and we will also
know better how to use this as a probe of the large-scale distribution of DM
in the GW domain, as discussed in this Letter.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
under Grants Nos. 12021003, 11690023, and 11920101003; the National Key R&D
Program of China No. 2017YFA0402600; the Strategic Priority Research Program
of the Chinese Academy of Sciences, Grant No. XDB23000000; the Interdiscipline
Research Funds of Beijing Normal University; the China Manned Space Project
(Nos. CMS-CSST-2021-B01 and CMS-CSST-2021-A01); and the CAS Project for Young
Scientists in Basic Research under Grant No. YSBR-006. M.B. was supported by
Foreign Talent Introducing Project and Special Fund Support of Foreign
Knowledge Introducing Project in China (No. G2021111001L).

## References

  * Abbott et al. (2016) Abbott, B. P., et al. 2016, PRL, 116, 061102
  * Abbott et al. (2017) Abbott, B. P., et al. 2017, PRL, 119, 161101
  * Abbott et al. (2021) Abbott, B. P., et al. 2021, arXiv:2111.03606
  * Aghanim et al. (2020) Aghanim, N., Akrami, Y., Ashdown, M., et al. 2020, A&A, 641, A6 [arXiv:1807.06209]
  * Atreya et al. (2018) Atreya, A., Bhatt, J.R. & Mishra, A. 2018, JCAP, 02, 024
  * Auger et al. (2010) Auger, M. W., et al. 2010, ApJ, 724, 511
  * Belgacem et al. (2018a) Belgacem, E., Dirian, Y., Foffa, D., & Maggiore, M. 2018a, PRD, 97, 104066
  * Belgacem et al. (2018b) Belgacem, E., Dirian, Y., Foffa, D., & Maggiore, M. 2018b, PRD, 98, 023510
  * Bonvin et al. (2018) Bonvin, V., et al. 2018, A&A, 616, A183
  * Cai et al. (2016) Cai, R.-G., et al. 2016, PRD, 93, 043517
  * Cai & Yang (2017) Cai, R.-G., & Yang, T. 2017, PRD, 95, 044024
  * Cao et al. (2011) Cao, S., Zhu, Z.-H. & Liang, N. 2011, A&A, 529, A61
  * Cao et al. (2015) Cao, S., et al. 2015, ApJ, 806, 185
  * Cao et al. (2016) Cao, S., et al. 2016, MNRAS, 457, 281
  * Cao et al. (2019) Cao, S., et al. 2019, NatSR, 9, 11608
  * Cao et al. (2021) Cao, S. et al. 2021, MNRAS Letters, 502, L16
  * Cao et al. (2022) Cao, S. et al. 2022, ApJ, in press [arXiv:2112.00237]
  * Chen et al. (2018) Chen, G. C.-F., et al. 2018, MNRAS, 481, 1115
  * Cheng & Bian (2019) Cheng, W. & Bian, L. 2019, PRD, 99, 035038
  * Choi et al. (2007) Choi, Y.-Y., et al. 2007, ApJ, 658, 884
  * Collett (2015) Collett, T. E. 2015, ApJ, 811, 20
  * Collett & Cunnington (2016) Collett, T. E. & Cunnington, S. D. 2016, MNRAS, 462, 3255
  * Collett & Bacon (2017) Collett, T. E., Bacon, D. 2017, PRL, 118, 091101
  * Cutler & Holz (2009) Cutler, C., & Holz, D. E. 2009, PRD, 80, 104009
  * de Martino et al. (2020) de Martino, I., et al. 2020, Universe, 6, 107
  * Dexter & Agol (2011) Dexter, J. & Agol, E. 2011, ApJL, 727, L24
  * Dh et al. (2011) Oh, S.-H., et al. 2011, ApJ, 141, 193
  * Ding, Biesiada & Zhu (2015) Ding, X., Biesiada, M. & Zhu, Z-H., 2015, JCAP, 12, 006
  * Dooley et al. (2005) Dooley, G. A., et al. 2016, MNRAS, 461, 710
  * Ehlers et al. (1996) Ehlers, J. & Prasanna, A. R., 1996, CQG, 13, 2231
  * Etherington (1933) Etherington, I. M. H. 1933, Philosophical Magazine, 15, 761 (reprinted in 2007, Gen. Rel. Grav. 39, 1055)
  * Foxley-Marrable et al. (2018) Foxley-Marrable, M., et al. 2018, MNRAS, 478, 5081
  * Geng et al. (2020) Geng, S., et al. 2020, ApJ, 905, 54
  * Goobar et al. (2017) Goobar, A., Amanullah, R., Kulkarni, S. R., et al. 2017, Science, 356, 291
  * Goswami et al. (2017) Goswami, G., et al. 2017, PRD, 95, 103509
  * Hawking (1966) Hawking, S. W., 1966, ApJ, 145, 544
  * Hayashi et al. (2020) Hayashi, K., et al. 2021, PRD, 103, 023017
  * Hlozek et al. (2019) Hlozek, R. A., et al. arXiv:1903.09324
  * Kaplinghat et al. (2016) Kaplinghat, M., et al. 2016, PRL, 116, 041302
  * Kawamura (2006) Kawamura, S., et al. 2006, CQG, 23, 125
  * Koopmans et al. (2005) Koopmans L.V.E. 2005, Proceedings of XXIst IAP Colloquium, “Mass Profiles & Shapes of Cosmological Structures” (Paris, 4-9 July 2005), eds G. A. Mamon, F. Combes, C. Deffayet, B. Fort (Paris: EDP Sciences) [astro-ph/0511121]
  * Kuzio et al. (2008) Kuzio, R., et al. 2008, ApJ, 676, 920
  * Li et al. (2018) Li, S. S., Mao, S. & Lu, Y., 2018, MNRAS, 476, 2220
  * Li (2015) Li, T. G. Extracting Physics from Gravitational Waves: Testing the Strong-field Dynamics of General Relativity and Inferring the Large-scale Structure of the Universe (Springer, 2015)
  * Liao et al. (2017) Liao, K., et al. 2017, Nature Communications, 8, 1148
  * Liao et al. (2019) Liao, K. 2019, ApJ, 871, 113
  * Lu et al. (2018) Lu, B. Q., et al. arXiv:1803.11397
  * Markevitch et al. (2004) Markevitch, M., et al. 2004, ApJ, 606, 819
  * Mastrogiovanni, Steer & Barsuglia (2020) Mastrogiovanni, S., Steer, D. & Barsuglia, M. 2020, PRD, 102, 044009
  * Nissanke et al. (2010) Nissanke et al. 2010, ApJ, 725, 496
  * Oguri (2010) Oguri, M. 2010, PASJ, 62, 1017
  * Oguri & Marshall (2010) Oguri, M. & Marshall, P. J. 2010, MNRAS, 405, 2579
  * Pan et al. (2021) Pan, Y., et al. 2021, ApJ, 911, 135
  * Piórkowska et al. (2021) Piórkowska, A., et al. 2021, ApJ, 908, 196
  * Qi et al. (2019) Qi, J. Z., et al. 2019, PRD, 99, 063507
  * Read et al. (2017) Read, J., et al. 2017, MNRAS, 467, 2019
  * Ruff et al. (2011) Ruff, A., et al. 2011, ApJ, 727, 96
  * Schaye et al. (2015) Schaye, J., et al. 2015, MNRAS, 446, 521
  * Seljak (1994) Seljak, U. 1994, ApJ, 436, 509
  * Sonnenfeld et al. (2012) Sonnenfeld, A., et al. 2012, ApJ, 752, 163
  * Sonnenfeld et al. (2013) Sonnenfeld, A., Gavazzi, R., Suyu, S.H., Treu, T., & Marshall, P. J. 2013, ApJ, 777, 97
  * Springel et al. (2018) Springel, V., et al. 2018, MNRAS, 475, 676
  * Suyu et al. (2010) Suyu, S. H., et al. 2010, ApJ, 711, 201
  * Suyu et al. (2012) Suyu, S. H., et al. 2012, ApJ, 750, 10
  * Suyu et al. (2017) Suyu, S. H., et al. 2017, MNRAS, 468, 2590
  * Tulin et al. (2013a) Tulin, S., Yu, H.-B., & Zurek, K. M. 2013, PRD, 87, 115007
  * Tulin et al. (2013b) Tulin, S., Yu, H.-B., & Zurek, K. M. 2013, PRL, 110, 111301
  * Tulin & Yu (2018) Tulin, S. & Yu, H. B. 2018, Phys. Rept. 730, 1
  * Yahalomi et al. (2017) Yahalomi, D., et al. 2017, arXiv:1711.07919
  * Zhao et al. (2011) Zhao, W., et al. 2011, PRD, 83, 023005
  * Zheng et al. (2021) Zheng, X., et al. 2021, EPJC, 81, 14
  * Zimdahl et al. (2001) Zimdahl, W., et al. 2001, PRD, 64, 3501

## Appendix

We summarize the uncertainties of multiple measurements in strongly lensed NS-
NS systems, that is, the precision of lens reconstruction, the precision of
time delays, and the precision of lensed standard sirens. The error strategy
below is used to produce the measurements of $\beta(z_{l},z_{s})$ expected to
be derived by future GW detectors (see Fig. 1 for details).

### Precision of lens reconstruction

Three sources of uncertainty are included in our simulation of lens
reconstruction. With the observation of host arcs by deep imaging of the HST,
the Einstein radius of multiple images and the lens mass slope can be measured
precisely and accurately (at the level of 1% ) based on advanced lens
modelling techniques (Suyu et al., 2010, 2012) and kinematic modelling methods
(Auger et al., 2010; Sonnenfeld et al., 2012). This error strategy has been
extensively used in the simulation of the LSST lens sample, with high-quality
(sub-arcsecond) imaging data in general (Collett & Cunnington, 2016). Although
the joint lensing and stellar-dynamical analysis could generate 5% accuracy on
the total mass-density slope (Ruff et al., 2011), the lessons learnen from
lensed quasar results showed that this accuracy is expected to increase to 1%,
with supplementary information about time delays obtained by the COSMOGRAIL
(Bonvin et al., 2018) or H0LiCOW collaboration (Suyu et al., 2017). Meanwhile,
we take 5% as the fractional uncertainty of the observed velocity dispersion,
with detailed follow-up spectroscopic information from other ground-based
facilities, that is, multi-object and single-object spectroscopy to enhance
dark energy science from LSST (Hlozek et al., 2019).

### Precision of time delays

We consider three sources of uncertainty in the simulation of time-delay
measurements. It is well recognized that coalescing NS-NS systems manifested
as short gamma-ray bursts can have many advantages as precise lensing time-
delay indicators. However, lensed GW signals could provide extremely accurate
$\Delta t$ measurements, with a detailed analysis of the waveform from the GW
detection pipeline (Liao et al., 2017) (accuracy ranging from $10^{-4}\;ms$ to
$0.1\;s$). Therefore, the fractional uncertainty of $\Delta t$ could be
neglected. On the other hand, the detection of the EM counterpart would be
very beneficial. The host galaxy image followed up with a quality typical to
HST observations with dedicated lens modelling techniques would result in a
precision of 0.5% in lens modelling; specifically, the Fermat potential
difference reconstruction for a well-measured lensed GW system (Liao et al.,
2017). The analogous uncertainty is $3\%$ for almost all quasar-galaxy
systems, with lensed quasar images quality typical to HST observations (Suyu
et al., 2017). In this analysis, we applied a strategy to calculate the
uncertainty of $\Delta\phi$, from the measurements of the lens mass profile,
the Einstein radius, and their corresponding uncertainties (Cao et al., 2019).
Finally, the effect of the light-of-sight density fluctuation on the
gravitational time delays should be taken into account. For each lens system,
we added an additional 1% relative uncertainty to the lens modelling, which
naturally appears due to mass along the line of sight to the source (Seljak,
1994). This additional line-of-sight systematic uncertainty (at the 1% level)
is also characteristic for that derived in typical strongly lensed quasar
systems (Suyu et al., 2017).

### Precision of lensed standard sirens

We account for two sources of uncertainties in the simulation of lensed GWs:
effective luminosity distance, and lensing magnification. Firstly, to obtain
the uncertainty of $D_{L,eff}(z)$ from the GW signal of NS-NS standard sirens
observed by different detectors, the distance precision was taken as
$\sigma^{2}_{\rm tot}=\sigma^{2}_{\rm ins}+\sigma^{2}_{\rm len}$, with the
instrumental measurement uncertainty $\sigma_{inst}\simeq
2D^{L}_{s}/\rho_{net}$, where $\rho_{net}$ denotes the combined signal-to-
noise ratio (S/N) for the network of independent interferometers, that is, the
square root of the inner product of the strain in Fourier space (Cai et al.,
2016). A detailed description of the noise power spectral density (PSD), which
characterizes the sensitivity of the detector, could be found in Cai et al.
(2016) and Kawamura (2006). Note that we assumed an ideal situation by setting
the inclination angle of the binary orbital angular momentum $\iota=0$.
However, Nissanke et al. (2010) demonstrated the strong correlation between
the luminosity distance and the inclination angle of the binary
($\iota=[0,180^{\circ}]$). This distance-inclination degeneracy, which
decreases the sensitivity to both parameters, might be broken by the strongly
beamed observations ($\iota<20^{\circ}$) of short gamma-ray bursts as the
electromagnetic counterparts of GWs (Li, 2015). Therefore, considering the
maximum effect of the inclination on the combined S/N, we chose to double the
instrumental uncertainty of $D_{L}$ as the upper limit, as was proposed in the
recent study of the third-generation GW detector (ET) based on a Fisher
information matrix (Cai & Yang, 2017). This strategy has been extensively
applied in subsequent works focusing on the constraint ability of standard
siren GWs on different cosmological parameters (Geng et al., 2020; Pan et al.,
2021; Zheng et al., 2021; Cao et al., 2022). Moreover, the measurements of
luminosity distance are also affected by the weak-lensing effect, which is
modeled as $\sigma_{lens}/D^{L}_{s}=0.05z$ for ET (Zhao et al., 2011) and
$\sigma_{lens}/D^{L}_{s}=0.044z$ for DECIGO (Cutler & Holz, 2009). Secondly,
the amplification factor of the GW signal along with the corresponding
uncertainty is related to the magnification ($\mu$) and its uncertainty of the
EM counterpart, following the recent discussion of gravitational lensing
statistics and the magnification by galaxies (Suyu et al., 2010, 2012). More
specifically, the uncertainty of $\mu$ is calculated from the observables and
their uncertainties (i.e. the Einstein radius and the mass profile slope) (Cao
et al., 2019) by solving the lens equation with _glafic_ in the framework of a
specific lens model (Oguri, 2010; Oguri & Marshall, 2010). In our simulated
data, the image magnification is expected to be determined at a level of $\sim
10\%$. Meanwhile, recent works have studied the microlensing (ML) effect
generated by stars in lensing galaxies, which potentially biases the
magnification map of strongly lensed sources (Foxley-Marrable et al., 2018).
Although the community is mainly of the opinion that this effect needs to be
incorporated in a strong-lensing analysis, much uncertain input for the ML
priors on the AGN accretion disk models remains (Dexter & Agol, 2011), as well
as on the local environments for images and mass function for the stars (Chen
et al., 2018). Therefore, to assess the impact of ML on luminosity distance
estimates, we assumed an additional systematic 5% uncertainty on the
magnification measurements for all lensed gravitational waves. A similar error
strategy can be found in the recent discussion of Liao et al. (2019), focusing
on ML time delays from the light curves of strongly lensed quasars. We
emphasize that for galactic-scale lenses with point-like background sources,
the model predictions in the generic case for image magnifications could be
very large (Yahalomi et al., 2017). However, for strongly lensed standard
sirens in the GW domain, the lensing magnification could be more precisely
measured in a lens-model-independent way by comparing the observed flux with
that of other unlensed GWs within a narrow redshift bin (Goobar et al., 2017).