We present numerical simulations of electron–positron pair cascades in the magnetospheres of magnetic neutron stars for a wide range of surface fields (Bp= 1012–1015 G), rotation periods (0.1–10 s) and field geometries. This has been motivated by the discovery in recent years of a number of radio pulsars with inferred magnetic fields comparable to those of magnetars. Evolving the cascade generated by a primary electron or positron after it has been accelerated in the inner gap of the magnetosphere, we follow the spatial development of the cascade until the secondary photons and electron–positron pairs leave the magnetosphere, and we obtain the pair multiplicity and the energy spectra of the cascade pairs and photons under various conditions. Going beyond previous works, which were restricted to weaker fields (B≲ a few × 1012 G), we have incorporated in our simulations detailed treatments of physical processes that are potentially important (especially in the high‐field regime) but were either neglected or crudely treated before, including photon splitting with the correct selection rules for photon polarization modes, one‐photon pair production into low Landau levels for the e±, and resonant inverse Compton scattering from polar cap hotspots. We find that even for B≫BQ= 4 × 1013 G, photon splitting has a small effect on the multiplicity of the cascade since a majority of the photons in the cascade cannot split. One‐photon decay into e+ e− pairs at low Landau levels, however, becomes the dominant pair production channel when B≳ 3 × 1012 G; this tends to suppress synchrotron radiation so that the cascade can develop only at a larger distance from the stellar surface. Nevertheless, we find that the total number of pairs and their energy spectrum produced in the cascade depend mainly on the polar cap voltage BpP−2, and are weakly dependent on Bp (and P) alone. We discuss the implications of our results for the radio pulsar death line and for the hard X‐ray emission from magnetized neutron stars.