  A long-standing issue in peculiar velocity cosmology is whether the
halo/galaxy velocity bias $b_v=1$ at large scale. The resolution of this
important issue must resort to high precision cosmological simulations.
However, this is hampered by another long-standing `sampling artifact' problem
in volume weighted velocity measurement. We circumvent this problem with a
hybrid approach. We first measure statistics free of sampling artifact, then
link them to volume weighted statistics in theory, finally solve for the
velocity bias. $b_v$ determined by our method is not only free of sampling
artifact, but also free of cosmic variance. We apply this method to a
$\Lambda$CDM N-body simulation of $3072^3$ particles and $1200 Mpc/{\rm h}$ box
size. For the first time, we determine the halo velocity bias to $0.1\%$-$1\%$
accuracy. Our major findings are as follows: (1) $b_v\neq 1$ at $k>0.1 h/{\rm
Mpc}$. The deviation from unity ($|b_v-1|$) increases with $k$. Depending on
halo mass and redshift, it may reach $\mathcal{O}(0.01)$ at $k=0.2 h/{\rm Mpc}$
and $\mathcal{O}(0.05)$ at $k\sim 0.3 h/{\rm Mpc}$. The discovered $b_v\neq 1$
has statistically significant impact on structure growth rate measurement by
spectroscopic redshift surveys, including DESI, Euclid and SKA. (2) Both the
sign and the amplitude of $b_v-1$ depend on mass and redshift. These results
disagree with the peak model prediction in that $b_v$ has much weaker deviation
from unity, varies with redshift, and can be bigger than unity. (3) Most of the
mass and redshift dependences can be compressed into a single dependence on the
halo density bias. Based on this finding, we provide an approximate
two-parameter fitting formula.
