  We compute the axion-pion scattering $a \pi \to \pi \pi$, relevant for the
axion thermalization rate in the early universe, within unitarized NLO chiral
perturbation theory. The latter extends the range of validity of the chiral
expansion of axion-pion scattering and thus provides a crucial ingredient for
the reliable determination of the relic density of thermal axions, whenever the
axion decoupling temperature is below that of the QCD phase transition.
Implications for cosmological observables are briefly discussed.
