This paper first introduces a high-order approximate formula for the Caputo–Hadamard fractional derivative using interpolation theory, and analyzes its convergence and accuracy. Although the formula bears resemblance to the logarithmic L2-1σ formula in form, it differs in the selection of mesh points, making it suitable for more complex models in numerical computations. Subsequently, a numerical method for the Allen–Cahn equation with time Caputo–Hadamard derivative is proposed. The temporal direction is discretized using the newly derived difference formula, while the spatial direction is approximated by the anisotropic nonconforming quasi-Wilson finite element method (FEM). It is proven that the derived scheme has optimal convergence accuracy in the L2-norm and global superconvergence property in the H1-norm. Numerical examples are provided to corroborate the theoretical claims.