In this paper, a numerical method to solve the multi-dimensional spatial fractional Allen–Cahn equations has been investigated. After semi-discretizating the equations, a system of nonlinear ordinary differential equations with a Toeplitz structure is induced. We propose to split the Toeplitz matrix into the sum of a circulant matrix and a skew-circulant matrix, and apply the Strang splitting method. Such a two-level Strang splitting method will reduce the computational complexity to O(qlog q). Moreover, it preserves not only the discrete maximum principle unconditionally but also second-order convergence as well. By introducing a new modified energy formula, the energy dissipation property can be guaranteed. Finally, some numerical experiments are conducted to confirm the theories we put forward.