Gravitational lensing of gravitational waves: A statistical perspective

Introduction

The four signals of gravitational waves (GWs) from binary black hole systems, GW150914 (Abbott et al. 2016c), GW151226 (Abbott et al. 2016d) , GW170104 (Abbott et al. 2017b), and GW170608 (Abbott et al. 2017a) detected by Advanced Laser Interferometer Gravitational Wave Observatory (aLIGO) during its first and second observing runs (O1, O2), marked the commencement of GW astronomy. More recently, with the Advanced Virgo detector becoming operational, we had the first joint detection GW170814 (Abbott et al. 2017c) and the first binary neutron star (BNS) signal GW170817 (Abbott et al. 2017d). These observations provide us a new opportunity to study astrophysics and cosmology.

Since Wang et al. 1996 proposed the possibility of observing several strongly lensed GW events in the context of aLIGO type detectors, gravitational lensing of GWs has been widely discussed over the past two decades. Such discussions involve diffraction effects in lensed GW events (Nakamura 1998; Takahashi & Nakamura 2003), the waveform distortion caused by the gravitational lensing (Cao et al. 2014; Dai & Venumadhav 2017), the influence on the statistical signatures of black hole mergers (Dai et al. 2017) as well as the potential for studying fundamental physics (Collett & Bacon 2017; Fan et al. 2017) and cosmology (Sereno et al. 2011; Liao et al. 2017; Wei & Wu 2017). Nevertheless, in spite of the broad range of topics discussed so far, the field of gravitational lensing of GWs is still worth an extensive exploration in order to fully understand the phenomenon and how to employ it to investigate the Universe.

One crucial question we have to answer before a further exploration of gravitational lensing of GW occurs is ‘how many lensed GW events are expected to be observed?’ Indeed, several discussions on this aspect already exist in the literature. For example, Sereno et al. 2010 studied lensed GW events from the merging of massive black hole binaries in the context of the LISA mission; Biesiada et al. 2014 considered the observational context for the Einstein Telescope (ET); and more recently, Ng et al. 2017 revisited the LIGO lensing rate. However, all the studies mentioned so far adopt the simplest lens model, which treats the lens mass distribution as axisymmetric.

In this paper, we present some extensions to the calculation of the lensed GW rate, making allowances for the ellipticity of the lens, the lens environment (as an external shear), and for magnification bias. This treatment not only provides a more precise prediction about the lensing rate, including more statistical properties, but also can serve as a useful tool for cosmological study (e.g. Chae 2003). We concentrate on the ground-based GW detectors, specifically aLIGO and the proposed ET. Nevertheless, the strategy developed here is general and can be easily extended to address other similar GW surveys as long as the geometrical optics approximation to GW propagation is valid.

The estimate of source rate dominates the prediction for the lensed event rate. Here we consider GWs from the coalescence of stellar binary black holes as the only sources, since they are the main signals received by ground-based detectors (Dominik et al. 2013; Abbott et al. 2016a). In order to obtain the source rate, we use a similar approach as in Cao et al. 2017 to estimate the merger rate of stellar binary black holes, and then use the GW detection theory developed by Finn 1996 to translate the intrinsic merger rate into the detectable source rate.

Another essential factor that can affect the observation of lensed events is the lensing time delay, as an image with time delay comparable to the survey’s duration has a high probability of being missed by the detector. We assess this selection bias by computing the distribution function of time delays corresponding to the lens properties adopted in this paper.

Our paper is organized as follows. In Section 2, we describe the approach to lensing rate calculation and the assumption of lens properties. We present our results in Section 3 and summary in Section 4. Throughout this paper, we adopt geometric units with $G=c=1$ and assume a Lambda cold dark matter universe with $(\Omega_{M},\Omega_{\Lambda})=(0.3,0.7)$ and a Hubble parameter $H_{0}=70~\text{km}~\text{s}^{-1}~\text{Mpc}^{-1}$.

Theoretical Model

In this section, we present our lens model (Section 2.1) and our GW detection model (Section 2.2). With the theory of lensing statistics (Section 2.3), we then derive the formulae to calculate the expected lensing rate in Section 2.4. The theory developed in this section is general and can be used to estimate strong gravitational lensing rates in any ground-based GW surveys so long as the geometrical optics approximation (see below) is valid.

When the lens mass is larger than $\sim 10^{5}M_{\odot}(f/Hz)^{-1}$ where $f$ is the frequency of the incident waves, the propagation of GWs is analogous to that of light. This is known as the geometrical optics approximation to GW propagation (Takahashi & Nakamura 2003). Since we here concentrate on the macrolensing by galaxies ($M\gtrsim 10^{10}M_{\odot}$) of high frequency GWs ($f\gtrsim 10$Hz), this condition is always satisfied. Hence, it is a reasonable approximation in the context of this paper to neglect the wave effect and adopt the standard optical gravitational lens theory to study the gravitational lensing of GWs.

As it is broadly reckoned that the strong lensing probability is dominated by early-type galaxies (Turner et al. 1984; Möller et al. 2007, and references therein), we only consider early-type galaxies as lensing objects. The singular isothermal ellipsoid (SIE) is adopted to model the mass distributions of the lensing galaxies. For the SIE convergence in Cartesian coordinates ($x,y$), we adopt the form developed by Keeton & Kochanek 1998:

where $q$ is the projected minor-to-major axis ratio, and $\lambda(q)$, the so-called ‘dynamical normalization’, depends on the three-dimensional shape of lensing galaxies (Chae 2003).

Furthermore, we consider the influence from the lens environment as an external shear $\bm{\gamma}$ whose lens potential is given by (Kochanek 1991; Witt & Mao 1997, and references therein)

where $(\gamma_{1},\gamma_{2})$ are the two components of the shear in Cartesian coordinates, and $(\gamma,\theta_{\bm{\gamma}})$ are the corresponding amplitude and direction components in polar coordinates. The connections between these two coordinate systems are: $\gamma_{1}=\gamma\cos 2\theta_{\bm{\gamma}}~,~\gamma_{2}=\gamma\sin 2\theta_{\bm{\gamma}}$.

More detailed discussions of the lens model can be found in Appendix A.

2 GW modelling

An estimate of the GW event rate density is required for calculating the expected number of lensed events. This involves the theory of GW detection, which has been discussed by many authors (Finn & Chernoff 1993; Finn 1996; Flanagan & Hughes 1998; Taylor & Gair 2012). Here we mainly follow the framework developed by Finn 1996.

For Gaussian and stationary noise, the optimal matched filtering signal-to-noise ratio (S/N) $\rho$ is defined as (e.g. Flanagan & Hughes 1998)

where $S_{n}(f)$ is the one-sided power spectral density of the detector’s noise, and $h(f)$ is the Fourier transform of the detector’s response to the GWs.

The GW generated by an inspiralling binary system can be approximately described by a quadrupolar formula (Newtonian order) with the frequency twice the binary’s orbital frequency. This waveform model does not meet the empirical requirement coming from the analysis of GW data, but is accurate enough for our statistical purpose. The amplitude given by the quadrupolar formula can be written as (Taylor & Gair 2012)

where $D_{L}$ is the luminosity distance and

is the observed (redshifted) chirp mass with $\mathcal{M}_{0}$ the intrinsic chirp mass, and $\Theta$ is the orientation function:

describing the detector’s responses to the different GW polarizations. Obviously, $\Theta$ depends only on the sky position and relative orientation of the source to the detector ($\theta,\phi,i,\psi$) $\theta$ and $\phi$ correspond to the usual spherical coordinates that describe the direction to the source, while $i$ and $\psi$ give the source’s orientation with respect to the detector (see Finn 1996 for a detailed discussion)., which are uncorrelated and uniformly distributed. A reasonable approximation of the probability distribution of $\Theta$ is given by (Finn 1996)

Combining equation (3) with equation (4), the S/N can be written as (Finn 1996; Taylor & Gair 2012)

is the detector’s characteristic distance parameter, with

is the dimensionless function reflecting the overlap between the GW signal generated by the inspiral stage and the detector’s effective bandwidth: $\zeta(f_{\rm max})$ is unity if $2f_{\rm max}$ is larger than the upper bound frequency of the detector’s bandwidth (i.e., the GW signal from the inspiral stage completely covers the detector’s effective bandwidth), and is less than unity if the inspiral terminates within the detector’s bandwidth.

The argument $f_{\rm max}$ is the redshifted orbital frequency at which the quadrupolar formula is no longer applicable (the binary finishes the inspiral and starts to merge). It is plausible to choose the entering of the innermost circular orbital (ICO) as the end of the inspiral stage. For binaries with equal-mass, this can be described as (Taylor & Gair 2012)

where $M$ is the total mass of the binary. For binaries with unequal mass, $f_{\text{ICO}}$ also depends on the mass asymmetry. In our simulation, we ignore this small correction, and make exclusive use of equation (13). We calculate $\zeta(f_{\rm max})$ for the typical total mass in our source sample ($M=10M_{\odot}$) and find it to be close to unity ($\sim 0.98$). Hence for simplicity, we adopt $\zeta(f_{\rm max})=1$ in the following calculations.

The distribution of the GW event rate in the observer’s frame with $z,\mathcal{M}_{0}$ and $\rho$ is given by (Finn 1996)

where $dV^{c}$ is the differential comoving volume and the factor $1/(1+z)$ accounts for the time dilation. $R_{\rm mrg}(\mathcal{M}_{0};z)$ is the intrinsic merger rate density with respect to the chirp mass $\mathcal{M}_{0}$ at redshift $z$. Our model to estimate this density is presented in Appendix B. The distribution $P_{\rho}(\rho|z,\mathcal{M}_{0})$ can be calculated by combining equations (8) and (9):

where $\Theta_{\rho}$ is rearranged from equation (9):

By marginalizing over $\mathcal{M}_{0}$ in equation (14), we can obtain the differential GW event rate with S/N $\rho$ at redshift $z$:

The GW event rate for a particular detector of threshold $\rho_{0}$ is given by

and the corresponding differential rate is

3 Lensing statistics

In the context of lensing statistics, the most important parameter is the so-called optical depth, or the differential lensing probability (e.g. Turner et al. 1984; Chae 2003; Huterer et al. 2005):

which describes the differential probability for a given source with S/N $\rho$ at redshift $z_{s}$ to be lensed.

The first integral takes into account the comoving volume between the observer and the source. It is required for calculating the total number of lensing galaxies.

The second integral gives the number density of lensing galaxies in comoving volume, where $\Psi(\sigma_{v})$ is the velocity distribution function of lensing galaxies. In the context of lensing statistics, the modified Schechter function (Choi et al. 2007)

is often used to fit the velocity distribution function, where $(\phi_{\ast},\sigma_{\ast},\alpha,\beta)$=$(8.0\times 10^{-3}h^{3}~\text{Mpc}^{-3},161\text{ km s}^{-1},2.32,2.67)$.

The third integral is over the distribution $p_{q}(q)$ of the projected axis ratio $q$. We adopt a Gaussian distribution to describe $p_{q}(q)$, with a mean of $0.7$, and standard deviation of $0.16$. The distribution is truncated at $q=0.2$ and $1.0$ . This is consistent with the observations (Jorgensen et al. 1995; Sheth et al. 2003).

The fourth integral is two-dimensional, where $p_{\bm{\gamma}}(\gamma,\theta_{\bm{\gamma}})$ denotes the distribution of external shear $\bm{\gamma}$. Following Huterer et al. 2005, we assume the amplitude $\gamma$ follows a log-normal distribution with mean $\ln 0.05$ and standard deviation $0.2$ (note: the mean and standard deviation are not the values for $\gamma$ itself, but of the underlying normal distribution it is derived from). The direction $\theta_{\bm{\gamma}}$ is assumed to be random.

The bias factor $B(\rho;z_{s})$ describes an enhancement of the representation of events due to the magnification caused by the lens (magnification bias): In gravitational lensing of GWs, the amplification in S/N is $\sqrt{\mu}$, since we directly observe the waveform instead of intensity.

Note that $\int_{0}^{\infty}d\mu~p_{\mu}(\mu)=1$ is required, in order to combine the bias factor into optical depth naturally.

The choice of the magnification factor $\mu$ is demanded for multiple-image systems in calculation of the bias factor. For double (two-image) lenses, we adopt the magnification factor of the fainter image as $\mu$, so that both images are magnified above the threshold. For quadruple (four-image) lenses, we adopt the magnification factor of the third brightest image, hence at least three images are magnified above the threshold. Also, this choice ensures the detection of the first lensed image to arrive, since the third brightest image is generally expected to arrive first (Oguri & Marshall 2010).

The two parts separated by square brackets are independent, and can be integrated separately.

As for the problem of determining the region of cross-section, we handle double and quadruple lenses separately. The naked cusp lenses are ignored in this paper, since they seldom happen at galaxy-scale lenses (Oguri & Marshall 2010). This treatment gives us the fraction of quadruple lenses. For double lenses, the condition that the fainter image is magnified above threshold $\rho_{0}$ is taken to define the region of cross-section. This treatment guarantees the theoretical detectability of multiple images.

4 Expected lensing rate

The expected lensing rate for a particular detector of threshold $\rho_{0}$ can be calculated as

In theory, by substituting equations (17) and (23) into equation (24), we can obtain the lensing rate as a function of threshold $\rho_{0}$. In practice, it is numerically more friendly if some rearrangements or reductions are undertaken. Hence we introduce a more practical form for calculation of the expected lensing rate:

is the integral which combines the velocity distribution function $\Psi(\sigma_{v})$ and the angular Einstein radius $\theta_{E}$.

For the differential rate, it is convenient to write the function as

Results

In this section, we present our prediction of the strongly lensed GW event rate (Section 3.1). We take aLIGO operating at its design sensitivity and the ET utilising its ‘xylophone’ configuration as illustrations. We calculate the lensing rate as a function of the characteristic distance $R_{0}$ (see below) to give a more general prediction for ground-based detectors with arbitrary sensitivity. In Section 3.2, we illustrate the probability distribution of lensing time delays to assess the detectability of multiple images in a finite duration.

The lensing rate is strongly dependent on the estimate of the GW event rate, which, in turn, depends on the estimate of the merger rate of stellar binary black holes. As an illustration, we use a simple recipe analogous to that in Cao et al. 2017 to compute the merger rate (see Appendix B for further details). Fig. 1 shows our results on the merger rate density distribution as a function of cosmic time (redshift). The two different lines represent estimates obtained by using observationally determined star formation rate (SFR) functions from Strolger et al. 2004 (red solid line) and from Madau & Dickinson 2014 (black dashed line). The discrepancy between these two results is noticeable at high redshift. This contradiction accounts for all the disparities in the following results. Despite the simplicity of our model, our results are comparable to those estimated through more sophisticated population synthesis models (see the comparison in Cao et al. 2017).

The GW detector’s sensitivity is described by a characteristic distance $R_{0}$, which depends only on the detector’s noise power spectral density $S_{n}(f)$ (see equation 10). Generally speaking, the larger $R_{0}$, the farther a detector can observe. For aLIGO, we use the data ‘ZERO_DET_high_P.txt’ from Shoemaker 2010 as the $S_{n}(f)$ for each interferometer operating at the design sensitivity. Since aLIGO consists of two interferometers with equal configurations and closely parallel orientations (one at Hanford, WA, and the other at Livingston, LA), we can treat aLIGO as a whole with the characteristic distance $\sqrt{2}$ times larger than that of each signal interferometer (Finn 1996). For the ET, which uses 3rd-generation technology and the ‘xylophone’ configuration, we adopt $R_{0}=1591$ Mpc (Taylor & Gair 2012).

Once the intrinsic merger rate of stellar binary black holes and the characteristic distance of the detector are determined, we can obtain the unlensed and lensed GW event rates through equations (18) and (25), respectively. The threshold $\rho_{0}$ is set to be eight, that means a signal is identified as detected when its S/N is above eight. Table 1 summarizes our results for the unlensed and lensed GW event rates in various detectors.

We predict that the unlensed GW event rates are $\sim 10^{3}\text{ yr}^{-1}$ for aLIGO at its design sensitivity ($R_{0}=155.4$ Mpc) and $\sim 10^{5}\text{ yr}^{-1}$ for ET ($R_{0}=1591$ Mpc). For comparison, we use the same strategy to compute the GW event rate at aLIGO’s O2 run ($R_{0}=63.7$ Mpc, $\rho_{0}=13$, $\zeta=0.4\sim 1.0$) The threshold $\rho_{0}$ is set according to GW170104, which has the lowest S/N among aLIGO’s O2 detections. The lower limit of $\zeta$ factor (see equation 12) is calculated using the total mass of GW170814 which has the largest mass among aLIGO’s O2 detections. and obtain $15\sim 75\text{ yr}^{-1}$, which is consistent with the current aLIGO detection rate. The improved sensitivity and the lower threshold of S/N account for the much higher expected source rate at aLIGO’s design sensitivity compared with that of the O2 run.

Based on these estimates of the GW event rate, we find that gravitational lensing of GWs is promising for both aLIGO at its design sensitivity and the proposed ET. More specifically, when the SFR function from Strolger et al. 2004 is adopted, both detectors have the largest expected numbers of lensed events (aLIGO $\sim 1\text{ yr}^{-1}$ and ET $\sim 80\text{ yr}^{-1}$). For the SFR function adopted from Madau & Dickinson 2014, the number in ET declines dramatically to $\sim 40\text{ yr}^{-1}$ due to the lower source rate expected at high redshift. The number in aLIGO drops only slightly and is still close to $1\text{ yr}^{-1}$, since aLIGO is insensitive to the event rate at high redshift. Furthermore, we compute the fraction of quadruple lenses in each survey. Our calculation indicates that the quadruple fraction is approximately $30$ per cent for aLIGO events and $6$ per cent for ET events. The higher quadruple fraction in aLIGO corresponds to the larger magnification bias.

We also calculate the most probable redshifts of the lensed sources (from $16$ per cent to $84$ per cent) and find it ranges from $\sim 1.1$ to $\sim 2.7$ for aLIGO events and from $\sim 1.5$ to $\sim 3.7$ for ET events based on the SFR function from Madau & Dickinson 2014. If the SFR function from Strolger et al. 2004 is adopted, the redshifts are slightly higher due to the higher estimates of source rates at high redshift (see Fig. 1), ranging from $\sim 1.2$ to $\sim 3.3$ for aLIGO events and from $\sim 1.8$ to $\sim 5.7$ for ET events.

We demonstrate the rate of lensed GW events as a function of the characteristic distance in Fig. 4. The notation for the different lines is the same as above. As expected, the larger $R_{0}$, the larger number of lensed events a detector can observe. This result indicates that any detectors more sensitive to aLIGO are expected to observe several strongly lensed events per year. The declining tendency of the quadruple fraction in the bottom panel is again due to the decrease in the magnification bias.

2 Distribution of time delays

A prediction for the time delay distribution is required in order to assess the detectability of multiple images during a finite duration GW survey. We achieve this goal through a semi-analytic technique based on Monte Carlo sampling (see Mao 1992 for a similar calculation for gamma-ray bursts).

The specific procedure is as follows. First, we randomly generate a sample of $10^{7}$ lens systems at a given source redshift. The lens objects are considered to be uniformly distributed on the sky, and the lens properties are distributed as described in Section 2.3. Then, we solve each lens system to see if it has multiple images, and for those with multiple images, we calculate their time delays through equation (36). By grouping these lens systems according to their time delays, we obtain the distribution of time delays. Since we do not set a threshold of S/N in this calculation, the distribution derived here considers all the lens systems satisfying the lens properties described in Section 2.3, not just those observable by a particular survey.

Fig. 5 shows the cumulative distribution function of the time delay for four representative source redshifts, $z_{s}=0.5$, $1.5$, $3.5$, and $10.5$, respectively. For double lenses (top left-hand panel) with a typical source redshift ($z_{s}=1.5$), $90$ per cent of the systems have time delays less than $\sim 1$ month. Even for the systems with a high source redshift ($z_{s}=10.5$), nearly $80$ per cent have time delays less than $1$ month. Almost all the systems have time delays less than $10$ months. This result indicates that the selection bias raised by the lensing time delay is insignificant, since the data-taking phases of GW detectors in the future will have durations well beyond most lens systems’ time delays. For example, the first and second runs (O1, O2) of aLIGO lasted for approximately $4$ months and $9$ months, respectively.

For quadruple lenses, we calculate time delays for three independent image pairs, in order of the arrival time: between the first and the second images [top right-hand panel; quad(12)], between the first and the third images [bottom left-hand panel; quad(13)], and between the first and the fourth images [bottom right-hand panel; quad(14)]. The result shows that for a typical time delay ($z_{s}=1.5$), $90$ percent of the systems have time delays between the first and the last images shorter than $\sim 0.4$ month, which implies missing any images due to the finite observation duration is unlikely. It is worth pointing out that the time delays between image pairs in quadruple lenses are typically shorter than those in double lenses. This feature implies a possible bias with quadruple lenses being over-represented in a finite survey.

Summary and discussion

In this paper, we have investigated the statistical properties of the strong gravitational lensing of GWs from stellar binary black hole coalescences in the context of ground-based detectors. By taking more realistic lens and source properties into account, we make a prediction for the rate of lensed GW events. Moreover, we calculate the probability distribution of lensing time delays to assess the selection bias due to the finite duration of a survey. Our main results can be summarized as follows.

We predict that aLIGO operating at its design sensitivity is expected to detect several lensed GW events (approximately $1$ event per year). The ET prediction is much higher (approximately $40\sim 80$ events per year) due to its much-higher sensitivity. The results are dominated by double lenses, with an expected quadruple fraction of $\sim 30$ per cent for aLIGO events and $\sim 6$ per cent for ET events. According to the SFR function from Madau & Dickinson 2014, the most probable redshifts of the lensed GW sources range from $\sim 1.1$ to $\sim 2.7$ for aLIGO events and from $\sim 1.5$ to $\sim 3.7$ for ET events. We emphasize the strong dependence between the predicted lensing rate and the source rate. This dependence leaves space for further improvement of the lensing rate prediction.

Specifically, the estimate of the merger rate density is calibrated to the current observations of stellar binary black hole GW sources by aLIGO and VIRGO, i.e. a mean rate density of $\sim 103~{\rm Gpc^{-3}\,yr^{-1}}$ in the local Universe. However, the current constraint on this mean rate density has a large uncertainty as shown in Abbott et al. 2017b, and it could range from $40$ to $213~{\rm Gyr^{-1}\,yr^{-1}}$ assuming a power-law distribution for the primary black hole masses. Considering this uncertainty, the strongly lensed GW event rate should be in the range of about a factor of $0.4$ to $2.1$ of the estimates listed above (a factor of $\sim 5$ uncertainty). Note also that the merger rate density estimated from the simple model presented in this paper seems to be smaller than some estimates by using binary population synthesis models (Cao et al. 2017, see discussion in). This may suggest that the strongly lensed GW event rate, especially for ET, may be even larger than the estimates obtained here.

Furthermore, we have developed a general calculation formalism of the lensing rate, not restricted to any specific detectors. The result indicates that any ground-based detectors more sensitive than aLIGO are anticipated to observe several lensed GW events per year (see Fig. 4). Detectors need not be a single instrument with unprecedented sensitivity such as ET but can be a network of interferometers such as aLIGO together with Virgo. As networks of GW detectors become routine in the near future (e.g. Abbott et al. 2016b, for a review of the commissioning roadmap), the detection of lensed GW events is expected even before ET becomes operational.

We have evaluated the chance of missing some images in a finite observation period, by examining the probability distribution of the lensing time delays. We find most lens systems involved in this study have time delays less than $\sim 1$ month (see Fig. 5). Since GW surveys in the future will have a duration much longer than a month, we expect the selection bias raised by the finite observation time should be small. Nevertheless, we emphasize that the time delays of quadruple lenses are systematically smaller than those of double lenses due to the smaller impact parameters in quadruple systems (the source is closer to the centre of the lens galaxy). This feature may result in a slightly higher fraction of quadruple lenses in a finite observation period.

In a real GW survey, other factors besides the finite duration may cause the absence of some images from detection, such as unexpected glitches in the detector, detector downtime for improvement and so on. Most of these effects can be eliminated by building up a network of several detectors (see Abbott et al. 2017d, for a treatment of the glitch in a real GW observation). This implies another advantage of joint detection in GW astronomy.

There are also some systematic errors due to the uncertainty of the velocity distribution function of lensing galaxies (equation 21). In this work, we adopt the modified Schechter function with parameters from Choi et al. 2007 based on the SDSS DR3 data, while other authors using different data bases obtain somewhat different parameters (see Montero-Dorta et al. 2017, for a recent comparison). Also, different strategies for sample selection and function modelling can affect the shape of the velocity distribution function (see e.g. Sohn et al. 2017). Furthermore, the velocity distribution function is expected to evolve with time at high redshift, though the details of this evolution are somewhat uncertain (see e.g. Bezanson et al. 2011). All these factors may introduce uncertainties to the results.

In this paper, we only consider GWs arising from stellar binary black hole coalescences. Although these sources as a whole dominate the high-band GW events, there are other types of double compact objects that can generate GWs, such as the inspiral of neutron star-neutron star or black hole-neutron star binaries. These sources are especially intriguing in multimessenger observations, as these systems are believed to be associated with kilonovae and can produce electromagnetic counterparts (see e.g. Metzger & Berger 2012, for a theoretical study and Abbott et al. 2017e for a real observation). In consideration of these systems’ enormous potential for physical and cosmological research (see e.g. Collett & Bacon 2017; Wei & Wu 2017), a further statistical study involving these systems is warranted.

Acknowledgements

We thank Richard Long for many constructive comments that improved the paper. This work was supported by the National Natural Science Foundation of China (Grant No. 11333003, 11390372 to SM; and 11690024 and 11390372 to YL). YL was also partly supported by the Strategic Priority Program of the Chinese Academy of Sciences (Grant No. XDB 23040100), and the National Key Program for Science and Technology Research and Development (Grant No. 2016YFA0400704).

References

Appendix A Lens Theory

In this appendix, we present further details of the lens theory based on the SIE model with external shear which is used in this paper. We refer the interested reader to Schneider et al. 1992; Schramm 1990; Kochanek 1991; Keeton & Kochanek 1998, and references therein for thorough discussions.

Bearing in mind the two-dimensional nature of the lensing calculation, we adopt $(x,y)$ and $(x_{s},y_{s})$ as position vectors in the lens plane (the “thin lens” approximation) and source plane, respectively. Using equations (1) and (2), the first derivatives of the lens potential $\phi(=\phi^{\rm SIE}+\phi^{\rm shear})$ are shown as (e.g. Keeton & Kochanek 1998)

where $b_{I}(q)=\lambda(q)\sqrt{q}$, $\psi=\sqrt{x^{2}+q^{2}y^{2}}$, and $e=\sqrt{1-q^{2}}$ is the eccentricity of lensing galaxies. The second derivatives are

is a system of nonlinear equations. Directly employing numerical calculation to solve equation (32) could be time-consuming. Introducing the polar coordinates $x=r\cos\alpha$, $y=r\sin\alpha$, (32) becomes

Now we can numerically solve the one-dimensional equation (33) to obtain the polar angle $\alpha$, then substitute $\alpha$ into equation (34) to obtain the radius $r$.

The time delay $\tau$ and magnification $\mu$ are given by (e.g. Schneider et al. 1992)

Appendix B Binary Black Hole Merger Rate

In this appendix, we describe the approach we use in computing the merger rate of binary black holes. We consider only stellar binary black boles formed from isolated massive binary stars in galaxies. These are the most promising sources of GWs that can be detected by ground-based GW surveys. The approach is similar to that presented in Cao et al. 2017 and Dvorkin et al. 2016.

Generally, the birth rate per unit volume of single black holes with mass $M_{\bullet}$ at the cosmic time $t$ is given by

Here $\phi(m_{\star})$ is the initial mass function of the star with the Chabrier initial mass function (Chabrier 2003) being adopted, $\dot{\psi}(Z;t)$ is the star formation rate (SFR) per unit volume with metallicity $Z$ at cosmic time $t$, and $\delta$ is Dirac-$\delta$ function. The relation between the mass of a stellar remnant black hole and the mass of its progenitor star is given by $M_{\bullet}=g(m_{\star},Z)$. We adopt the version obtained by Spera et al. 2015.

We assume that $\dot{\psi}(Z;t)$ can be separated into two independent functions, one is the total SFR function at redshift $z$ and the other is the metallicity distribution function at that redshift. For the total SFR function, we adopt the observationally determined functions from Madau & Dickinson 2014 and from Strolger et al. 2004. For the metallicity distribution function, we adopt the mean metallicity given by Belczynski et al. 2016.

Assuming that a fraction ($f_{\rm eff}$) of black holes exist as the primary components The primary component of a binary has mass ($M_{\bullet,1}$) larger than that ($M_{\bullet,2}$) of the secondary one. of binaries which can merge within the Hubble time, the merger rate density of stellar binary black holes is then given by

Here $P_{t}(t_{\rm d})$ is the probability distribution of the time delays $t_{\rm d}$ between the formation of stellar binary black holes and merger. We adopt the form $P_{t}(t_{\rm d})\propto t_{\rm d}^{-1}$ (O’Shaughnessy et al. 2010; Belczynski et al. 2016; Lamberts et al. 2016) and assume the minimum and maximum values of $t_{\rm d}$ are $50$ Myr and the Hubble time, respectively. $P_{q}(q)$ is the probability distribution of the mass ratio $q=M_{\bullet,2}/M_{\bullet,1}$ and is assumed to be independent of the black hole mass. We assume $P_{q}(q)\propto q$ over the range from $0.5$ to $1$, which seems to be consistent with binary population synthesis results (Belczynski et al. 2016; Cao et al. 2017). The parameter $f_{\rm eff}$ is determined by adopting the constraint on the mean detection rate of $103~{\rm Gpc}^{-3}\,{\rm yr}^{-1}$ given by the current aLIGO detection (Abbott et al. 2017b) to calibrate the merger rate density at $z\sim 0$ obtained from the model.

The merger rate density with respect to chirp mass $\mathcal{M}_{0}$ at redshift $z$ can be obtained as

where $\mathcal{M}_{q,M_{\bullet,1}}=q^{3/5}M_{\bullet,1}/(1+q)^{1/5}$ is the chirp mass of a black hole binary with primary black hole mass $M_{\bullet,1}$ and the mass ratio $q$.

By marginalizing over $\mathcal{M}_{0}$ in equation (40), we can obtain the merger rate density at redshift $z$: