  For over a decade there have been contradictory claims in the literature
about whether the local bulk flow motion of galaxies is consistent or in
tension with the $\Lambda$CDM model. While it has become evident that
systematics affect bulk flow measurements, systematics in the estimators have
not been widely investigated. In this work, we thoroughly evaluate the
performance of four estimator variants, including the Kaiser maximum likelihood
estimator (MLE) and the minimum variance estimator (MVE). We find that these
estimators are unbiased, however their precision may be strongly correlated
with the survey geometry. Small biases in the estimators can be present leading
to underestimated bulk flows, which we suspect are due to the presence of
non-linear peculiar velocities. The uncertainty assigned to the bulk flows from
these estimators is typically underestimated, which leads to an overestimate of
the tension with $\Lambda$CDM. We estimate the bulk flow for the CosmicFlows-4
data and use mocks to ensure the uncertainties are appropriately accounted for.
Using the MLE we find a bulk flow amplitude of $408\pm165 \mathrm{km s}^{-1}$
at a depth of $49\, \mathrm{Mpc} h^{-1}$, in reasonable agreement with
$\Lambda$CDM. However using the MVE which can probe greater effective depths,
we find an amplitude of $428\pm108 \mathrm{km s}^{-1}$ at a depth of $173\,
\mathrm{Mpc} h^{-1}$, in tension with the model, having only a 0.11%
probability of obtaining a larger $\chi^2$. These measurements appear directed
towards the Great Attractor region where more data may be needed to resolve
tensions.
