## 2 Basic Method  
The first-order velocity-stress formulation of the acoustic wave equation is as follows:  
\[\frac{\partial P}{\partial t}=-v^{2}\bigg{(}\frac{\partial v_{z}}{\partial z}+ \frac{\partial v_{x}}{\partial x}\bigg{)} \tag{1}\]  
\[\frac{\partial v_{x}}{\partial t}=-\frac{\partial P}{\partial x} \tag{2}\]  
\[\frac{\partial v_{z}}{\partial t}=-\frac{\partial P}{\partial z} \tag{3}\]  
where \(P\) denotes the pressure, \(v\) is the acoustic wave propagation velocity, \(v_{z}\) and \(v_{x}\) are its propagation velocities in the \(x\) and \(z\) directions, and \(x\) is a spatial coordinate. The first spatial derivative SGFD operator of (2) can be expressed as the following equation [7; 23]:  
\[\frac{\partial P}{\partial x}\approx\frac{1}{h}\sum_{n=1}^{N}a_{n}[P(x+nh-1/2 h)-P(x-nh+1/2h)] \tag{4}\]  
where \(h\) is the spatial step size, \(a_{n}\)(\(n=1\), \(\ldots\), \(N\)) is the SGFD coefficient, and \(N\) is the operator length. According to the plane wave theory, it is possible to make the following equation:  
\[P(x)=P_{0}e^{ikx} \tag{5}\]  
where \(P_{0}\) is a constant, \(i^{2}=-1\), and \(k\) denotes the wavenumber. By substituting (5) into (4), the following equation can be derived [23]:  
\[\sum_{n=1}^{N}a_{n}sin[(2n-1)kh/2]\approx kh/2 \tag{6}\]  
For simplification, (6) can be expressed as follows:  
\[\sum_{n=1}^{N}a_{n}\Psi(\beta)]\approx\psi(\beta) \tag{7}\]  
where  
\[\Psi_{n}(\beta)=sin[(2n-1)\beta],\ \beta=kh/2,\ \psi(\beta)=\beta \tag{8}\]  
and \(\beta\in[0,\mu]\), \(\mu\) is a given numerical value since, under the Nyquist wavenumber condition, \(\beta=\pi/2\), the range of \(\beta\) and \(\mu\) is from 0 to \(\pi/2\). The left side of (7) can be expanded using TEM, which yields the following expression [23]:  
\[\sum_{m=1}^{\infty}\big{[}\frac{(-1)^{m-1}}{(2m-1)!}\sum_{n=1}^{N}(2n-1)^{2m-1 }a_{n}\beta^{2m-1}\big{]}\approx\beta \tag{9}\]By comparing the polynomial coefficients of \(\beta\) on both sides of (9), we can obtain the following equation [23]:  
\[\frac{\left(-1\right)^{m-1}}{\left(2m-1\right)! }\sum_{n=1}^{N}\left[\left(2n-1 \right)^{2m-1}\right]a_{n}\,=\left\{ \begin{array}{cc}1,\ m=1\\ 0,\ m=2,3\ldots N\end{array} \right. \tag{10}\]  
Equation (10) can be transformed into a matrix equation as follows [23]:  
\[ \begin{bmatrix}1&3&\cdots&2N-1\\ 1^{3}&3^{3}&\cdots&\left(2N-1\right)^{3}\\ \vdots&\vdots&\vdots&\vdots\\ 1^{2N-1}&3^{2N-1}&\cdots&\left(2N-1\right)^{2N-1}\end{bmatrix}\begin{bmatrix} a_{1}\\ a_{2}\\ \vdots\\ a_{N}\end{bmatrix} = \begin{bmatrix}1\\ 0\\ \vdots\\ 0\end{bmatrix}  \tag{11}\]  
The formula above can be solved for the SGFD coefficients using the Gaussian elimination method. The Gaussian elimination method is a very classic calculation method in linear algebra that can obtain the SGFD coefficients of (11) by taking the inverse of the coefficient matrix. While the LSM has demonstrated robust capabilities in suppressing numerical dispersion across a broad wavenumber range, as noted by [PERSON] [28], its efficacy diminishes in the small wavenumber range. Conversely, the TEM exhibits exceptional numerical accuracy within the small wavenumber range, which is predicated on the condition that the parameter \(\beta\) approaches zero. In pursuit of developing optimal SGFD coefficients that effectively minimize numerical dispersion across small, medium, and large wavenumber ranges, a novel approach has been employed. This approach involves the integration of TEM with LSM, resulting in the formulation of a combined methodology referred to as the TLM-based SGFD method. This innovative method seeks to harness the strengths of both LSM and TEM, thereby ensuring enhanced precision in numerical simulations across the entire wavenumber range. The relationship delineated in (10), where the exponent of \(\beta\) is established at 1, elucidates a specific mathematical interaction within the SGFD framework.