11institutetext: Department of Physics, University of California, Berkeley,
California 94720, USA22institutetext: Theoretical Physics Group, Lawrence
Berkeley National Laboratory, Berkeley, California 94720, USA33institutetext:
Theoretical Physics Department, CERN, Geneva, Switzerland

#  Dark Radiation Constraints on Heavy QCD Axions

David I. Dunsky 1,2 Lawrence J. Hall 3 Keisuke Harigaya

###### Abstract

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.

##  1 Introduction

The smallest dimensionless parameter of the Standard Model is the strong CP
parameter $\bar{\theta}\lesssim 10^{-10}$. This small parameter can be
understood as resulting from a discrete spacetime symmetry, $CP$ Nelson:1983zb
; Barr:1984qx  or $P$ Beg:1978mt ; Mohapatra:1978fy ; Babu:1989rb , or from a
global Abelian Peccei-Quinn (PQ) symmetry Peccei:1977hh ; Peccei:1977ur .
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 Georgi:1981pu . In such scenarios, the PQ symmetry is expected to
be explicitly broken by higher-dimensional operators, typically preventing
sufficient dynamical relaxation of $\bar{\theta}$ towards zero even if they
are suppressed by powers of the Planck mass, $M_{\rm Pl}$.

For example, if $\phi$ is the field that spontaneously breaks the PQ symmetry
at scale $f_{a}$, in standard axion theories the interactions
$\lambda_{n}\phi^{n+4}/M_{\rm Pl}^{n}$ must be suppressed to solve the strong
CP problem,

$\displaystyle\lambda_{n}\;<\;10^{(-46+10n)}\left(\frac{10^{8}\,{\rm
GeV}}{f_{a}}\right)^{4+n}.$ (1)

Even for the lowest values of $f_{a}$ allowed by observations, of order
$10^{8}\,{\rm 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 Harlow:2018tng ; Banks:2010zn , reintroducing the quality
problem Holman:1992us ; Barr:1992qq ; Kamionkowski:1992mf ; Dine:1992vx .
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 Choi:1992xp ;
Dine:1992ya  making them attractive avenues for the strong CP problem.

The severity of this PQ quality problem, shown in (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,111In models where the spontaneous PQ breaking occurs by hidden quark
condensation Choi:1985cb , this only requires the removal of hidden quark mass
terms. which may underlie the accidental PQ symmetry, the shift of
$\bar{\theta}$ from zero by a dimension 6 operator is sufficiently small if

$\displaystyle m_{a}\;\mathop{}_{\textstyle\sim}^{\textstyle>}\;{\rm
MeV}\left(\frac{f_{a}}{3\times 10^{4}\,{\rm GeV}}\right)^{2}.$ (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 Cadamuro:2011fd . There is a powerful bound on
this dark radiation from measurements of the Cosmic Microwave Background
radiation (CMB) by the Planck Collaboration Planck:2018vyg , $N_{\rm
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 CMB-S4:2016ple . 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. Chang:1993gm ; Ferreira:2018vjj ; Arias-Aragon:2020qtn ; Arias-
Aragon:2020shv ; Ferreira:2020bpb ; DEramo:2021psx ; DEramo:2021lgb ;
DEramo:2021usm  for studies on axions that are light and stable and directly
contribute to dark radiation, and Ref. Hannestad:2003ye ; Giare:2020vzo ;
Caloni:2022uya  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
Dimopoulos:1979pp ; Tye:1981zy ; Holdom:1982ex ; Holdom:1985vx ; Flynn:1987rs
; Agrawal:2017ksf . It is important that these instantons do not probe new CP
violating phases, so that the new contribution to the potential aligns
$\bar{\theta}$ to be sufficiently small. The growth in the QCD coupling in the
UV could also arise from extra spatial dimensions Gherghetta:2020keg .

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^{\prime}\gg 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, $\Lambda^{\prime}\gg\Lambda$.
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
$(\Lambda^{\prime}/\Lambda)^{2}$. The first implementations of this idea
Rubakov:1997vp ; Berezhiani:2000gh  used a Weinberg-Wilczek axion
Weinberg:1977ma ; Wilczek:1977pj , with the PQ symmetry spontaneously broken
by Higgs vevs. In this case the axion decay constant $f_{a}$ is large, of
order $v^{\prime}$, and while these theories ameliorate the quality problem of
(1, 2), solving the problem requires contrived arrangements. On the other
hand, in the theories considered in Fukuda:2015ana ; Hook:2019qoh ;
Kelly:2020dda  using a KSVZ axion Kim:1979if ; Shifman:1979if , the heavy QCD
axion mass is

$\displaystyle m_{a}\;\sim\;100\,{\rm MeV}\left(\frac{v^{\prime}}{10^{8}\,{\rm
GeV}}\right)^{8/11}\left(\frac{10^{4}\,{\rm GeV}}{f_{a}}\right)$ (3)

so that the quality problem is solved by taking $f_{a}\ll v^{\prime}$. As
always, one still needs to understand PQ in operators of dimension $\leq 4$ as
an accidental consequence of gauge symmetries Georgi:1981pu . 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. Cadamuro:2011fd ; Millea:2015qra ;
Depta:2020wmr  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.
Fukuda:2015ana  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_{\rm 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_{\rm
eff}$ constraint is not applicable. We conclude the paper with Sec. 6.

##  2 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

$\displaystyle{\cal
L}=\frac{1}{2}\partial^{\mu}a\partial_{\mu}a-\frac{1}{2}m_{a}^{2}a^{2}+\frac{g_{3}^{2}}{32\pi^{2}}\frac{a}{f_{a}}G_{\mu\nu}\tilde{G}^{\mu\nu}+\frac{e^{2}}{32\pi^{2}}\frac{E}{N}\frac{a}{f_{a}}F_{\mu\nu}\tilde{F}^{\mu\nu}.$
(4)

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_{\rm 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

$\displaystyle{\cal L}_{\rm chiral}$ $\displaystyle=\frac{f_{\pi}^{2}}{4}{\rm
Tr}\\{D^{\mu}U^{\dagger}D_{\mu}U\\}+\frac{f_{\pi}^{2}}{4}{\rm
Tr}\\{2BM_{q}U+\text{h.c.}\\}+\frac{1}{2}\frac{\partial^{\mu}a}{f_{a}}{\rm
Tr}\\{Q_{A}\lambda^{a}\\}J{{}_{\mu}^{a}}$
$\displaystyle+\frac{e^{2}}{32\pi^{2}}(\frac{E}{N}-\frac{8}{3}Q_{u}-\frac{2}{3}Q_{d}-\frac{2}{3}Q_{s})\frac{a}{f_{a}}F_{\mu\nu}\tilde{F}^{\mu\nu},$
(5)

where $f_{\pi}=93$ MeV, $B$ is a strong interaction parameter of order the QCD
scale, $J_{\mu}^{a}$ is the $SU(3)$ axial current given by

$\displaystyle J{{}_{\mu}^{a}}$ $\displaystyle=\frac{i}{4}f_{\pi}^{2}{\rm
Tr}\\{\lambda^{a}(UD^{\mu}U^{\dagger}-U^{\dagger}D^{\mu}U)\\},$ (6)

and

$\displaystyle Q_{A}=\begin{pmatrix}Q_{u}&&\\\ &Q_{d}&\\\
&&Q_{s}\end{pmatrix}$ (7)

is the quark charge matrix for the transformation that eliminates the axion-
gluon coupling with ${\rm Tr}\;Q_{A}=1$.

The meson nonet and quark mass matrices are given by

$\displaystyle U$
$\displaystyle=\exp{i\frac{\sqrt{2}}{f_{\pi}}\begin{pmatrix}\frac{\pi_{0}}{\sqrt{2}}+\frac{\eta_{8}}{\sqrt{6}}+\frac{\eta_{0}}{\sqrt{3}}&\pi^{+}&K^{+}\\\
\pi^{-}&-\frac{\pi_{0}}{\sqrt{2}}+\frac{\eta_{8}}{\sqrt{6}}+\frac{\eta_{0}}{\sqrt{3}}&K^{0}\\\
K^{-}&\bar{K}^{0}&-2\frac{\eta_{8}}{\sqrt{6}}+\frac{\eta_{0}}{\sqrt{3}}\end{pmatrix}},$
(8) $\displaystyle M_{q}$
$\displaystyle=\begin{pmatrix}m_{u}e^{i\frac{a}{f_{a}}Q_{u}}&&\\\
&m_{d}e^{i\frac{a}{f_{a}}Q_{d}}&\\\
&&m_{s}e^{i\frac{a}{f_{a}}Q_{s}}\end{pmatrix}.$ (9)

When calculating the effect of axion-meson scattering at energies below
$\Lambda_{\rm QCD}\approx 150$ MeV, we limit ourselves to the relevant two-
dimensional subspace, $SU(2)_{L}\times SU(2)_{R}\to SU(2)_{V}$, since the only
active QCD degrees of freedom are pions. For axion masses above $\Lambda_{\rm
QCD}$, the chiral perturbation based on the symmetry $SU(2)_{L}\times
SU(2)_{R}\to SU(2)_{V}$ breaks down. In this case, the axion may decay to
heavy mesons like $\eta$ and $K$, or for sufficiently large $m_{a}$, directly
to gluons. Consequently, we use the results of Aloni:2018vki  for the axion
decay rate into mesons, gluons, and photons (including the enhancement in the
axion-photon coupling from $a-\eta(\eta^{\prime})$ mixing) for
$m_{a}>m_{\pi}$. 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) generates the following axion-pion mass matrix
and interactions

$\displaystyle{\cal L}_{a,\pi}$
$\displaystyle=-\frac{1}{2}\begin{pmatrix}\pi_{0}&&a\end{pmatrix}\begin{pmatrix}B(m_{u}+m_{d})&&B(f_{\pi}/f_{a})\;(Q_{u}m_{u}-Q_{d}m_{d})\\\
B(f_{\pi}/f_{a})\;(Q_{u}m_{u}-Q_{d}m_{d})&&m_{a}^{2}+B(f_{\pi}/f_{a})^{2}\;(Q_{u}^{2}m_{u}+Q_{d}^{2}m_{d})\end{pmatrix}\begin{pmatrix}\pi_{0}\\\
a\end{pmatrix}$
$\displaystyle+\frac{B(Q_{u}m_{u}+Q_{d}m_{d})}{24f_{\pi}^{2}}(\pi_{0}^{4}+4\pi_{-}\pi_{+}\pi_{0}^{2})+\frac{B(Q_{u}m_{u}-Q_{d}m_{d})}{6f_{\pi}}\frac{a}{f_{a}}(\pi_{0}^{3}+2\pi_{0}\pi_{+}\pi_{-})$
$\displaystyle+\frac{1}{3f_{\pi}^{2}}(\pi_{+}\pi_{0}\,\partial^{\mu}\pi_{-}\partial_{\mu}\pi_{0}+\pi_{-}\pi_{0}\,\partial^{\mu}\pi_{+}\partial_{\mu}\pi_{0}-\pi_{0}^{2}\,\partial^{\mu}\pi_{-}\partial_{\mu}\pi_{+}-\pi_{-}\pi_{+}\,\partial^{\mu}\pi_{0}\partial_{\mu}\pi_{0}).$
(10)

From the $\pi_{0}-a$ mass matrix in the limit $f_{a}\gg f_{\pi}$, we can
identify $B=m_{\pi}^{2}/(m_{u}+m_{d})$. Moreover, since we study axions
heavier than the standard QCD axion, $m_{a}f_{a}\gg m_{\pi}f_{\pi}$, and we
can thus drop the second term in the bottom right entry of the mass matrix.

The $\pi_{0}-a$ mass matrix is diagonalized by the rotation

$\displaystyle\begin{pmatrix}\pi_{0}\\\
a\end{pmatrix}=\begin{pmatrix}\cos\theta&&\sin\theta\\\
-\sin\theta&&\cos\theta\\\ \end{pmatrix}\begin{pmatrix}\hat{\pi}_{0}\\\
\hat{a}\end{pmatrix}$ (11)

with

$\displaystyle\tan
2\theta=\frac{f_{\pi}}{f_{a}}\frac{1-z}{1+z}\frac{1}{1-r^{2}},$ (12)

where $r\equiv m_{a}/m_{\pi}$ and $z\equiv m_{u}/m_{d}\simeq 0.47$
ParticleDataGroup:2020ssz . In terms of the mass eigenstates
$\hat{\pi}_{0},\hat{a}$, the interaction of the axion with three pions is
described by

$\displaystyle{\cal L}_{\hat{a},\hat{\pi}}$
$\displaystyle\supset\frac{A}{f_{a}f_{\pi}}\frac{1}{1-r^{2}}\left[\partial^{\mu}\hat{a}(\pi_{-}\hat{\pi}_{0}\,\partial_{\mu}\pi_{+}+\pi_{+}\hat{\pi}_{0}\,\partial_{\mu}\pi_{-}-2\pi_{+}\pi_{-}\,\partial_{\mu}\hat{\pi}_{0})+\frac{r^{2}}{4}m_{\pi}^{2}\hat{a}(2\hat{\pi}_{0}\pi_{-}\pi_{+}+\hat{\pi}_{0}^{3})\right]$
(13)

where $A=\frac{1}{3}(1-z)/(1+z)\simeq 0.12$. For the remainder of this paper,
we drop the hats and refer to the mass eigenstates as $a$ and $\pi_{0}$. Note
that we do not consider the case where $m_{a}$ and $m_{\pi}$ are so highly
degenerate that the axion-pion mixing angle becomes of order unity since this
does not occur as long as

$\displaystyle\left|\frac{m_{\pi}-m_{a}}{m_{\pi}}\right|\gtrsim\frac{f_{\pi}}{f_{a}}\simeq
10^{-4}\left(\frac{10^{3}{\,\rm
GeV}}{f_{a}}\right)\qquad\Big{(}\parbox{132.30513pt}{\centering Axion-Pion \\\
Non-Degeneracy Condition\@add@centering}\Big{)}\,,$ (14)

which is only violated for axions that are extremely degenerate with pions.

The coupling of the axion with photons is

$\displaystyle{\cal L}_{a,\gamma}$
$\displaystyle=\frac{g_{\gamma}}{4}\frac{a}{f_{a}}F_{\mu\nu}\tilde{F}^{\mu\nu},$
(15) $\displaystyle
g_{\gamma}=\frac{e^{2}}{8\pi^{2}}\left(\frac{E}{N}-\frac{5}{3}-\mathcal{F}_{\theta}(m_{a})\right)$
(16)

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 $\mathcal{F}_{\theta}$ from axion-meson mixing. For
$m_{a}\ll m_{\eta}$, $\mathcal{F}_{\theta}$ reduces to $2\sin\theta
f_{a}/f_{\pi}$, where $\theta$ is the axion-pion mixing angle, (12). Further,
the value of $g_{\gamma}$ in the massless axion case is recovered in the limit
$m_{a}\ll m_{\pi}$ in which the term in parenthesis reduces to the standard
result, $\frac{E}{N}-\frac{2}{3}\frac{4+z}{1+z}\simeq\frac{E}{N}-2.03$
ParticleDataGroup:2020ssz . For $m_{a}>m_{\pi}$, we extract
$\mathcal{F}_{\theta}$ from the calculations of Aloni:2018vki  which include
$a-\eta$ and $a-\eta^{\prime}$ mixing. Note that previous considerations of
heavy-axion cosmological constraints, except for Kelly:2020dda , neglect the
effect of axion-meson mixing on $\mathcal{F}_{\theta}$ which leads to
significantly different values of $g_{\gamma}$ for $m_{a}>m_{\pi}$. 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$.

##  3 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}\gtrsim 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
$\Delta N_{\rm 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}(\mathbf{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 Salvio:2013iaa  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.

###  3.1 Axion Initial Conditions

At temperatures above the QCD phase transition temperature, $T_{\rm
QCD}\approx 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_{\rm 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_{\rm 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_{\chi\rm PT}\equiv 100$ MeV DiLuzio:2021vjd
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.

####  3.1.1 Equilibrium from Scatterings

When the temperature $T$ is much greater than $m_{a}$ and $T_{\rm QCD}$, the
axion-gluon interaction $a+g\leftrightarrow g+g$ dominates axion production.
The thermally averaged rate of axion-gluon scattering is given by
Salvio:2013iaa ; DEramo:2021lgb

$\displaystyle\Gamma_{ag\leftrightarrow
gg}\simeq\frac{16}{\pi}\left(\frac{g_{3}^{2}}{32\pi^{2}}\right)^{2}\frac{T^{3}}{f_{a}^{2}}\mathcal{F}_{g}(T)\,,$
(17)

where $\mathcal{F}_{g}(T)$ is a temperature-dependent function that captures
the axion production enhancement in the plasma from thermal gluon decays.
$\mathcal{F}_{g}(T)$ is numerically computed in Salvio:2013iaa ;
DEramo:2021lgb , and for $g_{3}\lesssim 1$, takes the approximate analytic
form $\mathcal{F}_{g}\approx 2g_{3}^{2}\ln 1.5/g_{3}$ Graf:2010tv ;
Salvio:2013iaa .

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_{\chi\rm PT}$. In
the yellow region, meson or gluon decays and inverse decays keep the axion in
equilibrium below $T_{\chi\rm 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_{\rm FO,g}$. The red
region indicates where the axions freeze out at temperatures above the
validity of chiral perturbation theory ($T_{\chi\rm PT}\equiv 100$ MeV) but
below the validity of gluon peturbation theory ($T_{\rm gPT}\equiv 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_{\chi{\rm PT}}$; see Eq. 23 for details.

On the other hand, when $T\lesssim T_{\rm QCD}$, the interaction
$a+\pi\leftrightarrow\pi+\pi$ dominates axion production. The thermally
averaged rate of axion-pion scattering is given by

$\displaystyle\Gamma_{a\pi\leftrightarrow\pi\pi}=\frac{T^{5}}{f_{a}^{2}f_{\pi}^{2}}\frac{A^{2}}{(1-r^{2})^{2}}\mathcal{F}_{\pi}(m_{a},T)$
(18)

where $\mathcal{F}_{\pi}(m_{a},T)$ is a temperature and axion mass dependent
function that we compute numerically in Appendix A. For reference, Fig. 16
shows $\mathcal{F}_{\pi}$ as a function of $m_{\pi}/T$ for a variety of
$m_{a}$.

We define the axion decoupling temperature, $T_{\rm FO}$, when
$\Gamma_{ag\leftrightarrow gg}=3H(T_{\rm FO,g})$ for $T\gg T_{\rm QCD}$, or
$\Gamma_{a\pi\leftrightarrow\pi\pi}=3H(T_{\rm FO,\pi})$ for $T\ll T_{\rm
QCD}$, where $H(T_{\rm FO})$ is the Hubble rate at $T_{\rm FO}$. Contours of
$T_{{\rm 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_{\rm FO,\pi}$ is less than
$T_{\chi\rm PT}$. The estimation of $T_{\rm FO,g}$ based on Eq. (17) breaks
down if $m_{a}>T_{\rm FO,g}$, but we find that it anyway occurs in the red
parameter region which possesses greater uncertainty: In the red-shaded
region, $T_{\rm FO}$ occurs above the temperature at which chiral perturbation
theory breaks down ($T_{\chi\rm PT}$) but below the temperature where the
strong coupling constant, $g_{3}$, becomes non-perturbative ($T_{g\rm PT}$).
We take $T_{\chi\rm PT}\simeq 100$ MeV, the temperature above which one-loop
corrections in chiral perturbation theory become comparable to tree-level
results DiLuzio:2021vjd . Similarly, $T_{g\rm PT}\simeq 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 Aloni:2018vki . In the red-shaded region, we cannot precisely
determine $T_{\rm FO}$. In Sec. 3.1.3, we evaluate the uncertainty in $N_{\rm
eff}$ arising from this uncertainty in $T_{\rm FO}$.

####  3.1.2 Equilibrium from Decays and Inverse Decays

Even when $T_{\chi\rm PT}<T_{\rm FO}<T_{g\rm PT}$, and hence the freeze-out
temperature is uncertain, it is still possible to infer the axion abundance at
$T_{\chi\rm PT}$ for sufficiently large $m_{a}$. Specifically, if the axion
decay rate is greater than Hubble at $T_{\chi\rm PT}$, then the axion
possesses a thermal distribution at $T_{\chi\rm PT}$. For example, when
$m_{a}\gg T$, the decay $a\rightarrow g+g$ can dominate over axion-gluon
scattering when $T>T_{g\rm PT}$, or, for example, $a\rightarrow 3\pi$ can
dominate at $T=T_{\chi\rm PT}$. The axion decay rate to strongly coupled
particles is given by

$\displaystyle\Gamma_{a\rightarrow\text{QCD}}=\frac{m_{a}^{3}}{f_{a}^{2}}\mathcal{F}_{c}(m_{a})\,,$
(19)

where $\mathcal{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 $\mathcal{F}_{c}(m_{a})$ as numerically computed in
Aloni:2018vki , which includes the following meson decay channels:
$a\rightarrow 3\pi,\,\pi\pi\gamma,\,\eta\pi\pi,\,KK\pi,\,\eta^{\prime}\pi\pi,$
$\,\rho\rho,\,\omega\omega,\,K^{*}\overline{K}^{*},$ and $\phi\phi$. For
$m_{a}\geq T_{g\rm PT}$, the axion to gluon decay rate dominates and
$\mathcal{F}_{c}$ smoothly interpolates to the perturbative result
Aloni:2018vki

$\displaystyle\mathcal{F}_{c}\simeq\frac{2}{\pi}\left(\frac{g_{3}^{2}}{32\pi^{2}}\right)^{2}\left(1+\frac{83g_{3}^{2}}{16\pi^{2}}\right)\,.$
(20)

Last, the axion to photon decay rate is given by

$\displaystyle\Gamma_{a\rightarrow\gamma\gamma}=\frac{g_{\gamma}(m_{a})^{2}}{64\pi}\frac{m_{a}^{3}}{f_{a}^{2}}\,,$
(21)

where $g_{\gamma}(m_{a})$ is given in (16). Note that $g_{\gamma}(m_{a})$ is a
function of the axion mass due to the effects of axion-meson mixing as encoded
in the mixing function $\mathcal{F}_{\theta}(m_{a})$.

We define the axion decay temperature, $T_{\rm decay}$, when
$\Gamma_{a\rightarrow\rm QCD}+\Gamma_{a\rightarrow\gamma\gamma}=3H(T_{\rm
decay})$. The yellow region of Fig. 1 shows the region where $T_{\rm
decay}>T_{\chi\rm PT}$. In this region, the axion possesses a thermal
distribution when we begin our Boltzmann code at $T_{\chi\rm PT}$, even if
$T_{\rm FO}$ is uncertain.

Note that if the axion decays far before or after neutrino decoupling, the
initial condition of the axion at $T_{\chi\rm PT}$ becomes insensitive to the
calculation of $\Delta N_{\rm eff}$. In particular, below the lower dashed
line, $\Gamma_{a}=\Gamma_{a\rightarrow\rm
QCD}+\Gamma_{a\rightarrow\gamma\gamma}$, 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 $\Delta N_{\rm eff}$ from the axion is erased,
regardless of the initial axion abundance at $T_{\chi\rm PT}$. Likewise, above
the upper dashed line. $\Gamma_{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
$\Delta N_{\rm eff}\ll-0.3$, which is already excluded by the observations of
CMB Planck:2018vyg .

####  3.1.3 Initial Condition and Its Uncertainty

Figure 2: When $T_{\chi\rm PT}<T_{\rm FO}<T_{\rm 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_{\rm FO}\in(T_{\chi\rm
PT},T_{\rm 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_{\rm FO}$ around $0.3$ GeV arises from the balance between
Boltzmann suppression and dilution by $g_{*}$ across the QCD phase transition.

For $(m_{a},f_{a})$ in the purple or yellow regions of Fig. 1, the axion is in
thermal equilbrium at $T_{\chi\rm PT}$. Consequently, in these regions, we
take the initial distribution function of the axion at $T_{\chi{\rm PT}}$ to
be a Bose-Einstein distribution of temperature $T_{\chi{\rm PT}}$. For
sufficiently large $f_{a}$, however, the axion decouples from the bath at
$T_{\rm FO}>T_{\chi{\rm PT}}$ (blue region), and we take the initial axion
distribution to be a thermal one at $T_{\rm FO}$, red-shifted down to
$T_{\chi{\rm PT}}$. Inside the red-shaded region, however, $T_{\rm FO}$ is
uncertain. Despite this uncertainty, we can still bound $\Delta N_{\rm eff}$
by running our axion Boltzmann code with both the maximum and minimum possible
axion abundance at $T_{\chi\rm PT}$. We scan over possible freeze-out
temperatures between $T_{\chi\rm PT}$ and $T_{g\rm PT}$ for all axion masses
in the red-shaded region to determine the smallest and largest abundance at
$T_{\chi\rm PT}$ as a function of $m_{a}$. In Fig. 2, the left panel shows the
$T_{\rm 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_{\chi\rm PT}$ and the largest that of a frozen-out abundance at $T_{\rm
FO}=T_{g\rm PT}$. This is because large changes in the degrees of freedom of
the thermal bath between $T_{\chi\rm PT}$ and $T_{g\rm PT}$ can dilute the
previously frozen-out axion. For all future plots, we show the results of
$\Delta N_{\rm eff}$ arising from these two possible initial conditions.

In summary, we take the initial axion distribution function at $T_{\chi\rm
PT}$ (time $t_{\chi\rm PT}$) to be

$\displaystyle f_{a}(\mathbf{p}_{a},t=t_{\chi\rm
PT})=\left(\exp{\dfrac{\sqrt{|\mathbf{p_{*}}|^{2}+m_{a}^{2}}}{T_{*}}}-1\right)^{-1},\quad|\mathbf{p_{*}}|=|\mathbf{p}_{a}|\frac{a(T_{*})}{a(T_{\chi\rm
PT})}$ (22)

where $f_{a}(\mathbf{p}_{a})$ is a Bose-Einstein distribution with momentum
$\mathbf{p}_{a}$ and effective temperature $T_{*}$ given by

$\displaystyle T_{*}=\begin{cases}T_{{\rm FO},g}\quad&\text{if }T_{\rm FO}\geq
T_{g\rm PT}\\\ T_{\chi\rm PT}\quad&\text{if }T_{\rm FO}\leq T_{\chi\rm
PT}\text{ or }T_{\rm decay}\geq T_{\chi\rm PT}\\\ \text{Max \& Min }T_{\rm
FO}\text{ (see Fig.~{}\ref{fig:TFOUncertainty})}\quad&\text{if }T_{\chi\rm
PT}<T_{\rm FO}<T_{g\rm PT}\text{ and }T_{\rm decay}<T_{\chi\rm PT}\end{cases}$
(23)

Note that if $T_{\rm RH}<T_{*}$ or if there is a source of dilution in the
universe after the axion freezes-out, then the initial abundance of axions can
be small and the bounds on $N_{\rm eff}$ discussed in this work are weakened.

###  3.2 Axion Abundance Below $\Lambda_{\rm 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}(\mathbf{p}_{a})$, is

$\displaystyle\frac{\partial f_{a}}{\partial t}-p_{a}H\frac{\partial
f_{a}}{\partial p_{a}}$
$\displaystyle=(C_{\gamma}+C_{P}+C_{\pi}+C_{\Gamma})(f_{a,\rm eq}-f_{a}),$
(24)

where $p_{a}=|\mathbf{p}_{a}|$ is the magnitude of the axion momentum, $H$ is
the Hubble expansion rate, and $C_{\gamma}$, $C_{P}$, $C_{\pi}$, and
$C_{\Gamma}$ 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+A+B+...\leftrightarrow I+J+...$, is

$\displaystyle C(f_{a,\rm eq}-f_{a})$ $\displaystyle=\frac{1}{2E_{a}}\int
d\Pi_{A}d\Pi_{B}...d\Pi_{I}...d\Pi_{J}...S\,|\mathcal{M}|^{2}\Lambda\,(2\pi)^{4}\delta^{4}(p_{a}+p_{A}+p_{B}...-p_{I}-p_{J}...)$
(25)

where $E_{a}=\sqrt{p_{a}^{2}+m_{a}^{2}}$ is the axion energy,
$d\Pi=d^{3}p/(2\pi)^{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

$\displaystyle\Lambda$ $\displaystyle=\left[(1\pm f_{A})(1\pm f_{B})...(1\pm
f_{a})f_{I}f_{J}...-f_{A}f_{B}...f_{a}(1\pm f_{I})(1\pm f_{J})...\right]$ (26)
$\displaystyle\simeq(f_{a,\rm
eq}-f_{a})\exp(-E_{A}/T)\exp(-E_{B}/T)...\quad\text{(Kinetic equilibrium
limit)}$

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. (26)
shows $\Lambda$ in the limit where $1\pm f\simeq 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_{\gamma}$, 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 Cadamuro:2010cz

$\displaystyle C_{\gamma}$
$\displaystyle\approx\frac{m_{a}^{2}-4m_{\gamma}^{2}}{m_{a}^{2}}\frac{m_{a}}{E_{a}}\left[1+\frac{2T}{p_{a}}\log\frac{1-e^{-(E_{a}+p_{a})/2T}}{1-e^{-(E_{a}-p_{a})/2T}}\right]\Gamma_{a\rightarrow\gamma\gamma}$
(27) $\displaystyle C_{P}$
$\displaystyle\approx\frac{\Gamma_{a\rightarrow\gamma\gamma}}{m_{a}^{3}}\sum_{i=e,\mu,\pi_{\pm}}n_{i}e_{i}^{2}\log\left[1+\frac{16E_{a}^{2}(m_{i}+3T)^{2}}{m_{\gamma}^{2}(m_{i}^{2}+(m_{i}+3T)^{2})}\right].$
(28)

Here, $\Gamma_{a\rightarrow\gamma\gamma}$ is the axion-to-two photon decay
rate, (21), $m_{\gamma}$ is the photon plasma mass, and $n_{i}$ is the number
density of the $i$th electromagnetically charged particle of charge $e_{i}$ in
the thermal bath. When the electron is relativistic, $m_{\gamma}\simeq eT/3$
weldon1982covariant , but when the electron becomes non-relativistic at
$T\lesssim m_{e}$, $m_{\gamma}$ reduces to the classical plasma frequency of
$m_{\gamma}\simeq e\sqrt{n_{e}/m_{e}}$, which is exponentially suppressed. For
simplicity, we piecewise-connect the two regimes for $m_{\gamma}$ when they
intersect, which occurs roughly at $T\simeq 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_{\pi}$ takes the form

$\displaystyle
C_{\pi}=\left(\frac{A}{f_{a}f_{\pi}}\frac{1}{1-r^{2}}\right)^{2}\frac{T^{6}}{2E_{a}}\times\mathcal{F}_{\rm
PS}\left(\frac{m_{a}}{T},\frac{p_{a}}{T}\right)\,,$ (29)

where $A=\frac{1}{3}(1-z)/(1+z)\simeq 0.12$ and $r\equiv m_{a}/m_{\pi}$ as
before. The function $\mathcal{F}_{\rm 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_{\pi}$ agrees with the massless axion result found in
literature Hannestad_2005 ; DiLuzio:2021vjd . Note that integration of
$C_{\pi}$ over the axion phase space defines $\Gamma_{a\pi\leftrightarrow
a\pi}$ given in Eq. (18).

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}\gtrsim 3\times 10^{4}$ GeV, the axion-pion
scattering dominates over Primakoff scattering, keeping the axions in thermal
equilbrium until $T\sim 10$ MeV for $f_{a}\lesssim 10^{6}$ GeV. Figure 4:
(Right) Total axion decay rate ($\tau^{-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_{\pi^{0}}$, $m_{\eta}$, and
$m_{\eta}^{\prime}$. Troughs arise from cancellations in $g_{\gamma}$ between
the anomaly, axial rotation, and meson-mixing contributions. For $m_{a}\gtrsim
2$ GeV, the decay rate is set by gluons. The dashed contours show the axion-
photon decay rate without meson mixing.

The axion-meson decay collision term has also never been considered in the
literature. Under the Boltzmann approximation in the kinetic equilbrium limit,
$C_{\Gamma}$ is simply

$\displaystyle C_{\Gamma}\approx\frac{m_{a}}{E_{a}}\Gamma_{\rm QCD}\,,$ (30)

where $\Gamma_{\rm QCD}$ is the total axion decay rate to all final states
containing mesons (19). In (30), 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. 4 show the axion-pion decoupling temperature,
$T_{\rm FO,\pi}$ vs $f_{a}$ over a range of axion masses. For $f_{a}\lesssim
10^{6}\,{\rm GeV}$, the axion remains in thermal equilibrium from pionic
interactions until a few $10$s of MeV which is typically far lower than the
Primakoff decoupling temperature Cadamuro:2011fd ,

$\displaystyle T_{\rm FO,P}\approx 91\,\frac{\sqrt{g_{*}(T_{\rm
FO,P})}}{g_{c}(T_{\rm FO,P})}\left(\frac{f_{a}}{10^{6}\,\rm
GeV}\right)^{2}\left(\frac{g_{\gamma}}{\alpha/2\pi}\right)^{-2}\,{\rm GeV}\,,$
(31)

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. 4. For $f_{a}\gtrsim 5\times 10^{4}\,{\rm GeV}$,
$T_{\rm FO,P}\gg T_{\rm FO,\pi}$, 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_{\rm FO,\pi}$, 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_{\chi\rm 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 $\Gamma_{a\rightarrow\gamma\gamma}^{-1}$, 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
$\Gamma_{a\rightarrow\gamma\gamma}$ 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_{\rm 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, (24), 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

$\displaystyle N_{\rm
eff}=\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\frac{\rho_{\nu}}{\rho_{\gamma}}\,.$
(32)

Here, $\rho_{\nu}$ and $\rho_{\gamma}$ are the neutrino and photon energy
densities. Note that $N_{\rm 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_{\rm\chi 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, $\rho_{\delta_{\rm QCD}}$, is approximately $5\%$ of the
total energy density at $T_{\rm\chi 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

$\displaystyle\sum_{i=\gamma,e,\mu,\pi,\delta_{\rm
QCD}}\frac{\partial\rho_{i}}{\partial t}+3H(\rho_{i}+P_{i})$
$\displaystyle=\frac{\Gamma_{\nu
e}}{T^{4}}(\rho_{\nu_{e}}^{2}-\rho_{\nu_{e},\rm
eq}^{2})+\frac{\Gamma_{\nu_{\mu\tau}}}{T^{4}}(\rho_{\nu_{\mu\tau}}^{2}-\rho_{\nu_{\mu\tau},\rm
eq}^{2})$ (33)
$\displaystyle+\int\frac{d^{3}p_{a}}{(2\pi)^{3}}\sqrt{p_{a}^{2}+m_{a}^{2}}(C_{\gamma}+C_{P}+C_{\pi}+C_{\Gamma})(f_{a}-f_{a,\rm
eq}),$

where

$\displaystyle\rho_{i}=\int\frac{d^{3}p_{i}}{(2\pi)^{3}}\sqrt{p_{i}^{2}+m_{i}^{2}}\,f_{i}(p_{i}),\hskip
21.68121ptP_{i}=\int\frac{d^{3}p_{i}}{(2\pi)^{3}}\frac{1}{3}\frac{p_{i}^{2}}{\sqrt{p_{i}^{2}+m_{i}^{2}}}\,f_{i}(p_{i})$
(34)

are the energy densities and pressures of the $i$th 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 $\rho_{\delta_{\rm QCD}}$ and
$P_{\delta_{\rm QCD}}$ from the work of Saikawa:2018rcs , 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 (33) follow the Boltzmann equations

$\displaystyle\frac{\partial\rho_{\nu_{e}}}{\partial t}+4H\rho_{\nu_{e}}$
$\displaystyle=-\frac{\Gamma_{\nu
e}}{T^{4}}(\rho_{\nu_{e}}^{2}-\rho_{\nu_{e},\rm eq}^{2})$ (35)
$\displaystyle\frac{\partial\rho_{\nu_{\mu\tau}}}{\partial
t}+4H\rho_{\nu_{\mu\tau}}$
$\displaystyle=-\frac{\Gamma_{\nu_{\mu\tau}}}{T^{4}}(\rho_{\nu_{\mu\tau}}^{2}-\rho_{\nu_{\mu\tau},\rm
eq}^{2}).$ (36)

Here, $\Gamma_{\nu_{e}}\simeq 0.68\,G_{F}^{2}T^{5}$ and
$\Gamma_{\nu_{\mu\tau}}\simeq 0.15\,G_{F}^{2}T^{5}$ Cadamuro:2011fd  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,

$\displaystyle H^{2}=\left(\frac{\dot{R}}{R}\right)^{2}$
$\displaystyle=\frac{8\pi
G}{3}\left(\rho_{\gamma}+\rho_{e}+\rho_{\mu}+\rho_{\pi}+\rho_{\delta_{\rm
QCD}}+\rho_{\nu_{e}}+\rho_{\nu_{\mu\tau}}+\rho_{a}\right).$ (37)

We numerically solve the system of equations (24)-(37) using the method of
lines methodOfLines . 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 $\\{|\mathbf{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.

###  3.3 $\Delta N_{\rm eff}$

In this section, we present the numerical results of $N_{\rm eff}$ as
determined from the Boltzmann equations of Sec. 3.2 describing the
cosmological evolution of the heavy QCD axion below $T_{\chi\rm PT}$.

Figure 5:  Contours of $N_{\rm eff}$ in the $m_{a}-\tau_{\gamma}$ (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_{\gamma}$, 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 axion-lifetime 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.

First, to compare with past literature and highlight the importance of the new
effects discussed in this work, we show $N_{\rm 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_{\gamma}$, the proper frozen-out initial axion abundance, and the QCD
contributions to the background energy density as described by $\rho_{\pi}$
and $\delta_{\rm QCD}$. Fig. 5 shows the numerical results of $N_{\rm eff}$
when these terms are neglected. That is, including only the photon
($C_{\gamma}$) and Primakoff ($C_{P}$) collision terms in Eq. (24) and Eq.
(33); including only $\gamma,e$ and $\mu$ in the sum of the thermally coupled
species of Eq. (33); including only
$\rho_{\gamma},\rho_{e},\rho_{\mu},\rho_{\nu_{e}},\rho_{\nu_{\mu\tau}}$, and
$\rho_{a}$ in Hubble (37); 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_{*}={\rm Max}(T_{\rm FO,P},T_{\chi\rm PT})$; and
finally, in Eq. (16), setting $E/N=0$ and $\mathcal{F_{\theta}}=(1-z)/(1+z)$,
which is the axion-meson mixing contribution in the inapplicable $m_{a}\ll
m_{\pi}$ limit. The left panel of Fig. 5 shows contours of $\Delta N_{\rm
eff}$ in the $m_{a}-\tau_{\gamma}$ plane, where $\tau_{\gamma}$ is the axion
to photon lifetime when axion-meson mixing is neglected. Note that
$\tau_{\gamma}$ is the total lifetime of the axion since QCD decay channels
are neglected in this particular case.

For $\tau_{\gamma}\gtrsim 10^{-1}$ s the axion decays after neutrino
decoupling, heating up the photons relative to the neutrinos and giving rise
to $\Delta N_{\rm eff}<0$ as can be seen from the enhanced denominator of Eq.
(32). For $\tau_{\gamma}\lesssim 10^{-1}$ s and $m_{a}\gtrsim 10$ MeV, the
axion decays sufficiently early that the photons and neutrinos rethermalize
before the neutrino decouples. In this scenario, the $N_{\rm eff}$ signal of
the heavy axion is absent and $\Delta N_{\rm eff}\simeq 0$. For
$\tau_{\gamma}\lesssim 10^{-1}$ s and $m_{a}\lesssim 10$ MeV, the axion
remains in thermal equilibrium past neutrino decoupling, heating up the
photons and again giving rise to negative $\Delta N_{\rm eff}$. The right
panel of Fig. 5 shows the same $\Delta N_{\rm 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_{\rm eff}$ for heavy axion cosmologies incorporating the
new effects included in this work: axion-pion scattering, axion-meson decay,
axion-meson mixing in $g_{\gamma}$, the proper (and occasionally uncertain)
frozen-out initial axion abundance, the QCD contributions to the background
energy density as described by $\rho_{\pi}$ and $\delta_{\rm QCD}$, as well
results for the KSVZ $E/N=0$ and the GUT motivated $E/N=8/3$. Figs. 7 and 7
show the contours of $\Delta N_{\rm 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},\tau)$ and $(m_{a},f_{a})$ planes, respectively. Note that here,
$\tau$ is the total lifetime of the axion, which begins differing from
$\tau_{\gamma}$ for $m_{a}\gtrsim m_{\pi}$ where axion-meson mixing becomes
important and then becomes even more disparate when axion-meson channels open
for $m_{a}\gtrsim 3m_{\pi}$, 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_{\rm FO}$ lies in the
uncertain region between $T_{\chi\rm PT}$ and $T_{g\rm PT}$, as indicated in
Fig. 2. The separation between the solid and dashed blue contours indicates
the uncertainty in $N_{\rm eff}$ arising from the uncertainty in $T_{\rm FO}$
in this region. As can be seen, this region is localized roughly between
$250\,{\rm MeV}\lesssim m_{a}\lesssim 800\,{\rm MeV}$ and, for any value of
$m_{a}$ in this region, the uncertainty in the value of $f_{a}$ for any
$N_{\rm eff}$ contour is typically only a several 10s of percent and always
less than a factor 3.

Figure 6: Contours of $N_{\rm eff}$ in the $m_{a}-\tau$ (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_{\rm FO}$ lies between
$T_{\chi\rm PT}$ and $T_{\rm gPT}$. The dashed contours bound this uncertainty
by showing the value of $N_{\rm eff}$ taking the minimum initial axion yield
while the solid contours show the value of $N_{\rm eff}$ taking the maximum
initial axion yield as given in Fig. 2.

Figure 7:  Same as Fig. 7, but for E/N = 8/3.

Figs. 7 and 7 demonstrate three important differences in the $N_{\rm 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_{\pi}^{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}\gtrsim 3m_{\pi}$, $|\Delta N_{\rm 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_{\pi},\,m_{\eta}$, and $m_{\eta^{\prime}}$. 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_{\gamma}$ (21). The increased
$|\Delta N_{\rm eff}|$ for fixed $(m_{a},f_{a})$ due to the increased axion
lifetime is important for $m_{a}<3m_{\pi}$ when only the axion-photon decay
channel is open. Third, for $m_{a}\lesssim m_{\pi}$, $|\Delta N_{\rm 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 $|\Delta N_{\rm 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\times 10^{6}$
GeV) and ($1000$ MeV, $1.0\times 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}=\rho_{i}R^{4}$ of the $i$th 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, $\gamma,e,\mu,\pi,\delta_{\rm
QCD}$ is shown in blue, $\nu_{\mu,\tau}$ in green, $\nu_{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.

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},\tau)=(100\,{\rm MeV}$, $0.1{\rm s})$, left,
$(m_{a},\tau)=(1000\,{\rm MeV}$, $0.1{\rm s})$, 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. (32) minus the Standard Model result, $N_{\rm
eff}^{\rm SM}=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_{\chi\rm 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_{\chi\rm PT}$. In the left panel, the axion lifetime
with or without mesons is approximately the same ($\tau\approx 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 ($\tau\approx 0.01$ s with meson interactions and $24$ s
without). This leads to a reduced $\Delta N_{\rm eff}$, as shown on the bottom
right panel. In both panels, we take $E/N=0$.

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 leads to a relative suppression in the non-relativistic abundance of the
axion prior to decaying around $2$ MeV. Consequently, including axion-pion
scattering $\Delta N_{\rm eff}$ is not as negative as previous results in the
literature. The temporal evolution of $\Delta N_{\rm eff}(t)=N_{\rm
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_{\pi}$,
the axion follows the dashed thermal abundance for a longer duration which
generates the large triangular peak in the allowed $N_{\rm eff}$ plot of Fig.
7.

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\rightarrow\eta\pi\pi,\pi\pi\gamma$ 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 $\Delta N_{\rm eff}$
as shown explicitly by the evolution of $\Delta N_{\rm 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\sigma$ limit on $N_{\rm eff}$ as constrained by Planck approximately
follows the contour $\tau\approx 0.05$ s around $m_{a}\sim$ GeV and slowly
drops with increasing axion mass. The slight decrease in the maximum allowed
$\tau$ 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,
$\rho_{a}\propto e^{-t/\tau}$, leads to a logarithmic decrease in the maximum
allowed $\tau$ given by the semi-analytic function

$\displaystyle\tau_{\rm max}(m_{a})<\frac{4.3\times 10^{-2}\,\rm
s}{1+0.25\ln{\frac{m_{a}}{2\rm\,GeV}}}\qquad(m_{a}\geq 2\,\rm GeV).$ (38)

Eq. (38) can also be written in terms of $f_{a}$ by equating $\tau_{\rm
max}^{-1}$ with the analytic decay rate into gluons given by Eqns. (19) and
(20),

$\displaystyle f_{a}\lesssim 1.2\times 10^{10}\,{\rm
GeV}\left(\frac{m_{a}}{2\,\rm
GeV}\right)^{\scalebox{1.01}{$\frac{3}{2}$}}\left(\frac{\alpha_{3}}{0.3}\right)\left(\frac{1+\frac{83}{4\pi}\alpha_{3}}{1+\frac{83}{4\pi}0.3}\right)^{\scalebox{1.01}{$\frac{1}{2}$}}\left(1+0.25\ln\frac{m_{a}}{2\,\rm
GeV}\right)^{\scalebox{1.01}{$-\frac{1}{2}$}}\quad(m_{a}\geq 2\,\rm GeV).$
(39)

Thus, for $m_{a}\gg 2$ GeV and $f_{a}$ small enough to be probed by
accelerator experiments Ertas:2020xcc ; Dobrich:2015jyk ; Dolan:2017osp ;
NA64:2020qwq ; FASER:2018eoc ; Gori:2020xvq ; Aloni:2018vki ; Mariotti:2017vtv
; Kelly:2020dda ; Chakraborty:2021wda ; Bertholet:2021hjl , the dark radiation
constraint is absent.

##  4 Including a Mirror Photon

In this section we add a mirror photon $\gamma^{\prime}$ 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}\times 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

$\displaystyle{\cal L}_{a,\gamma^{\prime}}$
$\displaystyle=\frac{g_{\gamma^{\prime}}}{4}\frac{a}{f_{a}}F^{\prime}_{\mu\nu}\tilde{F^{\prime}}^{\mu\nu}.$
(40)

Above the mirror QCD scale, $\Lambda_{\rm QCD}^{{}^{\prime}}$, the $Z_{2}$
symmetry ensures that

$\displaystyle g_{\gamma^{\prime}}=\frac{e^{\prime
2}}{8\pi^{2}}\left(\frac{E}{N}\right)\,,$ (41)

with $e^{\prime}$ differing from $e$ only by renormalization group scaling,
which we ignore.

Note that unlike $g_{\gamma}$ in (21), $g_{\gamma^{\prime}}$ 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 $\Lambda_{\rm QCD}^{{}^{\prime}}$ and thus irrelevant to the
theory below $\Lambda_{\rm QCD}^{{}^{\prime}}$. If, however, any mirror quark
is lighter than the mirror QCD scale, then there is an additional contribution
to $g_{\gamma^{\prime}}$ analogous to the second and third terms in (16) for
$g_{\gamma}$. 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^{\prime}\gg v$, all mirror quarks are
heavier than $\Lambda_{\rm QCD}^{{}^{\prime}}$ for $f_{a}m_{a}\gtrsim 25\;{\rm
GeV}^{2}$ Dunsky:2019upk . 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_{\gamma^{\prime}}$ is thus uncorrected and given by (41).

Numerical results in this section are calculated taking $e^{\prime}=e$ (that
is, neglecting the small running of $e^{\prime}$ below $v^{\prime}$) and with
a non-zero $E/N$ so that $g_{\gamma}^{\prime}$ is non-zero. As a result, the
axion-mirror photon decay rate is

$\displaystyle\Gamma_{a\rightarrow\gamma^{\prime}\gamma^{\prime}}=\frac{g_{\gamma^{\prime}}^{2}}{64\pi}\frac{m_{a}^{3}}{f_{a}^{2}}\simeq\frac{1}{64\pi}\left(\frac{e^{2}}{8\pi^{2}}\frac{E}{N}\right)^{2}\frac{m_{a}^{3}}{f_{a}^{2}}\,.$
(42)

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.

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.

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 $\gamma^{\prime}$, $\gamma$, and QCD degrees of freedom,
respectively. Because there is no cancellation of terms in
$g_{\gamma^{\prime}}$ as compared to $g_{\gamma}$, the
$a\rightarrow\gamma^{\prime}\gamma^{\prime}$ decay rate (blue) is roughly an
order of magnitude greater than the $a\rightarrow\gamma\gamma$ decay rate
(orange) until $m_{a}\gtrsim 3m_{\pi}$: 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_{\gamma}$ is negligible which leaves the now smaller
$a\rightarrow\gamma^{\prime}\gamma^{\prime}$ decay rate roughly an order of
magnitude weaker than the $a\rightarrow\gamma\gamma$ decay rate.

We highlight these two representative values of $E/N$ since they generate
substantially different decay branching ratios into mirror photons.222The case
for $E/N=2/3$, which can be achieved by a KSVZ quark that possesses the same
gauge quantum numbers as down quarks, is intermediate between these two cases
and like the case $E/N=1/3$, also yields a region of $\Delta N_{\rm
eff}\approx 0$ at low $m_{a}$. This disparity is important since the parameter
space where the axion branching ratio into dark photons is $\mathcal{O}(1)$
can be cosmologically dangerous as the mirror photon decay mode (42) increases
$N_{\rm eff}$, by directly generating dark radiation in the form of
$\gamma^{\prime}$, and reduces the heating of the Standard Model bath as fewer
axions decay into $\gamma$.

Quantitatively, $N_{\rm eff}$ with a mirror photon is

$\displaystyle N_{\rm
eff}=\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\frac{\rho_{\nu}+\rho_{\gamma^{\prime}}}{\rho_{\gamma}}\,,$
(43)

where $\rho_{\nu}$, $\rho_{\gamma^{\prime}}$, and $\rho_{\gamma}$, are the
relic energy densities of neutrinos, mirror photons, and photons,
respectively. According to Eq. (43), an additional mirror photon in thermal
equilibrium significantly increases $N_{\rm eff}$ and is generally excluded by
current $\Delta N_{\rm eff}$ limits Planck:2018vyg . 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}(\mathbf{p})$ in
the freeze-in picture

$\displaystyle\frac{\partial f_{a}}{\partial t}-p_{a}H\frac{\partial
f_{a}}{\partial p_{a}}$
$\displaystyle=(C_{\gamma}+C_{P}+C_{\pi}+C_{\Gamma})(f_{a,\rm
eq}-f_{a})-C_{\gamma^{\prime}}f_{a}\,,$ (44)

where the mirror photon collision term, $C_{\gamma^{\prime}}$, is given by

$\displaystyle
C_{\gamma^{\prime}}\simeq\frac{m_{a}}{E_{a}}\Gamma_{a\rightarrow\gamma^{\prime}\gamma^{\prime}}\,.$
(45)

Eq. (44) replaces Eq. (24) when including mirror photons. In addition, while
the energy density evolution of the Standard Model bath remains as given in
Eqns. (33) and (35), 333The axion-mediated interaction
$\gamma^{\prime}+\gamma^{\prime}\leftrightarrow\gamma+\gamma$ can contribute
to the energy transfer to the Standard Model thermal bath and hence to
additional terms on the right-hand-side of Eq .(33). However, this interaction
is $\mathcal{O}(1/f_{a}^{4})$ and generally negligible. Similarly, the decay
and inverse decay $a\leftrightarrow\gamma+\gamma^{\prime}$ can also modify
(33), but this requires $O(1)$ kinetic mixing or $aF^{\prime}\tilde{F}$. The
former is constrained by searches for relic mirror charged particle
Dunsky:2018mqs , and the latter requires $U(1)\times U(1)^{\prime}$ charged
particles around the mass scale $f_{a}$, whose relic is also constrained. the
energy density evolution of ${\gamma^{\prime}}$ is described by

$\displaystyle\frac{\partial\rho_{\gamma^{\prime}}}{\partial
t}+4H\rho_{\gamma}^{\prime}=\int\frac{d^{3}p}{(2\pi)^{3}}m_{a}\Gamma_{a\rightarrow\gamma^{\prime}\gamma^{\prime}}f_{a}=m_{a}\Gamma_{a\rightarrow\gamma^{\prime}\gamma^{\prime}}\,n_{a}\,.$
(46)

Figure 11: Contours of $N_{\rm eff}$ in the $m_{a}-\tau$ (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_{\rm FO}$ lies between $T_{\chi\rm PT}$ and $T_{\rm gPT}$. The dashed
contours bound this uncertainty by showing the value of $N_{\rm eff}$ taking
the minimum initial axion yield while the solid contours show the value of
$N_{\rm 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 12:  Same as Fig. 12, but for $E/N=1/3$.

Conservatively, we take the initial $\gamma^{\prime}$ density at $T_{\chi\rm
PT}$ to be zero. Due to the substantial change in Standard Model degrees of
freedom across $T_{\rm QCD}$, freeze-in production of $\gamma^{\prime}$ much
earlier than $T_{\chi\rm PT}$ is diluted and this conservative estimate is a
fairly good approximation to the true initial abundance of $\gamma^{\prime}$.
Last, the Hubble expansion rate, (37), is modified to include the additional
mirror photon energy density,

$\displaystyle H^{2}=\left(\frac{\dot{R}}{R}\right)^{2}$
$\displaystyle=\frac{8\pi
G}{3}\left(\rho_{\gamma}+\rho_{e}+\rho_{\mu}+\rho_{\pi}+\rho_{\delta_{\rm
QCD}}+\rho_{\nu_{e}}+\rho_{\nu_{\mu\tau}}+\rho_{a}+\rho_{\gamma^{\prime}}\right).$
(47)

Figures 12 and 12 show the contours of $N_{\rm 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_{\rm eff}$ is
excluded by current CMB measurements at the 2$\sigma$ level. Note that whereas
$\Delta N_{\rm eff}$ is strictly negative in the case without the mirror
photon (Figs. 7 and 7), the case with the mirror photon gives positive $\Delta
N_{\rm 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
$\Gamma_{a\gamma^{\prime}\gamma^{\prime}}\gtrsim H(T=m_{a})$, or equivalently,
roughly when

$\displaystyle f_{a}\lesssim 10^{5}{\,\rm GeV}\left(\frac{m_{a}}{125\,{\rm
MeV}}\right)^{\scalebox{1.01}{$\frac{1}{2}$}}\left(\frac{E/N}{8/3}\right).\qquad(\text{Mirror
Photon reaches thermal equil.})$ (48)

Within this green region, the contour values for $N_{\rm eff}$ are
artificially high because the mirror photon acquires a greater than thermal
abundance due to the lack of a back reaction in Eq. (44). In principle,
capping the mirror photon abundance at a thermal abundance suggests that
realistic $N_{\rm eff}$ contours within the green region are roughly fixed at
the value of $N_{\rm eff}$ on the boundary of the green region; that is, the
value of $N_{\rm 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_{\rm eff}\gtrsim 3.8$ in the green
region when $E/N=8/3$. Such a large $\Delta N_{\rm 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
$\rho_{\gamma^{\prime}}$ energy deposit (positive $\Delta N_{\rm eff}$
contribution) and the heating of $\rho_{\gamma}$ relative to neutrinos
(negative $\Delta N_{\rm 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 $\tau\gtrsim 1$ sec is excluded by BBN.
This is because to cancel the positive $\Delta N_{\rm eff}$, the axion decays
before the proton-neutron conversion completes, and the Helium abundance will
be affected.

##  5 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 $\Delta N_{\rm
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.

Figure 13: Cosmological constraints on the low mass axion region, $m_{a}<1$
MeV. The axion decays prior to recombination ($t_{\rm CMB}\approx 370,000$ yr)
in the orange region, which is excluded by the $N_{\rm eff}$ constraint. The
boundary denoting $\tau=t_{\rm 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 Planck:2018vyg , assuming the axion
froze-out in the early Universe. For $f_{a}\lesssim 10^{8}$ GeV, $T_{\rm FO}$
occurs below the electroweak scale, and hence the axion energy density at
$t_{\rm 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.

As can be seen from Fig. 7 and 7, the CMB limit on $\Delta N_{\rm eff}$
excludes 1 MeV $\lesssim m_{a}\lesssim 3$ MeV for any $f_{a}$. This exclusion
from $\Delta N_{\rm eff}$ continues for $m_{a}<1$ MeV, until the axions decay
after recombination.444The parameter space where the axions decay before
recombination but are in thermal equilibrium at BBN and generate a positive
$\Delta N_{\rm eff}$ is discussed in Salvio:2013iaa ; DEramo:2021lgb ;
DEramo:2021psx , but is weaker than the bounds in this work. Hence, the blue
region of Fig. 13, where $\tau<t_{\rm CMB}\equiv 370,000$ yrs ($z\approx
1100$) ParticleDataGroup:2020ssz , is excluded by $\Delta N_{\rm eff}$. For
$\tau\gtrsim t_{\rm CMB}$ and $m_{a}>0.13(\Omega_{\rm
DM}h^{2}/0.12)Y_{a}(T_{\rm 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_{\rm FO,g})$, typically
between $2\times 10^{-3}$ to $2\times 10^{-2}$. Relaxing this assumption, by
taking $T_{\rm RH}$ below $T_{\rm 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 Xu:2021rwg
for our framework and exclude the green-shaded region in Fig. 13. Here we
conservatively impose the bound only for $\tau>t_{0}$, where $t_{0}$ is the
present age of the Universe, but we expect that the bound is also applicable
as long as $\tau>t_{\rm 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 ciaran_o_hare_2020_3932430  for an overview),
but they are generically weaker.

Finally, the misalignment mechanism Preskill:1982cy ; Abbott:1982af ;
Dine:1982ah  overproduces axion dark matter above the dotted-dashed line in
Fig. 13. Here we assume an $O(1)$ misalignment angle.

##  6 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}\sim(100\;{\rm 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_{\rm eff}$, which has been
precisely measured by the Planck Collaboration Planck:2018vyg  and will be
significantly improved by CMB Stage 4 experiments CMB-S4:2016ple . 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_{\rm 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_{\rm eff}$ and axion decays to mirror photons directly
enhancing $N_{\rm eff}$. We have also provided a detailed analysis of the
$N_{\rm eff}$ bound in theories with a light mirror photon.

Our analysis of the $N_{\rm 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 axion-meson 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.

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\sigma$ exclusion region on $N_{\rm eff}$ as determined by
Planck for the heavy QCD axion with $E/N=0$ (see Fig. 7). The dashed orange
contour shows the future reach of $N_{\rm 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_{\rm 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
Ertas:2020xcc ; Dobrich:2015jyk ; Dolan:2017osp ; NA64:2020qwq ; FASER:2018eoc
; Gori:2020xvq ; Aloni:2018vki ; Mariotti:2017vtv ; Kelly:2020dda ;
Chakraborty:2021wda ; Bertholet:2021hjl . The brown region shows the bound
from supernova 1987a on axions with hadronic interactions, Chang:2018rso ,
green from horizontal-branch cooling Ayala:2014pea , red from CAST
CAST:2007jps ; CAST:2017uph , and pink from solar neutrinos Vinyoles:2015aba .
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_{\rm 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.

Our results for the CMB $N_{\rm eff}$ constraints on the heavy QCD axion, in
the absence of a mirror photon, are shown in Figs. 7 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. 7 and 7. First, resonances
occur when the axion mass is around the $\pi_{0}$, $\eta$, and $\eta^{\prime}$
masses, greatly affecting $N_{\rm 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 $\sim
1$ GeV; this point is apparent in the right panels of Figs. 7 and 7 where
regions with higher $f_{a}$ open up.

In the presence of a light mirror photon, our results for the CMB $N_{\rm
eff}$ constraints are shown in Figs. 12 and 12 for $E/N$ = 8/3 and 1/3
respectively. For $E/N=8/3$, $m_{a}<100\;{\rm MeV}$ is excluded for all values
of $f_{a}$. A substantial fraction of the allowed region with $100\;{\rm
MeV}<m_{a}<500\;{\rm MeV}$ will be probed by CMB-S4 via a positive signal for
$\Delta N_{\rm eff}$. For $E/N=1/3$, the CMB $N_{\rm eff}$ bound is
considerably weaker. A new allowed region opens up at lower axion masses,
$1\;{\rm MeV}<m_{a}<100\;{\rm 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 $\Delta N_{\rm 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_{\rm
eff}$ in this work are shown in comparison to other cosmological and
astrophysical constraints in Fig. 14. The limits on $N_{\rm eff}$ in this work
provide the strongest constraints on heavy QCD axions for $m_{a}\gtrsim 1$ MeV
and $f_{a}\gtrsim 10^{5}$ GeV. Complementary constraints at small $f_{a}$
arise from direct heavy axion searches at accelerators Kelly:2020dda ;
Bertholet:2021hjl  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 allow the region enclosed by the dashed purple contour to be
probed Kelly:2020dda . 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 Bertholet:2021hjl . 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
Chang:2018rso ; we show a conservative case. Constraints on the decay of the
axion from extragalactic background light or CMB spectral distortions are
derived in Cadamuro:2011fd , 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_{\rm 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_{\rm RH}>4$ MeV Kawasaki:1999na ;
Kawasaki:2000en ; Hasegawa:2019jsa  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. (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 $\Delta N_{\rm eff}$
signal DEramo:2021psx , heavy QCD axions can generate substantial $\Delta
N_{\rm 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 $\Delta
N_{\rm eff}$. Such a measurement would determine a correlation between the
axion mass and decay constant.

## Acknowledgement

We thank Jeffrey Anderson, Raymond Co, and Jonathan Wurtele for providing
useful computational resources. This work was supported in part by the
Director, Office of Science, Office of High Energy and Nuclear Physics, of the
US Department of Energy under Contracts DE-AC02-05CH11231 (LJH) and by the
National Science Foundation under grant PHY-1915314 (LJH).

##  Appendix A Calculation of Axion-Pion Collision Term

In this section, we compute the axion-pion scattering collision term,
$C_{\pi}$, used in the axion-Boltzmann equation, (33). The interaction between
an axion, neutral pion, and two charged pions includes the following three
interactions,
$a(p_{a})+\pi_{0}(p_{1})\leftrightarrow\pi_{+}(p_{2})+\pi_{-}(p_{3})$,
$a(p_{a})+\pi_{-}(p_{1})\leftrightarrow\pi_{0}(p_{2})+\pi_{-}(p_{3})$, and
$a(p_{a})+\pi_{+}(p_{1})\leftrightarrow\pi_{+}(p_{2})+\pi_{0}(p_{3})$ as shown
in Fig. 15. According to the chiral Lagrangian (2), the matrix element for
$a+\pi_{0}\rightarrow\pi_{+}+\pi_{-}$ is

$\displaystyle\mathcal{M}_{a+\pi_{0}\rightarrow\pi_{+}+\pi_{-}}=\frac{3}{2}\frac{A}{f_{a}f_{\pi}}\frac{1}{1-r^{2}}(s-m_{\pi}^{2}),$
(49)

where $s=(p_{a}+p_{1})^{2}$. The matrix elements for
$a+\pi_{-}\rightarrow\pi_{-}+\pi_{0}$ and
$a+\pi_{+}\rightarrow\pi_{+}+\pi_{0}$ scattering are obtained by the four
momentum mapping $p_{1}\rightarrow-p_{2}$ and $p_{1}\rightarrow-p_{3}$,
respectively. As a result, the total squared amplitude for scatterings
involving the axion and charged pions is

$\displaystyle|\mathcal{M}|^{2}_{a+\pi_{i}\rightarrow\pi_{j}^{\dagger}+\pi_{k}^{\dagger}}$
$\displaystyle=|\mathcal{M}|^{2}_{a+\pi_{0}\rightarrow\pi_{+}+\pi_{-}}+|\mathcal{M}|^{2}_{a+\pi_{-}\rightarrow\pi_{-}+\pi_{0}}+|\mathcal{M}|^{2}_{a+\pi_{+}\rightarrow\pi_{+}+\pi_{0}}$
(50)
$\displaystyle=\left(\frac{3}{2}\frac{A}{f_{a}f_{\pi}}\frac{1}{1-r^{2}}\right)^{2}\left[(s-m_{\pi}^{2})^{2}+(t-m_{\pi}^{2})^{2}+(u-m_{\pi}^{2})^{2}\right]$
$\displaystyle=\left(\frac{3}{2}\frac{A}{f_{a}f_{\pi}}\frac{1}{1-r^{2}}\right)^{2}\left[s^{2}+t^{2}+u^{2}-3m_{\pi}^{4}-2m_{a}^{2}m_{\pi}^{2}\right],$

where $t\equiv(p_{a}-p_{2})^{2},\,u\equiv(p_{a}-p_{3})^{2}$ and
$\pi_{i},\pi_{j},\pi_{k}$ in the subscript of (50) refer to three pions of
different charge $\pi_{0}$, $\pi_{+}$, $\pi_{-}$. Likewise, the interaction
between an axion and three neutral pions generates the squared scattering
amplitude

$\displaystyle|\mathcal{M}|^{2}_{a+\pi_{0}\rightarrow\pi_{0}+\pi_{0}}=\left(\frac{3}{2}\frac{A}{f_{a}f_{\pi}}\frac{1}{1-r^{2}}\right)^{2}m_{a}^{4}.$
(51)

Inserting the sum of the squared scattering matrix elements, (50) and (51),
into Eq. (25) and integrating over the phase spaces of the three pions gives
the product of the axion-pion collision term and $f_{a,\rm eq}-f_{a}$

$\displaystyle C_{\pi}(f_{a,\rm eq}-f_{a})=\frac{1}{2E_{a}}\int\Bigl{[}$
$\displaystyle
d\Pi_{1}\,d\Pi_{2}\,d\Pi_{3}(|\mathcal{M}|^{2}_{a+\pi_{i}\rightarrow\pi_{j}^{\dagger}+\pi_{k}^{\dagger}}+\frac{1}{2!}|\mathcal{M}|^{2}_{a+\pi_{0}\rightarrow\pi_{0}+\pi_{0}})$
$\displaystyle\times\Lambda(2\pi)^{4}\delta^{4}(p_{a}+p_{1}-p_{2}-p_{3})\Bigr{]},$
(52)

Figure 15: Feynman diagrams associated with axion-pion scattering involving
charged pions.

where

$\displaystyle\Lambda$
$\displaystyle=\left(1+f_{a}\right)\left(1+f_{1}\right)f_{2}f_{3}-f_{a}f_{1}\left(1+f_{2}\right)\left(1+f_{3}\right)$
(53) $\displaystyle\simeq(f_{a,\rm
eq}-f_{a})f_{1}\left(1+f_{2}+f_{3}\right)\,,$ (54)

and $f_{i}$ is the distribution function of particle $i$ possessing momentum
$p_{i}$ in accordance with Fig. 15. In going from (53) to (54), 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$:

$\displaystyle\sum_{i=\gamma,e,\mu,\pi,\delta_{\rm
QCD}}\rho_{i}=\rho_{\gamma}(T)+\rho_{e}(T)+\rho_{\mu}(T)+\rho_{\pi}(T)+\rho_{\delta\rm
QCD}(T)$ (55)

where the left-hand side of (55) is solved from Eq. (33) and the right-hand
side is calculated from Eq. 34 for $\gamma,\mu,e,\pi$ and from Saikawa:2018rcs
for $\rho_{\delta_{\rm QCD}}$, as described more in Appendix C.

To calculate $C_{\pi}$, we first introduce another $\delta$-function in (A) by
writing
$d\Pi_{3}(2\pi)^{3}=\frac{d^{3}\mathbf{p}_{3}}{2E_{3}}=d^{4}p_{3}\delta(p_{4}^{2}-m_{4}^{2})\Theta(p_{4}^{0})$.
By integrating $p_{3}$ over the other delta function
$\delta^{4}(p_{1}+p_{2}-p_{3}-p_{4})$, Eq. (A) simplifies to

$\displaystyle C_{\pi}(f_{a,\rm eq}-f_{a})=\frac{1}{2E_{a}}\int\Bigl{[}$
$\displaystyle
d\Pi_{1}\,d\Pi_{2}(|\mathcal{M}|^{2}_{a+\pi_{i}\rightarrow\pi_{j}^{\dagger}+\pi_{k}^{\dagger}}+\frac{1}{2!}|\mathcal{M}|^{2}_{a+\pi_{0}\rightarrow\pi_{0}+\pi_{0}})\Lambda
2\pi\delta(p_{4}^{2}-m_{\pi 4}^{2})\Theta(p_{4}^{0})\Bigr{]}\,,$ (56)

with the understanding that $p_{3}=p_{a}+p_{1}-p_{2}$. Note that the argument
of the remaining delta-function, $p_{3}^{2}-m_{3}^{2}$, can be written as

$\displaystyle
Q+2\left((E_{a}E_{1}-|\mathbf{p}_{a}||\mathbf{p}_{1}|\cos\alpha)-(E_{1}E_{2}-|\mathbf{p}_{1}||\mathbf{p}_{2}|\cos\gamma)-(E_{a}E_{2}-|\mathbf{p}_{a}||\mathbf{p}_{2}|\cos\theta)\right)$
(57)

where, in the notation of Hannestad:1995rs ,
$Q=m_{3}^{2}-m_{a}^{2}+m_{1}^{2}+m_{2}^{2}$, and $\alpha,\theta,\gamma$ are
the angles between $\mathbf{p}_{a}$ and $\mathbf{p}_{1}$, $\mathbf{p}_{a}$ and
$\mathbf{p}_{2}$, and $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$, respectively. It
is convenient to express the latter angle in terms of the former two by
$\cos\gamma=\cos\alpha\cos\theta+\sin\alpha\sin\theta\cos\beta$.

In the massless axion limit, the argument of the remaining delta function can
easily be expressed in terms of $E_{a}$ as done in DiLuzio:2021vjd . 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 Hannestad:1995rs , which we follow. In particular, the choice of the
angle $\beta$ is most convenient as it only occurs once in (57). The resulting
integral for $C_{\pi}$ is

$\displaystyle C_{\pi}(f_{a,\rm
eq}-f_{a})=\frac{1}{2E_{a}}\int\Bigl{[}\left(\frac{d|\mathbf{p}_{1}|}{(2\pi)^{3}}\frac{\mathbf{p}_{1}^{2}}{2E_{1}}d\cos\alpha\,d\beta\right)\left(\frac{d|\mathbf{p}_{2}|}{(2\pi)^{3}}\frac{\mathbf{p}_{2}^{2}}{2E_{2}}d\cos\theta\,d\phi\right)\left(\frac{2\pi\delta(g(\beta))}{\left|\frac{dg}{d\beta}\right|}\right)$
$\displaystyle\times(|\mathcal{M}|^{2}_{a+\pi_{i}\rightarrow\pi_{j}^{\dagger}+\pi_{k}^{\dagger}}+\frac{1}{2!}|\mathcal{M}|^{2}_{a+\pi_{0}\rightarrow\pi_{0}+\pi_{0}})\Lambda
2\pi\delta(p_{4}^{2}-m_{\pi 4}^{2})\Theta(p_{4}^{0})\Bigr{]}\ $ (58)

where $g(\beta)=p_{3}(\beta)^{2}-m_{3}^{2}$. The integrals over the azimuthal
angles $\beta$ and $\phi$ can be done analytically due to the delta function
and the lack of $\phi$ dependence in the integrand. The remaining integrals
over the pion 3-momenta $|\mathbf{p}_{1}|\in[0,\infty)$ and
$|\mathbf{p}_{2}|\in[0,\infty)$ and the polar angles $\cos\alpha\in[-1,1]$ and
$\cos\theta\in[-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 $\Theta(1-\cos^{2}\beta_{i})$ in (A),
where $\beta_{i}$ is the location of the two (equal and opposite) roots of
$g(\beta)$ Hannestad:1995rs . Note that $\beta_{i}$ are functions of the other
four integration variables.

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. (18),

$\displaystyle\Gamma_{a\pi\leftrightarrow\pi\pi}=\frac{1}{n_{a,\rm
eq}}\int\frac{d^{3}\mathbf{p}_{a}}{(2\pi)^{3}}\mathcal{C}_{\pi}f_{a,\rm
eq}\,\equiv\frac{T^{5}}{f_{a}^{2}f_{\pi}^{2}}\frac{A^{2}}{(1-r^{2})^{2}}\mathcal{F}_{\pi}(m_{a},T),$
(59)

where $n_{a,\rm eq}$ is the thermal number density of axions with mass $m_{a}$
at temperature $T$. As before, $A\equiv\frac{1}{3}(1-z)/(1+z)$, $z\equiv
m_{u}/m_{d}$, and $r\equiv m_{a}/m_{\pi}$.

Figure 16: Numerical evaluation of $\mathcal{F}_{\pi}$ as a function of
$m_{\pi}/T$ for a variety of axion masses. $\mathcal{F}_{\pi}$ becomes
suppressed for $m_{a}>m_{\pi}$ and $T<m_{\pi}$ 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 Hannestad:2005df ;
DiLuzio:2021vjd , corresponds to the $m_{a}/m_{\pi}=0$ contour.

Fig. 16 shows $\mathcal{F}_{\pi}$ as a function of $m_{\pi}/T$ for a variety
of axion masses. For $m_{a}\ll m_{\pi}$, $\mathcal{F}_{\pi}$ reduces to
previous results in the literature for massless axions Hannestad_2005 ;
DiLuzio:2021vjd , with $\mathcal{F}_{\pi}$ related to the $h_{\rm LO}$
function defined in Hannestad_2005 ; DiLuzio:2021vjd  by the mapping
$\mathcal{F}_{\pi}(m_{\pi}/T,m_{a}=0)\equiv 0.212\,h_{\rm LO}(m_{\pi}/T)$.
According to Fig. 16, for $T<m_{\pi}$, $\mathcal{F}_{\pi}$ 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.

##  Appendix B $N_{\rm eff}$ for Large Axion Masses

In the main text, we show the $N_{\rm 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. 7, 7, 12 and 12, the $N_{\rm eff}$ contours become
approximately horizontal above $m_{a}\approx 1$ GeV, signifying that $N_{\rm
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 $\tau$, 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_{\rm 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 $\rho_{\rm max}$ be the maximum energy density that can be
deposited at a certain time $t_{*}$ so that $N_{\rm eff}$ does not drop below
an arbitrary contour, $N_{\rm eff,0}$, which we will take to be $\simeq 2.62$,
the $2\sigma$ limit on $N_{\rm eff}$ allowed by Planck. Choose a point
$(m_{0},{\tau_{0}})$ that lies on this $N_{\rm eff,0}$ contour in the
$m_{a}\gtrsim 1$ GeV region. Analytically, the energy density of this non-
relativistic axion at time $t_{*}$ is

$\displaystyle\rho(t_{*})=\rho_{\rm init}\left(\frac{a_{\rm
init}}{a(t_{*})}\right)^{3}e^{-t_{*}/\tau_{0}}=m_{0}Y_{0}s(t_{*})e^{-t_{*}/\tau_{0}}\equiv\rho_{\rm
max}$ (60)

where $Y_{0}$ is the axion yield, and $s(t_{*})$ the entropy density at time
$t_{*}$. Note Eq. (60) defines $\rho_{\rm max}$. It follows that for axions of
different $(m_{a},\tau)$ to possess the same energy density as $\rho_{\rm
max}$ at time $t_{*}$, requires

$\displaystyle\frac{m_{0}Y_{0}}{m_{a}Y_{a}}=\exp
t_{*}\left(\frac{1}{\tau_{0}}-\frac{1}{\tau}\right)\,$ (61)

or equivalently,

$\displaystyle\tau=\frac{\tau_{0}}{1+\frac{\tau_{0}}{t_{*}}\ln\left(\frac{m_{a}Y_{a}}{m_{0}Y_{0}}\right)},$
(62)

We perform a numerical fit of the $N_{\rm eff,0}=2.62$ contour with the anchor
point $(m_{0},\tau_{0})=(2.0\,{\rm GeV},\,4.3\times 10^{-2}\,{\rm s})$ and
find $t_{*}\simeq 0.17\,{\rm s}$. Other anchor points give similar $t_{*}$.
Eq. (38) follows from this fit.

##  Appendix 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, (24), into a partition of
$N$ ordinary differential equations, $\\{f_{\mathbf{\tilde{p},i}}(t)\\}$ of
time, with each ODE corresponding to the time evolution of the phase space
density at a fixed comoving momentum,
$|\tilde{\mathbf{p}}_{i}|=|\mathbf{p}_{i}|a(t)$, with $i\in\\{1,...,N\\}$. In
our numerical setup, we split (24) into $N=24$ ODEs of logarithmically
equidistant $|\tilde{\mathbf{p}}_{i}|$ where $i=1$ corresponds to the fixed
comoving momentum $|\tilde{\mathbf{p}_{1}}|=e^{-8}T_{\chi\rm PT}$ and $i=N$
corresponds to $|\tilde{\mathbf{p}}_{N}|=e^{7/2}T_{\chi\rm 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, $|\mathbf{p}|/T$. We have verified the convergence of our results
by checking that $N_{\rm eff}$ changes by less than 1% when using larger $N$.
In addition, our results for $N_{\rm 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_{\rm eff}^{\rm SM}=3.044$ Bennett:2020zkv ; Froustey:2020mcq
; Zyla:2020zbs  by $\sim 0.2\%$ because we do not include effects from QED
corrections or neutrino oscillations.

The dynamical timescale ($\sim 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_{\chi\rm PT})$, where $t_{\chi\rm
PT}=1/2H(T_{\chi\rm PT})$ is the starting time of the Boltzmann code.

Last, we determine the extra QCD contributions to the energy density
($\rho_{\delta_{\rm QCD}}$) and pressure ($P_{\delta_{\rm QCD}}$) arising from
heavy mesons and from ideal gas law deviations using the calculations of
Saikawa:2018rcs , which tabulated $g_{*}(T)$ and $g_{*S}(T)$ for the Standard
Model across the QCD phase transition.

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 Saikawa:2018rcs . Top
right: $\rho_{\delta\rm 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}\equiv P/T^{4}$.

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\sim 100$ MeV, is the
Standard Model value of $g_{*}$ from Saikawa:2018rcs , which includes
contributions, $\delta_{\rm QCD}$, from heavy mesons and from ideal gas law
deviations. We extract $\rho_{\delta_{\rm QCD}}$ by taking the difference
between the blue and orange contours, $\Delta g_{*}$, and multiplying this
result by $\pi^{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\sim 100$ MeV, is the Standard Model value of $g_{*P}$ from
Saikawa:2018rcs , which we infer through the relationship
$g_{*P}=\frac{\pi^{2}}{30}(\frac{4}{3}g_{*S}-g_{*})$. We extract
$P_{\delta_{\rm QCD}}$ by taking the difference between the blue and orange
contours, $\Delta g_{*P}$, and multiplying this result by $T^{4}$, as shown by
the bottom right panel.

## References

  * (1) A. E. Nelson, _Naturally Weak CP Violation_ , _Phys. Lett. B_ 136 (1984) 387.
  * (2) S. M. Barr, _Solving the Strong CP Problem Without the Peccei-Quinn Symmetry_ , _Phys. Rev. Lett._ 53 (1984) 329.
  * (3) M. A. B. Beg and H. S. Tsao, _Strong P, T Noninvariances in a Superweak Theory_ , _Phys. Rev. Lett._ 41 (1978) 278.
  * (4) R. N. Mohapatra and G. Senjanovic, _Natural Suppression of Strong p and t Noninvariance_ , _Phys. Lett. B_ 79 (1978) 283.
  * (5) K. S. Babu and R. N. Mohapatra, _A Solution to the Strong CP Problem Without an Axion_ , _Phys. Rev. D_ 41 (1990) 1286.
  * (6) R. D. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * (7) R. D. Peccei and H. R. Quinn, _Constraints Imposed by CP Conservation in the Presence of Instantons_ , _Phys. Rev._ D16 (1977) 1791.
  * (8) H. M. Georgi, L. J. Hall and M. B. Wise, _Grand Unified Models With an Automatic Peccei-Quinn Symmetry_ , _Nucl. Phys. B_ 192 (1981) 409.
  * (9) D. Harlow and H. Ooguri, _Symmetries in quantum field theory and quantum gravity_ , _Commun. Math. Phys._ 383 (2021) 1669 [1810.05338].
  * (10) T. Banks and N. Seiberg, _Symmetries and Strings in Field Theory and Gravity_ , _Phys. Rev. D_ 83 (2011) 084019 [1011.5120].
  * (11) R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, R. Watkins and L. M. Widrow, _Solutions to the strong CP problem in a world with gravity_ , _Phys. Lett. B_ 282 (1992) 132 [hep-ph/9203206].
  * (12) S. M. Barr and D. Seckel, _Planck scale corrections to axion models_ , _Phys. Rev. D_ 46 (1992) 539.
  * (13) M. Kamionkowski and J. March-Russell, _Planck scale physics and the Peccei-Quinn mechanism_ , _Phys. Lett. B_ 282 (1992) 137 [hep-th/9202003].
  * (14) M. Dine, _Problems of naturalness: Some lessons from string theory_ , in _Conference on Topics in Quantum Gravity_ , 7, 1992, hep-th/9207045.
  * (15) K.-w. Choi, D. B. Kaplan and A. E. Nelson, _Is CP a gauge symmetry?_ , _Nucl. Phys. B_ 391 (1993) 515 [hep-ph/9205202].
  * (16) M. Dine, R. G. Leigh and D. A. MacIntire, _Of CP and other gauge symmetries in string theory_ , _Phys. Rev. Lett._ 69 (1992) 2030 [hep-th/9205011].
  * (17) K. Choi and J. E. Kim, _DYNAMICAL AXION_ , _Phys. Rev. D_ 32 (1985) 1828.
  * (18) D. Cadamuro and J. Redondo, _Cosmological bounds on pseudo Nambu-Goldstone bosons_ , _JCAP_ 02 (2012) 032 [1110.2895].
  * (19) Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * (20) CMB-S4 collaboration, _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * (21) S. Chang and K. Choi, _Hadronic axion window and the big bang nucleosynthesis_ , _Phys. Lett. B_ 316 (1993) 51 [hep-ph/9306216].
  * (22) R. Z. Ferreira and A. Notari, _Observable Windows for the QCD Axion Through the Number of Relativistic Species_ , _Phys. Rev. Lett._ 120 (2018) 191301 [1801.06090].
  * (23) F. Arias-Aragón, F. D’eramo, R. Z. Ferreira, L. Merlo and A. Notari, _Cosmic Imprints of XENON1T Axions_ , _JCAP_ 11 (2020) 025 [2007.06579].
  * (24) F. Arias-Aragón, F. D’Eramo, R. Z. Ferreira, L. Merlo and A. Notari, _Production of Thermal Axions across the ElectroWeak Phase Transition_ , _JCAP_ 03 (2021) 090 [2012.04736].
  * (25) R. Z. Ferreira, A. Notari and F. Rompineve, _Dine-Fischler-Srednicki-Zhitnitsky axion in the CMB_ , _Phys. Rev. D_ 103 (2021) 063524 [2012.06566].
  * (26) F. D’Eramo, F. Hajkarim and S. Yun, _Thermal axion production at low temperatures: a smooth treatment of the QCD phase transition_ , 2108.04259.
  * (27) F. D’Eramo, F. Hajkarim and S. Yun, _Thermal QCD Axions across Thresholds_ , 2108.05371.
  * (28) F. D’Eramo and S. Yun, _Flavor violating axions in the early Universe_ , _Phys. Rev. D_ 105 (2022) 075002 [2111.12108].
  * (29) S. Hannestad and G. Raffelt, _Cosmological mass limits on neutrinos, axions, and other light particles_ , _JCAP_ 04 (2004) 008 [hep-ph/0312154].
  * (30) W. Giarè, E. Di Valentino, A. Melchiorri and O. Mena, _New cosmological bounds on hot relics: axions and neutrinos_ , _Mon. Not. Roy. Astron. Soc._ 505 (2021) 2703 [2011.14704].
  * (31) L. Caloni, M. Gerbino, M. Lattanzi and L. Visinelli, _Novel cosmological bounds on thermally-produced axion-like particles_ , 2205.01637.
  * (32) S. Dimopoulos, _A Solution of the Strong CP Problem in Models With Scalars_ , _Phys. Lett. B_ 84 (1979) 435.
  * (33) S. H. H. Tye, _A Superstrong Force With a Heavy Axion_ , _Phys. Rev. Lett._ 47 (1981) 1035.
  * (34) B. Holdom and M. E. Peskin, _Raising the Axion Mass_ , _Nucl. Phys. B_ 208 (1982) 397.
  * (35) B. Holdom, _Strong QCD at High-energies and a Heavy Axion_ , _Phys. Lett. B_ 154 (1985) 316.
  * (36) J. M. Flynn and L. Randall, _A Computation of the Small Instanton Contribution to the Axion Potential_ , _Nucl. Phys. B_ 293 (1987) 731.
  * (37) P. Agrawal and K. Howe, _Factoring the Strong CP Problem_ , _JHEP_ 12 (2018) 029 [1710.04213].
  * (38) T. Gherghetta, V. V. Khoze, A. Pomarol and Y. Shirman, _The Axion Mass from 5D Small Instantons_ , _JHEP_ 03 (2020) 063 [2001.05610].
  * (39) V. A. Rubakov, _Grand unification and heavy axion_ , _JETP Lett._ 65 (1997) 621 [hep-ph/9703409].
  * (40) Z. Berezhiani, L. Gianfagna and M. Giannotti, _Strong CP problem and mirror world: The Weinberg-Wilczek axion revisited_ , _Phys. Lett. B_ 500 (2001) 286 [hep-ph/0009290].
  * (41) S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * (42) F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * (43) H. Fukuda, K. Harigaya, M. Ibe and T. T. Yanagida, _Model of visible QCD axion_ , _Phys. Rev. D_ 92 (2015) 015021 [1504.06084].
  * (44) A. Hook, S. Kumar, Z. Liu and R. Sundrum, _High Quality QCD Axion and the LHC_ , _Phys. Rev. Lett._ 124 (2020) 221801 [1911.12364].
  * (45) K. J. Kelly, S. Kumar and Z. Liu, _Heavy axion opportunities at the DUNE near detector_ , _Phys. Rev. D_ 103 (2021) 095002 [2011.05995].
  * (46) J. E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , _Phys. Rev. Lett._ 43 (1979) 103.
  * (47) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493.
  * (48) M. Millea, L. Knox and B. Fields, _New Bounds for Axions and Axion-Like Particles with keV-GeV Masses_ , _Phys. Rev. D_ 92 (2015) 023010 [1501.04097].
  * (49) P. F. Depta, M. Hufnagel and K. Schmidt-Hoberg, _Robust cosmological constraints on axion-like particles_ , _JCAP_ 05 (2020) 009 [2002.08370].
  * (50) D. Aloni, Y. Soreq and M. Williams, _Coupling QCD-Scale Axionlike Particles to Gluons_ , _Phys. Rev. Lett._ 123 (2019) 031803 [1811.03474].
  * (51) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (52) A. Salvio, A. Strumia and W. Xue, _Thermal axion production_ , _JCAP_ 01 (2014) 011 [1310.6982].
  * (53) L. Di Luzio, G. Martinelli and G. Piazza, _Breakdown of chiral perturbation theory for the axion hot dark matter bound_ , _Phys. Rev. Lett._ 126 (2021) 241801 [2101.10330].
  * (54) P. Graf and F. D. Steffen, _Thermal axion production in the primordial quark-gluon plasma_ , _Phys. Rev. D_ 83 (2011) 075011 [1008.4528].
  * (55) D. Cadamuro, S. Hannestad, G. Raffelt and J. Redondo, _Cosmological bounds on sub-MeV mass axions_ , _JCAP_ 02 (2011) 003 [1011.3694].
  * (56) H. A. Weldon, _Covariant calculations at finite temperature: the relativistic plasma_ , _Physical Review D_ 26 (1982) 1394.
  * (57) S. Hannestad, A. Mirizzi and G. Raffelt, _A new cosmological mass limit on thermal relic axions_ , _Journal of Cosmology and Astroparticle Physics_ 2005 (2005) 002.
  * (58) K. Saikawa and S. Shirai, _Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model_ , _JCAP_ 05 (2018) 035 [1803.01038].
  * (59) “The Numerical Method of Lines.” https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html.
  * (60) F. Ertas and F. Kahlhoefer, _On the interplay between astrophysical and laboratory probes of MeV-scale axion-like particles_ , _JHEP_ 07 (2020) 050 [2004.01193].
  * (61) B. Döbrich, J. Jaeckel, F. Kahlhoefer, A. Ringwald and K. Schmidt-Hoberg, _ALPtraum: ALP production in proton beam dump experiments_ , _JHEP_ 02 (2016) 018 [1512.03069].
  * (62) M. J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer and K. Schmidt-Hoberg, _Revised constraints and Belle II sensitivity for visible and invisible axion-like particles_ , _JHEP_ 12 (2017) 094 [1709.00009].
  * (63) NA64 collaboration, _Search for Axionlike and Scalar Particles with the NA64 Experiment_ , _Phys. Rev. Lett._ 125 (2020) 081801 [2005.02710].
  * (64) FASER collaboration, _FASER’s physics reach for long-lived particles_ , _Phys. Rev. D_ 99 (2019) 095011 [1811.12522].
  * (65) S. Gori, G. Perez and K. Tobioka, _KOTO vs. NA62 Dark Scalar Searches_ , _JHEP_ 08 (2020) 110 [2005.05170].
  * (66) A. Mariotti, D. Redigolo, F. Sala and K. Tobioka, _New LHC bound on low-mass diphoton resonances_ , _Phys. Lett. B_ 783 (2018) 13 [1710.01743].
  * (67) S. Chakraborty, M. Kraus, V. Loladze, T. Okui and K. Tobioka, _Heavy QCD axion in b→s transition: Enhanced limits and projections_ , _Phys. Rev. D_ 104 (2021) 055036 [2102.04474].
  * (68) E. Bertholet, S. Chakraborty, V. Loladze, T. Okui, A. Soffer and K. Tobioka, _Heavy QCD axion at Belle II: Displaced and prompt signals_ , _Phys. Rev. D_ 105 (2022) L071701 [2108.10331].
  * (69) D. Dunsky, L. J. Hall and K. Harigaya, _Dark Matter, Dark Radiation and Gravitational Waves from Mirror Higgs Parity_ , _JHEP_ 02 (2020) 078 [1908.02756].
  * (70) D. Dunsky, L. J. Hall and K. Harigaya, _CHAMP Cosmic Rays_ , _JCAP_ 07 (2019) 015 [1812.11116].
  * (71) W. L. Xu, J. B. Muñoz and C. Dvorkin, _Cosmological Constraints on Light (but Massive) Relics_ , 2107.09664.
  * (72) C. O’HARE, _cajohare/axionlimits: Axionlimits_ , July, 2020.  10.5281/zenodo.3932430.
  * (73) J. Preskill, M. B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127.
  * (74) L. F. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133.
  * (75) M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * (76) J. H. Chang, R. Essig and S. D. McDermott, _Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle_ , _JHEP_ 09 (2018) 051 [1803.00993].
  * (77) A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi and O. Straniero, _Revisiting the bound on axion-photon coupling from Globular Clusters_ , _Phys. Rev. Lett._ 113 (2014) 191302 [1406.6053].
  * (78) CAST collaboration, _An Improved limit on the axion-photon coupling from the CAST experiment_ , _JCAP_ 04 (2007) 010 [hep-ex/0702006].
  * (79) CAST collaboration, _New CAST Limit on the Axion-Photon Interaction_ , _Nature Phys._ 13 (2017) 584 [1705.02290].
  * (80) N. Vinyoles, A. Serenelli, F. L. Villante, S. Basu, J. Redondo and J. Isern, _New axion and hidden photon constraints from a solar data global fit_ , _JCAP_ 10 (2015) 015 [1501.01639].
  * (81) M. Kawasaki, K. Kohri and N. Sugiyama, _Cosmological constraints on late time entropy production_ , _Phys. Rev. Lett._ 82 (1999) 4168 [astro-ph/9811437].
  * (82) M. Kawasaki, K. Kohri and N. Sugiyama, _MeV scale reheating temperature and thermalization of neutrino background_ , _Phys. Rev. D_ 62 (2000) 023506 [astro-ph/0002127].
  * (83) T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram and S. Hannestad, _MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative and hadronic decays of massive particles_ , _JCAP_ 12 (2019) 012 [1908.10189].
  * (84) S. Hannestad and J. Madsen, _Neutrino decoupling in the early universe_ , _Phys. Rev. D_ 52 (1995) 1764 [astro-ph/9506015].
  * (85) S. Hannestad, A. Mirizzi and G. Raffelt, _New cosmological mass limit on thermal relic axions_ , _JCAP_ 07 (2005) 002 [hep-ph/0504059].
  * (86) J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor et al., _Towards a precision calculation of $N_{\rm eff}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED_, _JCAP_ 04 (2021) 073 [2012.02726].
  * (87) J. Froustey, C. Pitrou and M. C. Volpe, _Neutrino decoupling including flavour oscillations and primordial nucleosynthesis_ , _JCAP_ 12 (2020) 015 [2008.01074].
  * (88) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.