{'fraction_non_alphanumeric': 0.056076467910787436, 'fraction_numerical': 0.02266727355484752, 'mean_word_length': 4.232772308669419, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Frozen Gaussian Approximation for 3-D Elastic Wave Equation and Seismic Tomography

The purpose of this work is to generalize the frozen Gaussian approximation (FGA) theory to solve the 3-D elastic wave equation and use it as the forward modeling tool for seismic tomography with high-frequency data. FGA has been previously developed and verified as an efficient solver for high-frequency acoustic wave propagation (P-wave). The main contribution of this paper consists of three aspects: 1. We derive the FGA formulation for the 3-D elastic wave equation. Rather than standard ray-based methods (e.g. geometric optics and Gaussian beam method), the derivation requires to do asymptotic expansion in the week sense (integral form) so that one is able to perform integration by parts. Compared to the FGA theory for acoustic wave equation, the calculations in the derivation are much more technically involved due to the existence of both P- and S-waves, and the coupling of the polarized directions for SH- and SV-waves. In particular, we obtain the diabatic coupling terms for SH- and SV-waves, with the form closely connecting to the concept of Berry phase that is intensively studied in quantum mechanics and topology (Chern number). The accuracy and parallelizability of the FGA algorithm is illustrated by comparing to the spectral element method for 3-D elastic wave equation in homogeneous media; 2. We derive the interface conditions of FGA for 3-D elastic wave equation based on an Eulerian formulation and the Snell's law. We verify these conditions by simulating high-frequency elastic wave propagation in a 1-D layered Earth model. In this example, we also show that it is natural to apply the FGA algorithm to geometries with non-Cartesian coordinates; 3. We apply the developed FGA algorithm for 3-D seismic { wave-equation-based traveltime tomography and full waveform inversion, respectively


INTRODUCTION
Images computed by seismic tomography can provide crucial information for the subsurface structures of Earth at different scales, and the understanding of tectonics, volcanism, and geodynamics (e.g. Aki & Lee 1976;Romanowicz 1991;Rawlinson et al. 2010;Zhao 2012). Wave-equation-based seismic tomography solves the nonlinear optimization problem iteratively for velocity models by computing seismograms and sensitivity kernels in 3-D complex models Liu & Gu 2012;Tong et al. 2014). This leads to successful applications including imaging the velocity models of the southern California crust (Tape et al. 2009(Tape et al. , 2010, the European upper mantle (Zhu et al. 2012), the North Atlantic region (Rickers et al. 2013), and the Japan islands (Simute et al. 2016). The performance of seismic tomography is restricted by how accurate one can solve the 3-D wave equation for synthetic seismograms and sensitivity kernels. The dominant frequency of a lowfrequency earthquake is around 1 Hz, while the one of a typical earthquake is around 5 Hz (e.g. Nakamichi et al. 2003), which leads to demanding and even unaffordable computational cost. Real applications usually chose significantly low-frequency data in order to be computationally feasible (e.g. Tape et al. 2009;Zhu et al. 2012;Simute et al. 2016), since simulating high-frequency seismic waves requires much more powerful computational resources on both memories and CPU time than the low-frequency waves. This imposes a demand on improving the efficiency and accuracy of numerical methods for computing high-frequency waves in order to make use of real seismic data around their dominant frequencies.
In previous works (Chai et al. 2017(Chai et al. , 2018, the authors have developed and verified frozen Gaussian approximation (FGA) as an efficient solver for computing high-frequency acoustic wave propagation (P-wave), wave-equation-based traveltime tomography and full waveform inversion (FWI). The core idea of FGA is to approximate seismic wavefields by fixed-width Gaussian wave packets, whose dynamics follow ray paths with the prefactor amplitude equation derived delicately from an asymptotic expansion on phase plane. With a multicore processors computer station, the property that the FGA algorithm is embarrassingly parallel makes possible the application of FGA to compute 3-D high frequency sensitivity kernels, and further used for 3-D traveltime tomography and FWI. Compared to other ray-based methods including WKBJ (e.g. Chapman 1976;Chapman & Drummond 1982), WKM (e.g. Ni et al. 2000;Helmberger & Ni 2005), generalized ray theory (e.g. Helmberger 1968;Vidale & Helmberger 1988), seismic traveltime tomography (e.g. Aki et al. 1977;Tong et al. 2017), Kirchhoff migration (e.g. Gray 1986;Keho & Beydoun 1988) and Gaussian beam migration (e.g. Hill 1990Hill , 2001Nowack et al. 2003;Gray 2005;Gray & Bleistein 2009;Popov et al. 2010), FGA does not need to solve ray paths by shooting to reach the receivers, and provides accurate solutions at the presence of caustics and multipathing, with no requirement on tuning beam width parameters to achieve a good resolution (Cerveny et al. 1982;Hill 1990;Fomel & Tanushev 2009;Qian & Ying 2010;Lu & Yang 2011).
In this paper, we first generalize the FGA theory to solve the 3-D elastic wave propagation. Different from the WKBJ theory or Gaussian beam method which can be derived by direct asymptotic expansion, the derivation of FGA formulation requires to do the asymptotic expansion in an integral form (weak sense) so that one is able to perform integration by parts to eliminate the extra constraints yielded by direct asymptotic expansion. Compared to the previous works on FGA (Lu & Yang 2012a;Chai et al. 2017Chai et al. , 2018, the calculations for the derivation are much more technically involved due to the existence of both P-and S-waves. Particularly, the diabatic coupling of the polarized directions for SH-and SV-waves leads to a term closely connected to the concept of Berry phase which is intensively studied in quantum mechanics and topology (Chern number) (e.g. Berry 1984;Simon 1983). We prove the accuracy and parallelizability of the FGA algorithm by comparing to the spectral element method (e.g. Komatitsch & Tromp 1999;Komatitsch et al. 2005;Tromp et al. 2008) for the 3-D elastic wave equation in homogeneous media, where one has the analytical solution as a benchmark. We derive the interface conditions of FGA for the 3-D elastic wave equation based on an Eulerian formulation and the Snell's law. These interface conditions are verified by simulating high-frequency elastic wave propagation in a 1-D layered Earth model. We also show how natural to apply the FGA algorithm to geometries with non-Cartesian coordinates in the simulation of this model.
The second part of this paper contributes to the study of seismic tomography with high-frequency data. Wave-equation-based traveltime tomography and FWI are able to generate higher-resolution images of the Earth's interior than the conventional ray-based tomography method (Virieux & Operto 2009). Especially, some recent novel studies on FWI using optimal transport distance (e.g. Wasserstein metric) are able to capture traveltime differences between seismic signals, and thus overcomes the cycling effect in standard FWI method (e.g. Engquist & Froese 2014;Métivier et al. 2016;Yang et al. 2018;Chen et al. 2017). However, seismic tomography using high-frequency data increases computational cost drastically, which restricts the application of these wave-equation-based inversion methods, and only low-frequency data are modeled and inverted in some real applications (e.g. Tape et al. 2007;Zhu et al. 2012;Simute et al. 2016). To address the computation challenge, we use the developed FGA algorithm to compute the 3-D high-frequency sensitivity kernels for wave-equationbased traveltime tomography and FWI, respectively. We apply FGA to both wave-equation-based traveltime tomography and FWI on synthetic crosswell seismic data with dominant frequencies as high as those of real crosswell data, and use a hierarchical approach which first uses traveltime tomography to create a macro-scale model and then adopts FWI to generate a high-resolution micro-scale model.

FROZEN GAUSSIAN APPROXIMATION
We derive the frozen Gaussian approximation (FGA) formulation for the following inhomogeneous, isotropic elastic wave equation: where u(t, x) is the wavefield, ρ(x) is the density, λ(x) and µ(x) are the Lamé parameters, and F(t, x) is the body force. We shall assume F = 0 in this section as a sake of simplicity for presenting the FGA formula. Remark that eq. (1) can be rewritten as, with the P-wave, and S-wave speeds defined as, We consider the elastic wave eq. (1), with the following initial conditions Note that all the bolded variables and parameters will be referred to vectors in R 3 without any further clarification.

Solution Ansatz
The FGA approximates the wavefield u(t, x) in eq. (1) by a summation of dynamic frozen Gaussian wave packets, with the weight functions α k p,s (y, q, p) = 1 2kc p,s |p| 3 kf k (y)c p,s |p| ± ig k (y) ·n p,s .
In eq. (5), i = √ −1 is the imaginary unit, and we use unbolded subscripts/superscripts "p" and "s" to indicate P-and S-waves respectively, and superscripts k to indicate quantities depending on wave number k. The quantities that have both p and s as subscripts/superscripts mean that they can be referred to both P-and S-waves. Here G p,s ± refers to the initial sets of Gaussian center q and propagation vector p for P-and S-waves respectively, and ± indicates the two-way wave propagation directions correspondingly.N p,s (t) are unit vectors indicating the polarized directions of P-and Swaves, i.e. if (q, p) ∈ G p ± , thenN p P p ; and if (q, p) ∈ G s ± , thenN s ⊥ P s . In (7),n p,s are the initial directions of P-and S-waves, i.e.n p,s =N p,s (0), and the "±" on the right-hand-side of the equation indicate that the α k p,s correspond to (q, p) ∈ G p,s ± . Remark that e ikp·(y−q)− k 2 |y−q| 2 can be understood as complex localized wave-packet centered at q with propagation vector p, and ψ k p,s is the projection of the initial wavefield onto each wave-packet computed by eq. (6) with y serving as the dummy variable in the integration. Associated with each frozen Gaussian wave-packet, the timedependent quantities are the position center Q p,s (t, q, p), momentum center P p,s (t, q, p), amplitude a p,s (t, q, p) and unit direction vectorsN p,s (t). Note that all the S-waves discussed in the formulation will include both SH-and SV-waves.

Formulation and Algorithm
The derivation of the FGA formulation involves with the asymptotic expansion in the integral form (weak sense), with proper integration by parts performed to convert powers of distance to the Gaussian center x − Q to orders of the wavelength k −1 . It is quite lengthy and technically involved, and thus we leave it to Appendix A for the readers who are interested in the mathematical details, and only present the FGA formulation as below.
For the "+" wave propagation, i.e., (q, p) ∈ G p,s + , the Gaussian center Q p,s (t, q, p) and propagation vector P p,s (t, q, p) follow the ray dynamics with initial conditions Q p,s (0, q, p) = q and P p,s (0, q, p) = p.
For the "-" wave propagation, i.e., (q, p) ∈ G p,s − , the Gaussian center Q p,s (t, q, p) and propagation vector P p,s (t, q, p) follow the ray dynamics with initial conditions Q p,s (0, q, p) = q and P p,s (0, q, p) = p.
Remark that, the equations for "-" wave propagation (10) have opposite signs of the right-hand side to the equations for "+" wave propagation (8), and other than that, they are the same. Actually they can be both viewed as the Hamiltonian system with H(Q, P ) = ±c p,s (Q)|P | as the Hamiltonian function for "+" and "-" wave propagation, respectively.
The prefactor amplitudes a p,s (t, q, p) satisfy the following equations, where S-waves have been decomposed into SH-and SV-waves, with the initial conditions a p,sv,sh = 2 3/2 , andN sv andN sh are the two unit directions perpendicular to P s , referring to the polarized directions of SV-and SH-waves, respectively. Here "±" corresponds to the two-way wave propagation directions, and we have used the short-hand notations, Note that, the prefactor equation (12) is consistent with the one for acoustic wave equation with c 2 = (λ + 2µ)/ρ (Chai et al. 2017) , and the last terms on the right-hand-side of (13)-(14) indicate the diabatic coupling of the polarized directions for SH-and SV-waves, which are closely connected to i.e., to choose the proper sets G p,s ± by the fast FBI transform algorithm (Yang et al. 2013, Page 2) to decompose the initial wavefields into Gaussian wave packets, and calculate the corresponding weight function ψ k p,s defined in eq. (6); Step 2, Time propagation, i.e., to solve by some numerical solver eqs.  Fig. 1 to describe the FGA algorithm, with the technique details discussed in the figure caption.

Accuracy and Parallelizability
We check the accuracy and parallelizability of the FGA method by simulating the elastic wave propagation at a homogeneous media with absorbing boundary conditions, for which one can solve the analytical solution as a benchmark. Specifically, we consider T as the source time function at x 0 , and ρ, λ and µ as constants. Eq. (16) has the following analytical solution where we have used the Einstein's index summation convention, and δ ij is the Kronecker delta function. Here we use x = (x 1 , x 2 , x 3 ) T for a simple formulation of the analytical solution, which is actually x = (x, y, z) T in standard notations.
In the numerical tests, we simulate the elastic wave propagation in a domain of the size 128 km × 128 km × 128 km, and take the source time function as with f as the frequency, T 0 = 0.1768 s and σ = 0.8660 s. We choose P-and S-wave speeds as shows that the speed-up ratio for FGA is almost 2 when one doubles the number of processors, which indicates that the FGA algorithm is embarrassingly parallel. Remark that, in the simulation for the https://geodynamics.org/cig/software/specfem3d/  8) and (12) where Q and P are 3-D real vectors and the prefactor amplitude a is a complex number; while in SPECFEM3D, we use 128 elements in each direction with 5 3 nodes in each element. One needs to roughly to compute a total number of 3 × 128 3 × 5 3 ≈ 800 millions of variables in the simulation, where the prefactor 3 is due to u is a 3-D real vector in eq. (1). In addition, the stability conditions on time step for solving the ODE systems (8) and (12) by RK4 is better than solving the elastic wave equation (1) by SPECFEM, since the CFL condition is restricted by small wavelength in solving eq. (1).

INTERFACE CONDITIONS AND SIMULATION OF A 1-D LAYERED EARTH MODEL
In this section, we derive the interface conditions of FGA for 3-D elastic wave propagation in heterogeneous media with strong discontinuities, e.g., the Moho surface and Core-mantle boundary. We verify these interface conditions by simulating a 1-D layered Earth model, which actually requires to revise the FGA solution in spherical coordinates (non-Cartesian coordinates). This can be easily done by an observation that the only step of the FGA algorithm depending on the x-coordinate is the reconstruction part given by eq. (5) and all the other steps including the initial wavefield decomposition by eqs. (6) and (7) and propagation of ray equations eqs. (8), (10), (12), (13) and (14) are x-coordinatefree. Therefore, to generalize the FGA method in spherical coordinates, one simply needs to change x = (x, y, z) T in eq. (5) to x = (r cos θ cos φ, r sin θ cos φ, r sin φ) with r as the distance to earth center, θ and φ as the longitude and latitude degrees.

Transmission and Reflection Interface Conditions
For a sake of clarity and simplicity, we only give the interface conditions for a flat interface located at z = z 0 , and for a reflecting geometry of general shape, one needs to apply the formulation in the local  tangent-normal coordinates by treating the tangential direction as the local flat horizontal interface.
The wave speeds of the two layers are assumed to be, In Fig. 7, we only consider an incident Gaussian wave packet for P-wave hitting the interface at z = z 0 , and then reflected and transmitted as Gaussian wave packets for P-and SV-waves, respectively.
The other cases including an incident Gaussian wave packet for SV-and SH-waves can be handled similarly, although there is no interaction between P-and S-waves for the case of SH-wave.
Associated with each Gaussian wave packet, one needs to provide the reflection and transmission conditions for a p,s , Q p,s , and P p,s , which follow the Snell's Law and the Zoeppritz equations (Yilmaz 2001). However, for FGA, what is different from standard ray theory on the flat interface is that, one also needs to derive the interface conditions for Z p,s which will change after the Gaussian wave packet hits the interface and affect the dynamics of a p,s given by eqs. (12)-(14). Eq. (15) implies that to derive  1.9393, which is slightly smaller than 2 indicating an almost perfectly parallel efficiency. On the other hand, as a comparison, SPECFEM3D is used with 128 elements in each spatial direction to achieve a comparable accuracy to FGA. The speed-up ratio for the SPECFEM3D solver is around 1.4828, which is smaller that those of FGA.
This is because SPECFEM3D solves eq. (1) on a parallel computer with N processors by partitioning the whole domain into N slabs with each processor solving the equation in each slab. Therefore, for each time step, each processor needs to communicate with its neighbors to get necessary boundary information, which decreases the speed-up ratio.  Figure 7. Cartoon illustration of an incident Gaussian wave packet for P-wave hitting the interface at z = z 0 , and then reflected and transmitted as Gaussian wave packets for P-and SV-waves. Here the G in,re,tr p,s stands for the Gaussian wave packet for the incident, reflected and transmitted P-and SV-waves, respectively. We denote θ i , θ r , θ t to be the incident, reflection and transmission angles of P-waves, and φ r , φ t to be the reflection and transmission angles of SV-waves, respectively. the interface conditions of Z p,s will be equivalent to derive the interface conditions for ∂ z Q p,s and ∂ z P p,s , which requires to use the conservation of level set functions designed in the Eulerian frozen Gaussian approximation formula (Lu & Yang 2012b;Wei & Yang 2012). The mathematical details of the derivation are lengthy and technical, and thus we leave them to Appendix B for the interested readers, and only present the results here: where (Q in,re,tr , P in,re,tr ) corresponds to the center and propagation vector of incident, reflected and transmitted Gaussian wave packet for either P-or S-waves, respectively, and P in,re,tr = (p x , p y , p in,re,tr z ).
Here e 3 = (0, 0, 1) is a row vector, F and W are two 3 × 3 matrices, F T = W −1 , and Let us explain the formulation of reflection and transmission interface conditions with more detials for the incident P-wave as illustrated in Fig. 7 to be the center, propagation vector and prefactor amplitude of the corresponding incident, reflected and transmitted Gaussian wave packets for P-and S-waves respectively, and P in p = (p x , p y , p z ), then for P re,tr p = (p x , p y , p re,tr z,p ) and P re,tr s = (p x , p y , p re,tr z,s ), p tr z,s = (p z n 3 ) 2 + (n 2 3 − 1)(p 2 x + p 2 y ), Moreover, if one denotes θ i , θ r , θ t to be the P-wave incident, reflection and transmission angles, and φ r , φ t to be the SV-wave reflection and transmission angles, respectively, then the Zoeppritz equations (Yilmaz 2001 with the matrix M as where ρ 1,2 are the densities for the layers 1 and 2, respectively. More explicitly in eq. (19) corresponding to the case in Fig. 7

Waveguide example in a 1-D layered Earth model
We verify the interface conditions (19) and (21) by simulating a waveguide example in a 1-D layered Earth model, with the layered P-wave velocity given in Fig. 8(a) following the data in the IASP91 model (Kennett 1991;IRIS DMC 2010). We are particularly interested in choosing the 410-km discontinuity for the numerical proof of the conditions (19) and (21), which presents to a 5 − 6% increase on P-wave velocity and calibrates the mantle transition zone. We consider a radially symmetric surface source as shown in Fig. 8(b), so that the elastic wave equation (1)  i.e., We choose the initial condition for the elastic wave equation (1) as where r = |x|, r 0 = 600 km, and φ 0 (r) = exp − r 2 2σ 2 cos 2πr , with σ = 3.3146 km and = 7.3631 km. The dominant frequency is around 1.36 Hz. We solve eq. (23) using 1-D finite difference method with fine enough grid points as the reference solution, to which we compare the full elastic wave solution of (1) computed by the FGA algorithm. Fig. 9 shows the seismic signals of P-waves received at stations of depth 480 km, 420 km and 360 km, respectively, where one can see a good agreement of FGA simulation with the reference solution for the P-waves.

WAVE-EQUATION-BASED TRAVELTIME TOMOGRAPHY AND FWI
In this section, we apply the developed FGA algorithm for 3-D wave-equation second layer and homogeneity in the y-direction, We start with the following initial velocity model then use FGA to compute the forward and adjoint wavefields and construct 3-D kernels of different phases by the methods of wave-equation-based traveltime tomography Liu & Gu 2012) and FWI (Pratt & Shipp 1999;Virieux & Operto 2009). Since FWI requires a more sophisticated initial velocity model for the convergence than wave-equation-based traveltime tomography, we use a hierarchical strategy as in our previous work (Chai et al. 2018), which first uses waveequation-based traveltime tomography to create a macro-scale model and then adopts FWI to generate a high-resolution micro-scale model. For the signals received at the top station (Fig. 11), we show the corresponding kernels of different phases in Fig. 12. The inversion results are shown in Fig. 13. The damping and smoothing parameters are chosen empirically by forcing the variation in each iteration less than 3%. Note that artifacts are visible in the final images, which are unavoidable and mainly caused by the uneven data coverages.

DISCUSSION AND CONCLUSIONS
This is the last part of our theoretical study on the application of the frozen Gaussian approximation (FGA) method to seismic tomography using high-frequency seismic data. Together with Chai et al. (2017Chai et al. ( , 2018, we systematically derive the FGA formulation for acoustic and elastic wave equations, derive the transmission and reflection conditions of FGA for sharp interfaces (e.g., Moho surface and Core-mantle boundary), reformulate the equations of the forward and adjoint wavefields for a convenient application of the FGA algorithm, propose a fast Gaussian summation algorithm for the benchmark solution is available. We also test the performance of FGA in a 1-D layered Earth model for checking the accuracy of the interface conditions, and apply FGA to the cross-well example for the 3-D seismic tomography using high-frequency data. We remark that, in the original mathematical work on FGA for strictly linear hyperbolic systems (Lu & Yang 2012a), the FGA formulation is rather complicated, and done for a purpose of rigorously proving its accuracy in the high-frequency regime.
In this paper, we weakly expand the solution ansatz of FGA with a consideration of decomposition into P-and S-waves, and reach simple prefactor equations (12), (13) and (14) after lengthy and tricky simplifications. This also generalizes the results in Lu & Yang (2012a) to some extent in the sense that the elastic wave equation is actually not strictly hyperbolic, which has an eigenfunction space of multiplicity-2 (corresponding to S-waves). By these efforts, we hope to bridge a gap between the semiclassical approximation theory in mathematics and the application of 3-D seismic tomography in geophysics. We also remark that, it will be interesting to study the performance of the FGA method in the optimal transport theory-based seismic tomography, with the normalization strategies proposed in, e.g., Engquist   train neural networks for seismic inversion; 3. to develop the FGA algorithm for the optimal transport theory-based seismic tomography using high-frequency data.

ACKNOWLEDGMENTS
We are grateful to the support from the Center for Scientific Computing from the CNSI, MRL: an NSF MRSEC (DMR-1121053) and NSF CNS-0960316. JH, LC and XY were partially supported by the NSF grants DMS-1818592, DMS-1418936 and DMS-1107291. JH was also partially supported by the graduate summer fellowship of Department of Mathematics at University of California, Santa Barbara. PT was partially supported by the MOE AcRF Tier-1 grant (M4011899.110). XY also thanks Professor Qinya Liu for useful discussions on the WKM and GRT methods. We also thank the anonymous referee for bringing us the interesting questions on the optimal transport theory-based seismic tomography.

APPENDIX A: DERIVATION OF THE FGA FORMULATION
This appendix is devoted to the mathematical details of deriving the FGA formulation. It is lengthy and technically involved, making use of a few mathematical concepts including weak asymptotic expansion, Schwartz class functions, canonical transform and symplecticity. We shall do our best to make the details as straightforward as possible so that they are accessible for interested readers with different background.
For a sake of clarity and simplicity, we only consider the case of λ and µ being constants in eq. (1) with only "+" wave propagation direction, and the derivation of FGA for general λ(x) and µ(x) with two-way wave propagation directions will be essentially the same but just with more calculations. We plug eq. (6) into eq. (5), combine the phase functions, and define the total phase Φ p,s as We call two functions f (y, q, p) and g(y, q, p) are weakly equal to each other, denoted by f ∼ g, if f (y, q, p) exp(ikΦ p,s ) dy dq dp = g(y, q, p) exp(ikΦ p,s ) dy dq dp.
Note that both f and g can be scalar-, vector-, matrix-or tensor-valued functions.
A map: (q, p) → Q p,s (q, p), P p,s (q, p) is called canonical transformation if the associated Jacobian matrix is symplectic, i.e., for any (q, p), where Id 3 is a 3×3 identity matrix. It is easy to verify that the Hamiltonian flow given by either eq. (8) or (10) is a canonical transform by observing that with the Hamiltonian function H p,s (Q, P ) = c p,s (Q p,s )|P p,s |, thus J p,s (q, p) satisfies eq. (A.3). This property guarantees that the matrix Z p,s defined in eq. (15) is always invertible for all t > 0 by the following argument, with the subscripts p, s omitted for convenience, where Z * is the conjugate transpose of Z, and in the last equality, we have used eq. (A.3). Therefore, the Z −1 p,s used in eqs. (12), (13) and (14) is a well behaved quantity. With the theoretical preparation above, we are ready to derive the formulation of FGA. We only consider the case F = 0 in eq. (1), and refer to Chai et al. (2018, Section 3.1) for the treatment of a general source time function. We start with the following form of Gaussian wave packet as the solution ansatz, u p,s (t, x, y, q, p) = A p,s (t, q, p) exp ikΦ p,s (t, x, y, q, p) , (A.4) with Φ given by eq. (A.1). Remark that, for the readers who are familiar with Gaussian beam method (e.g. Hill 1990Hill , 2001Nowack et al. 2003;Gray 2005;Gray & Bleistein 2009;Popov et al. 2010), it is easy to see that, a direct plugin of eq. (A.4) into eq. (1) will yield more equations than the number of variables due to the missing of a time-dependent Hessian function in the total phase Φ. Therefore, one can only derive the FGA formulation using the weak asymptotic expansion with a sense of equality defined in eq. (A.2), where the following integration by parts lemma will help to eliminate the extra constraints by converting the powers of x − Q to the powers of k −1 .
Lemma of integration by parts: For any vector a(y, q, p) = (a j ) and matrix M (y, q, p) = (M ij ) in Schwartz class, i.e., with decay properties at infinity so that the integrals in eq. (A.2) are well defined, one has the following integration by parts formula in the componentwise form, with ∂ z = where the Einstein's index summation convention is used, the matrix Z is given in eq. (15), and we refer to Lu & Yang (2011, 2012a for the detailed proofs. Note that eq. (1) with F = 0 can be rewritten as the following "curl" form, Notice that eq. (A.6) is linear, and thus one can derive the prefactor equations for P-and S-waves individually by assuming A p P or A s ⊥ P in eq. (A.4), respectively. Without loss of generality, we first consider the prefactor equation for the P-wave, with the governing equations for Q p and P p given by eq. (8). For an ease of notations, we shall omit the subscript p in calculating the prefactor equation for the P-wave.
Plugging eq. (A.4) into eq. (A.6) and expanding the asymptotics in the weak sense of (A.2) yield The spatial and temporal derivatives of Φ are given by Notice that the terms containing k(x − Q) will be of O(1) by the lemma of integration by parts, and for P-waves, P × A = 0 and ∇× (x − Q) × A = −2A. Plugging the derivatives of Φ in eq. (A.8) into eq. (A.7) produces, after neglecting the O(1) and lower order terms, where ⊗ means the tensor product, e.g., (A ⊗ P ) jl = A j P l .
Expanding ρ(x) around Q and truncating at order third order, and noticing that ρ(Q)c 2 (Q) = λ + 2µ is constant, one has Taking the second derivative for ∂ QQ ρ and substituting eq. (A.11) bring Plugging eqs, (8) and (A.12) into eq. (A.9), and dividing by ρ yield, in componentwise form, where M jl and N i are given as follows, (A.14) Assuming that A = a pN p , withN p = P |P | , then Plugging this into eq. (A.13) with the ray equations (8), using the P-wave velocity eq. (3) and grouping in powers of (x − Q) produce Denoting ∂ z = (∂ 1 , ∂ 2 , ∂ 3 ) for an ease of notations, applying eq. (A.5) and dropping the lower order terms yield To derive an ordinary differential equation (ODE) instead of a partial differential equation (PDE) for a p , one needs to simplify the terms containing ∂ k a p as Recall that eq. (A.16) holds in the sense of integral form (A.2), and now we shall consider a strong form of eq. (A.16), i.e., equate the integrands of the integrals on both sides. After taking the dot product of integrands with P , one has where the terms (A.17) actually become zero since P · (|P | 2 − P ⊗ P ) = 0, and we have used the Since Z = ∂ z (Q + iP ) by eq. (15), ∂ t Z = ∂ t ∂ z Q + i∂ t ∂ z P . Then eq. (8) implies (A.18) Using eq. (A.18) for further simplifications give where the last two terms can be grouped as which implies a clean ODE for a p as in eq. (12), Similarly, one can derive the prefactor equations for SV-and SH-waves by assuming A = a svN sv + a shN sh withN sv ⊥ P ,N sh ⊥ P andN sv ⊥N sh in eq. (A.4). The calculations will be essentially the same as the prefactor equation for P-waves except that one will have the diabatic coupling terms ofN sv andN sh as shown below, where the interaction terms are given by Also, note that byN sv ⊥N sh , one has that dN sh dt ·N sv + dN sv dt ·N sh = 0.
is symmetric, and then m sh→sv = m sv→sh = 0, which brings eqs. (13) and (14)  The derivation of the interface conditions for ∂ z Q p,s and ∂ z P p,s requires to use the conservation of level set functions designed in the Eulerian frozen Gaussian approximation formula (Lu & Yang 2012b;Wei & Yang 2012), where the idea is to use the following Liouville operator to describe the dyanmics of Gaussian wave packet on phase plane, L p,s = ∂ t + ∂ p H p,s · ∂ q − ∂ q H p,s · ∂ p (B.1) whose corresponding characteristic equations are given by the Hamiltonian systems eqs. (8) and (10) with H p,s = ±c p,s (q)|p|. For example, the prefactor equation of the P-wave for the "+" wave propagation direction is given by, in the Eulerian formulation, We shall derive the interface conditions of ∂ z Q and ∂ z P for the transmitted P-wave and the "+" wave propagation direction, and omit the subscript "p" in the following derivation for a sake of simplicity. Consider a level set function φ(t, q, p) = (φ 1 , φ 2 , φ 3 ) which satisfies Lφ = 0, with φ(0, q, p) = p + iq, (B.2) then the Eulerian formulation of FGA (Lu & Yang 2012b) shows that We will follow the strategy described in Chai et al. (2018); Wei & Yang (2012), and consider the case illustrated in Fig. 7, where the level set functions φ re,tr for the transmitted P-waves satisfy the same evolution as φ in eq. (B.2) with the following interface conditions φ tr (t, q tr , p tr ) = φ tr (t, q in , p in ).

(B.4)
Differentiating eq. (B.4) by the definition of partial derivatives and chain rule, and making use of