Time delay of timelike particles in gravitational lensing of Schwarzschild spacetime

TL;DR

This paper analyzes the time delays of timelike particles in Schwarzschild spacetime's gravitational lensing, deriving formulas for different scenarios and exploring implications for neutrino mass measurements and gravitational wave velocity constraints.

Contribution

It provides exact and approximate formulas for time delays of signals with arbitrary velocities in Schwarzschild spacetime and applies these to neutrinos and gravitational waves.

Findings

Time delay formulas depend on signal velocity and distances.

Neutrino time delays can help determine neutrino mass ordering.

Gravitational wave delay differences are too small for current measurement accuracy.

Abstract

Time delay in Schwarzschild spacetime for null and timelike signals with arbitrary velocity $v$ is studied. The total travel time $t_\mathrm{if}$ is evaluated both exactly and approximately in the weak field limit, with the result given as functions of signal velocity, source-lens and lens-observer distances, angular position of the source and lens mass. Two time delays, $\Delta t_v$ between signals with different velocities but coming from same side of the lens and $\Delta t_\mathrm{p}$ between signals from different sides of the lens, as well as the difference $\Delta t_{\mathrm{p}v}$ between two $\Delta t_\mathrm{p}$'s are calculated. These time delays are applied to the gravitational-lensed supernova neutrinos and gravitational waves (GW). It is shown that the $\Delta t_v$ between different mass eigenstates of supernova neutrinos can be related to the mass square difference of these eigenstates and therefore could potentially be used to discriminate neutrino mass orderings, while the difference $\Delta t_{\mathrm{p}v}$ between neutrino and optical signals can be correlated with the absolute mass of neutrinos. The formula for time delay in a general lens mass profile is derived and the result is applied to the singular isothermal sphere case. For GWs, it is found that the difference $\Delta t_{\mathrm{p}v}$ between GW and GRB can only reach $1.45\times 10^{-5}$ second for very large source distance ($2\times 10^4$ [Mpc]) and source angle (10 [as]) if $v_{GM}=(1-3\times 10^{-15})c$. This time difference is at least three order smaller than uncertainties in time measurement of the recently observed GW/GRB signals and thus calls for improvement if $\Delta t_{\mathrm{p}v}$ is to be used to further constrain the GW velocity.

Equations

$$\mathrm{d}{s}^2=A(r)\mathrm{d}{t}^2-B(r)\mathrm{d}{r}^2-C(r)\left(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\phi }^2\right)$$

$$\frac{\mathrm{d}t}{\mathrm{d}r}=\frac{E\sqrt{B(r)C(r)}}{\sqrt{A(r)}\sqrt{E^{2}C(r)-L^{2}A(r)-\kappa A(r)C(r)}},$$

$$A(r)\frac{\mathrm{d}t}{\mathrm{d}\lambda }=E,C(r)\frac{\mathrm{d}\phi }{\mathrm{d}\lambda }=L.$$

$$E=\frac{1}{\sqrt{1-v^2}},L=|\mathbf{p}\times \mathbf{r}|=\frac{v}{\sqrt{1-v^2}}b.$$

$$b=\frac{L}{Ev}$$

$$(bEv{)}^2=\frac{[E^2-\kappa A(r_0)]C(r_0)}{A(r_0)},$$

$$A(r)=1-\frac{2M}{r},B(r)={\left(1-\frac{2M}{r}\right)}^{-1},C(r)=r^2$$

$$\frac{\mathrm{d}t}{\mathrm{d}r}=\frac{E{r}^2}{(r-2M)\sqrt{E^2{r}^2-L^2(1-\frac{2M}{r})-r(r-2M)}},$$

$$b^2(E^2-1)=\frac{r_0^2\left[(E^2-1)r_0+2M\right]}{r_0-2M}.$$

$$\frac{\mathrm{d}t}{\mathrm{d}u}=\frac{E{r}_0^3}{u^2(2Mu-r_0)\sqrt{r_0\left(2Mu+\left(E^2-1\right)r_0-\frac{u^2(2Mu-r_0)[2M+(E^2-1)r_0]}{2M-r_0}\right)}}\equiv g(u,r_0,E),$$

$$t_{\mathrm{if}}={\int }_1^{\frac{r_0}{r_{\mathrm{i}}}}g(u,r_0,E)\mathrm{d}u+{\int }_1^{\frac{r_0}{r_{\mathrm{f}}}}g(u,r_0,E)\mathrm{d}u\equiv t_{\mathrm{i}}+t_{\mathrm{f}}.$$

$$g(u,r_0,E)=\frac{c_g}{u^2(u-u_3)\sqrt{(u-1)(u-u_1)(u-u_2)}},$$

$$u_{\begin{array}{c}1\\ 2\end{array}}=$$

$$u_3=$$

$$c_g=$$

$$t_{i/f}=c_g\sum _{j=0}^7{C}_j\left[f_j\left(\frac{r_0}{r_{i/f}}\right)-f_j(1)\right],$$

$$C_0$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_{\begin{array}{c}4\\ 5\end{array}}$$

$$C_{\begin{array}{c}6\\ 7\end{array}}$$

$$f_0(x)=$$

$$f_1(x)=$$

$$f_2(x)=$$

$$f_3(x)=$$

$$f_{\begin{array}{c}4\\ 5\end{array}}(x)=$$

$$f_{\begin{array}{c}6\\ 7\end{array}}(x)=$$

$$h(x)=\frac{\sqrt{u_1-1}+\sqrt{x-1}}{\sqrt{u_1-1}-\sqrt{x-1}}.$$

$$\underset{u\to 1}{\lim }{C}_0{f}_0(u)+C_1{f}_1(u)$$

$$t_{\mathrm{if}}$$

$$r_0=\frac{1}{3}\left(-\frac{2M}{E^2-1}+\frac{\frac{4M^2}{(E^2-1{)}^2}+3b^2}{f(E,b)}+f(E,b)\right),$$

$$f(E,b)=$$

$$\displaystyle+3\sqrt{3}\sqrt{\frac{16b^{2}M^{4}}{(E^{2}-1)^{3}}+\frac{(8-36E^{2}+27E^{4})b^{4}M^{2}}{(E^{2}-1)^{2}}-b^{6}}\bigg{\}}^{1/3}.$$

$$r_0=\frac{b^2}{\sqrt{3}{\left[b^2\left(-3\sqrt{3}M+\sqrt{27M^2-b^2}\right)\right]}^{1/3}}+\frac{{\left[b^2\left(-3\sqrt{3}M+\sqrt{27M^2-b^2}\right)\right]}^{1/3}}{\sqrt{3}}.$$

$$r_{\mathrm{i}},r_{\mathrm{f}}\gg r_0\gg M,$$

$$\beta =\theta -\frac{r_{\mathrm{i}}}{r_{\mathrm{i}}+r_{\mathrm{f}}}\varphi ,$$

Full text

Gravitational lensing; Time delay; Schwarzschild spacetime; Timelike particles

I Introduction

With the discovery of supernova neutrino (SNN) from SN1987A Hirata:1987hu ; Bionta:1987qt , the recent observation of gravitational wave (GW) signal Abbott:2016blz ; Abbott:2016nmj ; Abbott:2017oio ; TheLIGOScientific:2017qsa and more recent confirmation of neutrino emission from blazer TXS 0506+056 IceCube:2018cha ; IceCube:2018dnn , astronomy has certainly entered the multi-messenger era. In particular, the observation of the binary neutron star merger GW170817 and GRB 170817A TheLIGOScientific:2017qsa ; GBM:2017lvd ; Monitor:2017mdv is a simultaneous observation of the Gamma-ray burst (GRB) and GW signals. The 1.74 [s] time difference between the GRB and GW signals constrains the speed difference of gravity and light, the Equivalence Principle and the physical properties of the central engine of the GRB GBM:2017lvd . It can also put stringent constraints on the parameter space of general scalar-tensor Sakstein:2017xjx , vector-tensor gravity Baker:2017hug , dark energy Creminelli:2017sry ; Ezquiaga:2017ekz and dark matter models Boran:2017rdn . This time difference however, can originate from one or more of the three sources: (i) the generation region of the GW-GRB, (ii) the propagation path from source to local galaxy, and (iii) near the Galaxy until received by detectors. For signals from sources of high redshift such as GW170817, it is generally expected that the gravitational potential along propagation is more important than the potential at the originating or receiving sites.

If there exist large mass along the propagation path, then the signal will experience gravitational lensing (GL), regardless whether the signal is neutrino, GRB or GW, massive or massless. Many observables of the GL such as the apparent angles, the time delay between different images and the magnification can be used to deduce properties of the signal source, the signal itself and spacetime it went through Walsh:1979nx ; Aubourg:1993wb ; Alcock:2000ph ; Gaudi:2008zq ; Gould:2010bk ; Oguri:2010ns ; Treu:2010uj . The time delay between different images of the same kind of signal has a special advantage over the time delay between signals of different particles in the same event: the former do not suffer the uncertainty of emission time since these different images are from the same emission. Even if the emission times of this signal is not known exactly, the time delay produced during traveling are still valuable to deduce properties of the sources, the signal particle/wave or the spacetime transmitted.

However, in the computation of the time delay of timelike particles, the formulas for null ones are usually used barrow1987lensing ; Eiroa:2008ks ; Fan:2016swi ; Wei:2017emo ; Yang:2018bdf , i.e., the timelike nature of the massive particle/wave was not fully accounted. Although the neutrinos from supernova and GWs from mergers are usually relativistic, the time delay itself indeed is the difference between the total travel times, which are large quantities too. Therefore for high accuracy calculations, especially when the timelike particle is not that relativistic or the lensing is strong, the timelike nature of the particles shall be fully addressed.

To fulfill this purpose, in this paper we propose to study the time delays between signals with different velocities and between different images of same or different kind of signals, especially those with nonzero masses. For simplicity, we will concentrate on the Schwarzschild spacetime. Previously, the time delay of light has been extensively studied using different approaches in the weak and strong field limits in various spacetimes or gravities Richter:1982zz ; Edery:1997hu ; Bozza:2003cp ; Sereno:2004 ; Jacob:2008bw ; Bailey:2009me ; Eiroa:2013 ; Sahu:2013 ; Wang:2014yya ; Zhao:2016kft ; He:2016cya ; Zhao:2017cwk ; Zhao:2017jmv ; Deng:2017umx . For massive particles, the time delay of massive photon was investigated in Ref. Glicenstein:2017lrm to constrain photon mass. Ref. Baker:2016reh studied the difference between time delay of massive neutrino and that of massless GW to place bounds on total neutrino mass and some cosmological parameters. We emphasize that our work is different from theirs in a few ways. First, unlike Ref. Baker:2016reh ; Glicenstein:2017lrm which works only in the relativistic limit of the signal and weak field limit of the lens, we computed the exact total travel time that works for arbitrary velocity, i.e., velocities not very close to $c$. Secondly, even in the relativistic limit, the total time and time delay in these works are only to the order ${\cal O}(c-v)^{1}$ while our approximation formulas can works to higher orders. Thirdly, we have computed two time delays, $\Delta t_{v}$ due to velocity difference and $\Delta t_{\mathrm{p}}$ due to path difference, instead of only $\Delta t_{\mathrm{p}}$. These time delays and their differences are used to constrain the mass ordering and absolute mass of neutrinos and GW speed. Finally, the velocity correction to the time delay of signal with arbitrary velocity in general mass profiles is found.

The paper is organized as follows. In Sec. II, we set up the general framework for the calculation of the total travel time $t_{\mathrm{if}}$ of signal with arbitrary velocity in Schwarzschild spacetime. In Sec. III, the $t_{\mathrm{if}}$ is first evaluated exactly and the result is found as a combination of several elliptic functions. The same $t_{\mathrm{if}}$ was then computed in the weak field limit and a much simpler expression is found. Both the exact and approximate $t_{\mathrm{if}}$ are expressed as functions of signal velocity, source-lens distance, lens-observer distance, angular position of the source and lens mass. In Sec. IV, time delay $\Delta t_{v}$ between signals with different velocities but coming from same side of the lens, and time delay $\Delta t_{\mathrm{p}}$ between signals from different sides of the lens, as well as the difference $\Delta t_{\mathrm{p}v}$ between two $\Delta t_{\mathrm{p}}$ are found. These time delays are then applied to the cases of supernova neutrino and GW in Sec. V. It is shown that the $\Delta t_{v}$ might be related to the mass square difference between neutrino mass eigenstates and therefore be used to discriminate neutrino mass orderings, while the difference $\Delta t_{\mathrm{p}v}$ between neutrino and optical signals can be correlated with the absolute mass of neutrinos. For GW, the formula for time delay in a general gravitational potential is derived first. Then the difference $\Delta t_{\mathrm{p}v}$ between GW and GRB is evaluated. It is shown that even for distance as large as $2\times 10^{4}$ [Mpc] and $v_{\mathrm{GW}}=(1-10^{-15})c$, $\Delta t_{\mathrm{p}v}$ can only reach to the value of about $1.45\times 10^{-5}$ [s]. To further constrain $v_{\mathrm{GW}}$ therefore calls for the improvement in uncertainty of GW and GRB time measurements.

II Geodesic equations and the total travel time

We consider the time delay of signal particle/wave to be of different kinds and have different velocities. When passing by a gravitational center, they will experience different travel time. To calculate this travel time, we start from the general spherically symmetric spacetime with metric

[TABLE]

where $(t,~{}r,~{}\theta,~{}\phi)$ are the coordinates. Using the geodesic and normalization equations, we can obtain the following equation of motion in the equatorial plane for coordinate $t$ and $r$

[TABLE]

where $\kappa=0,1$ for null and timelike particles respectively. Here $E$ and $L$ are the first integrals of the geodesic equations of $t$ and $\phi$, satisfying

[TABLE]

They can be interpreted as the energy and orbital angular momentum of the particle per unit mass at infinity. For timelike particles, we have

[TABLE]

Here $v$ is the speed of the particle at infinity and $b$ is the impact parameter. For null particles, $E$ would approach infinity but the relation

[TABLE]

holds for both null and timelike particles.

In order to calculate the travel time, in principle we should integrate Eq. (2) from the source coordinate $r_{\mathrm{i}}$ to the observer coordinate $r_{\mathrm{f}}$. This is usually done by integrating from $r_{\mathrm{i}}$ to the closest radius $r_{0}$ first and then adding the same integral from $r_{0}$ to $r_{\mathrm{f}}$. The closest radius $r_{0}$ is defined as the maximal radial coordinate $r$ satisfying $\displaystyle\frac{\mathrm{d}r}{\mathrm{d}t}=0$. Using Eq. (2) and then substituting $L$ in Eq. (5), this is equivalent to

[TABLE]

from which $r_{0}$ can be formally solved in terms of $E$ and $b$.

To facilitate the relevant computation in Schwarzschild spacetime, we now substitute

[TABLE]

into Eq. (2) and (6), and find for timelike particles

[TABLE]

Note that we use the natural unit $G=c=1$ throughout the paper. The corresponding equations of null particles can be obtained by taking the infinite $E$ limit in the above equations.

For the purpose of later integration, it is convenient to make the change of variable in Eq. (8) from $r$ to $\displaystyle u=\frac{r_{0}}{r}$, so that this equation becomes

[TABLE]

where $L$ was replaced using Eq. (4) and for simplicity we used $g(u,r_{0},E)$ to denote the right hand side. Then the total travel time $t_{\mathrm{if}}$ from the source at $r_{\mathrm{i}}$ to the observer at $r_{\mathrm{f}}$ is given by the following integral

[TABLE]

Here the first and second integrals $t_{\mathrm{i}}$ and $t_{\mathrm{f}}$ are the times from $r_{\mathrm{i}}$ to $r_{0}$ and time from $r_{0}$ to $r_{\mathrm{f}}$ respectively. In next section, we will calculate this travel time using exact integration and the approximation method.

III Exact travel time and its approximation

III.1 Analytical integration

In order to integrate the Eq. (11), we first simplify the integrand $g(u,~{}r_{0},~{}E)$ to the form

[TABLE]

where $u_{1},\ u_{2},\ u_{3}$ are the roots of the denominator and $c_{g}$ is a coefficient, given by

[TABLE]

The time $t_{\mathrm{i}}$ and $t_{\mathrm{f}}$ in Eq. (11) then can be integrated analytically and the result is a linear combination of seven elliptic functions

[TABLE]

where $c_{g}$ is defined in Eq. (15) and the coefficients $C_{j}$’s are

[TABLE]

and the functions $f_{j}(u)$ are

[TABLE]

Here “$\mathrm{i}$” is the imaginary unit, $F,~{}E,~{}\Pi$ are respectively the elliptical integral of the first, second kind and the incomplete elliptic integral defined in Appendix A, and $h(x)$ is an axillary function defined as

[TABLE]

The second term inside the bracket in Eq. (16) contains the functions $f_{j}(u)$ evaluated at $u=1$. Evaluation of $f_{0}(1)$ and $f_{1}(1)$ however demands special care because they are separately divergent but their divergences cancel each other exactly and therefore their combination is still finite. For these two terms, we find

[TABLE]

Substituting Eq. (30) into (16) and then (11), then the total travel time is found as

[TABLE]

Eq. (31) is a function of five parameters $E$ (or $v$), $r_{\mathrm{i}}$, $r_{\mathrm{f}}$, $r_{0}$ and $M$, which once are known, the total travel time would immediately follow. From the definitions of $c_{i},~{}C_{i}$ and $F_{i}$ in Eqs. (13)-(15) and (17)-(28) and dimension counting, one can recognize that if all distance variables $r_{\mathrm{i}}$, $r_{\mathrm{f}}$ and $r_{0}$ are measured in the unit of $M$, then $t_{\mathrm{if}}$ would be linear to $M$, i.e., $t_{\mathrm{if}}=M\cdot q(E,r_{\mathrm{i}}/M,r_{\mathrm{f}}/M)$ for some function $q$. Therefore effectively, the dependence of $t_{\mathrm{if}}$ and the time delays that will be discussed later on $M$ is simple and we can concentrate on dependance on other parameters. Among all parameters $(E,r_{\mathrm{i}},r_{\mathrm{f}},r_{0},M)$, $r_{\mathrm{i}},~{}r_{\mathrm{f}}$ and $M$ are usually deducible using other astrophysical observations or theoretical tools. $E$ (or $v$) can be measured at the observatory. Therefore, there is only one last obstacle: the closest radius $r_{0}$, that is not practically known or easily measurable. Although in Eq. (9) we can solve $r_{0}$ in terms of $E$ and $b$, the impact parameter $b$ is not known explicitly either. One therefore has to find a way to further link $b$ to some measurable quantities in the GL setup. In the remaining part of this section, we will first solve Eq. (9) for $r_{0}$ and then show that $b$ can be tied to the apparent angle $\theta$ in the lens equation which is further solvable in terms of $r_{\mathrm{i}},~{}r_{\mathrm{f}}$ and the angular position $\beta$ of the source (see Fig. 1).

From Eq. (9), which is a cubic polynomial, $r_{0}$ can be solved as the only positive solution that is accessible by particles coming from and going back to infinity,

[TABLE]

where

[TABLE]

For null rays, this $r_{0}$ can be simplified to

[TABLE]

Eq. (32) is equivalent to Eq. (14) in Ref. Jia:2015zon while Eq. (34) agrees with Eq. (6.3.37) in Ref. Wald:1984 .

Now to connect $b$ to other measurable quantities, we have to use the lens equation. We will only consider this equation in the weak field limit, i.e.,

[TABLE]

although the lens equations in the strong field limit and their solutions are also known Jia:2015zon . The reason is that GL in the strong field limit for both the neutrinos and GW are beyond observational capability in the near future. Technically, the lens equation in the strong field limit are also much more involved to solve. The lens equation in the weak field limit is given by (see Fig. 1)

[TABLE]

where $\beta$ is the angular position of the source, $\theta$ is the apparent angle and $\phi$ is the deflect angle of the geodesic trajectory. Both $\theta$ and $\phi$ can be linked to the impact parameter $b$. For $\theta$, its relation with $b$ under the weak field limit is given by

[TABLE]

And for $\phi$ in the weak field limit, its value to the order ${\cal O}(b)^{-1}$ for particles with arbitrary $E$ or $v$ is Accioly:2002 ; Pang:2018

[TABLE]

Substituting Eqs. (37) and (38) into (36), one obtains a quadratic equation of $b$, whose solutions and the corresponding apparent angles of the images are

[TABLE]

The positive $b_{+}$ corresponds to the particle trajectory along the path on the same side of the lens-observer axis as the source, while the negative $b_{-}$ corresponds to the path on the other side. Their size satisfies $|b_{+}|\geq|b_{-}|$.

Substituting Eq. (39) back into Eq. (33) and further into Eq. (32), $r_{0}$ can be expressed as a function of parameters $E,~{}\beta,~{}r_{\mathrm{i}},~{}r_{\mathrm{f}}$ and $M$

[TABLE]

The explicit formula of Eq. (40) is elementary but too long to show here. It enables the computation of $t_{\mathrm{if}}$ in Eq. (31) in terms of measurables $E~{}(\mbox{or}~{}v),~{}\beta,~{}r_{\mathrm{i}},~{}r_{\mathrm{f}}$ and $M$.

Because there are two solutions of the impact parameter $b_{\pm}$ for one set of parameters $(E,~{}\beta,~{}r_{\mathrm{i}},~{}r_{\mathrm{f}},~{}M)$, there are two trajectories connecting the source and the observer and correspondingly two $t_{\mathrm{if}}$’s. In turn, there will be two basic types of time delay in the lensing of particles with different velocities: (1) the time delay between the total travel times of particles with different velocities along path on the same side of the lens, and (2) the time delay between total travel times of particles along paths on different sides of the lens. We will denote these two types time delay as $\Delta t_{v}$ and $\Delta t_{\mathrm{p}}$ respectively. Of course, if both the particle velocities and path sides are different, the time delay will be a mixed of these two types.

We emphasis that Eqs. (31) and (32) are exact formulas for all kinds of lensing including weak, strong or retro- lensings. More importantly, these results are valid for all particle velocity and therefore allow us to study the time delay in lensing of neutrinos and (potentially) massive gravitons. Eqs. (31) do have a drawback that it is expressed using complex elliptical functions which might hinder a simple and clear understanding of the physics, e.g., the effect of various parameters $E$ (or $v$), $\beta$, $r_{\mathrm{i}}$ and $r_{\mathrm{f}}$, on the time delay. Therefore, it is desirable to consider an approximation of these results.

III.2 Approximation in weak field limit

In the derivation of Eq. (39), we have used the weak field limit (35). In this subsection, we will extend the application of this limit to the integration of Eq. (10) and to the solution (32). The key is to note that in this limit, $r_{0}$ is much larger than $M$. Therefore, making an asymptotic expansion of small quantity $M/r_{0}$ in Eq. (10), it is transformed into

[TABLE]

Integrating this using the limits in Eq. (11) and dropping the terms of order ${\cal O}(M/r_{0})$ and higher, the total travel time in the weak field limit becomes

[TABLE]

Beside $M/r_{0}$, there exists another small ratio $r_{0}/r_{\mathrm{i/f}}$ in the weak field limit. Expanding Eq. (42) to the order ${\cal O}(M/r_{0})^{0}$ and ${\cal O}(r_{0}/r_{i/f})^{2}$ respectively, one finds

[TABLE]

Eq. (III.2) and other previous formulas for the total travel time were for particles with energy $E$ of unit mass at infinite radius. Parameter $E$ however is not very convenient for comparing the total travel time for particles with different rest masses. Therefore we replace $E$ by velocity $v$ at infinity using Eq. (4) in various formulas. In particular, Eq. (42) becomes

[TABLE]

The first term is of geometrical origin and represents the propagation time for particle with general velocity $v$ along the bent path. The second term represent the effect of the general relativistic gravitational potential to the total travel time. When $v=c$, this becomes the well-known total travel time for null particles Hartle:2003yu . For Eq. (III.2), after substitution of $v$, the total travel time is approximated by

[TABLE]

The first term is the time cost for travel if the spacetime is Minkovski and the second term is the correction because the bending of the geodesic trajectory causes extra distance, and both these two terms originate from the geometric propagation time term in Eq. (44). The third term in Eq. (45), when setting $v=c$, is half of the conventional Shapiro time delay for a returning light signal. The last term is the high order term from the general relativistic potential term in Eq. (44). Note that although we have used the weak field limit (35) to various orders in the derivation of Eqs. (44) and (45), no assumption on velocity $v$ was used and therefore Eqs. (44) and (45) should be valid for any velocity.

The $r_{0}$ in Eq. (III.2), which was given by Eq. (40), can also be further simplified in the weak field limit. The key is that this limit not only implies $r_{0}\gg M$ but also $|b_{\pm}|\gg M$. Therefore expanding the right side of solution (40) for large $b$, one obtains to the ${\cal O}(b_{\pm})^{0}$ order the following result

[TABLE]

which after substituting Eq. (39) for $b_{\pm}$ and $v$ for $E$ becomes

[TABLE]

Substituting this into Eq. (44) or (45), one then can obtain the total time for each set of parameters in a much simpler way than Eqs. (31) and (40). To check the correctness of these results, we have numerically calculated the $t_{\mathrm{if}}$’s using both Eqs. (31) and (45) with variables within the parameter ranges used in Sec. IV and V and excellent agreement was found.

IV Two types of time delay

As we pointed out in the previous section, there will be two basic types of time delay if signals with different velocities are lensed: the time delay $\Delta t_{v}$ between particles of different velocities, and the time delay $\Delta t_{\mathrm{p}}$ between particles traveled along different paths. In principle we can use both the exact result Eq. (31) and the weak field limit result Eq. (44) to find all these time delays. However, for the purpose of later usage and more intuitive understanding of relevant results, in this section we will derive perturbative results for these two types of time delay starting from Eq. (45). Note that the total travel time $t_{\mathrm{if,w}}$ is dependent on variables $v,~{}r_{\mathrm{i}},~{}r_{\mathrm{f}},~{}M$ explicitly and on $\beta$ and the path choice ($\pm$ sign) implicitly through $r_{0\mathrm{w}\pm}$ in Eq. (47). In other words,

[TABLE]

In the computations below, for the simplicity of the notation, we will only keep necessary variables and suppress the rest.

IV.1 Time delay $\Delta t_{v}$

We consider the situation that two signals with velocities $v$ and $v^{\prime}~{}(v^{\prime}>v)$ traveling on same side of the lens first. If $v$ and $v^{\prime}$ are very different, then the time delay should be evaluated directly using Eq. (45)

[TABLE]

Given that for all practical $\beta$, the weak field limit (35) implies that the main contribution to $\Delta t_{v}$ comes from the first term in Eq. (45). In other words, we should roughly have

[TABLE]

Since the dependance of $\Delta t_{v}$ on $\beta$ is hidden in $r_{0}$, this approximation not only fixes the main dependance of $\Delta t_{v}$ on the coordinates and velocities, but also suggest that the time delay in this case is largely insensitive to the angular position $\beta$ of the source. For the high order terms in Eq. (50), from Eq. (45) one can see that they will be maximal when $r_{0}$ is large while $r_{\mathrm{i}}$ and $r_{\mathrm{f}}$ is relatively small. Expansion (47) further suggests that $r_{0}$ is large only when $\beta$ is large. Therefore the dependance of $\Delta t_{v}$ on $\beta$ is stronger when $\beta$ is large and less so when it is small.

If $v$ and $v^{\prime}$ are very close so that their difference $\Delta v=v^{\prime}-v$ is much smaller than $v^{\prime}$ and $v$, e.g. for a null ray and an ultra-relativistic ray or two ultra-relativistic rays, then a further expansion of Eq. (49) around $v$ can be carried out to find

[TABLE]

where

[TABLE]

and $r_{0\pm\mathrm{w}}$ was still given by Eq. (47) and consequently

[TABLE]

To reveal more physical insights from these results, we plot in Fig. 2 the time delay $\Delta t_{v}$ using the definition (49) and total time formula Eq. (45). Note the $\Delta t_{v}$ depends on four variables nontrivially: $\beta,~{}v$ and $r_{\mathrm{f}}/M$ and $r_{\mathrm{i}}/M$. The typical range of $\beta$ is from $10^{-6}$ [as] to 10 [as] and $v/c$ is between 0 and 1. For $r_{\mathrm{i}}$ and $r_{\mathrm{f}}$, for simplicity we fix them to be equal and denote them collectively by $D$. It can range from $(1\sim 10~{}\mathrm{[kpc]})/(10^{-1}\sim 10~{}[M_{\odot}])$ in the typical microlensing case to $(1\sim 10^{2}~{}\mathrm{[Mpc]})/(10^{6}\sim 10^{10}~{}M_{\odot})$ if the lensing is due to supermassive black holes (except the Galatic one), to even $(10^{2}\sim 10^{5}~{}\mathrm{[Mpc]})/(10^{10}\sim 10^{12}~{}M_{\odot})$ for typical lensing by galaxies. Therefore we will use a range of $(10^{-7}\sim 10^{2}~{}\mathrm{[kpc]})/M_{\odot}$ for $D/M$. The value of $M$ is implicitly fixed at $M=4.12\times 10^{6}M_{\odot}$ in this plot. $\Delta t_{v}$ at other values of $M$ can be obtained by scalings.

In Fig. 2 (a), the time delay $\Delta t_{v}$ between signals with velocity $v=(1-10^{-15})c$ and light as a function of $\beta$ and $D$ is plotted. This velocity $v$ is chosen to represent an allowed value of the GW velocity GBM:2017lvd ; Monitor:2017mdv . It is seen that as discussed in Eq. (50), the dependance of the time delay on $D$ is linear, while its dependance on $\beta$ is not noticeable in this plot. For $D=8.02$ [kpc], a source located at the same distance from the galaxy center as Earth and then lensed by the galactic supermassive BH, the $\Delta t_{v}$ for signal with this velocity is about 0.0020 [s]. This is more than one order smaller than the 0.054 [s] uncertainty in the measurement of time delay between GW and GRB GBM:2017lvd ; Monitor:2017mdv and therefore calling for improvement in GW and GRB detection time accuracy if $\Delta t_{v}$ in lensing of galactic central BH for galactic merger event is ever used. In Fig. 2 (b), the time delay $\Delta t_{v}$ between light and signal with arbitrary velocity $v$ is plotted as a function of $v$ and $D$ for $\beta=1$ [as]. This more clearly verified the observation in Eq. (50) that $\Delta t_{v}$ is not only linear to the distance but also linear to $\left(\frac{1}{v}-\frac{1}{v^{\prime}}\right)=\left(\frac{1}{1-\Delta v}-1\right)\approx\Delta v$. It is seen that when the velocity difference decreases from $v=(1-10^{-14})c$ to $v=(1-10^{-16})c$, the time delay also decreases linearly by two orders for all $D$. We also numerically verified that, if $\beta$ is changed to larger value (e.g. 10 [as]) or smaller value (e.g. $10^{-6}$ [as]), the change in plot Fig. 2 (b) is indeed unnoticeable, in agreement with previous observation.

IV.2 Time delay $\Delta t_{\mathrm{p}}$

For particles traveling on two sides of the lens along the $b_{-}$ and $b_{+}$ paths respectively, the time delay can be expressed using Eq. (45) as

[TABLE]

where $v_{\pm}$ are the velocity of the signal on two sides of the lens and $r_{0}(v_{\pm},|b_{\pm}(v_{\pm})|)$ are given by Eq. (40). This equation is valid for comparing any two kinds of signal with arbitrary velocities.

In astrophysical observation, it is often the case that total travel time of signals with same velocities from different lensing images are compared. In this situation, $v_{+}=v_{-}=v$ and the time delay Eq. (54) becomes after using Eq. (45)

[TABLE]

where we have ignored the last term in Eq. (45), which is valid when $v$ is not extremely small: $(v^{2}/c^{2}>M/r_{0\pm\mathrm{w}})$.

When $\beta$ is small, the $|r_{\mathrm{0+w}}|-|r_{\mathrm{0-w}}|\propto\beta^{1}$ will also be small and therefore the perturbative expansion of Eq. (55) in powers of $\beta$ and then in powers of $M/r_{\mathrm{i/f}}$ can be done. To the $\beta^{3}$ order and leading order of $M/r_{\mathrm{i/f}}$ in each order of $\beta^{i}$, one finds

[TABLE]

In order for this expansion to converge, then we should demand that

[TABLE]

If $\beta$ is in this range, then clearly the first term will dominate and therefore it is expected that the time delay $\Delta t_{\mathrm{p}}$ will be proportional to $\beta^{1}$ and $r_{\mathrm{i/f}}^{1/2}$. Note however for some gravitational lensing with large $\beta$ and/or $r_{\mathrm{i/f}}$, this condition is violated (see Ref. Treu:2010uj for ranges of parameters $\beta$ and $r_{\mathrm{i/f}}$) and therefore for those cases the expansion (56) will not be accurate. In those cases, one can easily show that the logarithmic term of Eq. (55) will be much smaller than its first term, which after substituting Eq. (47) for $r_{0\pm\mathrm{w}}$ becomes

[TABLE]

Note that when the relativistic limit of the velocity $v$ is taken, the ${\cal O}(1-v)^{1}$ order term in Eq. (58) agrees with Eq. (20) of Ref. Glicenstein:2017lrm . If we further take the limit of large $r_{\mathrm{i}}$ and $r_{\mathrm{f}}$, this becomes

[TABLE]

Unlike the situation in Eq. (56), this time delay is proportional to $\beta^{2}$ and $r_{\mathrm{i}}^{1}$.

Regarding the dependance of $\Delta t_{\mathrm{p}}$ on the signal velocity $v$ that is close to $c$, in both cases of Eq. (56) and (59), a further expansion of their first term implies that the dominate part of $\Delta t_{\mathrm{p}}$ is always proportional to $1+\Delta v$. This is very close to 1 for relativistic particles, suggesting that in this case the variation of the time delay due to velocity change is much smaller than the time delay itself.

In Fig. 3, we plot the time delay $\Delta t_{\mathrm{p}}$ as functions of various parameters using Eq. (55). In Fig. 3 (a), $\Delta t_{\mathrm{p}}$ as a function of $\beta$ and $r_{\mathrm{i}}=r_{\mathrm{f}}=D$ is plotted for light signal $(v=c)$ using logarithmic scale. It is seen that for fixed $D$ and increasing $\beta$, $\Delta t_{\mathrm{p}}$ increases linearly with slope 1 in the log-scale plot when $\beta$ is small. Similarly, for fixed $\beta$ and increasing $D$, $\Delta t_{\mathrm{p}}$ also increases linearly but with slope 1/2 when $D$ is small. These features are in agreement with Eq. (56). When Eq. (57) is about to be violated as $\beta$ and $D$ increase, i.e., when

[TABLE]

a transition from Eq. (56) to (59) happens. This can be seen from the coincidence of the red line representing Eq. (60) and the bending region in Fig. 3 (a). Beyond this line, the slope of the plot in both the $\beta$ and $D$ directions are doubled, reflecting that Eq. (59) now takes place.

In Fig. 3 (b), $\Delta t_{\mathrm{p}}$ as a function of $\beta$ and $v$ is present. It is seen that for the range of $\beta$, the same transition from small $\beta$ expansion Eq. (56) to large $\beta$ expansion Eq. (59) happens. Moreover, for the entire parameter range, the dependance of $\Delta t_{\mathrm{p}}$ on $v$ is not noticeable in this log-scale plot, as argued below Eq. (59). Indeed, one can plot $\Delta t_{\mathrm{p}}$ as function of $v$ and $D$ too and also find the very weak dependance on $v$.

In Sec. V, we will be interested in the difference of two $\Delta t_{\mathrm{p}}$’s of different velocities $v$ and $v^{\prime}$, i.e., $\Delta t_{\mathrm{p}v}\equiv\Delta t_{\mathrm{p}}(v)-\Delta t_{\mathrm{p}}(v^{\prime})$, which can be calculated using Eq. (55). In the entire parameter ranges of $M$, $\beta$ and $r_{\mathrm{i/f}}$ considered in Fig. 3 (a), one can verify that when both $v^{\prime}$ and $v$ are close to $c$, the contribution from the $M/v^{2}$ term in $r_{0\pm\mathrm{w}}$ in Eq. (47) to this difference can be ignored. Further expand this difference to the first order of $(v^{\prime}-v)$, the result is found as

[TABLE]

Similar to the situation in Eqs. (56) and (59), the coefficient of $(v^{\prime}-v)$ in Eq. (61) also depends on $\beta$ and $r_{\mathrm{i/f}}$ as $\beta^{1}r_{\mathrm{i/f}}^{1/2}$ when they are small and as $\beta^{2}r_{\mathrm{i/f}}^{1}$ when they are large. Later in Sec. V we will use this to find the time delay difference between two neutrino mass eigenstates in the neutrino lensing case and between GRB and GW signals in the binary neutron star merger case.

V Time delay of neutrinos and GWs

V.1 SNN time delay

The neutrino mass ordering and the absolute value of neutrino masses are important problems for not only particle physics but also cosmology and astrophysics. The latest constraints on the neutrino mass square differences are Tanabashi:2018oca

[TABLE]

and the cosmological bound on the sum of the neutrino masses is Ade:2015xua

[TABLE]

Using the Eqs. (62) to (65), we can estimate the masses of neutrinos for both the normal and inverted orderings. Assuming the lightest neutrino is massless, Eq. (62) to (64) suggests that for normal order

[TABLE]

and for inverted order

[TABLE]

If we assume that the bound (65) is saturated, then the corresponding neutrino masses for normal order are

[TABLE]

and for inverted order

[TABLE]

We first consider the time delay of SNN signals from the same side of the lens. In this work, we focus on the SNNs because their properties are better understood comparing to neutrinos of other astrophysical origin Aartsen:2016ngq ; IceCube:2018cha ; IceCube:2018dnn . Because neutrinos have three mass eigenstates, for any given SNN spectrum that usually last a few seconds the three mass eigenstates will decouple from each other during propagation from supernova to observer which costs long time. Therefore it is expected that three separate signals corresponding to the $|\nu_{1}\rangle$, $|\nu_{2}\rangle$ and $|\nu_{3}\rangle$ eigenstates will be received from the same side of the lens. Here we will show that these three signals will have a time delay that might be used to resolve the neutrino mass ordering problem.

We assume that the SNN has a fixed energy of 10 MeV which is about the average of the spectrum Hirata:1987hu ; Bionta:1987qt . Using the masses in Eqs. (66) to (69), we then can calculate the velocity $v_{i},~{}v_{j}$ of each mass eigenstate $|\nu_{i}\rangle$ and $|\nu_{j}\rangle$ respectively and use Eq. (49) to find the time delay $\Delta t_{v,ij}$ between them, i.e., $\Delta t_{v,ij}\equiv t_{\mathrm{if,w}}(v_{i})-t_{\mathrm{if,w}}(v_{j})$, for both mass orderings. For simplicity, we assume that the lens is located at 2 [Mpc] away from both the observer and the source and the source angular position $\beta=1$ [as]. In Fig. 4, we show the time delays of both the normal and inverted orderings with the masses given by Eqs. (66)-(67). It is seen that for the normal ordering, there is a delay of about 0.1 [ms] for $|\nu_{2}\rangle$ signal comparing to $|\nu_{1}\rangle$ and about 2.6 [ms] for $|\nu_{3}\rangle$ comparing to $|\nu_{2}\rangle$. For the inverted ordering, $|\nu_{1}\rangle$ appears 2.6 [ms] after the $|\nu_{3}\rangle$ and $|\nu_{2}\rangle$ appears 0.1 [ms] after $|\nu_{1}\rangle$. In both the mass orderings, if $|\nu_{1}\rangle$ and $|\nu_{2}\rangle$ signal are to be resolved, then these signals should have a characteristic time that is narrower than 0.1 [ms]. Fortunately for supernovae that collapse into BHs, it is known that the SNN spectrum tail has a characteristic termination time that last usually about $2R/c\approx 0.1$ [ms], where $R\approx 10$ [km] is the size of neutrino emission region in supernova. If the location of the lens and source are 30 time larger, then it can be seen through Eq. (50) that the time intervals between the mass eigenstates will be 30 larger too. That is, $|\nu_{1}\rangle$ and $|\nu_{2}\rangle$ are separated by about 3 [ms]. This is about the minimal time duration of another feature that is widely believed to exist in SNN spectrum, the neutronization burst peak Hempel:2011mk ; Lentz:2011aa ; Bruenn:2012mj ; Couch:2013kma . Therefore in this case the mass eigenstates from the neutronization burst peaks might also be split in time in different ways in these two mass orderings. These all suggest that the two mass orderings will appear as two different sequence of events separated by different time intervals if the distance of the lens and supernova are large enough.

We also varied the mass from Eqs. (66)-(67) to (68)-(69) and repeated the calculation in Fig. 4. It is found that $\Delta t_{v,ij}$ in both mass ordering scenarios are independent of the masses in the given range. Indeed, for $\Delta t_{v}$ between ultra-relativistic particles with same energy but different masses, one can replace the velocity in the first term of Eq. (50) by $v_{i}=\sqrt{1-m_{i}^{2}/E^{2}}$ and then expand around the small $m_{i}/E$. One finds that $\Delta t_{v,ij}$ is equivalent to the time delay found for a single gravitational potential in Ref. Zatsepin:1968kt , which is given by

[TABLE]

where $D_{os}$ is the total distance from the SN to the observer. Clearly, the leading term of the time delay is only sensible to mass square difference but not the absolute value of the masses of the neutrinos.

Indeed, the time delay between different mass eigenstates of neutrinos originating form a gravitational potential has been considered in Ref. Jia:2017oar using Eq. (70) with more realistic SNN spectrum and neutrino-matter interaction cross-sections. The findings there was similar to what was observed in Fig. 4 that the mass eigenstates might be separated by different time intervals in different mass orderings. However, it was pointed out that in order to have enough statistics, the distance $D_{os}$ cannot be too large even for a gigaton water or liquid scintillator detector. Limited $D_{os}$ in Eq. (70) then implies that in order to have enough temporal resolution, only sharp features with very short time duration ($\leq 0.1$ ms) in the SNN spectrum can be used for the purpose of discriminating the mass orderings. Comparing to the case without GL however, the time delay in GL case fortunately do have an important advantage that the neutrino flux is significantly magnified, up to 100 times Rubin:2017ipu . This will strongly increase the detection event rates for sources from the same $D_{os}$ and therefore makes the method more practical.

For SNN signal from different paths, the time delay $\Delta t_{\mathrm{p}}(v_{i})$ between same mass eigenstate $|\nu_{i}\rangle$ with velocity $v_{i}$ is given by Eq. (55). Comparing this with Eq. (70) we can find that for $\beta>10^{-3}$ [as], $\Delta t_{\mathrm{p}}(v_{i})$ is much larger than $\Delta t_{v,ij}$. Therefore the two series of neutrino signals will not have any overlap for most of the ranges of $\beta$. If we consider the difference of two time delays, $\Delta t_{\mathrm{p}}(v_{i})$ of the ultra-relativistic $|\nu_{i}\rangle$ signal and $\Delta t_{\mathrm{p}}(c)$ of the optical signal, then Eq. (61) should be used. Using $v_{i}=\sqrt{1-m_{i}^{2}/E^{2}}$, it becomes

[TABLE]

Formally, measuring $\Delta t_{\mathrm{p}v}$ and $E$ allows the determination of the absolute mass of $|\nu_{i}\rangle$ for given $\beta,~{}r_{\mathrm{i/f}}$ and $M$. This was not possible when using solely the time difference $\Delta t_{v}$ because the difference in emission times of neutrino and optical signal in SN cannot be determined very precisely (e.g. SN1987A), while using the difference of two time delays can avoid this uncertainty. Practically however, one can find using reasonable $r_{\mathrm{i}}$ and $r_{\mathrm{f}}$, and typical $\beta$ and $E$ that this time difference is too small to be experimentally resolved if the neutrino masses are in the range specified by Eq. (66) to (69). For example, in the large $\beta$ and $r_{\mathrm{i/f}}$ limit, Eq. (71) becomes after restoring all units

[TABLE]

Therefore even for features in SNN spectrum that is as narrow as the neutrino observatory uncertainty ($\sim$1 [ns]), to resolve the peaks form different mass eigenstate would require an extremely large $r_{f}$. This in turn requires extremely large detectors to reach high enough statistics. Therefore until such detectors are built, this practically will not put any constraint on the neutrino absolute mass.

V.2 Time delay in a general mass profile and time delay of GW

In the analysis of the GW170817 and GRB 170817A signal Monitor:2017mdv , it was deduced from the $+1.74\pm 0.05$ [s] time delay of GRB that the speed of GW is constrained to the range

[TABLE]

The main uncertainty in Eq. (73) comes from the fact that this time delay is not necessarily the difference of the total travel times of GW and GRB, because their emission time could also be different. To avoid this problem, Refs. Fan:2016swi ; Wei:2017emo ; Yang:2018bdf proposed to use the difference between time delays of GW images and time delay of GRB images in GL to accurately determine the GW velocity. However, the time delay used in these works (Eq. (15) of Ref. Fan:2016swi , Eq. (1) of Ref. Wei:2017emo and Eq. (2) of Ref. Yang:2018bdf ) was derived for a singular isothermal sphere profile from the the time delay between lensed and unlensed rays of light Biesiada:2007rk

[TABLE]

but not the time delay of signal of arbitrary velocity. Therefore it requires a revision. Note in Eq. (74), redshift $z$ was set to zero because we are in a Schwarzschild spacetime, $\psi(\mathbf{\theta})$ is the effective lensing potential at angle $\mathbf{\theta}$ and $D_{ol}$ and $D_{ls}$ are the distance from observer to lens and lens to source respectively. In this work we will update this equation, and show that the first (geometric) term in the bracket of Eq. (74) receives a factor of $1/v$, and the second (potential) term gets a factor of $(3v^{2}-1)/(2v^{3})$ where $v$ is the speed of the lensed signal. In other words, the time delay formula becomes

[TABLE]

In order to illustrate this, we only need to consider the time delay caused by a point mass, i.e., a Schwarzschild spacetime. The time delay formula for an light ray in this case was known to be Keeton:2005jd

[TABLE]

Now for timelike particles, their motion in Schwarzschild metric satisfy the following normalization condition

[TABLE]

where $\lambda$ is the proper time and $l$ is the length parameter along the path $\mathrm{d}l^{2}=\mathrm{d}r^{2}+r^{2}(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2})$ and the approximation is valid because $r\gg M$. Using the definition of energy $E$ in Eq. (3) to replace 1 on the left hand side of Eq. (77), it can be rewritten as

[TABLE]

Introducing Newtonian potential $U(r)=-\frac{M}{r}$ and rearranging the equation, the travel time can be expressed as

[TABLE]

Since $U(r)$ is small, expanding this equation to the first order of $U(r)$ yields

[TABLE]

where $E$ was replaced by $v$ using Eq. (4).

To carry out the integral in Eq. (81), we do a change of variables from length parameter $l$ to $x$ in Fig. 1 by using the geometric relation

[TABLE]

so that the time from the source to lens plane becomes

[TABLE]

Because $\xi,~{}\xi^{\prime}\ll D_{ls},~{}D_{ol}$, this integral can be carried out to the leading orders of $\xi/D_{ls}$ or $\xi^{\prime}/D_{ls}$ to find

[TABLE]

Similarly, the time from the lens plane to the observer is found to be

[TABLE]

and the total travel time becomes

[TABLE]

The time delay between the travel time for the lensed ray and unlensed ray then is

[TABLE]

Since we are in the weak field limit, the angle $\beta$ and $\theta$ are small. From Fig. 1 we have

[TABLE]

which lead to a solution of $\xi$ and $\xi^{\prime}$ in terms of other observables

[TABLE]

Substituting Eqs. (89) into Eq. (87), one finally finds

[TABLE]

Note that if we expand $\Delta t_{g}$ for relativistic velocity, this result is in accordance with Eq. (7) of Ref. Glicenstein:2017lrm (although it seems a factor 2 difference occur at the $(1-v)^{1}$ order).

Comparing to formula (76), one can immediately read off the extra factors due to velocity in the geometric term and the potential term. Since the time delay for general potential is a convolution of point mass potential, then it is apparent that for a general potential $\psi(\mathbf{\theta})$, the time delay of an relativistic particle with velocity $v$ becomes

[TABLE]

which is exactly the desired Eq. (75). Note that if the relativistic limit is taken, then the first term agrees with Eq. (19) of Ref. Baker:2016reh at the ${\cal O}(1-v)^{1}$ order. The potential term there however has the same dependance on velocity as the first term in the bracket at order ${\cal O}(1-v)^{1}$, which is at odds with our result and Ref. Glicenstein:2017lrm .

Eq. (91) is applicable to general mass profiles to study the difference between two time delays of different signals. However, here we will only not further pursue along this direction but concentrate on the simple point mass result, Eq. (61). For GW and GRB signal, their time delay difference becomes

[TABLE]

Clearly, this difference is linear to $1-v_{\mathrm{GW}}$. Its dependance on $\beta$ and the distance is the same as in Eq. (61) and similar to Eq. (58) .

Since previously $v_{\mathrm{GW}}$ has been constrained to the range of Eq. (73), using the maximum deviation of $v_{\mathrm{GW}}$ from $c$ that is still allowed, i.e., $1-v_{\mathrm{GW}}/c=3\times 10^{-15}$, we can estimate the difference $\Delta t_{\mathrm{p}v,\mathrm{GW}}$. In Fig. 5 we plot $\Delta t_{\mathrm{p}v,\mathrm{GW}}$ as a function of $\beta$ and $D$ for velocity $\Delta v/c=3\times 10^{-15}$. It is seen that the parameters considered take their maximal values that are considered, i.e., $D_{ol}=D_{ls}=2\times 10^{4}$ [Mpc] and $\beta=10$ [as], the difference in GW and GRB time delay also reach maximum, $\Delta t_{\mathrm{p}v,\mathrm{GW}}=1.45\times 10^{-5}$ [s]. This distance corresponds to $z=4.3$ and is already 10 times the maximal distance of currently detected GW candidate gwevent1 and slightly larger than the largest distance of detected strong optical GLs opgllist . Similarly, the value of $\beta$ we used here is also larger than known optical GL values opgllist . On the other hand, it is also known that current uncertainties in the measurement of GW event time is 0.002 [s] and that of GRB event time is 0.05 [s] Monitor:2017mdv . These numbers are 2 and 3 orders larger than the above difference of $1.45\times 10^{-5}$ [s] at largest $D$ and $\beta$. Therefore, in order to further constraint GW velocity, either these uncertainties have to be improved by 2 and 3 orders respectively for GW and GRB measurements, or some very rare gravitationally lensed binary neutron merger event from a distance larger by a factor of 3 orders, i.e., $\sim 10^{7}$ [Mpc], has to be detected.

To evaluate the best constraining power to velocity difference $\Delta v\equiv 1-v_{\mathrm{GW}}/c$ by a given measured $\Delta t_{\mathrm{p}v,\mathrm{GW}}$, in Fig. 5 (b) we plot $\Delta v$ as a function of $\beta$ and $D$ assuming that a 1 [ms] time delay difference was measured by the detectors. The upper surface is for $M=4.12\times 10^{11}M_{\odot}$ and the lower surface is for $M=4.12\times 10^{-1}M_{\odot}$. It is seen that the larger the $D$ and $\beta$, the smaller the $1-v_{\mathrm{GW}}/c$ can be constrained. Moreover, at large $D$ and $\beta$, the two surface due to different $M$ converge. Indeed, taking the large $D$ and $\beta$ limit, one can obtain the velocity difference $1-v_{\mathrm{GW}}/c$ from Eq. (92) as

[TABLE]

This shows clearly that constraining power of lens and source at different distance and angular positions.

VI Conclusion and Discussions

The time delay of timelike particles in GL is important for astrophysical applications such as constraining GW velocity and neutrino mass/mass orderings. In this work, we studied the time delay of signals with arbitrary velocity in GL of the Schwarzschild spacetime. Exact formula of the total time $t_{\mathrm{if}}$ is obtained in Eq. (31) as elliptical function and the approximation result of $t_{if}$ is found in Eq. (45) under the weak field limit. Both the exact and approximate $t_{\mathrm{if}}$ are functions of the gravitational center mass, the source and observer distances, the particle velocity and minimal radius $r_{0}$, the last of which is linked to the angular position of the source using lens equation.

The time delay $\Delta t_{v}$ for signals of different velocities but on the same side of the lens and various limits of $\Delta t_{v}$ are obtained in Sec. (IV.1). The time delay $\Delta t_{\mathrm{p}}$ for signals on opposite sides of the lens and its approximations are obtained in Sec. (IV.2). The dependance of $\Delta t_{v}$ and $\Delta t_{\mathrm{p}}$ on various parameters including $\beta$, $v$ and $r_{\mathrm{i/f}}$ are discussed carefully.

These time delays are applied to the time delay of SNNs and GW. It is shown that the time delay $\Delta t_{v}$ between relativistic neutrino mass eigenstates $|\nu_{i}\rangle$ and $|\nu_{j}\rangle$ is proportional to $D_{os}\Delta m^{2}_{ij}/E$. Therefore the three mass eigenstates will yield a sequence of neutrino signals that is different for normal and inverted orderings. This implies a possibility to discriminate these orderings, although a very large detector is required to have enough statistics.

For GW application, we first updated the formula of time delay in GL of a general mass profile for signal with arbitrary velocity. Then in Schwarzschild spacetime, it was shown that for distance as large as $2\times 10^{4}$ [Mpc] and source angle of 1 [as], the difference in time delays of GW signal with velocity $v_{\mathrm{GW}}=(1-3\times 10^{-15})c$ and GRB signal can only reach $1.45\times 10^{-5}$ [s]. To utilize this difference to further constraint GW velocity, the measurement uncertainty of GW and GRB has to be improved or very high redshift GW/GRB event has to be observed.

It is instructive to comment on the future extensions of the current work. The first is to apply the time delay formula for general mass profile to more specific lens mass distributions and study the corresponding implications on GW properties. The perturbative methods used in the work can also be extended to other spacetimes, preferably those with spherical/axial symmetries, such as Kerr spacetime. It would be interesting to know how the spin angular momentum of the spacetime affects the time delay of massive particles from different images in GL. We are currently working along these directions.

Acknowledgements.

The authors appreciate the discussion with Dr. Yungui Gong and Xilong Fan. This work is supported by the National Nature Science Foundation of China No. 11504276.

Appendix A Definitions of the elliptic functions

And the special functions appear in the equations are elliptic integral. Then we introduce the definitions of the elliptic integrals. The elliptic integral of the first kind is

[TABLE]

The elliptic integral of the second kind is

[TABLE]

The incomplete elliptic integral is

[TABLE]