Superpositions of Gaussian beams and column Gaussian packets in heterogeneous anisotropic media

S U M M A R Y Based on the integral superposition of Gaussian packets, we derive the equations for the integral superposition of Gaussian beams and for the integral superposition of column Gaussian packets in smoothly heterogeneous media. Whereas Gaussian beams extend along their central rays, column Gaussian packets extend along an arbitrary system of lines. The equations are applicable to both the anisotropic ray theory and the coupling ray theory in anisotropic media, or to the isotropic ray theory in isotropic media. The superpositions corresponding to the coupling ray theory can be applied to various kinds of reference rays. The equations can be used in both Cartesian and curvilinear coordinates.


I N T RO D U C T I O N
A Gaussian beam or packet is a localized solution of a wave equation under consideration, expressed in terms of a complex-valued vectorial amplitude and complex-valued traveltime.Gaussian beams and packets may serve as building blocks of a wavefield (Ralston 1983).A Gaussian beam is a frequency-domain solution localized in a vicinity of its central ray.A Gaussian packet is a time-domain solution localized in a vicinity of its central point moving along the central ray.The integral superposition of Gaussian beams is performed over all central rays of Gaussian beams.The integral superposition of Gaussian packets is performed over all central points of Gaussian packets and yields the frequency-domain solution corresponding to the reference frequency of the superposed time-domain Gaussian packets.The integral superposition of Gaussian beams or packets overcomes the problems of the standard ray theory with caustics.The integral superposition of Gaussian beams and packets is considerably comprehensive and flexible.It may be formulated in many ways, including the Maslov method and its various generalizations as special cases.The form of the integral superposition depends primarily on the specification of the wavefield, also comprising the system of Gaussian packets scattered from Gabor functions forming medium perturbations.In this paper, we concentrate on the representation of non-directional wavefields without pronounced diffractions rather than on the decomposition of directional beams or of diffracted waves.
The equations for the integral superposition of Gaussian packets in heterogeneous isotropic media were derived using the Maslov asymptotic theory applied to a general 3-D subspace of the 6-D complex phase space by Klimeš (1984b), and demonstrated in the computation of seismic wavefields by Klimeš (1989).These equations were rederived for heterogeneous anisotropic media by Klimeš (2014Klimeš ( , 2018)).Červený & Pšenčík (2015) generalized the integral superposition of Gaussian beams to the integral superposition of column Gaussian packets, which infinitely extend along the reference lines different from central rays.The two main objectives of this paper are to rederive the equations for the integral superposition of Gaussian beams for heterogeneous anisotropic media in both general and ray-centred coordinates, and to compare the integral superposition of column Gaussian packets with the integral superposition of Gaussian beams.
In this paper, we start with the three-parametric superposition of Gaussian packets in a smoothly heterogeneous generally anisotropic medium by Klimeš (2018).We consider the endpoints of central rays along a reference surface and a system of reference lines intersecting the reference surface.We then asymptotically calculate the one-parametric superpositions along the reference lines and in this way convert the three-parametric superposition of Gaussian packets into a two-parametric superposition of column Gaussian packets, which infinitely extend along the reference lines.The two-parametric superposition of column Gaussian packets is performed along the reference surface, and represents a generalization of the two-parametric superposition of paraxial Gaussian beams.
We then demonstrate that the column Gaussian packets are not singular at caustics only if the reference lines are tangent to the central rays at the reference surface.If we choose the reference lines along the central rays, the two-parametric superposition of column Gaussian packets reduces to the two-parametric superposition of paraxial Gaussian beams along the reference surface.
There are two different high-frequency asymptotic ray theories with frequency-independent amplitudes: the isotropic ray theory based on the assumption of equal velocities of both elastic S waves, and the anisotropic ray theory assuming both waves are strictly decoupled.
Here the term 'different' means that the isotropic ray theory is not a special case of the anisotropic ray theory for decreasing anisotropy, and that both theories yield different S waves in equal media.In the isotropic ray theory, the polarization vectors of waves do not rotate about the third eigenvector of the Christoffel matrix, whereas in the anisotropic ray theory they coincide with the first two eigenvectors of the Christoffel matrix which may rotate rapidly.In 'weakly anisotropic' media, at moderate frequencies, the actual wave polarization tends to remain unrotated about the third eigenvector of the Christoffel matrix, but is partly attracted by the rotation of the first two eigenvectors of the Christoffel matrix.The intensity of the attraction increases with frequency.This behaviour of the actual wave polarization is described by the coupling ray theory proposed, for example, by Kravtsov (1968) or Naida (1977Naida ( , 1979) ) for electromagnetic waves, and by Coates & Chapman (1990) for elastic S waves.The frequency-dependent coupling ray theory is the generalization of both the zero-order isotropic and anisotropic ray theories and provides continuous transition between them.The coupling ray theory is applicable to coupled waves at all degrees of anisotropy, from isotropic to considerably anisotropic media.The numerical algorithm for calculating the frequency-dependent coupling-ray-theory tensor Green function has been designed by Bulant & Klimeš (2002).The coupling-ray-theory tensor Green function is frequency-dependent, and is usually calculated for many frequencies.This frequency dependence can be overcome by the prevailing-frequency approximation of the coupling ray theory with frequency-independent amplitudes according to Klimeš & Bulant (2012, 2016).
The presented integral superpositions of column Gaussian packets and paraxial Gaussian beams may correspond to the anisotropic ray theory, to the frequency-dependent coupling ray theory for S waves or to the prevailing-frequency approximation of the coupling ray theory in anisotropic media, or to the isotropic ray theory in isotropic media.The equations can be used in both Cartesian and curvilinear coordinates.
The lower-case indices i, j, ... take the values of 1, 2, 3.The upper-case indices I, J, ... take the values of 1, 2. The Einstein summation convention over repeated indices is used.

I N T E G R A L S U P E R P O S I T I O N O F G AU S S I A N PA C K E T S
We consider Cartesian or curvilinear coordinates x i in 3-D space.We consider an orthonomic system x i = xi (γ a ) of real-valued central rays parametrized by real-valued ray coordinates γ a , where γ 1 and γ 2 are the ray parameters, and γ 3 is the parameter along central rays determined by the form of the Hamiltonian function.In the case of the coupling ray theory, the orthonomic system of central rays is represented by the system of reference rays along which the coupling equations are solved.The reference rays are calculated using the reference Hamiltonian function, and correspond to the reference traveltime and the reference slowness vector.In the case of attenuation, we choose the reference Hamiltonian function for tracing the real-valued central rays according to Klimeš & Klimeš (2011, eq. 7).The time-harmonic integral superposition of Gaussian packets, Gaussian beams or column Gaussian packets will be considered with respect to the central rays.
Along the central rays, we calculate the ray-theory traveltime and the corresponding ray-theory amplitude using an arbitrary ray method.The ray-theory traveltime may differ from the reference traveltime, but will be assumed to represent just a perturbation of the reference traveltime.In this way, we shall substitute the first-order and second-order spatial derivatives of the reference traveltime for the corresponding spatial derivatives of the ray-theory traveltime.
The integral superposition of paraxial Gaussian packets in heterogeneous isotropic media was derived using the Maslov asymptotic theory applied to a general 3-D subspace of the 6-D complex phase space by Klimeš (1984b, eq. 51).In a heterogeneous anisotropic medium, the time-harmonic integral superposition of paraxial Gaussian packets reads (eq.11; Klimeš 2014Klimeš , 2018) where function √ det M ab is the product of the square roots of the eigenvalues of matrix M ab .The individual square roots are taken with positive real parts.Here ω is the circular frequency, τ ( xn ) is the traveltime at point xi , p k ( xn ) is the slowness vector corresponding to the orthonomic system of central rays at point xi , f kl ( xn ) is the complex-valued matrix with positive imaginary part describing the shape of the paraxial Gaussian packet centred at point xi , N kl ( xm ) is the matrix of the second-order partial derivatives of the reference traveltime at point xi , and A i[ j] ( xm , ω) is the complex-valued vectorial or tensorial ray-theory amplitude.
The amplitudes of Gaussian packets corresponding to a particular source are vectorial without optional subscript [j] .The amplitudes of Gaussian packets designed to compose the Green tensor (Klimeš 2012) are tensorial with optional subscript [j] .
In anisotropic media, vectorial or tensorial amplitude A i[ j] ( xm , ω) may represent either the anisotropic ray-theory vectorial or tensorial amplitude, or the frequency-dependent coupling ray-theory vectorial or tensorial amplitude calculated using the scalar reference amplitude Downloaded from https://academic.oup.com/gji/article-abstract/215/3/1739/5071951 by guest on 11 March 2019 corresponding to the orthonomic system of central rays, or the vectorial or tensorial amplitude of the prevailing-frequency approximation of the coupling ray theory (Klimeš & Bulant 2012, 2016).In isotropic media, vectorial or tensorial amplitude A i[ j] ( xm , ω) represents the isotropic-ray-theory vectorial or tensorial amplitude.
Matrix f kl ( xn ) may be chosen arbitrarily, but must smoothly vary with coordinates xn .Matrix f kl ( xn ) need not satisfy the equations for Gaussian packets propagating along the central rays, because Gaussian packets arriving at different points of the same central ray may correspond to different initial conditions for the shape of Gaussian packets.Matrix f kl ( xn ) may also depend on frequency ω.

G E N E R A L I N T E G R A L S U P E R P O S I T I O N O F C O L U M N G AU S S I A N PA C K E T S
We consider the points of intersection of central rays determined by ray parameters γ 1 and γ 2 with the reference surface.In the integral superposition, we introduce curvilinear coordinates with ξ 3 = 0 along the reference surface which becomes the coordinate surface.We choose ξ I = γ I along the reference surface.The ξ 3 coordinate lines corresponding to constant ξ 1 and ξ 2 will represent the reference lines of column Gaussian packets which will infinitely extend along these reference lines.The reference lines for one-parametric analytical asymptotic integration in superposition (1) are thus represented by the ξ 3 coordinate lines.
We change the coordinates in integral superposition (1), where For fixed ξ I , we consider the dependence of quantities in the integral superposition on ξ 3 .The quadratic Taylor expansion of coordinates with respect to ξ 3 reads xi (ξ 3 ) xi (0) where and Note that contravariant vector Z i 3 is tangent to the reference lines.The quadratic Taylor expansion of the reference traveltime with respect to ξ 3 reads Since we assume that traveltime τ ( xn ) in superposition ( 1) is a perturbation of the reference traveltime, we apply expansion (8) approximately also to the traveltime in superposition (1).The perturbation is contained in term τ (0).The linear Taylor expansion of the reference slowness vector with respect to ξ 3 reads The mixed quadratic Taylor expansion of with respect to ξ 3 and x i − xi (0) then reads where We perform the summations in eq. ( 11) and obtain the following expansion: Downloaded from https://academic.oup.com/gji/article-abstract/215/3/1739/5071951 by guest on 11 March 2019 which does not contain Z i 33 (0) and is thus independent of the curvature of the reference lines x i = xi (ξ 3 ).Expansion (13) may be expressed as We now integrate superposition (3) with respect to ξ 3 and obtain the two-parametric integral superposition where ξ 1 = γ 1 and ξ 2 = γ 2 along the reference surface of integration.Complex-valued traveltime can be expressed as with where all quantities are taken at the reference surface.Note that which means that the column Gaussian packets extend infinitely along the reference lines corresponding to constant ξ 1 and ξ 2 .We define the derivatives of coordinates x i along the reference surface with respect to ray parameters ξ M = γ M .Contravariant vectors Z i 1 and Z i 2 are tangent to the reference surface.
We introduce 2 × 2 matrix We define transformation matrix which represents the inverse matrix to matrix Z i m given by definitions ( 5) and (20), Covariant vectors Z 1i and Z 2i are perpendicular to the reference ξ 3 lines.Covariant vector Z 3i is perpendicular to the reference surface.We may consider decomposition f i j = ( f i j − N i j ) + N i j , apply identity (19) and definition (21) to term ( f i j − N i j ), and express matrix (18) as Considering identity following from definition (23), we express matrix (24) as which reads where all quantities are taken at the reference surface.
We covariantly transform matrices N ij and f ij from general coordinates x i to reference coordinates ξ m using the transformation matrix given by definitions ( 5) and ( 20), Two-parametric integral superposition (15) then reads and matrix ( 21) reads We calculate the determinant of the 2 × 2 matrix F M N − M M N using expression (31) and obtain det ( We insert relation (32) into integral superposition (30) and obtain the two-parametric integral superposition of column Gaussian packets which infinitely extend along the reference lines corresponding to constant ξ 1 and ξ 2 .Complex-valued paraxial traveltime θ is given by expansion ( 17) with the second-order derivatives (27).

I N T E G R A L S U P E R P O S I T I O N S O F G AU S S I A N B E A M S
Since matrix f ik is always finite only if N ik Z k 3 is always finite.To keep N ik Z k 3 finite at caustics, we must choose contravariant vector Z k 3 tangent to the central ray.Curvilinear coordinate ξ 3 then represents a local parameter along the central ray.Covariant vectors Z 1i and Z 2i are then perpendicular to the central rays, and covariant vector Z 3i is perpendicular to the reference surface.
Hamilton's equations of rays (ray tracing equations) yield For the reference lines coinciding with central rays, we may thus express Considering identity we express matrix (27) for the superposition of Gaussian beams as For the special case of the reference surface coinciding with a wavefront, relation ( 38) is analogous to the transformation of the second-order derivatives of the traveltime from ray-centred coordinates to general coordinates ( Červený & Klimeš 2010, eq. 36).
For matrix (38) of the second-order traveltime derivatives, integral superposition (33) with complex-valued traveltime (17) represents the integral superposition of Gaussian beams.

G AU S S I A N B E A M S I N R AY-C E N T R E D C O O R D I N AT E S
Since the integral superposition of Gaussian beams is frequently expressed in terms of the quantities specified in ray-centred coordinates, we shall now transform integral superposition (33) with complex-valued traveltime (17) corresponding to the second-order traveltime derivatives (38) from general coordinates to ray-centred coordinates.
We define transformation matrix from ray-centred coordinates q a to Cartesian or curvilinear coordinates x i , and transformation matrix from Cartesian or curvilinear coordinates x i to ray-centred coordinates q a .We define covariant derivatives of the traveltime in ray-centred coordinates.Derivatives in the wavefront tangent plane coincide with the partial derivatives of the traveltime in ray-centred coordinates ( Červený & Klimeš 2010, eq. 38).
We insert traveltime derivatives (38) into transformation relation ( 42) and obtain the 2 × 2 complex-valued matrix of the second-order derivatives of the complex-valued traveltime of a Gaussian beam in ray-centred coordinates.
We finally express traveltime derivatives (38) in terms of matrix (43) using relation by Červený & Klimeš (2010, eq. 36), We now inspect expression (44) by comparing it with expression (38).We calculate the third column of covariant derivatives (41) from expression (44), We now compare this expression with the third column of covariant derivatives (41) calculated directly from expression (38), We insert relation into expression (47) and obtain the following expression: identical with expression (46).Using relation (38), we express matrix F M N appearing in superposition (33) in terms of matrix f i j as Relations ( 28) and ( 50) yield Analogously to transformation relation (44), we express traveltime derivatives N ij in terms of the 2 × 2 matrix M (q) AB of the second-order derivatives of the ray-theory traveltime in ray-centred coordinates ( Červený & Klimeš 2010, eq. 36) We insert transformations (44) and (52) into relation ( 51) and obtain the following relation: is identical to the matrix Q K L of geometrical spreading in ray-centred coordinates.

C O N C L U S I O N S
Two-parametric integral superposition (33) with complex-valued paraxial traveltime (17) corresponding to the second-order derivatives ( 27) represents the superposition of column Gaussian packets which infinitely extend along the reference lines corresponding to constant ξ 1 and ξ 2 .The integral superposition of Gaussian beams is a special case of the integral superposition of column Gaussian packets.
The column Gaussian packets in the integral superposition are always finite only if the reference lines coincide with the central rays, that is, if the column Gaussian packets coincide with Gaussian beams.To avoid problems with caustics, we thus have to consider the integral superposition of Gaussian beams instead of the integral superposition of general column Gaussian packets.
In general Cartesian or curvilinear coordinates, the integral superposition of Gaussian beams is given by the two-parametric integral superposition (33) with complex-valued paraxial traveltime (17) corresponding to the second-order derivatives (38).If we consider the quantities specified in ray-centred coordinates, the integral superposition of Gaussian beams is given by the two-parametric integral superposition (56) with complex-valued paraxial traveltime (17) corresponding to the second-order derivatives (44).

A C K N O W L E D G E M E N T S
I am indebted to Vlastislav Červený for formulating the problem and for many invaluable discussions on the topic of this paper.The suggestions by Einar Iversen and an anonymous reviewer made it possible for me to improve the paper considerably.
The research has been supported by the Grant Agency of the Czech Republic under contracts 16-01312S and 16-05237S, and by the members of the consortium 'Seismic Waves in Complex 3-D Structures' (see 'http://sw3d.cz').