# Accurate Pseudo-Spectral Acoustic Wave Modelling with Time Dispersion Elimination
## 2 Methodology  
The system of a coupled first-order velocity-stress acoustic equation without source term in a 2D case has the form  
\[ \begin{split}&\rho\partial_{t}v_{x}=-\partial_{x}p\\ &\rho\partial_{t}v_{y}=-\partial_{y}p\\ &\partial_{t}\sigma=-\rho\big{(}\partial_{x}v_{x}+\partial_{y}v_ {y}\big{)}\\ & p=c^{2}\sigma\end{split}  \tag{1}\]  
here \((v_{x},v_{y})\) are \(x\) and \(y\) components of acoustic particle velocity, respectively. \(p\) is acoustic pressure, \(\rho\) is density, and \(c\) is acoustic velocity. \(\sigma\) is an intermediate variable. Following the standard staggered grid [20], the time-stepping scheme of Equation (1) can be written as  
\[ \begin{split}& v_{x}^{n+\frac{1}{2}}=v_{x}^{n-\frac{1}{2}}-\Delta t \partial_{x}p^{n}/\rho\\ & v_{y}^{n+\frac{1}{2}}=v_{y}^{n-\frac{1}{2}}-\Delta t\partial_{y }p^{n}/\rho\\ &\sigma^{n+1}=\sigma^{n}-\Delta t\rho\bigg{(}\partial_{x}v_{x}^{n +\frac{1}{2}}+\partial_{y}v_{y}^{n+\frac{1}{2}}\bigg{)}\\ & p^{n+1}=c^{2}\sigma^{n+1}\end{split}  \tag{2}\]  
the superscript \(n\) indicates the integer time step, \(n+1/2\) and \(n-1/2\) are the indices for the forward and backward half time steps based on time step \(n\). For the spatial staggered grid discretization, the illustration is shown in Figure 1.  
Note that the spatial derivatives defined on the spatial staggered grid in Equation (2) should be computed using Fourier transform, together with spatial shift properties of Fourier transform. The staggered grid-based pseudo-spectral scheme of Equation (2) has the form  
\[ \begin{split} v_{x}^{n+\frac{1}{2}}&=v_{x}^{n-\frac{1} {2}}-\Delta t/\rho F^{-1}\left(ik_{x}e^{ik_{x}\Delta x/2}F(p^{n})\right)\\ v_{y}^{n+\frac{1}{2}}&=v_{y}^{n-\frac{1}{2}}-\Delta t /\rho F^{-1}\left(ik_{y}e^{ik_{y}\Delta y/2}F(p^{n})\right)\\ \sigma^{n+1}&=\sigma^{n}-\Delta t\rho F^{-1}\left( ik_{x}e^{-ik_{x}\Delta x/2}F\bigg{(}v_{x}^{n+\frac{1}{2}}\bigg{)}+ik_{y}e^{-ik_{y} \Delta y/2}F\bigg{(}v_{y}^{n+\frac{1}{2}}\bigg{)}\right)\\ p^{n+1}&=c^{2}\sigma^{n+1}\end{split}  \tag{3}\]  
where symbols \(F\) and \(F^{-1}\) represent forward and inverse Fourier transform, respectively. \(k_{x}\) and \(k_{y}\) are wavenumbers in the \(x\) and \(z\) direction. \(\Delta x\) and \(\Delta y\) are spatial discretization intervals in the \(x\) and \(z\) directions. To reduce boundary reflection, the perfectly matched layer (PML) is applied to enclose the computational domain [21]. Note that the wave wrapping around effect can also be largely eliminated by adopting PML [22]. Following [23], the discrete Equation (3) including PML can be written as  
\[ \begin{split} v_{x}^{n+\frac{1}{2}}&=e^{-\alpha_{x} \Delta t/2}\bigg{(}e^{-\alpha_{x}\Delta t/2}v_{x}^{n-\frac{1}{2}}-\Delta t/ \rho F^{-1}\Big{(}ik_{x}e^{ik_{x}\Delta x/2}F(p^{n})\Big{)}\bigg{)}\\ v_{y}^{n+\frac{1}{2}}&=e^{-\alpha_{y}\Delta t/2} \bigg{(}e^{-\alpha_{y}\Delta t/2}v_{y}^{n-\frac{1}{2}}-\Delta t/\rho F^{-1} \Big{(}ik_{y}e^{ik_{y}\Delta y/2}F(p^{n})\Big{)}\bigg{)}\\ \sigma_{x}^{n+1}&=e^{-\alpha_{x}\Delta t/2}\bigg{(} e^{-\alpha_{x}\Delta t/2}\sigma_{x}^{n}-\Delta t\rho F^{-1}\Big{(}ik_{x}e^{-ik_{x} \Delta x/2}F\bigg{(}v_{x}^{n+\frac{1}{2}}\bigg{)}\bigg{)}\bigg{)}\\ \sigma_{y}^{n+1}&=e^{-\alpha_{y}\Delta t/2}\bigg{(} e^{-\alpha_{y}\Delta t/2}\sigma_{y}^{n}-\Delta t\rho F^{-1}\Big{(}ik_{y}e^{-ik_{y} \Delta y/2}F\bigg{(}v_{y}^{n+\frac{1}{2}}\bigg{)}\bigg{)}\bigg{)}\\ p^{n+1}&=c^{2}\Big{(}\sigma_{x}^{n+1}+\sigma_{y}^{n +1}\Big{)}\end{split}  \tag{4}\]  
here in Equation (4), we split the acoustic pressure field \(\sigma\) into its \(x\) and \(z\) counterparts to implement the PML boundary condition. \(\alpha_{x}\) and \(\alpha_{y}\) are attenuation coefficients along the \(x\) and \(z\) directions, and the specific choice of these coefficients has been well investigated and can be found in many published studies such as those by [23] or [24]. Figure 1: Spatial staggered grid discretization scheme.