# Lensing Magnification Seen by Gravitational Wave Detectors

Giulia Cusin, Ruth Durrer, Irina Dvorkin

► To cite this version:

Giulia Cusin, Ruth Durrer, Irina Dvorkin. Lensing Magnification Seen by Gravitational Wave Detectors. Universe, 2021, 8 (1), pp.19. 10.3390/universe8010019. hal-03559991

HAL Id: hal-03559991

https://hal.science/hal-03559991

Submitted on 2 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Lensing Magnification Seen by Gravitational Wave Detectors

Giulia Cusin<fnref n="1" />,*, Ruth Durrer<fnref n="1" /> and Irina Dvorkin<fnref n="2" />

**Abstract:** In this paper, we studied the gravitational lensing of gravitational wave events. The probability that an observed gravitational wave source has been (de-)amplified by a given amount is a detector-dependent quantity which depends on different ingredients: the lens distribution, the underlying distribution of sources and the detector sensitivity. The main objective of the present work was to introduce a semi-analytic approach to study the distribution of the magnification of a given source population observed with a given detector. The advantage of this approach is that each ingredient can be individually varied and tested. We computed the expected magnification as both a function of redshift and of the *observed source luminosity distance*, which is the only quantity one can access via observation in the absence of an electromagnetic counterpart. As a case study, we then focus on the LIGO/Virgo network and on strong lensing ($\mu > 1$).

**Keywords:** gravitational waves; gravitational lensing; detector sensitivity; magnification

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Cusin, G.; Durrer, R.; Dvorkin, I. Lensing Magnification Seen by Gravitational Wave Detectors. *Universe* **2022**, *8*, 19. https://doi.org/10.3390/universe8010019

Academic Editor: B. P. Bonga
(Béatrice)

Received: 26 November 2021

Accepted: 20 December 2021

Published: 30 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The direct detection of gravitational waves (GWs) from inspiraling binary black holes and neutron stars by the LIGO and Virgo experiments [1–11] has been one of the greatest successes in experimental physics and observational astronomy in recent years. Not only did it lead to unprecedented insights into black holes and especially neutron stars, but it also allowed us to test general relativity (GR) and cosmology.

An interesting aspect is the fact that the clustered matter between the gravitational wave source and the observer can amplify the observed signal. This may lead to an underestimation of the distance to the source and therefore of its redshift which in turn leads to an overestimation of the chirp mass [12]. The effects of gravitational lensing on the apparent properties of detected sources and the possible identification of unique lensing features were extensively studied [13–28]. The prospects of the lensing of gravitational waves alone or in combination with other data as a probe of dark matter or of general relativity have also been investigated [29–32].

In this paper, we sought to determine the probability that an event of a given strain amplitude, or rather of a given signal-to-noise ratio (SNR) for fixed experimental characteristics, is being significantly magnified by lensing. In particular, we studied the dependence of the lensing probability and of the mean magnification on redshift and on the strain amplitude (or SNR). Similar studies were already present in the literature [12,15,18,25,26,33–35]. The main difference with the present work is that, instead of full numerical simulations with ray tracing, we presented semi-analytic formulae which can be applied to an arbitrary lens and source distribution, maintaining full control of the modeling and transparency of all physical effects. This allowed us to study the dependence of the result on the lens and/or source distribution which is otherwise quite obscure. Our analysis can also be easily adapted to an arbitrary experiment given the expected number of sources as a function of their SNR and redshift. Moreover, unlike previous studies on this subject, we compute the expected magnification as a function of redshift, but also as a function of the *observed*

<footnote n="1">¹ Centre for Astroparticle Physics, Département de Physique Théorique, Université de Genève, 24 quai Ernest-Ansermet, 1211 Geneva, Switzerland; giulia.cusin@unige.ch</footnote>

<footnote n="2">² Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis Boulevard Arago, 75014 Paris, France; dvorkin@iap.fr</footnote>

<footnote>* Correspondence: ruth.durrer@unige.ch</footnote>
luminosity distance, which is the only quantity one can access observationally in the absence of an electromagnetic counterpart. While in principle it is possible to determine the redshift of an inspiraling binary black hole from the light curve, the LIGO/Virgo observations have by far too large error bars to do this. Here, we assume the Planck concordance cosmology to determine the redshift from the luminosity distance. If lensing is relevant, this redshift is, however, not correct, as the inferred distance is affected by lensing. The advantage of having a semi-analytic analysis is that it allows us to answer the question of whether it is feasible that some or most of the observed LIGO/Virgo events are strongly lensed, as proposed in [16,17,27,36]. We show that, in order to achieve a significant lensing probability for more than one of these events, one has to tweak the lens and source distribution in a way that is in contradiction with the present observational constraints and with results from simulations.

In this work, we focus on strong lensing. We are indeed interested in working out the distribution of magnification $\mu > 1$ for sources visible by LIGO/Virgo/Kagra. The magnification induced by a cosmological large scale-structure, which can also lead to de-magnification, $\mu < 1$, is negligible for this study. For the distances of the LIGO/Virgo sources, it only leads to small corrections to the results presented herein. We modeled the lenses as singular isothermal spheres (SISs) which is reasonable for galaxies. While SISs are not realistic when considering lensing by individual clusters, they are roughly sufficient for statistical purposes, as can be seen in the detailed study of [37,38]. The main advantage of SISs is of course that they allow for analytic results. Note that we are only interested in the magnification of the strongest image and not in the more detailed structure of the lensing event such as the time delays or the relative strength of images.

We concentrate on transient sources, such as compact binary coalescences which can be observed with LIGO/Virgo/Kagra [1–6,8–11]. The degeneracy between the intrinsic chirp mass and redshift of such systems may pose interesting questions. Many of the binary black holes (BBHs) detected by LIGO/Virgo are much more massive than those known from X-ray binaries in the galaxy (the latter have masses only up to ~25$M_{\odot}$) and in some cases, BBHs detected in GW have masses well above 60$M_{\odot}$ [6,8]. Black hole (BH) masses above ~25$M_{\odot}$ are expected in low-metallicity environments and are indeed predicted by various stellar evolution models (e.g., [39–48]). However, standard stellar evolution models struggle to explain masses in the so-called ‘pair-instability mass gap’ in the range between ~60$M_{\odot}$ and ~130$M_{\odot}$, where no BHs are expected to form [49]. An intriguing possibility could be that the apparent high masses are the result of redshifting, implying that at least some of the observed BBHs are high-redshift sources with low intrinsic masses that underwent significant gravitational lensing [16].

A detailed analysis of the individual events in Ref. [33] concluded that the probability of lensing for the sources observed by LIGO/Virgo runs $O(1)$ and $O(2)$ is vanishingly small, with the consequence that BH population models can be reasonably informed by GW observations to date [8,11,50]. In this work, we confirmed this finding with a statistical analysis. We present predictions for the probability that an event observed by LIGO/Virgo at a given distance has been magnified by a given amount: while the probability for more than one event is very small (the probability that two events are magnified by more than a factor 5 is approximately $10^{-10}$ for the typical distances of LIGO/Virgo events), it slowly decays with increasing magnification, i.e., there is a non-negligible probability of having extremely highly magnified events (the probability that an event at a distance larger than 2500 Mpc is magnified by a factor of more than 100 is approximately $10^{-7}$).

The remainder of this paper is structured as follows. In the next section, we present the lens model and we derive expressions for the cross-section, optical depth and lensing probability for a distribution of lenses. The results presented in this section are not new. The optical depth for strong lensing was already derived in the pioneering works of Turner et al. and later by Fukugita et al. and Oguri et al. [51–53]. Probability distributions for different lensing properties were derived from numerical N-body simulations in [54,55]. Our results are in agreement with this work. In Section 3, we derive the expression for
the probability that a source seen from a given *observed* distance was magnified by a given factor and we compute the mean magnification and its variance both as a function of the (unobservable) cosmological redshift and the observed distance which can be converted into an 'apparent' redshift under the assumption of ΛCDM cosmology with Planck values for the cosmological parameters. In Section 4, we introduce a numerical study for the population of lenses and sources and we present numerical results for the magnification as a function of both cosmological and observed redshift, for the LIGO/Virgo detector network. We also compute the probability of magnification for a source with given observed properties. In Section 5, we discuss our results and conclude. In a series of appendices, we discuss a simple analytical model which is useful to understand the physics and we present results for the mean magnification for different choices of source distribution.

**Notation:** In this paper, we set $c = 1$ but keep the gravitational constant $G$. We define $\theta$ and $\beta$ as the lensed and the unlensed angular position of the source. The optical axis connects the observer to the center of mass of the lens. We introduce a vector in the source plane from the optical axis to the source, $\mathbf{b} = d_s\beta$ and a vector $\zeta$ in the lens plane defined as $\zeta = \theta d_\ell$. This fixes the geometry. Here, $d_\ell = d(z_\ell)$, $d_\ell = d(z_\ell, z_s)$ and $d_s = d(z_s)$ are the angular diameter distances from the observer to the lens, from the lens to the source and from the observer to the source, respectively. These quantities are illustrated in Figure 1. A glossary with a list of the main symbols introduced in this article and their corresponding definitions can be found in Table 1.

<table><caption>**Table 1.** Glossary of the main quantities which we introduce and use in this paper.</caption><thead><tr><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td>$z_s$</td><td>(cosmological) source redshift</td></tr><tr><td>$z_\ell$</td><td>(cosmological) lens redshift</td></tr><tr><td>$D_s$</td><td>luminosity distance to the source</td></tr><tr><td>$D_\ell$</td><td>luminosity distance to the lens</td></tr><tr><td>$d_s$</td><td>angular diameter distance to the source</td></tr><tr><td>$d_\ell$</td><td>angular diameter distance to the lens</td></tr><tr><td>$d_{\ell s}$</td><td>angular diameter distance from the lens to the source</td></tr><tr><td>$z_{\text{obs}}$</td><td>observed source redshift (inferred from standard ΛCDM assuming no magnification)</td></tr><tr><td>$\sigma_v$</td><td>galaxy velocity dispersion</td></tr><tr><td>$\mu$</td><td>magnification</td></tr><tr><td>$\Delta t$</td><td>time-delay between multiple images</td></tr><tr><td>$\sigma$</td><td>lensing cross-section</td></tr><tr><td>$\tau(\mu, z_s)$</td><td>optical depth for magnification by at least $\mu$ for source at $z_s$</td></tr><tr><td>$p(\mu, z_s)$</td><td>probability density of magnification for source at $z_s$</td></tr><tr><td>$P(> \mu, z_s)$</td><td>probability of magnification by at least $\mu$ for source at $z_s$</td></tr><tr><td>$F_{\text{lim}}$</td><td>limiting flux of a given observatory</td></tr><tr><td>$\rho_{\text{lim}}$</td><td>limiting signal-to-noise ratio (corresponding to $F_{\text{lim}}$)</td></tr><tr><td>$d\mathcal{N}_{\text{obs}}(F_{\text{lim}}, z_s)$</td><td># sources observed with liming flux $F_{\text{lim}}$ in a bin around $z_s$</td></tr><tr><td>$d\mathcal{N}_{\text{obs}}(F_{\text{lim}}, D_{\text{obs}})$</td><td># sources observed with liming flux $F_{\text{lim}}$ in a bin around $D_{\text{obs}}$</td></tr><tr><td>$d\mathcal{N}(\mu, F_{\text{lim}}, z_s)$</td><td># sources that we see with liming flux $F_{\text{lim}}$ in a bin around $z_s$, if the magnification is $\mu$</td></tr><tr><td>$d\mathcal{N}(\mu, F_{\text{lim}}, D_{\text{obs}})$</td><td># sources that we see with liming flux $F_{\text{lim}}$ in a bin around $D_{\text{obs}}$, if the magnification is $\mu$</td></tr><tr><td>$P_z(\mu)$</td><td>probability density that an observed source at (cosmological) redshift $z$ is magnified by $\mu$</td></tr><tr><td>$P_{\text{obs}}(\mu)$</td><td>probability density that an observed source at observed redshift $z_{\text{obs}}$ is magnified by $\mu$</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 1. Schematic representation of a gravitational lens: introduction of the relevant lengths and angles.</figcaption></figure>

## 2. Lens Modeling

We describe galaxies acting as lenses as singular isothermal spheres and we derive expressions for the magnification between multiple images as a function of the geometrical quantities entering the picture. Most of this section follows Refs. [51,52,56]. SISs are not a very realistic lens model but since we are interested in statistical averages and not in the detailed modeling of a single lens, the fact that individual lenses will in general be triaxial is expected to average out. Furthermore, the neglected substructure may actually lead to some additional magnification. In this sense, our model is rather conservative. However, a recent numerical study has shown that the SIS model is surprisingly close to the observed lensing probability [38]. Note also that our results are only sensitive to the magnification of the most highly magnified image but not to the geometry or the number of images.

### 2.1. Basic Strong Lensing Quantities

The matter density of a singular isothermal sphere is given by

$$ \rho(r) = \frac{\sigma_v^2}{2\pi Gr^2}, \qquad (1) $$

where $\sigma_v$ denotes the velocity dispersion of the lens. Despite the singularity in the center and the infinite total mass, this can be considered as a rather realistic mass distribution for lensing by a galaxy (see [34,56]), where $\sigma_v$ is the velocity dispersion inside the galaxy. Integrating along the line of sight, we obtain the surface density:

$$ \Sigma(\theta) = \frac{\sigma_v^2}{2G} \frac{1}{d_{\ell}\theta}. \qquad (2) $$

We introduce $\Sigma_c = (4\pi G)^{-1}d_s/(d_{\ell s}d_{\ell})$ so that the convergence $\kappa$ is given by

$$ \kappa(\theta) = \frac{\Sigma(\theta)}{\Sigma_c} = 2\pi\sigma_v^2 \frac{d_{\ell s}}{d_s} \frac{1}{\theta}. \qquad (3) $$

An SIS has a constant deflection angle given by

$$ \alpha(\theta) = 4\pi\sigma_v^2 \frac{d_{\ell s}}{d_s} = 2\theta\kappa(\theta) \equiv \alpha_0. \qquad (4) $$

We now rescale our variables by $\alpha_0$, defining the rescaled image and source positions as $x = \theta/\alpha_0$ and $y = \beta/\alpha_0$, respectively. In terms of these rescaled variables, the lens equation for the SIS simply becomes $y = x - x/|x|$. For $y = 0$, the solution is the Einstein ring $|x| = 1$, and for $y = |y| < 1$, there are two solutions: one with $x_1 = |x_1| = 1 + y$ on the same side of the line of sight as the source (positive parity); and one with $x_2 = |x_2| = 1 - y$ on the opposite side (negative parity). For $y > 1$, the second solution no longer exists. The Jacobian of the lens map is:
$$A_{ij} = \delta_{ij} \left(1 - \frac{1}{|\mathbf{x}|}\right) + \frac{x_i x_j}{|\mathbf{x}|^3}, \\ \det A = 1 - \frac{1}{|\mathbf{x}|}, \qquad \mu = \frac{1}{|\det A|} = \frac{|\mathbf{x}|}{|1 - |\mathbf{x}||}, \qquad (5)$$

where we introduced the magnification $\mu$. Expressing the total magnification of a point source at position $y$ in terms of $y$, we find:

$$\mu(y) = \begin{cases} \mu(\mathbf{x}_1) + \mu(\mathbf{x}_2) = \frac{y+1}{y} + \frac{1-y}{y} = \mu_1 + \mu_2 = \frac{2}{y}, & y \le 1, \\ \frac{y+1}{y} = 1 + \frac{1}{y}, & y \ge 1. \end{cases} \qquad (6)$$

## 2.2. Optical Depth and Lensing Probability

The impact parameter of the source (in the source plane) is given by $|b| = b = \beta d_s = y\alpha_0 d_s$. A source with an impact parameter smaller or equal to $b$ was amplified by a factor of at least $\mu(y)$. Hence, the cross-section for amplification by more than a factor of $\mu_1$ of the stronger image from a SIS lens with velocity dispersion $\sigma_v$ is:

$$\sigma(\mu_1, z_\ell, z_s, \sigma_v) = \pi b^2 = \pi (y\alpha_0 d_s)^2 = \frac{\pi (4\pi)^2 \sigma_v^4 d_{ls}^2}{(\mu_1 - 1)^2}. \qquad (7)$$

Note that this cross-section gives the area, centered along the line of sight of the lens, within which a source at $z_s$ must lie so that it is amplified by a factor $\mu_1$ or larger by the lens at $z_\ell$. Expression (7) also remains valid for $y \ge 1$, where we have only one image with magnification $\mu_1$ which tends to 1 when $y \to \infty$.

In our study of the strong lensing of gravitational waves, we only considered one image and not the sum of both since we expect to see a short burst of GWs which only comes from one image. The second image is delayed in time, with a typical time delay of the order of a few months (as can be seen in, e.g., [12]), which is much longer than the GW burst. Since we are interested in magnification, we shall compute the cross-section for the stronger image. This point has been neglected in the previous literature: in Equation (7), $y^{-1} = \mu_1 - 1$ is usually replaced by $y^{-1} = \mu/2 = (\mu_1 + \mu_2)/2$ which is the correct expression for a static situation in which both images are seen together. For strong amplification, $\mu_1 \sim \mu_2 \gg 1$—a difference which reduces the cross-section by a factor of 4.

To compute the optical depth to the source of magnification by at least $\mu$, denoted by $\tau(\mu, z_s)$, we need to know the physical density $n(\sigma_v, z_\ell)$ of lenses (galaxies) with a given velocity dispersion $\sigma_v$ at redshift $z_\ell$. We define the density function $n(\sigma_v, z_\ell)$ such that its integral $\int n(\sigma_v, z_\ell) d\sigma_v$ simply gives the total density of lenses at redshift $z_\ell$.

The optical depth for lensing with magnification $\ge \mu$ (of the most strongly magnified image in the case of two images) for a source at redshift $z_s$ is:

$$\tau(\mu, z_s) = \int_0^{z_s} dz_\ell \frac{dr}{dz_\ell} \left(\frac{d_\ell}{d_s}\right)^2 \int \sigma(\mu, z_\ell, z_s, \sigma_v) n(\sigma_v, z_\ell) d\sigma_v, \qquad (8)$$

where $dr$ is the physical length element at redshift $z_\ell$. Note also that we rescaled the cross-section to the lens redshift, $\sigma \to \sigma(d_\ell/d_s)^2$, since in Equation (8), we multiply by the lens density at $z_\ell$. Inserting this and (7) for the cross-section in Equation (8), we obtain:

$$\tau(\mu, z_s) = \frac{\pi (4\pi)^2}{(\mu - 1)^2} \int_0^{z_s} dz \frac{\chi^2(z, z_s) \chi^2(z)}{\chi^2(z_s) (1+z)^3 H(z)} \int_0^\infty d\sigma_v \sigma_v^4 n(\sigma_v, z), \qquad (9)$$

where $\chi(z_1, z_2)$ denotes the co-moving distance from the redshift $z_1$ to $z_2$ and $\chi(z) = \chi(0, z)$. Equation (9) is the optical depth for magnification larger than $\mu$.

Often, only the strong lensing case is considered and the magnification from the two images is added to give the total magnification. To do this, one has to replace $1/(\mu - 1)^2$
by $4/\mu^2$. As already mentioned above, we cannot do this herein since in the case of strong magnification and double images, we expect a considerable time delay, so that typically, we observe only one image at one given time. Here, we assume this to be the stronger image. In the strong magnification case, $\mu \gg 2$, this difference is roughly a factor 4, while in the limit $y \to 1$ where $\mu \to 2$ and the second image disappears, the two expressions agree.

The probability of having an event from redshift $z_s$ magnified by more than $\mu$ is:

$$P(> \mu, z_s) = 1 - \exp(-\tau(\mu, z_s)) = \int_{\mu}^{\infty} p(\mu, z_s) d\mu, \qquad (10)$$

where $p(\mu, z_s)$ is the probability density for magnification by $\mu$:

$$p(\mu, z_s) = -\frac{d}{d\mu} P(> \mu, z_s) = -\frac{d\tau}{d\mu} \exp(-\tau(\mu, z_s)). \qquad (11)$$

To understand Equation (10), note that $d\tau/dz$ is the scattering rate leading to a magnification larger than $\mu$ (per redshift). Hence the probability for magnification larger than $\mu$ satisfies the differential equation $dP(> \mu, z)/dz = (1-P)d\tau/dz$ with solution (10). The factor $(1-P)$ is essential to keep the probability normalized also when $\tau$ becomes large. In the limit of small optical depth, $P(> \mu, z_s) = \tau(\mu, z_s)$.

Note that, while our formula also holds if there is only one image, $\mu < 2$, it cannot describe the de-magnification which happens when a beam (from a GW or a light source) crosses an under-density. An SIS is always an over-density. Therefore, our formulas are restricted to $\mu > 1$.

The probability density and probability defined in this section only depend on the lens model and lens distribution but not on the distribution of sources or on the sensitivity of the instrument. They are detector-independent functions which depend on the cosmological source redshift, which is a quantity that cannot be accessed observationally (without assuming the existence of an optical counterpart) as it is degenerate with the (unknown) magnification. The physically interesting quantity is rather the probability that a source that we observe from a given observed distance has been magnified by a given amount. We introduced this quantity in the next section and we will see that it depends on both the distribution of sources and the specifications of a given instrument.

### 3. Inclusion of Lensing Selection Effects

In this section, we derive expressions for the probability that a given source observed from a given distance has a magnification $\mu$, as a function of both redshift and observed source distance. Several ingredients enter the computation: the probability distribution of magnification given a lens population $p(\mu, z)$ given by Equations (9) and (11), the underlying source population and the sensitivity curve of a given instrument. As an instrument has finite sensitivity, highly magnified sources are easier to detect than non-magnified ones: hereafter, we dub these effects *lensing selection effects*. These results are valid for a generic distribution of magnification $p(\mu, z)$, in particular, they also stay valid when one adds weak lensing contribution to the distribution of strong lensing studied in this work. We will come back to this point in the discussion session.

#### 3.1. Magnification Distribution as Function of Redshift

Compact binary coalescences can be detected using matched filtering techniques, where a model of the signal, which depends on the physical parameters of the source, is matched against the observed signal (see [9] for a description and further references). For a given detector sensitivity and a given template bank used to perform the matched filtering, the obtained SNR $\rho$ can be related to a false alarm rate (FAR), which in turn sets a threshold for detection. For example, the O3a LIGO/Virgo catalog [9] includes sources with the FAR
below $2 y^{-1}$. For simplicity, we will work directly with a limiting SNR $\rho_{\lim}$ and assume that the signal is detected if $\rho \ge \rho_{\lim}$, where the SNR $\rho$ is defined by [57]

$$ \rho^2 = 4 \int \frac{|h(f)|^2}{S_n(f)} df = \int \rho^2(f) df. \qquad (12) $$

Here, $S_n(f)$ is the noise of the detector at frequency $f$, as can be seen in [58].

If an event is magnified by $\mu > 1$, we observe it even if without magnification, its SNR would have been $\rho_{\lim}/\sqrt{\mu} < \rho_{\lim}$. It follows that in the presence of magnification, the total number of objects which we expect to see from a redshift bin $dz$ is given by the following convolution:

$$ d\mathcal{N}_{\text{obs}}(z, \rho_{\lim}) = \int_{0}^{\infty} d\mu p(\mu, z)\mathcal{N}(\mu, z), \qquad (13) $$

where $p(\mu, z)$ is the distribution of lensing magnification as a function of cosmological redshift and:

$$ d\mathcal{N}(\mu, z) = dz \int_{\rho_{\lim}/\sqrt{\mu}}^{\infty} N(\rho, z) d\rho, \qquad (14) $$

is the number of sources that we see from redshift $z$ in the presence of magnification $\mu$. In Equation (14) $N(\rho, z)$ is the number density of sources as a function of redshift and $\rho$. We denote by $\rho$ the SNR of a source computed assuming no lensing effects. The mean amplification of a source at redshift $z$ is:

$$ \langle \mu \rangle (z) = \frac{\int_{0}^{\infty} d\mu \, \mu \, p(\mu, z) d\mathcal{N}(\mu, z)}{\int_{0}^{\infty} d\mu \, p(\mu, z) d\mathcal{N}(\mu, z)}. \qquad (15) $$

Note that (15) is the first moment of the probability distribution of magnification:

$$ \mathcal{P}_z(\mu, z) \equiv \mathcal{C} \, p(\mu, z) \frac{d\mathcal{N}(\mu, z)}{dz}, \qquad (16) $$

where $\mathcal{C}$ is a normalization constant. Equation (16) depends on three basic ingredients: (i) distribution of magnification $p(\mu, z)$; (ii) distribution of sources as a function of redshift and luminosity (or analogously $\rho$); (iii) sensitivity curve of a given experiment which determines $\rho_{\lim}$.

For sufficiently high $\mu$, $p(\mu, z_s) \simeq -d\tau/d\mu \propto \mu^{-3}$, while the number of sources which we see with magnification $\mu$ or more from a redshift bin around $z_s$, defined in Equation (14), becomes independent of $\mu$ for large values of $\mu$. Indeed, in the presence of very high magnification, we simply see all sources that are present. This last statement is general and does not depend on the details of the distribution of sources. It follows from the physical assumption that there exists a minimum value for the intrinsic luminosity emitted, $L_{\min}$, and once $\mu$ is such that $L_{\lim}/\mu < L_{\min}$, we see all sources and the luminosity integral no longer depends on $\mu$. The probability of having magnification bigger than $\mu$ for a given source redshift is given by the integral of the probability distribution of magnification, Equation (16), from a minimum value $\mu$. We denote this quantity as $\mathcal{P}_z(> \mu)$. For large values of magnification $\mathcal{P}_z(> \mu) \propto \mu^{-2}$, the variance of the distribution, defined by $\text{Var}(\mu) \equiv \langle \mu^2 \rangle - \langle \mu \rangle^2$, is log-divergent. We stress that this is a general property of the probability distribution of magnification in strong lensing. It does not depend on the details of the distribution of sources, nor on the details of the lens distribution. Moreover, it also applies to the case of transient electromagnetic sources lensed by galaxies. Of course, this divergence is not physical. When the magnification becomes very large, $y$ is very small and the geometric optics treatment adopted here breaks down. Nevertheless, this indicates that the distribution has a significant tail with large magnification. Even though, as we shall see, the probability for a magnification of, e.g., 10, is very small; it is not exponentially small as it would be for a Gaussian distribution, but it only decays as a mild power law.
## 3.2. Magnification Distribution as Function of Observed Distance

The functions determined thus far, however, are not observable as we typically do not know the intrinsic cosmological redshift of a GW event. We usually infer it by assuming that the observed luminosity distance is the one of the background universe without magnification. This is of course not true when magnification is relevant. We therefore rewrite (13) and (15) as functions of the observed luminosity distance which is directly extracted from observations. It is given by

$$D_{\text{obs}}(z_s, \mu) \equiv \frac{D_s(z_s)}{\sqrt{\mu}}. \qquad (17)$$

The number of sources that we can see from an observed luminosity distance $D_{\text{obs}} \pm dD_{\text{obs}}/2$ is obtained from (13), rewriting all quantities as functions of the observed distance and replacing the redshift bin $dz_s$ by $(dD_{\text{obs}}/dz_s)^{-1}dD_{\text{obs}} = \sqrt{\mu}(dD_s/dz_s)^{-1}dD_{\text{obs}}$. This gives:

$$d\mathcal{N}_{\text{obs}}(D_{\text{obs}}, \rho_{\text{lim}}) = \int_{0}^{\infty} d\mu p(\mu, z_s(D_{\text{obs}}, \mu)) d\mathcal{N}(\mu, D_{\text{obs}}), \qquad (18)$$

where:

$$d\mathcal{N}(\mu, D_{\text{obs}}) = dD_{\text{obs}} \frac{\sqrt{\mu}}{D'_{s}(z_s)} \int_{\rho_{\text{lim}}/\sqrt{\mu}}^{\infty} N(\rho, z_s(D_{\text{obs}}, \mu)) d\rho, \qquad (19)$$

is the number of sources which we see from a luminosity distance $D_{\text{obs}} \pm dD_{\text{obs}}/2$, with magnification $\mu$. In these expressions, $z_s(D_{\text{obs}}, \mu) \equiv z_s(D_s = \sqrt{\mu}D_{\text{obs}})$ is the cosmological redshift that corresponds to the observed distance $D_{\text{obs}} = \sqrt{\mu}D_s(z_s)$. Not to be confused with the redshift inferred by an observer who assumes that $D_{\text{obs}}$ is the luminosity distance of the background universe (without magnification). This latter, which we denote by $z_{\text{obs}}$, is given by $D_s(z_{\text{obs}}) = D_{\text{obs}}$. In Equation (19) $D'_{s}(z_s) \equiv dD_s/dz_s$. The average amplification as a function of the observed luminosity distance is therefore:

$$\langle \mu \rangle (D_{\text{obs}}) = \frac{\int_{0}^{\infty} d\mu \mu p(\mu, D_{\text{obs}}) d\mathcal{N}(\mu, D_{\text{obs}}) / dD_{\text{obs}}}{\int_{0}^{\infty} d\mu p(\mu, D_{\text{obs}}) d\mathcal{N}(\mu, D_{\text{obs}}) / dD_{\text{obs}}}, \qquad (20)$$

which of course depends on the threshold $\rho_{\text{lim}}$, i.e., on the detector sensitivity. In full analogy with Equation (16), we introduce the probability distribution of magnification as a function of the observed luminosity distance of the source:

$$\mathcal{P}_{\text{obs}}(\mu, D_{\text{obs}}) \equiv C p(\mu, D_{\text{obs}}) \frac{d\mathcal{N}(\mu, D_{\text{obs}})}{dD_{\text{obs}}}, \qquad (21)$$

where $C$ is a constant of normalization which is given by the denominator of (20). Equations (13)–(16) and Equations (18)–(21) are the main results of this section.

From Equation (19), it follows that $d\mathcal{N}(\mu, D_{\text{obs}})/dD_{\text{obs}} \propto \sqrt{\mu}$ for $\mu \gg 1$ as both the lower integration bound and $D'_{s}$ tend towards constant asymptotic values for large magnification $\mu$. The probability density $p(\mu, D_{\text{obs}})$ in (21) also depends on $\mu$ through the observed luminosity as we defined $p(\mu, D_{\text{obs}}) \equiv p(\mu, z_s(\sqrt{\mu}D_{\text{obs}}))$. Recalling that $p(\mu, z_s) \propto \mu^{-3}$, it follows that $p(\mu, D_{\text{obs}})$ will have a milder decay with magnification. Indeed, for a fixed $D_{\text{obs}}$ and $\mu > 1$, $p(\mu, D_{\text{obs}}) \equiv p(\mu, z_s(\sqrt{\mu}D_{\text{obs}}))$ also receives contributions from redshifts higher than $z_s(D_{\text{obs}})$ and therefore the lensing probability is higher. The exact scaling of $p(\mu, D_{\text{obs}})$ with $\mu$ depends on the details of the lens distribution.

We anticipate here some results that we will obtain in Section 4.1 when considering a realistic model for the galaxy distribution and evolution. For a realistic distribution with co-moving number density of lenses which reproduces observations (or hydrodynamical simulations calibrated to observations), we find that $p(\mu, D_{\text{obs}}) \propto \mu^{-\alpha}$. The value of $\alpha$ increases from 2 to 3 as we increase the observed luminosity distance from ~300 Mpc to 5000 Mpc, where it saturates to the asymptotic value $\alpha = 3$. It follows that the probability
distribution (21) scales as $P_{\text{obs}}(\mu) \propto \mu^{-\alpha+1/2}$ and its integral from $\mu$ to infinity behaves as $P_{\text{obs}}(> \mu) \propto \mu^{-\alpha+3/2}$ for large values of $\mu$. It follows that, for any fixed value of $D_{\text{obs}}$, the variance of the distribution diverges as $\mu^{-\alpha+7/2}$. We stress that this is a general property of the probability distribution of magnification (for different observed redshift), independent from the details of the distribution of sources. Moreover, in that case, the variance of the distribution of magnification (for a given observed redshift or luminosity distance) is divergent. We observed that for realistic models of lens distribution, the mean of (21) also diverges for small observed distances where $\alpha \sim 2$. The fact that the variance of the distribution diverges is due to the presence of extremely highly magnified sources, i.e., sources aligned along the lens–observer direction for which the magnification $\mu \to \infty$, as can be seen in Equation (6). This is also at the origin of the divergence of the mean of the magnification at small observed distances. Physically, the magnification tends to a finite value and as mentioned above, the presence of the divergence is due to the geometric optics approximation, which breaks down in the vicinity of caustics, where wave effects need to be taken into account. The fact that the mean of the magnification diverges only for small observed distances is due to the fact that, if a given source is very highly magnified, it appears to be at a small observed distance $D_{\text{obs}}$. We will get back to this point in Sections 4 and 5.

# 4. Numerical Results

In this section, we model the redshift distribution and redshift evolution of lenses (galaxies) using fits to hydrodynamical simulations. We also model the distribution of GW sources using a phenomenological model of the binary BH population. We then compute the distribution of magnification as a function of the cosmological source redshift and of the observed one. We focus on the case of the two LIGO/Virgo detectors.

## 4.1. The Lens Distribution

We model the number density of galaxies (lenses) as a function of the velocity dispersion $\sigma_v$, taking into account the evolution of galaxies with redshift. We used the results of Torrey et al. [59], based on the Illustris hydrodynamical simulation, to describe the number density of galaxies as a function of their velocity dispersion. The Illustris simulation [60–62] is a cosmological hydrodynamical simulation run in a periodic box of size $L = 106.5$ Mpc, dark matter mass resolution of $m_{\text{DM}} = 6.3 \times 10^6 M_{\odot}$ and baryon mass resolution of $m_{\text{bar}} = 1.3 \times 10^6 M_{\odot}$. Illustris includes prescriptions for gas cooling, star formation and feedback, chemical evolution of the interstellar medium and supermassive black hole growth and feedback. While many of these physical processes are described with sub-grid models, the simulation was calibrated to match the present day stellar mass function and the cosmic star formation rate history, and was shown to reproduce the observed stellar mass function out to a high redshift. Torrey et al. [59] calculated the cumulative velocity dispersion function of the Illustris galaxies $N(> \sigma_v, z)$, where the velocity dispersion $\sigma_v$ is the three-dimensional standard deviation of the stellar velocity within the stellar half-mass radius of the galaxy, and provided fitting functions for this distribution. We use the values in Table 6 in the ArXiv version of Ref. [59], and the fit is [63]:

$$ \log_{10} N(> \sigma_v, z) = A(z) + \alpha(z) (\log_{10} \sigma_v - \gamma(z)) + \\ + \beta(z) (\log_{10} \sigma_v - \gamma(z))^2 - (\sigma_v \times 10^{-\gamma(z)})^{1/\ln(10)}, \quad (22) $$

where the numerical value of $\sigma_v$ is to be taken in units of km/sec. We stress that, here, $\sigma_v$ represents the redshift-independent velocity dispersion. The functions $A$, $\alpha$, $\beta$ and $\gamma$ are polynomial functions of redshift given in Appendix A. While these simulations cannot resolve sub-structures smaller than $10^6 M_{\odot}$, we expect these to have a relatively small effect on the total magnification (in Ref. [38]; it was found that a resolution of approximately $10^{11} M_{\odot}$ is sufficient for a good estimation of the total magnification). From the comparison
with the simulations, we expect the fit to be accurate within a factor of at least 2. We define the co-moving number density of galaxies as

$$n^{\text{com}}(\sigma_v, z) = -\frac{dN}{d\sigma_v}(> \sigma_v, z). \qquad (23)$$

The result for different redshifts is shown in Figure 2.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Co-moving number density of galaxies as a function of the velocity dispersion $\sigma_v$ for different redshifts. The color code is: $z = 0, 1, 2, 3, 4, 5, 6$ from black to red. The figure is obtained from fits to the Illustris simulation of [59].</figcaption></figure>

The analysis of the Illustris simulation [59] which we use here revealed a very mild evolution in the galaxy velocity dispersion function for redshifts $0 \le z \le 2$, and a significant evolution beyond $z = 2$. Interestingly, the level of redshift evolution of the velocity dispersion function was found to be much lower than for the stellar mass function.

It is important to stress that the galaxy velocity dispersion distribution function that we use in this work depends on galaxy build-up and therefore the details of the baryonic physics implementation in the Illustris simulation. However, the analysis in [59] found that the median number density evolution of galaxies is nearly identical when traced using the stellar velocity dispersion of dark matter halo mass. It is therefore reasonable to assume that the results are relatively insensitive to the specific implementation of baryonic physics, and mostly depends on the dark matter halos, whose evolution is a rather robust outcome of the simulation.

## 4.2. Results for Magnification Probability Density

In Figure 3, we plot the (normalized) lensing probability density as a function of the magnification for different redshifts. On the left hand panel, we present curves for fixed-source redshifts, while in the right hand panel, we fix the observer distance $D_{\text{obs}} = D(z_s)/\sqrt{\mu}$. Note that the probability distribution is always strongly peaked at $\mu \simeq 1$. We stress that both optical depth and probability density only depend on the lens model and on the model for the distribution of lenses, but they are observatory-independent and they do not depend on the distribution of sources. Similar results have been obtained in the past, as can be seen in [38,54,55,64] and are in good quantitative agreement with our findings.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** (Left) Lensing probability density $p(\mu, z)$ as a function of magnification for $z = 0.6, 1, 2, 3$ (from bottom to top) using the Illustris result (22). (Right) Lensing probability density $p(\mu, D_{\text{obs}})$ as a function of magnification for $D_{\text{obs}} = 500, 1000, 2500, 4500$ Mpc (from bottom to top) using the Illustris result (22). For a sufficiently high redshift, $p(\mu, z)$ scales as $\mu^{-3}$ as mentioned in the text. The thin solid line on the left panel illustrates this scaling.</figcaption></figure>

### 4.3. Compact Binary Population

We focused below on BBH events and assumed that their masses are in the range of $[5, 50] M_{\odot}$ [8,50]. We further assume that their formation rate is proportional to the cosmic star formation rate (SFR), for which we use the following parametric form [65]:

$$ \Psi(t(z)) = A \frac{ae^{b(z-z_m)}}{a-b+b \cdot e^{a(z-z_m)}}, \qquad (24) $$

with $A = 0.24 M_{\odot}/\text{yr}/\text{Mpc}^3$, $z_m = 2.3$, $a = 2.2$ and $b = 1.4$ which peaks at $z_m$. Note that we use this parameterization in the redshift range $[0, 20]$ even though it was derived using observations up to $z = 9$. We define an efficiency parameter $\epsilon$ to describe the fraction of mass in BBH that merges within the age of the Universe that form out of a given stellar population. Moreover, we assume a delay time $t_{\text{delay}}$ between the formation of the progenitor star and the coalescence of the BBH, where the latter is distributed as

$$ P_d(t_{\text{delay}}) \propto t_{\text{delay}}^{-1}, \qquad (25) $$

in the range $[t_{\text{min}}, t_{\text{max}}]$ with $t_{\text{min}} = 50$ Myr and $t_{\text{max}}$ equal to the age of the Universe. We also assume that the mass of the primary BH is distributed as

$$ P(m_1) \propto m_1^{-\alpha}, \qquad (26) $$

in the range $[5, 50] M_{\odot}$ with $\alpha = 2.35$, while the mass of the secondary BH is uniformly distributed in the range $[5, m_1]$ so that $p(m_2) = \text{const.}$ in this range and $p(m_1, m_2) = p(m_1)p(m_2)$ [50]. Note that this mass distribution is independent of redshift. We assume zero spins for all BHs. The merger rate is then given by the following integral over the delay time distribution:

$$ \begin{aligned} \frac{dR}{dm_1 dm_2 dz} &= \\ &= \epsilon T_{\text{obs}} \int \Psi(t(z) - t_{\text{delay}}) p(m_1, m_2) P_d(t_{\text{delay}}) \frac{dV}{dz dt_{\text{delay}}}. \end{aligned} \qquad (27) $$

We calculate the SNR of each merger using Equation (12). We use the PhenomB inspiral-merger-ringdown waveforms to calculate the strain [66] and we use the noise power spectral density publicly available from [67,68]. We consider the event observable if it has $\rho \ge \rho_{\text{lim}} = 10$. The normalization parameter $\epsilon$ was obtained by comparing the number of observable events in our model to the number of events observed during the first two LIGO/Virgo observational runs O1 + O2; note that the most recent merger rate from the O3 LIGO/Virgo catalog is on the low end but consistent with the O1 + O2 rate:
$$N_{\text{obs}} = \int_{\rho(m_1, m_2, z) \ge 10} \frac{dR}{dm_1 dm_2 dz} dm_1 dm_2 dz \quad (28)$$

The merger rate per unit redshift and SNR is then calculated from Equation (27).

In Figure 4, we show the number of GW events per unit of SNR in different redshift bins. At low redshift, the slope and the width of the curves are constant; while at high redshift, the curves develop a low-SNR tail. The reason for this can be understood from Figure 5 where we plot the characteristic strain of BBH mergers with two different values of chirp mass (left and right panel) and for different redshifts, together with the O2 LIGO-Hanford detector noise (black solid) and the ET detector noise (black dashed). At low redshift, the SNR of high-mass events is always bigger than the one of low-mass events due to their higher amplitude. At high redshift, however, the sharp increase in the noise curve at low frequencies reduces the SNR of high masses more than the one of low masses which have higher frequencies so that the situation becomes actually reversed and low masses have a higher SNR. This leads to a spread of the SNR for different events over a large range. Of course, for all masses at high redshift, the SNR is very small and they can only be observed if they are highly magnified:

<figure><img src="image_8.png" /><figcaption>**Figure 4.** (Left) Number of events per unit of $\log_{10}\text{SNR}$. Each line corresponds to the number of events in a redshift bin $z \pm 0.025$ and the values of redshift are listed in the legend on the right. (Right) Number of events per unit of redshift. Each line corresponds to the number of events in a bin of signal-to-noise $\log_{10}\text{SNR} \pm 0.015$ and the values of $\rho$ are listed in the legend on the right of the figure.</figcaption></figure>

### 4.4. Probability of Strong Lensing for LIGO/Virgo Events

In Figure 6, we present the number density of magnified objects that we observe from a distance $D_{\text{obs}}$ as a function of magnification. We plot this quantity for three different values of $D_{\text{obs}}$. In the right-hand panel of Figure 6, we multiply each line of the left-hand panel by the probability density of lensing as a function of magnification to obtain the (normalized) probability distribution of magnification for objects observed from different distances, $D_{\text{obs}}$, i.e., Equation (21), as a function of magnification. Note that, while the number of objects which can be seen by the given survey (with fixed magnification) is somewhat lower for $D_{\text{obs}} = 3000$ Mpc than for $D_{\text{obs}} = 500$ Mpc, the probability of seeing such objects is higher for $D_{\text{obs}} = 3000$ Mpc than for $D_{\text{obs}} = 500$ Mpc. This is due to the fact that $p(\mu, D_{\text{obs}})$ is monotonically increases with $D_{\text{obs}}$, as can be seen in Figure 3. For the product of $p(\mu, D_{\text{obs}})dN/dD_{\text{obs}}$, the increase in $p$ wins over the decrease in $N$ for $D_{\text{obs}} = 3000$ Mpc.

In Figure 7, we plot the same quantities for different values of the cosmological redshift of the source. At a low redshift, $dN/dz$ tends towards a constant at high magnification since we see all events (black curve, left panel). The product $p(z, \mu)dN/dz$ nearly tends towards a constant at high redshift which indicates that $dN/dz$ grows similarly to $\mu^3$ at high redshift (blue dashed curve, right panel).
<figure><img src="image_9.png" /><figcaption>**Figure 5.** We fix in each panel a value of the chirp mass and we plot the strain for sources with that mass and in different redshift bins. The solid black curve is the O2 strain sensitivity of the LIGO-Hanford detector [67], and the dashed-dotted curve is the strain sensitivity of the Einstein Telescope [69].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** (Left) Number of events per unit of Mpc⁻¹ visible if the magnification is µ as a function of µ, for different observed luminosity distances: Dobs = (500, 3000, 4500) Mpc, indicated by red, black and orange curves, respectively. This quantity is given by dN/dDobs in Equation (19). (Right) The same quantity is multiplied by p(µ, Dobs) to obtain the probability distribution of magnification in Equation (21) per unit of observed distance. The little peaks in the red line come from the fact that when fixing Dobs = Ds(zs)/√µ, we pick up contributions from different zs, as can be seen in Figure 4. These are numerical artifacts from our binning of (dN/d log10 ρ dzs)(ρ, zs).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** (Left) Number of events visible if the magnification is µ per unit of redshift as a function of µ, for different cosmological redshift: z = (0.6, 1, 2, 3) indicated by black, red, green, blue lines, respectively. This quantity is given by dN/dzs in Equation (14). (Right) The same quantity multiplied by p(µ, zs) to obtain the probability distribution of magnification (per unit of redshift) in Equation (16).</figcaption></figure>

In Figure 8, we present the results for the average amplification factor as a function of the cosmological redshift and luminosity distance of the sources (left panel) and as a function of observed redshift and luminosity distance of the sources (right panel). No sources are observable beyond 5500 Mpc. (This can be understood by computing the
probability of seeing a source from a given observed distance, i.e., integrating $p(\mu, D_{\text{obs}})$ over $\mu$. The probability of seeing an object from a distance bigger than 5000 Mpc is smaller than $10^{-7}$.) The minimum amplification needed to detect sources from a redshift $z_s$ and with a given $\rho_{\text{max}}(z_s) > \rho_{\text{lim}}$ is given by

$$ \mu_{\text{min}}(z_s) = \left( \frac{\rho_{\text{lim}}}{\rho_{\text{max}}} \right)^2, \qquad (29) $$

hence, we can see sources from a maximum observed distance given by

$$ D_{\text{obs}}^{\text{max}}(z_s) = \frac{D(z_s)}{\sqrt{\mu_{\text{min}}}} = \frac{D(z_s)\rho_{\text{max}}(z_s)}{\rho_{\text{lim}}}. \qquad (30) $$

For low redshifts, $z_s \lesssim 2$, and we can work out the scaling of this quantity with redshift using a simple analytical approach, as can be seen in Appendix C. The probability $P_{\text{obs}}(> \mu)$ that an object that LIGO/Virgo observes from a distance $D_{\text{obs}}$ was magnified more than $\mu$ is presented in Figure 9 for different values of $D_{\text{obs}}$ as a function of $\mu$.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** (Left) Amplification factor for events observable by LIGO/Virgo, as a function of the cosmological redshift (and luminosity distance) of a source. (Right) Amplification factor for events observable by LIGO/Virgo as a function of the observed redshift (luminosity distance) of a source. As explained in the text, no events are observed after 5500 Mpc, independently from the value of magnification.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Probability that a given GW event observed by LIGO/Virgo has a bigger magnification than a given value $\mu$, as a function of $\mu$, for different observed distances.</figcaption></figure>

## 5. Discussion and Conclusions

Due to cosmological expansion, the frequencies of gravitational waves are redshifted. Since the frequencies of the GW signal from BBH mergers are degenerate with the mass, the quantity that we directly measure from the signal is the redshifted chirp mass $M_z = (1+z)M$, where $M$ is the intrinsic (true) chirp mass of the binary at a redshift $z$. The estimated luminosity distance can be converted into a redshift estimate (by choosing a cosmological model), which can, in turn, be used to estimate the intrinsic chirp mass $M$ of the binary. The unknown lensing magnification will bias our estimation of the intrinsic
mass and the distance (or equivalently, the redshift) of the binary. Therefore, lensed binaries will appear as a population at a lower redshift and with higher masses. The question is whether this effect is relevant for a given population of GW sources or whether it can be neglected.

In this article, we proposed a semi-analytical approach to study the lensing of gravitational waves by galaxies. The framework we propose can be applied to any detector or detector network once the strain noise of the observatory is known. We presented results for the probability of strong lensing for BBH visible by LIGO/Virgo, as a function of the redshift and magnification and for the average amplification as a function of cosmological redshift and observed distances. The numerical results are in good qualitative agreement with a simple model (valid for low redshift sources) that we present in a series of appendices.

Results for the probability that a source with cosmological redshift $z_s$ is magnified by $\mu$ are presented in Figure 3 as a function of both the cosmological and observed redshift (or rather distance) of the source. We stress that this result does not depend on the properties of the distribution of sources, nor on the specifics of a given observatory. The probability is small, but it slowly decays as a function of magnification. More specifically, for large values of magnification, the probability as a function of cosmological redshift scales as $p(\mu, z_s) \propto \mu^{-3}$. This scaling is independent of the lens distribution model. The probability density, as a function of the observed distance, $p(\mu, D_{\text{obs}}) \equiv p(\mu, z_s(\sqrt{\mu}D_{\text{obs}}))$, has an even slower decay with magnification, as can be seen in the right panel of Figure 3. Indeed, for a fixed $D_{\text{obs}}$ and $\mu > 1$, $p(\mu, D_{\text{obs}}) \equiv p(\mu, z_s(\sqrt{\mu}D_{\text{obs}}))$ takes contributions from a higher redshift than $z_s(D_{\text{obs}})$ and is therefore larger. The exact scaling of $p(\mu, D_{\text{obs}})$ with $\mu$ depends on the details of the lens distribution. For a realistic distribution with a co-moving number density of lenses which fits numerical simulations (Illustris), we find that $p(\mu, D_{\text{obs}}) \propto \mu^{-\alpha}$. The value of $\alpha$ increases from 2 to 3 as the observed luminosity distance varies from $\sim 300$ Mpc to 5000 Mpc where it saturates to the asymptotic value $\alpha = 3$.

Selection effects related to the sensitivity curve (and hence to the horizon) of a given detector network enter into the game when computing the probability that a source *that we can observe* with a given observatory has a given magnification. This probability distribution $P_{\text{obs}}(\mu)$ is defined in Equations (16) and (21) as a function of the cosmological redshift of a source and as a function of the observed luminosity distance, respectively. We focused on the LIGO/Virgo network in our analysis. Plots of the distribution of magnification for different cosmological redshifts and observed luminosity distances are presented in Figures 6 and 7 (right panels). The distribution highly peaks at $\mu = 1$ for low observed distances but it has more and more power at higher magnification as we move towards higher redshift or observed distances. The integral of this distribution multiplied by $\mu$ gives the average magnification expected for a given observed distance and the given experiment. For the LIGO/Virgo observatory, this is shown in Figure 8.

The probability distribution of magnification $P_{\text{obs}}(\mu)$ for a given observed distance scales as $P_{\text{obs}}(\mu) \propto \mu^{-\alpha+1/2}$ and its integral from $\mu$ to infinity behaves as $P_{\text{obs}}(> \mu) \propto \mu^{-\alpha+3/2}$ for large values of $\mu$ with $\alpha \in [2, 3]$. The probability $P_{\text{obs}}(> \mu)$ is presented in Figure 9 for different values of $D_{\text{obs}}$. From this scaling, it follows that the mean of the magnification at small observed distances, which have $\alpha \le 2.5$, diverges at large $\mu$. This is due to the presence of rare but extremely highly magnified sources, i.e., sources aligned along the lens–observer direction for which the magnification $\mu \to \infty$, see Equation (6). It also follows that, for any fixed value of $D_{\text{obs}}$, the variance of the distribution diverges as $\mu^{-\alpha+7/2}$. We stress that this is a general property of the probability distribution of magnification, which is independent from the details of the distribution of sources. Moreover, it also applies to the case of strongly lensed sources of electromagnetic radiation. Furthermore, in that case, the variance of the distribution of magnification (for different observed redshift) is divergent. Physically, the magnification tends to a finite value and as mentioned above, the presence of the divergence is due to the geometric optics approximation, which breaks down in the vicinity of caustics, where wave effects need to be taken into account. Nevertheless, this indicates that the distribution has a significant tail with large magnification. The fact that
the mean of the magnification only diverges for small observed distances is due to the fact that highly magnified sources appear to have small distances, $D_{\text{obs}} = D_s / \sqrt{\mu}$. In Figure 8, we present the results for the average amplification as a function of the cosmological and observed redshift of a source. A cut-off in $\mu$ was introduced to compute the mean for sources at low distances $D_{\text{obs}} < 1000$ Mpc. We checked that our cut-off choice is such that increasing its value does not significantly affect the result. The average magnification, as a function of the observed redshift, is very close to 1 for a low redshift and it starts growing from $z_{\text{obs}} = 0.5$, reaching the maximum value of 40 at a redshift 0.8 where we stop seeing sources, i.e., where the horizon of the observatory (computed to allow for the presence of magnification) lies.

In a series of appendices, we present a semi-analytic model which allows us to study the impact of changing the underlying astrophysical model for both BBH and lens distribution on the resulting magnification distribution. Using this model, we verified that the average magnification (mildly) depends on the details of the source population, but that the qualitative picture of Figure 8 is expected to be valid—independently from the details of the chosen astrophysical model.

Let us finally derive some consequences for the LIGO/Virgo O1, O2 and O3 observations. The 74 published observed BBH (or binary black hole-Neutron star) events have distances of $D_{\text{obs}} \in [300, 5300]$ Mpc. Assuming that the sensitivity has not significantly changed between O2 and O3, from the right panel of Figure 9, it follows that (for our modeling) the probability that one of them was magnified with a magnification of 10 or bigger is $P_{\text{obs}}(> 10) \sim 0.5$. (In Figure 9, the steep drop of the curves to the flatter slope happens between $1 < \mu < 1.01$.) However, the probability that the magnification of one of them is bigger than 100 is not much smaller, namely $P_{\text{obs}}(> 100) \sim 0.03$. In order to infer that most or at least several of the observed events originated from cosmological redshifts $z_s \sim 1-2$, these would have to be magnified by a factor of approximately 100. While this is not improbable, for one of the events to be substantially magnified, it would require very significant changes in the lens population, which are in contradiction with both the observations and numerical simulations, indicating that this cannot be the case for several events, details of which can be found in Appendix B. Another possibility would be a very significant change in the source distribution which would have to be approximately 1000 times larger at $z_s \sim 1-2$ and correspondingly smaller at $z \sim 0$. This, however, is at very strong odds with stellar evolution models. We therefore conclude that even though the lensing of GW observations will be relevant in the future, e.g., for the analysis of standard sirens, we cannot reasonably 'tweak' galaxy population models or stellar evolution models so that we can conclude that most of the observed LIGO/Virgo detections are strongly lensed, as proposed in [16,17,27,36].

The framework used in this analysis is flexible and can be adapted to study the strong and weak lensing of sources in any frequency band, and it can be applied to any detector network. In this work, we have focused on the case of the LIGO/Virgo network. We stress that while for the case of LIGO/Virgo, the contribution of weak lensing is negligible, this is not the case for instruments with a sensitivity such that the peak of the redshift distribution of sources is typically within the instrument horizon. In this case, the weak lensing contribution to the distribution of $\mu$ has to be added using, e.g., the results of ray tracing simulations. Then, a $p(\mu, z)$ which describes both weak and strong lensing can be used to derive the corresponding distribution of $P(\mu)$ and $P_{\text{obs}}(\mu)$ (including selection effects) following the procedure illustrated in this work, as can be seen in, e.g., [70], for a recent analysis of the weak and strong lensing of massive black hole binaries visible by LISA.

**Author Contributions:** Conceptualization, R.D. and G.C. All three authors contributed to the required theoretical calculations and to the writing of the paper. All authors have read and agreed to the published version of the manuscript.
**Funding:** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 693024) and from the Swiss National Science Foundation.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** It is a pleasure to thank Michele Maggiore and Nicola Tamanini for helpful discussions. We also thank Paul Torrey for providing us with the analytic fitting formula of [59]. We are grateful also to XiKai Shan who pointed out a relevant error in the first version of this paper.

**Conflicts of Interest:** The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## Appendix A. Fit to the Illustris Simulation

The functions $A$, $\alpha$, $\beta$ and $\gamma$ of the fitting formula used in Section 4.1 to describe the evolution of galaxy as a function of redshift and velocity dispersion [63] are modeled as

$$A = A_0 + A_1z + A_2z^2, \qquad (A1)$$

$$\alpha = \alpha_0 + \alpha_1z + \alpha_2z^2, \qquad (A2)$$

$$\beta = \beta_0 + \beta_1z + \beta_2z^2, \qquad (A3)$$

$$\gamma = \gamma_0 + \gamma_1z + \gamma_2z^2. \qquad (A4)$$

In the following table, we list the values of the coefficients.

<table><thead><tr><th>$i$</th><th>$i = 0$</th><th>$i = 1$</th><th>$i = 2$</th></tr></thead><tbody><tr><td>$A_i$</td><td>7.391498</td><td>5.729400</td><td>-1.120552</td></tr><tr><td>$\alpha_i$</td><td>-6.863393</td><td>-5.273271</td><td>1.104114</td></tr><tr><td>$\beta_i$</td><td>2.852083</td><td>1.255696</td><td>-0.286638</td></tr><tr><td>$\gamma_i$</td><td>0.067032</td><td>-0.048683</td><td>0.007648</td></tr></tbody></table>

## Appendix B. Analytic Description of Lens Distribution

If there was no evolution of lenses, we would expect the physical number density of galaxies to scale with redshift as $n(\sigma_v, z) \propto (1+z)^3$, hence:

$$\int_{0}^{\infty} d\sigma_{v} \sigma_{v}^{4} n(\sigma_{v}, z) = (1 + z)^{3} N \langle \sigma_{v}^{4} \rangle, \qquad (A5)$$

where $N$ is the present galaxy density:

$$N = \int_{0}^{\infty} d\sigma_{v} n(\sigma_{v}, z = 0), \qquad (A6)$$

and $\langle \sigma_v^4 \rangle$ is the mean of the velocity dispersion to power 4. Let us now assume that, in the considered redshift range, it is legitimate to neglect evolution, and hence the physical density scales with the expansion of the Universe. We postulate this law of evolution, not because it is very realistic, but because it allows us to analytically integrate the optical depth. However, as we shall see in the next section, for low redshifts to which present gravitational wave experiments are sensitive, this model is actually quite good.
Neglecting evolution, the factor $1/(1+z)^3$ in (9) cancels out and using $\chi(z, z_s)' = -H(z)^{-1}$, we obtain:

$$ \begin{aligned} \tau(\mu, z_s) &= \frac{\pi(4\pi)^2 \langle\sigma_v^4\rangle N}{\chi_s^2 (\mu-1)^2} \int_0^{\chi_s} d\chi (\chi_s - \chi)^2 \chi^2 \\ &= \frac{\pi(4\pi)^2 \langle\sigma_v^4\rangle N \chi_s^3}{30(\mu-1)^2}. \end{aligned} \quad (A7) $$

Crude estimates for $N$ and $\sigma_v$ are (our value of $\sigma_v$ is one of large galaxies and $N$ is obtained by simply dividing $\rho_m = \Omega_m\rho_c$ by the mass of a large galaxy, $M_{\text{gal}} \sim 10^{12} M_{\odot}$.)

$$ N = 10^9 H_0^3, \quad \langle\sigma_v^4\rangle = (150 \text{ km/s})^4 \approx 5 \times 10^{-14}. \quad (A8) $$

As we set $c=1$, we have $H_0 = h/(3000 \text{ Mpc})$. With this, we obtain:

$$ \tau(\mu, z_s) \approx \frac{0.001}{(\mu-1)^2} (H_0 \chi(z_s))^3 \left[ \frac{N}{10^9 H_0^3} \frac{\langle\sigma_v^4\rangle}{5 \times 10^{-14}} \right]. \quad (A9) $$

for the optical depth.

This analytic model gives a result for optical depth in good agreement with the one obtained in Section 4.2 considering a more realistic distribution of lenses, which evolves with redshift (fractional deviations of a few percent). For our simple model, we find for the probability density of magnification as a function of the cosmological source redshift:

$$ p(\mu, z_s) = \frac{2p_1(z_s)}{(\mu-1)^3} \exp\left(-\frac{p_1(z_s)}{(\mu-1)^2}\right), \quad \text{where} \quad (A10) $$

$$ p_1(z_s) = 0.001(H_0\chi(z_s))^3 \left[ \frac{N}{10^9 H_0^3} \frac{\langle\sigma_v^4\rangle}{5 \times 10^{-14}} \right]. \quad (A11) $$

Already from this expression, it is clear that we cannot achieve a significant probability of amplification at low redshifts. With this probability distribution, the mean magnification for a source from redshift $z_s$ is $\langle\mu(z_s)\rangle = 1 + \sqrt{\pi p_1(z_s)}$. To obtain a non-negligible probability for magnification by a factor of 50 or more, which is needed to explain the observed GW events as highly magnified events from redshift $z_s \sim 1-2$, we require that $p_1 \gtrsim 1000$. The factor $H_0\chi(z_s)$ in Equation (A11) becomes of order unity at redshift $z_s \sim 1.5$ and then very slowly grows up to $H_0\chi(z_s) \sim 2$ at $z_s \sim 10$. This requires a galaxy density which is more than $10^6$ times higher than $10^9 H_0^3$ or $\sigma_v \gtrsim 4500 \text{ km/s}$ or a combination of both such that:

$$ \frac{N}{10^9 H_0^3} \frac{\langle\sigma_v^4\rangle}{(150 \text{ km/s})^4} > 10^6. \quad (A12) $$

These possibilities are very strongly ruled out by present astronomical observations and by numerical simulations, as can be seen, e.g., in Section 4.

To have some intuition regarding the role of the redshift evolution of the galaxy velocity dispersion function, we compare the results obtained in Section 4 using the Illustris simulation to a non-evolving lens distribution. For this, we use the distribution function (22) for $z=0$, assuming no redshift evolution. We then have:

$$ n^{\text{no evo}}(\sigma_v, z) = (1+z)^3 n_{\text{com}}(\sigma_v, z=0). \quad (A13) $$

The comparison is shown in Figure A1. The result is valid for arbitrary magnification $\mu$.
<figure><img src="image_14.png" /><figcaption>Figure A1. Relative difference between the probability with and without redshift evolution. The result is valid for arbitrary magnification $\mu$.</figcaption></figure>

At low redshift, $0 < z < 2$, the non-evolving distribution actually has somewhat less galaxies with low-velocity dispersion, which are the most numerous ones, and therefore the evolving distribution leads to a higher probability for magnification than the non-evolving one, $\delta P > 0$. At higher redshift $z > 4$, the loss of lenses in the evolving distribution becomes more relevant and $\delta P$ changes sign. However, the difference between the evolving and non-evolving distribution is never larger than 10% for redshift $z \le 6$. This tells us that our results are quite insensitive to the redshift evolution of the lenses.

## Appendix C. Analytical Description of BBH Population

In this appendix, we work out an analytic model to describe the number density of sources as a function of the emitted gravitational-wave strain amplitude (or the SNR for a given detector configuration) and as a function of redshift. As a first step, we review the derivation of the SNR of a binary system of compact objects at fixed redshift as a function of the component masses and of the noise power spectral density of the instrument(s). We introduce the SNR per unit frequency $\rho$ for a binary that merges at redshift $z_s$ and is observed in a narrow frequency band around $f$ as [57]

$$ \rho^2(f, z_s) = 4 \frac{|h(f, z_s)|^2}{S_n(f)} \propto \frac{F(f, z_s)}{S_n(f)}, \qquad (A14) $$

where $S_n(f)$ is the noise power spectrum of a given detector network, with a dimension of $Hz^{-1}$ and $h(f, z_s)$ is the Fourier transform of the strain, with dimensions $Hz^{-1}$. The total SNR in the frequency band of a given detector is simply the integral of Equation (A14) over that band. In Equation (A14), we defined the received flux per unit of frequency as $F(f, z_s)$. Hence, there is a one-to-one mapping between flux and SNR. For a single binary system, the SNR from the inspiral phase integrated over observed frequency is given by [57] (neglecting the merger and ringdown phases and writing the quantities upon which the SNR explicitly depends:)

$$ \rho^2(z_s, m, \mathcal{M}) = \frac{5}{96\pi^{4/3}} \frac{\Theta^2}{D^2(z_s)} (GMz)^{5/3} f_{7/3}(z_s, m), \qquad (A15) $$

where $M_z = M(1 + z_s)$ is the redshifted chirp mass of the system, $m = m_1 + m_2$ is the total mass and the integral over frequencies can be approximated by

$$ f_{7/3}(z_s, m) \equiv \int_0^{f_{\text{ISCO}}/(1+z_s)} df [f^{7/3}S_n(f)]^{-1}, \qquad (A16) $$

with dimension $[f_{7/3}] = [Hz]^{-1/3}$, while $\Theta^2$ (with $0 \le \Theta^2 \le 16$) is a geometrical factor that depends on the inclination of the binary and on the antenna pattern of the detector. Averaging over all angles, one finds that $\Theta^2 = 64/25$. We do not need this explicit expression in what follows and the interested reader is referred to [57] for details. In Equation (A16),
the upper integration bound is given by the (redshifted) frequency corresponding to the innermost stable circular orbit (ISCO) of the system, i.e., the frequency at which we consider the inspiraling phase of the system to end in our approximation. It is defined as [58]

$$f_{\text{ISCO}} \equiv \frac{1}{6\sqrt{6}(2\pi)} \frac{c^3}{Gm} \approx 2.2 \, \text{kHz} \left(\frac{M_{\odot}}{m}\right), \qquad (\text{A17})$$

where *m* is again the total mass of the binary system. This frequency is not the frequency where the signal is maximal, but it is the frequency beyond which we can no longer trust the quadrupole formula used for the SNR estimate in Equation (A15). It is therefore a conservative estimate of the SNR which has contributions from higher frequencies, e.g., the coalescence frequency is typically approximately 3$f_{\text{ISCO}}$. However, we follow here the standard procedure as detailed in [57]. We rewrite the parameterization (A15) in a way which will be useful for what follows:

$$\begin{aligned} \rho^2(z_s, \mathcal{M}, m) &= \\ &= 7.7 \cdot 10^{-39} \Theta^2 \left(\frac{\mathcal{M}}{M_{\odot}}\right)^{5/3} \left(\frac{\text{Mpc}}{D(z_s)}\right)^2 \frac{f_{7/3}(z_s, m)}{\text{Hz}^{-1/3}} (1+z_s)^{5/3}. \end{aligned} \qquad (\text{A18})$$

Let us understand the redshift dependence: from Equation (A16), we see that if $f_{\text{ISCO}}/(1+z_s) > f_{\text{max}}$, then $f_{7/3}(z_s)$ does not depend on redshift. Here, with $f_{\text{max}}$, we denote the maximum frequency a given detector is sensitive to (i.e., the frequency above which the noise becomes very large) and it is fixed by the form of the spectral noise in Equation (A16). This is determined by the condition:

$$(1+z_s) < \frac{f_{\text{ISCO}}}{f_{\text{max}}} \approx \left(\frac{2.2 \cdot \text{kHz}}{f_{\text{max}}}\right) \left(\frac{M_{\odot}}{m}\right). \qquad (\text{A19})$$

If this condition is verified for a given source, then (A16) is redshift independent.

From a comparison with numerical models for source formation and evolution (see Section 4), we find that the maximum SNR that we receive from a given redshift $z_s$ can be effectively written as

$$\rho_{\text{max}}(z_s) = \mathcal{A} \left( \frac{M_{\text{max}}}{M_{\odot}} \right)^{5/6} \frac{\text{Mpc}}{D(z_s)} (1+z_s)^{5/6} \bar{f}_{7/3}(z_s)^{1/2}, \qquad (\text{A20})$$

and similarly for the minimum value:

$$\rho_{\text{min}}(z_s) = \mathcal{A} \left( \frac{M_{\text{min}}}{M_{\odot}} \right)^{5/6} \frac{\text{Mpc}}{D(z_s)} (1+z_s)^{5/6} \bar{f}_{7/3}(z_s)^{1/2}, \qquad (\text{A21})$$

where $\bar{f}_{7/3}(z_s)^{1/2}$ is obtained using the average ISCO frequency in Hz, rescaled with redshift, as the upper cutoff in the integral (A16), i.e., the ISCO frequency for the average total mass *m*. We introduced a constant $\mathcal{A} = 1.4 \times 10^{-19}$ which corresponds to the square root of the numerical factor in Equation (A18) multiplied by the square root of $\bar{\Theta}^2 = 64/25$. We denoted by $M_{\text{min}}$ and $M_{\text{max}}$ the minimum and maximum value of the chirp mass (for a given mass distribution), respectively. We observed that for redshifts that satisfy the condition (A19), Equations (A20) and (A21) simply scale with redshift as $(1+z_s)^{5/6}D(z_s)^{-1}$.

As a second step, we need a parameterization for the number of sources per unit of SNR as a function of SNR and of redshift. The number density of sources per unit of chirp mass (and as a function of redshift) is given by

$$N(z_s, \mathcal{M}_c) \propto \mathcal{R}(z_s, \mathcal{M}_c) \propto p(\mathcal{M}_c), \qquad (\text{A22})$$

where $\mathcal{R}$ is the merger rate for a binary system with chirp mass $\mathcal{M}_c$ at redshift $z_s$, which in turn is proportional to the distribution of chirp masses (with a redshift-dependent
pre-factor). For simplicity, we assume it to be monomial, $p(M_c) \propto M_c^{-\beta_M}$. We now convert the number density (A22) into a number density per unit of SNR. Neglecting the mass dependence of theISCO frequency, as we did in (A20) and (A21), we can simply invert (A18) and we obtain for the distribution of sources as a function of SNR:

$$N(\rho, z_s) = \begin{cases} N_{\rho}(z_s)\rho^{-\gamma} & \text{for } \rho_{\min}(z_s) < \rho < \rho_{\max}(z_s), \\ 0 & \text{otherwise,} \end{cases} \qquad (A23)$$

with $\gamma = (6\beta_M - 1)/5$ and $N_\rho(z_s)$ which is a redshift-dependent function. For example, the distribution of masses in the astrophysical model discussed in Section 4 is such that the resulting distribution of chirp masses is (independent on redshift) $p(M_c) \propto M_c^{-3.2}$. Then, $N(\rho, z_s) \propto \rho^{-3.6}$, which is the scaling of the curves in Figure A2. Note that when comparing this figure to the scaling derived here, one needs to keep in mind that $dN/d\rho = \rho^{-1}dN/d\log\rho$.

<figure><img src="image_15.png" /><figcaption>Figure A2. Comparison between the analytic model for source distribution presented in Appendix C and the numerical simulations of Section 4.3, for different redshifts $0.0125 \le z \le 2.5075$. The dashed line represents the scaling given in Equation (A23) with $\gamma = 3.6$.</figcaption></figure>

## Appendix D. Analytic Derivation of the Average Magnification

The number of sources which we see from redshift $z_s$ with our limiting SNR $\rho_{\lim}$, if magnified by $\mu$, is:

$$N(\rho_{\lim}, z_s, \mu) = N_{\rho}(z_s) \int_{\alpha_{\lim}}^{\rho_2(z_s)} d\rho \rho^{-\gamma} \\ = \frac{N_{\rho}(z_s)}{\gamma - 1} \left[ \left(\frac{1}{\alpha_{\lim}}\right)^{\gamma - 1} - \left(\frac{1}{\rho_2(z_s)}\right)^{\gamma - 1} \right], \qquad (A24)$$

where $\alpha_{\lim} = \max(\rho_{\lim}/\sqrt{\mu}, \rho_1(z_s))$ and we set $\rho_1 = \rho_{\min}$, $\rho_2 = \rho_{\max}$. We assume that $\rho_{\lim}/\sqrt{\mu} < \rho_2(z_s)$. If this is not the case, then $N(\rho_{\lim}, z_s, \mu) = 0$. We assume here that $\gamma \neq 1$ but the case $\gamma = 1$ can be easily worked out following steps similar to the ones performed below. The result without magnification is simply $N(\rho_{\lim}, z_s, 1)$. The total number of sources which we see is therefore:

$$N(\rho_{\lim}, z_s) = \int_{1}^{\infty} d\mu p(\mu, z_s) N(\rho_{\lim}, z_s, \mu). \qquad (A25)$$

With a bit of algebra, this can be rewritten as

$$N(\rho_{\lim}, z_s) = \frac{N_{\rho}}{\gamma-1} (\rho_1^{-\gamma+1} - \rho_2^{-\gamma+1}) + \frac{N_{\rho}}{\gamma-1} \left[ \rho_2^{-\gamma+1} \int_1^{(\rho_{\lim}/\rho_2)^2} d\mu p(\mu, z_s) - \rho_1^{-\gamma+1} \int_1^{(\rho_{\lim}/\rho_1)^2} d\mu p(\mu, z_s) + \int_{(\rho_{\lim}/\rho_2)^2}^{(\rho_{\lim}/\rho_1)^2} d\mu p(\mu, z_s) \left(\frac{\sqrt{\mu}}{\rho_{\lim}}\right)^{\gamma-1} \right], \qquad (A26)$$
where the redshift dependence in $\rho_{1,2}$ and $N_\rho$ is understood. We eventually want to compute the average magnification. For this, we need to compute:

$$M(\rho_{\text{lim}}, z_s) \equiv \int_{1}^{\infty} d\mu \, \mu \, p(\mu, z_s) N(\rho_{\text{lim}}, z_s, \mu), \qquad (\text{A27})$$

which can be done with similar steps as above, yielding:

$$M(\rho_{\text{lim}}, z_s) = \frac{N_\rho}{\gamma-1} (\rho_1^{-\gamma+1} - \rho_2^{-\gamma+1}) \bar{\mu} + \frac{N_\rho}{\gamma-1} \left[ \rho_2^{-\gamma+1} \int_1^{(\rho_{\text{lim}}/\rho_2)^2} d\mu \, \mu \, p(\mu, z_s) - \rho_1^{-\gamma+1} \int_1^{(\rho_{\text{lim}}/\rho_1)^2} d\mu \, \mu \, p(\mu, z_s) + \int_{(\rho_{\text{lim}}/\rho_2)^2}^{(\rho_{\text{lim}}/\rho_1)^2} d\mu \, \mu \, p(\mu, z_s) \left(\frac{\sqrt{\bar{\mu}}}{\rho_{\text{lim}}}\right)^{\gamma-1} \right], \qquad (\text{A28})$$

where $\bar{\mu}$ depends on the source redshift and is defined as

$$\bar{\mu}(z_s) = \int_{1}^{\infty} d\mu \, \mu p(\mu, z_s). \qquad (\text{A29})$$

The average magnification is now given by the ratio of (A28) over (A24):

$$\langle \mu \rangle (z_s) = M(\rho_{\text{lim}}, z_s) / N(\rho_{\text{lim}}, z_s). \qquad (\text{A30})$$

The computation of the average magnification presented in this section relies on the assumption that we can replace the mass dependence of the ISCO frequency appearing in the expression of the SNR (A18) by an average. From a comparison with the numerical models of Section 4, it turns out that this approximation is well justified for low redshift sources ($z_s < 2$ when considering the sensitivity curve of the LIGO/Virgo detectors). The comparison is shown in Figure A2. As already explained, the average magnification does not depend on the scaling of the number of events with redshift, i.e., on $N_\rho(z_s)$ in Equation (A23). To compare with a realistic population model, in Figure A2, we reconstructed this dependence fitting the numerical model, i.e., $N_\rho(z_s) \propto \rho^{-3}$. For high redshift sources, a simple analytic derivation of the average magnification is not possible, but the relation (A18) can be numerically inverted to find the chirp mass and the average magnification can be numerically computed as a function of redshift.

## Appendix E. Analytical Results for the LIGO Detector(s)

Our analytic model for the distribution of sources has four free parameters: ($\gamma, M_{\text{min}}, M_{\text{max}}$) defining the distribution of chirp masses and $\rho_{\text{lim}}$ (or equivalently $S_n$) which gives the minimum SNR for detection, and depends on the detector network chosen. We focus now on the case of the LIGO detectors and we compute the average magnification for different choices of the parameters ($\gamma, M_{\text{min}}, M_{\text{max}}$). For definiteness, we consider the noise curve of the LIGO-Livingstone detector for the O2 run. The sensitivity for O3 is slightly better. In Figure A3, we plot the spectral noise of the LIGO-Livingstone detector that we find to be well described by the fitting function $S_{\text{fit}}(f)$ given by

$$S_{\text{fit}}(f) = A_s \left[ \left( \frac{f}{30 \text{ Hz}} \right)^{\alpha_s} + 1 \right] \exp \left[ \beta_s \left( \log \left( \frac{f}{150 \text{ Hz}} \right) \right)^2 \right], \qquad (\text{A31})$$

with $\log_{10} A_s = -46.3$, $\alpha_s = -7.8$ and $\beta_s = 0.4$.
<figure><img src="image_16.png" /><figcaption>**Figure A3.** Spectral noise curve for the LIGO-Livingstone detector [71]. The red curve is the fit given in Equation (A31).</figcaption></figure>

We analyze the dependence of the average amplification on the details of the astro-
physical model. In particular, we want to explore the dependence on the distribution
of the component masses and on the upper cut-off in the distribution. We assume that
the primary BH mass is distributed as $p(m_1) \propto m_1^{-\alpha}$ in the range [5, X] $M_{\odot}$ where both $\alpha$
and X are parameters that we vary. We consider two cases: $\alpha = 2.35$ and $\alpha = 1.5$. The
corresponding normalized distribution of chirp masses is polynomial in the chirp mass, i.e.,
$p(M_c) \propto M_c^{-\beta_M}$ with $\beta_M = 3.2$ and $\beta_M = 2.1$, respectively. As explained in the previous
section, for a sufficiently small redshift, the scaling of the number density of sources per unit
of signal to noise is given by (A23) with $\gamma = -1/5(1 - 6\beta_M)$; hence, $\gamma = 3.6$ and 2.3 for the
two values of $\alpha$, respectively. In Figure A4, we show results for the average magnification
as a function of redshift for the two models considered and for two different values of the
upper cut-off in the mass distribution. The result for the average amplification in Figure A4
is in excellent agreement with the one we found using a realistic model for source formation
and evolution for low redshift sources (as can be seen in Figure 8, left panel).

<figure><img src="image_17.png" /><figcaption>**Figure A4.** (Left) Average magnification for $\alpha = 2.35$ ($\gamma = 3.6$) for two different upper cut-offs in the distribution of the primary mass. (Right) Average magnification for $\alpha = 1.5$ ($\gamma = 2.3$) for two different upper cut-offs in the distribution of the primary mass.</figcaption></figure>

References

1. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration] Observation of Gravitational Waves from a Binary Black Hole Merger. Phys. Rev. Lett. **2016**, *116*, 061102. [CrossRef] [PubMed]

2. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration] GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence. Phys. Rev. Lett. **2016**, *116*, 241103. [CrossRef] [PubMed]

3. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration] Binary Black Hole Mergers in the first Advanced LIGO Observing Run. Phys. Rev. **2016**, *X6*, 041015.

4. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration] GW170608: Observation of a 19-solar-mass Binary Black Hole Coalescence. Astrophys. J. **2017**, *851*, L35. [CrossRef]

5. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration] GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence. Phys. Rev. Lett. **2017**, *119*, 141101. [CrossRef] [PubMed]

6. Abbott, B.P.; et al. [LIGO Scientific and Virgo Collaboration] GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2. Phys. Rev. Lett. **2017**, *118*, 221101.
7. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration] GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral. *Phys. Rev. Lett.* **2017**, *119*, 161101. [CrossRef] [PubMed]

8. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration] GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs. *Phys. Rev.* **2019**, *X9*, 031040. [CrossRef]

9. Abbott, R.; et al. [LIGO Scientific Collaboration and Virgo Collaboration] GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run. *Phys. Rev. X* **2021**, *11*, 021053. [CrossRef]

10. Abbott, R.; et al. [LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration] Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences. *Astrophys. J. Lett.* **2021**, *915*, L5. [CrossRef]

11. Abbott, R.; et al. [LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration] GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run. *arXiv* **2021**, arXiv:2111.03606.

12. Oguri, M. Effect of gravitational lensing on the distribution of gravitational waves from distant binary black hole mergers. *Mon. Not. R. Astron. Soc.* **2018**, *480*, 3842–3855. [CrossRef]

13. Wang, Y.; Stebbins, A.; Turner, E.L. Gravitational Lensing of Gravitational Waves from Merging Neutron Star Binaries. *Phys. Rev. Lett.* **1996**, *77*, 2875–2878. [CrossRef] [PubMed]

14. Sereno, M.; Sesana, A.; Bleuler, A.; Jetzer, P.; Volonteri, M.; Begelman, M.C. Strong Lensing of Gravitational Waves as Seen by LISA. *Phys. Rev. Lett.* **2010**, *105*, 251101. [CrossRef] [PubMed]

15. Ng, K.K.Y.; Wong, K.W.K.; Broadhurst, T.; Li, T.G.F. Precise LIGO Lensing Rate Predictions for Binary Black Holes. *Phys. Rev.* **2018**, D97, 023012, [CrossRef]

16. Broadhurst, T.; Diego, J.M.; Smoot, G. Reinterpreting Low Frequency LIGO/Virgo Events as Magnified Stellar-Mass Black Holes at Cosmological Distances. *arXiv* **2018**, arXiv:1802.05273.

17. Broadhurst, T.; Diego, J.M.; Smoot, G.F. Twin LIGO/Virgo Detections of a Viable Gravitationally-Lensed Black Hole Merger. *arXiv* **2019**, arXiv:1901.03190.

18. Haris, K.; Mehta, A.K.; Kumar, S.; Venumadhav, T.; Ajith, P. Identifying strongly lensed gravitational wave signals from binary black hole mergers. *arXiv* **2018**, arXiv:1807.07062.

19. Cao, Z.; Li, L.F.; Wang, Y. Gravitational lensing effects on parameter estimation in gravitational wave detection with advanced detectors. *Phys. Rev. D* **2014**, *90*, 062003. [CrossRef]

20. Dai, L.; Venumadhav, T. On the waveforms of gravitationally lensed gravitational waves. *arXiv* **2017**, arXiv:1702.04724.

21. Dai, L.; Li, S.S.; Zackay, B.; Mao, S.; Lu, Y. Detecting Lensing-Induced Diffraction in Astrophysical Gravitational Waves. *Phys. Rev. D* **2018**, *98*, 104029, [CrossRef]

22. Dai, L.; Venumadhav, T.; Sigurdson, K. Effect of lensing magnification on the apparent distribution of black hole mergers. *Phys. Rev. D* **2017**, *95*, 044011, [CrossRef]

23. Jung, S.; Shin, C.S. Gravitational-Wave Fringes at LIGO: Detecting Compact Dark Matter by Gravitational Lensing. *Phys. Rev. Lett.* **2019**, *122*, 041103, [CrossRef] [PubMed]

24. Lai, K.H.; Hannuksela, O.A.; Herrera-Martín, A.; Diego, J.M.; Broadhurst, T.; Li, T.G.F. Discovering intermediate-mass black hole lenses through gravitational wave lensing. *Phys. Rev. D* **2018**, *98*, 083005, [CrossRef]

25. Smith, G.P.; Jauzac, M.; Veitch, J.; Farr, W.M.; Massey, R.; Richard, J. What if LIGO's gravitational wave detections are strongly lensed by massive galaxy clusters? *Mon. Not. R. Astron. Soc.* **2018**, *475*, 3823–3828. [CrossRef]

26. Shan, X.; Wei, C.; Hu, B. Lensing magnification: gravitational waves from coalescing stellar-mass binary black holes. *arXiv* **2020**, arXiv:2012.08381.

27. Diego, J.M.; Broadhurst, T.; Smoot, G. Evidence for lensing of gravitational waves from LIGO-Virgo. *arXiv* **2021**, arXiv:2106.06545.

28. Pagano, G.; Hannuksela, O.A.; Li, T.G.F. lensingGW: A Python package for lensing of gravitational waves. *Astron. Astrophys.* **2020**, *643*, A167, [CrossRef]

29. Mukherjee, S.; Wandelt, B.D.; Silk, J. Probing the theory of gravity with gravitational lensing of gravitational waves and galaxy surveys. *Mon. Not. R. Astron. Soc.* **2020**, *494*, 1956–1970. [CrossRef]

30. Chung, K.W.; Li, T.G.F. Lensing of Gravitational Waves as a Novel Probe of Graviton Mass. *arXiv* **2021**, arXiv:2106.09630.

31. Urrutia, J.; Vaskonen, V. Lensing of gravitational waves as a probe of compact dark matter. *arXiv* **2021**, arXiv:2109.03213.

32. Cremonese, P.; Mota, D.F.; Salzano, V. Characteristic features of gravitational wave lensing as probe of lens mass model. *arXiv* **2021**, arXiv:2111.01163.

33. Hannuksela, O.A.; Haris, K.; Ng, K.K.Y.; Kumar, S.; Mehta, A.K.; Keitel, D.; Li, T.G.F.; Ajith, P. Search for gravitational lensing signatures in LIGO-Virgo binary black hole events. *Astrophys. J.* **2019**, *874*, L2.

34. Oguri, M. Strong gravitational lensing of explosive transients. *Rept. Prog. Phys.* **2019**, *82*, 126901, [CrossRef] [PubMed]

35. Diego, J.M. Constraining the abundance of primordial black holes with gravitational lensing of gravitational waves at LIGO frequencies. *arXiv* **2019**, arXiv:1911.05736.

36. Broadhurst, T.; Diego, J.M.; Smoot, G.F. Interpreting LIGO/Virgo “Mass-Gap” events as lensed Neutron Star-Black Hole binaries. *arXiv* **2020**, arXiv:2006.13219.

37. Gavazzi, R.; Treu, T.; Rhodes, J.D.; Koopmans, L.V.; Bolton, A.S.; Burles, S.; Massey, R.; Moustakas, L.A. The Sloan Lens ACS Survey. 4. The mass density profile of early-type galaxies out to 100 effective radii. *Astrophys. J.* **2007**, *667*, 176–190. [CrossRef]
38. Robertson, A.; Smith, G.P.; Massey, R.; Eke, V.; Jauzac, M.; Bianconi, M.; Ryczanowski, D. What does strong gravitational lensing? The mass and redshift distribution of high-magnification lenses. *arXiv* 2020, arXiv:2002.01479.

39. Belczynski, K.; Holz, D.E.; Bulik, T.; O'Shaughnessy, R. The first gravitational-wave source from the isolated evolution of two 40–100 Msun stars. *Nature* **2016**, *534*, 512. [CrossRef]

40. Eldridge, J.J.; Stanway, E.R. BPASS predictions for Binary Black-Hole Mergers. *Mon. Not. R. Astron. Soc.* **2016**, *462*, 3302–3313. [CrossRef]

41. Wysocki, D.; Gerosa, D.; O'Shaughnessy, R.; Belczynski, K.; Gladysz, W.; Berti, E.; Kesden, M.; Holz, D.E. Explaining LIGO's observations via isolated binary evolution with natal kicks. *Phys. Rev. D* **2018**, *97*, 043014. [CrossRef]

42. Dvorkin, I.; Uzan, J.P.; Vangioni, E.; Silk, J. Exploring stellar evolution with gravitational-wave observations. *Mon. Not. R. Astron. Soc.* **2018**, *479*, 121–129. [CrossRef]

43. Kruckow, M.U.; Tauris, T.M.; Langer, N.; Kramer, M.; Izzard, R.G. Progenitors of gravitational wave mergers: Binary evolution with the stellar grid-based code ComBinE. *Mon. Not. R. Astron. Soc.* **2018**, *481*, 1908–1949. [CrossRef]

44. Spera, M.; Mapelli, M.; Giacobbo, N.; Trani, A.A.; Bressan, A.; Costa, G. Merging black hole binaries with the SEVN code. *Mon. Not. R. Astron. Soc.* **2019**, *485*, 889–907. [CrossRef]

45. Mapelli, M.; Giacobbo, N.; Santoliquido, F.; Artale, M.C. The properties of merging black holes and neutron stars across cosmic time. *Mon. Not. R. Astron. Soc.* **2019**, *487*, 2–13. [CrossRef]

46. Rodriguez, C.L.; Zevin, M.; Amaro-Seoane, P.; Chatterjee, S.; Kremer, K.; Rasio, F.A.; Ye, C.S. Black holes: The next generation—Repeated mergers in dense star clusters and their gravitational-wave properties. *Phys. Rev. D* **2019**, *100*, 043027. [CrossRef]

47. Stevenson, S.; Sampson, M.; Powell, J.; Vigna-Gómez, A.; Neijssel, C.J.; Szécsi, D.; Mandel, I. The impact of pair-instability mass loss on the binary black hole mass distribution. *Astrophys. J.* **2019**, *882*, 121. [CrossRef]

48. Fragione, G.; Leigh, N.; Perna, R. Black hole and neutron star mergers in Galactic Nuclei: the role of triples. *Mon. Not. R. Astron. Soc.* **2019**, *488*, 2825–2835. [CrossRef]

49. Woosley, S.E.; Heger, A. The Pair-Instability Mass Gap for Black Holes. *Astrophys. J. Lett.* **2021**, *912*, L31. [CrossRef]

50. Abbott, B.P. et al. [LIGO Scientific Collaboration and Virgo Collaboration]. Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo. *Astrophys. J.* **2019**, *882*, L24. [CrossRef]

51. Turner, E.L.; Ostriker, J.P.; Gott, J. Richard, I. The Statistics of gravitational lenses: The Distributions of image angular separations and lens redshifts. *Astrophys. J.* **1984**, *284*, 1–22. [CrossRef]

52. Fukugita, M.; Futamase, T.; Kasai, M.; Turner, E.L. Statistical Properties of Gravitational Lenses with a Nonzero Cosmological Constant. *Astrophys. J.* **1992**, *393*, 3. [CrossRef]

53. Oguri, M.; Marshall, P.J. Gravitationally lensed quasars and supernovae in future wide-field optical imaging surveys. *Mon. Not. R. Astron. Soc.* **2010**, *405*, 2579–2593. [CrossRef]

54. Hilbert, S.; White, S.D.; Hartlap, J.; Schneider, P. Strong lensing optical depths in a LambdaCDM universe. *Mon. Not. R. Astron. Soc.* **2007**, *382*, 121–132. [CrossRef]

55. Hilbert, S.; White, S.D.; Hartlap, J.; Schneider, P. Strong lensing optical depths in a LCDM Universe. 2. The influence of the stellar mass in galaxies. *Mon. Not. R. Astron. Soc.* **2008**, *386*, 1845–1854. [CrossRef]

56. Schneider, P.; Ehlers, J.; Falco, E.E. *Gravitational Lenses*; Springer: Berlin/Heidelberg, Germany, 1992. [CrossRef]

57. Finn, L.S.; Chernoff, D.F. Observing binary inspiral in gravitational radiation: One interferometer. *Phys. Rev.* **1993**, *D47*, 2198–2219. [CrossRef] [PubMed]

58. Maggiore, M. *Gravitational Waves. Vol. 1: Theory and Experiments*; Oxford Master Series in Physics; Oxford University Press: Oxford, UK, 2007.

59. Torrey, P.; Wellons, S.; Machado, F.; Griffen, B.; Nelson, D.; Rodriguez-Gomez, V.; McKinnon, R.; Pillepich, A.; Ma, C.P.; Vogelsberger, M.; Springel, V.; Hernquist, L. An analysis of the evolving comoving number density of galaxies in hydrodynamical simulations. *Mon. Not. R. Astron. Soc.* **2015**, *454*, 2770–2786. [CrossRef]

60. Vogelsberger, M.; Genel, S.; Springel, V.; Torrey, P.; Sijacki, D.; Xu, D.; Snyder, G.F.; Nelson, D.; Hernquist, L. Introducing the Illustris Project: Simulating the coevolution of dark and visible matter in the Universe. *Mon. Not. R. Astron. Soc.* **2014**, *444*, 1518–1547. [CrossRef]

61. Vogelsberger, M.; Genel, S.; Springel, V.; Torrey, P.; Sijacki, D.; Xu, D.; Snyder, G.F.; Bird, S.; Nelson, D.; Hernquist, L. Properties of galaxies reproduced by a hydrodynamic simulation. *Nature* **2014**, *509*, 177–182. [CrossRef]

62. Genel, S.; Vogelsberger, M.; Springel, V.; Sijacki, D.; Nelson, D.; Snyder, G.; Rodriguez-Gomez, V.; Torrey, P.; Hernquist, L. Introducing the Illustris Project: the evolution of galaxy populations across cosmic time. *Mon. Not. R. Astron. Soc.* **2014**, *445*, 175–200. [CrossRef]

63. Torrey, P. (Department of Astronomy, University of Florida, Gainesville, FL, USA). Personal communication, 2019.

64. Takahashi, R.; Sato, M.; Nishimichi, T.; Taruya, A.; Oguri, M. Revising the Halofit Model for the Nonlinear Matter Power Spectrum. *Astrophys. J.* **2012**, *761*, 152. [CrossRef]

65. Vangioni, E.; Olive, K.A.; Prestegard, T.; Silk, J.; Petitjean, P.; Mandic, V. The Impact of Star Formation and Gamma-Ray Burst Rates at High Redshift on Cosmic Chemical Evolution and Reionization. *Mon. Not. R. Astron. Soc.* **2015**, *447*, 2575. [CrossRef]

66. Ajith, P.; Hannam, M.; Husa, S.; Chen, Y.; Brügmann, B.; Dorband, N.; Müller, D.; Ohme, F.; Pollney, D.; Reisswig, C.; et al. *Inspiral-merger-ringdown waveforms for black-hole binaries with non-precessing spins*. *Phys. Rev. Lett.* **2011**, *106*, 241101. [CrossRef] [PubMed]
67. Collaboration, L.S. *H1 Calibrated Sensitivity Spectra June 10 2017*; Technical Report; LIGO Scientific Collaboration: Cascina, Italy, 2017.

68. Collaboration, L.S. *L1 Calibrated Sensitivity Spectra Aug 06 2017*; Technical Report; LIGO Scientific Collaboration: Cascina, Italy, 2017.

69. Collaboration, L.S. *Unofficial Sensitivity Curves for Various Detector Configurations Used in ISWP 2016*; Technical Report; LIGO Scientific Collaboration: Cascina, Italy, 2016.

70. Cusin, G.; Tamanini, N. Characterization of lensing selection effects for LISA massive black hole binary mergers. *Mon. Not. R. Astron. Soc.* **2021**, *504*, 3610–3618. [CrossRef]

71. LIGO-Virgo. 2019. Available online: https://dcc.ligo.org/LIGO-P1900011/public (accessed on 16 November 2019).