In this paper, we propose two new exponential integrator Fourier pseudo-spectral schemes for nonlinear Dirac (NLD) equation. The proposed schemes are time symmetric, unconditionally stable and preserve the total energy in the discrete level. We give rigorously error analysis and establish error bounds in the general H m -norm for the numerical solutions of the new schemes applied to the NLD equation. In more details, the proposed schemes have the second-order temporal accuracy and spectral spatial accuracy, respectively, without any CFL-type condition constraint. The error analysis techniques include the energy method and the techniques of either the cut-off of the nonlinearity to bound the numerical approximate solutions or the mathematical induction. Extensive numerical results are reported to confirm our error bounds and theoretical analysis.