Abstract We study the probability distribution function (PDF) of relative velocity between two different dark matter halos (i.e., pairwise velocity) with high-resolution cosmological N-body simulations. We revisit a non-Gaussian framework to predict pairwise velocity statistics developed in Tinker. We investigate the pairwise velocity PDFs over a wide range of halo masses of 1012.5 ≲ M [h −1 M ⊙] ≲ 1015 and redshifts of 0 < z < 1. At a given set of masses, redshift, and separation length between two halos, our model requires three parameters to set the pairwise velocity PDF, whereas previous non-Gaussian models in the literature assumed four or more free parameters. At length scales of 5 < r [h −1 Mpc] < 40, our model predicts the mean and dispersion of the pairwise velocity for dark matter halos with masses of 1012.5 ≲ M [h −1 M ⊙] ≲ 1013.5 at 0.3 < z < 1 with a 5%-level precision. We demonstrate that our model of the pairwise velocity PDF provides an accurate mapping of the two-point clustering of massive-galaxy-sized halos at scales of O(10)h −1 Mpc between redshift and real space for a given real-space correlation function. For a mass-limited halo sample with masses greater than 1013.5 h −1 M ⊙ at z = 0.55, our model can explain the monopole and quadrupole moments of the redshift-space two-point correlations with a precision better than 5% at the scales of 5−40 and 10–30 h −1 Mpc, respectively. Our model of the pairwise velocity PDF will give a detailed explanation of the statistics of massive galaxies at intermediate scales in redshift surveys.