Numerical estimation of the curvature of a light wavefront in a weak gravitational field
NNumerical estimation of the curvature of a light wavefront in aweak gravitational field
A. San Miguel, F. Vicente and J.-F. Pascual-S´anchez
Dept. de Matem´atica Aplicada,Facultad de Ciencias.Universidad de Valladolid,47005 Valladolid, Spain
Abstract
The geometry of a light wavefront evolving in the 3–space associated with a post-Newtonianrelativistic spacetime from a flat wavefront is studied numerically by means of the ray tracingmethod. For a discretization of the bidimensional wavefront the surface fitting technique is usedto determine the curvature of this surface at each vertex of the mesh. The relationship betweenthe curvature of a wavefront and the change of the arrival time at different points on the Earth isalso numerically discussed. a r X i v : . [ g r- q c ] O c t . INTRODUCTION The description of the propagation of light in a gravitational field is even today a centralproblem in the general theory of relativity. The deflection of light rays and time delays ofelectromagnetic signals due to the presence of a gravitational field are phenomena detectablewith current experimental techniques which allow design new tests for general relativity. Inthis line, Samuel [1] recently proposed a method for the direct measurement of the curva-ture of a light wavefront initially flat, curved when light crosses regions where the gravita-tional field is non vanishing. He found a relationship between the differences of arrival timerecorded at four points on the Earth, measured by employing techniques of very long baseinterferometry and the volume of a parallelepided determined by four points in the curvedwavefront surface. This surface is described by means of a polynomial approximation of theeikonal in a Schwarzschild gravitational field. For more complex gravitational models, suchas those considered by Klioner and Peip [2], de Felice et al. [3] or Kopeikin and Sch¨afer [4]in studies of light propagation in the solar system, the use of numerical methods for thedetermination of the geometry of the wavefront surface would also be required. An ana-lytical approach to the relativistic modeling of light propagation has also been developedrecently by Le Poncin-Lafitte et al. [5] and Teyssandier and Le Poncin-Lafitte [6], wherethey present methods based on Synge’s world function and the perturbative series of powersof the Newtonian gravitational constant, to determine the post-Minkowskian expansions ofthe time transfer functions.Nowadays there are numerous techniques in computational differential geometry whichallow to analyze geometric properties of surfaces embedded in the ordinary Euclidean space.Techniques of this type are widely applied in different areas such as Computational Ge-ometry, Computer Vision or Seismology. In one of these methods, developed in works byGarimella and Swartz [7] and Cazals and Pouget [8], the estimation of differential quantitiesis established using a fitting of the local representation of the surface by means of a heightfunction given by a Taylor polynomial. A survey of methods for the extraction of quadricsurfaces from triangular meshes is found in Petitjean [9].In this work we consider a discretization of the wavefront surface, replacing this surfaceby a polyhedral whose faces are equilateral triangles. At initial time, the surface is assumedto be flat and far enough from a gravitational source (say, the Sun) and moving towards2his source. We study the deformation of the instantaneous polyhedral representing thewavefront when crossing a region in the relativistic 3–space near the Sun due to the bendingof light rays by the gravitational field. In this study, we apply the ray tracing method withinitial values on the vertices of the triangular mesh to obtain the corresponding discretesurface at each instant of time. Then we apply the techniques given in [7] and [8] to describethe wavefront as a surface embedded in the Riemannian 3–space of the post-Newtonianformalism of general relativity. For each vertex in the instantaneous mesh we obtain aquadric which represents locally the surface by applying the least-squares method to theimmediate neighboring vertex around the considered point which is represented in normalcoordinates adapted to the light rays.The structure of the paper is as follows: In Section 2 we briefly introduce the basic modelfor the wavefront propagation in the post-Newtonian formalism. In Section 3, we establisha discretized model of the wavefront surface by means of a regular triangulation and wedescribe the method employed in this work for the study of the curvature of this surface.In Section 4, a numerical estimation of the curvature of the surface is derived using theray tracing method. For the numerical integration of the light ray equation, we use the
Taylor algorithm implemented by Jorba and Zou [10] which is based on the Taylor seriesmethod for the integration of ordinary differential equations and which allow the use of highorder numerical integrators and arbitrary arithmetic accuracy, as is required to describe theinfluence of weak gravitational perturbations on the bending of light rays. Finally, in Section5, we apply the method discussed above to study the effect of the wavefront curvature onthe variation of the arrival time of the light at points on the Earth surface, following themodel proposed in Samuel’s test [1]. The paper concludes by giving another approach to theestimation of the curvature of the wavefront, derived from an approximation of the Waldcurvature [12] associated with a quadruple of points in the wavefront.
II. LIGHT PROPAGATION IN A GRAVITATIONAL FIELD
Let us consider a spacetime ( M , g ) corresponding to a weak gravitational field and choosea coordinate system { ( z , ct ) } such that the coordinate representation of the metric tensor is g αβ = η αβ + h αβ , with η αβ = diag (1 , , , − . (1)3here the coordinate components of the metric deviation h αβ are given by: h ab = 2 c − κ (cid:107) z (cid:107) − δ ab , h a = − c − κ (cid:107) z (cid:107) − ˙ Z a , h = 2 c − κ (cid:107) z (cid:107) − (2)(Greek indices run from 1 to 4 and Latin indices from 1 to 3) where κ := GM represents thegravitational constant of a monopolar distribution of matter (say the Sun) located at Z a ( t )and c represents the light speed.In the post-Newtonian framework one may consider a simultaneity space Σ t at eachcoordinate time t . From the fundamental equation of the geometrical optics for the phase ψ ( z, t ) of an electromagnetic wave [11]: g αβ ∂ψ∂z α ∂ψ∂z β = 0 , (3)and Cauchy data ψ ( z,
0) = const, given on a spacelike surface D := { ( z, | φ ( z ) = 0 } ⊂ Σ one obtains the characteristic hypersurface (light cone) Ω := { ( z, t ) | ψ ( z, t ) = const } ⊂ M of the light propagation. The intersection of the characteristic hypersurface and thecorresponding simultaneity space is the spacelike wavefront at a time t which will be denotedby S t := Ω ∩ Σ t .An alternative formulation of the problem of light propagation in the spacetime ( M , g )is based on the determination of the bicharacteristics generated by the isotropic vectors k := grad ψ . From both the equation for the null geodesics, z ( t ) = (cid:0) z ( t ) , t (cid:1) expressed interms of the coordinate time and the isotropy condition g ( ˙ z, ˙ z ) = 0, one obtains in thepost-Newtonian approach to general relativity that, neglecting terms of order O ( c − ), thenull geodesics of ( M , g ) must satisfy the equations¨ z a = ϕ a ( z , ˙ z , t ) (4)0 = g αβ ˙ z α ˙ z β . (5)In (4) the components of the acceleration ϕ a ( z , ˙ z , t ) are given by (see [13]) ϕ a ( z , ˙ z , t ) = c h ,a − [ h ,t δ ak + h ak,t + c ( h a,k − h k,a )] ˙ z k − ( h ,k δ al + h ak,l − h kl,a ) ˙ z k ˙ z l − ( c − h k,j − c − h jk,t ) ˙ z j ˙ z k ˙ z a , (6)where the first and third terms in the right hand side of (6) are of order O (1), while theremaining terms are O ( c − ). 4or initial values ( z , ˙ z ), with z ∈ S and ( ˙ z /c, , T ] allows todetermine the spacelike wavefront S T . This surface is embedded in the Riemannian 3–space(Σ T , ˜ γ ) where the components of the metric tensor are given by˜ γ ab := g ab − g a g b g , (7)and at the time T the tangent vectors ˙ z ( T ) to the integral curves z ( t ) are ˜ γ –orthogonal to S T . III. NUMERICAL DESCRIPTION OF A SPACELIKE BIDIMENSIONAL WAVE-FRONT
Hereafter, we consider the simplest gravitational model generated by a static point mass.Let E denote the quotient space of M by the global timelike vector field ∂ t associated withthe global coordinate system used in the post-Newtonian formalism. We will consider aregion of the wavefront in E described by a coordinate chart { z } . Further, we assume thatthe Riemannian manifold ( E , ˜ γ ) is almost flat and that the metric corresponding to ˜ γ isquasi-Cartesian in the chosen coordinates. A. Discretization of the initial wavefront
Given an asymptotically Cartesian coordinate system, we consider a set S formed bypoints with coordinates ( z , z , − ζ ) where ζ > S maybe considered as a flat surface. The direction determined by the point O ∗ : (0 , , − ζ ) in S and the center of the Sun O : (0 , ,
0) is perpendicular to S . We will study the geometry ofa region C ⊂ S determined by points P whose Euclidean distances d ( P, O ∗ ) to the point O ∗ satisfy R (cid:12) (cid:54) d ( P, O ∗ ) (cid:54) R (cid:12) , where R (cid:12) is the radius of the Sun.For the discretization of the problem we consider, in the first place, the set of points( a , a , a ) ∈ A ∗ , where A = { , N , N + 1 , . . . , N } and N < N are two natural numbers(see Figure 1(a)). In A ∗ the point (0 , ,
0) is excluded. Next, we construct in the complexplane a regular triangular mesh whose vertices are located between two hexagons as shownin figure 1(b), and the edges have length (cid:96) . The inner hexagon has sides of length (cid:96)N and5 IG. 1: Triangulation and enumeration of the initial wavefront. (a) Triplets ( a , a , a ) ∈ A ∗ usedto label the vertex of a hexagonal mesh ¯ V . (b) Regular triangulation V ∗ of a hexagonal annularregion in the wavefront. Points (3 , , , (0 , ,
4) and (3 , ,
0) in ¯ V correspond with points 36 , V ∗ respectively. the outer hexagon sides of length (cid:96)N . In this triangulation each vertex is represented by acomplex number of the set (see [14])¯ V := { z = a + a ω + a ω | a , a , a ∈ A , ω := exp(2 π i / } , (8)where i := √−
1. The vertices ( a , a , a ) with some of their components equal to N (resp. N ) are located on the inner (resp. outer) boundary of the mesh ¯ V . We establish anenumeration of the vertices as shown in Figure 1(b), in such a way that the inner vertices z j have subscripts j = 1 , . . . , J . Finally, we apply the change of scale:¯ V → V ∗ , z (cid:55)→ zrN , (9)so that the inner boundary is a hexagon of radius r . In consequence the length of each edgein this triangulation is equal to r/N .The complex plane and the plane S may be identified by means of the mapping ι : z (cid:55)→ ( (cid:60) ( z ) , (cid:61) ( z ) , − ζ ). Thus one obtains a discretization of the initial region of S which we areconsidering here, triangulated with vertices given by ι ( z j ) whose corresponding mesh on S will also be denoted by V ∗ . 6 . Normal coordinates around a point Now we assume that at each vertex in V ∗ a photon with velocity ˙ z := (0 , , c ) is located.The null geodesics equation (4) may be written as a first order differential system ˙ u = F ( u , t )in phase space u := ( z , ˙ z ) which determines a flow in E : z ( t ) = Φ t ( z , ˙ z ) , (10)in terms of the initial values z := z (0) ∈ V ∗ , ˙ z := ˙ z (0). Then, for each time t thereis a surface S t image of S under the flow (10). At a point z ∈ S t the normal vector τ ( z ) coincides with the tangent vector to the curve z ( t ) at that point. Furthermore, beforereaching the focal points of the beam of light, the triangulation V ∗ induces a triangulation V on the wavefront S t whose vertices we enumerate using the same labels used for thecorresponding vertices in S .In order to simplify the description of the geometry of the surface S t on a neighborhoodof a point P ∈ S t we use a ˜ γ –orthonormal reference frame { e i } i =1 centered on that point,where one of its vectors, say e , is parallel to the vector n ( P ) := τ ( P ) / (cid:112) ˜ γ ( τ , τ ) tangentto the ray passing through that point. Let { y j } j =1 be a normal coordinate system with poleat the point P and associated normal reference frame { e i } i =1 .By using the classic formulae of Riemannian geometry (see [15] §
18, and [16]), the coor-dinate transformation y i (cid:55)→ z i , from normal to post-Newtonian coordinates, is determinedby y i = (Λ − ) ia (cid:0) z a − z a + (cid:0) ˜Γ abc (cid:1) ( z b − z b )( z c − z c ) (cid:1) , (11)where Λ is a non-singular constant matrix and (cid:0) ˜Γ abc (cid:1) := ˜Γ abc ( P ) are the Christoffel symbolsat the point P . By neglecting terms of order higher than ˜Γ abc , the inverse transformationmay be approximated by z a = a a + Λ ai (cid:0) y i − (cid:0) ˜Γ ijk (cid:1) y j y k (cid:1) , (12)whose corresponding Jacobian determinant is ∂z a ∂y i = Λ ak (cid:0) δ ki − (cid:0) ˜Γ kil (cid:1) y j (cid:1) . (13)In normal coordinates a metric tensor γ on the space E is determined from ˜ γ by γ ij = ∂z a ∂y i ∂z b ∂y j ˜ γ ab . (14)7nd, as it is well known, at point P the tensor γ ij is reduced to δ ij and the associatedChristoffel symbols at this point are Γ ijk = 0. C. Local approximation of the wavefront
From the triangulation V ∗ of the initial wavefront the ray tracing method furnishes adiscrete surface determined by the mesh V . To compute differential magnitudes of thewavefront surface corresponding to this mesh one needs to define a discrete neighborhood ofeach vertex in V . For this, we consider firstly at each inner vertex z ∗ j ∈ V ∗ , j = 1 , . . . , J aneighborhood (named hereafter in the terminology of computational geometry, e.g. [18]) formed by the six vertices z ∗ j k closest to z ∗ j :[ z ∗ j ; z ∗ j k ] j =1 ,...,J,k =0 ,..., , where z ∗ j k := z j + exp( kπ i / , (15)Then, for each 1–ring in (15) one may determine on the mesh V a corresponding 1–ringformed by the image of the points z ∗ j , z ∗ j k ∈ S under the flow (10)[ z j ; z j k ] := [ Φ t ( z ∗ j ); Φ t ( z ∗ j k )] . (16)In a neighborhood of a point z j the wavefront can be approximated by a height functionon the orthogonal plane to e determined by means of a least-squares fitting of the data(16) expressed in normal coordinates { y j } i =1 . As a model for this surface we chose a quadricpassing through the coordinate origin whose gradient at this point is parallel to e , y = f ( y , y ) := a ( y ) + a y y + a ( y ) , (17)where a , a and a are the indeterminate coefficients to be obtained by the least squaresmethod.The quadric (17) provides a surface ˜ S j that approximates the surface S t in a neigh-borhood of the point z j and is defined in parametric form y i = y i ( x A ), with A = 1 , y = x , y = x , y = f ( x , x ) . (18)In coordinates ( x , x ), the metric γ on E induces a metric on ˜ S j whose associated metrictensor g has the form g AB := γ ij ∂y i ∂x A ∂y j ∂x B . (19)8oreover, on the tangent plane to S t at z j one defines a tensor B associated with the normal n at that point as: B : T z j S × T z j S → ( T z j S ) ⊥ , B ( ∂ A , ∂ B ) = ∂ y ∂x A ∂x B n (20)Therefore, if { v , v } represents an orthonormal basis for the vector space T z j S consistingof eigenvectors of B with associated eigenvalues λ and λ , and II denotes the secondfundamental form, then the difference of sectional curvatures K and ¯ K associated withthe plane generated by { v , v } , in S and E respectively, is given by a generalized Gaussformula (see [19], p.131) K rel := K ( v , v ) − ¯ K ( v , v ) = λ λ (21)which is named (see [15]) relative sectional curvature K rel , whereas the mean curvature isdetermined by half the trace of II : H = ( λ + λ ) (22) IV. NUMERICAL ESTIMATION OF THE CURVATURE OF A WAVEFRONTA. Integration of the equations of light rays
Here we deal with the problem of the numerical integration of the initial value problemcorresponding to (4)-(5). The equation (4) for the light propagation in the gravitational fieldgenerated by a static material point, may be rewritten in terms of new variables defined as: u := ( u , . . . , u ) , with u i := z i and u i +3 := ˙ z i , ( i = 1 , , u = u , ˙ u = u , ˙ u = u , ˙ u = − κc c u − u u − u u u − u u u + u u + u u ( u + u + u ) / , ˙ u = − κc c z − u u u − u u − u u u + u u + u u ( u + u + u ) / , ˙ u = − κc c z − u u u − u u u − u u + u u + u u ( u + u + u ) / . (24)9he solution u ( t ) must satisfy the constraint (5) which, in the notation (23), may be ex-pressed as F ( u ) := u + u + u − c + 2 κc (cid:112) u + u + u ( u + u + u + c ) = 0 . (25)To obtain a numerical solution of the initial value problem given by (24) and the initialdata u (0) = (cid:0) z (0) , ˙ z (0) (cid:1) we use the Taylor integrator developed by Jorba and Zhou [10],based on the classic Taylor series method for ordinary differential equations. In this methoda Taylor expansion of the vector field ˙ u defined in (24) is made using techniques of automaticdifferentiation to obtain the corresponding Taylor coefficients. The Taylor integrator allowsthe control of both the order and the step size employed in the method. Furthermore, the
Taylor integrator is implemented so that one may use extended precision arithmetic for thehighly accurate computation required in this problem.From now on we use normalized units taking the radius and the mass of the Sun as unitsof length and mass, respectively. Then, the initial values corresponding to a photon initiallylocated at 100 astronomical units from the Sun are given by u := u (0) = (1 . , . , − . , . , . , . . (26)For the integration over a period of time of 50400 seconds, using a precision of 120 binarydigits and a tolerance Tol = 1 . E −
20, the
Taylor algorithm gives the solution in termsof Taylor polynomials of degree twenty-four. As a test of the accuracy of the method, inFigure 2 it is shown the behavior of the function F ( u ) appearing in the isotropy constraint(5). In this figure we see that at the arrival point, the isotropy constraint on the tangentvector to the light ray is not satisfied with the same degree of accuracy required at theinitial position, reaching this deviation its maximum value at the instant when the photonis nearest to the Sun.In order to obtain from the numerical solution ˜ u n given by Taylor another value u n satisfying the condition (25) and such that the function (cid:107) u n − ˜ u n (cid:107) reaches a minimum, wewill apply at each step in the Taylor integrator the method of standard projection [20] toproject ˜ u n on the manifold F ( u ) = 0. This leads to a constrained extremum problem witha Lagrangian function L ( u , λ ) := (cid:107) u − ˜ u (cid:107) − F ( u n ) T λ , where λ is a Lagrange multiplier.10 IG. 2: Semilogarithmic representation of the difference F (˜ u n ) − F ( u ) along the numerical solutionof the light ray equation obtained by means of the Taylor algorithm. The maximum value of ∆ F is reached at the time when the photon is nearest to the Sun. The necessary condition of extremum and the constraint condition leads to u n = ˜ u n + ∇ F ( u n ) T λ (27)0 = F ( u n ) (28)and replacing (27) in (28) the following nonlinear equation for λ is obtained F (cid:18) ˜ u n + ∇ u F (cid:0) u n ( ˜ u n , λ ) (cid:1) T λ (cid:19) = 0 . (29)This equation may be solved by applying the simplified Newton method. The projectionstage in the algorithm spent a 3% of the 0 .
06 seconds-CPU employed by an Intel Core 2 Quadprocessor to determine the trajectory of a photon in the time interval we are considering.
B. Deformation of a wavefront by a static gravitational field
Now we apply the method described in Section 3 to that region of a wavefront propagatingalong a tubular neighborhood around the Oz -axis and whose radius is 2 R (cid:12) . Initially, thewave surface S is flat, perpendicular to the Oz -axis and with position and velocity givenin (26). The wavefront S T is then determined after a trip of 101 astronomical units.
1. Curvature at a point
Firstly, we obtain an estimation of the curvature at a point in S T by using a sequenceof 1–rings with decreasing radii until a stable value of the curvature wavefront at that11 IG. 3: 1-ring with vertices [ V , V , . . . , V ] and corresponding least-squares quadric fitting of thesevertices represented in normal coordinates centered at V . point is reached. Let us consider the 1–ring of radius r and centered at the point V [0]0 :=(1 . , . , − . V [ k +1]0 := (cid:0) r cos( kπ/ , sin( kπ/ , − . (cid:1) , ( k = 0 , . . . . (30)Now, we apply the Taylor algorithm for the values of tolerance and arithmetical precisionpointed out in Subsection IV A, to determine the images V [ k ] of the vertices V [ k ]0 , k = 0 , . . . , . , . , . T = 54400 seconds, the 1–ring[ V [0] ; V [1] , . . . , V [6] ] provides a discretized neighborhood of the point V [0] ∈ S T (see Figure 3).Applying the scheme developed in Subsection III C we obtain that the mean and therelative total curvatures at point V [0] take, for both the different values of the radius r and the tolerances in the Taylor algorithm, the values shown in Table IV B 1, where oneobserves that for values
Tol = 1 . E −
20 and r = R (cid:12) /
50 the three first significant decimalnumbers are correct.Table: Variation of the relative total curvature and mean curvature with respect to boththe radius of the 1-ring and the tolerance used in the numerical integration of the rayequation. 12OL=1 E −
10 TOL=1 E − K rel H K rel HR (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . − R (cid:12) / − . −
10 0 . − − . −
10 0 . −
2. Wavefront surface
To obtain an estimation of the curvature of a region of the wavefront propagating alongthe Oz –axis, we apply the ray tracing method with initial values on the wavefront surfacedescribed at the beginning of Subsection IV B. We use a regular discretization of the wave-front surface as that described in Subsection III A; Furthermore, taking into account theresults shown in Table IV B 1 we choose the length of the edges in the corresponding meshequal to R (cid:12) / R (cid:12) /
25 and 2 R (cid:12) respec-tively. Therefore a number of N := 28044 vertices is required. The Taylor algorithm withprojection in a time interval [0 , T ] (the time required to run a path of length equal to 101AU) is applied to each photon located at an initial vertex; the CPU-time employed to carryout this computation is of 1544 seconds.For the numerical solution u n , n = 1 , . . . , N , of (24), both the mean and relative-totalcurvatures at each inner vertex of the mesh on the surface S T are computed by applyingthe method described in Section III, schematized in pseudocode in the next Table.13 ata: u ∗ n := ( z ∗ n , ˙ z ∗ n ) , n = 1 , . . . N for n = 1 . . . N do u n := Taylor ( t, u ∗ n ) // rays tracing y n := NormalCoordinates ( z n ) // see Eq. (11) for i = 0 . . . do y n i := Ring ( y n ) // see Eq. (16) end ( a , a , a ) = LeastSquares ( y n i ) // see Eq. (18) γ AB ( x n ) := Metric ( a , a , a , x n ) // see Eq. (19) B = SecondFundamentalForm ( x n ) // see Eq. (20)( λ , λ ) = Diagonalize ( B )( K rel , H ) = Curvature ( λ , λ ) // see Eqs. (21,22) end In Figure 4, the surface S T at the time when the wavefront arrives at the Earth is shownusing a gray-scale to represent the mean curvature (note we have used a different scale onthe Oz –axis). One sees in this figure that the absolute value of the mean curvature functiondefined on S T increases as the distance between the photon and the Oz –axis decreases. FIG. 4: Wavefront surface and relative total curvature (gray scale) deformed by a spherical gravi-tational field (a different scale is used for the vertical axis). . Curvature of a wavefront in the PPN formalism To derive the deflection angle for light rays, instead of assuming the validity of generalrelativity, one may consider a more general expression for the metric generated by a sphericalcentral body that is valid for different gravitational metric theories. In the parameterizedpost-Newtonian formalism, the expression of a spherically symmetric metric, written atthe linearized order we are considering here, contains one parameter γ which is usuallyinterpreted as a modification of the curvature of the space. In the parameterized post-Newtonian formalism the total relativistic deflection angle of light rays passing near thelimb of the Sun is given approximately as (see [21]):∆ φ (cid:39) γ )1 (cid:48)(cid:48) . . (31)Using very long baseline interferometry techniques, one may obtain high precision generalrelativistic measurements for the deflection of radio signals from distant radio sources withan accuracy at the 0 .
02 percent level. In [22] an estimation of 0 . ± . γ .Here we consider a gravitational field depending on the Eddington parameter γ andadapt the numerical method described above to determine the dependence of the wavefrontcurvature at a point, corresponding to a light ray grazing the Sun with respect to thisparameter. For the numerical discussion, we take a gravitational field in which only thedominant terms in the post-Newtonian metric are included, so that the metric deviationsmay be written as h ab = 2 c − κγ (cid:107) z (cid:107) − δ ab , h a = 0 , h = 2 c − κ (cid:107) z (cid:107) − . (32)We take values of γ in the interval [0 . , . γ . For these values of γ thevariations of K rel and H are in the second and third significant decimal digit respectively.In consequence, the numerical treatment we applied gives account of the variations of thegeometric model when a change of the PPN Eddington parameter is made.15 IG. 5: Relative total curvature K rel (a) and mean curvature H (b) versus post-Newtonian pa-rameter γ for a wavefront numerically determined. V. VARIATION OF THE TIME OF ARRIVAL AND CURVATURE OF THEWAVEFRONT
In this section we consider the problem treated by Samuel [1] where the curvature of awavefront in a gravitational field is detected by measuring the arrival times T a , a = 0 , , , E a on an Earth hemisphere. The arrival timedifferences between these stations depend on the curvature of the wavefront and are relatedto the volume of a parallelepiped formed by the vectors −−−→ E E j − c ( T j − T ) n j , j = 1 , , n := (0 , , A. Volume of a tetrahedron with vertices on the wavefront surface
Here we compute the volume and the arrival time differences associated with a simplexdetermined by four points on the wavefront numerically determined above. For this, weconsider a light ray reaching the reference station E (see Figure 6) at the instant of time T and whose tangent vector in this point is n . Let P ∗ be the position of this photon ata previous instant T ∗ such that c ( T − T ∗ ) = 101 AU. This position may be obtained bya backward integration of the equations of motion (24) with initial values ( P , − c n ). On16he plane determined by ( P ∗ , n ∗ ), where n ∗ is the tangent vector to the ray at the point P ∗ , we consider the l–ring [ P ∗ ; P ∗ j ], j = 1 , . . . ,
6, centered at P ∗ . Then, from (24) and theinitial values ( P ∗ j , c n ∗ ) we determine the 1–ring [ Q ; Q j ] image of [ P ∗ ; P ∗ j ] by the flow (10)integrating again over the interval [0 , T − T ∗ ].In the numerical model we are considering, the coordinates of P are (2 R (cid:12) , , n = (0 , , P ∗ and direction n ∗ of the photon are then (expressedhere using only three decimal digits): P ∗ := (1 . , . , − . , n ∗ = (0 .
183 E − , . , . . (33)We choose as vertices P ∗ j the points P ∗ j = R ( P ∗ , e ,ϑ ) (cid:0) r cos( πj/ , r sen( πj/ , − AU (cid:1) j = 1 , . . . , R ( P ∗ , e ,ϑ ) represents a rotation around axis ( P ∗ , e ) through an angle ϑ := e · n ∗ / (cid:107) n ∗ (cid:107) .Following the same method employed in Subsection IV B 2 we perform a least-squaresfitting of the data { Q a } a =0 to get an approximation of the wavefront surface in a neighbor-hood of Q given by a quadratic surface Q which, in the normal coordinates correspondingto ray ( P , n ), takes the form Q ( z ) : z = − . × − z + 0 . × − z z + 0 . × − z . (35)Now, on the Earth surface we choose four points E a with geocentric coordinates(0 , , − R ⊕ ) , ( R ⊕ , , , ( √ R ⊕ , R ⊕ ,
0) and ( √ R ⊕ , − R ⊕ ,
0) respectively, which will betransformed to normal coordinates. Assuming that the geometry of the 3–space in thevicinity of the Earth is Euclidean, one may determine the points Q a ∈ Q whose distancesto the corresponding stations E a are minima.Once the points Q a ∈ Q corresponding to the images of the stations E a on the Earthare determined, we may compute the differences τ a = dist( Q a , Q (cid:48) a ) /c between the arrivaltimes measured at each station under the assumptions that the wavefront surface is either acurved surface Q or a plane P determined by the pair ( P , n ). These differences are given,for the data provided above, by τ = 0 . , τ = 0 . , τ = τ = 0 . , (36)17hereas the distances l ij := dist( Q i , Q j ), i, j = 0 , , ,
3, between the projected points Q a on Q are (in the normalized units we are using) l = 0 . − , l = 0 . − , l = 0 . − ,l = 0 . − , l = 0 . − l = 0 . − A, B ) denotes the Euclidian distance between two arbitrary points A and B and Q (cid:48) a represents the points of minimum distance from the station E a to the plane P . Thevolume of the tetrahedron determined by the points Q a is proportional to the Cayley-Mengerdeterminant [12] defined in terms of the lengths l ij of the edges and it is given by:Vol = √ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) l l l l l l l l l l l l (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (37)which for the lengths (37) takes the value Vol = 0 . −
7. Therefore the metric quadruple( Q a , l ab ) determines a non degenerate simplex in R (see for instance Saucan [17]). This isequivalent to say that the metric quadruple ( Q a , l ab ) is not congruent to any quadruple ofpoints in the Euclidean plane and consequently the curvature of the wavefront surface atthe point Q is non vanishing. B. Estimation of the wavefront curvature from arrival time measurements
Now we assume that the points { Q a } a =0 on the wavefront are directly determined throughmeasurements of arrival times. Here we consider the problem of determining an approxima-tion of the wavefront curvature in a region far enough from the Sun (say the Earth), withoutresorting to the ray tracing method.An estimation of the Gaussian curvature of the wavefront surface can be obtained usingthe notion of the Wald curvature of a metric space established in the Distance Geometry (see[12]), that in the case of 2-dimensional manifolds agrees with the Gaussian curvature. TheWald curvature is determined as the limit of the embedding curvatures of metric quadruplesisometrically embedded in surfaces of constant curvature (the Euclidean plane R , the 2–sphere S √ κ or the hyperbolic space H √− κ ). 18 IG. 6: Schematic diagram of a wavefront plane (line P P ) and curve (line Q Q ). Continuouslines E i Q i represent the light rays associated with Q Q while the dashed lines E i P i represent theassociated light rays to P P . In the hyperbolic plane H r of curvature − /r , represented by the Blumenthal model([12], p. 19) we consider the metric quadruple ( Q j , l ij ) defined in Subsection V A. Thecurvature of a hyperbolic plane on which there exists a quadruple congruent with Q j mustfulfill both of the following conditions: A ( r ) := l /r ) cosh( l /r ) cosh( l /r )cosh( l /r ) 1 cosh( l /r ) cosh( l /r )cosh( l /r ) cosh( l /r ) 1 cosh( l /r )cosh( l /r ) cosh( l /r ) cosh( l /r ) 1 = 0 , (38)and each non-zero principal minor of A ( r ) of order m +1 has the sign ( − m (see [12], p. 274).For small values of the curvature 1 /r the determinant det (cid:0) A ( r ) (cid:1) can be approximated by aTaylor polynomial. Using the symbolic processor Maple and employing numerical precisionof Digits =50 to perform the Taylor expansion we obtain that (38) may be approximatedby the following algebraic equation for ρ := 1 /r . × − ρ + 0 . × − ρ − . × − ρ = 0 (39)having only one positive real root. This approximated solution is taken as an initial guessedsolution to solve the transcendental equation (38) by means of the command solve in19aple. On the other hand, the conditions imposed on the sign of the principal minors arealso satisfied. Therefore the Wald curvature K W associated with the chosen quadruple { Q j } is K W := − r = − .
16 E − . (40)This result gives an approximation of the total curvature of the wavefront surface under theassumption that locally this surface may be identified with a hyperbolic plane in which thequadruple considered is isometrically embedded. C. Scheme of the method
In this subsection we present an outline of the construction of the method we used aboveto obtain the Wald curvature corresponding to four receiving stations located at points E a and the arrival time differences τ a :– Define on the wavefront the points P a corresponding to the stations E a : −−→ OQ a = −−→ OE a + cτ a n , (41)where O is the coordinate origin and n is the unit vector in the direction of the lightrays.– Determine the relative distances l ij between the points Q a : l ij = L ij + 2 c ( τ j − τ i )( −−→ OE j − −−→ OE i ) · n + c ( τ j − τ i ) , (42)where L ij are the Euclidean distances between the stations E i and E j .– For the point Q a so obtained, establish the nonlinear equation (38) in terms of thecurvature of the surface.– Solve the nonlinear equation (38) for the unknown r to obtain an estimation of thecurvature of the wavefront in a neighborhood of the station E by means of the Waldcurvature K W of this surface modeled as a hyperbolic plane.20 I. CONCLUSIONS
The ray tracing numerical method provides a useful tool for the description of spacelikebidimensional wavefronts within the framework of the general relativity. We have studieda method, based on techniques of computational geometry, that allows to estimate thecurvature properties of the surface by making a least-squares fitting of the wavefront surfaceby a quadric surface in the neighborhood of each point of this surface. The computation ofthe light rays is carried out using an algorithm based on the Taylor method for the solutionof differential equations and employing high arithmetic precision. Further, we have applieda projection at each step of the numerical integration process that allows to guarantee thefulfillment (at machine precision) of the isotropy condition for the tangent vector to thelight ray. We have also studied numerically the dependence of the curvature properties ofthe wavefront surface on the value of the Eddington parameter γ . On the other hand, wehave employed a geometric computational approach to the study of the model proposed bySamuel as a new general relativity test, by determining a numerical approximation of thevolume corresponding to a tetrahedron formed by four points on the wavefront that reachesfour receiving stations on the Earth surface. Finally, we have obtained an estimation of theWald curvature for the wavefront in a vicinity of the Earth by using the differences of arrivaltime recorded at four receiving stations on the Earth. VII. ACKNOWLEDGEMENTS
This research was supported by the Spanish Ministerio de Educaci´on y Ciencia, MEC-FEDER Project ESP2006-01263. [1] Samuel J 2004
Class. Quantum Grav. , , L83–L88.[2] Klioner S A and Peip M 2003 Astron. Astrophys. , 1063–1074.[3] de Felice F, Crosta M T, Vecchiato A, Lattanzi M G and Bucciarelli B 2004
The AstronomicalJournal , 580–595.[4] Kopeikin S and Sch¨afer G 1999
Phys. Rev. D , 124002.[5] Le Poncin-Lafitte C, Linet B and Teyssandier P 2004 Class. Quantum Grav. , , 4463-4483.
6] Teyssandier P and Le Poncin-Lafitte C 2008
Class. Quantum Grav. , , 145020.[7] Garimella R V and Swartz BK 2003 Technical Report, LA-UR-03-8240, Los Alamos NationalLaboratory.[8] Cazals F and Pouget M 2003 SGP ’03: Proceedings of the 2003 Eurographics/ACM SIG-GRAPH symposium on Geometry processing (Aachen, Germany) , (Aire-la-Ville: Eurograph-ics Association) p. 177.[9] Petitjean S 2002
ACM Computing Surveys , 1–6.[10] Jorba `A and Zou M 2005 Experimental Mathematics , 99-117.[11] Landau L D and Lifshitz E M 1962 The Classical Theory of Fields (Oxford: Pergamon Press).[12] Blumenthal L M 1970
Theory and Applications of Distance Geometry . (New York: ChelseaPublishing Company, 2nd edition).[13] Brumberg V A 1991
Essential Relativistic Celestial Mechanics , (Bristol: Adam Hilger).[14] Bobenko A I and Hoffmann T 2003
Duke Math. J.
Riemannian Geometry (Princeton: Princeton University Press).[16] Brewin L 1998
Class. Quantum Grav . What is Geometry? , edG Sica, ( Monza/Italy: Polimetrica International Scientific Publisher).[18] Meyer M, Desbrum M, Schr¨oder P and Barr AH (2002)
Proceedings of VisMath ’02 (Berlin,Germany).[19] do Carmo M P 1992
Riemannian Geometry (Boston: Birkh¨auser).[20] Hairer E 2000
BIT , , 726-734.[21] Misner C W, Thorne K S and Wheeler J A 1973. Gravitation , (San Francisco: Freeman).[22] Shapiro S S, Davis J L, Lebach D E and Gregory J S 2004
Phys. Rev. Lett , 121101., 121101.