Motion of charged particles in electrostatic plasma turbulence is studied analytically and numerically. Analytical study is based on a stochastic differential equation, the generalized Langevin equation, which is derived by the projection operator method with the assumption that the random fluctuation field is a Gaussian and wide sense stationary process. The equation describes the particle velocity as a stochastic process driven by field fluctuations, and is applied to evaluate the velocity diffusion coefficient. The non‐Markovian velocity diffusivity is characterized by the retarded turbulent collision term which contributes to the memory effect of the wave–particle interaction. Test‐particle numerical experiments are carried out by following trajectories of charged particles in Langmuir turbulence. The observed diffusion coefficient in a relatively strong turbulence deviates from predicted values of the quasilinear and the extended resonance broadening theories. The present non‐Markovian formulation explains well such a departure of the diffusion rate.