Multi-axial perfectly matched layer(M-PML),known to have lost the perfect-matching property owing to multi-axial coordinate stretching,has been numerically validated to be long-time stable and it is thus used extensively in linear anisotropic wave simulation and in isotropic cases where the PML becomes unstable.We are concerned with the construction of the M-PML for anisotropic wave simulation based on a second order wave equation implemented with the displacement-based numerical method.We address the benefit of the incorrect chain rule,which is implicitly adopted in the previous derivation of the M-PML.We show that using the frequency-shifted stretching function improves the absorbing efficiency of the M-PML for near-grazing incident waves.Then,through multi-axial complex-coordinate stretching the second order anisotropic wave equation in a weak form,we derive a time-domain multi-axial unsplit frequency-shifted PML(M-UFSPML)using the frequency-shifted stretching function and the incorrect chain rule.A new approach is provided to reduce the number of memory variables needed for computing convolution terms in the M-UFSPML.The obtained M-UFSPML is well suited for implementation with a finite element or the spectral element method.By providing several typical examples,we numerically verify the accuracy and long-time stability of the implementation of our M-UFSPML by utilizing the Legendre spectral element method.