In laboratory fusion devices radio frequency electromagnetic waves are routinely used for heating plasmas and for controlling current profiles. The evolution of particle distribution function in the presence of electromagnetic waves is derived from fundamental equations using the action-angle variables of the dynamical Hamiltonian. Unlike conventional quasilinear theories (QLTs), the distribution function is evolved concurrently with the particle motion. Since the particle dynamics is time reversal invariant, the master equation for the evolution of the distribution function is also time reversal invariant. A sequential averaging of the master equation over the angles leads to a hierarchy of diffusion equations. The diffusion operator in the equation obtained after averaging over all angles is time dependent, in direct contrast to time independent diffusion operator in QLTs. The evolution of the distribution function with time-dependent diffusion operator is markedly different from quasilinear evolution and is illustrated for current drive by a spectrum of coherent electrostatic waves. A proper description of wave–particle interactions is important for fusion plasmas since the velocity space gradients of the distribution function decisively affect collisional relaxation and the associated transport processes.