  We present 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
catalogue on galaxy number counts relaxes the bounds on the thermal axion mass,
mildly favouring a $\sim 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.
