Ion–electron two-stream instabilities in high intensity heavy ion fusion beams, described self-consistently by the nonlinear Vlasov–Maxwell equations, are studied using a three-dimensional multispecies perturbative particle simulation method. Large-scale parallel particle simulations are carried out using the recently developed Beam Equilibrium, Stability, and Transport (BEST) code. For a parameter regime characteristic of heavy ion fusion drivers, simulation results show that the most unstable mode of the ion–electron two-stream instability has a dipole-mode structure, and the linear growth rate decreases with increasing axial momentum spread of the beam particles due to Landau damping by the axial momentum spread of the beam ions in the longitudinal direction.