Abstract: We present the cosmological bounds on the thermal axion mass in an extended cosmological scenario in which the primordial power spectrum of scalar perturbations differs from the usual power-law shape predicted by the simplest inflationary models. The power spectrum is instead modeled by means of a “piecewise cubic Hermite interpolating polynomial” (PCHIP). When using cosmic microwave background measurements combined with other cosmological data sets, the thermal axion mass constraints are degraded only slightly. The addition of the measurements of sigma(8) and Omega(m) from the 2013 Planck cluster catalog on galaxy number counts relaxes the bounds on the thermal axion mass, mildly favoring a similar to 1 eV axion mass, regardless of the model adopted for the primordial power spectrum. However, in general, such a preference disappears if the sum of the three active neutrino masses is also considered as a free parameter in our numerical analyses, due to the strong correlation between the masses of these two hot thermal relics.