How does the Earth's rotation affect predictions of gravitational wave strong lensing rates?
Lilan Yang Xuheng Ding Marek Biesiada Kai Liao Zong-Hong Zhu







School of Physics and Technology, Wuhan University, Wuhan 430072, China; yang_lilan@whu.edu.cn; zhuzh@whu.edu.cn

School of Physics and Technology, Wuhan University, Wuhan 430072, China; yang_lilan@whu.edu.cn; zhuzh@whu.edu.cn


Department of Astronomy, Beijing Normal University, Beijing, 100875, China; zhuzh@bnu.edu.cn
Department of Astrophysics and Cosmology, Institute of Physics, University of Silesia, 75 Pułku Piechoty 1, 41-500, Chorzów , Poland



School of Science, Wuhan University of Technology, Wuhan 430070, China;



School of Physics and Technology, Wuhan University, Wuhan 430072, China; yang_lilan@whu.edu.cn; zhuzh@whu.edu.cn
Department of Astronomy, Beijing Normal University, Beijing, 100875, China; zhuzh@bnu.edu.cn






The next generation of ground-based gravitational wave (GW) detectors, e.g. the Einstein Telescope, is expected to observe a significant number of  strongly lensed GW events as predicted in many previous papers.
However, all these works ignored the impact of the Earth's rotation on this prediction.
Multiple lensed images arrive at the Earth at different time, thus the ground-based detector has different responses to the lensed images due to different orientations of the detector relative to the GW source direction. Therefore the amplitudes of the GW signal from different images are modulated appropriately, in addition to the lensing magnification.
In order to assess this effect, we performed Monte Carlo simulations to calculate the event rate of lensed GW signals.

Our conclusion is that the Earth's rotation has a non-negligible impact on the event rate of lensed GW image. The updated event rates decrease by factors of ∼40%, ∼20%, ∼10%, for NS-NS, BH-NS, BH-BH systems respectively.









§ INTRODUCTION


After the first detection of gravitational wave (GW) event GW150914, which was produced by merging binary black holes (BBH)<cit.> and other subsequent detections <cit.>, GW astronomy came into being and a brand new window on the Universe was opened.
Especially the event GW170817 registered by the LIGO-Virgo gravitational wave detector network <cit.> accompanied by the detection of electromagnetic counterpart was a breakthrough that commenced a new era in the multi-messenger astronomy creating unique opportunities to deepen our understanding of the Universe. One may expect that LIGO-Virgo network will keep providing new detections and the next generation of ground-based detectors such as the Einstein Telescope (ET) operating with increased sensitivities will eventually yield 10^3-10^7 inspiral events per year reaching the redshift z=17 <cit.>.
Reaching so deep, one may expect that a significant number of such signals could be lensed by intervening galaxies.

Gravitational lensing of GWs has been extensively studied in many works since the pioneering paper <cit.>.
In particular, the effect of lensing on the parameter extraction of gravitational wave signals was discussed by <cit.>.
It was also suggested that cosmological parameters can be significantly constrained
using time delays measurements of strongly lensed GW events  <cit.>.
Moreover, strongly lensed GW signals can be used to test fundamental physics. For example, the speed of gravity can be tested with strongly lensed GW events accompanied by electromagnetic counterparts <cit.>.
Admittedly, constraints on the speed of GWs (expressed as bounds on the graviton mass or equivalently on its Compton wavelength) obtained by the LIGO Collaboration <cit.>  with un-lensed events are already very strigent. However, they strongly rely on the PPN waveform templates fitted to the data. The idea here is that if the Compton wavelength of the graviton was finite,
lower frequencies would propagate slower than higher frequencies, leading to the dispersion which would modify phasing of the coalescing signal. Similarly, the bound obatined from the coalescing NS-NS system  <cit.> was a conservative one assuming that GW and EM signals were emitted simultaneously and the observed delay was attributed solely to the difference of propagation speeds. Lensed GW-EM signals are free from such pre-assumptions. 
Therefore, the accurate prediction of lensed GW event rate becomes an important issue.


Lensed event rates for the ET detector were studied in <cit.>. The intrinsic merger rates of the whole class of double compact objects – DCOs thereafer – i.e. (NS-NS,BH-NS,BH-BH systems) located at different redshifts were taken from the StarTrack population synthesis evolutionary code <cit.>. Optical depth for lensing
was calculated using the singular isothermal sphere (SIS) lens model and Schechter-like velocity dispersion distribution in the population of lenses according to <cit.>.
The general conclusion was that the ET would register about 50 – 100 strongly lensed inspiral events per year. These would be dominated by BH-BH events contributing 91% – 95%  to the total rate, depending on details of evolutionary scenarios considered. Recently, <cit.> extended these predictions to more realistic lens properties allowing for the ellipticity of the lens (and thus the quadruple lensed images), the lens  environment (modeled as an external shear), and magnification bias.
Meanwhile, the prediction of the lensing rate for Advanced LIGO was revisited by <cit.>, where it had been noticed that the detection is correlated with source position and detector's orientation.
Comparing to ground based detectors,  lensing rate for the space-borne detector LISA were also discussed in  <cit.>.


All these works done so far, ignored the impact of the Earth's rotation. The orientation between the source and the detector was considered as  fixed. However, since different lensed signals travel along different paths and probe different depths of the potential of the lens (Shapiro effect), they actually arrive at the Earth at different times. The typical time delay of multiple signals varies from days to hundred of days. After a time delay, the rotation of the Earth changes the orientation of the ground-based detector with respect to the direction of the GW source.  As will see, this will affect the strength of the observed signal, in addition to the relative magnification due to lensing. If not taken into account, it would bias the predictions.

In this paper, we are filling this gap by accounting for the Earth rotation in estimating the rates of lensed GW signals. We base our forecasts on the Monte Carlo simulation. In section <ref>, we describe our methodology. In section  <ref>, we present and discuss the results. Finally, the conclusions are summarized in section <ref>.



§ METHODOLOGY


Our predictions are focused on the ET, which will consist of three nested detectors placed underground at depth of 100 - 200 m, arranged in a triangular pattern. Initial design assumed that each detector would be built from a single interferometer, where the high power needed to achieve good high-frequency performance compromises the low-frequency performance. The next step in the ET design is the so called “xylophone configuration”, where each detector is split into two interferometers, one specialised for detecting low-frequency gravitational waves and the other one for the high-frequency part. 

We briefly review the detection rate for unlensed events in section <ref> and lensing statics in section  <ref>. Then, we describe the details of our Monte Carlo simulation and update the event rate of lensed GW sources that could be  detected by the ET in section <ref>. 





 §.§ Detection rate of unlensed events


Since both the rudiments and the details of GW detection theory have been introduced many times <cit.>, we only recap the main points for clarity. The matched filtering is a standard technique applied in GW data-analysis to efficiently search for GW signals with known characteristics (templates) hidden in noisy data. Consequently, the strength of the signal is measured by the signal-to-noise ratio (SNR).
The optimal matched-filtering SNR of an inspiraling DCO system at the redshift z_s registered by a single detector is
 
    ρ = 8 Θr_0/d_L(z_s)( M_z/1.2 M_⊙)^5/6√(ζ(f_max)),

where: d_L is the luminosity distance to the inspiralling DCO, r_0 is detector's characteristic distance parameter. For two configurations considered here, i.e. the ET initial design and the  advanced “xylophone” design, r_0 = 1527 Mpc  and r_0 = 1918 Mpc,  respectively. M_z is the observed (redshifted)  chirp mass (i.e. M_z= M_0(1+z)) . We have assumed that M_0 are 1.2 M_⊙ for NS-NS, 3.2 M_⊙ for BH-NS and 6.7 M_⊙ for BH-BH systems respectively. According to <cit.>, these values represent average chirp mass for each category of DCO simulated by population synthesis. They were also used in <cit.>. ζ(f_max) is the dimensionless function reflecting the overlap between the GW signal and the detector's effective bandwidth. 
 Θ is the orientation factor determined by four angles,
according to:
 
    Θ = 2 [ F_+^2(1 + cos^2ι )^2 + 4 F_×^2 cos^2ι ]^1/2

where: F_+ = 1/2 (1 + cos^2θ) cos2ϕcos2 ψ - cosθsin2 ϕsin 2 ψ, and
F_× = 1/2 (1 + cos^2θ)  cos2ϕsin2 ψ + cosθsin2 ϕcos 2 ψ are called antenna patterns.
Two angles (θ ,ϕ) describe the direction to the DCO binary relative to the detector; another two angles ( ψ,ι) describe the binary orientation relative to the line-of-sight between it and the detector.
These four angles (θ ,ϕ, ψ,ι) are independent and one can assume that (cosθ, ϕ/π, ψ/π, cosι) are distributed uniformly over the range [-1, 1].
Due to rotation of the Earth, (θ ,ϕ) of a given source would change appropriately, while ( ψ,ι) would be the same.  Note that different values of (θ ,ϕ) could modify the value of Θ up to an order of magnitude, thus the effect of rotation of the Earth is non-trivial. For this reason, ρ of lensed images is not constant but is a function of time ρ(t). Single inspiral GW signal crosses the detector at one well defined moment t_0. SNR of such event can be considered as a single value ρ(t_0).
GW signal is detectable, if its SNR exceeds the threshold, i.e. ρ > ρ_0 = 8.
Then, the yearly detection rate of DCO sources can be expressed as:
 
    Ṅ(>ρ_0|z_s) = ∫_0^z_sd Ṅ(>ρ_0)/dz dz

where dṄ(>ρ_0)/dz_s is the yearly merging rate of DCO sources in the redshift interval [z_s, z_s+dz_s].
In one such redshift interval, the rate of DCO events is
 
    dṄ=4π(c/H_0)^3ṅ_0(z_s)/1+z_sr̃^2(z_s)/E(z_s)dz_s

where ṅ_0(z_s) denotes intrinsic inspiral rate at redshift z_s, r̃(z_s) is dimensionless comoving distance to the source and E(z_s) - dimensionless expansion rate of the Universe at redshift z_s. Concerning cosmological model, we assume flat ΛCDM model with Ω_m=0.3 and H_0 = 70  km s^-1 Mpc^-1 in order to comply with <cit.>. 
More importantly, we use the intrinsic merger rates ṅ_0(z_s) predicted by the population synthesis model (using StarTrack code) in <cit.>. We have taken the
data from the website http://www.syntheticuniverse.org, more specifically the so called “rest frame rates” in cosmological scenario. In this code  binary systems were evolved from ZAMS until the compact binary formation (after supernova explosions) making a number of assumptions about star formation rate, galaxy mass distribution, stellar populations, their metallicities and galaxy metallicity evolution with redshift . 
In order to investigate the
uncertainties of the chemical evolution of the Universe <cit.> employed two distinct scenarios for metallicity evolution with redshift, called “low-end” and “high-end”scenarios.
These assumptions result with different DCO merger rates predictions and consequently to distinct lensed GW event rates. 
Because the compact object formation depends critically on the physics of common envelope (CE) phase of evolution and on SN explosion mechanism and both of them are to some degree uncertain, <cit.> considered four scenarios: standard one and three of its modifications — Optimistic Common Envelope (OCE), delayed SN explosion and high BH kicks scenario. Some of the underlying assumptions could be constrained in light of existing GW detections <cit.>. However, the aim of this work is to evaluate the effect of the Earth's rotation, hence we focus only on the standard scenario to facilitate the comparison with previous works.




 §.§ Lensing statistics


Similarly as in <cit.>, we assume the SIS model of the lens. Of course, adopting more sophisticated lens models (such as SIE or power-law profile) could make the prediction more realistic and would enable to study lensing systems with quadruple images. In the recent paper <cit.>, whose authors used the SIE model, the results turned out to be consistent with the ones by SIS model (<cit.>), thus we conclude that the SIS model is sufficient for the purpose of our study, which is the prediction of the event rates.
From the physical point of view, mass distribution of lensing galaxies is the most relevant parameter for gravitational lensing.  
In the SIS model assumed here, the same information is conveyed by the  stellar velocity dispersion, which is much easier to assess.

Einstein radius θ_E gives the characteristic angular scale of lensing phenomenon. In the SIS model, it can be expressed  as θ_E=4π(σ/c)^2d_A(z_l,z_s)/d_A(z_s), where σ denotes velocity dispersion of lensing galaxy, d_A(z_s) and  d_A(z_l,z_s) are angular diameter distances to the source and between the source and the lens, respectively. Introducing the angle β, which is the angle between directions to the source location and to the center of the lens, one
concludes that if β<θ_E, two strong lensed images (brighter one I_+ and fainter one I_-) will appear at locations θ_±=θ_E±β. Introducing dimensionless quantities:  x=θ/θ_E and y=β/θ_E, the strong lensing condition becomes y<1 and positions of the images are x_±=1± y with magnifications μ_±=1/y±1. Gravitationally lensed GW signals, corresponding to these two images would have SNRs: ρ_±=ρ_intr._±√(μ_±)=ρ_intr._±√(1/y±1). Let us emphasize that the ρ_intr._- and ρ_intr._+ are now different due to the rotation of the Earth.

In previous works, only the I_- image was considered necessary for the detection in order to establish the lensed nature of the signal.
It was because its SNR ρ_- = ρ_intr._-√(1/y-1) was mistakenly assumed as being always lower than the SNR of I_+ i.e.  ρ_+ =ρ_intr._+√(1/y+1), which lead to a misconception that detecting the I_- ensures the ability to detect I_+. However, due to  the rotation of the Earth during the time delay between images it need not be the case. Therefore, we should treat the I_- and I_+ images  separately and we propose that the lensed GW event registered by the ET requires both ρ_- and ρ_+ exceeding the threshold ρ_0 = 8. This means that √(1/y±1)ρ_intr._±>ρ_0 is required, which leads to the condition:

    y_±≤ y_±,max = [ ( 8/ρ_intr.±)^2 ∓ 1 ]^-1

Thus, the elementary cross section for lensing reads:

    S_cr, ±(σ, z_l, z_s, ρ_±) = πθ_E^2 y_±,max^2  = 16 π^3 ( σ/c)^4 ( r̃_ls/r̃_s)^2 y_±, max^2 .

Finally, total optical depth τ describes probability that GW source at redshift z_s would be lensed  and detected by the ET.
Since the images I_+ and I_- are considered separately, we split τ into two terms, i.e. τ_± corresponding to I_±, respectively:

    τ_±(z_s, ρ_±) = 1/4 π∫_0^z_s  dz_l  ∫^∞_0   d σ 
    4 π( c/H_0)^3 r̃_l^2/E(z_l) S_cr, ±(σ, z_l, z_s, ρ_±) d n/d σ

where z_l is the redshift of lens (we assume that lenses are distributed homogeneously in redshift). Moreover, 
we model the velocity dispersion distribution in the population of lensing galaxies as a modified Schechter function d n/d σ = n_*( σ/σ_*)^αexp( - ( σ/σ_*)^β)β/Γ (α/β) 1/σ, with the parameters n_*,σ_*,α and β  taken after <cit.>. This is the same assumption as made in <cit.> - we reproduced it for comparison of results. 
Admittedly, there exist more recent data for the velocity distributions of galaxies like e.g. <cit.>. However,  these fits were made for galaxies of all types, unlike  <cit.> who used early type galaxies only.




 §.§ Monte Carlo method




Our goal is to update the event rates of lensed GW signals, explicitly considering variation of the
detector's orientation due to Earth's rotation, which will modulate the SNRs of I_-and I_+. This can hardly be done analytically, hence we need to perform the appropriate Monte Carlo simulation.
We describe the details of the simulation in this section.

We build up a mock universe by creating a large sample of DCO events to represent the overall DCO events, and randomly generate the values of their key parameters. For a DCO event at redshift z_s, its ρ at a specific moment can be randomly assigned with Eq.(<ref>). This is done by randomly assigning the orientation factor Θ based on sampling the four angles (θ ,ϕ, ψ,ι)[(cosθ, ϕ/π, ψ/π, cosι) are distributed uniformly in the  range [-1, 1].]. Let us note, that alternatively to the Monte Carlo method, the common practice is to define the Θ by taking its averaged numerical probability density as (see formula 3.11 in <cit.>):

    P_Θ(Θ)    =    5 Θ (4 - Θ)^3 /256,      if   
    0< Θ < 4  
    
    P_Θ(Θ)    =    0,      otherwise.


It is instructive to recall the origin of the Eq(<ref>). First,  <cit.> numerically estimated P(Θ^2) using Monte Carlo simulation – see Table 1
in <cit.>, where cumulative probability distribution of Θ^2 is reported. Later, 
<cit.> noticed that “to an excelent approximation” P_Θ(Θ) can be described by the Eq(<ref>). Since then this equation has been used abundantly due to its simple analytical form. 
In Figure <ref> we compare our random sampling of Θ based on 10^7 sample points (the same number as used by <cit.>) to this numerical probability density. 
One can see, that the numerical probability distribution Eq(<ref>) is not an excellent approximation.  

Identifying lensed DCO system by the ET requires that both lensed images  I_- and I_+ are detected.
We break this identification in two steps.
We start with answering how many lensed DCO events have I_- available to the ET and not considering I_+ for a moment, and the corresponding optical depth is τ_-(z_s,ρ_intr._-) (see Eq(<ref>)).
Knowing this probability, the total number of I_- detected  can be sampled by Monte Carlo simulation by accumulating the events which meet the requirement though all the space, i.e., ∫_0^z_max∫_0^∞ dρτ(z_s,ρ) ∂Ṅ/∂ z∂ρ  dz.
Let us remark that quantitatively this number should be equivalent to numerical values reported in previous papers <cit.>.

Once the set of “I_- detected” events is sampled, we can randomly assign the position y for each of them[The square of source position, i.e., y^2 follows a uniform distribution between (0, y_max^2).]. The SNR for corresponding I_+ images can thus be calculated as ρ_+ = ρ_intr._+√(1/y+1). As mentioned above, ρ_intr._+ is different from ρ_intr._- due to time delay. To take this difference into account, we keep the (ψ,ι) values unchanged and re-generate the values of (θ ,ϕ) to derive the ρ_intr._+.  This re-generating depends on the location of the ET on Earth and on the time delay between  “I_-” and  “I_+”.

In order to re-generate the values of (θ_+, ϕ_+) at t_+, knowing their initial values (θ_-, ϕ_-) at t_-, where t_+ and t_- are the times of arrival for I_+ and I_-, respectively[Let us remind that t_+ = t_- - Δ t, where Δ t is lensing time delay.], we proceed in the following way. The angles (θ, ϕ) are defined in a local detector's coordinate system, see Figure <ref>. For simplicity we assume that it is a cardinal coordinate system of the ET located at latitude 45^∘ N, longitude 0^∘ <cit.>. Then we transform the coordinate of the GW signal to the center of the Earth with z axis coincident with the Earth's rotation axis. 
In this new non-rotating coordinate system, the direction of the GW
signal is fixed [We neglect the orbital motion of the Earth
and corresponding change
of direction to the GW source.].
The angles (β, α) are used to describe the location of the ET: 
β=const.=45^∘, while α = α_0 + Ω_r t, where Ω_r = 2 π/T is the rotational angular velocity of the Earth. 



Consequently, 
α_+ = α_- - Ω_r Δ t,
and knowing the time delay Δ t, we can transform the coordinate of the
GW signal back to the detector's frame. Furthermore, we assume the time delay Δ t follows a uniform distribution as [0, 24hr]. Of course what matters here is Δ t  24hr.

Knowing ρ_intr._+, we obtain the optical depth τ_+(z_s,ρ_intr._+) for “I_+" which enables one to count the number of events having ρ_+>ρ_0 = 8 across all the redshift bins. Finally, we count the events for which both the I_- and I_+ could be detected by the ET.
Note that our procedure of considering the lensed image of I_- first  and then I_+ does not mean that the I_- arrives earlier than I_+. Actually, the I_+ arrives first. However, the identification of lensed GW requires that both I_- and I_+ are detected. In this work, we consider detection of I_- as the prior condition.




We performed 10^5 realizations of Monte Carlo simulation for all DCO systems and averaged the results.










§ RESULTS AND DISCUSSION

Table <ref>
shows the  expected yearly rates of lensed GW, based on the standard DCO scenario. Cases when only I_- image is detected and when both I_- and I_+ are detected by the ET, are shown.
One can see that the event rate of “I_-& I_+" is smaller than the event rate of “only I_-". This is the result of rotation of the Earth affecting the prediction of lensed events rates.
Table <ref> demonstrates that, compared to previous works, our updated event rates decrease by factors of ∼40%, ∼20%, ∼10%, for NS-NS, BH-NS and BH-BH systems, respectively. 
The reason why different DCO systems are affected in a distinct way can be understood from the distribution of the intrinsic SNR for different types of binaries (see Fig. 1. of <cit.> ). In the case of NS-NS systems, detectable lensed events (i.e. with ρ >8) correspond to the upper tail of the distribution, while the events with ρ < 8 dominate in NS-NS systems. Therefore, these systems are more affected by the Earth's rotation than BH-BH systems which are dominated by ρ > 8 cases. Threshold value ρ=8 splits the distribution of BH-NS systems into two approximately equal parts, which explains the intermediate value of their lensing rate reduction due to rotation of the Earth.
However, since the event rate is dominated by BH-BH systems, total GW lensing rate decreases by a factor of ∼10%. Despite the above described reduction of yearly rates, this means that if the I_- image of lensed GW event is detected, then the I_+ image is very likely to have been detected as well.
One can expect that since the effect of the Earth's rotation has also been neglected in <cit.>, their prediction of the total lensed event rate are overestimated by ∼10%.
Figure <ref> shows the relationship between yearly lensed event rate of “I_-& I_+" and the source redshift.
This relationship looks very similar to previous results <cit.> which indicates that our Monte Carlo simulation achieved  results consistent with numerical calculations.


In this work, we conservatively considered the SIS model, although considering the SIE model (or the power-law profile) would be more realistic and would enable one to study the quadruple lensing systems. However, for the purpose of predicting the lensed GW event rate, SIS model is sufficient, considering fact that the quad fraction is only about ∼10% of all strong lensing systems <cit.>. Of course, the predictions taking into account the Earth's rotation would be more difficult and challenging for the SIE model since there would be three mutual time delays between images. However, the Monte Carlo simulation similar to used in this paper, would then be the only reasonable approach.

The results obtained in this work are applicable to the ground-based detectors only. They are not valid for the next generation space-based detectors, like LISA which will be placed on an Earth-trailing heliocentric orbit <cit.> with three satellites forming an equilateral triangle. Since the satellites follow their own orbits the triangle will rotate as well, but the rotation period would be about 1 year. This means that the rotation of the LISA detector would have less effect on the lensed GW detection rate than the rotation of the Earth has in the case of ground-based detectors.




§ CONCLUSION

In this work, we used the Monte Carlo simulation technique to evaluate the effect of the Earth's rotation on the prediction of GWs strong lensing rates for the third generation ground-based detector – the ET. All previous works concerning GW lensing  ignored this effect and assumed that detecting the fainter image I_-  guarantees that the I_+ image could be also be observed. This was a  wrong, unjustified assumption. Our results show that the rate of detecting GW signals from both I_- and I_+ images is less than the rate of detecting I_- only.
In particular,  for NS-NS systems, the event rate of both I_- and I_+ images is ∼40%  less than the rate of detecting only one of them. For BH-BH systems, the impact of the Earth's rotation is the smallest (∼10%).
Furthermore, we have shown that the Monte Carlo method provides reliable and more accurate results compared D to previous work <cit.> using numerical calculation.

Since the total event rate is dominated by BH-BH systems, the impact of the Earth's rotation on the total rate is at the level of ∼10%. Therefore, one should not worry much about this effect making  cosmological inference, such as determination of the Hubble constant or cosmic equation of state, using catalogs of inspiral events. One should keep in mind, however, that this effect is the highest (∼40%) for  NS-NS systems, which are accompanied by electromagnetic counterparts.
Therefore, the Earth's rotation should be taken into account in all considerations concerning lensed GW signals from coalescing NS binaries.













This work was supported by the National Natural Science Foundation of China under Grants Nos. 11633001 and 11373014, the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB23000000 and the Interdiscipline Research Funds of Beijing Normal University.
X. Ding acknowledges support by China Postdoctoral Science Foundation Funded Project (No. 2017M622501).
M.B. expresses his gratitude for hospitality of the Wuhan University where part of this work was done.
M.B. was supported by Foreign Talent
Introducing Project and Special Fund Support of Foreign Knowledge
Introducing Project in China.
K. Liao was supported by the National Natural Science Foundation of China (NSFC) No. 11603015.




