  We revisit the joint constraints in the mixed hot dark matter scenario in
which both thermally produced QCD axions and relic neutrinos are present. Upon
recomputing the cosmological axion abundance via recent advances in the
literature, we improve the state-of-the-art analyses and provide updated bounds
on axion and neutrino masses. By avoiding approximate methods, such as the
instantaneous decoupling approximation, and limitations due to the limited
validity of the perturbative approach in QCD that forced to artificially divide
the constraints from the axion-pion and the axion-gluon production channels, we
find robust and self-consistent limits. We investigate the two most popular
axion frameworks: KSVZ and DFSZ. From Big Bang Nucleosynthesis (BBN) light
element abundances data we find for the KSVZ axion $\Delta N_{\rm eff}<0.31$
and an axion mass bound $m_a < 0.53 $ eV (i.e., a bound on the axion decay
constant $f_a > 1.07 \times 10^7$ GeV) both at $95\%$ CL. These BBN bounds are
improved to $\Delta N_{\rm eff}<0.14$ and $m_a< 0.16$ eV ($f_a > 3.56 \times
10^7$ GeV) if a prior on the baryon energy density from Cosmic Microwave
Background (CMB) data is assumed. When instead considering cosmological
observations from the CMB temperature, polarization and lensing from the Planck
satellite combined with large scale structure data we find $\Delta N_{\rm
eff}<0.23$, $m_a< 0.28$ eV ($f_a > 2.02 \times 10^7$ GeV) and $\sum m_\nu <
0.16$ eV at $95\%$ CL. This corresponds approximately to a factor of $5$
improvement in the axion mass bound with respect to the existing limits. Very
similar results are obtained for the DFSZ axion. We also forecast upcoming
observations from future CMB and galaxy surveys, showing that they could reach
percent level errors for $m_a\sim 1$ eV.
