  The naturalness problem of PQ symmetry motivates study of the heavy QCD
axion, with masses $m_a >$ 1 MeV generated at scales above the QCD scale, and
low values of the PQ symmetry breaking scale, $f_a$. We compute the abundance
of such axions in a model-independent way, assuming only that they freeze-out
after reheating from inflation, and are not subsequently diluted by new
physics. If these axions decay between neutrino decoupling and the last scatter
era of the Cosmic Microwave Background (CMB), they dilute the neutrinos and
their abundance is constrained by CMB measurements of the energy density in
dark radiation, $N_{\rm eff}$. We accurately compute this bound using a
numerical code to evolve the axion momentum distribution, including many key
processes and effects previously ignored. We assume that the only relevant
axion decays are to final states involving Standard Model particles. We
determine regions of $(m_a, f_a)$ that will give a signal in $N_{\rm eff}$ at
CMB Stage 4 experiments. We similarly compute the $N_{\rm eff}$ bound and CMB
Stage 4 signal for heavy axions that can decay to light mirror photons.
Finally, we compute the bounds on heavy axions with mass below 1 MeV that decay
after the era of CMB last scatter, from their contribution to cold or hot dark
matter or $N_{\rm eff}$ at this era.

High Energy Physics - Phenomenology