Recursive integral time extrapolation of elastic waves using low-rank symbol approximation

Conventional solutions of elastic wave equations rely on inaccurate ﬁnite-difference approximations of the time derivative, which result in strict limit on dispersion and stability conditions. In this work, we derive a general solution to the elastic anisotropic wave equation, in the form of a Fourier integral operator (FIO). The proposed method is a generalization of the previously developed recursive integral time extrapolation (RITE) operators from acoustic to elastic media, and can accurately propagate waves in time using the form of the analytical solution in homogeneous media. The formulation is closely connected to elastic wave-mode decomposition, and can be applied to the most general anisotropic medium. The numerical calculation of the FIO makes use of a low-rank approximation to enable highly accurate and stable wave extrapolations. We present numerical examples including wave propagation in 3-D heterogeneous orthorhombic and triclinic models. and S arrivals are identical, while reﬂected P and S arrivals show obvious difference in both amplitude and phase due to the missing of stiffness gradient terms. To validate our observation, we also carry out wave extrapolation using the pseudospectral method by solving ﬁrst- and


I N T RO D U C T I O N
Elastic wave extrapolation honours elastic effects such as wave-mode conversions, and provides reliable amplitude information which is crucial in seismic imaging of the subsurface. With advances of high-performance computing, seismic processing algorithms such as reverse-time migration (RTM) and full-waveform inversion (FWI) based on elastic kernels are becoming more affordable (Lu et al. 2010;Vigh et al. 2014). In exploration seismology, the most widely used method to solve elastic wave equations involves finite-difference approximations of both spatial and temporal derivatives (Virieux 1984(Virieux , 1986Etgen 1987;Levander 1988;Bernth & Chapman 2010. The pseudospectral method (Kosloff et al. 1984;Reshef et al. 1988;Fornberg 1996) provides accurate calculation of spatial derivatives, but still requires small time steps to avoid temporal dispersion. The numerical accuracy of temporal differentiation can be improved by employing higher order Taylor expansions in the explicit case (Dablain 1986;Crase 1990) and Padé expansions in the implicit case (Chu 2009;Liu & Sen 2009). When implemented using coupled first-order particle velocity-stress systems, these methods usually require to use staggered grids to correctly centre first-order differences of different model parameters (Özdenvar & McMechan 1996;Corrêa et al. 2002;Bernth & Chapman 2010). In global seismology, the finite-element method is a popular choice due to its accurate representation of complex geometry of interfaces (Virieux et al. 2011). In particular, the spectral element method has seen rapid development and wide adoption mainly due to its efficiency (Komatitsch & Vilotte 1998;Komatitsch & Tromp 1999;Chaljub et al. 2007). They are also subject to temporal dispersion and stability constraints when used in explicit time stepping. Recently, several methods have been introduced for stable and dispersion-free time extrapolation of scalar wavefields based on the form of analytical solutions of the acoustic and viscoacoustic wave equations in homogeneous media (Tal-Ezer 1986;Tal-Ezer et al. 1987;Tabei et al. 2002;Etgen & Brandsberg-Dahl 2009;Zhang & Zhang 2009;Pestana & Stoffa 2010;Chu & Stoffa 2010;Fomel et al. 2013;Song et al. 2013;Fang et al. 2014;Sun et al. 2015Sun et al. , 2016. Du et al. (2014) provide a review of existing operators of such nature and refer to them as recursive integral time extrapolation (RITE). Chu & Stoffa (2011) and Firouzi et al. (2012) extend this approach to elastic wave extrapolation in isotropic media.
To mitigate cross-talk between P and S waves, it is often necessary to decouple different wave modes prior to imaging. In isotropic media, wave-mode separation can be achieved using the divergence and curl operators (Aki & Richards 1980). Dellinger & Etgen (1990) implement wave-mode separation in anisotropic media by projecting the vector wavefield onto the polarization directions defined by the Christoffel equation. Yan & Sava (2009 implement wave-mode separation in vertical transversely isotropic (VTI) and tilted transverse isotropic Low-rank elastic wave extrapolation 1479 (TTI) media by introducing space-domain non-stationary filters to handle spatial heterogeneity, and improve the efficiency using the idea of phase shift plus interpolation (Yan & Sava 2011). Zhang & McMechan (2010) further investigate wavefield vector decomposition method in the wavenumber domain and apply it to VTI media. Cheng & Fomel (2014) formulate the wave-mode separation and decomposition operators in heterogeneous media as Fourier integral operators (FIOs) and efficiently apply them using the low-rank approximation . Sripanich et al. (2017) extend the low-rank decomposition operator to wave-mode decomposition in orthorhombic media.
Conventionally, wave-mode decomposition and wave extrapolation are considered as two separated steps. Hou et al. (2014) and Cheng et al. ( , 2016 combine these two steps into a single FIO, which can be implemented by low-rank approximation. These methods are based on the assumption that the medium properties are sufficiently smooth so that their spatial derivatives can be neglected. However, the Earth model can be strongly heterogeneous and contain discontinuities, for example, at salt/sediment boundaries. In such cases, the assumption about the smoothness of the Earth model is no longer valid and could lead to inaccurate calculation of polarization directions. More importantly, simultaneous wave extrapolation and wave-mode separation based on such an assumption may fail to provide reliable phase and amplitude information. In this paper, we introduce a general framework for elastic wave extrapolation in strongly heterogeneous and anisotropic media without the assumption of the smoothness of the medium. The proposed method uses FIOs which allow accurate and stable wave extrapolation to be performed without explicit wave-mode separation. The proposed formulation reveals a simple connection between wave-mode decomposition and wave extrapolation through matrix exponentials. We also show that it is not necessary to explicitly decompose the wavefield into separate wave modes in order to apply recursive integral operators. We first derive one-step elastic wave extrapolation in homogeneous and smooth media motivated by one-step acoustic wave extrapolation (Zhang & Zhang 2009;Sun et al. 2016). In one-step extrapolation of elastic waves, only positive or negative frequency components are propagated, naturally providing the directional information that is useful for efficient wavefield up-down separation and angle gathers computation during imaging (Shen & Albertin 2015;Hu et al. 2016;Sun et al. 2016). We also construct the corresponding two-step scheme which uses only a real-valued vector wavefield. We draw connections of the proposed method to simultaneous propagation of decoupled elastic wave modes. Next, we show that to accurately model wave propagation in strongly heterogeneous media, spatial gradients of stiffnesses need to be included in the Christoffel matrix, leading to complex eigenvalues and polarization directions. Efficient calculation of the proposed FIOs in heterogeneous media is enabled by approximating the wave extrapolation matrix symbol with a low-rank decomposition. Our numerical examples demonstrate that the proposed method is stable and free of dispersion artefacts and therefore is suitable for accurate elastic RTM and FWI in (strongly) heterogeneous anisotropic media.

A generic wave equation
Following the notation of Du et al. (2014), a generic linear second order in time wave equation can be expressed in the following form where u is the wavefield, t is time and A is the matrix operator containing material parameters and spatial derivative operators. Eq. (1) can also be expressed using the first-order system ∂ ∂t where u t ≡ ∂u/∂t. The solution of eq. (2) can be formulated using the definition of the matrix exponential: Defining ≡ √ A, the eigenvalue decomposition of B can be written as ) where The solution to the first-order system of eq. (3) can now be written as To simplify the system, we can define the analytical wavefield The solution to eq. (1) finally takes the form Selecting the first of the two decoupled solutions of eq. (8) leads to a one-step time extrapolation operator Using the form of eq. (9), we can also define the backward time extrapolator If we require the data to be real-valued at every time step, and sum the forward and backward extrapolators, we then arrive at the two-step formulation: The wavefieldû is an analytical signal (denoted by ∧ ), with its imaginary part being the Hilbert transform of its real part (Zhang & Zhang 2009). To see this, we can perform Hilbert transform to the real-valued wavefield u in the frequency domain, and use the dispersion relation ω 2 = 2 and derivative property of Fourier Transform where F and F −1 denote forward and inverse Fourier transform in time, and H denotes Hilbert transform. The output corresponds to the imaginary part of analytical wavefieldû.

Fourier integral operator
The FIO is a type of linear operator that is frequently encountered in the solution of differential and pseudo-differential equations, and is widely used in physics. These operators take the form (Candès et al. 2007) where (x, ξ ) is a smooth phase function, a(x, ξ ) is a smooth amplitude term and F(ξ ) is the Fourier transform of the input function f(x): In the case of acoustic isotropic constant-density wave equations for pressure waves, the matrix operator A of eq. (1) is equal to v 2 |k| 2 , where v is velocity and |k| 2 is the Laplacian operator. One (positive frequency) branch of the solution to this partial differential equation corresponds to the one-step extrapolation method (Zhang & Zhang 2009): where k = 2πξ. Eq. (15) is simple to compute when velocity is constant, but becomes a form of FIO when velocity is heterogeneous and/or anisotropic. In such cases, straightforward application of the FIO amounts to a matrix-vector multiplication with a quadratic cost, which is prohibitive for real-data applications. A low-rank approximation to the operator matrix defined by the wave extrapolation symbol is usually required for practical application of such FIOs . For seismic wave extrapolation, the low-rank property is related to the smoothness of the phase function that describes phase velocity of seismic waves.

Homogeneous media
We first assume homogeneous model properties and investigate the form of time-stepping solutions in general elastic media. For elastic anisotropic wave equations, the matrix operator A corresponds to the spatial Fourier transform of the matrix − /ρ = −DCD /ρ, where is the Christoffel matrix, ρ is density, C is the elastic stiffness tensor expressed in Voigt notation as a 6 × 6 matrix and D is the derivative matrix operator given by In the Fourier (wavenumber) domain, the first-order spatial derivative corresponds to ik i , where i can be x, y or z. Therefore, the homogeneous A takes the form˜ /ρ after i 2 cancels the negative sign, that is, − F −→˜ , where ∼ indicates Fourier domain. For example, the Christoffel matrix˜ in the case of orthorhombic anisotropy takes the form: Since the matrix˜ is symmetric positive definite (SPD), it can be diagonalized with its eigenvalues corresponding to the square of phase velocity of separate wave modes and its orthogonal eigenvectors corresponding to the polarization directions (Dellinger & Etgen 1990): Because Q is orthogonal, Q projects the input vector to its column space. The square root matrix ≡ √ A is analogous to the phase function in the acoustic case, and corresponds to the angular frequency ω according to the dispersion relation. It can be found by taking the square root of the eigenvalues in the diagonal matrix. Analogously, the wave extrapolation operator, e i t , can be computed as: An analogous idea of using matrix exponentials for wave propagation was studied previously by Kosloff & Kessler (1987) in application to the one-way wave equation. Physically, the operator in eq. (19) first decomposes the input vector wavefield into three wave modes, phase shifts them using the corresponding phase velocities, and then aligns them in the polarization directions of the decomposed wave modes. The operator defined in eq. (19) can be expressed as a summation of rank-one matrices Note that the a i a i term in eq. (20) is the wave-mode decomposition operator (Dellinger 1991;Zhang & McMechan 2010).
In the wavenumber domain, we can write the wavefield vector asû To apply the Fourier integral operator, we omit the pair of forward and backward Fourier transforms and formally write where s xx = e ivp k t a px a px + e iv s1 k t a s1x a s1x + e iv s2 k t a s2x a s2x , s xy = e ivp k t a px a py + e iv s1 k t a s1x a s1y + e iv s2 k t a s2x a s2y , s xz = e ivp k t a px a pz + e iv s1 k t a s1x a s1z + e iv s2 k t a s2x a s2z , s yy = e ivp k t a py a py + e iv s1 k t a s1y a s1y + e iv s2 k t a s2y a s2y , s yz = e ivp k t a py a pz + e iv s1 k t a s1y a s1z + e iv s2 k t a s2y a s2z , s zz = e ivp k t a pz a pz + e iv s1 k t a s1z a s1z + e iv s2 k t a s2z a s2z , In a two-step formulation (eq. 11), using eq. (20), the cosine term has the following expression: and can be calculated as where c yy = cos(v p k t)a py a py + cos(v s1 k t)a s1y a s1y + cos(v s2 k t)a s2y a s2y , Note that, because the elements of the 3 × 3 matrix cos ( t) in eq. (24) are real-valued, the wavefield will stay real-valued as long as data are real-valued. The matrix cos ( t) is also closely connected with the k-space adjustment of the Christoffel matrix (Liu 1995;Firouzi et al. 2012;Cheng et al. 2016).
Expanding the wave extrapolation operator in the form of eqs (22) and (25) reveals a simple relationship between wave propagation and wave-mode decomposition. In the one-step formulation, each individual term contained in each element of e i t is essentially a sequence of wave-mode decomposition, phase shift and recomposition. For example, the action of the first term in s xx , e ivp k t a px a px , can be interpreted as projecting the x component of a vector elastic wavefield onto the P-wave mode, phase shifting it using the P-wave phase velocity and then aligning it with the x component of the P-wave polarization direction. If individual wave modes are needed to perform imaging, the decoupled operators can be separately applied. The operators concerning a specific wave mode are indicated by the corresponding phase velocity. The first columns in eqs (22) and (25) are P-wave propagators, while the second and third columns are S1 and S2 wave propagators, respectively. In general anisotropic media beyond TTI symmetry, the two S-wave modes do not decouple easily. Therefore, in order to avoid S-wave singularities in wave propagation, the two-coupled S waves should be computed together (Cheng et al.2016;Sun et al. 2016).
The proposed framework, however, does not require explicit wave-mode decomposition for wave extrapolation. Therefore, it can significantly reduce the computational cost of wave extrapolation, yet still obtain waves free of instability and dispersion artefacts. We also emphasize that the proposed method is capable of handling general anisotropic media, including the case of triclinic anisotropy.

Heterogeneous media
In heterogeneous media, the elastic wave equation can be expressed using the Einstein notation: where dot on top denotes time derivative, comma in subscript denotes spatial derivative and repeated indices imply summation. Analogously, the Christoffel matrix in the case of orthorhombic symmetry can be expanded using the chain rule After spatial Fourier transform, that is, − F −→˜ , the Christoffel matrix˜ in the case of orthorhombic anisotropy takes the form: When the model is smoothly varying, which is the underlying assumption of wave-mode separation , the gradients of stiffnesses are insignificant. In such cases, the imaginary matrix can be dropped, which leads to the conventional real-valued Christoffel matrix in eq. (17). However, in the case of strong heterogeneity, such as at medium interfaces, the gradients of stiffnesses become significant and the imaginary part needs to be taken into account. Because˜ becomes complex and non-Hermitian, the eigendecomposition of the generic matrix A is expressed as where the superscript * denotes conjugate transpose. Note that the eigenvalues ν 2 1 , ν 2 2 and ν 2 3 , as well as their corresponding eigenvectors, are complex-valued. Since Q is square and Q −1 Q = I, it follows that Physically, this means that, in strongly heterogeneous media, wave-mode separation cannot be clearly defined due to wave-mode conversion, which is expressed as the imaginary part of eigenvalues. The polarization directions are no longer orthogonal in the real sense-they become orthogonal complex vectors.
Comparing eq. (30) with eq. (21), we can see that because the eigenvalues that appear in the exponent in eq. (30) are no longer real-valued, the corresponding two-step operators to that in eq. (23) cannot be constructed using Euler's formula, as done in eqs (23)-(25).

Low-rank approximation
So far, we have laid out our basic theory of recursive integral time extrapolation of elastic waves. In mildly heterogeneous media, the Christoffel matrix is SPD. In strongly heterogeneous media, the Christoffel matrix becomes complex-valued and non-Hermitian. However, in both cases, the eigenvalues and eigenvectors of the Christoffel matrix become dependent on both spatial location and propagation direction, in other words, they are functions of both space x and wavenumber k. If these operators are implemented straightforwardly, one is faced with the daunting task of computing and storing the complete eigenvalue decomposition of the Christoffel matrix using all the combinations of x and k, leading to O(N 2 x ) computational and memory complexity, where N x refers to the total number of mesh points in 3-D. To perform wave extrapolation in the form of integral operators, one would have to multiply matrices with vectors in dimension of N x , leading to a computational complexity of O(N 2 x ). This is simply infeasible for practical applications. In this work, to efficiently apply the derived FIOs, we proposed to apply the low-rank decomposition ) on the mixed-domain wave extrapolation matrices. Take the wave extrapolation operator in eq. (30), e i t , as an example. We propose to apply low-rank approximation on each individual element of its expansion. For instance, theŝ xx element, which operates on the x component of the input vector wavefield and outputs to the x component of the output vector wavefield, can be approximated as whereŝ xx (x, k m ) = U andŝ xx (x n , k) = V * are sampled representative columns and rows from the original matrixŝ xx (x, k) = W, M and N are the numerical ranks of matrix W and the matrix a mn = A is obtained from minimizing the objective function defined in the Frobenius norm Similarly,ŝ xy andŝ xz can be approximated aŝ The computation ofû x (x, t + t) then becomeŝ b mn e ixkŝ xy (x n , k)û y (k, t)dk is the complexity of one forward or inverse fast Fourier transform (FFT), and N is the numerical rank of the low-rank approximation, which is 1 for homogeneous media and O(1) for heterogeneous media. It is worth noting that the low-rank decomposition is a pre-computation step, that is, it is only performed prior to time-stepping wave extrapolation. Additionally, only the decomposed matrices, instead of model parameters, need to be stored into memory for repeated application. The time cost by low-rank decomposition, along with the on-the-fly eigenvalue decompositions, is marginal compared with the cost of wave extrapolation, since the decomposition only accesses a subsample of the entire matrix. The computational cost of the application of each FIO is mainly dependent on the numerical ranks for the corresponding low-rank approximation. In homogeneous media, the rank is one. For heterogeneous anisotropic media, the rank is largely determined by model complexity. Examples of approximation ranks are shown in Section 3.

N U M E R I C A L E X A M P L E S
In this section, we use various synthetic models to test the low-rank one-step RITE method. The two-step formulations are not discussed here due to its limitation to smooth media.
We first test the accuracy and stability of the proposed method in a two-layer orthorhombic model, and compare it with the conventional finite-difference and pseudospectral methods. We construct a two-layer orthorhombic model on a 100 × 100 × 100 grid with the densitynormalized stiffness tensor (in km 2 s −2 ) of the first layer given as (Schoenberg & Helbig 1997 , and the second layer being a scaled version of the first layer by a factor of 1.8. We apply triangle smoothing in the vertical direction with a radius of 4. The spatial sampling rate in all directions of the grid is 10 m. A displacement source is oriented at 45 • tilt and 45 • azimuth and injected at x = 0.5 km, y = 0.5 km and z = 0.4 km. The source wavelet has a peak frequency of 35 Hz. A wavefield snapshot shown in Fig. 1 was taken at 0.18 s. The wavefield that was modeled by eighth-order finite difference in space and second-order finite difference in time using a time-step size of 1 ms suffers from strong dispersion artefacts (Fig. 1a), while the wavefield modeled by pseudospectral method in space and second-order finite difference in time is free of spatial dispersions. However, increasing the time-step size to 2 ms leads to numerical (b) the low-rank method with (in red) and without (in blue) including the stiffness gradient terms; (c) the pseudospectral method with (in red) and without (in blue) including the stiffness gradient terms and (d) the pseudospectral method (in red) and the low-rank method (in blue), both including the stiffness gradient terms.
instability of both methods, and thus the results are not shown. In comparison, the proposed low-rank RITE method shows an accurate result free from dispersion and instability using increasing time-step sizes of 1, 2, 4 and 8 ms (Fig. 2). The remarkable numerical stability of the proposed method has been observed and theoretically proved in a preceding study of acoustic wave extrapolation (Sun et al. 2016). This experiment confirms that the method remains stable for applications to elastic wave extrapolation. Unlike finite-difference and finite-element methods, the low-rank RITE method is not subject to the Courant-Friedrichs-Lewy (CFL) condition and does not suffer from numerical dispersion artefacts. Next, to test wave extrapolation accuracy at a medium interface, where strong heterogeneity occurs, we use the same two-layer orthorhombic model without smoothing in the vertical direction. Fig. 3 shows wavefield snapshots taken at t = 0.18 s. Because of the strong contrast at the medium interface, the transmitted and reflected waves calculated by the proposed method including the stiffness gradient terms (Fig. 3a) demonstrate notable amplitude and phase differences compared with the waves calculated by ignoring the gradient terms (Fig. 3b). Fig. 4(a) illustrates a shot gather extracted at Z = 100 m modeled by a reference pseudospectral method, with different arrivals marked on the plot. Fig. 4(b) shows overlaid shot gathers extracted at the same locations, which compares the modeled data by the low-rank method, with and without accounting for the stiffness gradient terms. We can observe that for the two low-rank RITE methods, modeled direct P and S arrivals are identical, while reflected P and S arrivals show obvious difference in both amplitude and phase due to the missing of stiffness gradient terms. To validate our observation, we also carry out wave extrapolation using the pseudospectral method by solving first-and  second-order elastic wave equations. Their shot gather, shown in Fig. 4(c), demonstrates the same phenomena. Finally, in Fig. 4(d), we overlay the amplitude-normalized shot gathers from the low-rank RITE solution including the stiffness gradients, and the pseudospectral solution of the first-order equations. They demonstrate almost identical amplitude and phase behaviour. In terms of computational cost, the numerical ranks of the wave extrapolation matrices are 2 when the stiffness gradient terms are not considered, while the ranks increase to 5 when the gradient terms are included.
In the third example, we test the ability of the proposed method to model decomposed wave-mode propagation in a more general anisotropic medium, a triclinic medium. We use the lab measurements from Mah & Schmitt (2003) as the background model, with the  which is then normalized by a mass density of 1.395 kg m −3 . We then add spherical perturbations centred at Z = 400 m, X = 500 m and Y = 500 m to all the density-normalized stiffness tensors to make the models mildly heterogeneous. Using the same mesh size as the previous example, and a time-step size of 1 ms, the proposed method is capable of accurately and separately propagating P and S waves in the triclinic medium (Fig. 5) using eq. (21). The two S-wave modes are propagated together, since they do not decouple easily. Both the P and coupled S waves are free of numerical dispersion and instability. The numerical ranks of the wave extrapolation matrices for this test are 3.
In the final example, we demonstrate one advantage of the proposed method for seismic imaging applications, the accessibility of directional information of wavefield. For a real-valued wavefield, the propagation direction cannot be uniquely defined in the wavenumber domain because of its symmetric wavenumber spectrum about zero (Hu et al. 2016). This is due to the fact that the real-valued wavefield contains both positive and negative frequency components. The analytical wavefield, on the other hand, has decoupled positive and negative frequency components, so its wavenumber spectrum reveals the direction of wavefield propagation. To demonstrate this fact, we construct an orthorhombic model based on the French model (French 1974), with its C11 component shown in Fig. 6. We use eq. (21) to calculate a decoupled P wavefield, whose x component at t = 0.39 s is shown in Fig. 7(a). We select a portion of the wavefield centred at Z = 0.4 km, X = 0.58 m and Y = 0.58 m propagating in the direction negative along all three axes, as shown in Fig. 7(b). Fig. 8(a) shows the wavenumber amplitude spectrum of the wavefield. Its energy is all distributed in the positive octant. One the other hand, the wavenumber amplitude spectrum of the real part of the wavefield shows symmetric distribution of energy about the origin, as demonstrated by Fig. 8(b). Finally, we decouple the upward-and downward-traveling wavefield based on the sign k z , with the upward-traveling P and S waves shown in Figs 9(a) and (b), and downward-traveling P and S waves shown in Figs 9(c) and (d). The horizontal traveling plane waves are the artefact due the abrupt truncation in the wavenumber domain, and can be removed by tapering the amplitude at the truncation angle. It is worth noting that, with the flexibility of the proposed method, wavefield decomposition is not limited to up/down separation. In fact, the wavefield can be decomposed into arbitrary angle bins.

D I S C U S S I O N
Compared with the finite-element method, which is the state-of-the-art method for modeling wave propagation in global seismology, the proposed method has certain advantages and disadvantages, which dictates its application space.
As is demonstrated in the numerical example, the low-rank RITE method is free of dispersion artefacts and numerical instability, even for large time steps beyond the Nyquist sampling limit. This is not to say that the wave extrapolation is free of phase errors for large time steps, but for similar step sizes, the proposed method is able to produce a cleaner wavefield. Other numerical properties possessed by the proposed method that is not available in finite-element method includes: access to separate wave modes during wave propagation, explicit control of wave propagation direction and free of pseudo-shear artefacts for modeling pseudoacoustic wave propagation. For this study, we have employed the tapering absorbing boundary condition (Cerjan et al. 1985). A direction-dependent absorbing boundary condition can also be implemented to reduce artificial reflections from large incident angles (Sun et al. 2016). A perfectly-match-layers implementation has not been investigated in this work.
However, since the application of the proposed method relies on FFTs, it shares certain constraints with regular grid methods, in particular the pseudo-spectral method. First, the free surface boundary condition can be cumbersome to implement, while for the continuous Galerkin finite-element method such as the spectral element method, the free surface boundary condition is naturally satisfied. Additionally, unlike the finite-element method, the regular grid setup will need very fine discretization to account for complex geometry of the interfaces, which may hinder its application when topography plays an important role in recorded data. Another advantage of the finite-element method is its ability to use finer meshes in low-velocity zones, which can help speed up computation. Fortunately, the stability of the proposed method allows for large step sizes even in the presence of low-velocity zones, which helps to ameliorate the issue with regular grid.
In this paper, we have considered second-order displacement equations for elastic wave extrapolation. Nevertheless, the framework presented here is not limited to such a formulation. For example, second-order stress equations can be used, which has the advantage of not calculating spatial derivative of stiffnesses. Other first-order formulations should also be investigated.
Our formulation can be viewed as an extension of the one-step extrapolation method by Zhang & Zhang (2009) to elastic anisotropic media. It has been shown by Bleistein et al. (2008) that the one-step solution is asymptotically true amplitude, that is, it provides the same traveltime and leading order amplitude as conventional acoustic wave propagation. It is reasonable to expect that the new formulation for elastic waves has similar behaviour. However, more rigorous theoretical proof is required to arrive at a definitive conclusion.
In the last numerical example, we have demonstrated that our method is advantageous in providing directional information about the wavefield. This allows for efficient computation of wavefield up-down separation, RTM angle gathers and absorbing boundary conditions (Shen & Albertin 2015;Hu et al. 2016;Sun et al. 2016) in the context of elastic imaging.

C O N C L U S I O N S
We have presented a recursive integral time extrapolation method for modeling elastic wave propagation in general heterogeneous anisotropic media. The one-step formulation involves analytical wavefields that contain either positive or negative frequencies, which provides crucial information about the direction of wave propagation. The two-step formulation, on the other hand, involves only real-valued wavefields. The proposed method employs a low-rank approximation of the wave extrapolation symbol to efficiently apply the FIOs defined by the mixed-domain components of the modified Christoffel matrix. In practice, this reduces the computational cost to a small number of spatial FFTs per time step. The low-rank decomposition only needs to be computed once prior to wave extrapolation. Numerical examples show that the proposed method has superior accuracy and stability compared with conventional finite-difference and pseudospectral methods.