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 ΛCDM model. We utilize a set of paired quasi-N-body FastPM simulations that were designed to mitigate the sample variance effect on the BAO feature and evaluated the BAO systematics as precisely as ∼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.