Scattering and decay processes of thermal bath particles in the early universe can dump relativistic axions in the primordial plasma. If produced with a significant abundance, their presence can leave observable signatures in cosmological observables probing both the early and the late universe. We focus on the QCD axion and present recent and significant improvements for the calculation of the axion production rate across the different energy scales during the expansion of the universe. We apply these rates to predict the abundance of produced axions and to derive the latest cosmological bounds on the axion mass and couplings.