# An Efficient Finite-Difference Stencil with High-Order Temporal Accuracy for Scalar Wave Modeling
## 2 Method
### Determining FD Coefficients of the New Stencil through Taylor-Series Expansion  
Assuming the plane wave propagates in the grid, we let  
\[p_{i,j}^{l}=e^{i[k_{x}(x+ih)+k_{z}(z+jh)-\omega(t+l\Delta t)]}. \tag{5}\]  
Here, \(k_{x}=k\cos(\theta)\), \(k_{z}=\sin(\theta)\), and \(\  
When comparing the weights of \(k_{x}^{2\xi}k_{z}^{2j-2\xi}\), we obtain  
\[\left\{ \begin{array}{ll}\sum\limits_{m=1}^{M}m^{2j}a_{m}+2\sum \limits_{n=1}^{N-1}n^{n2j}b_{n}=r^{2j-2}\ \ \ (\xi=0\ or\ \xi=j)\\ 2\sum\limits_{n=1}^{N-1}\frac{n^{2j}}{(2j-2\xi)!(2\xi)!}b_{n}= \frac{j!r^{2j-2}}{(2j)!\xi!(j-\xi)!}\ \ (\xi=1,2,\ldots,j-1).\end{array} \right. \tag{11}\]  
We only need \(M+N-1\) independent equations to solve the FD coefficients \(a_{m}\) and \(b_{n}\). Through redundancy analysis, we simplify Equation (11) as  
\[\left\{ \begin{array}{ll}\sum\limits_{m=1}^{M}m^{2j}a_{m}+2\sum \limits_{n=1}^{N-1}n^{2j}b_{n}=r^{2j-2}\ \ \ (j=1,2,\ldots,M);\\ \sum\limits_{n=1}^{N-1}n^{2j}b_{n}=\frac{j!(2\xi)!(2j-2\xi)!}{2(2j)!\xi!(j- \xi)!}r^{2j-2}\ \ \ (j=2,3,\ldots,N);\\ \mbox{and}\ \ \xi=int(j/2)).\end{array} \right. \tag{12}\]  
where \(int\) represents a function to find the integer part of the real value. In Appendix A, we prove that the FD coefficients obtained by Equation (12) can reach at least sixth-order temporal accuracy. Aside from that, we give the solution of Equation (12) in an analytic form when \(N=2\), and the details are shown in Appendix B. Next, we analyze the relationship between the new stencil and some existing methods:  
(1) When \(M>0\) and \(N>0\) (for example, \(M=3\) and \(N=3\)), this is a general situation, and Equation (12) can be expressed as  
\[\left[ \begin{array}{cccccc}1&4&9&2&8&18\\ 1&16&81&2&32&162\\ 1&16&81&2&128&1458\\ 0&0&0&1&16&81\\ 0&0&0&1&64&729\\ 0&0&0&1&256&6561\end{array} \right]\left[ \begin{array}{c}a_{1}\\ a_{2}\\ a_{3}\\ b_{1}\\ b_{2}\\ b_{3}\end{array} \right]=\left[ \begin{array}{c}1\\ r^{2}\\ r^{4}\\ r^{2}/6\\ r^{4}/10\\ 3r^{6}/10\end{array} \right]. \tag{13}\]  
The FD coefficients \(a_{m}\) and \(b_{n}\) can be determined by solving the above linear equation, and \(a_{0}\) can be determined by Equation (8). (2) When \(N=0\), there is no rotated cross-stencil. Equation (12) can be simplified as follows:  
\[\sum\limits_{m=1}^{M}m^{2j}a_{m}=r^{2j-2}\ \ \ \ \ \ (j=1,2,\ldots,M). \tag{14}\]  
The new FD scheme is simplified to the conventional time-space domain FD scheme with a cross-stencil [19]. (3) When \(N=1\), Equation (12) becomes the fourth-order temporal accuracy FD scheme with a cross-rhombus stencil [14; 15]. (4) When \(r=0\), the FD coefficients \(b_{n}=0\), Equation (12) can be expressed as  
\[\left\{ \begin{array}{ll}\sum\limits_{m=1}^{M}m^{2}a_{m}=1&\ \ \ (j=1)\\ \sum\limits_{m=1}^{M}m^{2j}a_{m}=0&\ \ \ (j=2,\ldots,M)\end{array} \right.. \tag{15}\]  
In this case, the new FD scheme is simplified to the space domain FD scheme with a cross-stencil. Determining the FD Coefficients of the New Stencil with the Dispersion Relationship-Preserving Method  
TE-based FD coefficients are prone to dispersion errors for the high-wavenumber components [17; 22]. As an alternative, we present the dispersion relationship-preserving method to obtain FD coefficients satisfying more elevated wavenumbers in the dispersion relationship [6; 7; 20; 21]. The DRP-based FD coefficients can effectively mitigate the dispersion errors from the high-wavenumber components. Following our previous work [7], we define a new function \(\psi_{m,\beta,\theta}\) to represent the weights of the FD coefficients \(a_{m}\) in the dispersion relation (Equation (6)). Then, \(\psi_{m,\beta,\theta}\) can be defined as  
\[\psi_{m,\beta,\theta}=2[\cos(m\beta\cos(\theta))+\cos(m\beta\sin(\theta))]. \tag{16}\]  
Here, \(\beta=kl\), and \(\theta\) represents the propagation angle. Similarly, we define another function \(\varphi_{n,\beta,\theta}\) to represent the weights of \(b_{n}\).