We compute the axion-pion scattering a π →π π , relevant for the axion thermalization rate in the early Universe, within unitarized next-to-leading-order 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.