Constraining the abundance of primordial black holes with gravitational lensing of gravitational waves at LIGO frequencies

Gravitational waves from binary black holes that are gravitationally lensed can be distorted by small microlenses along the line of sight. Microlenses with masses of a few tens of solar masses, and that are close to a critical curve in the lens plane, can introduce a time delay of a few millisecond. Such time delay would result in distinctive interference patterns in the gravitational wave that can be measured with current experiments such as LIGO/Virgo. We consider the particular case of primordial black holes with masses between 5 and 50 solar masses acting as microlenses. We study the effect of a population of primordial black holes constituting a fraction of the dark matter, and contained in a macrolens (galaxy or cluster), over gravitational waves that are being lensed by the combined effect of the macrolens plus microlenses. We find that at the typical magnifications expected for observed GW events, the fraction of dark matter in the form of compact microlenses, such as primordial black holes, can be constrained to percent level. Similarly, if a small percentage of the dark matter is in the form of microlenses with a few tens of solar masses, at sufficiently large magnification factors, all gravitational waves will show interference effects. These effects could have an impact on the inferred parameters. The effect is more important for macroimages with negative parity, which usually arrive after the macroimages with positive parity.


I. INTRODUCTION
Gravitational waves (GW hereafter) open new opportunities to study the universe. Due to the all-sky sensitivity of an array of GW detectors, an event with an amplitude above the sensitivity of the detector can be observed independently of its location in the sky (although in reality a detector will see the signal modulated by the geometric factor that does depend on the sky location). Rare events can then be observed provided their amplitude is large enough. Strong gravitational lensing of distant background objects can be considered as a rare event since the probability that a distant source (z > 1) is strongly lensed is less than 0.1 %. Observations of transients (like supernovae) that are strongly lensed are extremely rare since detecting transients has traditionally required to be observing that particular region of the sky at the time of the transient and only a small fraction of these transients are expected to be strongly lensed. In the case of GWs, one does not need to be pointing a telescope to a particular region in the sky. Then, if a distant GW is being strongly lensed, and its magnified amplitude surpasses the detector threshold sensitivity, it will be detected independently of its sky position. If the rate of GWs is sufficiently high at z > 1, the number of GWs per redshift interval at large redshifts (that scales to first order as z 3 ) may be large enough to compensate for the small probability of lensing. Then, some of the events that could not be observed at z > 1 because they are too distant, could be promoted beyond the sensitivity threshold of the detector if they are magnified by some factor. Strong lensing of GWs is different than standard * jdiego@ifca.unican.es lensing of galaxies in different ways. Owing to the very low frequency of the GWs, and depending on the mass of the lens, diffraction effects may be important. This will be studied in more detail in the sections below. On the other hand, due to the very small size of the volume emitting GWs, magnification factors of many thousands are possible. This is not true in the case of strongly lensed galaxies, where magnification factors above a few tens are rare. This is due to the fact that the maximum magnification is proportional to √ R −1 , where R is the radius of the object being lensed. Allowing for magnification factors larger than a few hundred has important implications on the detectability of distant objects. GWs that are emitted at a redshift z ≈ 3, if magnified by a factor ≈ 300 can compensate the increase in the luminosity distance and appear as luminous as a similar event at z ≈ 0.3. On the contrary, the galaxy hosting this GW would have its total flux magnified only by a factor of a few tens at most.
A high rate of events at high redshift may require the existence of an exotic candidate for dark matter; the primordial black holes (or PBHs hereafter). PBHs have been hypothesized as a possible candidate for dark matter [1][2][3][4][5]. Rates of events (with masses ≈ 30 M ) as high as 10 4 events per Gpc 3 and year (at z ≈ 1-2) are possible if the fraction of dark matter in the form of PBH is of order 10% [6]. If a fraction of the dark matter is composed of PBHs, they will not only have an impact on the rate of merger events. PBHs can also act as microlenses, and split an incoming GW into multiple images with a time delay between them. If the mass of the microlens is in a certain range, diffraction effects may take place on the GW [7,8]. The GW may interfere with itself if the time delay between the lensed microimages is of the order of the period of the GW. Lensing of GWs by compact astrophysical bodies has been considered in the past, although normally earlier work considers masses above 100 M [9][10][11][12][13]. In [14], the authors consider lensing by PBHs with masses between ∼ 10 and ∼ 10 5 M . In that work, the lensing probability is computed assuming a simple model in which all microlenses are isolated, that is, the halo (galaxy, group or cluster) hosting the microlens is being ignored. As shown by [15], the role of the macrolens can not be ignored since it can modify significantly the time delays, and consequently the detectability of the microlensing event. The calculation of the optical depth of lensing is also affected by the macrolens. This is particularly important for the most distant detected GWs, most of which will be strongly lensed. At large magnification factors, the probability of microlensing grows as the magnification of the macrolens [16]. As discussed in [15], at sufficiently large magnifications, all strongly lensed GWs will be affected by microlensing. This happens when the effective optical depth is ≈ 1 In this work we focus on lensed BBH, although the results presented in this paper are also relevant for other GWs at higher frequencies like BNS or NSBH. However, since these alternative candidates have a smaller chirp mass, the probabilities of observing these events at higher redshifts is smaller. As demonstrated in earlier work, and also shown in the sections below, the probability of strong lensing is sufficiently high only at redshifts above 0.5. BBH are luminous enough that, after being magnified by some factor, can be observed even if they are produced at z > 1.
Under the hypothesis that dark matter is composed of PBHs with masses around 30 solar masses, GWs that are being strongly lensed by massive haloes are useful probes of dark matter. Events that are strongly lensed necessarily must travel through areas with a relatively large surface mass density. As strong lensing events must take place near the critical curves of halos, the surface mass density near critical curves is close to the critical surface mass density for lensing, which typically is in the range of a few thousand solar masses per parsec square. This means that, in projection, a parsec square must contain dozens of PBHs in the intervening macrolens alone (and more when one accounts for additional matter along the line of sight) We adopt a cosmological model with Ω m = 0.3, Λ = 0.7, and h = 0.7. A few terms will be used throughout the paper. We refer to the massive haloes (galaxies, groups and clusters) as macrolenses. The stars, remnants and PBHs in these macrolenses that perturb the macrolens on the small scale are referred to as microlenses. The images produced by the macrolens are referred to as macroimages. In the presence of microlenses in the lens plane, macroimages of small sources can break up into smaller microimages around the positions of the microlenses. The source frame is the reference frame where the GW is originated while the observer frame is Earth. The structure of this paper is as follows. In section II we give a brief introduction to lensing of GWs in the regime of geometric optics. This regime is the appropriate one to address the probability of lensing of a GW by a macrolens. Section III introduces the chirp mass function of GWs. Section IV discusses the probability of lensing by macrolenses. In section V we describe two models for the evolution of the rate of GWs as a function of time. Section VI estimates the expected number of strongly lensed events for the two models assumed. In section VII we discuss the effects that microlenses have on GWs, including wave effects that, at LIGO frequencies, become relevant for microlenses with masses below a few hundred M . In section VIII we describe the simulation used to test the sensitivity of microlensing of GWs to the abundance of PBHs. Finally, in section IX we present the main results of this work and discuss its implications.

II. LENSING OF GRAVITATIONAL WAVES
Lensing of GWs has been studied in detail in earlier work [9][10][11] and more recently in [14,[17][18][19][20][21][22][23]. In general, wave optics is the appropriate regime for studying the lensing of GWs when the Schwarzschild radius of the lens is comparable to the wavelength of the GW [8,10,24]. For the frequencies probed by experiments such as LIGO/Virgo, wave optics is relevant when the lens has a mass smaller than a few thousand solar masses. For lenses with masses above 10 4 solar masses, one can rely on geometric optics. We discuss the regime of wave optics in section VII.This section presents only a brief discussion of lensing of GWs in the regime of geometric optics.
The most basic quantity that can be measured from a GW is the observed (redshifted) chirp mass, where m 1 , and m 2 are masses of the two objects before coalescence (in the source frame), and M chirp is the chirp mass in the source frame. Only GWs that are above the detection threshold of the detector can be observed. This threshold is determined by the signal-to-noise-ratio, ρ, of the GW. To first order, the value of ρ depends on the luminosity distance of the GW, D l (z), and the observed chirp mass, M. More specifically, ρ is obtained after integrating the square of the Fourier transform of the observed signal over a frequency range, and weighted by the inverse of the power spectral density of the detector (see for instance [25,26]).
The term F(Θ, s, p, θ) accounts for the geometric configuration, including the angular position in the sky, orientation of the detector, spin, polarization and orbital inclination of the binary. The term ζ encapsulates the detector response through the noise spectral density and the maximum frequency of the GW, f max . This maximum frequency (twice the frequency of the last orbit) is determined by the innermost stable circular orbit (ISCO). The ISCO frequency scales as the inverse of the chirp mass. Larger chirp masses correspond to smaller f max , and consequently to smaller ζ, partially compensating the dependency of ρ with M 5/6 . More specifically, the term ζ scales as fmax fmin df f −7/3 /S(f ) [26], where S(f ) is the power spectral density of the detector. For the range of frequencies and sensitivity in LIGO, one can use the published S(f ) for the observed LIGO events to estimate ζ. For the typical BBH events detected by LIGO, ζ has a mild dependency with f max . In particular, ζ scales roughly as f 0.1 max (or as M −0.1 ) for f max between 300 Hz and 700 Hz, which is the typical range of the observed BBH in LIGO. This modest dependency with f max is due to the fact that most of the power of the signal is contained in the lower frequencies, so above 200 Hz, the contribution to ρ is relatively small. The focus of this paper is in the often neglected first term in equation 2. This is the square root of the magnification, µ, which is usually assumed to be µ = 1, and hence ignored [27]. If a GW is being amplified by strong lensing, by ignoring µ, the inferred luminosity distance will be smaller than the true distance by a factor √ µ, and the inferred chirp mass (in the rest frame) larger in order to maintain M constant. This mechanism has been invoked in the past to explain the unusually high chirp masses observed for some of the events detected in LIGO/Virgo [18]. A fundamental ingredient in this work is the distribution of chirp masses. Little is known about this distribution. The observed events by LIGO/Virgo are consistent with standard power law models (M −2.3 ) but also with a bimodal function having a peak at high masses (at ≈ 40 M ) [28]. In [18], the authors make the interesting suggestion that the peak in the bimodal mass function could be the consequence of lensing of the GWs. Since a lensed GW can be misinterpreted (if lensing is ignored) as a closer GW with a larger chirp mass, lensed GWs could naturally produce observed shallow, or even bimodal mass functions. In this work we assume a simple model for the Chirp mass function. We start from a standard stellar model where heavy stars follow the standard power law M −2.3 , appropriate for masses above 1 M [29], and assume that the remnants (BHs) left by these massive stars have a mass that is proportional to the parent star. Individual masses (m 1 and m 2 in equation 1) are drawn from this distribution (with m 1 and m 2 in the range 5M < m 1 , m 2 < 50M ) and the Chirp mass is computed from the randomly generated m 1 and m 2 . This approach ignores possible correlations between m 1 and m 2 . The resulting distribution is shown in Fig. 1  The resulting mass function resembles the low mass peak of Model C considered in [28], although with a lower mean mass (see Fig.1 in that reference). It is important to note that, as demonstrated below (and originally proposed by [18]), the high mass peak of model C in [28] can be naturally produced by lensing, provided the rate of events at high redshift is sufficiently high.

IV. PROBABILITY OF OBSERVING A LENSED GW
The probability of lensing is referred to as the optical depth of lensing, τ , and in its most basic form depends on the magnification, µ and the redshift of the source, z s , i.e τ (µ, z s ). Different authors have estimated τ (µ, z s ) adopting different assumptions. Most estimates of τ are based on analytic models which have the advantage of flexibility, speed and achieving high spatial resolution (needed to resolve small areas with larger magnifications). On the negative side, analytic calculations do not account for corrections such as projection effects or substructure that can increase the optical depth. N-body simulations can account for both projection effects and substructure but at the expense of resolution. In this work we follow [30] and use an analytic model based on the mass function of [31] and an elliptical halo model to compute τ (µ, z s ). We improve on the work of [30] and modify the elliptical profile by adding a steeper component in the central region. This component accounts for the baryonic contribution that is important in smaller haloes and was neglected in [30]. We also extend the mass range down to 10 11 M from the minimum mass of 10 12 M considered in [30]. The baryonic contribution makes smaller haloes more relevant, and even though their contribution to τ is still subdominant, they can not be ignored. An additional improvement involves computing the optical depth in the image plane rather than in the source plane. In [30], the author used the optical depth computed in the source plane and later corrected for the multiplicity of lensed images. The optical depth computed in the source plane accounts for the total magnification and it is the best choice when estimating rates of unresolved images. Since GWs can be resolved in time (two counterimages from the same event will arrive separated by a time interval that can range between a few hours to a few days), the optical depth computed in the source plane has the advantage that no corrections are needed to account for multiple images. It also has the advantage of achieving higher spatial resolution (since there is no need to map the image plane into the smaller source plane). The improvement in spatial resolution translates into an improvement into the maximum magnification that can be computed before limited resolution affects the computation. Our calculation can reach magnification factors of 100 before being affected by resolution effects (even for the smallest haloes). Above magnification 100, one can safely extrapolate the optical depth with the standard µ −3 law as τ (µ >100 , z s ) = τ (µ 100 , z s )(100/µ) 3 . In Fig. 2 we show examples of the optical depth as a function of z s and for different values of the magnification factor. Note how between µ = 10 and µ = 100 the optical depth has decreased by two orders of magnitude, as expected from the scaling P (µ) ∝ µ −3 in strong lensing. This result is similar to the one the one derived from N-body simulations and ray tracing [32,33].
For magnification values µ > 1010, and for redshifts z s > 1, the probability of lensing is ≈ 10 −4 . This number can be confronted with the number of observed strongly lensed QSOs (few dozens) out of the total number of known QSOs (> 10 5 ). Between z = 1 and z = 2, one expects that 1 in ≈ 4 × 10 5 events will be magnified with factors µ > 30. Since the volume in this redshift interval is ≈ 450 Gpc 3 , if the rate of GWs above z = 1 is above 1000 events per Gpc 3 and year, one would expect to observe at least 1 strongly lensed event above z = 1 per year and with µ > 30. The rate of events as a function of redshift is discussed in the next section.

V. RATE EVOLUTION OF BBH
The rate of GWs as a function of redshift is unknown and depends on the formation mechanism of BBH. The simplest models rely on the assumption that the rate of GW events trace the star formation rate. We adopt the model of [34] as a conservative one (and we denote it as the SFR case). This model predicts that the evolution of the rate of GWs is mild in redshift and peaks at z ≈ 2. We normalize this model at redshift z = 0 to a rate of a few tens of events per year and Gpc 3 , consistent with the inferred rate from the LIGO collaboration ( [28], under the assumption of no lensing). This model is shown as a dashed line in Fig. 3. As an alternative case [inspired by 18, and that we refer to as the BDS FIG. 3. Rate of events as a function of (lookback) time. The solid line represents a model that grows as the star formation rate at high redshift but below redshift z ≈ 1.8) it falls of exponentially with a half-life time of 1 Gyr (in the main text, we refer to this model as the BDS model). The dashed line is a model that traces the star formation rate history. For comparison, the volumetric rate of SNe (of all types) at z ≈ 2 is approximately 10 6 per year and Gpc 3 model], we consider also a model in which lensed GWs occur more frequently at high redshift, and has the attractive feature of explaining the apparent bi-modality of the mass function of BHs, as well as the unexpectedly high number density of very massive events (with chirp masses of 30 M or above). In this model, the rate of events at high redshift (where the lensing probability is highest) is increased by ≈ 2 orders of magnitude (but still comfortably below the rate of SNe at all redshifts) in order to compensate for the low probability of lensing. To keep the observed rate of events comparable to the ones predicted by the conservative model, the rate at the lowest redshifts must be smaller. We attain this by allowing the rate to decay exponentially between the maximum of the star formation rate at z ≈ 2 and z = 0. Using a half-life time of 1.5 Gyr results in similar predicted number of observed events for both models so we adopt this value. The fast decline in the rate can be justified by the rapid evolution of BBH if these are formed in dense environments like globular clusters [see e.g 35].
In this case, the BHs are formed after the SN explosion of the parent star (typically a short-lived massive star), and through mass segregation they sink towards the center of the cluster (since after the massive stars undergo SN, the most massive objects in the globular cluster are still their remnants). Typical relaxation times are 1 Gyr or less [36][37][38]. The duration of the phase between the BBH formation and the final coalescence is more uncertain and depends on mechanisms like dynamical friction and three-body interactions. Most models agree that, once the binary is formed, only a small fraction of the binaries will merge within a Hubble time due to radiation of GWs [39, see for instance]. Other mechanisms may need to be invoked to explain a rapid evolution in the merger rate between the peak of star formation and redshift z = 0.
Also, it is important to consider the fact that massive BHs are typically produced from massive stars with low metallicity [40]. Below z ≈ 2, massive stars with low metallicity are expected to be exceedingly rare. Hence, one may expect that BBH with elevated chirp masses are mostly produced at high redshifts and, if the time of BBH formation and coalescence is relatively short, a rapid evolution between redshift 2 and redshift 0 is expected. Above z ≈ 2, the rate decays with redshift in a similar fashion as in the SFR model.
The BDS model is shown as a solid line in Fig. 3. Note that we choose to express the rate as a function of lookback time in order to accentuate the extended period up to the epoch of maximum star formation rate. Interestingly, the LIGO collaboration finds evidence for a rapid evolution in the rate (R ∝ (1 + z) 6.5 , [28]) when GW170729 is included in their analysis. On the contrary, if GW170729 is excluded, the rate is consistent with a non-evolving law, reflecting the large degree of uncertainty on the redshift dependency of the merger rate. Note that the BDS model is consistent with the model C in the above reference when the 10 BBHs from O1 and O2 are included in the analysis [See the redshift evolution model in Figure 6 in 28]. On the other hand, it is worth noticing that the relatively high rate of events at high redshift of this model may be in tension with upper limits on the stochastic background of GWs [41]. However, the sensitivity of LIGO to the stochastic background is still ≈ 2 orders of magnitude above the current estimate of this background [42] so a model such a BDS is still marginally consistent with the lack of detection of this background. Moreover, since GWs at low redshift (with a lower rate in the BDS model at z < 2) contribute more to the stochastic background, a model such as the BDS model is expected to produce a background less than 2 orders of magnitude above the model used in [42] VI. EXPECTED NUMBER OF STRONGLY LENSED EVENTS With the ingredients presented in the previous sections, we can now compute the expected number of strongly lensed events. This number is basically given by the integral, is the rate of events discussed in the previous section, dN /dM is the chirp mass function discussed in section III (that to first order we assume it does not evolve with redshift). The integration mass limits in the integral are fixed to M a = 5 M and M b = 50M . P (Θ) is the probability distribution for the geometric factor, encoding all possible orientations of the detector and GW. For this term we follow [26]. Finally, µ min is the minimum magnification required for a GW at redshift z, with chirp mass M and geometric factor Θ to be above the detection threshold of the experiment (see Eq. 2). We set a threshold in the signal-to-noise ratio of ρ = 8 for noise properties similar to those in O1 and O2 runs in LIGO/Virgo (see section II). By setting µ min = 1 we compute also the number of events that are not being lensed. The resulting number of observed events it is shown in Fig. 4 as a function of the redshift of the GW and for the two rate models discussed in section V. The dashed line (SFR model) shows results consistent with earlier estimates that predict that a small probability of lensed events is expected. The peak in lensed events at z ≈ 0.3 corresponds to events with a mild magnification of ≈ 2 for which the optical depth of lensing is higher (see Fig. 5). On the contrary, the BDS model (solid line) predicts similar number of lensed and not lensed events. The lensed events originate mostly at z > 1. As in the SFR model, a peak of low magnification can be appreciated also at z ≈ 0.3 Fig. 5 shows the distribution of magnifications that result in observed events. For the SFR model, most of the events are magnified by moderate factors which makes them difficult to be recognized as lensed events. However, for the BDS model, most of the lensed events have magnifications of a few tens, and a non-negligible fraction can have magnifications above 100.
Also interesting is the distribution of inferred chirp masses. Figure 6 shows the chirp masses inferred for the not lensed and lensed events. As expected the inferred masses for the not lensed events resembles the underlying chirp mass (with a bias towards higher masses because lower masses are less likely to exceed the detection threshold). However, for the lensed events (for which the mass has been inferred assuming the wrong magnification, µ = 1), the peak of the distribution is around masses of 20 or 30 solar masses, in striking resemblance to the current observed LIGO masses. [18] originally suggested that this mechanism was responsible for the apparently high chirp masses of most of the LIGO events.

VII. LENSING BY A MACROLENS PLUS MICROLENSES
This section briefly reviews the lensing formalism for microlenses in a macrolens. This topic has been extensively covered in the literature [43][44][45][46]. Our model involves a macrolens and a population of microlenses. The microlenses include stellar microlenses (stars and remnants) and a population of PBHs. The PBHs account for a fraction f of the total dark matter where f = 1 would imply that all dark matter is made of PBHs. Since in the range of ≈ 5 − 50M , the fraction f = 1 is already excluded by other observations, we consider only values consistent with current constraints. The surface mass density of dark matter is given by the convergence, κ of the macrolens. Then, the surface mass density of PBHs is simply f κ.
For the macrolens, we follow [15] and define the macrolens with just two parameters, the magnification factors in the radial and tangential directions, or µ r and µ t respectively. This is a valid approach for the macrolens since we are dealing with very small regions of the sky. Without loss of generality, we assume that µ t µ r and that the main direction of the shear, γ, is oriented in the horizontal direction, that is γ 2 = 0 and γ = γ 2 1 + γ 2 2 = γ 1 . Given µ t and µ r , the corresponding values of κ (convergence) and γ (shear) can be found from the relation between κ, γ, µ r and µ t . For a given choice of κ, and γ, the lens equation ( β = θ − α( θ)) of the macrolens can be expressed as FIG. 6. Inferred observed chirp masses in the source frame assuming there is no magnification. Note how the lensed events appear to have a higher chirp mass and the apparent bi-modality of the observed mass function.
where the positions in the source plane are given by the coordinates β = (β x , β y ) and the positions in the image plane are given by the coordinates θ = (θ x , θ y ). The lensing potential of the macrolens, ψ, is given by where we remind the reader that we adopt a reference system where γ 2 = 0, θ x and θ y are given in radians, and we ignore a constant additive term (i.e., the potential is identically zero at the origin of coordinates of θ).
Since both the deflection field and lensing potential are linear with the addition of new masses, if a population of N point masses are present, the deflection, α P S ( θ), and potential, ψ P S ( θ), from the distribution of point masses can be simply added to the above equations with; and, where and D s (z s ) the angular diameter distances between the lens and the source, between the observer and the lens, and between the observer and the source respectively.
A quantity of interest, that will be relevant in sections VIII and IX, is the effective optical depth, τ ef f introduced by [16], where the total magnification (µ) is the product of the tangential and radial magnifications (i.e., µ = µ t × µ r ), and Σ (expressed in units of M /pc 2 in the expression above) is the microlens surface mass density. When τ ef f ≈ 1, the saturation regime is reached. In this regime, caustics always overlap in the source plane, and any source moving across a field with τ ef f > 1 will be constantly experiencing microlensing [16,30]. Since typical values for Σ(M /pc 2 ) range between a few to a few tens, and assuming a typical value for µ r ∼ 1, it is clear from the expression above, that the saturation regime is reached when the macrolens magnification is in the range of a few hundred to a few thousand. Similar values for the macrolens magnification have been observed already at the position of the microlensing events of the Icarus and Warhol high redshift stars [47,48]. Given the fact that GW experiments have access to the entire sky at any given moment (as opposed to the aforementioned examples of Icarus and Warhol where a significant luck-factor had to be involved) events with similar, or even more extreme, macrolens magnifications are expected [see 30, for a detailed estimation of the probability of these events].
In those scenarios, we expect microlensing to play a significant role. Finally, the time delay is given by Where we have assumed that all point masses are at the same redshift in the lens plane so the factors D i (z l , z s ) are the same for all of them. The time delay can be expressed in dimensionless units by re-scaling both the angular positions and potential by the Einstein radius, θ 2 E = (4GM/c 2 )D(z l , z s ). This redefinition of the time delay expression makes most sense when one is dealing with a single microlens since in this case the Einstein radius can be defined without ambiguity, but in general one can still set the Einstein radius to any arbitrary mass, or scale, and still redefine the time delay equation.
where x = θ/θ E , y = β/θ E ,ψ = ψ/θ 2 E , and R s (z l ) is the redshifted Schwarzschild radius of the lens. The first term sets the scale of the time delay for a given lens mass, 2R s (z l )/c = 1.97 × 10 −5 (1 + z l )(M/M ) seconds. This simple scaling, shows that for GW with frequencies ν ∼ 100 Hz, interference and diffraction effects are expected for isolated microlenses with masses > 500/(1 + z l ) M [7,10]. Below this mass, the Schwarzschild radius of the lens is smaller than the wavelength of the GW (with ν ∼ 100 Hz) and diffraction effects are not important resulting in the GW not seeing the microlens. However, as discussed in the section below, this minimum mass can be lowered substantially when the microlens is near a macrolens critical curve.

A. The impact of microlenses on gravitational waves
In earlier work, [16] discussed how a microlens with mass M embedded in a region of a macrolens where the macrolens magnification is µ, behaves as a microlens with an effective mass µM . Also, at larger magnifications, an area A in the lens plane, maps into a smaller area A/µ in the source plane. Hence the density of microcaustics increases by a factor µ. Since the surface mas density of microlenses is roughly a constant fraction of the convergence (which is of order 1 near the critical curves of macrolenses), it is then easy to realize that a GW that is being magnified by a macrolens with a factor µ, will inevitably be affected by microlensing (and its associated wave effects), if the magnification factor of the macrolens is large enough.
When the wavelength of a GW is comparable to the Schwarzschild radius of the microlens, one needs to consider wave optics instead of geometric optics. The magnification factor in wave optics is given by the diffraction integral [see for instance 49].
where ν is the frequency of the GW in Hz and the normalization A o guarantees that at very high frequencies one recovers the geometric optics magnification when averaged over a frequency range. ∆T is the time delay between lensed images expressed in seconds. The total magnification, µ, and phase shift, φ, are given by [see for For isolated microlenses and simple macrolenses, analytic solutions can be found for the magnification. The magnification depends on the frequency and that at low frequencies the magnification tends to 1 (i.e the microlens has a negligible effect over the GW). In realistic situations, where a GW is strongly lensed, it may intersect not only one microlens but different microlenses. In this case, the integral in Eq. 11 can be solved numerically as described in [8,15,50]. The case of a pair of overlapping microlenses was studied in [30]. In this paper we study the more realistic situation where a distribution of microlenses intersects the GW. The probability that a GW intersects a microlens and suffers a distortion in its waveform depends on the density of microlenses and the size of the caustic region. As discussed at the beginning of this section, this probability increases with µ. Hence GWs with large magnification factors are more likely to suffer interference than GWs with modest magnification factors. We consider two scenarios, the first one where the macrolens magnification is 30 and a second scenario where the macrolens magnification is 5 times larger (that is, 150). As shown in section VI, the most likely magnification factor for models such as the BDS one is µ ≈ 30.
Despite the rapid decline of the optical depth with magnification, observed events with µ = 150 are only ≈ 5 times less likely than observed events with µ = 30 (for the BDS model) since the reduction in lensing probability is partially compensated by the increase in the accessible volume. For the SFR model, the rate of probabilities is similar but the probability of observing events with µ = 30 is two orders of magnitude smaller, so events like the ones described in this work would be observed only after the sensitivity of the detectors improve.

VIII. PRIMORDIAL BLACK HOLES AS DARK MATTER
In the previous sections we discussed the lensing probability for the SFR and the BDS models. In both cases, lensing events are expected to take place although in the case of the BDS model, these events are ≈ two orders of magnitude more likely, and should have already been observed by LIGO. As mentioned above, the BDS model is appealing because it predicts a bi-modal mass function in the observed chirp mass, consistent with the observed distribution of chirp masses. On the contrary, if the BDS model is correct, lensed events with inferred large chirp masses should come in pairs with time delays between a few hours to several days, and from the same location in the sky. To this date, no clear multiply lensed GW has been observed [51][52], although several candidates have been proposed [53]. More recently, two events during the 03 run were observed on August 28th 2019, which are separated in time by ≈ 19 minutes, and seem to originate from similar luminosity distances and positions in the sky. These positions are separated by only ≈ 10 degrees. To first order, the probability of having two events so close in time and space is comparable to the probability of having a strongly lensed event. How-ever, based on the sky location of these two events, [54] concluded that the two events can not be interpreted as a pair of strongly lensed images. As shown below, if GWs are strongly lensed and microlenses are present, the observed strains can be perturbed, specially at the largest frequencies. The impact of this perturbation on the estimation of the parameters of the GW is something that needs to be addressed carefully.
In this section we assume that strongly lensed GWs have been observed already (as suggested by the BDS model), or that they will be observed in the near future, and investigate the possibility of constraining the abundance of PBHs. We take advantage of the magnifying power of the macrolens hosting the PBHs, which increase the concentration of caustics in the source plane, resulting in complex time delay distributions and interference patterns of the lensed GWs. We create a set of simulations where we vary the fraction of dark matter contained in PBHs, and compute the distortion in the magnification as a function of frequency for different positions of the GWs in the source plane. Sampling the source plane allows us to perform a statistical analysis, where we estimate the probability that a GW is distorted by some fraction at a given frequency. The simulations involve three components, a macrolens (galaxy or cluster), stellar microlenses from the macrolens, and PBHs that constitute a fraction, f , of the dark matter of the macrolens. By construction f < 1, and we consider fractions, f , that are consistent with current upper limits. We neglect projection effects and only consider deflectors that are in the macrolens plane. Projection effects are expected to contribute at the percent level since the adopted value for the convergence of the macrolens is of order 1. For the macrolens, we follow [15] and model the macrolens with just two parameters, the convergence κ, and the shear γ. The specific values of κ and γ are derived from predetermined values of the magnification in the tangential, µ t , and radial directions, µ r . We consider two scenarios for the magnification, µ, of the macrolens. In the first scenario, we adopt a relatively modest magnification of 30, resulting from the product of a tangential magnification of 10 and a radial magnification of 3, that is µ = µ t µ r = 10 × 3 = 30. In the second scenario, we consider a more extreme value of the magnification µ = µ t µ r = 50 × 3 = 150. The two scenarios can be considered as two nearby positions in the image (or source) plane where the second scenario corresponds to the position that is closest to the critical curve (or caustic). For typical macrolenses with the mass of a massive galaxy, the magnification can change from 30 to 150 in image plane distances smaller than one arcsecond. In addition to the two considered values for the macrolens magnification, we also consider the two possible parities of the macroimages. A source located near a caustic produces typically two bright macroimages (in general additional macroimages are also formed, but with significantly less magnification and time delays much larger than the duration of a GW event, so they can be neglected). One of the two macroimages has positive parity (i.e, it has the same orientation of the orbital momentum as the source) and the other has negative parity (i.e it resembles a mirror image of the source). Other than the change in symmetry, in the presence of microlenses in the lens plane, the parity of the macroimage is important when dealing with very small sources. As shown in earlier work (see for instance [15] for a discussion in the context of GW), macroimages of very small sources with negative parity are significantly more likely to have smaller magnifications than the corresponding macroimage with positive parity, if microlenses intersect the line of sight. The effect is even more dramatic for GWs since the microlenses introduce a time delay that can be more pronounced between microimages forming around a macroimage with negative parity. The time delay between microimages results in constructive and destructive interference patterns affecting the GW. This effect depends on the frequency of the GW and can resemble (or can be confused with) spin misalignment. Since the results depend on the parity of the macroimage, this allows to confirm that a GW is being strongly lensed, since one would expect (on average) that the image with positive parity arrives first and a few hours (or days) later the second image with negative parity would arrive. This second image is expected to show stronger distortions with frequency, on average, if microlensing is involved. If µ is sufficiently large (µ > larger than a few hundred), this type of behavior is always expected since, even in the absence of PBHs, microlenses from the intergalactic or intracluster medium can saturate the source plane with microcaustics. Events with extreme magnification can be recognized because they may originate at very high redshifts (z > 2) and be misinterpreted (if lensing is ignored) as events with unusually high masses m 1 and m 2 , that could even fall within the mass gap above 50 M .
In all simulations we include the stellar component using the same model described in [15]. For the two macrolens magnifications considered (30 and 150), we simulate a circular area in the image plane of 3 mas in diameter, but consider only the central horizontal rectangular area of 3×0.5 mas 2 . The resolution is 100 nanoarcsec per pixel, sufficient to properly resolve the caustic regions around microlenses with stellar masses. The area in the source plane is reduced by the corresponding magnification factor of the macrolens. The contribution of the stellar component is modest and only the (rare) most massive stars and remnants are expected to produce significant effect on the GWs. The contribution to κ of the stars and remnants is typically at the ∼ 1% level, and most of the mass is contained in objects with masses ∼ 1 M or less. At masses above ∼ 5 M , where the interference effects are more prominent, stars and remnants contribute to the total mass a fraction much smaller than the contribution of the PBHs considered in this work. Because of this, the impact of the particular modeling of the stellar component is expected to be minor in our conclusions. Nevertheless, they are included assuming a surface mass density of 11.7 M pc −2 which is the same value used in [15], and a typical value found around the critical curves of macrolenses for background sources with z > 1 [47,55].
For the mass function of the PBHs, we assume that their mass function is not necessarily the same as the one assumed for the sources (see section III). Instead, we assume two different models. First, we consider a Gaussian distribution. This is a simple extension to the monochromatic model with a fixed mass, that allows for a broadening of the mass function. A log-normal mass function could be used instead, but we expect negligible differences in our conclusions, provided both Gaussian and log-normal have comparable dispersion. Similar models have been considered in the past to describe a possible population of PBHs with masses comparable to those found by LIGO. This type of model is already constrained by current observations which rule out fractions larger than f ≈ 0.1. In this work we consider fractions that do not violate these limits. In particular we consider values of f = 0.033, 0.01, 0.003, that is, 3.3%, 1%, and 0.3% the amount of dark matter respectively. We refer to these three models as Gaussian models and more specific details above these models are provided below. Interestingly, a fraction ≈ 1% may be sufficient to explain the rate of mergers observed by LIGO [6,56]. As a second model for the population of PBH, we consider a much broader mass function, which is still consistent with current constraints [57], but can allow for a larger fraction of the dark matter to be explained by PBHs. This is possible as long as the PBH span a wide range of masses, and the fraction of mass contained in each range does not violate current limits [57][58][59]. We model this second model as a power law between M min and M max and refer to this type of model as power-law model.
where κ is the convergence, Σ crit is the critical surface mass density, and Ψ(M ) is the mass function of the PBH per unit area. We assume that it follows a power law, Ψ(M ) = dN/dM ∝ M γp−1 . The index, γ p , is a free parameter and we consider two values, γ p = 0 and γ p = 1. When γ p = 0, f is a constant as a function of mass. When γ p = 1, the number of PBHs per mass interval is constant and hence more massive PBHs contribute more to f than lower mass PBHs. This mass function was also considered in [60] that argues that "a mass function of this form arises naturally if the PBHs form from scale-invariant density fluctuations or from the collapse of cosmic strings". For the limits of the integral, we consider M min = 5 M and M max = 50 M . This range is justified because wave effects are not very important for masses below 5 M . On the other extreme, above 50 M the fraction of PBH is heavily constrained. Outside this mass range, there may still be PBHs that contribute to f . To account for this possibility, we consider values of f smaller than 1, but still larger than the values of f adopted for the Gaussian model. Note however that in [61] the authors propose a mass function that, they claim, is consistent with current microlensing observations and can explain all dark matter as PBHs. In particular, in this model the mass function of PBHs is multi-modal, with the peaks of the mass function concentrated around 10 −6 , 1, 30, and 10 6 M .
For a given choice of the mass function, the number of PBHs per unit area inserted in the simulation is determined by the specific values of f and κ. The convergence is determined by the values of µ t and µ r as described in [16].

A. Gaussian Model
First we consider the case of the Gaussian model. We assume a Gaussian mass function for the PBHs with mean 30 M and dispersion 5 M . This model is inspired by the high-mass peak found in the bi-modal mass function in [28]. Since the dispersion is relatively small, the results presented in this section would be close to those obtained with a monochromatic mass function with masses equal to the mean value of the Gaussian. In Figure 7 we show the magnification in the image plane for the simulations in the Gaussian model case and for the case with µ = 150. The magnification diverges at the critical curves (white lines). The top 3 panels show the magnification on the side of the image plane where the parity is positive and for the three values of f that are consistent with current constraints [62]. Form top to bottom f = 0.3%,f = 1%, and f = 3.3%. The bottom three panels are the corresponding magnification maps in the image plane for the side of the image plane with negative parity and the same values of f . In both cases (positive and negative parities), the dark-gray areas have magnification factors of a few. The largest critical curves are around the PBHs. Smaller critical curves correspond to stars and remnants in the macrolens. Note how the orientation (and shape) of the critical curves depend on the sign of the parity. The case with µ = 30 is not shown but would look similar with the exception that the critical curves around the microlenses would be smaller by a factor 150/30. By using the lens equation, we map the magnification in the image plane into the magnification in the source plane, where the critical curves map into caustics. Since the image plane gets compressed by a factor µ in the source plane, it is often the case that, at large macrolens magnifications, the caustics overlap. In this scenario, images can form around critical curves that are separated by relatively large distances in the image plane. For these images, the geometric time delay can be larger than the time delay that would be observed if the caustics did not overlap, resulting in interference effects at even lower frequencies. This overlapping effect in the source plane is shown in Figure 8 where the three top panels are for the case with positive parity and the three panels in the bottom are for the case with negative parity. In this case, the caustics look very different depending on the sign of the parity. On the side with positive parity, the minimum magnification (dark gray areas) is of the order of the macrolens magnification (i.e, µ min ≈ 150) while in the side with negative parity the dark-gray ar- eas have magnification of just a few. For both parities, at fractions f above a few percent and at this macrolens magnification, caustics start to overlap across the source plane. At even larger macrolens magnification factors, the overlap start to take place at smaller values of f .
We use the simulated source planes to estimate the probability that a GW is distorted by some amount at a given frequency. We place 300 GWs at random positions in the source plane shown in Fig. 8 and compute the magnification of the GW (relative to the macrolens magnification that does not depend on frequency) as a function of frequency following [15] for each position. Each of the 300 random positions (β in Eq. 11) results in a distortion curve of relative magnification as a function of frequency (w in Eq. 11). In Figure 9 we show the first 50 curves, out of the 300 examples for a particular choice of the PBH (that is described below in subsection VIII B). Let us remind the reader that these curves represent the relative magnification with respect to the macromodel and that the macromodel magnification is expected to be independent of frequency. For clarity purposes, we highlight the first three curves as colored solid lines (blue, orange, and red). As shown in the figure, some of the curves show significant deviations from unity at frequencies as low as 100 Hz. For instance, the dark blue solid curve shows a large decrement between ≈ 80 Hz and ≈ 150Hz. The same curve, shows an increase in the magnification (relative to the macrolens) at ≈ 400 Hz, followed by a large reduction in the relative magnification at ≈ 500 Hz. All curves tend to 1 at lower frequencies, as expected, since at these frequencies the GWs become insensitive to microlenses with the masses considered in this work.
With the 300 magnification curves, we compute the probability that a the microlensing distortion of GW is larger than some amount as a function of frequency. For this goal we compute the envelopes of the curves at different percentage levels (50%, 90% and 97%), and as a function of frequency. The envelopes at the 50% level indicate the minimum distortion (at a given frequency) in the magnification of 50% of the curves, that is 1 in 2 GWs will have a distortion at least as large as the one shown by the envelope at 50 % at that frequency. The envelopes at the 90% level indicate the minimum distortion in the magnification of 10% of the cases. That is, 1 in 10 GWs, will have a distortion in the magnification at a given frequency that is at least as large as the one shown by the envelope. Similarly, for the envelope at 97%, approximately 1 in 30 GWs, will have a distortion at a given frequency at least as large as the one shown by the envelope. Figure 10 shows the envelopes (or probabilities) for the case with macrolens magnification µ = 30. The cases where the macrolens parity is negative are shown in the left column while the cases where the parity is positive are on the right column. From top to bottom, we show the cases with f = 0.033, 0.01, 0.003 and f = 0. The case with f = 0 has only the stellar microlenses (stars and remnants with a surface mass density of 11.7 M /pc 2 following the model of [63]), and no PBHs. The serrated pattern observed in the envelopes with negative parity at low frequencies is an artifact due to numerical limitations (and related to the limited size of the simulation) when solving numerically the integral in equation 11. The envelope at this low frequencies is supposed to be a smooth version of the curves. As expected, the distortions are larger at higher frequencies. If one considers the case with f = 0.01 (or 1% of the dark matter in the form of PBHs) and positive parity (i.e., right column, second panel starting from the top) at 500 Hz, then 10% (90% envelope) of the strongly lensed GWs with µ = 30 are expected to show a strong distortion in the magnification where destructive interference demagnifies the GW by a factor ≈ 3 (relative magnification ≈ 0.3). Since the strain goes as the square root of the magnification, this translates into a reduction of a factor ≈ 2 in the amplitude of the strain at 500 Hz. If the fraction of dark matter in PBHs increases, at f = 3.3%, destructive interference with even greater intensity is evident at frequencies as low as ≈ 200Hz, well within range of detectors such as LIGO. The distortion is even more apparent for the counterimages with negative parity (that we remind the reader should arrive between a few hours to a few days after the GW with positive parity) as shown by the left column. In this case, at least 1 in 10 GWs with magnification above µ = 30 should show a reduction in the strain by a factor ≈ 2 at frequencies of 200 Hz. We should note that the envelope accounts only for the most extreme distortion at that frequency. This means, that the event mentioned above with a reduction of a factor ≈ 2 at frequencies of 200 Hz, at higher frequencies may show a smaller distortion (or even no distortion at all in a given frequency range).
The above result shows the potential of strongly lensed GWs to probe the abundance of PBHs but it requires observing several events in order for one of them to show significant distortions in the observed strain. However, if we focus on events with even larger macrolens magnification factors, the caustics will overlap even more and individual microlenses will behave as microlenses with larger masses (with an effective mass that scales as µM ). In Figure 11 we show the case for a macrolens magnification µ = 150. Note that although the probability of an event with µ > 150 is 5 2 times smaller than the probability with µ > 30, the probability of observing such event is actually only a few times smaller since at µ > 150 one is probing a larger volume than at µ > 30 (see for instance Fig. 5). For lensed events with µ > 150, and if one considers the case where 1% of the dark matter is in the form of PBHs, the 50% envelope is already at magnification factors 0.3 at frequencies ≈ 300 Hz. That is, a reduction in the strain of a factor ≈ 2 at 300 Hz is expected in 50% of the lensed events with µ > 150. One in ten of these events (90% envelope) would have a similar reduction of a factor 2 in the strain at frequencies as low as ≈ 130 Hz. The same fraction would show reductions of a factor 5 in the strain at 400 Hz.

B. Power Law Model
The second scenario we consider is the case where the mass function of PBHs between 5 M and 50 M is described by a power law. Given the current constraints, PBH can make up to ≈ 10% of the dark matter [57]. If one relaxes the constraints at very low masses (asteroids) from PBH evaporation into gamma rays and at masses above 50 M from the CMB (Planck), PBHs can make up to 100% of the dark matter. In [60], the authors explicitly consider extended mass functions with power laws similar to the ones considered in this work and find that in the mass regime we are interested in, a fraction as high as 10% is consistent with current constraints. We consider this value of f = 0.1, in the case of the power law and make a set of simulations similar to the ones in the previous section (µ = 30, µ = 150, positive, and negative parities) for the two power laws considered, γ p = 0, and γ p = 1.
The result is shown in Figure 12. In this case, and for simplicity, we consider only two envelopes at 50% and 90%. Not surprisingly, since there are more PBHs per unit area, the effect is larger than for the Gaussian model described in the previous subsection. As a general trend, we observe again that images with negative parity are more likely to show interference effects. The differences between the two power law models is relatively small. The difference is also small when comparing the µ = 30 and µ = 150 cases, except at low frequencies, where interference effects in the case with µ = 150 are more evident. This can be explained because at fractions f = 0.1, and for magnifications µ > 30, the saturation level is reached in the source plane [see 16, for a discussion of the saturation level], where microcaustics overlap over the entire source plane, and the statistical properties of the magnification are very similar at µ = 30 and mu = 150. Only at smaller fractions, f , the probability of interference becomes sensitive to the macrolens magnification for magnification factors µ > 30, as shown in Fig. 13. For smaller magnification factors, the saturation regime is not reached and the probability of intersecting a caustic becomes sensitive to the masses of the microlenses.
An interesting trend is observed in all cases with negative parity (left column in Figs. 12 and 13), where in more than 50% of the cases the magnification decreases with increasing frequency. In the 90% envelope, both for the images with positive and negative parity, GWs are heavily suppressed above 100 Hz. For the images with negative parity, this happens at even lower frequencies. In particular, 10% of the GWs with magnification µ = 150 would be heavily suppressed below 60 Hz. Since most of the signal to noise in the events observed by LIGO is concentrated below 100 Hz, this means that these events would basically be unobserved. In 50% of the cases, at µ = 150 the suppression factor is typically a factor > 2 at frequencies > 100 Hz for images with negative parity and a bit smaller than 2 for the images with positive parity, and at the same frequencies.

IX. DISCUSSION AND CONCLUSIONS
The results in the previous sections show how, if a fraction of the dark matter is composed of microlenses with masses above a few M , a strongly lensed GW may reveal this population of microlenses through interference effects. In general, the distortion is larger the larger is the fraction of dark matter in the form of compact microlenses, but at sufficiently high magnifications, the saturation level is reached (when τ ef f ≈ 1, see equation 8) and the microcaustics overlap in the source plane. Once this saturation level is reached, the results are less sensi- tive to the number and masses of the microlenses but the point where the saturation level is reached depends on the surface mass density of microlenses. Since the critical surface mass density for typical redshifts of lenses and background sources is Σ crit ≈ 2000 M pc −2 and near the critical curves κ ≈ 1, if the fraction of dark matter in microlenses is ≈ 3%, then τ ef f ≈ 1 when the magnification of the macrolens is µ ∼ 100 (assuming the radial magnification factor is ≈ 2). At smaller magnification factors, such as µ ≈ 30, one expects to be more sensitive to the particular mass function of microlenses but at the expense of having fewer GWs with a significant distortion in the magnification, since the area in the source plane which is not close to a microcaustic is larger in this case.
In the previous section, the envelope curves shown in Figures fig:FigStatMu1-12 account for the probability of distortion of certain magnitude at a given frequency. Al-ternatively, it is interesting to ask the question of what is the probability of having a magnification of at least a factor µ * (a strain distortion of a factor √ µ * ) below some frequency ν * . In Figure 14 we show the case for µ * = 4 and ν * = 200 Hz. Below this frequency, one can expect that a distortion in the strain of a factor > 2 at some frequency ν < 200 Hz could be picked by optimally filtering the raw data. Each curve represents all the models discussed in the previous section. The Gaussian models are one the left side of the x-axis while the power law models are on the right side. Dashed lines are for macroimages with negative curvature and solid lines are for macroimages with positive curvature. In general, and as expected, it can be appreciated how models with a larger fraction f are more likely to produce distortions in the strain. The exception is model C in the blue solid line which has fewer cases than model B in the same curve but this is purely a fluctuation due to small number statistics. The power law models are also more likely to produce distortions than the Gaussian models except when comparing models D (Gaussian with 3.3% dark matter) and model E (power law with 3.3% dark matter and γ = 0) which predict almost identical results, although with the Gaussian model predicting a few more cases with this magnitude in the distortion. Note how models G and H are close to the saturation regime and produce virtually identical results. The plot also shows the gain in probability of distortion when one considers events magnified by larger factors (red curves). If one considers for instance the Gaussian model with 1% dark matter at µ = 150 (red curves, model C), between ≈ 20% and ≈ 30% of the cases show distortions in the magnification of at least a factor 4 below 200 Hz. At even larger magnification factors, PBH models with smaller fractions of dark matter would show similar results. Note how the Gaussian model with 0.3% dark matter as PBHs (model B) is clearly more sensitive than model A (only stars and remnants) despite the fact that model A has a larger surface mass density than model B. However, most of the mass in model A is contained in low-mass stars, for which diffraction effects are small below 200 Hz. It is important to stress out that the mass function of the PBHs assumed in this work is not the same as the one assumed for the two components producing the GWs (see section III). In other words, we assume that the merger rate of BBHs is fixed and independent of the PBHs. For the GWs we assumed a chirp mass that results from drawing m 1 and m 2 from a power law which in our notation would have an exponent γ p = −1.3. Such a steep mass function would produce, on average, diffraction effects which are smaller than the ones presented in this work for γ p = 0, 1, since most of the mass would concentrate in smaller objects.
Although all the information of the distortion due to microlensing is encapsulated in the relative magnification as a function of frequency shown in the previous section, it is interesting to show the effect over the observed strain. In Fig. 15 we show examples for the case with magnification µ = 150, negative parity, and for the power FIG. 14. Probability of having a distortion in the magnification larger than a factor 4 at frequencies below 200 Hz. The solid lines are for images with positive parity while the dashed lines are for images with negative parity. Blue lines are for macroimages with magnification 30 and red lines are for macroimages with magnification 150. Model A is the one with no PBHs (i.e only stars and remnants with Σ = 11.7 M pc −2 ). Models B, C and D are the three Gaussian models with f =0.3%, 1%, and 3.3% respectively. Models E and F are power law models with f =3.3% and with γ = 0 (E) and γ = 1 (F). Models G and H are power law models with f =10% and with γ = 0 (G) and γ = 1 (H).  Fig. 9 in section VIII A. The time of coalescence takes place at t = t o . Note how in some cases, the amplitude in the last cycles is affected substantially or even suppressed almost completely. Also, it is interesting to see how on average (black-dashed line), the amplitude in the last cycles is reduced by some factor that increases with frequency. This trend was already observed in Figs. 12 and 13 but it is observed more clearly in this plot. The suppression of the last cycles is more prominent in the macroimages with negative parity than in macroimages with positive parity [see 30, where it is also shown how the mean of the magnification remains unchanged when microlenses are present]. This follows the expected behavior of lensing in the wave optics regime, where at low frequencies one is expected to be insensitive to the presence of microlenses (i.e, one recovers the mean of the magnification) while at high frequencies, the GWs can see the microlenses (i.e, the typical magnification tends towards the median of the magnification). If microlens-ing is predominant in strongly lensed GWs, and as shown by Figure 15, the last cycle before the ringdown is more likely to be damped (specially in macroimages with negative parity), it is possible that the estimated parameters of the GW (specially the one with negative parity) are biased if microlensing effects are ignored.
If one assumes that the BDS model is correct (i.e, all the low-frequency events detected by LIGO are strongly lensed with magnification factors up to several hundreds), this would imply that the fraction of dark matter in the form of PBHs between 5 M and 50 M must be below the 10% level. Otherwise, about half of the observed lensed GWs would show clear signs of microlensing in the strains. Similarly, if the fraction of dark matter in the form of PBHs in the same mass regime is of the order, or above, 10%, this would imply that the BDS model can not be correct, and that the observed low frequency (or high chirp mass) GWs are not strongly lensed. More data, and with higher sensitivity (specially at high frequencies) is needed to settle this question.
Future detectors of GWs, such as the Einstein telescope, will soon increase the sensitivity and extend the range in frequencies where events can be observed. Extending the range towards lower frequencies is interesting, not only for early warning reasons (events can be detected hundreds of seconds before the merger takes place), but also because at these low frequencies, microlensing effects are expected to be very small, and the parameters estimated at low frequencies can be later used to unveil signs of microlensing at high frequencies. Increasing the sensitivity at high frequencies would allow to measure, with higher significance, the ring-down (very sensitive to microlensing) and also detect events with smaller chirp masses at cosmological distances. Binary neutron stars (or BNS) are expected to be more common than BBHs (the inferred volumetric rate of BNS at z = 0 is ≈ 10 3 per year and Gpc 3 [64]). However, since the chirp mass is smaller, the signal-to-noise ratio is also smaller than for BBHs and hence they can only be observed at smaller distances (where the optical depth of lensing is small). Once the sensitivity of future GW detectors increase, it will allow to observe BNS at larger cosmological distances. If magnified by factors of a few tens to a few hundreds, future observations can reveal lensed GWs from distant BNS at z > 1. At higher redshifts, this rate could be even higher, offering a chance of observing these events through lensing. Interestingly, these events, when redshifted (assuming z ∼ 1−2), would appear to be falling within the forbidden mass gap at low mass. Also, the relatively smaller chirp mass results in a larger f max and it is at these higher frequencies where wave effects are more prominent.