  A multiscale time integrator Fourier pseudospectral (MTI-FP) method is
proposed and analyzed for solving the Klein-Gordon-Schr\"{o}dinger (KGS)
equations in the nonrelativistic limit regime with a dimensionless parameter
$0<\varepsilon\le1$ which is inversely proportional to the speed of light. In
fact, the solution to the KGS equations propagates waves with wavelength at
$O(\varepsilon^2)$ and $O(1)$ in time and space, respectively, when
$0<\varepsilon\ll 1$, which brings significantly numerical burdens in practical
computation. The MTI-FP method is designed by adapting a multiscale
decomposition by frequency to the solution at each time step and applying the
Fourier pseudospectral discretization and exponential wave integrators for
spatial and temporal derivatives, respectively. We rigorously establish two
independent error bounds for the MTI-FP at $O(\tau^2/\varepsilon^2+h^{m_0})$
and $O(\varepsilon^2+h^{m_0})$ for $\varepsilon\in(0,1]$ with $\tau$ time step
size, $h$ mesh size and $m_0\ge 4$ an integer depending on the regularity of
the solution, which imply that the MTI-FP converges uniformly and optimally in
space with exponential convergence rate if the solution is smooth, and
uniformly in time with linear convergence rate at $O(\tau)$ for
$\varepsilon\in(0,1]$ and optimally with quadratic convergence rate at
$O(\tau^2)$ in the regime when either $\varepsilon=O(1)$ or $0<\varepsilon\le
\tau$. Thus the meshing strategy requirement (or $\varepsilon$-scalability) of
the MTI-FP is $\tau=O(1)$ and $h=O(1)$ for $0<\varepsilon\ll 1$, which is
significantly better than classical methods. Numerical results demonstrate that
our error bounds are optimal and sharp. Finally, the MTI-FP method is applied
to study numerically convergence rates of the KGS equations to the limiting
models in the nonrelativistic limit regime.
