  A group of high order Gautschi-type exponential wave integrators (EWIs)
Fourier pseudospectral method are proposed and analyzed for solving the
nonlinear Klein-Gordon equation (KGE) in the nonrelativistic limit regime,
where a parameter $0<\varepsilon\ll1$ which is inversely proportional to the
speed of light, makes the solution propagate waves with wavelength
$O(\varepsilon^2)$ in time and $O(1)$ in space. With the Fourier pseudospectral
method to discretize the KGE in space, we propose a group of EWIs with designed
Gautschi's type quadratures for the temporal integrations, which can offer any
intended even order of accuracy provided that the solution is smooth enough,
while all the current existing EWIs offer at most second order accuracy. The
scheme is explicit, time symmetric and rigorous error estimates show the
meshing strategy of the proposed method is time step $\tau=O(\varepsilon^2)$
and mesh size $h=O(1)$ as $0<\varepsilon\ll1$, which is `optimal' among all
classical numerical methods towards solving the KGE directly in the limit
regime, and which also distinguish our methods from other high order approaches
such as Runge-Kutta methods which require $\tau=O(\varepsilon^3)$. Numerical
experiments with comparisons are done to confirm the error bound and show the
superiority of the proposed methods over existing classical numerical methods.
