  The velocity distribution function (VDF) of dark matter (DM) halos in
$\Lambda$CDM dissipationless cosmological simulations, which must be
non-separable in its radial and tangential components, is still poorly known.
We present the first single-parameter, non-separable, anisotropic model for the
VDF in $\Lambda$CDM halos, built from an isotropic $q$-Gaussian (Tsallis) VDF
of the isotropic set of dimensionless spherical velocity components (after
subtraction of streaming motions), normalized by the respective velocity
dispersions. We test our VDF on 90 cluster-mass halos of a dissipationless
cosmological simulation.
  Beyond the virial radius, $r_{\rm vir}$, our model VDF adequately reproduces
that measured in the simulated halos, but no $q$-Gaussian model can adequately
represent the VDF within $r_{\rm vir}$, as the speed distribution function is
then flatter-topped than any $q$-Gaussian can allow. Nevertheless, our VDF fits
significantly better the simulations than the commonly used Maxwellian
(Gaussian) distribution, at virtually all radii within $5\,r_{\rm vir}$. Within
0.4 (1) $r_{\rm vir}$, the non-Gaussianity index $q$ is (roughly) linearly
related to the slope of the density profile and also to the velocity anisotropy
profile. We provide a parametrization of the modulation of $q$ with radius for
both the median fits and the fit of the stacked halo. At radii of a few percent
of $r_{\rm vir}$, corresponding to the Solar position in the Milky Way, our
best-fit VDF, although fitting better the simulations than the Gaussian one,
overproduces significantly the fraction of high velocity objects, indicating
that one should not blindly use these $q$-Gaussian fits to make predictions on
the direct detection rate of DM particles.
