Abstract Electromagnetic waves due to electron–positron clouds (bunches), created by cascading processes in pulsar magnetospheres, have been proposed to explain the pulsar radio emission. In order to verify this hypothesis, we utilized for the first time Particle-in-Cell (PIC) code simulations to study the nonlinear evolution of electron–positron bunches dependant on the initial relative drift speeds of electrons and positrons, plasma temperature, and distance between the bunches. For this sake, we utilized the PIC code ACRONYM with a high-order field solver and particle weighting factor, appropriate to describe relativistic pair plasmas. We found that the bunch expansion is mainly determined by the relative electron–positron drift speed. Finite drift speeds were found to cause the generation of strong electrostatic superluminal waves at the bunch density gradients that reach up to E ∼ 7.5 × 10⁵ V cm⁻¹ (E/(m _e c ω _p e ⁻¹) ∼ 4.4) and strong plasma heating. As a result, up to 15% of the initial kinetic energy is transformed into the electric field energy. Assuming the same electron and positron distributions, we found that the fastest (in the bunch reference frame) particles of consecutively emitted bunches eventually overlap in momentum (velocity) space. This overlap causes two-stream instabilities that generate electrostatic subluminal waves with electric field amplitudes reaching up to E ∼ 1.9 × 10⁴ V cm⁻¹ (E/(m _e c ω _p e ⁻¹) ∼ 0.11). We found that in all simulations the evolution of electron–positron bunches may lead to the generation of electrostatic superluminal or subluminal waves, which, in principle, can be behind the observed electromagnetic emissions of pulsars in the radio wave range.