## 1 Introduction  
Regarding the numerical modeling of wave propagation in anisotropic media, numerous studies have been implemented using finite difference and spectral methods in the time domain [1, 2, 3] to obtain the numerical solution in discrete geological media. These media may contain rocks with various properties, defined as acoustic, elastic, viscoelastic isotropic, and anisotropic. With the flourishing development of computational technology and numerical algorithms, wavefield simulation methods have played a crucial role in imaging subsurface characteristics. Various imaging methods, such as tomography [4, 5], reverse time migration [6, 7, 8], and full-waveform inversion [9, 10, 11], require an efficient wavefield simulation method to generate synthetic waveform data to match the observed data. Accurately performing wavefield modeling with high efficiency in complex media, which contains anisotropy and attenuation, is critical. In contrast to time-domain wavefield modeling, although it is challenging to solve the sparse linearized matrix, multi-source modeling can be performed efficiently, and attenuation can be easily introduced [12]. Inaddition, the physical parameters of the reservoir can be estimated by seismic inversion and imaging methods within a limited number of frequency components [13]. [PERSON] et al. [14] developed frequency-domain finite-difference viscoacoustic wavefield modeling in tilted transversely isotropic (TTI) media. [PERSON] et al. [15] used the finite element method to simulate the frequency-domain wavefield of vertically transversely isotropic (VTI) media. [PERSON] et al. [16] implemented seismic wave modeling in elastic anisotropic media using the generalized stiffness reduction method. [PERSON] et al. [17] developed a new version of the generalized stiffness reduction method for anisotropic media and incorporated it into the 2.5D spectral element method. [PERSON] et al. [3] derived attenuating anisotropic wave equations involving fractional time derivatives and solved them using the pseudospectral method. Reducing artificial reflections is still essential despite numerous frequency-domain modeling methods. Figure 1 shows a sketch of the absorbing boundary condition (ABC). The middle area is the computational area, and the dashed-dot area is the truncated boundary. When the wavefield propagates to the boundary, the ABC needs to attenuate the reflections inwards. Since [PERSON] [18] developed a perfectly matched layer (PML) to reduce artificial reflections of the electromagnetic wavefield, many studies have shown success in its applications to acoustic and elastic wave modeling [19; 20], and improvements have also been made to the perfectly matched layer method [21]. The PML method has become popular for wave modeling problems because it is applicable to first-order wave equations in the time domain [22] and frequency domain [23]. However, it is not straightforward how to implement it for second-order wave equations [24]. In addition, the most severe problem of PML may exist in absorbing artificial reflections in elastic anisotropic media. This is because PML can only function effectively when all slowness vector components of the waves are parallel to their corresponding group velocity vectors [25]. As is known, an elastic anisotropic wavefield generates multiple modes of waves and becomes more complex when the PML becomes unstable in dealing with the edge-effect problem. Thus, we need to replace PML in wavefield simulations of elastic anisotropic media. A multi-axial perfectly matched layer (M-PML) boundary condition was developed in the time domain [26], and its stability has been verified in time-domain isotropic elastic modeling [27; 28]. In this paper, the M-PML boundary condition is applied to frequency-domain finite-difference modeling in several elastic anisotropic media.