One intriguing BSM particle is the QCD axion, which could simultaneously provide a solution to the Strong CP problem and account for some, if not all, of the dark matter density in the universe. This particle is a pNGB of the conjectured Peccei-Quinn (PQ) symmetry of the Standard Model. Its mass and interactions are suppressed by a heavy symmetry breaking scale, $f_a$, whose value is roughly greater than $10^{9}$ GeV (or, conversely, the axion mass, $m_a$, is roughly less than $10^4\ \mu \text{eV}$). The density of axions in the universe, which cannot exceed the relic dark matter density and is a quantity of great interest in axion experiments like ADMX, is a result of the early-universe interplay between cosmological evolution and the axion mass as a function of temperature. The latter quantity is proportional to the second derivative of the QCD free energy with respect to the CP-violating phase, $\theta$. However, this quantity is generically non-perturbative and previous calculations have only employed instanton models at the high temperatures of interest (roughly 1 GeV). In this and future works, we aim to calculate the temperature-dependent axion mass at small $\theta$ from first-principle lattice calculations, with controlled statistical and systematic errors. Once calculated, this temperature-dependent axion mass is input for the classical evolution equations of the axion density of the universe. Due to a variety of lattice systematic effects at the very high temperatures required, we perform a calculation of the leading small-$\theta$ cumulant of the theta vacua on large volume lattices for SU(3) Yang-Mills with high statistics as a first proof of concept, before attempting a full QCD calculation in the future. From these pure glue results, the misalignment mechanism yields the axion mass bound $m_a \geq (14.6\pm0.1) \ \mu \text{eV}$ when PQ-breaking occurs after inflation.