  This paper provides a comprehensive overview of how fitting of Baryon
Acoustic Oscillations (BAO) is carried out within the upcoming Dark Energy
Spectroscopic Instrument's (DESI) 2024 results using its DR1 dataset, and the
associated systematic error budget from theory and modelling of the BAO. We
derive new results showing how non-linearities in the clustering of galaxies
can cause potential biases in measurements of the isotropic
($\alpha_{\mathrm{iso}}$) and anisotropic ($\alpha_{\mathrm{ap}}$) BAO distance
scales, and how these can be effectively removed with an appropriate choice of
reconstruction algorithm. We then demonstrate how theory leads to a clear
choice for how to model the BAO and develop, implement and validate a new model
for the remaining smooth-broadband (i.e., without BAO) component of the galaxy
clustering. Finally, we explore the impact of all remaining modelling choices
on the BAO constraints from DESI using a suite of high-precision simulations,
arriving at a set of best-practices for DESI BAO fits, and an associated theory
and modelling systematic error. Overall, our results demonstrate the remarkable
robustness of the BAO to all our modelling choices and motivate a combined
theory and modelling systematic error contribution to the post-reconstruction
DESI BAO measurements of no more than $0.1\%$ ($0.2\%$) for its isotropic
(anisotropic) distance measurements. We expect the theory and best-practices
laid out to here to be applicable to other BAO experiments in the era of DESI
and beyond.
