Dark Radiation Constraints on Heavy QCD Axions

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.


Introduction
The smallest dimensionless parameter of the Standard Model is the strong CP parameter θ 10 −10 .This small parameter can be understood as resulting from a discrete spacetime symmetry, CP [1,2] or P [3][4][5], or from a global Abelian Peccei-Quinn (PQ) symmetry [6,7].Imposing a PQ symmetry appears odd, as the symmetry is necessarily broken by the QCD anomaly; but it may appear as an approximate accidental symmetry at low energies [8].In such scenarios, the PQ symmetry is expected to be explicitly broken by higher-dimensional operators, typically preventing sufficient dynamical relaxation of θ towards zero even if they are suppressed by powers of the Planck mass, M Pl .
For example, if φ is the field that spontaneously breaks the PQ symmetry at scale f a , in standard axion theories the interactions λ n φ n+4 /M n Pl must be suppressed to solve the strong CP problem, λ n < 10 (−46+10n) 10 8 GeV f a 4+n . (1.1) Even for the lowest values of f a allowed by observations, of order 10 8 GeV, operators of dimension 5 through 8 are highly constrained; and the problem gets worse rapidly as f a is increased.It is non-trivial to find theories where an accidental symmetry is protected to such high order, typically requiring significant additions to the theory.Even if one simply imposes the PQ symmetry as a classical symmetry, it may be broken by quantum gravity [9,10], reintroducing the quality problem [11][12][13][14].Since the PQ symmetry must have a QCD anomaly, this PQ quality problem cannot be avoided by promoting it to a gauge symmetry.On the other hand, P and CP can be embedded in higher-dimensional gauge symmetries [15,16] making them attractive avenues for the strong CP problem.The severity of this PQ quality problem, shown in (1.1), applies to the standard QCD axion, where its mass arises from non-perturbative QCD physics at the Fermi scale.It motivates theories with a heavy QCD axion, where the axion mass arises from physics at higher energy scales and is much larger.While the conventional QCD axion mass is less than the eV scale, these theories allow the axion mass to be larger than the MeV scale, removing constraints from stellar cooling and/or beam-dump experiments and allowing greatly reduced symmetry breaking scales.Removing dimension 5 or lower operators by a gauge symmetry, 1which may underlie the accidental PQ symmetry, the shift of θ from zero by a dimension 6 operator is sufficiently small if (1.2) In heavy axion theories there is a limit to how heavy the axion can be, and solving the quality problem then motivates low values of f a .
For f a of order (10 4 -10 7 ) GeV, a strong cosmological limit on the heavy QCD axion arises for masses in the (MeV -GeV) range.Such axions may decay after neutrino decoupling, diluting the neutrino abundance, as found for axion-like particles in [18].There is a powerful bound on this dark radiation from measurements of the Cosmic Microwave Background radiation (CMB) by the Planck Collaboration [19], N eff = 2.96 +0. 34 −0.33 at 95% c.l., and a significantly more accurate determination, with uncertainties smaller by almost an order of magnitude, is a key objective of CMB Stage 4 experiments [20].In this paper we study this bound on (m a , f a ) in a model-independent way, including many effects previously ignored, several arising from axion-meson interactions.See Refs.[21][22][23][24][25][26][27][28] for studies on axions that are light and stable and directly contribute to dark radiation, and Ref. [29][30][31] for axions as hot dark matter.
There is a long history of theories with a heavy QCD axion, motivated by both the quality problem and the interest in reducing f a so that the axion is more visible.One simple possibility is that QCD, or part of the gauge group in which it is embedded, becomes strong in the UV, so that there is an important contribution to the axion potential from short distance instantons [32][33][34][35][36][37].It is important that these instantons do not probe new CP violating phases, so that the new contribution to the potential aligns θ to be sufficiently small.The growth in the QCD coupling in the UV could also arise from extra spatial dimensions [38].
Another simple way to make the QCD axion heavy is to introduce a Z 2 symmetry that transforms the Standard Model (SM) into a mirror sector.The Z 2 symmetry is spontaneously or softly broken so that the mirror electroweak scale is much larger than the SM weak scale, v v.The mirror quarks are then much heavier than the SM quarks, so that below the mirror quark masses the QCD coupling runs faster than the QCD coupling and confines at a scale much above the QCD scale, Λ Λ.When introducing a PQ field that is Z 2 even, the resulting axion couples with the same strength to SM and mirror gluons, and hence its mass is larger than the conventional QCD axion by a factor of roughly (Λ /Λ) 2 .The first implementations of this idea [39,40] used a Weinberg-Wilczek axion [41,42], with the PQ symmetry spontaneously broken by Higgs vevs.In this case the axion decay constant f a is large, of order v , and while these theories ameliorate the quality problem of (1.1, 1.2), solving the problem requires contrived arrangements.On the other hand, in the theories considered in [43][44][45] using a KSVZ axion [46,47], the heavy QCD axion mass is m a ∼ 100 MeV v 10 8 GeV 8/11 10 4 GeV f a (1.3) so that the quality problem is solved by taking f a v .As always, one still needs to understand PQ in operators of dimension ≤ 4 as an accidental consequence of gauge symmetries [8].If the mirror photon in these theories is light, the CMB constraints from dark radiation are modified, which we also study.
The constraints from dark radiation on the axion mass and its couplings have been studied in the literature.Refs.[18,48,49] consider an axion-like particle that couples only to photons, and do not consider axion-gluon couplings.As we will see, the axion-gluon coupling, which leads to axion-meson couplings, helps to keep the axion in thermal equilibrium, so that the Boltzmann suppression of the axion abundance is more effective, relaxing the constraint on (m a , f a ).Ref. [43] studies the mirror QCD case with an axion-mirror photon coupling but, while the decay of the axion into three pions is taken into account, axion-pion scattering and other axion-meson interactions are not included.
This paper is organized as follows.Sec. 2 shows the Lagrangian of the theory above and below the QCD scale.Sec. 3 describes the computation of the dark radiation abundance with a set of Boltzmann equations and shows the resultant N eff .Sec. 4 discusses the case with a mirror photon.Sec. 5 shows a complementary constraint from dark matter overproduction for a sufficiently light axion that decays after the matter-radiation equality, where the N eff constraint is not applicable.We conclude the paper with Sec. 6.

The Effective Theory Above and Below the QCD Scale
In this paper, we study a heavy axion with interactions above the QCD scale given by We assume that the axion couplings with up, down, and electron axial currents are negligible.E/N is the ratio of the electromagnetic and QCD anomalies of the PQ symmetry.For a KSVZ axion with electrically neutral heavy quarks, E/N = 0.For complete representations of SU (5), E/N = 8/3.In models with a large axion mass from mirror QCD, the axion may also couple to mirror photons.We include the effect of mirror photons on N eff in Sec. 4. Note that mirror-photon effects are negligible if the PQ symmetry does not have an electromagnetic anomaly and the mirror quarks are much heavier than the mirror QCD scale, which causes the mixing between the axion and mirror mesons composed of mirror quarks to be small.The analysis of Sec. 2 and 3 is also applicable to the case where the mirror photon is massive and decouples by the QCD phase transition.
After an axial rotation to remove the coupling of the axions to gluons, below the QCD scale, the interactions of the axion with mesons and photons are described by the chiral Lagrangian, which, to leading order in p 2 , is where f π = 93 MeV, B is a strong interaction parameter of order the QCD scale, J a µ is the SU (3) axial current given by and is the quark charge matrix for the transformation that eliminates the axion-gluon coupling with Tr Q A = 1.
The meson nonet and quark mass matrices are given by When calculating the effect of axion-meson scattering at energies below Λ QCD ≈ 150 MeV, we limit ourselves to the relevant two-dimensional subspace, SU (2 since the only active QCD degrees of freedom are pions.For axion masses above Λ QCD , the chiral perturbation based on the symmetry SU (2) L × SU (2) R → SU (2) V breaks down.In this case, the axion may decay to heavy mesons like η and K, or for sufficiently large m a , directly to gluons.Consequently, we use the results of [50] for the axion decay rate into mesons, gluons, and photons (including the enhancement in the axion-photon coupling from a − η(η ) mixing) for m a > m π .We discuss this further in the following section, but for now focus on the SU (2) V subspace which is sufficient for inferring axion-meson scattering in the early Universe.
In this two-dimensional subspace, we take the Q matrix proportional to the identity in isospin space with Q u = Q d = 1/2, so that kinetic mixing between the axion and pion is absent, though mass mixing is present.Expanding out the chiral Lagrangian (2.2) generates the following axion-pion mass matrix and interactions From the π 0 − a mass matrix in the limit f a f π , we can identify B = m 2 π /(m u + m d ).Moreover, since we study axions heavier than the standard QCD axion, m a f a m π f π , and we can thus drop the second term in the bottom right entry of the mass matrix.
The π 0 − a mass matrix is diagonalized by the rotation where r ≡ m a /m π and z ≡ m u /m d 0.47 [51].In terms of the mass eigenstates π0 , â, the interaction of the axion with three pions is described by

10)
where For the remainder of this paper, we drop the hats and refer to the mass eigenstates as a and π 0 .Note that we do not consider the case where m a and m π are so highly degenerate that the axion-pion mixing angle becomes of order unity since this does not occur as long as which is only violated for axions that are extremely degenerate with pions.
The coupling of the axion with photons is ) where E/N arises from the UV contribution associated with the anomalies of the PQ symmetry, 5/3 from the axial rotation that removes the axion couplings to gluons, and F θ from axion-meson mixing.For m a m η , F θ reduces to 2 sin θf a /f π , where θ is the axion-pion mixing angle, (2.9).Further, the value of g γ in the massless axion case is recovered in the limit m a m π in which the term in parenthesis reduces to the standard result, E N − 2 3 4+z 1+z E N −2.03 [51].For m a > m π , we extract F θ from the calculations of [50] which include a − η and a − η mixing.Note that previous considerations of heavy-axion cosmological constraints, except for [45], neglect the effect of axion-meson mixing on F θ which leads to significantly different values of g γ for m a > m π .Finally, we shall consider two reference values: E/N = 8/3, which is the case for the KSVZ model with SU (5) unification, and E/N = 0.

Computation of the Dark Radiation Density
In this section, we present the numerical results of the Boltzmann equations describing the cosmological evolution of the heavy QCD axion in the early universe.Unlike the standard QCD axion which is light and very long-lived, the heavy QCD axion, with m a 1 MeV and f a limited by the quality problem, is cosmologically unstable.When the heavy axion decays during or after neutrino decoupling, photons are subsequently heated relative to neutrinos, producing a potentially observable negative contribution to ∆N eff .Since the heavy axion can be out of thermal equilibrium around neutrino decoupling and contain a population of axions which have large momenta that decay (dangerously) late, it is crucial to track the momentum space distribution function of the axion, f a (p), throughout neutrino decoupling.
While neutrino decoupling occurs around the MeV era, the axion often decouples at earlier times and hence its abundance must be traced back to temperatures far above the MeV scale.The dominant interactions between the axion and thermal bath change as the universe cools.For temperatures above the QCD scale, axion-gluon scattering dominates and ensures the axions reach a thermal distribution for sufficiently high temperatures [52] as discussed in Sec.3.1.For temperatures below the QCD scale, axion-meson scattering and axion-photon scattering can be effective as discussed in Sec.3.2.

Axion Initial Conditions
At temperatures above the QCD phase transition temperature, T QCD ≈ 150 MeV, axion-gluon interactions may be strong enough to keep the axion in thermal equilibrium.Generally, this process is UV dominated so that for sufficiently high temperatures, the axion reaches thermal equilibrium.As the universe cools, depending on f a , the axion-gluon interactions may decouple.Likewise, below T QCD , axion-pion interactions may be strong enough to keep the axion in thermal equilibrium.In this subsection, we compute the temperature at which the axion scattering rate with strongly coupled particles decouples, T FO , and in the following subsection we use this freeze-out temperature to set the initial conditions of the axion distribution function of our Boltzmann code, which evolves the axion phase space distribution from temperatures where first-order chiral perturbation theory is valid, T χPT ≡ 100 MeV [53] to temperatures past neutrino and electron decoupling.We further discuss the region of parameter space in the (m a , f a ) plane where perturbation theory in gluon and pion descriptions breaks down at axion decoupling, and quantify the resulting uncertainty in the initial axion abundance.

Equilibrium from Scatterings
When the temperature T is much greater than m a and T QCD , the axion-gluon interaction a + g ↔ g + g dominates axion production.The thermally averaged rate of axion-gluon scattering is given by [27,52] Γ ag↔gg 16 π Figure 1.Overview of the parameter space in the m a −f a plane where axion interactions with mesons or gluons thermalize.In the purple region, scatterings with pions keep the axion in equilbrium below T χPT .In the yellow region, meson or gluon decays and inverse decays keep the axion in equilibrium below T χPT .In the blue region, scatterings with gluons once kept the axion in equilbrium.In this region, the horizontal contours denote the gluon freeze-out temperature, T FO,g .The red region indicates where the axions freeze out at temperatures above the validity of chiral perturbation theory (T χPT ≡ 100 MeV) but below the validity of gluon peturbation theory (T gPT ≡ 2 GeV).In this region, the freeze-out temperature is uncertain.In each region, we choose the appropriate initial condition of the axion distribution at T χPT ; see Eq. 3.7 for details.
where F g (T ) is a temperature-dependent function that captures the axion production enhancement in the plasma from thermal gluon decays.F g (T ) is numerically computed in [27,52], and for g 3 1, takes the approximate analytic form F g ≈ 2g 2 3 ln 1.5/g 3 [52,54].On the other hand, when T T QCD , the interaction a + π ↔ π + π dominates axion production.The thermally averaged rate of axion-pion scattering is given by where F π (m a , T ) is a temperature and axion mass dependent function that we compute numerically in Appendix A. For reference, Fig. 16 shows F π as a function of m π /T for a variety of m a .We define the axion decoupling temperature, T FO , when Γ ag↔gg = 3H(T FO,g ) for T T QCD , or Γ aπ↔ππ = 3H(T FO,π ) for T T QCD , where H(T FO ) is the Hubble rate at T FO .Contours of T FO,g in the m a − f a plane are shown by the horizontal lines in the blue-shaded region of Fig. 1.Likewise, the purple region indicates where the axion-pion freeze-out temperature, T FO,π is less than T χPT .The estimation of T FO,g based on Eq. (3.1) breaks down if m a > T FO,g , but we find that it anyway occurs in the red parameter region which possesses greater uncertainty: In the red-shaded region, T FO occurs above the temperature at which chiral perturbation theory breaks down (T χPT ) but below the temperature where the strong coupling constant, g 3 , becomes non-perturbative (T gPT ).We take T χPT 100 MeV, the temperature above which one-loop corrections in chiral perturbation theory become comparable to tree-level results [53].Similarly, T gPT 2 GeV is conservatively associated with the energy scale below which g 3 becomes non-perturbative and one-loop corrections become comparable to tree-level results [50].In the red-shaded region, we cannot precisely determine T FO .In Sec.3.1.3,we evaluate the uncertainty in N eff arising from this uncertainty in T FO .

Equilibrium from Decays and Inverse Decays
Even when T χPT < T FO < T gPT , and hence the freeze-out temperature is uncertain, it is still possible to infer the axion abundance at T χPT for sufficiently large m a .Specifically, if the axion decay rate is greater than Hubble at T χPT , then the axion possesses a thermal distribution at T χPT .For example, when m a T , the decay a → g + g can dominate over axion-gluon scattering when T > T gPT , or, for example, a → 3π can dominate at T = T χPT .The axion decay rate to strongly coupled particles is given by where F c (m a ) is an axion mass dependent function that captures the variety of strongly coupled degrees of freedom the axion can decay to.We use F c (m a ) as numerically computed in [50], which includes the following meson decay channels: a → 3π, ππγ, ηππ, KKπ, η ππ, ρρ, ωω, K * K * , and φφ.For m a ≥ T gPT , the axion to gluon decay rate dominates and F c smoothly interpolates to the perturbative result [50] Last, the axion to photon decay rate is given by where g γ (m a ) is given in (2.13).Note that g γ (m a ) is a function of the axion mass due to the effects of axion-meson mixing as encoded in the mixing function F θ (m a ).We define the axion decay temperature, T decay , when Γ a→QCD + Γ a→γγ = 3H(T decay ).The yellow region of Fig. 1 shows the region where T decay > T χPT .In this region, the axion possesses a thermal distribution when we begin our Boltzmann code at T χPT , even if T FO is uncertain.
Note that if the axion decays far before or after neutrino decoupling, the initial condition of the axion at T χPT becomes insensitive to the calculation of ∆N eff .In particular, below the lower dashed line, Γ a = Γ a→QCD + Γ a→γγ , is so large that the axion always decays before the universe is 10 −2 seconds old.In this region, the neutrinos are still strongly coupled to the thermal bath when the axion decays so that any effect to ∆N eff from the axion is erased, regardless of the initial axion abundance at T χPT .Likewise, above the upper dashed line.Γ a is so small that the axion always decays after the universe is 1 second old.In this region, neutrinos have long since decoupled from the thermal bath when the axion decays, leading to ∆N eff −0.3, which is already excluded by the observations of CMB [19].

Initial Condition and Its Uncertainty
For (m a , f a ) in the purple or yellow regions of Fig. 1, the axion is in thermal equilbrium at T χPT .Consequently, in these regions, we take the initial distribution function of the axion at T χPT to be a Bose-Einstein distribution of temperature T χPT .For sufficiently large f a , however, the axion decouples from the bath at T FO > T χPT (blue region), and we take the initial axion distribution to be a thermal one at T FO , red-shifted down to T χPT .Inside the red-shaded region, however, T FO is uncertain.Despite this uncertainty, we can still bound ∆N eff by running our axion Boltzmann code with both the maximum and minimum possible axion abundance at T χPT .We scan over possible freeze-out temperatures between T χPT and T gPT for all axion masses in the red-shaded region to determine the smallest and largest abundance at T χPT as a function of m a .In Fig. 2, the left panel shows the T FO that gives the maximum (blue) and minimum (orange) axion yield, Y a = n a /s, as a function of m a , and the right panel shows these maximum and minimum yields.Note that the smallest possible axion number density is not necessarily that of an axion in thermal equilibrium at T χPT and the largest that of a frozen-out abundance at T FO = T gPT .This is because large changes in the degrees of freedom of the thermal bath between T χPT and T gPT can dilute the previously frozen-out axion.For all future plots, we show the results of ∆N eff arising from these two possible initial conditions.In summary, we take the initial axion distribution function at T χPT (time t χPT ) to be where f a (p a ) is a Bose-Einstein distribution with momentum p a and effective temperature T * given by Note that if T RH < T * or if there is a source of dilution in the universe after the axion freezesout, then the initial abundance of axions can be small and the bounds on N eff discussed in this work are weakened.

Axion Abundance Below Λ QCD : Boltzmann Equations
Accurately capturing the effect of heavy axion decoupling and decay on our cosmology requires understanding the phase space evolution of the axion in the primordial thermal bath.The Boltzmann equation describing the evolution of the axion phase space density, f a (p a ), is where p a = |p a | is the magnitude of the axion momentum, H is the Hubble expansion rate, and C γ , C P , C π , and C Γ are the collision terms for axion-two photon scattering, axion-Primakoff scattering, axion-pion scattering, and axion-meson decay, respectively.Generally, the collision term, C, corresponding to the axion interaction a where E a = p 2 a + m 2 a is the axion energy, dΠ = d 3 p/(2π) 3 2E is the phase space measure per particle, |M | 2 is the matrix element of the interaction, S = 1/m! is the symmetry factor for every m identical particles in the initial or final states, and is the phase space density factor for all the incoming and outgoing particles interacting with the axion, where the plus sign refers to stimulated emission (boson) and the minus to Fermi blocking (fermion).The second line of Eq. (3.10) shows Λ in the limit where 1 ± f 1 with f = exp(−E/T ) the Boltzmann distribution, which is an excellent approximation for particles in kinetic equilibrium.
In past literature, only the axion-two photon, C γ , and the axion-Primakoff, C P , collision terms have been considered in calculations involving the axion Boltzmann equation.Consequently, these terms have already been computed, and their values are [55] 1 − e −(Ea−pa)/2T Γ a→γγ (3.11) Here, Γ a→γγ is the axion-to-two photon decay rate, (3.5), m γ is the photon plasma mass, and n i is the number density of the ith electromagnetically charged particle of charge e i in the thermal bath.When the electron is relativistic, m γ eT /3 [56], but when the electron becomes non-relativistic at T m e , m γ reduces to the classical plasma frequency of m γ e n e /m e , which is exponentially suppressed.For simplicity, we piecewise-connect the two regimes for m γ when they intersect, which occurs roughly at T m e /2.
The axion-pion scattering collision term for a massive axion has not been computed in the literature and we do so for the first time in Appendix A. C π takes the form where A = 1 3 (1 − z)/(1 + z) 0.12 and r ≡ m a /m π as before.The function F PS contains the phase space integration over the axion-pion scattering matrix.In Appendix A, we determine this axion-pion scattering matrix, numerically perform the phase integration, and show how C π agrees with the massless axion result found in literature [53,57].Note that integration of C π over the axion phase space defines Γ aπ↔aπ given in Eq. (3.2).
The axion-meson decay collision term has also never been considered in the literature.Under the Boltzmann approximation in the kinetic equilbrium limit, C Γ is simply where Γ QCD is the total axion decay rate to all final states containing mesons (3.3).In (3.14), we use the Boltzmann kinetic equilibrium approximation.
It is fruitful to estimate the impact of the axion-pion scattering and axion-meson decays compared with the standard Boltzmann calculations in literature which only include Primakoff scattering and axion-photon decays.For example, the solid contours in Fig. 3 show the axionpion decoupling temperature, T FO,π vs f a over a range of axion masses.For f a 10 6 GeV,  the axion remains in thermal equilibrium from pionic interactions until a few 10s of MeV which is typically far lower than the Primakoff decoupling temperature [18], where g c is the the sum of the charged relativistic degrees of freedom in the bath.The Primakoff decoupling temperature is shown by the dashed contours in the left panel of Fig. 3.For f a 5 × 10 4 GeV, T FO,P T FO,π , demonstrating the importance of the pions at maintaining thermal equilibrium with the axion all the way to O(10) MeV temperatures.The lower pionic decoupling temperature compared to standard Primakoff decoupling leads to two important effects: (1) it can reduce the abundance of axions with masses above T FO,π , as they now follow an exponentially suppressed distribution relative to the standard, non-Boltzmann suppressed distribution; and (2) it can enhance the initial abundance of axions with masses below T χPT as they are not diluted by g * S like axions that decouple from Primakoff interactions in the early universe.
Similarly, axion-meson decays and inverse decays can keep the axion in thermal equilbrium at low temperatures.Most importantly, for a fixed f a , the additional QCD decay channels can significantly decrease the axion lifetime relative to Γ −1 a→γγ , which is the standard axion lifetime taken in previous Boltzmann calculations.For example, the solid blue and orange contours of Fig. 4 show the total axion decay rate including QCD decay channels for E/N = 0 and 8/3, respectively.Compared to the dashed contours used in past literature, which show Γ a→γγ when meson mixing and QCD channels are absent, the realistic total axion decay rate can be significantly different.Moreover, because the axion lifetime relative to neutrino decoupling dominantly sets the N eff signal, we expect that incorporating axion-to-QCD channels will significantly alter the allowed parameter space in the (m a , f a ) plane.
To precisely quantify these new effects, the phase space evolution of the axion, (3.8), together with the evolution of the Standard Model particles in the thermal bath must be computed to determine the effect of axion decoupling and decay on the relative temperature differences between photons and neutrinos in our present Universe, as typically characterized by the effective number of neutrino species Here, ρ ν and ρ γ are the neutrino and photon energy densities.Note that N eff is most sensitive to heavy axions that decay at temperatures near neutrino decoupling, which occurs around the MeV scale.The particles in the thermal bath from T χPT through neutrino decoupling are photons, neutrinos, electrons, muons, pions and a small density of heavier mesons.The contribution to the energy density from these heavier mesons and from deviations of pions from the ideal gas law from self-interactions, ρ δ QCD , is approximately 5% of the total energy density at T χPT and quickly drops far below 1% by T = 60 MeV (see Fig. 17).
The evolution of the energy density of species in the thermal bath that are tightly thermally coupled follows the energy density Boltzmann equation i=γ,e,µ,π,δ QCD where are the energy densities and pressures of the ith tightly coupled species in the bath.We assume that the photons and pions follow Bose-Einstein distributions and the electrons and muons follow Fermi-Dirac distributions.As discussed more in Appendix C, we infer ρ δ QCD and P δ QCD from the work of [58], which computes the Standard Model equation of state across the QCD phase transition and takes into account the deviations from the ideal gas law arising from the strongly coupled QCD bath.
The energy densities of particles not strongly thermally coupled electromagnetically, namely the axion and neutrinos, must be solved for numerically.Specifically, the rate of change of the neutrino energy densities on the right side of (3.17) follow the Boltzmann equations Here, Γ νe 0.68 G 2 F T 5 and Γ νµτ 0.15 G 2 F T 5 [18] are thermally averaged neutrino interaction rates with the thermal bath for electron neutrinos and for muon and tau neutrinos, respectively.
Last, the Hubble rate, H, quantifies the expansion rate of the universe and sets the decoupling time of all interactions.The squared Hubble rate is set by the sum of all energy densities, We numerically solve the system of equations (3.8)-(3.21)using the method of lines [59].The method of lines is a numerical technique for solving a system of partial differential equations by discretizing one independent variable direction (comoving momentum in our case) while keeping the other independent variable continuous (logarithmic time in our case).The main advantage of the method of lines technique is the conversion of the Boltzmann system of partial differential equations in {|p a |, t} into a system of many ordinary differential equations in {t} which is computationally easier to solve.Moreover, by keeping the time-like variable continuous, useful techniques such as dynamical step-sizes can be employed to speed up the computation by automatically taking large temporal time-steps when changes in the interactions are small (such as at thermal equilibrium) while taking small temporal time-steps when changes are sudden (such as decays, decouplings, or re-thermalizations).See Appendix C for more details of our numerical setup.

∆N eff
In this section, we present the numerical results of N eff as determined from the Boltzmann equations of Sec.3.2 describing the cosmological evolution of the heavy QCD axion below T χPT .
First, to compare with past literature and highlight the importance of the new effects discussed in this work, we show N eff for heavy axion cosmologies without incorporating the following crucial elements in the Boltzmann code: axion-pion scattering, axion-meson decay, axion-meson mixing in g γ , the proper frozen-out initial axion abundance, and the QCD contributions to the background energy density as described by ρ π and δ QCD .Fig. 5 shows the numerical results of N eff when these terms are neglected.That is, including only the photon (C γ ) and Primakoff (C P ) collision terms in Eq. (3.8) and Eq.(3.17); including only γ, e and µ in the sum of the thermally coupled species of Eq. (3.17); including only ρ γ , ρ e , ρ µ , ρ νe , ρ νµτ , and ρ a in Hubble (3.21); setting the initial abundance of axions at the start of the Boltzmann code to that of a frozen-out abundance set by Primakoff scatterings such that T * = Max(T FO,P , T χPT ); and finally, in Eq. (2.13), setting E/N = 0 and F θ = (1 − z)/(1 + z), which is the axionmeson mixing contribution in the inapplicable m a m π limit.The left panel of Fig. 5 shows contours of ∆N eff in the m a − τ γ plane, where τ γ is the axion to photon lifetime when axionmeson mixing is neglected.Note that τ γ is the total lifetime of the axion since QCD decay channels are neglected in this particular case.
For τ γ 10 −1 s the axion decays after neutrino decoupling, heating up the photons relative to the neutrinos and giving rise to ∆N eff < 0 as can be seen from the enhanced denominator of Eq. (3.16).For τ γ 10 −1 s and m a 10 MeV, the axion decays sufficiently early that the photons and neutrinos rethermalize before the neutrino decouples.In this scenario, the N eff signal of the heavy axion is absent and ∆N eff 0. For τ γ 10 −1 s and m a 10 MeV, the axion remains in thermal equilibrium past neutrino decoupling, heating up the photons and again giving rise to negative ∆N eff .The right panel of Fig. 5 shows the same ∆N eff contours as the left panel but in the m a − f a plane.Both panels assume the usual hadronic axion with E/N = 0. Taking the GUT motivated value of E/N = 8/3 only slightly shifts the contours in the right panel vertically.
We now consider N eff for heavy axion cosmologies incorporating the new effects included in this work: axion-pion scattering, axion-meson decay, axion-meson mixing in g γ , the proper (and occasionally uncertain) frozen-out initial axion abundance, the QCD contributions to the background energy density as described by ρ π and δ QCD , as well results for the KSVZ E/N = 0 and the GUT motivated E/N = 8/3.Figs. 6 and 7 show the contours of ∆N eff for heavy axions with these additional contributions for E/N = 0 (KSVZ) and E/N = 8/3 (GUT), respectively.In both figures, the left and right panels show the parameter space in the (m a , τ ) and (m a , f a ) planes, respectively.Note that here, τ is the total lifetime of the axion, which begins differing from τ γ for m a m π where axion-meson mixing becomes important and then becomes even more disparate when axion-meson channels open for m a 3m π , as shown in Fig. 4. The solid and dashed blue contours in each panel correspond to taking the maximum possible Y a (blue) and minimum possible Y a (dashed) when T FO lies in the uncertain region between T χPT and T gPT , as indicated in Fig. 2. The separation between the solid and dashed blue contours indicates the uncertainty in N eff arising from the uncertainty in T FO in this region.As can be seen, this region is localized roughly between 250 MeV m a 800 MeV and, for any value of m a in this region, the uncertainty in the value of f a for any N eff contour is typically only a several 10s of percent and always less than a factor 3.
Figs. 6 and 7 demonstrate three important differences in the N eff signal from heavy axions as currently considered in literature (Fig. 5): First, the effect of axion-pion resonance on the mixing angle can be seen by the triangular shaped peaks near m a = m 0 π .In this regions, the axion is tightly coupled thermally to pions so that when the axion decays, its abundance is sufficiently exponentially suppressed that it does not heat up the photons even when decaying past neutrino decoupling.Second, for m a 3m π , |∆N eff | is reduced for fixed (m a , f a ) due to the meson-decay channels now open which cause the axion to decay earlier, especially near resonances in the mixing angle at m a = m π , m η , and m η .This can be seen more clearly in the (m a , f a ) planes.Note that for some m a , incorporating axion-meson mixing can increase the axion lifetime due to cancellations between contributions of g γ (3.5).The increased |∆N eff | for fixed (m a , f a ) due to the increased axion lifetime is important for m a < 3m π when only the axion-photon decay channel is open.Third, for m a m π , |∆N eff | is reduced for fixed (m a , f a ) because the axion is kept in thermal equilibrium by axion-pion scattering to lower temperatures compared to when the pions are absent.This leads to a reduced axion abundance at neutrino decoupling which reduces the |∆N eff | contribution from the axion.
Each of these effects can be seen more clearly in the top panels of Fig. 8 which show the evolution of the energy densities of the axion and other species in the thermal bath as a function of time and temperature (top horizontal axis) for fixed (m a , f a ) = (100 MeV, 2.5×10 6 GeV) and (1000 MeV, 1.0 × 10 9 GeV), in the top left and top right panels, respectively.The dark colored contours show the evolution of the comoving energy density, X i = ρ i R 4 of the ith species when including the axion-meson interactions while the light colored contours show the same evolution when the axion-meson interactions are absent.The electromagnetic component of the thermal bath, γ, e, µ, π, δ QCD is shown in blue, ν µ,τ in green, ν e in orange, and the axion in red.The dashed red contour shows the comoving energy density of the axion if it were to maintain a thermal distribution for all times.As can be seen from the m a = 100 MeV panel, the axion starts off in thermal equilbrium compared to the case without meson interactions in which the axion possesses a g * S suppressed abundance from  earlier Primakoff freeze-out.Moreover, the axion abundance with meson interactions follows the dashed thermal distribution to lower temperatures than without meson interactions.This Figure 8. Top panels show the comoving energy density evolution of the axion (dark red), electron neutrinos (dark orange), muon and tau neutrinos (dark green), and strongly coupled species in the Standard Model thermal bath (dark blue) for (m a , τ ) = (100 MeV, 0.1s), left, (m a , τ ) = (1000 MeV, 0.1s), right.The correspondingly lighter shaded contours show the same evolution but without pion or meson interactions.The dashed red contour shows the axion comoving energy density if it were to maintain a thermal density at all times.The bottom panels show the instaneous value of Eq. (3.16) minus the Standard Model result, N SM eff = 3.044.The left and right panels respectively highlight how including QCD interactions can keep the axion thermally coupled to the Standard Model bath for a longer time or enhance the axion decay rate at a fixed f a .The left panel demonstrates that for low m a , the axion can be in thermal equilbrium at T χPT when including pions-interactions but be out of equilbrium when including only Primakoff processes which freeze-out at much higher temperatures that lead the axion to possess a g * S diluted abundance at T χPT .In the left panel, the axion lifetime with or without mesons is approximately the same (τ ≈ 0.1 s) because the meson decay channels are forbidden and axion-pion mixing is not too appreciable yet.In contrast, the right panel demonstrates that for high m a , the axion can decay much earlier due to the kinematic availability of meson decay channels (τ ≈ 0.01 s with meson interactions and 24 s without).This leads to a reduced ∆N eff , as shown on the bottom right panel.In both panels, we take E/N = 0. leads to a relative suppression in the non-relativistic abundance of the axion prior to decaying around 2 MeV.Consequently, including axion-pion scattering ∆N eff is not as negative as previous results in the literature.The temporal evolution of ∆N eff (t) = N eff (t) − 3.044, as shown by the dark blue (with mesons) and light blue (without mesons) contours in the lower left panel of Fig. 8, demonstrates this difference explicitly.Note that for m a = 100 MeV, the axion mass is on the cusp of the axion-pion resonance for E/N = 0.For m a closer to m π , the axion follows the dashed thermal abundance for a longer duration which generates the large triangular peak in the allowed N eff plot of Fig. 6.
For m a = 1000 MeV, the axion with meson interactions again follows the dashed thermal distribution to slightly lower temperatures than the axion without meson interactions.More important though is the difference in decay time between the two cases.In particular, the axion with meson interactions (dark red) decays earlier than the axion without (light red) due to the a → ηππ, ππγ decay channels that are now kinematically open to the 1000 MeV axion that are absent from the 100 MeV axion.These extra decay channels lead to a much smaller ∆N eff as shown explicitly by the evolution of ∆N eff in the bottom right panel of Fig. 8.We note that for m a > 2 GeV, the axion decays dominantly into gluons.Here, the 2σ limit on N eff as constrained by Planck approximately follows the contour τ ≈ 0.05 s around m a ∼ GeV and slowly drops with increasing axion mass.The slight decrease in the maximum allowed τ in this region originates from the increase in the axion energy density at decay with axion mass: For such heavy axions with long lifetimes, f a is large and the axion freezes-out early, leading to freeze-out yield roughly independent of m a .Consequently, the heavier the axion, the earlier it must decay so that its energy density at neutrino decoupling is further exponentially suppressed to counter its larger frozen-out energy density.As shown in Appendix B, the exponentially decaying energy density of axions, ρ a ∝ e −t/τ , leads to a logarithmic decrease in the maximum allowed τ given by the semi-analytic function Thus, for m a 2 GeV and f a small enough to be probed by accelerator experiments [45,50,[60][61][62][63][64][65][66][67][68], the dark radiation constraint is absent.

Including a Mirror Photon
In this section we add a mirror photon γ to the theory, with a mass sufficiently small that it can be ignored in our analysis.A mirror photon is natural in theories with a Z 2 symmetry that not only doubles the SU (3) c sector of the Standard Model to achieve a heavy QCD axion, but also doubles the SU (2) L × U (1) Y sector.This complete mirroring of the Standard Model gauge group introduces another axion coupling relevant in computing the amount of dark radiation Figure 9.The axion total decay width (left) and branching ratios (right) as a function of axion mass, for E/N = 8/3.The orange (blue) contours refer to decays to photons (mirror photons); the green contour is for decays to states that include hadrons.In the left panel, the solid (dashed) black contours give the axion decay width in the theory with (without) mirror photons.
Above the mirror QCD scale, Λ QCD , the Z 2 symmetry ensures that with e differing from e only by renormalization group scaling, which we ignore.Note that unlike g γ in (3.5), g γ does not include contributions from axion-mirror meson mixing nor from the axial rotation onto mirror quarks because the masses of the lightest mirror quarks are typically much heavier than Λ QCD and thus irrelevant to the theory below Λ QCD .If, however, any mirror quark is lighter than the mirror QCD scale, then there is an additional contribution to g γ analogous to the second and third terms in (2.13) for g γ .In the minimal theory, where the Z 2 symmetry exchanges the Standard Model with its mirror and is spontaneously broken by a difference between the electroweak vevs with v v, all mirror quarks are heavier than Λ QCD for f a m a 25 GeV 2 [69].This relation is satisfied for nearly the entire parameter region of interest to us, so that all mirror quarks are well above the QCD scale; g γ is thus uncorrected and given by (4.2).
Numerical results in this section are calculated taking e = e (that is, neglecting the small running of e below v ) and with a non-zero E/N so that g γ is non-zero.As a result, the axion-mirror photon decay rate is In particular, we consider two values for E/N : 8/3, motivated by grand unification, and 1/3, which can be achieved by an appropriate choice of KSVZ fermions.For E/N = 8/3, the total axion decay rate with and without mirror photons is shown by the solid and dashed black contours in the left panel of Fig. 9.The light blue, orange, and green contours indicate the axion decay rates into γ , γ, and QCD degrees of freedom, respectively.Because there is no cancellation of terms in g γ as compared to g γ , the a → γ γ decay rate (blue) is roughly an order of magnitude greater than the a → γγ decay rate (orange) until m a 3m π : thus, axions below this mass dominantly decay into dark photons.This can be seen more clearly in the right panel of Fig. 9, which shows the branching ratios for the same three axion decay channels.The case for E/N = 1/3 is significantly different as demonstrated in Fig. 10.In particular, the cancellation between terms in g γ is negligible which leaves the now smaller a → γ γ decay rate roughly an order of magnitude weaker than the a → γγ decay rate.
We highlight these two representative values of E/N since they generate substantially different decay branching ratios into mirror photons. 2 This disparity is important since the parameter space where the axion branching ratio into dark photons is O(1) can be cosmologically dangerous as the mirror photon decay mode (4.3) increases N eff , by directly generating dark radiation in the form of γ , and reduces the heating of the Standard Model bath as fewer axions decay into γ.
Quantitatively, N eff with a mirror photon is where ρ ν , ρ γ , and ρ γ , are the relic energy densities of neutrinos, mirror photons, and photons, respectively.According to Eq. (4.4), an additional mirror photon in thermal equilibrium significantly increases N eff and is generally excluded by current ∆N eff limits [19].Nevertheless, a mirror photon can be allowed if it does not achieve a thermal abundance.
This leads us to consider the Boltzmann equation for f a (p) in the freeze-in picture where the mirror photon collision term, C γ , is given by Eq. (4.5) replaces Eq. (3.8) when including mirror photons.In addition, while the energy density evolution of the Standard Model bath remains as given in Eqns.(3.17) and (3.19), 3the energy density evolution of γ is described by Conservatively, we take the initial γ density at T χPT to be zero.Due to the substantial change in Standard Model degrees of freedom across T QCD , freeze-in production of γ much earlier than T χPT is diluted and this conservative estimate is a fairly good approximation to the true initial abundance of γ .Last, the Hubble expansion rate, (3.21), is modified to include the additional mirror photon energy density, Figures 11 and 12 show the contours of N eff when including a massless mirror photon in heavy axion cosmologies for E/N = 8/3 and 1/3, respectively.As before, the blue region indicates where N eff is excluded by current CMB measurements at the 2σ level.Note that whereas ∆N eff is strictly negative in the case without the mirror photon (Figs. 6 and 7), the case with the mirror photon gives positive ∆N eff for most of the parameter space where the mirror photon dominates the branching ratio.
The green region indicates where the mirror photon reaches equilibrium and the freeze-in picture breaks down.This occurs when Γ aγ γ H(T = m a ), or equivalently, roughly when . (Mirror Photon reaches thermal equil.)(4.9) Within this green region, the contour values for N eff are artificially high because the mirror photon acquires a greater than thermal abundance due to the lack of a back reaction in  Eq. (4.5).In principle, capping the mirror photon abundance at a thermal abundance suggests that realistic N eff contours within the green region are roughly fixed at the value of N eff on the boundary of the green region; that is, the value of N eff when the mirror photon just acquires a thermal abundance from the freeze-in picture.For axions decaying prior to neutrino decoupling, this argument suggests N eff 3.8 in the green region when E/N = 8/3.Such a large ∆N eff is already excluded by experiments and thus the freeze-in picture is generally valid within the experimentally allowed region.However, for axions decaying after neutrino decoupling, it is possible that a tuned cancellation between the ρ γ energy deposit (positive ∆N eff contribution) and the heating of ρ γ relative to neutrinos (negative ∆N eff contribution) can occur in the green region.We leave the calculation of such a tuned cancellation to future work, but we expect that the parameter region with τ 1 sec is excluded by BBN.This is because to cancel the positive ∆N eff , the axion decays before the proton-neutron conversion completes, and the Helium abundance will be affected.

Light Axion and Dark Matter Over-Production
In this section, we discuss constraints on the heavy QCD axion for m a < 1 MeV.As its mass is decreased the axion remains excluded by ∆N eff until it decays after the CMB era.However, at this point, the axions remains excluded from its contribution to dark matter at the CMB era, until a significant further reduction in its mass.Before entering the allowed light axion region, there is an excluded region from free-streaming effects on large scale structure.
As can be seen from Fig. 6 and 7, the CMB limit on ∆N eff excludes 1 MeV m a 3 MeV for any f a .This exclusion from ∆N eff continues for m a < 1 MeV, until the axions decay after recombination. 4 Hence, the blue region of Fig. 13, where τ < t CMB ≡ 370, 000 yrs (z ≈ 1100) [51], is excluded by ∆N eff .For τ t CMB and m a > 0.13(Ω DM h 2 /0.12)Y a (T FO,g ) −1 eV, the axion is sufficiently heavy and long-lived to exceed the observed dark matter density at the CMB era, as shown by the excluded orange region of Fig. 13.Here, we assume the reheat temperature is sufficiently high that axions undergo freeze-out, as shown in Fig. 1, giving an axion freeze-out yield, Y a (T FO,g ), typically between 2 × 10 −3 to 2 × 10 −2 .Relaxing this assumption, by taking T RH below T FO,g or by introducing dilution between axion freeze-out and BBN, reduces the orange excluded region.
At lower axion masses, the free-streaming of axions suppresses the matter spectrum (i.e., the axion is hot dark matter) and for even smaller masses, the axion works as dark radiation.We reinterpret the bound derived in [71] for our framework and exclude the green-shaded region in Fig. 13.Here we conservatively impose the bound only for τ > t 0 , where t 0 is the present age of the Universe, but we expect that the bound is also applicable as long as τ > t eq , since the suppression of the matter spectrum is dominated by the free-streaming of axions before the matter-radiation equality.A part of orange, blue, or green region is also excluded by other astrophysical constraints (see [72] for an overview), but they are generically weaker.
Figure 13.Cosmological constraints on the low mass axion region, m a < 1 MeV.The axion decays prior to recombination (t CMB ≈ 370, 000 yr) in the orange region, which is excluded by the N eff constraint.The boundary denoting τ = t CMB is given for E/N = 0 (solid) and 8/3 (dashed) without a mirror photon; and for 1/3 (dotted) and 8/3 (dot-dashed) with a mirror photon decay channel.In the blue region, the axion decays after recombination, but possesses a matter energy density during the CMB era exceeding that measured by Planck [19], assuming the axion froze-out in the early Universe.For f a 10 8 GeV, T FO occurs below the electroweak scale, and hence the axion energy density at t CMB grows because of the reduction in g * .In the green region, the axion suppresses the structure formation as hot dark matter or contributes to dark radiation.

Conclusions
The strong CP problem can be addressed in a wide variety of axion models.The minimal ones, where the QCD axion mass is solely given by strong QCD dynamics, predict m a f a ∼ (100 MeV) 2 , but are typically plagued by a quality problem.This quality problem can be ameliorated or solved in a range of "Heavy QCD Axion" theories, where m a f a is orders of magnitude larger than in the minimal models.The constraints and search strategies for these heavy axions are completely different from those for the conventional lighter axion.An important constraint from CMB data arises if the axion lifetime is in the range of 10 −1 s -10 −12 s, decaying after neutrino decoupling at the MeV era, but before last scattering of the CMB at the eV era.In this case, the energy density of neutrinos is diluted, affecting the dark radiation at the CMB era, N eff , which has been precisely measured by the Planck Collaboration [19] and will be significantly improved by CMB Stage 4 experiments [20].Thus, theory and experiment both strongly motivate a detailed study of this cosmological bound on heavy axions.For axion masses above 1 MeV, except for accelerator searches at low values of f a , N eff is the strongest bound on the heavy QCD axion, and is the focus of this work.
A well-motivated and predictive model involves a mirror copy of the SM with a large axion mass generated by the mirror QCD interaction.In this case there is a competition between axion decays to photons diluting the neutrino contribution to N eff and axion decays to mirror photons directly enhancing N eff .We have also provided a detailed analysis of the N eff bound in theories with a light mirror photon.
Our analysis of the N eff bound takes into account key pieces missing from previous studies of the heavy QCD axion by developing a Boltzmann code that follows the evolution of the momentum distribution for the axion.The mesons and gluons of QCD play a key role; we include axion-pion scattering, axion decay to final states involving mesons, and axionmeson mixing.In addition we follow a detailed cosmological evolution from the initial axion abundance from freeze-out to the non-trivial QCD contributions in the Friedmann equations.
Our results for the CMB N eff constraints on the heavy QCD axion, in the absence of a mirror photon, are shown in Figs. 6 and 7, and are very powerful.Planck excludes large areas of parameter space, especially at large f a , but large areas remain at low f a , where the quality problem is solved for operators of dimension 6 and larger.The discovery reach of CMB-S4 at larger values of m a is modest, but improves at lower m a : for example, if f a is of order 10 4 GeV, CMB-S4 will see a signal for m a in the range of (3-10) MeV.
We find two important differences from standard results, illustrated by comparing Fig. 5 with our results shown in Figs. 6 and 7. First, resonances occur when the axion mass is around the π 0 , η, and η masses, greatly affecting N eff for axion masses between (100 − 1000) MeV.Second, by including mesons and gluons, we correctly take account of the axion lifetime.This is a large effect, especially at large m a , increasing the decay rate by orders of magnitude as m a rises above ∼ 1 GeV; this point is apparent in the right panels of Figs. 6 and 7 where regions with higher f a open up.
In the presence of a light mirror photon, our results for the CMB N eff constraints are shown in Figs.11 and 12 for E/N = 8/3 and 1/3 respectively.For E/N = 8/3, m a < 100 MeV is excluded for all values of f a .A substantial fraction of the allowed region with 100 MeV < m a < 500 MeV will be probed by CMB-S4 via a positive signal for ∆N eff .For E/N = 1/3, the CMB N eff bound is considerably weaker.A new allowed region opens up at lower axion masses, 1 MeV < m a < 100 MeV, where dark radiation from the mirror photon compensates neutrino dilution from axion decays.A large fraction of this region gives a CMB-S4 signal, with ∆N eff positive (negative) for smaller (larger) values of f a .These allowed regions both solve the quality problem for operators of dimension 6 and larger.
The current and future 95% confidence limit on the axion mass from N eff in this work are shown in comparison to other cosmological and astrophysical constraints in Fig. 14.The limits on N eff in this work provide the strongest constraints on heavy QCD axions for m a 1 MeV and f a 10 5 GeV.Complementary constraints at small f a arise from direct heavy axion searches at accelerators [45,68] as shown by the purple shaded regions.Producing axions in a beam dump, such as the DUNE Near Detector, and discovering their subsequent decays, will  6).The dashed orange contour shows the future reach of N eff from CMB-S4.The left boundary of the orange region indicates where the axion decays after recombination which marks the parameter space where N eff contstraints become inapplicable (see Fig. 13).The blue region indicates where the axion energy density at recombination is greater than that of dark matter, assuming the reheat temperature of the universe is high enough that the axion froze-out with a thermal abundance.The solid and dashed purple regions shows the present and future bounds from accelerator searches, respectively [45,50,[60][61][62][63][64][65][66][67][68].The brown region shows the bound from supernova 1987a on axions with hadronic interactions, [76], green from horizontal-branch cooling [77], red from CAST [78,79], and pink from solar neutrinos [80].Below the dot-dashed contour, the axion energy density arising from the misalignment mechanism, with a misalignment angle of unity, is greater than the observed dark matter energy density.In the green region, the axion suppresses the structure formation or contributes to N eff .Below the dotted contours, the axions suffers a PQ quality problem arising from operators of the labeled dimension.The diagonal yellow strip indicates the standard QCD axion, which highlights the severity of the PQ quality problem for low mass axions.allow the region enclosed by the dashed purple contour to be probed [45].Furthermore, the dashed purple contour at higher m a and low f a can be probed by observing axions in B meson decays at Belle [68].The bound from axion cooling of Supernova 1987A has uncertainties arising from the temperature and density profiles of the supernova, and has been computed for a variety of such profiles in [76]; we show a conservative case.Constraints on the decay of the axion from extragalactic background light or CMB spectral distortions are derived in [18], but the constraints do not exclude the parameter region that is allowed in Fig. 14.
The bounds shown in Fig. 14, and elsewhere in the paper, are computed assuming that the reheat temperature of the universe T RH is above the axion freeze-out temperature, and that there is no subsequent dilution of the axion abundance, for example from late decaying particles.Removing this assumption relaxes the bounds, since axion production occurs via freeze-in rather than freeze-out, or is diluted after freeze-out.Since the freeze-out temperature decreases as f a drops, relaxing the bounds becomes harder at lower f a .For f a < 10 4 GeV, the axion is kept into thermal equilibrium even at T < 4 MeV, and the BBN bound T RH > 4 MeV [81][82][83] excludes the possibility of relaxing the bound.It may be plausible that the reheating temperature is below the freeze-out temperature for large f a , but solving the quality problem favors low f a , and it is typically harder to obtain a large enough m a for large f a ; see e.g., Eq. (1.3).
The next decade will yield exciting and important answers to axion physics.Heavy QCD axions provide a highly-motivated solution to the strong CP problem.Unlike the standard QCD axion, which induces a small ∆N eff signal [26], heavy QCD axions can generate substantial ∆N eff signals that can be probed by the exquisite sensitivity of current and near future CMB telescopes.Moreover, in theories without a mirror photon, this signal results from a depletion of the cosmic neutrino abundance, providing a less common fingerprint of a negative contribution to ∆N eff .Such a measurement would determine a correlation between the axion mass and decay constant.As a result, the total squared amplitude for scatterings involving the axion and charged pions is where t ≡ (p a − p 2 ) 2 , u ≡ (p a − p 3 ) 2 and π i , π j , π k in the subscript of (A.2) refer to three pions of different charge π 0 , π + , π − .Likewise, the interaction between an axion and three neutral pions generates the squared scattering amplitude Inserting the sum of the squared scattering matrix elements, (A.2) and (A.3), into Eq.(3.9) and integrating over the phase spaces of the three pions gives the product of the axion-pion collision term and f a,eq − f a and f i is the distribution function of particle i possessing momentum p i in accordance with Fig. 15.In going from (A.5) to (A.6), we take the pions to be in thermal equilibrium with the Standard Model thermal bath so that f 1 , f 2 , and f 3 follow a Bose-Einstein distribution of temperature T .The temperature T of the strongly coupled thermal bath is inferred at each numerical time step by solving the following equation for T : i=γ,e,µ,π,δ QCD where the left-hand side of (A.7) is solved from Eq. (3.17) and the right-hand side is calculated from Eq. 3.18 for γ, µ, e, π and from [58] for ρ δ QCD , as described more in Appendix C.
To calculate C π , we first introduce another δ-function in (A.4) by writing dΠ 3 (2π) 3 = ).By integrating p 3 over the other delta function δ 4 (p 1 + p 2 − p 3 − p 4 ), Eq. (A.4) simplifies to with the understanding that p 3 = p a + p 1 − p 2 .Note that the argument of the remaining delta-function, p 2 3 − m 2 3 , can be written as where, in the notation of [84], , and α, θ, γ are the angles between p a and p 1 , p a and p 2 , and p 1 and p 2 , respectively.It is convenient to express the latter angle in terms of the former two by cos γ = cos α cos θ + sin α sin θ cos β.
In the massless axion limit, the argument of the remaining delta function can easily be expressed in terms of E a as done in [53].However, in the massive axion limit, this is impossible and it is thus more useful to move the argument of the delta function onto one of the scattering angles, as done in [84], which we follow.In particular, the choice of the angle β is most convenient as it only occurs once in (A.9).The resulting integral for C π is where g(β) = p 3 (β) 2 − m 2 3 .The integrals over the azimuthal angles β and φ can be done analytically due to the delta function and the lack of φ dependence in the integrand.The remaining integrals over the pion 3-momenta |p 1 | ∈ [0, ∞) and |p 2 | ∈ [0, ∞) and the polar angles cos α ∈ [−1, 1] and cos θ ∈ [−1, 1] are done numerically using Monte Carlo integration.To ensure the integration region is performed only in the kinematically allowed region, we include a Heaviside function Θ(1−cos 2 β i ) in (A.10), where β i is the location of the two (equal and opposite) roots of g(β) [84].Note that β i are functions of the other four integration variables.
Figure 16.Numerical evaluation of F π as a function of m π /T for a variety of axion masses.F π becomes suppressed for m a > m π and T < m π since the production of axions must occur at the Boltzmann tail of the incoming pions.The massless axion limit, which has been computed previously in literature [53,85], corresponds to the m a /m π = 0 contour.
We can gain intuition for the cosmological effect of axion-pion scattering by calculating the thermally averaged pion-to-axion scattering rate as introduced in Eq. (3.2), where n a,eq is the thermal number density of axions with mass m a at temperature T .As before, A ≡ 1 3 (1 − z)/(1 + z), z ≡ m u /m d , and r ≡ m a /m π .Fig. 16 shows F π as a function of m π /T for a variety of axion masses.For m a m π , F π reduces to previous results in the literature for massless axions [53,57], with F π related to the h LO function defined in [53,57] by the mapping F π (m π /T, m a = 0) ≡ 0.212 h LO (m π /T ).According to Fig. 16, for T < m π , F π drops as m a increases.This is because in this regime, only the Boltzmann tail of pions with high energies m a can kinematically scatter to produce axions.

B N eff for Large Axion Masses
In the main text, we show the N eff constraint for m a < 2 GeV.In this appendix, we derive the constraint for m a > 2 GeV.
In the left panels of Fig. 6, 7, 11 and 12, the N eff contours become approximately horizontal above m a ≈ 1 GeV, signifying that N eff is dictated mainly by the lifetime of the axion in this region.However, careful inspection indicates that the slope of the contours is not quite flat, but slightly decreases as m a grows.The reason is, for fixed τ , the energy density of the axion at decay increases with increasing mass.This follows because axions in this region have such large f a that they decouple early and decay non-relativistically. Thus, what actually sets the N eff contours in the m a > 2 GeV region is how much energy density they deposit into the thermal bath right at neutrino decoupling.
For example, let ρ max be the maximum energy density that can be deposited at a certain time t * so that N eff does not drop below an arbitrary contour, N eff,0 , which we will take to be 2.62, the 2σ limit on N eff allowed by Planck.Choose a point (m 0 , τ 0 ) that lies on this N eff,0 contour in the m a 1 GeV region.Analytically, the energy density of this non-relativistic axion at time t We perform a numerical fit of the N eff,0 = 2.62 contour with the anchor point (m 0 , τ 0 ) = (2.0GeV, 4.3 × 10 −2 s) and find t * 0.17 s.Other anchor points give similar t * .Eq. (3.22) follows from this fit.

C Numerical Approaches
In this section, we discuss the numerical techniques used to solve the Boltzmann equation describing the cosmological evolution of the axion.As mentioned in Sec.3.2, we employ the method of lines technique to convert the Boltzmann system of partial differential equations into a system of ordinary differential equations.In particular, we discretize the partial differential equation governing the axion phase space density, (3.8), into a partition of N ordinary differential equations, {f p,i (t)} of time, with each ODE corresponding to the time evolution of the phase space density at a fixed comoving momentum, |p i | = |p i |a(t), with i ∈ {1, ..., N }.
In our numerical setup, we split (3.8) into N = 24 ODEs of logarithmically equidistant |p i | where i = 1 corresponds to the fixed comoving momentum | p1 | = e −8 T χPT and i = N corresponds to |p N | = e 7/2 T χPT .Note that because of entropy conservation -which holds except for when the axion dominates the energy density of the universe before decaying -each fixed comoving momentum equals the ratio of the physical momentum to the temperature, |p|/T .We have verified the convergence of our results by checking that N eff changes by less than 1% when using larger N .In addition, our results for N eff for just the Standard Model cosmology, or equivalently, when the axion decays far before neutrino decoupling, is 3.040.This value slightly differs from the true Standard Model value of N SM eff = 3.044 [86][87][88] by ∼ 0.2% because we do not include effects from QED corrections or neutrino oscillations.
The dynamical timescale (∼ 1/H, where H is Hubble) involved in the cosmological evolution of the Boltzmann equation spans many orders of magnitude from the end of the QCD phase transition to past neutrino decoupling.Hence, we solve the system of Boltzmann equations in terms of the logarithmic timescale y = ln(t/t χPT ), where t χPT = 1/2H(T χPT ) is the starting time of the Boltzmann code.
Last, we determine the extra QCD contributions to the energy density (ρ δ QCD ) and pressure (P δ QCD ) arising from heavy mesons and from ideal gas law deviations using the calculations of [58], which tabulated g * (T ) and g * S (T ) for the Standard Model across the QCD phase transition.The top left panel of Fig. 17 shows g * as a function of temperature T .The orange contour shows the value of g * assuming an ideal gas comprised of photons, neutrinos, electrons, muons, and pions in thermal equilibrium.The blue contour, which diverges from the orange above T ∼ 100 MeV, is the Standard Model value of g * from [58], which includes contributions, δ QCD , from heavy mesons and from ideal gas law deviations.We extract ρ δ QCD by taking the difference between the blue and orange contours, ∆g * , and multiplying this result by π 2 T 4 /30, as shown by the top right panel.
Similarly, the bottom left panel of Fig. 17 shows the ratio P/T 4 = g * P as a function of temperature T , where P is the total pressure of the relevant species.As before, the orange contour shows the value of g * P assuming an ideal gas comprised of photons, neutrinos, electrons, muons, and pions in thermal equilibrium.The blue contour, which diverges from the orange above T ∼ 100 MeV, is the Standard Model value of g * P from [58], which we infer through the relationship g * P = π 2 30 ( 4 3 g * S − g * ).We extract P δ QCD by taking the difference between the blue and orange contours, ∆g * P , and multiplying this result by T 4 , as shown by the bottom right panel.

Figure 2 .
Figure 2. When T χPT < T FO < T gPT , the initial axion density at the start of the Boltzmann code is uncertain.The left panel bounds this uncertainty by showing the two extreme T FO ∈ (T χPT , T gPT ) that give the largest (blue) and smallest (orange) axion densities at the starting temperature of the Boltzmann code.The right panel shows the corresponding axion yields.The cross-over between the maximum and minimum T FO around 0.3 GeV arises from the balance between Boltzmann suppression and dilution by g * across the QCD phase transition.

Figure 3 .
Figure 3. (Left) Axion freeze-out temperature from pion (solid) and Primakoff (dashed) scatterings as a function of f a for a variety of m a between 0 and 200 MeV.For f a 3 × 10 4 GeV, the axion-pion scattering dominates over Primakoff scattering, keeping the axions in thermal equilbrium until T ∼ 10 MeV for f a 10 6 GeV.

Figure 4 .
Figure 4. (Right) Total axion decay rate (τ −1 ) as a function of m a for E/N = 0 (blue) and E/N = 8/3 (orange).Resonant peaks in the decay rate arise from axion-meson mixing when m a is near m π 0 , m η , and m η .Troughs arise from cancellations in g γ between the anomaly, axial rotation, and mesonmixing contributions.For m a 2 GeV, the decay rate is set by gluons.The dashed contours show the axion-photon decay rate without meson mixing.

Figure 5 .
Figure5.Contours of N eff in the m a − τ γ (left) and m a − f a (right) planes neglecting the following crucial effects that we take into account in future figures: axion-pion scattering, axion-meson decays, axion-meson mixing in g γ , the proper frozen-out initial axion abundance, and the QCD contributions to the background evolution.Note that by neglecting these effects, as done in past works, the axionlifetime is incorrectly set by the axion-photon decay rate with no axion-meson mixing, as shown by the y-axis of the left panel.The dark orange region is excluded at 95% confidence by Planck, while the light orange shows the future reach of CMB-S4 experiment at 95% confidence.

Figure 6 .
Figure 6.Contours of N eff in the m a − τ (left) and m a − f a (right) planes for E/N = 0.The dark orange region is excluded at 95% confidence by Planck, while the light orange shows the future reach of CMB-S4 experiment at 95% confidence.The dashed contours indicate where the initial axion yield is uncertain because T FO lies between T χPT and T gPT .The dashed contours bound this uncertainty by showing the value of N eff taking the minimum initial axion yield while the solid contours show the value of N eff taking the maximum initial axion yield as given in Fig. 2.

. 22 )
Eq.(3.22)  can also be written in terms of f a by equating τ −1 max with the analytic decay rate into gluons given by Eqns.(3.3) and (3.4), f a 1.2 × 10 10 GeV m a

Figure 10 .
Figure 10.As in Fig.9, but for E/N = 1/3.In the left panel, the solid and dashed black contours, giving the axion decay width with and without mirror photons, are almost conincident, as the partial width to mirror photons is sub-dominant, as shown by the straight blue line.

Figure 11 .
Figure 11.Contours of N eff in the m a − τ (left) and m a − f a (right) planes for E/N = 8/3 when including a mirror photon.The dashed contours indicate where the initial axion yield is uncertain because T FO lies between T χPT and T gPT .The dashed contours bound this uncertainty by showing the value of N eff taking the minimum initial axion yield while the solid contours show the value of N eff taking the maximum initial axion yield as given in Fig.2.The green region shows where f a is sufficiently small that the dark photon reaches a thermal abundance and the freeze-in picture we employ breaks down.The dark orange region is excluded at 95% confidence by Planck, while the light orange shows the future reach of CMB-S4 experiment at 95% confidence.

Figure 14 .
Figure 14.Overview of the excluded parameter space as determined in this work (blue region) in relationship to other QCD axion bounds.The solid orange region shows the 2σ exclusion region on N eff as determined by Planck for the heavy QCD axion with E/N = 0 (see Fig.6).The dashed orange contour shows the future reach of N eff from CMB-S4.The left boundary of the orange region indicates where the axion decays after recombination which marks the parameter space where N eff contstraints become inapplicable (see Fig.13).The blue region indicates where the axion energy density at recombination is greater than that of dark matter, assuming the reheat temperature of the universe is high enough that the axion froze-out with a thermal abundance.The solid and dashed purple regions shows the present and future bounds from accelerator searches, respectively[45,50,[60][61][62][63][64][65][66][67][68].The brown region shows the bound from supernova 1987a on axions with hadronic interactions,[76], green from horizontal-branch cooling[77], red from CAST[78,79], and pink from solar neutrinos[80].Below the dot-dashed contour, the axion energy density arising from the misalignment mechanism, with a misalignment angle of unity, is greater than the observed dark matter energy density.In the green region, the axion suppresses the structure formation or contributes to N eff .Below the dotted contours, the axions suffers a PQ quality problem arising from operators of the labeled dimension.The diagonal yellow strip indicates the standard QCD axion, which highlights the severity of the PQ quality problem for low mass axions.

Figure 17 .
Figure 17.Top left: The orange contour shows g * as a function of temperature assuming an ideal gas of photons, neutrinos, electrons, muons, and pions, while the blue contour shows g * for the Standard Model including heavy mesons and deviations from the ideal gas law[58].Top right: ρ δQCD as a function of temperature, extracted from the difference in the blue and orange g * contours on the top left panel.Bottom panels: Same as top panels but for g * P ≡ P/T 4 .