## 2 Theoretical Foundations
### Numerical Implementation  
Introducing the finite differences and derivatives can be summarized as follows:\[ \begin{array}{c}\frac{\partial u_{x}}{\partial x}\approx\frac{1}{2\Delta x}\Big{[}u _{x}^{i+1,j}-u_{x}^{i-1,j}\Big{]},\\ \\ \frac{\partial u_{x}}{\partial z}\approx\frac{1}{2\Delta z}\Big{[}u_{x}^{i,j+1} -u_{x}^{i,j-1}\Big{]},\\ \\ \frac{\partial^{2}u_{x}}{\partial x^{2}}\approx\frac{1}{\Delta x^{2}}\Big{[}u _{x}^{i+1,j}-2u_{x}^{i,j}+u_{x}^{i-1,j}\Big{]},\\ \\ \frac{\partial^{2}u_{x}}{\partial x^{2}}\approx\frac{1}{\Delta x^{2}}\Big{[}u _{x}^{i,j+1}-2u_{x}^{i,j}+u_{x}^{i,j-1}\Big{]},\\ \\ \frac{\partial^{2}u_{x}}{\partial x\partial z}\approx\frac{1}{4\Delta x\Delta z }\Big{[}u_{x}^{i+1,j+1}-u_{x}^{i+1,j-1}-u_{x}^{i-1,j+1}+u_{x}^{i-1,j-1}\Big{]},\end{array}  \tag{11}\]  
where \(\Delta x\) and \(\Delta z\) are the horizontal and vertical grid sizes, respectively; and the value of grid sizes must be small enough so that there are enough grid points per wavelength [10]; \(i\) and \(j\) are the indices of the horizontal and vertical grids. The finite differences of \(u_{z}\) are similar to the way of the horizontal one. Using the finite-difference star stencil [10]:  
\[ \begin{array}{ccc}\mathbf{M}_{1}&\mathbf{M}_{4}&\mathbf{M}_{7}\\ \mathbf{M}_{2}&\mathbf{M}_{5}&\mathbf{M}_{8}\\ \mathbf{M}_{3}&\mathbf{M}_{6}&\mathbf{M}_{9}\end{array} . \tag{12}\]  
Substituting the finite difference of derivatives and each grid point of Equation (12) according to Figure 2 can be collected as follows:  
After using the finite difference method in the frequency domain, the final equation can be written as Equation (14). Thus, for each frequency component, the matrix equation can be given as:  
\[\mathbf{M}(\omega)\mathbf{U}(\omega)=\mathbf{F}(\omega). \tag{14}\]  
The impedance matrix \(\mathbf{M}(\omega)\) can be obtained by Equation (13), \(\mathbf{U}(\omega)\) is the monochromatic wavefield, and \(\mathbf{F}(\omega)\) is the source term. The frequency-domain Ricker wavelet is used in the following numerical examples, and LU decomposition is used to solve the matrix Equation (14) to get the monochromatic wavefield. Figure 2: The second-order finite-difference star stencil used to model the viscoelastic anisotropic wave equation. at 15 Hz was calculated using classical PML and M-PML, respectively. The asterisk in Figures 3 and 4 is the vertical source. The black dashed line in the figure shows the onset of the absorbing zone. From the results, it can be seen that the PML appears unstable at the boundary intersection in the upper right corner. This phenomenon is due to the complexity of the wavefield in the absorption boundary region, which leads to numerical instability of the PML. Since the PML is a uniaxially defined damping profile, absorption incompleteness is also observed in the other three corners. In the M-PML, however, since the damping profile is multi-axial, the instability that appears in the upper right corner of the PML is well resolved, and the attenuation of the waves is better in the corner regions. Snapshots of the time-domain wavefield corresponding to Figure 3 are given in Figure 4. The black dashed line in Figure 4 is the onset of the absorbing boundary condition, and the source is indicated by an asterisk. Figure 4a,c,e,g are the wavefields calculated by PML, and Figure 4b,d,f,h are the wavefields calculated by M-PML. The results show that the PML-calculated wavefields display a significant instability because at 0.0818 s, the wavefield has not yet propagated to the boundary, but strong energy reflections have already occurred, and these reflections result from the instability remain almost constant as the wavefield propagates. Therefore, PML is unstable in strongly anisotropic media. In contrast, the results of M-PML in the right panel have no unstable reflections, and the waves are successfully absorbed in the truncated boundary regions, so the M-PML with multi-axial damping profiles solves the problem of instability of PML in strongly anisotropic media.