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.

33 pages, 19 figures, 4 tables. Submitted to MNRAS. Matches version accepted to MNRAS. Moderate changes were made during revision including a comparison between TreePM and FastPM BAO feature