[en] A multi-species perturbative nonlinear (δf) electromagnetic particle simulation scheme has been developed for studying the propagation of intense charged particle beams in high-intensity accelerators and transport systems. The scheme is based on the Darwin approximation of Ampere's law, in which the transverse inductive electric field is neglected, resulting in the elimination of high-frequency transverse electromagnetic effects and, consequently, the associated numerical restrictions from the simulation. However, as noted previously, the presence of the time derivative of the vector potential in the equations of motion for the Darwin model can cause numerical instability. To circumvent this difficulty, we have adopted an approach by replacing the mechanical momentum, p_z, in the direction of beam propagation, by the canonical momentum, P_z=p_z+qA_z/c, as the phase-space variable. The resulting Ampere's law is then modified by the presence of an additional shielding term associated with the skin depth of the species. In order to minimize the numerical noise and to easily access both linear and nonlinear regimes for the physics of interest, we have also adopted the δf formalism for the Darwin model. The absence of unwanted high-frequency waves also enables us to use the adiabatic particle pushing scheme to compensate for the mass-ratio disparities for the various species of charge. The scheme is ideal for studying two-stream and filamentation instabilities, which may cause deterioration of the beam quality in the heavy ion fusion driver and fusion chamber