Context. Neutron stars are surrounded by ultra-relativistic particles efficiently accelerated by ultra-strong electromagnetic fields. These particles copiously emit high-energy photons through curvature, synchrotron and inverse Compton radiation. To date, however, no numerical simulations have been able to handle such extreme regimes of very high Lorentz factors and magnetic field strengths close to or even above the quantum critical limit of 4.4 × 109 T. Aims. It is the purpose of this paper to study particle acceleration and radiation reaction damping in a rotating magnetic dipole with realistic field strengths of 105 T–1010 T typical of millisecond and young pulsars and of magnetars. Methods. To this end, we implemented an exact analytical particle pusher including radiation reaction in the reduced Landau–Lifshitz approximation where the electromagnetic field is assumed constant in time and uniform in space during one time step integration. The position update is performed using a velocity Verlet method. We extensively tested our algorithm against time independent background electromagnetic fields like the electric drift in cross electric and magnetic fields and the magnetic drift and mirror motion in a dipole. Finally, we apply it to realistic neutron star environments. Results. We investigated particle acceleration and the impact of radiation reaction for electrons, protons, and iron nuclei inserted around millisecond pulsars, young pulsars, and magnetars, in comparison to situations without radiation reaction. We found that the maximum Lorentz factor depends on the particle species, but only weakly on the neutron star type. Electrons reach energies up to γe ≈ 108 − 109, whereas protons reach energies up to γp ≈ 105 − 106 and iron up to γ ≈ 104 − 105. While protons and iron are not affected by radiation reaction, electrons are drastically decelerated, reducing their maximum Lorentz factor by four orders of magnitude. We also found that the radiation reaction limit trajectories agree quite well with the reduced Landau–Lifshitz approximation in almost all cases.