  We study the probability distribution function (PDF) of relative velocity
between two different dark matter halos (i.e. pairwise velocity) with a set of
high-resolution cosmological $N$-body simulations. We investigate the pairwise
velocity PDFs over a wide range of halo masses of $10^{12.5-15}\,
h^{-1}M_{\odot}$ and redshifts of $0<z<1$. At a given set of masses, redshift
and the separation length between two halos, our model requires three
parameters to set the pairwise velocity PDF, whereas previous non-Gaussian
models in the literature assume four or more free parameters. At the length
scales of $r=5-40\, [h^{-1}\, \mathrm{Mpc}]$, our model predicts the mean and
dispersion of the pairwise velocity for dark matter halos with their masses of
$10^{12.5-13.5} \, [h^{-1}M_{\odot}]$ at $0.3 < z < 1$ with a 5%-level
precision, while the model precision reaches a 20% level (mostly a 10% level)
for other masses and redshifts explored in the simulations. 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 the scales of $O(10)\,
h^{-1}\mathrm{Mpc}$ between redshift and real space for a given real-space
correlation function. For a mass-limited halo sample with their masses greater
than $10^{13.5}\, h^{-1}M_{\odot}$ at $z=0.55$, our model can explain the
monopole and quadropole 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}\mathrm{Mpc}$, respectively. Our model of the pairwise velocity PDF will
give a detailed explanation of statistics of massive galaxies at the
intermediate scales in redshift surveys, including the non-linear
redshift-space distortion effect in two-point correlation functions and the
measurements of the kinematic Sunyaev-Zel'dovich effect.
