Time delay in the strong field limit for null and timelike signals and its simple interpretation

Gravitational lensing can happen not only for null signal but also timelike signals such as neutrinos and massive gravitational waves in some theories beyond GR. In this work we study the time delay between different relativistic images formed by signals with arbitrary asymptotic velocity $v$ in general static and spherically symmetric spacetimes. A perturbative method is used to calculate the total travel time in the strong field limit, which is found to be in quasi-series of the small parameter $a=1-b_c/b$ where $b$ is the impact parameter and $b_c$ is its critical value. The coefficients of the series are completely fixed by the behaviour of the metric functions near the particle sphere $r_c$ and only the first term of the series contains a weak logarithmic divergence. The time delay $\Delta t_{n,m}$ to the leading non-trivial order was shown to equal the particle sphere circumference divided by the local signal velocity and multiplied by the winding number and the redshift factor. By assuming the Sgr A* supermassive black hole is a Hayward one, we were able to validate the quasi-series form of the total time, and reveal the effects of the spacetime parameter $l$, the signal velocity $v$ and the source/detector coordinate difference $\Delta\phi_{sd}$ on the time delay. It is found that as $l$ increase from 0 to its critical value $l_c$, both $r_c$ and $\Delta t_{n,m}$ decrease. The variation of $\Delta t_{n+1,n}$ for $l$ from 0 to $l_c$ can be as large as $7.2\times 10^1$ [s], whose measurement then can be used to constrain the value of $l$. While for ultra-relativistic neutrino or gravitational wave, the variation of $\Delta t_{n,m}$ is too small to be resolved. The dependence of $\Delta t_{n,-n}$ on $\Delta \phi_{sd}$ shows that to temporally resolve the two sequences of images from opposite sides of the lens, $|\Delta \phi_{sd}-\pi|$ has to be larger than certain value.


I. INTRODUCTION
Deflection of light signal near massive celestial bodies was the most convincing phenomenon that helped establishing General Relativity as the correct description of gravity [1]. Nowdays deflection of signals by gravity is almost exclusively used in gravitational lensing (GL), which has become a powerful tool in astrophysics and cosmology. GL has been used to constrain not only properties of the lens, such as the mass distribution of galaxy [2][3][4], the structure of dark matter halos [5][6][7], accretion of materials [8,9], but also those of the source, such as supernova explosion mechanism [10,11]. Observables in GL can also correlate with properties of signal particles forming the GL images [12].
Traditionally, GL observations were always done using light signals of various wavelength.
However, with the discovery of neutrinos form SN1987A [13,14] and more recently from blazer TXS 0506 + 056 [15,16], and the discovery of gravitational waves (GW) due to binary black hole (BH)/neutron star mergers [17][18][19], especially the GBR+GW dual observation [20], it is clear that neutrinos and GW can also act as messengers to have GL effects.
In (strong) GL observations, the apparent angle of the images and the time delay between them are two most important observables that were widely used to reveal information about the lens, source and messenger. For both neutrino and GW observatories, currently their angular resolution are both too low to distinguish different images. However, the time measurement for these signals are usually very precise, reaching O(1) [ns] for neutrino events [21,22] and O(1) [ms] for GW observations [20]. Therefore theoretical study and corresponding observations on time delay in GLs of these signals might bear fruit earlier than resolving different images formed by them.
In this work, we study time delay of timelike and null signals with general asymptotic velocity in the strong field limit (SFL) in static and spherically symmetric (SSS) spacetimes.
In this limit, the signal's trajectory approaches the critical particle/photon sphere and the signal might loop around the central lens many time before reaching the detector, forming series of relativistic images from each side of the lens. The time delay of light signal in SSS spacetimes in this limit has been studied by Bozza in Ref. [23] as a distance estimator and then followed by many works in particular spacetimes or gravitational theories [24][25][26][27][28][29][29][30][31][32][33][34]. In this work, we not only extend it to arbitrary signal velocity, but develop a trackable way to calculate the total travel time and time delay to any desired order, which was never done before. Moreover, we also show that the time delay is given by a simple formula, Eq. (36), allowing a very simple and intuitive understanding, i.e., Eq. (37).
The work is organized as follows. In Sec. II, we develop the perturbative method used for the computation of the total travel time in the SFL. We will show that the total travel time takes a quasi-series form, which is then used in Sec. III to find the time delay between different images. It is show that the time delay equals to the circumference of the particle sphere divided by the local velocity of the signal and then multiplied by the redshift factor and the winding number. In Sec. IV, we then apply these results to the Hayward BH spacetime and study the dependance of the time delay on the spacetime charge parameter l and signal velocity v. The result reveals that the time delay can be used to constrain l quite well but not v since for both supernova neutrinos and GWs, their speeds have been well constrained to be extremely close to light speed.

II. THE PERTURBATION METHOD
The perturbative method we used here to calculate the total travel time is adopted from Ref. [35] which was to calculate the deflection angle in the SFL, and here this method is used to a different integral. Therefore in this section, we will first recap the essential steps to help to understand the method, and then apply the method to the suitable integral defining the total travel time of the signal.
We start from the most general SSS metric described by where (t, r, θ, ϕ) are coordinates and A, B, C are metric functions depending on r only.
The corresponding geodesic equations can always be transformed to the equatorial plane (θ = π/2), and then becomeṙ where κ = 0, 1 for massless and massive signals and L, E are the angular momentum and energy of the signal (per unit mass). L, E can relate to the impact parameter b and the asymptotic velocity v of the signal by The travel time from a source at r = r s to a detector at r = r d (see Fig. 1), after using Eqs. (2) and (3), is then where r 0 is the closest approach of the trajectory to the lens. According to Eq. (2), r 0 Together with Eq. (5), this provides a relation connecting the closest approach r 0 and the On the other hand, when r 0 or b is small enough, the signal will spiral into a compact sphere, the particle (or photon) sphere, and not escape back to infinity. The radius r c of this sphere is defined as the critical point of the denominator of Eq.
If r 0 approaches r c from above, we call the corresponding b the critical impact parameter and denote it by b c . Using Eq. (8), it is related to r c by The total travel time (6) in the SFL is difficult to compute analytically even for the simpler SSS spacetimes and therefore approximation methods are desired. In Ref. [35] we have developed a perturbative way to expand a similar integrand for the computation of signal's deflection angle. Here we adopt that method to the computation of the total travel time in the SFL. First, we define a function p(x) inspired by Eq. (8) as and define its inverse function as q(x). From Eq. (8), it is clear that so that taking its inverse function we have We then can do the following change of variable in Eq. (6) from r to ξ connected by Note that although the analytical form of p(x) is clear once the metric function is given, the inversion process to find q(x) is not always possible analytically. Fortunately, what is required in our later computation is the series expansion of q(x) and it can always be worked out using the Lagrange inverse theorem from the series form of p(x). Under this change of variable, and noticing Eqs. (11), (12), the integration limits and various factors in the integrand of Eq. (6) are changed according to where q = q ξ bc and q is its derivative. Here the η s,d are nothing but the sine value of the apparent angles of the signal at the source and detector respectively [cite a few of our works].
In the SFL and large r s,d limits, we have a → 0 + and η s,d → 1 − respectively. Grouping these terms together, we obtain the transformed total travel time as Note that this form of t depends on the impact parameter b through the parameter a, on the source/detector radius through η s,d . Its dependence on all other parameters of the spacetime is through the metric functions and the critical impact parameter b c , which also appears in A. Perturbative expansion of total travel time The beauty of the above change of variable (14) is that it transforms the infinite integration range to a finite range and allow the resultant integrand to be expanded in the small ξ limit, so that an perturbative integration can be carried out.
To see how the expansion is carried out, we first split the integrand into two factors, 1 √ ξ−a and y(ξ) with The factor 1 √ ξ−a will be directly integrated later while the factor y(ξ) should be further treated. Since in the SFL, the main part of the total time is contributed from the integration near small ξ, we can further split y(ξ) into two factors and expand them in the small ξ limit, i.e., where in Eq. (18b) the index starting from −1 because of the q term, and the f n are the expansion coefficients that can be worked out once the metric functions are known. We point out that it is in expansion (18b) that the series form of q(ξ/b c ) is needed and can be obtained using the Lagrange inverse theorem from the series expansion of p(x). We also emphasis that the coefficients f n will not depend on the initial/boundary conditions of the trajectory, such as the impact parameter b and r s,d , but only on the metric functions and spacetime parameters therein.
Collecting these expansions according to the powers of ξ, we see that y(ξ) becomes where a sum over a finite terms of a m appears because different powers of the denominator factor (2 − a) in Eq. (18a) mix into the coefficient of the same ξ n/2 power. The coefficients y n,m can be obtained from the coefficients f n and other factors in Eqs. (18) The integrability of this formula relies on the integration of the form Fortunately, this kind of integrals can always be worked out for integers n, and the results are given in Eq. (A1) in Appendix A. Using these results, the total travel time t becomes Here the first and second terms in the curl bracket are due to the integration of odd and even powers of ξ n/2 respectively. It is also important to notice that all the dependence of t on a and η s,d in this formula has been shown explicitly, and the coefficients y n,m only contain the signal kinetic parameter v and spacetime parameters through b c .
It is not difficult to notice that all the functions involved in t are quite elementary. Because of this, in the SFL (i.e, b → b + c , a → 0 + ), it can be further expanded into a quasi-series of small a, which yield a from where in the coefficient of each order of a k , there is only one term that contains ln a. The coefficients C k and D k can be worked out from Eq. (21) and it is seen that only the D k 's (but not the C k 's) depend on η s,d . This total time (22) resembles the same form as the deflection angle in the SFL in Ref. [35], although their coefficients C k and D k will be different. When from which we can read off the coefficients C 0 and D 0 in Eq. (22) as Later on, we will use the total time (23) to calculate the time delay between different relativistic images in the SFL. Different images correspond to different impact parameters and consequently, different a, but their r s,d and other spacetime parameters are exactly the same. This also implies that when subtracting two total times along the trajectories with slightly different b, the D 0 term will not contribute to the time delay at this order.

B. Computing coefficients y n,m
From the relation (13) and the change of variable (14) we knew that the ξ → 0 + limit is also the r 0 → r + c , b → b + c limit. Therefore the expansions (18) at small ξ or equivalently the coefficients y n,m should also be determined from the series expansion of the metric functions at r = r c . Assuming these expansions are where a i , b i and c i are the coefficients, then the very definition of r c in Eq. (9) becomes a constraint between the first few coefficients Using these metric expansions and going through the process from Eqs. (17) to (19), the y n,m 's can be computed. In particular, the first two y n,0 are found to be where Higher order y n,0 (n > 0) and y n,m (m > 0) can also be readily computed but are too long to present here. However from these higher orders, we are able to assert that for general n, y n,0 is determined by the metric expansion coefficients up to a n+3 , b n+1 and c n+3 . In time measurement in GL, what is measured is not the total travel time but the time delay between different images of the source. For GL in the SFL, a typical path is illustrated in Fig. 1. There will exist two basic kinds of time delay: the time delay between signals from the same side of the lens but with different winding numbers around the center, and the time delay between signals from different sides of the lens but with the same winding number around the center. A general time delay, which we denote as ∆t n,m (n, m ∈ Z) to represent the time delay between images winding around the center counter-clockwisely n and m times, should be a combination of the former two basic time delays. Using Eq. (24), this time difference can be written as

III. TIME DELAY IN THE SFL
where b n and b m are the impact parameters of the two images. Note the D 0 a 0 term in Eq.
(24) dose not affect the time delays since it is independent of b and therefore the same for all trajectories. Furthermore, here n, m could be negative integers if the winding is indeed clockwise (although b n and b m are always positive). Note that the case that winding dose not actually happen corresponds to the weak field limit of the time delay and was considered in our previous work [36].
Clearly, in order to compute the time delay, we shall find the corresponding impact parameters for the images in the SFL first. For this purpose, we will directly use the result from Eq. (38) of Ref. [35], where ∆φ sd = φ s − φ d is the difference between the angular coordinates φ s of the source and φ d of the detector (see Fig. 1). The coefficient while the coefficient D a0 is given in Eq. (28) where C 0 in Eq. (25), C a0 in Eq. (34) and b c in Eq. (29) are substituted and simplified. It is seen that the time delay depends only on the following parameters: the metric expansion coefficients a 0 , c 0 , the asymptotic signal velocity v and a angular factor determined by the number of loops n, m and angular coordinate difference ∆φ sd , but not on the finite distance r i and r f of the source and detector.
We emphasis that this time delay is a very general result: it applies to GL with general asymptotic velocity v, general source/detector angular coordinate difference ∆φ sd , general SSS spacetime with a particle sphere and arbitrary n and m. Setting v = 1 in Eq. (36) reduces it to previously known result in Ref. [23] (Eq. (40) and (41)), which concentrated on null signals.
Although the source and detector in GL in the SFL are usually far from the BH center, one would expect however when the winding number of two signals n and m are both large (in this case, n, m ≥ 1 are enough), the time delay between them, observed by a far away observer, should be equivalent to the circumference 2πr c of the particle sphere divided by local signal velocity v l and then multiplied by the difference of loops ∆l (not necessarily an integer) between the two paths, and finally the gravitational redshift factor γ from the particle sphere to the detector. That is, it is natural to expect that Here we show that indeed, the above is exactly the time delay result (36), which is found from more rigorous and lengthy calculations.
For arbitrary SSS spacetime, it is always possible to choose the metric function C(r) = r 2 and therefore the its expansion coefficient at r c is c 0 = r 2 c . That is, the particle sphere circumference 2πr c = 2π √ c 0 , i.e., the numerator of the first factor in Eq. (36). For the local velocity v l , in the SFL the signal circulates around the particle sphere, and therefore we only needs to consider the velocity due to angular motion. In an SSS spacetime, this is where γ sr = 1/ 1 − v 2 l is the special relativity gamma factor. Then using Eqs. (4), (5) (setting b to b c ), (29), and c 0 = r 2 c sequencically, it is just an elementary algebra to solve the local velocity as which is exactly the denominator of the first factor in Eq. (36). Thirdly, the difference in the number of loops ∆l, after properly taking into account the opposite direction case, is simply the last factor of Eq. (36). Lastly, in an asymptotically flat SSS spacetime described by metric (1), the gravitational redshift factor from r c to the detector which is located far away is simply 1/ A(r c ) = 1/ √ a 0 , the second factor of Eq. (36). Grouping these factors together, therefore it is verified that for general SSS spacetime and timelike or null signal, the time delay in the SFL, Eq. (36), has a very simple and intuitive understanding, Eq. (37).

IV. THE HAYWARD BH SPACETIME CASE
In this section, we apply our result to some particular spacetimes to check its validity, and examine the effect of v, ∆φ sd , number of loops n and m, and more importantly the spacetime parameters. The spacetime we study is the Hayward BH spacetime whose metric functions are [37] A where M is the spacetime mass and l is the charge parameter. |l| < 4M/(3 √ 3) ≡ l c in order for the spacetime to be a BH one. Using Eqs. (5) and (9), the equation determining the particle sphere radius r c of this spacetime becomes This is an six order polynomial of r c whose solution dose not have a closed algebraic form.
However, after formally or numerically solving it, then substituting r c into the metrics (40) and further into Eq. (10), the critical impact parameter b c is found as To solve the y n,m that are needed in the total travel time (23) and the time delay (36), then we should expand the metric functions at r = r c according to Eq. (27). The first few of these expansion coefficients are Substituting them into Eqs. (30) and (31), we can obtain the coefficients y n,0 for the Hayward spacetime. The first two of them, denoted as y −1,0,H and y 0,0,H , are where  To verify whether the time delay Eq. (46) is accurate, we can first truncate the sum in its second term to order m and define and then compare t m (r s , r d , b) with the total time t num (r s , r d , b) obtained directly from numerical integration of Eq. (6). As long as the numerical integration is done to high enough precision, t num can be thought as the true value of the total travel time. In Fig. 2, we plot t 1 , t 5 , t 10 , t 20 , t 30 as well as t num for b around b c , which is about 5.19M when we set l = 0.1M and r s = r d = 20M, v = 1. It is seen that as the truncation order increases, the total time approaches its true value for all b considered, and t 30 is non-distinguishable from t num in the plot. Moreover, in the inset we see that for any fixed truncation order (we took t 10 as an example), the smaller the b − b c , the better it approximates the true value of the travel time.
Although the truncated t m (r s , r d , b) in Eq. (48) approximate the true value pretty well in GL, a numerical study shows that to achieve the same accuracy in the total travel time as in Fig. 2, the truncation order of t m (r s , r d , b) would be formidably high, even in the SFL.
Fortunately, this will not affect the accuracy of the time delay (36) because in the SFL, the time delay between different trajectories mainly happens when the signal is very close to the particle sphere, for which part the total time (48) is a very good approximation. In other words, for any two trajectories, the travel times corresponding to the parts from large r d or r s to some small radius -below which the time delay happens -cancel out. This indeed leaves the time delay formula (36) very accurate even we truncated at a relatively low order, as can be seen from Fig. 3.
The effects of the spacetime charge l and signal velocity v, as well as ∆φ sd when sign(n) = sign(m) are shown in Fig. 3(b)-(d) by assuming that the Sgr A* supermassive BH is a Hayward BH. To help understanding these effects, we first plot the particle sphere radius r c as a function of l/M and v in Fig. 3(a). It is seen that for any particular signal velocity, as the charge l increases from 0 to l c = 4/(3 √ 3M ), r c decreases monotonically from its Schwarzschild value to a minimal value. In particular, for light signal, this is from 3M to about 2.65M . For fixed l and decreasing v on the other hand, r c increases from its light signal value to a larger but still finite radius. Comparing to the Reissner-Nordstrom spacetime, we see that the effects of l and v on r c are qualitatively similar to those of the electrostatic charge and particle velocity in that spacetime [38].
Then for the time delay ∆t n,m,H in Eq. (47), it is clear that as l increases to l c , r c decreases and changes of terms in both the denominator and numerator of the first factor cancel to a large extent. Therefore this factor has a very minimal variation. Similar trend happens for the second term and therefore the time delay depends very weakly on l, as can be seen from Fig. 3(b). In the entire range of l from 0 to l c , the time delay ∆t n,m,H for light signals from   [20]). Therefore using time delay caused by Sgr A* in the SFL to constrain speed of such ultra-relativistic signals seems not likely.
Finally, for the case of signals from opposite sides of the lens, i.e., sign(n) = sign(m), the dependence of the time delay on ∆φ sd is shown in Fig. 3(d). When n = −m and ∆φ sd = π, i.e., the source is perfectly aligned along the observer-lens axis and the trajectories from two sides have an mirror symmetry, then clearly we should have ∆t −m,m,H (∆φ sd = π) = 0, as shown in Fig. 3(d). Unlike GL in the weak field limit, angle ∆φ sd dose not need to be very close to π if GL in the SFL can really be observed in the future. As ∆φ sd deviates from π, the time delay ∆t −1,1,H becomes linear to (∆φ sd − π). The arrival time of the each series of images from one side of the lens will form an arithmetic sequence which is equivalent to ∆t n,1,H (n = 2, 3, · · · ) or ∆t n,−1,H (n = −2, −3, · · · ). The two sequences from two sides will have a relative shift ∆t −n,n,H that is linear to (∆φ sd − π) too, as shown in Fig. 33(d).
Because of this relation, it is seen that for a given characteristic time scale of the source event or observatory time resolution (whichever is larger), there exist a minimal ∆φ sd − π that allows the two sequences to be temporally separated. Taking 2 [ms] for the time resolution of GW signal as an example, then it is seen that for the two sequence to be separated by this interval, |∆φ sd − π| has to be larger than 2 [ ]. On the other hand, if GL of an event is to be both observed in the weak field limit and temporally resolved in the SFL, then since GL in the weak field limit are usually observed for β 10 [ ] [65] or roughly ∆φ sd 20 [ ], Fig. 3(d) implies that both the characteristic time scale of the event and the observatory time resolution have to be smaller than 20 [ms]. Therefore to temporally resolve the two sequences from two sides certainly imposes stringent requirement on the time measurement uncertainty of observatories, e.g., the GRB measurement uncertainty which is current about 50 [ms] [40].

V. CONCLUSION AND DISCUSSION
In this work we proposed a perturbative method to compute the total travel time t and time delay ∆t in the SFL in SSS spacetimes for signal with arbitrary asymptotic velocity.
The total travel time takes a simple quasi-series form where a = 1 − bc b and coefficients C k and D k can be expressed as rational functions of the metric expansion coefficients at the particle sphere radius. In the SFL, the leading contribution to ∆t is given by the ln a term. Using impact parameter corresponding to each relativistic image, we were able to show that ∆t is given by Eq. (36). This result allows an intuitive and yet quantitatively precise understanding: to the leading order of a, the time delay is given by the circumference of the particle sphere divided by the local velocity of the signal and then multiplied by the winding number difference and the redshift factor from the particle sphere to the far away detector.
We applied the results of t and ∆t to the Hayward BH spacetime. The correctness of the total travel time is verified by truncating the series to different orders. The time delay in this case is found in Eq. (46). To understand it, we first studied the dependence of the particle sphere radius r c on the spacetime charge parameter l and signal velocity v. It is found that as l increases or v decreases, r c decreases or increases respectively. Assuming the Sgr A* central BH is a Hayward BH, we were able to compute ∆t between images on the same side and opposite sides. It is found that as l increases from 0 to its critical value, the time delay per loop can vary by 7.2 × 10 1 [s], which is well within the time resolution of typical light, neutrino or GW observatories. Therefore measuring ∆t caused by the Sgr A* BH will constrain its charge l very well. On the other hand, for the supernova neutrinos or GW whose velocities deviate from that of the light by 10 −15 or less, measuring difference between ∆t's of different signals would not further constrain their velocities, because of the time measurement accuracy/characteristic time of the corresponding events are larger by a few orders.
Regarding directions to extend the current work, the first and most straightforward one is to extend the perturbative method to the case of the equatorial plane in the stationary and axialsymmetric spacetime. Based on the weak field limit experience [64], we expect that the spin parameter would play a non-trivial role in affecting the time delay between relativistic images with different winding directions. A more interesting extension is to apply the method to the time delay of asymptotically non-flat spacetimes. From the quasi-series (22) for the total time, we saw that the coefficients C k and D k are completely determined by the behavior of the metric functions around the particle sphere. Although GL in the weak field limit in these spacetimes is often problematic to study due to the difficulty to take the infinite radius limit, the metric functions behave normally at small radius and therefore the SFL can still be taken. We are pursuing along these directions.
The high order coefficient of y n,0 in Eq. (30) can be worked out with the help of a symbolic computation tool. Here we will only give one more order, i.e., where (A3)