Lensing Magnification Seen by Gravitational Wave Detectors

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 observedsource 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$).

1. Introduction

The direct detection of gravitational waves (GWs) from inspiraling binary black holes and neutron stars by the LIGO and Virgo experiments [1,2,3,4,5,6,7,8,9,10,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,14,15,16,17,18,19,20,21,22,23,24,25,26,27,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,30,31,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,34,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 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,2,3,4,5,6,8,9,10,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 ∼$25M_{\odot }$) and in some cases, BBHs detected in GW have masses well above $60M_{\odot }$ [6,8]. Black hole (BH) masses above ∼$25M_{\odot }$ are expected in low-metallicity environments and are indeed predicted by various stellar evolution models (e.g., [39,40,41,42,43,44,45,46,47,48]). However, standard stellar evolution models struggle to explain masses in the so-called ’pair-instability mass gap’ in the range between ∼$60M_{\odot }$ and ∼$130M_{\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,52,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 $\Lambda$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 $\mathbf{\theta }$ and $\mathbf{\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\mathbf{\beta }$ and a vector $\mathbf{\xi }$ in the lens plane defined as $\mathbf{\xi }=\mathbf{\theta }{d}_{\ell }$. This fixes the geometry. Here, $d_{\ell }=d\left(z_{\ell }\right),d_{\ell s}=d(z_{\ell },z_s)$ and $d_s=d\left(z_s\right)$ 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. Glossary of the main quantities which we introduce and use in this paper.

Symbol	Definition
$z_s$	(cosmological) source redshift
$z_{\ell }$	(cosmological) lens redshift
$D_s$	luminosity distance to the source
$D_{\ell }$	luminosity distance to the lens
$d_s$	angular diameter distance to the source
$d_{\ell }$	angular diameter distance to the lens
$d_{\ell s}$	angular diameter distance from the lens to the source
$z_{\mathrm{obs}}$	observed source redshift (inferred from standard $\Lambda$CDM assuming no magnification)
${\sigma }_v$	galaxy velocity dispersion
$\mu$	magnification
$\Delta t$	time-delay between multiple images
$\sigma$	lensing cross-section
$\tau (\mu ,z_s)$	optical depth for magnification by at least $\mu$ for source at $z_s$
$p(\mu ,z_s)$	probability density of magnification for source at $z_s$
$P(>\mu ,z_s)$	probability of magnification by at least $\mu$ for source at $z_s$
$F_{\lim }$	limiting flux of a given observatory
${\rho }_{\lim }$	limiting signal-to-noise ratio (corresponding to $F_{\lim }$)
$d{\mathcal{N}}_{\mathrm{obs}}(F_{\lim },z_s)$	# sources observed with liming flux $F_{\lim }$ in a bin around $z_s$
$d{\mathcal{N}}_{\mathrm{obs}}(F_{\lim },D_{\mathrm{obs}})$	# sources observed with liming flux $F_{\lim }$ in a bin around $D_{\mathrm{obs}}$
$d\mathcal{N}(\mu ,F_{\lim },z_s)$	# sources that we see with liming flux $F_{\lim }$ in a bin around $z_s$, if the magnification is $\mu$
$d\mathcal{N}(\mu ,F_{\lim },D_{\mathrm{obs}})$	# sources that we see with liming flux $F_{\lim }$ in a bin around $D_{\mathrm{obs}}$, if the magnification is $\mu$
${\mathcal{P}}_z\left(\mu \right)$	probability density that an observed source at (cosmological) redshift z is magnified by $\mu$
${\mathcal{P}}_{\mathrm{obs}}\left(\mu \right)$	probability density that an observed source at observed redshift $z_{\mathrm{obs}}$ is magnified by $\mu$

.

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 \left(r\right)=\frac{{\sigma }_v^2}{2\pi G{r}^2},$$

(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 \left(\theta \right)=\frac{{\sigma }_v^2}{2G}\frac{1}{d_{\ell }\theta }.$$

(2)

We introduce ${\Sigma }_c={\left(4\pi G\right)}^{-1}{d}_s/\left(d_{\ell s}{d}_{\ell }\right)$ so that the convergence $\kappa$ is given by

$$\kappa \left(\theta \right)=\frac{\Sigma \left(\theta \right)}{{\Sigma }_c}=2\pi {\sigma }_v^2\frac{d_{\ell s}}{d_s}\frac{1}{\theta }.$$

(3)

An SIS has a constant deflection angle given by

$$\alpha \left(\theta \right)=4\pi {\sigma }_v^2\frac{d_{\ell s}}{d_s}=2\theta \kappa \left(\theta \right)\equiv {\alpha }_0.$$

(4)

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

$$\begin{array}{cc}\hfill & A_{ij}={\delta }_{ij}\left(1-\frac{1}{\left|\mathbf{x}\right|}\right)+\frac{x_i{x}_j}{{\left|\mathbf{x}\right|}^3},\hfill \\ \hfill & \det A=1-\frac{1}{\left|\mathbf{x}\right|},\mu =\frac{1}{\left|\det A\right|}=\frac{\left|\mathbf{x}\right|}{|1-|\mathbf{x}\left|\right|},\hfill \end{array}$$

(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 \left(y\right)=\left\{\begin{array}{cc}\mu \left({\mathbf{x}}_1\right)+\mu \left({\mathbf{x}}_2\right)=\frac{y+1}{y}+\frac{1-y}{y}={\mu }_1+{\mu }_2=\frac{2}{y},\hfill & y\le 1,\hfill \\ \frac{y+1}{y}=1+\frac{1}{y},\hfill & y\ge 1.\hfill \end{array}\right.$$

(6)

2.2. Optical Depth and Lensing Probability

The impact parameter of the source (in the source plane) is given by $\left|\mathbf{b}\right|=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 \left(y\right)$. 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 {\left(y{\alpha }_0{d}_s\right)}^2=\frac{\pi {\left(4\pi \right)}^2{\sigma }_v^4{d}_{\ell s}^2}{{({\mu }_1-1)}^2}.$$

(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}d{z}_{\ell }\frac{dr}{d{z}_{\ell }}\begin{array}{c}\hfill {\left(\frac{d_{\ell }}{d_s}\right)}^2\end{array}\int \sigma (\mu ,z_{\ell },z_s,{\sigma }_v)n({\sigma }_v,z_{\ell })d{\sigma }_v,$$

(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:

$$\begin{array}{cc}\hfill & \tau (\mu ,z_s)=\hfill \\ \hfill & =\frac{\pi {\left(4\pi \right)}^2}{{(\mu -1)}^2}{\int }_0^{z_s}dz\frac{{\chi }^2(z,z_s){\chi }^2\left(z\right)}{{\chi }^2\left(z_s\right){(1+z)}^{3}H\left(z\right)}{\int }_0^{\infty }d{\sigma }_v{\sigma }_v^{4}n({\sigma }_v,z),\hfill \end{array}$$

(9)

where $\chi (z_1,z_2)$ denotes the co-moving distance from the redshift $z_1$ to $z_2$ and $\chi \left(z\right)=\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 ,$$

(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)).$$

(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{{\left|h\left(f\right)\right|}^2}{S_n\left(f\right)}df=\int {\rho }^2\left(f\right)df.$$

(12)

Here, $S_n\left(f\right)$ 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}}_{\mathrm{obs}}(z,{\rho }_{lim})={\int }_0^{\infty }d\mu p(\mu ,z)\mathcal{N}(\mu ,z),$$

(13)

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

$$d\mathcal{N}(\mu ,z)\equiv dz{\int }_{{\rho }_{lim}/\sqrt{\mu }}^{\infty }N(\rho ,z)d\rho ,$$

(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 \left(z\right)=\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)}.$$

(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},$$

(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 $\mathrm{Var}\left(\mu \right)\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_{\mathrm{obs}}(z_s,\mu )\equiv \frac{D_s\left(z_s\right)}{\sqrt{\mu }}.$$

(17)

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

$$\begin{array}{cc}\hfill & d{\mathcal{N}}_{\mathrm{obs}}(D_{\mathrm{obs}},{\rho }_{lim})={\int }_0^{\infty }d\mu p(\mu ,z_s(D_{\mathrm{obs}},\mu ))d\mathcal{N}(\mu ,D_{\mathrm{obs}}),\hfill \end{array}$$

(18)

where:

$$d\mathcal{N}(\mu ,D_{\mathrm{obs}})=d{D}_{\mathrm{obs}}\frac{\sqrt{\mu }}{D_s^{\prime }\left(z_s\right)}{\int }_{{\rho }_{lim}/\sqrt{\mu }}^{\infty }N(\rho ,z_s(D_{\mathrm{obs}},\mu ))d\rho ,$$

(19)

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

$$\langle \mu \rangle \left(D_{\mathrm{obs}}\right)=\frac{{\int }_0^{\infty }d\mu \mu p(\mu ,D_{\mathrm{obs}})d\mathcal{N}(\mu ,D_{\mathrm{obs}})/d{D}_{\mathrm{obs}}}{{\int }_0^{\infty }d\mu p(\mu ,D_{\mathrm{obs}})d\mathcal{N}(\mu ,D_{\mathrm{obs}})/d{D}_{\mathrm{obs}}},$$

(20)

which of course depends on the threshold ${\rho }_{\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}}_{\mathrm{obs}}(\mu ,D_{\mathrm{obs}})\equiv \mathcal{C}p(\mu ,D_{\mathrm{obs}})\frac{d\mathcal{N}(\mu ,D_{\mathrm{obs}})}{d{D}_{\mathrm{obs}}},$$

(21)

where $\mathcal{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_{\mathrm{obs}})/d{D}_{\mathrm{obs}}\propto \sqrt{\mu }$ for $\mu \gg 1$ as both the lower integration bound and $D_s^{\prime }$ tend towards constant asymptotic values for large magnification $\mu$. The probability density $p(\mu ,D_{\mathrm{obs}})$ in (21) also depends on $\mu$ through the observed luminosity as we defined $p(\mu ,D_{\mathrm{obs}})\equiv p(\mu ,z_s\left(\sqrt{\mu }{D}_{\mathrm{obs}}\right))$. Recalling that $p(\mu ,z_s)\propto {\mu }^{-3}$, it follows that $p(\mu ,D_{\mathrm{obs}})$ will have a milder decay with magnification. Indeed, for a fixed $D_{\mathrm{obs}}$ and $\mu >1$, $p(\mu ,D_{\mathrm{obs}})\equiv p(\mu ,z_s\left(\sqrt{\mu }{D}_{\mathrm{obs}}\right))$ also receives contributions from redshifts higher than $z_s\left(D_{\mathrm{obs}}\right)$ and therefore the lensing probability is higher. The exact scaling of $p(\mu ,D_{\mathrm{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_{\mathrm{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 ${\mathcal{P}}_{\mathrm{obs}}\left(\mu \right)\propto {\mu }^{-\alpha +1/2}$ and its integral from $\mu$ to infinity behaves as ${\mathcal{P}}_{\mathrm{obs}}(>\mu )\propto {\mu }^{-\alpha +3/2}$ for large values of $\mu$. It follows that, for any fixed value of $D_{\mathrm{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, 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 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_{\mathrm{obs}}$. We will get back to this point in Section 4 and Section 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,61,62] is a cosmological hydrodynamical simulation run in a periodic box of size $L=106.5$ Mpc, dark matter mass resolution of $m_{\mathrm{DM}}=6.3\times {10}^6{M}_{\odot }$ and baryon mass resolution of $m_{\mathrm{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]:

$$\begin{array}{cc}\hfill & {log}_{10}N(>{\sigma }_v,z)=A\left(z\right)+\alpha \left(z\right)\left({log}_{10}{\sigma }_v-\gamma \left(z\right)\right)+\hfill \\ \hfill & +\beta \left(z\right){\left({log}_{10}{\sigma }_v-\gamma \left(z\right)\right)}^2-{\left({\sigma }_v\times {10}^{-\gamma \left(z\right)}\right)}^{1/ln\left(10\right)},\hfill \end{array}$$

(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^{\mathrm{com}}({\sigma }_v,z)\equiv -\frac{dN}{d{\sigma }_v}(>{\sigma }_v,z).$$

(23)

The result for different redshifts is shown in Figure 2.

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_{\mathrm{obs}}=D\left(z_s\right)/\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.

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 \left(t\left(z\right)\right)=A\frac{a{e}^{b(z-z_m)}}{a-b+b·e^{a(z-z_m)}},$$

(24)

with $A=0.24M_{\odot }$/yr/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 $ϵ$ 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_{\mathrm{delay}}$ between the formation of the progenitor star and the coalescence of the BBH, where the latter is distributed as

$$P_d\left(t_{\mathrm{delay}}\right)\propto t_{\mathrm{delay}}^{-1},$$

(25)

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

$$P\left(m_1\right)\propto m_1^{-\alpha },$$

(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\left(m_2\right)=const.$ in this range and $p(m_1,m_2)=p\left(m_1\right)p\left(m_2\right)$ [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{array}{cc}\hfill & \frac{\mathrm{d}R}{\mathrm{d}{m}_1\mathrm{d}{m}_2\mathrm{d}z}=\hfill \\ \hfill & =ϵT_{\mathrm{obs}}\int \Psi (t\left(z\right)-t_{\mathrm{delay}})p(m_1,m_2)P_d\left(t_{\mathrm{delay}}\right)\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}{t}_{\mathrm{delay}}}.\hfill \end{array}$$

(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 }_{lim}=10$. The normalization parameter $ϵ$ 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_{\mathrm{obs}}={\int }_{\rho (m_1,m_2,z)\ge 10}\frac{\mathrm{d}R}{\mathrm{d}{m}_1\mathrm{d}{m}_2\mathrm{d}z}\mathrm{d}{m}_1\mathrm{d}{m}_2\mathrm{d}z$$

(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:

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_{\mathrm{obs}}$ as a function of magnification. We plot this quantity for three different values of $D_{\mathrm{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_{\mathrm{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_{\mathrm{obs}}=3000$ Mpc than for $D_{\mathrm{obs}}=500$ Mpc, the probability of seeing such objects is higher for $D_{\mathrm{obs}}=3000$ Mpc than for $D_{\mathrm{obs}}=500$ Mpc. This is due to the fact that $p(\mu ,D_{\mathrm{obs}})$ is monotonically increases with $D_{\mathrm{obs}}$, as can be seen in Figure 3. For the product of $p(\mu ,D_{\mathrm{obs}})d\mathcal{N}/d{D}_{\mathrm{obs}}$, the increase in p wins over the decrease in $\mathcal{N}$ for $D_{\mathrm{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, $d\mathcal{N}/dz$ tends towards a constant at high magnification since we see all events (black curve, left panel). The product $p(z,\mu )d\mathcal{N}/dz$ nearly tends towards a constant at high redshift which indicates that $d\mathcal{N}/dz$ grows similarly to ${\mu }^3$ at high redshift (blue dashed curve, right panel).

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_{\mathrm{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 }_{\max }\left(z_s\right)>{\rho }_{\lim }$ is given by

$${\mu }_{\min }\left(z_s\right)={\left(\frac{{\rho }_{\lim }}{{\rho }_{\max }}\right)}^2,$$

(29)

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

$$D_{\mathrm{obs}}^{\max }\left(z_s\right)=\frac{D\left(z_s\right)}{\sqrt{{\mu }_{\min }}}=\frac{D\left(z_s\right){\rho }_{\max }\left(z_s\right)}{{\rho }_{\lim }}.$$

(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 ${\mathcal{P}}_{\mathrm{obs}}(>\mu )$ that an object that LIGO/Virgo observes from a distance $D_{\mathrm{obs}}$ was magnified more than $\mu$ is presented in Figure 9 for different values of $D_{\mathrm{obs}}$ as a function of $\mu$.

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 ${\mathcal{M}}_z=(1+z)\mathcal{M}$, where $\mathcal{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 $\mathcal{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_{\mathrm{obs}})\equiv p(\mu ,z_s\left(\sqrt{\mu }{D}_{\mathrm{obs}}\right))$, has an even slower decay with magnification, as can be seen in the right panel of Figure 3. Indeed, for a fixed $D_{\mathrm{obs}}$ and $\mu >1$, $p(\mu ,D_{\mathrm{obs}})\equiv p(\mu ,z_s\left(\sqrt{\mu }{D}_{\mathrm{obs}}\right))$ takes contributions from a higher redshift than $z_s\left(D_{\mathrm{obs}}\right)$ and is therefore larger. The exact scaling of $p(\mu ,D_{\mathrm{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_{\mathrm{obs}})\propto {\mu }^{-\alpha }$. The value of $\alpha$ increases from 2 to 3 as the observed luminosity distance varies from ∼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 ${\mathcal{P}}_{\mathrm{obs}}\left(\mu \right)$ 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 Figure 6 and Figure 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 ${\mathcal{P}}_{\mathrm{obs}}\left(\mu \right)$ for a given observed distance scales as ${\mathcal{P}}_{\mathrm{obs}}\left(\mu \right)\propto {\mu }^{-\alpha +1/2}$ and its integral from $\mu$ to infinity behaves as ${\mathcal{P}}_{\mathrm{obs}}(>\mu )\propto {\mu }^{-\alpha +3/2}$ for large values of $\mu$ with $\alpha \in [2,3]$. The probability ${\mathcal{P}}_{\mathrm{obs}}(>\mu )$ is presented in Figure 9 for different values of $D_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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 ${\mathcal{P}}_{\mathrm{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 ${\mathcal{P}}_{\mathrm{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 $\mathcal{P}\left(\mu \right)$ and ${\mathcal{P}}_{\mathrm{obs}}\left(\mu \right)$ (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.