A kinematical approach to gravitational lensing using new formulae for refractive index and acceleration

This paper uses the Schwarzschild metric to derive an effective refractive index and acceleration vector that account for relativistic deflection of light rays, in an otherwise classical kinematic framework. The new refractive index and the known path equation are integrated to give accurate results for travel time and deflection angle, respectively. A new formula for coordinate acceleration is derived which describes the path of a massless test particle in the vicinity of a spherically symmetric mass density distribution. A standard ray-shooting technique is used to compare the deflection angle and time delay predicted by this new formula with the previously calculated values, and with standard first order approximations. Finally, the ray shooting method is used in theoretical examples of strong and weak lensing, reproducing known observer-plane caustic patterns for multiple masses.


INTRODUCTION
The deflection of light by massive bodies has been of interest to mathematicians and physicists from time to time since Newton suggested the possibility (Newton (1704)). This deflection was calculated in the late eighteenth and early nineteenth centuries, treating light as a classical particle. The deflection was again calculated by Einstein in the early twentieth century, using his new general theory of relativity, to be twice the previous classical result. The measurement of the deflection of light passing close by the sun was widely publicized as a dramatic confirmation of general relativity, in the now famous 1919 expeditions.
In the last three decades, gravitational lensing has become an important tool for astrophysicists, especially in searching for dark matter (first suggested in Paczynski (1986)) and exoplanets. By 1991, astrophysicists were suggesting that exoplanets could be found using microlensing, and since 2004 at least ten planets have been found in this way (Sumi et al (2010)). In microlensing, light from a background star passes close to the lensing system, and is deflected around the lens. Because of this, more light rays reach the observer, producing magnification of the background source. This magnification changes over time, as the source, lens and observer move into and out of alignment. The details of the magnification over time are plotted in a 'light curve', which is simply intensity versus time. A planet orbiting a lensing star can make changes in the magnification pattern at the observer's plane. Such changes show up as variations to the shape of the simple light curve.
Interpretation of these light curves is difficult, as this is an inverse problem, in which observers seek to determine the details of the lensing system which gave rise to the observed data. In particular, researchers are often trying to find planets in the lensing system, and to determine characteristics of such planets, primarily orbital radius and mass. Details on reproducing a model of such a multi-body lens were outlined at least as early as 1996 (Wambsganss (1997)). For a brief history see Schneider et al (1992), pages 1-2, Mao et al (2007), or the recent review by Ellis (2010). In this paper, microlensing will be considered, although the formulae derived here are just as valid for macrolensing situations.

Lensing Models
It is customary (Wambsganss (1997)) to use a 'thin lens' model, in which the effect of the lensing system is confined to the plane containing the lensing objects (the 'lens plane'), this plane being normal to the line from observer to lens (or alternatively, source to lens). A deflected light ray thus consists approximately of two straight lines, with an abrupt angle  Time to closest alignment Magnification Factor change in the lens plane. The magnitude of this change is given by a simple formula: ∆θ = 2rs/r0 where r0 is the point of closest approach of the light ray to the star or planet and rs is the Schwarzschild radius, rs = 2M G/c 2 . Here, M is the mass of the star, G is Newton's gravitational constant and c is the speed of light in a vacuum.
Current methods (Zabel & Peterson (2003) and Wambsganss (1997)) shoot rays from the observer to the lens plane, deflect by the angle as calculated above, and then draw the ray from there to the source plane. Equivalently, light rays may be followed from the source to the observer, mimicking more closely the actual physics. The density of rays at the source plane (alternatively, the observer's plane) is thus mapped. By tracing various linear paths across this 'magnification map', to simulate the relative movement of the source star, light curves are generated by plotting density of rays as a function of time. Fig. (1) shows schematically the method for mapping the amplification for an observer travelling at constant speed in the observer's plane. Simple light curves such as the one in the right of Fig. (1) provide information about the mass of the lensing object, provided that the distance to lens and distance to source can be estimated. If the lens is a binary system, there are deviations from this simple light curve, which may yield information about the orbital distance and relative masses in the system.
Thus the purpose of a microlensing model is to produce a magnification map due to light being deflected by the lensing system. Note that every such magnification map allows for arbitrarily many light curves, depending on the location of the observer's (or source's) path across the plane. Because of this, it is, in general, a difficult matter to find a model to fit an observed light curve for a lens involving more than one mass (that is, a binary or planetary system). The problem becomes much more difficult when the lens involves more than two objects (Mao et al (2007)).
The approach presented in this paper is to generate a magnification map by shooting rays through a lens of smoothly varying refractive index. The path will be curved rather than two straight line segments. Even so, the majority of the deflection occurs very close to the lens plane.

Kinematical Approach
It is the purpose of this paper to consider gravitational lensing using an alternative approach. This approach can be described as an equivalent classical kinematical formulation, which nevertheless replaces Newton's formula for gravitational acceleration with a formula derived from general relativity. Using this kinematical approach, each curved light ray is traced from the source to the plane of the observer, rather than using the simple angle change formula commonly used in lensing models.
We will start by considering the path predicted by general relativity for a massless test particle (a photon) near a massive body. The changes in angle along this path will be integrated to give the deflection of the light ray predicted by general relativity. As a check, we will compare this result with known formula for approximating the deflection to first order in the relevant small constant. Next we will derive an 'effective refractive index' due to the massive body. Integration will be used to predict the Shapiro delay, which will also be checked against the known first order prediction of general relativity. A new acceleration formula will then be derived and tested against these known values for the deflection and delay using a forward integration method. Finally, we will describe and exemplify a procedure for using the new acceleration formula to produce magnification maps for multi-body lenses.

LIGHT PATHS IN A SCHWARZSCHILD SYSTEM
The Schwarzschild metric describes spacetime outside a non-rotating, electrically neutral, spherically symmetric mass density distribution M and is taken as a valid approximation for local spacetime structure in the vicinity of any massive body (stellar systems, including black holes, but also planets) having negligible charge and angular momentum. In spherical coordinates the metric components are given in terms of the invariant interval dl 2 as (Misner et al. (1973)) where rs = 2M G/c 2 is the Schwarzschild radius, G is Newton's constant, and c is the constant speed of light in vacuum. Here r , θ, φ are coordinates which correspond to standard spherical coordinates in the reference frame of an observer at rest far from the system. Light travels on null geodesics with dl = 0, so equation (1) becomes cdt = r r − rs dr 2 + r(r − rs)(sin 2 θdφ 2 + dθ 2 ). ( Applying Fermat's principle, that light follows a path that extremizes travel time T , we can consider the functional: To illustrate the use of equation (3) we firstly derive a path equation which is equivalent to the usual trajectory equations for null geodesics, and numerically integrate it to give the deflection of a light ray according to general relativity.
Without loss of generality, for a single light ray, the coordinates can be oriented so that the ray is in the plane θ = π/2. Then dθ = 0 in equation (3), which consequently may be re-arranged to give adopting the notation φ ′ = dφ/dr . Let F be the integrand in equation (4), that is: Then the Euler-Lagrange equation is It can be seen in equation (5), that there is no explicit dependence of F upon φ so clearly ∂F/∂φ = 0. Therefore this term vanishes in equation (6) which then has the immediate first integral in which K is a constant. Now 1/φ ′ = dr/dφ, so rearrangement yields the first order separable ODE This constant K can be determined. At the point of closest approach to the mass (call this point r = r0), the radius is at a minimum, that is, dr/dφ = 0, so it follows that: Thus the path is defined by: It can be easily seen that substituting u = 1/r followed by differentiation gives the well known second order equation as (Capozziello et al (1997)) From equation (7), an integral can easily be written to evaluate the total change in φ of a light ray passing a massive object, which will be twice the change from the perihelion out to infinity: , This elliptic integral cannot be evaluated with a finite number of simple algebraic terms. Numerical integration could be used to evaluate the deflection. However, the integrand is infinite at r = r0, making the accuracy of any numerical evaluation questionable. By using a substitution r = 1/ cos(2ψ) , the singularity can be removed. After simplification we obtain: This integrand is well behaved over the interval, so it can be numerically integrated to any desired accuracy, for example, by gaussian quadrature. An undeflected ray will have △φ = π, so the deflection for a ray will be δ = △φ − π. The deflection is usually very small, so in order to avoid a 'loss of significance error', this subtraction should be performed in the integrand. For readability, let α = 2r0/rs, and let β = 2 − 4 cos 2 ψ − sec 2 ψ. In this notation, the deflection δ can therefore be expressed as The two terms in the integrand are combined, so that the subtraction can be performed explicitly. This yields: For the purposes of this paper, we take the solar radius to be 696000 kilometres, and the solar Schwarzschild radius to be 2.95 kilometres. Numerical integration of equation (9) for a ray passing near the surface of the sun (rs/r0 = 2.95/696000))gives a deflection of 1.74851634161261 arcseconds. As the path equation contains all of the information about the general relativistic path of the photon, this is the deflection predicted by general relativity to the level of precision shown. This deflection angle will be used later to confirm the accuracy of the kinematic approach. As a check, we can consider Einstein's estimate for the deflection (Schneider et al (1992), page 3) 2rs r0 = 1.74850913341648 arcseconds.
This estimate is found to correspond to our calculated value to first order in rs/r0, as expected.

NEW REFRACTIVE INDEX AND TRAVEL TIME DELAY
The approach presented here makes use of an expression for the refractive index, n. The Schwarzschild metric from section 2 will be used to derive an 'effective' refractive index. As suggested above, the functional (3) can be arranged in the form T = dt = (dt/ds)ds, where s is an arbitrary parametrization of the ray path. By choosing ds to be the element of arclength along the path, the speed is then v = ds/dt, and thus the refractive index is n = c/v = cdt/ds, and we obtain finally T = 1/c nds, with the refractive index then having the form where r ′ = dr/ds, θ ′ = dθ/ds and φ ′ = dφ/ds. (Note that there is no suggestion here of any physical effect on the speed of light -indeed, any local measurement of coordinate velocity is guaranteed to result in the usual speed c).
Taking again the two dimensional case, with θ = π/2, the refractive index has the form As an example, the delay can be calculated for light to travel from an object at earth radius, skim past the sun, and back out to earth radius. This is the calculation needed for the radar echo delay test of general relativity (Shapiro (1964)). (In fact, for that particular test, it would be necessary to consider another planet such as Venus or Mars, and the calculation would need to be performed for each leg of the journey. For the purposes of this paper, it is simpler to imagine a reflecting satellite at the same orbit as the earth). The problem can be pictured as in Fig. 2, not to scale. The deflection is small, so that the path appears as a straight line. In fact the curved path that the photon takes is derived above in equation (7). The path is symmetric about the point of closest approach (r0), so that exactly half the delay can be obtained by integrating from r0 to re, the radius of the earth's orbit. If the sun had no effect on the light path, the distance from perihelion to earth would be the straight line distance r 2 e − r 2 0 . The delay can be calculated using the new refractive index, equation (11), and multiplying by ds so that Again, (dφ/dr) 2 is given by the path equation (7). The time taken for the trip will be twice the time to go from the point of closest approach to the sun, r0 to the earth's orbit, re, that is dr.
The integrand is infinite at r = r0, so before integrating, we make the substitution cos(ψ) = r0/r. After rearrangement we get: This integrand is perfectly well behaved over the interval, so it can be integrated to arbitrary precision by using, for example, Gaussian Quadrature. As for the deflection, the delay is very small, so it is important to subtract the straight-line time before integrating. The straight-line time can be written as: The delay is obtained by subtracting equation (13) from equation (12) to give: where γ = 1 − (rs/r0) cos ψ 1 − rs cos 2 ψ (r0 − rs)(1 + cos ψ) This expression (14) is now in a form that minimizes errors caused by subtraction of large terms. Performing this calculation in Matlab gives a delay of 129.0896086 µs. As for the deflection, this calculation can be performed to arbitrary precision, and accurately describes the delay predicted by general relativity.
We now turn to a kinematic approach, and consider differential equations that relate position, velocity and acceleration. In a Newtonian system, the acceleration would be given by Newton's law of gravitational attraction, g = −c 2 rs/2r 2 er. This simple formula is not appropriate for the present application, and a new form will now be derived from the present relativistic approach.

A NEW ACCELERATION FORMULA
By combining the metric equation (1) with the path equation (7) for a photon in a Schwarzschild orbit, the velocity and acceleration of the photon due to the nearby mass are now derived. Here, the meaning is that of a coordinate acceleration. As a freely falling particle, the photon does not experience any locally measurable force. Beginning with the path equation (7), and making the substitutions the path equation becomes dr rdφ 2 = Ar 2 − µ.
Next, considering equation (1), setting dl 2 = 0, θ = π/2 and dividing both sides by dt 2 , we have whereṙ = vr is the radial velocity component and rφ = v φ is the tangential velocity component. Using the path equation (17), we can solve for vr and v φ in turn, to get: Thus the velocity vector of the photon along its path is To determine the acceleration vector, take the derivative with respect to time: a =vrer + vrėr +v φ eφ + v φėφ In polar coordinates, the derivatives of the unit vectors areėr =φeφ andėφ = −φer, and so the acceleration components ar and a φ are ar =vr −φv φ , a φ =v φ +φvr.
Differentiating vr in equation (19) yields: Equations (19) and (20) are now used to eliminate the square root terms, so thatvr simplifies tȯ The radial acceleration component is therefore The acceleration must be related directly to the Schwarzschild radius rs of the mass. There are two terms in equation (21) that do not have an rs coefficient. These two terms cancel if the positive sign is chosen. Thus, the correct form for the radial acceleration is: A similar treatment for tangential acceleration component yields Thus the acceleration vector for a photon near a Schwarzschild mass is: It is interesting to note that the radius for light to remain in a circular orbit about the mass can immediately be derived from this acceleration vector. In a circular orbit, there is no radial velocity, and so a = rs r 2 − 3v 2 2 er.
In addition, an object moving in a circular orbit in a classical kinematical framework has centripetal acceleration vector When equations (23) and (24) are equated, we obtain Thus, a photon at the '3/2' radius given in equation (25) is trapped in a circular orbit about the mass, in accordance with the known result predicted by general relativity (Carroll (2004), p 212).

Kinematic Ray Shooting
With known acceleration components, it is now possible to set up a standard system of differential equations for ray tracing in polar coordinates. The kinematical system is This new system (26) makes use of the new acceleration formula in equation (22). Such a system can be solved using forward integration. Initial conditions for the photon have initial position at perihelion of 696000 kilometres about a mass with Schwarzschild radius of 2.95 kilometres, zero radial velocity, and tangential velocity v φ = c/n = c √ µ to the right. The speed of light is taken to be c = 300000 kilometres/second. Using Matlab's ODE45 routine (an explicit Runge-Kutta 4-5 method) produces the path shown in Fig. (3) for the section shown; beyond this, the path is almost a straight line and so is not shown. The slope of the line between the last integration point before the photon reaches earth orbit, and the first point after is 1.74851634 ± 10 −8 arcseconds. The time delay is calculated in the same integration, and is found to be 129.089609 ± 10 −6 µseconds. The uncertainty is due to limitations on the precision of Matlab's ODE45 routine. Both of these values correspond well with the predictions from general relativity (as calculated above), more closely than the first order approximations, and use of a higher precision computation will allow a more accurate result, should such be required. Accuracy beyond first order is not commonly required, but these results give confidence that the acceleration vector presented here does accurately embody the effect on the photon due to a single Schwarzschild-type mass.
Using the values found for position and radial and tangential velocities when the photon reaches earth orbit, we can send the photon along the path from earth orbit past the sun and back out to earth orbit. The central section is shown in Fig. (4). Note that the scales differ by a factor of 10 6 .
The values for the deflection and the Shapiro delay calculated by this forward integration ray shooting are compared with the predictions from general relativity and the usual first order approximations(as calculated earlier) in Table 1, and demonstrate that the kinematic ray shooting method presented in this paper is an accurate representation of the effect of the gravitational field of a single Schwarzschild body on the motion of a photon, giving us some confidence in using this method in approximating more complicated systems. x (kilometres from perihelion) y (kilometres from perihelion)

Magnification Maps
When considering a multi-body system, such as a planetary system, it must be stressed that there is no known metric. That is, there is no known exact solution to Einstein's equations for such a system. Some sort of approximation is therefore required. Use of the 'weak-field metric' is one such approximation, as is the addition of the deflection angles due to each body (the method generally used in microlensing models). Here, we choose to approximate by adding the acceleration components due to each body in the system.
Having tested the radial and tangential acceleration components described above, it is a simple matter to set up a three dimensional ray tracing system for a planetary system. At each integration point along each ray, the acceleration components due to each massive body are calculated. This is done by a translation to put the massive body at the origin, followed by three rotations to place the photon's position vector and its velocity vector in the same plane as the massive body, with θ = π/2. The radial and tangential velocities are then used to calculate the radial and tangential acceleration components. The three rotations are then reversed, and the resultant cartesian acceleration components are added to the acceleration components due to any other masses in the system to determine the overall acceleration of the photon.
As a very simple (and artificial) example of this process, we first consider a two dimensional system, consisting of two very massive bodies in close proximity, and plot the path of the photon through this binary system. The smaller and larger stars represent bodies of 20 million and 50 million solar masses respectively (similar to the system described in Boroson and Lauer (2009), although in the present example we imagine that the system has decayed to the point where the black holes are only a billion kilometres apart). This example is purely to demonstrate the versatility of the present approach. We are making the gross simplification that the black holes are stationary throughout the period when the ray is passing through the system, and so the acceleration of the two masses towards each other is therefore being ignored. For the purpose of demonstrating the procedure used here, we are ignoring such limitations. The present model is clearly a coarse approximation in this extreme case. The light ray comes in from the far left, and is deflected by the summed acceleration components due to each mass. The path is shown in Fig. (5).
We now consider a more commonly modeled lensing system, the planetary lens. In describing a lensing system, it is common to use a parameter called the Einstein radius. This is the angular radius at which observers perfectly aligned with the point source and point lens would see a ring of light about the lens. Specifically, the system is designed as follows: a point source is at (−8000, 0, 0) and the observer is at (+8000, 0, 0). The lens star is placed at the origin, having Schwarzschild radius: rs = 99 * 10 −8 . An orbiting planet is placed at (0, 0.1208, 0) (in microlensing terms, it is at 1.35 times the Einstein radius), with rs = 1 * 10 −8 . For simplicity in this model, we ignore the motion of the planet. Rays are sent through the system, near the Einstein radius, and in the vicinity of the planet. Due to symmetry in the cases here, it is only necessary to calculate half the rays and plot the result both above and below the axis of symmetry. During the numerical integration, each ray is broken into several small sections, with the size of each section becoming smaller as the photon nears the lens star. This is important to ensure that the integration routine does not take too large a step and miss the strong deflection altogether. For the simulations presented in this paper, each ray is broken into 36 segments. The result of this process is a light density map, or magnification map (that is points in (y, z) where the rays cross the plane x = 8000). These simulations were run on MatLab version 6.1, under Windows XP, on an Athlon 64 3500+ processor with 1GB of ram. Running times are given as an indication, but no measures have been made to optimise the code for efficiency. Fig. (6) shows the magnification map produced when a rectangular array of rays (222 by 205) is sent through the lens system. The bending of the light towards the planet is clearly visible. This, combined with the bending caused by the lens star, produces the characteristic diamond shape for a system with a single planet outside the Einstein radius (Wambsganss (1997)). The running time was 26 hours. In Fig. (7), many more rays have been used (approximately 106,000 rays). In order to view the resulting density, it is necessary to colour or shade regions according to how many rays pass through each small area. The code used to do this was 'smoothhist2d', Perkins (2006). The running time was approximately 50 hours. This result clearly shows the caustic diamond structure.
This diamond structure is the expected shape for the magnification map, and suggests that the method used here can be considered as an alternative method for modeling a thin lens.

CONCLUSION AND DISCUSSION
We have considered the path of a photon near a Schwarzschild-type body. Using the Schwarzschild metric, a new refractive index has been derived. Integrating the angle along the path gives the total deflection angle, and integrating the new index along the path gives the travel time delay. These values for the delay time and deflection can be calculated, using the formulae here, to arbitrary precision. This is because these formulae are derived directly from the Schwarzschild metric. As a check, it was shown that the standard first order approximations used for the deflection and delay agree with these results to first order.
A new formula for the acceleration of a photon was derived by combining the Schwarzschild metric with the path equation, and differentiating. This new acceleration formula was tested with a ray shooting approach, using the new refractive index to provide the initial conditions for the velocity components. The deflection and delay values were found to be in excellent agreement with the precise values calculated earlier.
By making the approximation that the acceleration on the photon is the sum of the individual acceleration components due to each massive body, a simple microlensing model was developed to demonstrate a use of the new acceleration formula for a binary system. Sample light fields on an observer plane have been computed using this new approach, and reproduce the expected figures. No attempt has been made here to speed computations, since that was not the purpose of the present work, but future developments may address such issues of computational efficiency.
In summary, this work provides a 'classical' way of accurately describing the gravitational effect on a photon due to a single mass, and provides an alternative method for approximating the course of a photon through a complicated mass distribution. The authors believe that the approach presented here provides an insight into the effect of gravitating bodies on light rays that can be grasped without requiring a deep understanding of general relativity, and yet is still quantitatively accurate for a single mass, and can be used for approximating more complicated systems. As for the standard thin lens 'deflection angle' method, this approach to gravitational lensing may be used by applied mathematicians, computer modellers and others without requiring specialist knowledge of general relativity. Because this approach retains the 'delay' information as well as the deflection, it might conceivably be useful in analysing systems where the time delay plays a role, such as a pulsar source being lensed, should we observe such an event. It is also hoped that the formulae presented here will prove useful in producing models of more complicated mass distributions, such as galaxy clusters. Such models could be produced using the same method used here, simply by adding more bodies to the model.