Calculation of the amplitudes of elastic waves in anisotropic media in Cartesian or ray-centred coordinates

We derive various expressions for the amplitude of the ray-theory approximation of elastic waves in heterogeneous anisotropic media, and show their mutual relations. The amplitude of a wavefield with general initial conditions can be expressed in terms of two paraxial vectors of geometrical spreading in Cartesian coordinates, and in terms of the 2×2 matrix of geometrical spreading in ray-centred coordinates. The amplitude of the Green tensor can be expressed in six different ways: (a) in terms of the paraxial vectors corresponding to two ray parameters in Cartesian coordinates, (b) in terms of the 2×2 paraxial matrices corresponding to two ray parameters in ray-centred coordinates, (c) in terms of the 3×3 upper right submatrix of the 6×6 propagator matrix of geodesic deviation in Cartesian coordinates, (d) in terms of the 2×2 upper right submatrix of the 4×4 propagator matrix of geodesic deviation in ray-centred coordinates, (e) in terms of the 3×3 matrix of the mixed second-order spatial derivatives of the characteristic function with respect to the source and receiver Cartesian coordinates, and (f) in terms of the 2×2 matrix of the mixed second-order spatial derivatives of the characteristic function with respect to the source and receiver ray-centred coordinates. The step-by-step derivation of various equivalent expressions, both known or novel, elucidates the mutual relations between these expressions.


INTRODUCTION
The most important quantity used in ray methods is the travel time defined by means of the non-linear first-order partial differential equation called the Hamilton-Jacobi equation, which is also often referred to as the eikonal equation in wave propagation problems. The next most important quantity is the zero-order raytheory amplitude defined by means of the transport equation.
We may also calculate the propagator matrix of geodesic deviation. The propagator matrix may be defined and calculated in Cartesian coordinates or in raycentred coordinates. If we calculate the propagator matrix of geodesic deviation, we may obtain the paraxial vectors using the corresponding initial conditions and the propagator matrix. The paraxial vectors may then be used to determine the amplitude of a general wavefield or the amplitude of the Green tensor as mentioned above.
If we calculate the propagator matrix of geodesic deviation, we may also determine the amplitude of the Green tensor directly from the upper right submatrix of the propagator matrix of geodesic deviation, either from the 3×3 submatrix in Cartesian coordinates (Section 5.4) or from the 2×2 submatrix in ray-centred coordinates (Section 5.1).
We may also consider the mixed second-order spatial derivatives of the characteristic function with respect to the source and receiver coordinates. In this case, we may determine the amplitude of the Green tensor using the 3×3 matrix of these derivatives in Cartesian coordinates (Section 5.6) or the 2×2 matrix of these derivatives in ray-centred coordinates (Section 5.5).
We do not present the expressions for the amplitude in terms of the surface-tosurface paraxial vectors of geometrical spreading corresponding to two ray parameters, in terms of the 2×2 surface-to-surface paraxial matrices, or in terms of the upper right 2×2 submatrix of the 4×4 surface-to-surface propagator matrix of geodesic deviation, refer to Moser andČervený (2007). Nor do we demonstrate that the expression for the amplitude of the Green tensor in terms of the three different 2×2 matrices of the homogeneous second-order spatial derivatives of travel time and the characteristic function with respect to the source and receiver ray-centred coordinates byČervený (2001, Eq. 4.10.43) is applicable to anisotropic media.
The Einstein summation over repetitive lower-case Roman indices, corresponding to the 3 spatial coordinates, is used throughout the paper. Upper-case Roman indices correspond to the first two ray-centred coordinates or to two ray parameters.
where the function H(x i , p j ) of coordinates x i and of covariant vector p j from the cotangent space at point x i is referred to as the Hamiltonian function. The multivalued solution τ (x m ) of the Hamilton-Jacobi equation is determined by the initial conditions.

L. Klimeš
The characteristic function (two-point travel time, point-to-point distance) from pointx n to point x m , denoted by τ (x m ,x n ), is the solution of Hamilton-Jacobi equation (1) with initial conditions τ (x m ,x n ) = 0. Gradient of the travel time or of the characteristic function is referred to as the slowness vector.
When differentiating the Hamilton-Jacobi equation with respect to coordinates x j , we find that the multivalued solution τ of the Hamilton-Jacobi equation can be calculated along rays (geodesics). Hamilton's equations (equations of rays, ray tracing equations, equations of geodesics) read where d in the derivatives denotes differentiation along the ray. The meaning of independent parameter γ 3 along the rays depends on the particular form of the Hamiltonian function.
In addition to independent parameter γ 3 , corresponding to the Hamiltonian function, we may also parametrize the rays by travel time τ , and define the ray velocity vector as the derivative of coordinates x i of the ray with respect to travel time τ along the ray.
In 3-D space, the orthonomic system of rays corresponding to the given initial conditions consists of a two-parametric system of rays which are parametrized by two ray parameters γ 1 and γ 2 . The ray parameters together with the independent parameter along the rays form ray coordinates γ 1 , γ 2 , γ 3 .
Since the solution of Hamilton-Jacobi equation (1) calculated using Hamilton's equations (3) and (4) is multivalued, it is parametrized by ray coordinates γ a . Multivalued solution τ (x m ) or τ (x m ,x n ) is thus expressed as x m = x m (γ a ), τ = τ (γ a ).

. . R a y -c e n t r e d c o o r d i n a t e s
We may define ray-centred coordinates q a along a particular ray (Klimeš, 2006b). We parametrize the points along the ray by an arbitrary monotonic variable q 3 . At each point x i (q 3 ) of the ray, we choose two contravariant basis vectors h i 1 (q 3 ) and h i 2 (q 3 ) perpendicular to slowness vector p i , Contravariant basis vectors h i A should vary smoothly along the ray. The transformation from the ray-centred coordinates q a to Cartesian coordinates x i is defined by relation (7) Three contravariant basis vectors of the ray-centred coordinate system are In matrix notation, we shall denote the first two contravariant basis vectors h i 1 and h i 2 as h 1 and h 2 . The slowness vector in ray-centred coordinates at the central ray then reads p (q) a = p i h i a .
(9) Three covariant basis vectors of the ray-centred coordinate system arê In matrix notation, we shall denote the first two covariant basis vectorsĥ 1 i andĥ 2 i asĥ 1 andĥ 2 .
2 . 3 . P a r a x i a l r a y m a t r i c e s o f a n o r t h o n o m i c s y s t e m o f r a y s For an orthonomic system of rays corresponding to the given initial conditions, we define the 3×3 paraxial matrix of geometrical spreading in Cartesian coordinates. We analogously define the 3×3 paraxial matrix describing the paraxial slowness vectors in Cartesian coordinates. The third columns of the paraxial matrices can be obtained from the solution of Hamilton's equations (3) and (4). We shall refer to the first two columns X i 1 , X i 2 or Y i1 , Y i2 of the paraxial matrices as the paraxial vectors, and denote them as X 1 , X 2 or Y 1 , Y 2 in matrix notation. The paraxial vectors can be calculated using the Hamiltonian equations of geodesic deviation (paraxial ray equations, dynamic ray tracing equations), which can be obtained by differentiating Hamilton's equations (3) and (4) with respect to ray coordinates γ a (Červený, 1972 ). We also define the analogous paraxial matrices and in ray-centred coordinates.
In ray-centred coordinates, we can thus calculate just the 2×2 paraxial matrices Q I A and P IA using the Hamiltonian equations of geodesic deviation. In matrix notation, we shall denote these 2×2 paraxial matrices Q and P. We compare definitions (5) and (11), and obtain relation Definition (13) yields identity

. 4 . P r o p a g a t o r m a t r i c e s o f g e o d e s i c d e v i a t i o n
The propagator matrices of geodesic deviation represent the solution of the Hamiltonian equations of geodesic deviation with identity initial conditions, both in Cartesian coordinates and in ray-centred coordinates (Kendall et al., 1992 ;Klimeš, 1994 ).
The 6×6 propagator matrix of geodesic deviation in Cartesian coordinates may also be defined as the matrix of the derivatives of x i and p i with respect to their initial conditionsx j andp j for fixed γ 3 . In this paper, we shall use just the 3×3 upper right submatrix of this propagator matrix. The partial derivatives are calculated for fixed γ 3 . The 6×6 propagator matrix of geodesic deviation in ray-centred coordinates may also be defined as the matrix of the derivatives of q i and p (q) i with respect to their initial conditionsq j andp (q) j for fixed γ 3 . As the consequence of identities (15) and (16), we may also define the 4×4 propagator matrix of geodesic deviation in ray-centred coordinates as the matrix of the derivatives of q I and p (q) I with respect to their initial conditionsq J andp (q) J for fixed γ 3 . In this paper, we shall use just the upper right 2×2 submatrix of this propagator matrix. In matrix notation, we shall denote this 2×2 matrix as Q 2 .

. 1 . T r a n s p o r t e q u a t i o n
Multivalued zero-order ray-theory amplitude A = A(x m ) of a general elastic wavefield satisfies transport equation (Klimeš, 2006a, Eq. 10 ), where ray velocity vector V i is given by definition (5).
is a function parametrizing the transport equation. If A is the amplitude of the displacement of an elastic wavefield, ̺ is the density. From the point of view of differential geometry, amplitude A is a scalar if ̺ is the scalar density of weight −1 (scalar per volume).

. 2 . P h a s e s h i f t d u e t o c a u s t i c s
Transport equation (21) is a partial differential equation for the square A 2 of the amplitude, not for the amplitude itself. Even if the solution A 2 of transport equation (21) is real-valued, amplitude A becomes complex-valued if its square A 2 becomes negative. Amplitude A is thus complex-valued. Since the complex-valued square root has two branches, it is difficult to determine amplitude A from its square A 2 . We must determine which branch of the amplitude calculated along the ray is correct.
We thus separate square root A = √ A 2 into complex modulus |A| = |A 2 | and complex argument exp(i ϕ), Quantity ϕ in expression (22) is the phase shift due to caustics. If the rays are real-valued, there is a real-valued solution of transport equation (21) for the square A 2 of the amplitude. The phase shift due to caustics is then constant in the regions of equal sign of A 2 , and changes by an integer multiple of π/2 at the boundaries between the regions.
The calculation of amplitude A along the ray is thus composed of the calculation of its complex modulus |A| along the ray and of the determination of the phase shift ϕ due to caustics. This paper is devoted to the the calculation of the complex modulus |A| of the amplitude.
For an orthonomic system of rays corresponding to a wavefield with general initial conditions, the phase shift due to caustics may be determined using the 3×3 paraxial matrices in Cartesian coordinates with their derivatives (Klimeš, 2014, Secs 2.2 and 3.2 ), or using the 2×2 paraxial matrices in ray-centred coordinates with their derivatives (Klimeš, 2014, Secs 2.1 and 3.1 ).

L. Klimeš
The phase shift of the Green tensor due to caustics may be determined using the 3×3 right-hand submatrices of the 6×6 propagator matrix of geodesic deviation in Cartesian coordinates with their derivatives (Klimeš, 2010, Secs 2.2 and 3.2 ), or using the 2×2 right-hand submatrices of the 4×4 propagator matrix of geodesic deviation in ray-centred coordinates with their derivatives (Klimeš, 2010, Secs 2.1 and 3.1 ). The solution of transport equation (21) may be expressed in various forms. The square of the complex-valued amplitude may be expressed as (Klimeš, 2006a, Eqs 12, 34 ), where the 3×3 paraxial matrix X i a of geometrical spreading in Cartesian coordinates is given by definition (11).

AMPLITUDE OF
This equation represents the generalization of the analogous equation by Babich (1961, Eq. 3.7) from a homogeneous Hamiltonian function to a general Hamiltonian function.
The complex-valued amplitude then reads where ϕ is the phase shift due to caustics. Complex-valued factor C = C(γ 1 , γ 2 ) is constant along the ray in a smooth medium and is determined by the initial conditions. It is often referred to as the reduced amplitude (Červený et al., 1988, Eq. 5.19 ). Reduced amplitude C depends on the selection of ray parameters γ 1 and γ 2 .
For the special case of a homogeneous Hamiltonian function, paraxial vectors X 1 and X 2 are tangent to the wavefront, and their cross product is thus normal to the wavefront, In this case, (Červený, 1972, Eq. 29b;Kendall and Thomson, 1989, Eqs 29-30 ). This equation is not applicable to a general Hamiltonian function.
for the transformation from ray-centred coordinates to Cartesian coordinates. Considering definitions (8) and (13), relation (28) reads which implies relation for the determinants. The determinant of the transformation matrix (8) from raycentred to Cartesian coordinates is Since contravariant basis vectors h 1 and h 2 of the ray-centred coordinate system are tangent to the wavefront, their cross product is normal to the wavefront, and We insert relation following from definitions (7) and (8), into relation (31) with (32) and obtain relation Relation (30) with relation (34) reads Due to identities (15) and (18), we have relation We insert relation (35) with relation (36) into expression (24) and arrive at expression (Klimeš, 2012, Eqs 7, 9 ) for the amplitude of a general wavefield.

. 1 . A m p l i t u d e i n t e r m s o f t h e p r o p a g a t o r m a t r i x o f g e o d e s i c d e v i a t i o n i n r a y -c e n t r e d c o o r d i n a t e s
Using the representation theorem for elastic waves and relation (37) in raycentred coordinates, we can derive expression (Červený, 2001, Eq. 5.4.24 ;Klimeš, 2012, Eq. 55 ) for the amplitude of the Green tensor from pointx to point x in the frequency domain. Here (Klimeš, 2012, Eq. 13 ) is the relative geometrical spreading defined byČervený (2001, Eq. 4.14.45). Amplitude (38) corresponds to the Fourier transform of the Green tensor from time t to circular frequency ω. If the right-hand side of the Fourier transform includes multiplicative factor (2π) − 1 2 or (2π) −1 , the right-hand side of expression (38) should be multiplied by the same factor.
The determinant of the transformation matrix (10) from Cartesian to ray-centred coordinates is | det(ĥ i a )| = |ε abcĥ1 Since the covariant basis vectorsĥ 1 andĥ 2 of the ray-centred coordinate system are perpendicular to the ray, their cross product is tangent to the ray, and ε abcĥ1 analogous to relation (32) for the contravariant basis vectors. Here V is the ray velocity and V c is the ray-velocity vector. We insert definition (5) and part of definition (10) into relation (41) with (42) and obtain relation Inserting (42) into (41) and considering (44), we arrive at Since relations (34) and (45) yield identity which can be inserted into expression (39), both at pointx or point x, e.g.,

. 2 . A m p l i t u d e i n t e r m s o f t h e p a r a x i a l m a t r i c e s i n r a y -c e n t r e d c o o r d i n a t e s
Paraxial matrices Q(x) and P(x) corresponding to arbitrarily parametrized rays from a point source atx are related by equation Inserting relation (50) into expression (39), we obtain expression for the relative geometrical spreading. We may insert identity (47) into expression (51), both at pointx or point x, e.g., or

. 3 . A m p l i t u d e i n t e r m s o f t h e p a r a x i a l v e c t o r s i n C a r t e s i a n c o o r d i n a t e s
We supplement paraxial vectors Y i2 and Y i1 with slowness vector p i to create a 3×3 matrix. The transformation of this matrix from ray-centred coordinates to Cartesian coordinates reads It follows that the corresponding determinants are transformed as We insert relation (45) into relation (55) and obtain Stud. Geophys. Geod., 63 (2019)

. 4 . A m p l i t u d e i n t e r m s o f t h e p r o p a g a t o r m a t r i x o f g e o d e s i c d e v i a t i o n i n C a r t e s i a n c o o r d i n a t e s
Paraxial vectors X i A (x) and Y mA (x) corresponding to arbitrarily parametrized rays from a point source atx are related by equation where X im 2 (x,x) is the 3×3 upper right submatrix (19) of the 6×6 propagator matrix of geodesic deviation in Cartesian coordinates. Then We consider the skewness of Levi-Civita symbol ε ijk , and rewrite relation (64) to read We now insert identity ε mnl ε lrs = δ mr δ ns − δ ms δ nr (66) into relation (65), and arrive at (67) We insert relation (59) into relation (67), and obtain relation of the cofactors of matrix X im 2 (x,x), insert relation (68) with definition (69) into expression (61), and obtain relation (Kendall et al., 1992, Eq. 17b;Chapman, 2004, Eq. 5.4.23 ) for the relative geometrical spreading.
between the mixed second-order spatial derivatives of characteristic function τ (x, x ′ ) and the 3×3 upper right submatrix of the 6×6 propagator matrix of geodesic deviation in general coordinates including Cartesian coordinates. The meaning of functions γ(x,x) and Γ(x,x) is not significant here. Interested readers may refer to Klimeš (2013a,b).
If we define matrix [X −1 2 ] ji (x,x) inverse to matrix X ij 2 (x,x), we may express relation (71) as We now transform relation (72) into ray-centred coordinates. We transform the submatrix X ij 2 (x,x) of the propagator matrix of geodesic deviation from Cartesian to ray-centred coordinates, Stud. Geophys. Geod., 63 (2019)

L. Klimeš
The transformation of the inverse matrix then reads We transform the mixed second-order derivatives of the characteristic function from Cartesian to ray-centred coordinates, which can be expressed in terms of the contravariant basis vectors (8) of the raycentred coordinate system as Since (Hamilton, 1837, Eqs U, I ), and (Hamilton, 1837, Eqs Y, I ), we have identities Relation (72) in ray-centred coordinates then reads We define 2×2 matrix invert the 3×3 matrices in relation (80), and arrive at We see that the 2×2 matrix Q BA 2 (x,x) is a submatrix of the 3×3 matrix (73). Definition (81) may then be understood as the relation between the 2×2 submatrix Q BA 2 (x,x) of the 3×3 matrix (73) and the mixed second-order spatial derivatives of characteristic function τ (x, x ′ ) in ray-centred coordinates. This relation represents a generalization of the relation for the mixed second-order spatial derivatives of characteristic function τ (x, x ′ ) byČervený et al. (1984, Eq. 22) to anisotropic media.
We insert relation (83) into expression (39) and obtain expression for the relative geometrical spreading in terms of the mixed second-order spatial derivatives of the characteristic function in ray-centred coordinates. The special case of expression (84), corresponding to orthonormal contravariant basis vectors h 1 and h 2 , was presented by Schleicher et al. (2001, Eqs 5-6).
Expression (84) for the relative geometrical spreading may also be modified by inserting identity (47), both at pointx or point x, e.g., The transformation of the mixed second-order derivatives of the characteristic function from ray-centred to Cartesian coordinates reads Considering identities (79), we express relation (86) in terms of the covariant basis vectors (10) of the ray-centred coordinate system as We define the matrix of the cofactors of matrix (86). We insert relation (87) into definition (88) and arrive at relation We insert relation (42) into relation (89) and obtain relation We multiply relation (90) by p i (x) and p l (x), insert the product into expression (85), and obtain expression for the relative geometrical spreading in terms of the mixed second-order derivatives of the characteristic function in Cartesian coordinates. We insert expression (91) for the relative geometrical spreading into expression (38) and obtain new expression for the amplitude of the Green tensor.

CONCLUSIONS
The zero-order ray-theory amplitude satisfies the transport equation. There are two important kinds of the ray-theory amplitude: the amplitude corresponding to a wavefield with general initial conditions including finite sources, and the amplitude of the Green tensor. Thre already known equivalent expressions (24), (25) and (37) for the amplitude corresponding to the general initial conditions are presented in Section 4. Relation (38) for the amplitude of the Green tensor is presented in Section 5 together with twelve equivalent expressions for the relative geometrical spreading. Expressions (39), (58), (61) and (70) (91) for the relative geometrical spreading may be novel, although some special cases of some of these expressions have been published. The step-by-step derivation of various equivalent expressions starts with expressions (24) and (25), or expressions (38) and (39), respectively. The derivation thus elucidates the mutual relations between the equivalent expressions.
The Hamiltonian equations of geodesic deviation (paraxial ray equations, dynamic ray tracing equations) are sometimes expressed and solved in general Cartesian coordinates and sometimes in ray-centred coordinates connected with a particular ray. The expressions for the zero-order ray-theory amplitude in terms of the results of the Hamiltonian equations of geodesic deviation in general Cartesian coordinates are presented in Sections 4.1, 5.3, 5.4 and 5.6. The expressions for the zero-order ray-theory amplitude in terms of the results of the Hamiltonian equations of geodesic deviation in ray-centred coordinates are presented in Sections 4.2, 5.1, 5.2 and 5.5.
The zero-order ray-theory amplitude of the Green tensor can be expressed in terms of the propagator matrix of the Hamiltonian equations of geodesic deviation (Sections 5.1 and 5.4), in terms of the paraxial vectors corresponding to arbitrarily parametrized rays from a point source (Sections 5.2 and 5.3), or in terms of the second-order derivatives of the characteristic function (Sections 5.5 and 5.6).