  Future Baryon Acoustic Oscillation surveys aim at observing galaxy clustering
over a wide range of redshift and galaxy populations at great precision,
reaching tenths of a percent, in order to detect any deviation of dark energy
from the $\LCDM$ model. We utilize a set of paired quasi-\Nb\, FastPM
simulations that were designed to mitigate the sample variance effect on the
BAO feature and evaluated the BAO systematics as precisely as $\sim 0.01\%$. We
report anisotropic BAO scale shifts before and after density field
reconstruction in the presence of redshift-space distortions over a wide range
of redshift, galaxy/halo biases, and shot noise levels. We test different
reconstruction schemes and different smoothing filter scales, and introduce
physically-motivated BAO fitting models. For the first time, we derive a
Galilean-invariant infrared resummed model for halos in real and redshift
space. We test these models from the perspective of robust BAO measurements and
non-BAO information such as growth rate and nonlinear bias. We find that
pre-reconstruction BAO scale has moderate fitting-model dependence at the level
of $0.1\%-0.2\%$ for matter while the dependence is substantially reduced to
less than $0.07\%$ for halos. We find that post-reconstruction BAO shifts are
generally reduced to below $0.1\%$ in the presence of galaxy/halo bias and show
much smaller fitting model dependence. Different reconstruction conventions can
potentially make a much larger difference on the line-of-sight BAO scale, upto
$0.3\%$. Meanwhile, the precision (error) of the BAO measurements is quite
consistent regardless of the choice of the fitting model or reconstruction
convention.
