  Axions, if realized in nature, can be copiously produced in the early
universe via thermal processes, contributing to the mass-energy density of
thermal hot relics. In light of the most recent cosmological observations, we
analyze two different thermal processes within a realistic mixed
hot-dark-matter scenario which includes also massive neutrinos. Considering the
axion-gluon thermalization channel we derive our most constraining bounds on
the hot relic masses $m_a < 7.46$ eV and $\sum m_\nu< 0.114$ eV both at 95 per
cent CL; while studying the axion-pion scattering, without assuming any
specific model for the axion-pion interactions and remaining in the range of
validity of the chiral perturbation theory, our most constraining bounds are
improved to $m_a<0.91$ eV and $\sum m_\nu< 0.105$ eV, both at 95 per cent CL.
Interestingly, in both cases, the total neutrino mass lies very close to the
inverted neutrino mass ordering prediction. If future terrestrial double beta
decay and/or long baseline neutrino experiments find that the nature mass
ordering is the inverted one, this could rule out a wide region in the
currently allowed thermal axion window. Our results therefore strongly support
multi-messenger searches of axions and neutrino properties, together with joint
analyses of their expected sensitivities.
