# Axion hot dark matter bound, reliably

Luca Di Luzio<EMAIL_ADDRESS>DESY, Notkestraße 85, D-22607 Hamburg,
Germany Guido Martinelli<EMAIL_ADDRESS>Physics Department
and INFN Sezione di Roma La Sapienza, Piazzale Aldo Moro 5, 00185 Roma, Italy
Gioacchino Piazza<EMAIL_ADDRESS>IJCLab, Pôle Théorie
(Bât. 210), CNRS/IN2P3 et Université Paris-Saclay, 91405 Orsay, France

###### Abstract

We show that the commonly adopted hot dark matter (HDM) bound on the axion
mass $m_{a}\lesssim$ 1 eV is not reliable, since it is obtained by
extrapolating the chiral expansion in a region where the effective field
theory breaks down. This is explicitly shown via the calculation of the axion-
pion thermalization rate at the next-to-leading order in chiral perturbation
theory. We finally advocate a strategy for a sound extraction of the axion HDM
bound via lattice QCD techniques.

Introduction. The axion originally emerged as a low-energy remnant of the
Peccei Quinn solution to the strong CP problem Peccei and Quinn (1977a, b);
Wilczek (1978); Weinberg (1978), but it also unavoidably contributes to the
energy density of the Universe. There are two qualitatively different
populations of relic axions, a non-thermal one comprising cold dark matter
(DM) Preskill  _et al._ (1983); Abbott and Sikivie (1983); Dine and Fischler
(1983); Davis (1986), and a thermal axion population Turner (1987) which,
while still relativistic, would behave as extra dark radiation. Such hot dark
matter (HDM) component contributes to the effective number of extra
relativistic degrees of freedom Kolb and Turner (1990) $\Delta N_{\rm
eff}\simeq 4/7\left(43/[4g_{S}(T_{D})]\right)^{4/3}$, with $g_{S}(T_{D})$ the
number of entropy degrees of freedom at the axion decoupling temperature,
$T_{D}$. The value of $\Delta N_{\rm eff}$ is constrained by cosmic microwave
background (CMB) experiments, such as the Planck satellite Aghanim  _et al._
(2020a, b), while planned CMB Stage 4 (CMB-S4) experiments Abazajian  _et al._
(2016) will provide an observable window on the axion interactions.

There are several processes that can keep the axion in thermal equilibrium
with the Standard Model (SM) thermal bath. From the standpoint of the axion
solution to the strong CP problem, an unavoidable process arises from the
model-independent coupling to gluons,
$\frac{\alpha_{s}}{8\pi}\frac{a}{f_{a}}G\tilde{G}$.111Other thermalization
channels arise from model-dependent axion couplings to photons Turner (1987),
SM quarks Salvio  _et al._ (2014); Baumann  _et al._ (2016); Ferreira and
Notari (2018); Arias-Aragon  _et al._ (2020) and leptons D’Eramo  _et al._
(2018). For $T_{D}\gtrsim 1$ GeV thermal axion production proceeds via its
scatterings with gluons in the quark-gluon plasma Masso  _et al._ (2002); Graf
and Steffen (2011), while for $T_{D}\lesssim 1$ GeV processes involving pions
and nucleons must be considered Berezhiani  _et al._ (1992); Chang and Choi
(1993); Hannestad  _et al._ (2005). The latter, have the advantage of
occurring very late in the thermal history, so that it is unlikely that the
corresponding population of thermal axions could be diluted by inflation. The
transition between the two regimes depends on the strength of the axion
interactions set by $f_{a}$ or, equivalently, by $m_{a}\simeq
5.7\times(10^{6}\ \text{GeV}/f_{a})$ eV, and it encompasses the range
$m_{a}\in[0.01,0.1]$ eV (with heavier axions leading to lower decoupling
temperatures). Although the transition region cannot be precisely determined
due to the complications of the quark-hadron phase transition, for heavier
axions approaching the eV scale the main thermalization channel is
$a\pi\leftrightarrow\pi\pi$ Chang and Choi (1993); Hannestad  _et al._ (2005),
with $T_{D}\lesssim 200$ MeV. In this regime, scatterings off nucleons are
subdominant because of the exponential suppression in their number density.

The highest attainable axion mass from cosmological constraints on extra
relativistic degrees of freedom, also known as HDM bound, translate into
$m_{a}\lesssim$ 1 eV Zyla  _et al._ (2020). Based on a leading-order (LO)
axion-pion chiral effective field theory (EFT) analysis of the axion-pion
thermalization rate Chang and Choi (1993); Hannestad  _et al._ (2005), the
axion HDM bound has been reconsidered in Refs. Melchiorri  _et al._ (2007);
Hannestad  _et al._ (2008, 2010); Archidiacono  _et al._ (2013); Giusarma  _et
al._ (2014); Di Valentino  _et al._ (2015, 2016); Archidiacono  _et al._
(2015); Giarè  _et al._ (2020), also in correlation with relic neutrinos. The
most recent update Giarè  _et al._ (2020) quotes a 95$\%$ CL bound that ranges
from $m_{a}\lesssim$ 0.2 eV to 1 eV, depending on the used data set and
assumed cosmological model. Although the axion mass range relevant for the HDM
bound is in generic tension with astrophysical constraints, the latter can be
tamed in several respects.222Tree-level axion couplings to electrons are
absent in KSVZ models Kim (1979); Shifman  _et al._ (1980), thus relaxing the
constraints from Red Giants and White Dwarfs. The axion coupling to photons,
constrained by Horizontal Branch stars evolution, can be accidentally
suppressed in certain KSVZ-like models Kaplan (1985); Di Luzio  _et al._
(2017a, b). Finally, the SN1987A bound on the axion coupling to nucleons can
be considered less robust both from the astrophysical and experimental point
of view Raffelt (1990); Chang  _et al._ (2018); Carenza  _et al._ (2019); Bar
_et al._ (2020).

It is the purpose of this Letter to revisit the axion HDM bound in the context
of the next-to-LO (NLO) axion-pion chiral EFT. This is motivated by the simple
observation that the mean energy of pions (axions) in a heat bath of $T\simeq
100$ MeV is $\left\langle E\right\rangle\equiv\rho/n\simeq 350$ MeV ($270$
MeV), thus questioning the validity of the chiral expansion for the scattering
process $a\pi\leftrightarrow\pi\pi$. The latter is expected to fail for
$\sqrt{s}\sim\left\langle E_{\pi}\right\rangle+\left\langle
E_{a}\right\rangle\gtrsim 500$ MeV, corresponding to temperatures well below
that of QCD deconfinement, $T_{c}=154\pm 9$ MeV Bazavov  _et al._ (2012).

In this work, we provide for the first time the formulation of the full axion-
pion Lagrangian at NLO, including also derivative axion couplings to the
pionic current (previous NLO studies only considered non-derivative axion-pion
interactions Spalinski (1988); Grilli di Cortona  _et al._ (2016)), and paying
special attention to the issue of the axion-pion mixing. Next, we perform a
NLO calculation of the $a\pi\leftrightarrow\pi\pi$ thermalization rate (that
can be cast as an expansion in $T/f_{\pi}$, with $f_{\pi}\simeq 92$ MeV) and
show that the NLO correction saturates half of the LO contribution for
$T_{\chi}\simeq 62$ MeV. The latter can be considered as the maximal
temperature above which the chiral description breaks down for the process
under consideration. On the other hand, the region from $T_{\chi}$ up to
$T_{c}$, where chiral perturbation theory cannot be applied, turns out to be
crucial for the extraction of the HDM bound and for assessing the sensitivity
of future CMB experiments.

We conclude with a proposal for extracting the axion-pion thermalization rate
via a direct Lattice QCD calculation, in analogy to the well-studied case of
$\pi$-$\pi$ scattering.

Axion-pion scattering at LO. The construction of the LO axion-pion Lagrangian
was discussed long ago in Refs. Di Vecchia and Veneziano (1980); Georgi  _et
al._ (1986). We recall here its basic ingredients (see also Chang and Choi
(1993); Di Luzio  _et al._ (2020)), in view of the extension at NLO. Defining
the pion Goldstone matrix $U=e^{i\pi^{A}\sigma^{A}/f_{\pi}}$, with
$f_{\pi}\simeq 92$ MeV, $\pi^{A}$ and $\sigma^{A}$ ($A=1,2,3$) denoting
respectively the real pion fields and the Pauli matrices, the LO axion-pion
interactions stem from

$\mathscr{L}^{\rm LO}_{a\text{-}\pi}=\frac{f_{\pi}^{2}}{4}{\rm
Tr}\left[U\chi^{\dagger}_{a}+\chi_{a}U^{\dagger}\right]+\frac{\partial^{\mu}a}{2f_{a}}{\rm
Tr}\left[c_{q}\sigma^{A}\right]J^{A}_{\mu}\,,$ (1)

where $\chi_{a}=2B_{0}M_{a}$, in terms of the quark condensate $B_{0}$ and the
‘axion-dressed’ quark mass matrix
$M_{a}=e^{i\frac{a}{2f_{a}}Q_{a}}M_{q}e^{i\frac{a}{2f_{a}}Q_{a}}$, with
$M_{q}=\mbox{diag}\,(m_{u},m_{d})$ and $\mbox{Tr}\,Q_{a}=1$. The latter
condition ensures that the axion field is transferred from the operator
$\frac{\alpha_{s}}{8\pi}\frac{a}{f_{a}}G\tilde{G}$ to the phase of the quark
mass matrix, via the quark axial field redefinition
$q\to\exp(i\gamma_{5}\frac{a}{2f_{a}}Q_{a})q$. In the following, we set
$Q_{a}=M_{q}^{-1}/\mbox{Tr}\,M_{q}^{-1}$, so that terms linear in $a$
(including $a$-$\pi^{0}$ mass mixing) drop out from the first term in Eq. (1).
Hence, in this basis, the only linear axion interaction is the derivative one
with the conserved ${\rm SU}(2)_{A}$ pion current. The latter reads at LO

$J^{A}_{\mu}|^{\rm LO}=\frac{i}{4}f_{\pi}^{2}{\rm
Tr}\left[\sigma^{A}\left(U\partial_{\mu}U^{\dagger}-U^{\dagger}\partial_{\mu}U\right)\right]\,,$
(2)

while the derivative axion coupling in Eq. (1) is
$\mbox{Tr}\,\left[c_{q}\sigma^{A}\right]=(\frac{m_{u}-m_{d}}{m_{u}+m_{d}}+c^{0}_{u}-c^{0}_{d})\delta^{A3}$,
where the first term arises from the axial quark rotation that removed the
$aG\tilde{G}$ operator and the second one originates from the model-dependent
coefficient $c^{0}_{q}=\text{diag}(c^{0}_{u},c^{0}_{d})$, defined via the
Lagrangian term
$\frac{\partial^{\mu}a}{2f_{a}}\overline{q}c^{0}_{q}\gamma_{\mu}\gamma_{5}q$.
For instance, $c^{0}_{u,d}=0$ in the KSVZ model Kim (1979); Shifman  _et al._
(1980), while $c^{0}_{u}=\frac{1}{3}\cos^{2}\beta$ and
$c^{0}_{d}=\frac{1}{3}\sin^{2}\beta$ in the DFSZ model Zhitnitsky (1980); Dine
_et al._ (1981), with $\tan\beta$ the ratio between the vacuum expectation
values of two Higgs doublets. Expanding the pion matrix in Eq. (1) one obtains

$\mathscr{L}^{\rm
LO}_{a\text{-}\pi}\supset\epsilon\,\partial^{\mu}a\partial_{\mu}\pi_{0}+\frac{C_{a\pi}}{f_{a}f_{\pi}}\partial^{\mu}a[\partial\pi\pi\pi]_{\mu}\,,$
(3)

with the definitions
$[\partial\pi\pi\pi]_{\mu}=2\partial_{\mu}\pi_{0}\pi_{+}\pi_{-}-\pi_{0}\partial_{\mu}\pi_{+}\pi_{-}-\pi_{0}\pi_{+}\partial_{\mu}\pi_{-}$,
$\epsilon=-\frac{3f_{\pi}C_{a\pi}}{2f_{a}}$ and

$C_{a\pi}=\frac{1}{3}\left(\frac{m_{d}-m_{u}}{m_{u}+m_{d}}+c_{d}^{0}-c_{u}^{0}\right)\,.$
(4)

At the LO in $\epsilon$ the diagonalization of the $a$-$\pi^{0}$ term is
obtained by shifting $a\to a-\epsilon\pi^{0}$ and
$\pi^{0}\to\pi^{0}+\mathcal{O}(\epsilon^{3})a$, where we used the fact that
$m_{a}/m_{\pi}=\mathcal{O}(\epsilon)$. Hence, as long as we are interested in
effects that are linear in $a$ and neglect $\mathcal{O}(\epsilon^{3})$
corrections, the axion-pion interactions in Eq. (3) are already in the basis
with canonical propagators.

For temperatures below the QCD phase transition, the main processes relevant
for the axion thermalization rate are
$a(p_{1})\pi_{0}(p_{2})\rightarrow\pi_{+}(p_{3})\pi_{-}(p_{4})$, whose
amplitude at LO reads

$\mathcal{M}^{\rm
LO}_{a\pi_{0}\rightarrow\pi_{+}\pi_{-}}=\frac{C_{a\pi}}{f_{\pi}f_{a}}\frac{3}{2}\left[m_{\pi}^{2}-s\right]\,,$
(5)

with $s=(p_{1}+p_{2})^{2}$, together with the crossed channels
$a\pi_{-}\rightarrow\pi_{0}\pi_{-}$ and $a\pi_{+}\rightarrow\pi_{+}\pi_{0}$.
The amplitudes of the latter are obtained by replacing $s\to
t=(p_{1}-p_{3})^{2}$ and $s\to u=(p_{1}-p_{4})^{2}$, respectively. Taking
equal masses for the neutral and charged pions, one finds the squared matrix
element (summed over the three channels above) Hannestad  _et al._ (2005)

$\sum|\mathcal{M}|_{\rm
LO}^{2}=\left(\frac{C_{a\pi}}{f_{a}f_{\pi}}\right)^{2}\frac{9}{4}\left[s^{2}+t^{2}+u^{2}-3m_{\pi}^{4}\right]\,.$
(6)

Axion-pion scattering at NLO. To compute the axion thermalization process
beyond LO we need to consider the one-loop amplitudes from the LO Lagrangian
in Eq. (1) as well as the tree-level amplitudes stemming from the NLO axion-
pion Lagrangian, both contributing to $\mathcal{O}(p^{4})$ in the chiral
expansion. The NLO interactions include the derivative coupling of the axion
to the NLO axial current, which has been computed here for the first time.

We stick to the expression of the NLO chiral Lagrangian given in Ref. Gasser
and Leutwyler (1984) (see for example Appendix D in Scherer (2003) for trace
notation), which, considering only two flavours, depends on $10$ low-energy
constants (LECs) $\ell_{1},\ell_{2},\dots,\ell_{7},h_{1},h_{2},h_{3}$. The
axion field has been included in the phase of the quark mass matrix, as
described after Eq. (1). Note that since we are interested in $2\to 2$
scattering processes, we can neglect the $\mathcal{O}(p^{4})$ Wess-Zumino-
Witten term Wess and Zumino (1971); Witten (1983) since it contains operators
with an odd number of bosons.

To compute the axial current $J^{A}_{\mu}$ at NLO, we promote the ordinary
derivative to a covariant one, defined as
$D_{\mu}U=\partial_{\mu}U-ir_{\mu}U+iUl_{\mu}$, with
$r_{\mu}=r_{\mu}^{A}\sigma^{A}/2$ and $l_{\mu}=l_{\mu}^{A}\sigma^{A}/2$
external fields which can be used to include electromagnetic or weak effects.
The left and right SU(2) currents are obtained by differentiating the NLO
Lagrangian with respect to $l_{\mu}^{A}$ and $r_{\mu}^{A}$, respectively.
Taking the $R-L$ combination and switching off the external fields, the NLO
axial current reads

$\displaystyle J^{A}_{\mu}|^{\rm NLO}=\frac{i}{2}\ell_{1}{\rm
Tr}\left[\sigma^{A}\left\\{\partial_{\mu}U^{\dagger},U\right\\}\right]{\rm
Tr}\left[\partial_{\nu}U\partial^{\nu}U^{\dagger}\right]$
$\displaystyle+\frac{i}{4}\ell_{2}{\rm
Tr}\left[\sigma^{A}\left\\{\partial^{\nu}U^{\dagger},U\right\\}\right]{\rm
Tr}\left[\partial_{\mu}U\partial_{\nu}U^{\dagger}+\partial_{\nu}U\partial_{\mu}U^{\dagger}\right]$
$\displaystyle-\frac{i}{8}\ell_{4}{\rm
Tr}\big{[}\sigma^{A}\left\\{\partial_{\mu}U,\chi^{\dagger}_{a}\right\\}-\sigma^{A}\left\\{U,\partial_{\mu}\chi^{\dagger}_{a}\right\\}$
$\displaystyle\ \ \ \ \ \ \
+\sigma^{A}\left\\{\partial_{\mu}\chi_{a},U^{\dagger}\right\\}-\sigma^{A}\left\\{\chi_{a},\partial_{\mu}U^{\dagger}\right\\}\big{]}\,,$
(7)

where curly brackets indicate anti-commutators.

New $a\text{-}\pi_{0}$ mixings arise at NLO, both at tree level from the NLO
Lagrangian and at one loop from $\mathscr{L}^{\rm LO}_{a\text{-}\pi}$. These
mixings are explicitly taken into account in the Lehmann-Symanzik-Zimmermann
(LSZ) reduction formula Lehmann  _et al._ (1955) (focussing e.g. on the
$a\pi_{0}\to\pi_{+}\pi_{-}$ channel)

$\displaystyle\mathcal{M}_{a\pi_{0}\to\pi_{+}\pi_{-}}$
$\displaystyle=\frac{1}{\sqrt{Z_{a}Z_{\pi}^{3}}}\Pi_{i=1}^{4}\lim_{p_{i}^{2}\to
m_{i}^{2}}\left(p_{i}^{2}-m_{i}^{2}\right)$ $\displaystyle\times
G_{a\pi_{0}\pi_{+}\pi_{-}}(p_{1},p_{2},p_{3},p_{4})\,,$ (8)

where the index $i$ runs over the external particles, $Z_{a}$ ($Z_{\pi}$) is
the wave-function renormalization of the axion (pion) field and the full
4-point Green’s function is given by

$\displaystyle G_{a\pi_{0}\pi_{+}\pi_{-}}=\sum_{i,j=a,\pi_{0}}{\cal
G}_{ij\pi_{+}\pi_{-}}$ (9) $\displaystyle\times
G_{\pi_{+}\pi_{+}}(m^{2}_{\pi})G_{\pi_{-}\pi_{-}}(m^{2}_{\pi})G_{ai}(m^{2}_{a}=0)G_{\pi_{0}j}(m^{2}_{\pi})\,.$

The first term is the amputated 4-point function, multiplied by the 2-point
functions of the external legs with the axion mass to zero. Working with LO
diagonal propagators, the 2-point amplitude for the $a\text{-}\pi_{0}$ system
reads $\mathcal{P}_{ij}=\mbox{diag}\,(p^{2},p^{2}-m^{2}_{\pi})-\Sigma_{ij}$,
where $\Sigma_{ij}$ encodes NLO corrections including mixings. The 2-point
Green’s function $G_{ij}=(-i\mathcal{P})^{-1}_{ij}$ is hence

$G_{ij}=i\begin{pmatrix}\frac{1}{p^{2}}&\frac{\Sigma_{a\pi}}{p^{2}\left(p^{2}-m_{\pi}^{2}-\Sigma_{\pi\pi}\right)}\\\
\frac{\Sigma_{a\pi}}{p^{2}\left(p^{2}-m_{\pi}^{2}-\Sigma_{\pi\pi}\right)}&\frac{1}{p^{2}-m_{\pi}^{2}-\Sigma_{\pi\pi}}\end{pmatrix}\,.$
(10)

Plugging Eq. (9) and (10) into the LSZ formula for the scattering amplitude
and neglecting $\mathcal{O}(1/f_{a})^{2}$ terms, one finds (with $Z_{a}=1$,
$Z_{\pi}=1+\Sigma^{\prime}_{\pi\pi}(m^{2}_{\pi})$ and primes indicating
derivatives with respect to $p^{2}$)

$\displaystyle\mathcal{M}_{a\pi_{0}\to\pi_{+}\pi_{-}}=\left(1+\frac{3}{2}\Sigma^{\prime}_{\pi\pi}(m^{2}_{\pi})\right){\cal
G}_{a\pi_{0}\pi_{+}\pi_{-}}^{\rm LO}$
$\displaystyle-\frac{\Sigma_{a\pi}(m_{a}^{2}=0)}{m^{2}_{\pi}}{\cal
G}_{\pi_{0}\pi_{0}\pi_{+}\pi_{-}}^{\rm LO}+{\cal
G}_{a\pi_{0}\pi_{+}\pi_{-}}^{\rm NLO}\,,$ (11)

where the ${\cal G}$’s are evaluated at the physical masses of the external
particles. The one-loop amplitudes have been computed in dimensional
regularization. To carry out the renormalization procedure in the (modified)
$\overline{\text{MS}}$ scheme, we define the scale independent parameters
$\overline{\ell_{i}}$ as Gasser and Leutwyler (1984)

$\ell_{i}=\frac{\gamma_{i}}{32\pi^{2}}\left[\overline{\ell_{i}}+R+\ln\left(\frac{m_{\pi}^{2}}{\mu^{2}}\right)\right]\,,$
(12)

with $R=\frac{2}{d-4}-\log(4\pi)+\gamma_{E}-1$, in order to cancel the
divergent terms (in the limit $d=4$) with a suitable choice of the
$\gamma_{i}$. Eventually, only the terms proportional to $\ell_{1,2,7}$
contribute to the NLO amplitude, which is renormalized for $\gamma_{1}=1/3$,
$\gamma_{2}=2/3$ and $\gamma_{7}=0$. The latter coincide with the values
obtained in Ref. Gasser and Leutwyler (1984) for the standard chiral theory
without the axion.

The renormalized NLO amplitude for the $a\pi_{0}\rightarrow\pi_{+}\pi_{-}$
process (and its crossed channels) is given in Supplemetary Material. We have
also checked that the same analytical result is obtained via a direct NLO
diagonalization of the $a$ and $\pi^{0}$ propagators, without employing the
LSZ formalism with off-diagonal propagators. For consistency, we will only
consider the interference between the LO and NLO terms in the squared matrix
elements, $\sum|\mathcal{M}|^{2}\simeq\sum|\mathcal{M}|_{\rm LO}^{2}+\sum
2\mbox{Re}\,[\mathcal{M}_{\rm LO}\mathcal{M}^{*}_{\rm NLO}]$, since the NLO
squared correction is of the same order of the NNLO-LO interference, which we
neglect.

Breakdown of the chiral expansion at finite temperature. The crucial quantity
that is needed to extract the HDM bound is the axion decoupling temperature,
$T_{D}$, obtained via the freeze-out condition (following the same criterium
as in Hannestad  _et al._ (2005))

$\Gamma_{a}(T_{D})=H(T_{D})\,.$ (13)

Here, $H(T)=\sqrt{4\pi^{3}g_{\star}(T)/45}\,T^{2}/m_{\rm pl}$ is the Hubble
rate (assuming a radiation dominated Universe) in terms of the Planck mass
$m_{\rm pl}=1.22\times 10^{19}$ GeV and the effective number of relativistic
degrees of freedom, $g_{\star}(T)$, while $\Gamma_{a}$ is the axion
thermalization rate entering the Boltzmann equation

$\displaystyle\Gamma_{a}$ $\displaystyle=\frac{1}{n_{a}^{\rm
eq}}\int\frac{d^{3}\mathbf{p}_{1}}{(2\pi)^{3}2E_{1}}\frac{d^{3}\mathbf{p}_{2}}{(2\pi)^{3}2E_{2}}\frac{d^{3}\mathbf{p}_{3}}{(2\pi)^{3}2E_{3}}\frac{d^{3}\mathbf{p}_{4}}{(2\pi)^{3}2E_{4}}$
$\displaystyle\times\sum|\mathcal{M}|^{2}(2\pi)^{4}\delta^{4}\left(p_{1}+p_{2}-p_{3}-p_{4}\right)$
$\displaystyle\times f_{1}f_{2}(1+f_{3})(1+f_{4})\,,$ (14)

where $n_{a}^{\rm eq}=(\zeta_{3}/\pi^{2})T^{3}$ and $f_{i}=1/(e^{E_{i}/T}-1)$.

In the following, we will set the model-dependent axion couplings
$c^{0}_{u,\,d}=0$ (cf. Eq. (4)), to comply with the standard setup considered
in the literature Chang and Choi (1993); Hannestad  _et al._ (2005);
Melchiorri  _et al._ (2007); Hannestad  _et al._ (2008, 2010); Archidiacono
_et al._ (2013); Giusarma  _et al._ (2014); Di Valentino  _et al._ (2015,
2016); Archidiacono  _et al._ (2015); Giarè  _et al._ (2020) (see Ferreira
_et al._ (2020) for an exception). Moreover, we will neglect thermal
corrections to the scattering matrix element, since those are small for
$T\lesssim m_{\pi}$ Gasser and Leutwyler (1987a, b); Gerber and Leutwyler
(1989).

Figure 1: Numerical profile of the $h_{\rm LO}$ and $h_{\rm NLO}$ functions
entering the axion-pion thermalization rate in Eq. (Axion hot dark matter
bound, reliably).

By integrating numerically the phase space in Eq. (Axion hot dark matter
bound, reliably) and neglecting third order terms in the isospin breaking, we
find (see Supplementary Material for a useful intermediate analytical step, or
Hannestad and Madsen (1995) for a slightly different approach)

$\displaystyle\Gamma_{a}(T)$
$\displaystyle=\left(\frac{C_{a\pi}}{f_{a}f_{\pi}}\right)^{2}0.212\
T^{5}\Big{[}h_{\rm LO}(m_{\pi}/T)$
$\displaystyle-2.92\frac{T^{2}}{f_{\pi}^{2}}\ h_{\rm
NLO}(m_{\pi}/T)\Big{]}\,,$ (15)

where for the numerical evaluation we used the central values of the LECs
$\overline{\ell_{1}}=-0.36(59)$ Colangelo  _et al._ (2001),
$\overline{\ell_{2}}=4.31(11)$ Colangelo  _et al._ (2001),
$\overline{\ell_{3}}=3.53(26)$ Aoki  _et al._ (2020),
$\overline{\ell_{4}}=4.73(10)$ Aoki  _et al._ (2020) and $\ell_{7}=7(4)\times
10^{-3}$ Grilli di Cortona  _et al._ (2016), $m_{u}/m_{d}=0.50(2)$ Aoki  _et
al._ (2020), $f_{\pi}=92.1(8)$ MeV Zyla  _et al._ (2020) and $m_{\pi}=137$ MeV
(corresponding to the average neutral/charged pion mass). The $h$-functions
are normalized to $h_{\rm LO}(0)=h_{\rm NLO}(0)=1$ and plotted in Fig. 1.

We have checked that $h_{\rm LO}$ reproduces the result of Ref. Hannestad  _et
al._ (2005) within percent accuracy. It should be noted that Eq. (Axion hot
dark matter bound, reliably) is meaningful only for $m_{\pi}/T\gtrsim 1$,
since at higher temperatures above $T_{c}$ pions are deconfined and the axion
thermalization rate should be computed from the interactions with a quark-
gluon plasma. Nevertheless, we are interested in extrapolating the behaviour
of Eq. (Axion hot dark matter bound, reliably) from the low-temperature
regime, where the chiral approach is reliable.

In Fig. 2 we compare the LO and NLO rates contributing to
$\Gamma_{a}=\Gamma_{a}^{\rm LO}+\Gamma_{a}^{\rm NLO}$. In particular, the
$|\Gamma_{a}^{\rm NLO}/\Gamma_{a}^{\rm LO}|$ ratio does not depend on $f_{a}$.
Requiring as a loose criterium that the NLO correction is less than $50\%$ of
the LO one, yields $T_{\chi}\simeq 62$ MeV as the maximal temperature at which
the chiral description of the thermalization rate can be reliably extended.

Figure 2: Ratio between the NLO and the LO axion-pion thermalization rate.
$T_{\chi}\simeq 62$ MeV corresponds to a NLO correction of $50\%$.

Fig. 3 shows instead the extraction of the decoupling temperature (defined via
Eq. (13)) for two reference values of the axion mass (setting the strength of
the axion coupling via $f_{a}$), namely $m_{a}=1$ eV and 0.1 eV. Assuming a
standard analysis employing the LO axion thermalization rate Hannestad  _et
al._ (2005), the former benchmark (1 eV) corresponds to the most conservative
HDM bound Giarè  _et al._ (2020), while the latter (0.1 eV) saturates the most
stringent one Giarè  _et al._ (2020) and also represents the typical reach of
future CMB-S4 experiments Abazajian  _et al._ (2016). However, from Fig. 3 we
see that $T_{D}^{\rm LO}\simeq 59$ MeV for $m_{a}=1$ eV and $T_{D}^{\rm
LO}\sim 200$ MeV for $m_{a}=0.1$ eV. Hence, while in the former case the
decoupling temperature is at the boundary of validity of the chiral expansion,
set by $T_{\chi}\simeq 62$ MeV, in the latter is well above it. In particular,
the region where the chiral expansion fails, $T_{D}\gtrsim T_{\chi}$,
corresponds to $m_{a}\lesssim 1.2$ eV.

We hence conclude that in the mass range of interest, $m_{a}\in[0.1,1]$ eV,
the decoupling temperature and consequently the axion HDM bound cannot be
reliably extracted within the chiral approach. Note, also, that in the
presence of model-dependent axion couplings $c^{0}_{u,d}\gg 1$ (as in some
axion models Darmé  _et al._ (2020)), the same decoupling temperature as in
the $c^{0}_{u,d}=0$ case is obtained for larger $f_{a}$, thus shifting down
the mass window relevant for the axion HDM bound.

Figure 3: Axion-pion thermalization rate vs. Hubble rate for two reference
values of the axion mass, $m_{a}=1$ eV and 0.1 eV. The full $\Gamma_{a}$ has
been stopped at $T\simeq 85$ MeV, for which $|\Gamma^{\rm NLO}_{a}/\Gamma^{\rm
LO}_{a}|=90\%$.

Towards a reliable axion HDM bound. The failure of the chiral approach in the
calculation of the axion-pion thermalization rate can be traced back to the
fact that in a thermal bath with temperatures of the order of $T\simeq 100$
MeV the mean energy of pions is $\left\langle E_{\pi}\right\rangle\simeq 350$
MeV, so that $\pi$-$\pi$ scatterings happen at center of mass energies above
the validity of the 2-flavour chiral EFT. The latter can be related to the
scale of tree-level unitarity violation of $\pi$-$\pi$ scattering resulting in
$\sqrt{s}\lesssim\sqrt{8\pi}f_{\pi}\simeq 460$ MeV Weinberg (1966); Aydemir
_et al._ (2012). A possible strategy to extend the theoretical predictions at
higher energies is to compute the relevant $a\pi\to\pi\pi$ amplitudes using
lattice QCD simulations. To this end one may employ the standard techniques
used to compute weak non-leptonic matrix elements Blum  _et al._ (2015);
Abbott  _et al._ (2020) and $\pi$-$\pi$ scattering amplitudes as a function of
the energy at finite volume Luscher (1991); Rummukainen and Gottlieb (1995);
Kim  _et al._ (2005); Hansen and Sharpe (2012). Although this approach has
limitations with respect to the maximum attainable center of mass energy, we
believe that it can be used to compute the amplitudes up to values of
$\sqrt{s}\sim 600-900$ MeV or higher Briceno  _et al._ (2017).

We conclude by stressing the importance of obtaining a reliable determination
of the axion-pion thermalization rate, not only in view of the extraction of a
notable bound in axion physics, but also in order to set definite targets for
future CMB probes of the axion-pion coupling, which could represent a
‘discovery channel’ for the axion.

###### Acknowledgements.

Acknowledgments. We thank Enrico Nardi and Maurizio Giannotti for helpful
discussions. The work of LDL is supported by the Marie Skłodowska-Curie
Individual Fellowship grant AXIONRUSH (GA 840791) and the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy \- EXC 2121 Quantum
Universe - 390833306. The work of GP has received funding from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement N∘ 860881\.

## References

  * Peccei and Quinn (1977a) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn (1977b) R. D. Peccei and H. R. Quinn, Phys. Rev. D16, 1791 (1977b).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. 120B, 127 (1983).
  * Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Phys. Lett. 120B, 133 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. 120B, 137 (1983).
  * Davis (1986) R. L. Davis, Phys. Lett. B 180, 225 (1986).
  * Turner (1987) M. S. Turner, Phys. Rev. Lett. 59, 2489 (1987), [Erratum: Phys.Rev.Lett. 60, 1101 (1988)].
  * Kolb and Turner (1990) E. W. Kolb and M. S. Turner,  _The Early Universe_ , Vol. 69 (1990).
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A1 (2020a), arXiv:1807.06205 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020b), arXiv:1807.06209 [astro-ph.CO] .
  * Abazajian  _et al._ (2016) K. N. Abazajian _et al._ (CMB-S4), (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Salvio  _et al._ (2014) A. Salvio, A. Strumia, and W. Xue, JCAP 01, 011 (2014), arXiv:1310.6982 [hep-ph] .
  * Baumann  _et al._ (2016) D. Baumann, D. Green, and B. Wallisch, Phys. Rev. Lett. 117, 171301 (2016), arXiv:1604.08614 [astro-ph.CO] .
  * Ferreira and Notari (2018) R. Z. Ferreira and A. Notari, Phys. Rev. Lett. 120, 191301 (2018), arXiv:1801.06090 [hep-ph] .
  * Arias-Aragon  _et al._ (2020) F. Arias-Aragon, F. D’Eramo, R. Z. Ferreira, L. Merlo, and A. Notari, (2020), arXiv:2012.04736 [hep-ph] .
  * D’Eramo  _et al._ (2018) F. D’Eramo, R. Z. Ferreira, A. Notari, and J. L. Bernal, JCAP 11, 014 (2018), arXiv:1808.07430 [hep-ph] .
  * Masso  _et al._ (2002) E. Masso, F. Rota, and G. Zsembinszki, Phys. Rev. D 66, 023004 (2002), arXiv:hep-ph/0203221 .
  * Graf and Steffen (2011) P. Graf and F. D. Steffen, Phys. Rev. D 83, 075011 (2011), arXiv:1008.4528 [hep-ph] .
  * Berezhiani  _et al._ (1992) Z. Berezhiani, A. Sakharov, and M. Khlopov, Sov. J. Nucl. Phys. 55, 1063 (1992).
  * Chang and Choi (1993) S. Chang and K. Choi, Phys. Lett. B 316, 51 (1993), arXiv:hep-ph/9306216 .
  * Hannestad  _et al._ (2005) S. Hannestad, A. Mirizzi, and G. Raffelt, JCAP 07, 002 (2005), arXiv:hep-ph/0504059 .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Melchiorri  _et al._ (2007) A. Melchiorri, O. Mena, and A. Slosar, Phys. Rev. D 76, 041303 (2007), arXiv:0705.2695 [astro-ph] .
  * Hannestad  _et al._ (2008) S. Hannestad, A. Mirizzi, G. G. Raffelt, and Y. Y. Wong, JCAP 04, 019 (2008), arXiv:0803.1585 [astro-ph] .
  * Hannestad  _et al._ (2010) S. Hannestad, A. Mirizzi, G. G. Raffelt, and Y. Y. Wong, JCAP 08, 001 (2010), arXiv:1004.0695 [astro-ph.CO] .
  * Archidiacono  _et al._ (2013) M. Archidiacono, S. Hannestad, A. Mirizzi, G. Raffelt, and Y. Y. Wong, JCAP 10, 020 (2013), arXiv:1307.0615 [astro-ph.CO] .
  * Giusarma  _et al._ (2014) E. Giusarma, E. Di Valentino, M. Lattanzi, A. Melchiorri, and O. Mena, Phys. Rev. D 90, 043507 (2014), arXiv:1403.4852 [astro-ph.CO] .
  * Di Valentino  _et al._ (2015) E. Di Valentino, S. Gariazzo, E. Giusarma, and O. Mena, Phys. Rev. D 91, 123505 (2015), arXiv:1503.00911 [astro-ph.CO] .
  * Di Valentino  _et al._ (2016) E. Di Valentino, E. Giusarma, M. Lattanzi, O. Mena, A. Melchiorri, and J. Silk, Phys. Lett. B 752, 182 (2016), arXiv:1507.08665 [astro-ph.CO] .
  * Archidiacono  _et al._ (2015) M. Archidiacono, T. Basse, J. Hamann, S. Hannestad, G. Raffelt, and Y. Y. Wong, JCAP 05, 050 (2015), arXiv:1502.03325 [astro-ph.CO] .
  * Giarè  _et al._ (2020) W. Giarè, E. Di Valentino, A. Melchiorri, and O. Mena, (2020), arXiv:2011.14704 [astro-ph.CO] .
  * Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ (1980) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B166, 493 (1980).
  * Kaplan (1985) D. B. Kaplan, Nucl. Phys. B 260, 215 (1985).
  * Di Luzio  _et al._ (2017a) L. Di Luzio, F. Mescia, and E. Nardi, Phys. Rev. Lett. 118, 031801 (2017a), arXiv:1610.07593 [hep-ph] .
  * Di Luzio  _et al._ (2017b) L. Di Luzio, F. Mescia, and E. Nardi, Phys. Rev. D 96, 075003 (2017b), arXiv:1705.05370 [hep-ph] .
  * Raffelt (1990) G. G. Raffelt, Phys. Rept. 198, 1 (1990).
  * Chang  _et al._ (2018) J. H. Chang, R. Essig, and S. D. McDermott, JHEP 09, 051 (2018), arXiv:1803.00993 [hep-ph] .
  * Carenza  _et al._ (2019) P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo, and A. Mirizzi, JCAP 10, 016 (2019), [Erratum: JCAP 05, E01 (2020)], arXiv:1906.11844 [hep-ph] .
  * Bar  _et al._ (2020) N. Bar, K. Blum, and G. D’Amico, Phys. Rev. D 101, 123025 (2020), arXiv:1907.05020 [hep-ph] .
  * Bazavov  _et al._ (2012) A. Bazavov _et al._ , Phys. Rev. D 85, 054503 (2012), arXiv:1111.1710 [hep-lat] .
  * Spalinski (1988) M. Spalinski, Z. Phys. C 41, 87 (1988).
  * Grilli di Cortona  _et al._ (2016) G. Grilli di Cortona, E. Hardy, J. Pardo Vega, and G. Villadoro, JHEP 01, 034 (2016), arXiv:1511.02867 [hep-ph] .
  * Di Vecchia and Veneziano (1980) P. Di Vecchia and G. Veneziano, Nucl. Phys. B 171, 253 (1980).
  * Georgi  _et al._ (1986) H. Georgi, D. B. Kaplan, and L. Randall, Phys. Lett. B 169, 73 (1986).
  * Di Luzio  _et al._ (2020) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980), [Yad. Fiz.31,497(1980)].
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B104, 199 (1981).
  * Gasser and Leutwyler (1984) J. Gasser and H. Leutwyler, Annals Phys. 158, 142 (1984).
  * Scherer (2003) S. Scherer, Adv. Nucl. Phys. 27, 277 (2003), arXiv:hep-ph/0210398 .
  * Wess and Zumino (1971) J. Wess and B. Zumino, Phys. Lett. B 37, 95 (1971).
  * Witten (1983) E. Witten, Nucl. Phys. B 223, 422 (1983).
  * Lehmann  _et al._ (1955) H. Lehmann, K. Symanzik, and W. Zimmermann, Nuovo Cim. 1, 205 (1955).
  * Ferreira  _et al._ (2020) R. Z. Ferreira, A. Notari, and F. Rompineve, (2020), arXiv:2012.06566 [hep-ph] .
  * Gasser and Leutwyler (1987a) J. Gasser and H. Leutwyler, Phys. Lett. B 184, 83 (1987a).
  * Gasser and Leutwyler (1987b) J. Gasser and H. Leutwyler, Phys. Lett. B 188, 477 (1987b).
  * Gerber and Leutwyler (1989) P. Gerber and H. Leutwyler, Nucl. Phys. B 321, 387 (1989).
  * Hannestad and Madsen (1995) S. Hannestad and J. Madsen, Phys. Rev. D 52, 1764 (1995), arXiv:astro-ph/9506015 .
  * Colangelo  _et al._ (2001) G. Colangelo, J. Gasser, and H. Leutwyler, Nucl. Phys. B 603, 125 (2001), arXiv:hep-ph/0103088 .
  * Aoki  _et al._ (2020) S. Aoki _et al._ (Flavour Lattice Averaging Group), Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Darmé  _et al._ (2020) L. Darmé, L. Di Luzio, M. Giannotti, and E. Nardi, (2020), arXiv:2010.15846 [hep-ph] .
  * Weinberg (1966) S. Weinberg, Phys. Rev. Lett. 17, 616 (1966).
  * Aydemir  _et al._ (2012) U. Aydemir, M. M. Anber, and J. F. Donoghue, Phys. Rev. D 86, 014025 (2012), arXiv:1203.5153 [hep-ph] .
  * Blum  _et al._ (2015) T. Blum _et al._ , Phys. Rev. D 91, 074502 (2015), arXiv:1502.00263 [hep-lat] .
  * Abbott  _et al._ (2020) R. Abbott _et al._ (RBC, UKQCD), Phys. Rev. D 102, 054509 (2020), arXiv:2004.09440 [hep-lat] .
  * Luscher (1991) M. Luscher, Nucl. Phys. B 354, 531 (1991).
  * Rummukainen and Gottlieb (1995) K. Rummukainen and S. A. Gottlieb, Nucl. Phys. B 450, 397 (1995), arXiv:hep-lat/9503028 .
  * Kim  _et al._ (2005) C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl. Phys. B 727, 218 (2005), arXiv:hep-lat/0507006 .
  * Hansen and Sharpe (2012) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 86, 016007 (2012), arXiv:1204.0826 [hep-lat] .
  * Briceno  _et al._ (2017) R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. Lett. 118, 022002 (2017), arXiv:1607.05900 [hep-ph] .
  * Alloul  _et al._ (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. 185, 2250 (2014), arXiv:1310.1921 [hep-ph] .
  * Christensen and Duhr (2009) N. D. Christensen and C. Duhr, Comput. Phys. Commun. 180, 1614 (2009), arXiv:0806.4194 [hep-ph] .
  * Hahn (2001) T. Hahn, Comput. Phys. Commun. 140, 418 (2001), arXiv:hep-ph/0012260 .
  * Shtabovenko  _et al._ (2020) V. Shtabovenko, R. Mertig, and F. Orellana, Comput. Phys. Commun. 256, 107478 (2020), arXiv:2001.04407 [hep-ph] .
  * Shtabovenko  _et al._ (2016) V. Shtabovenko, R. Mertig, and F. Orellana, Comput. Phys. Commun. 207, 432 (2016), arXiv:1601.01167 [hep-ph] .
  * Mertig  _et al._ (1991) R. Mertig, M. Bohm, and A. Denner, Comput. Phys. Commun. 64, 345 (1991).
  * Patel (2015) H. H. Patel, Comput. Phys. Commun. 197, 276 (2015), arXiv:1503.01469 [hep-ph] .

##  I Supplementary Material

The calculation of the amplitudes was carried out using the computational
tools FeynRules Alloul  _et al._ (2014); Christensen and Duhr (2009), FeynArts
Hahn (2001), FeynCalc Shtabovenko  _et al._ (2020, 2016); Mertig  _et al._
(1991) and Package-X Patel (2015). The full analytical expression of the
renormalized NLO amplitude for the $a\pi_{0}\rightarrow\pi_{+}\pi_{-}$ process
reads

$\displaystyle\mathcal{M}^{\rm
NLO}_{a\pi_{0}\rightarrow\pi_{+}\pi_{-}}=\frac{C_{a\pi}}{192\pi^{2}f_{\pi}^{3}f_{a}}\Bigg{\\{}15m_{\pi}^{2}(u+t)-11u^{2}-8ut-11t^{2}-6\overline{\ell_{1}}\left(m_{\pi}^{2}-s\right)\left(2m_{\pi}^{2}-s\right)$
$\displaystyle-6\overline{\ell_{2}}\left(-3m_{\pi}^{2}(u+t)+4m_{\pi}^{4}+u^{2}+t^{2}\right)+9m_{\pi}^{4}\overline{\ell_{3}}+18\overline{\ell_{4}}m_{\pi}^{2}(m_{\pi}^{2}-s)+576\pi^{2}\ell_{7}m_{\pi}^{4}\left(\frac{m_{d}-m_{u}}{m_{d}+m_{u}}\right)^{2}$
$\displaystyle+3\left[3\sqrt{1-\frac{4m_{\pi}^{2}}{s}}s\left(m_{\pi}^{2}-s\right)\ln{\left(\frac{\sqrt{s\left(s-4m_{\pi}^{2}\right)}+2m_{\pi}^{2}-s}{2m_{\pi}^{2}}\right)}\right.$
$\displaystyle+\sqrt{1-\frac{4m_{\pi}^{2}}{t}}\left(m_{\pi}^{2}(t-4u)+3m_{\pi}^{4}+t(u-t)\right)\ln{\left(\frac{\sqrt{t\left(t-4m_{\pi}^{2}\right)}+2m_{\pi}^{2}-t}{2m_{\pi}^{2}}\right)}$
$\displaystyle\left.+\sqrt{1-\frac{4m_{\pi}^{2}}{u}}\left(m_{\pi}^{2}(u-4t)+3m_{\pi}^{4}+u(t-u)\right)\ln{\left(\frac{\sqrt{u\left(u-4m_{\pi}^{2}\right)}+2m_{\pi}^{2}-u}{2m_{\pi}^{2}}\right)}\right]\Bigg{\\}}$
$\displaystyle+\frac{4\ell_{7}m_{\pi}^{2}m_{d}\left(s-2m_{\pi}^{2}\right)m_{u}\left(m_{d}-m_{u}\right)}{f_{\pi}^{3}f_{a}\left(m_{d}+m_{u}\right){}^{3}}\,,$
(16)

where the terms proportional to $\overline{\ell_{3}}$, $\overline{\ell_{4}}$
and $\ell_{7}$ in the second row arise from the LO amplitude, via the NLO
corrections to $m_{\pi}$ and $f_{\pi}$ (see Ref. Gasser and Leutwyler (1984)).
The latter include the charged-neutral pion mass difference arising at second
order in the isospin breaking parameter $m_{d}-m_{u}$, which has been
neglected in the numerical integration of the rate. The amplitudes for the
crossed channels $a\pi_{-}\rightarrow\pi_{0}\pi_{-}$ and
$a\pi_{+}\rightarrow\pi_{+}\pi_{0}$ are obtained by cross symmetry through the
replacements $s\leftrightarrow t$ and $s\leftrightarrow u$, respectively.

Next, we describe our procedure to analytically reduce the 12-dimensional
phase space integral of Eq. (Axion hot dark matter bound, reliably) down to a
5-dimensional one. We first integrate the fourth-particle phase space in Eq.
(Axion hot dark matter bound, reliably) using the relation
$d^{3}\mathbf{p}_{4}/(2E_{4})=d^{4}p_{4}\delta\left(p_{4}^{2}-m_{4}^{2}\right)\theta(p_{4}^{0})$.
Therefore, defining the angles $\alpha$ and $\theta$ via
$\cos\alpha=\frac{\mathbf{p}_{1}\cdot\mathbf{p}_{2}}{|\mathbf{p}_{1}||\mathbf{p}_{2}|}$
and
$\cos\theta=\frac{\mathbf{p}_{1}\cdot\mathbf{p}_{3}}{|\mathbf{p}_{1}||\mathbf{p}_{3}|}$,
the thermalization rate becomes

$\displaystyle\Gamma_{a}$ $\displaystyle=\frac{1}{n_{a}^{\rm
eq}}\int\frac{dp_{1}|\mathbf{p}_{1}|^{2}}{2E_{1}}\frac{dp_{2}|\mathbf{p}_{2}|^{2}}{2E_{2}}\frac{dp_{3}|\mathbf{p}_{3}|^{2}}{2E_{3}}\int_{-1}^{1}d\cos\alpha\int_{-1}^{1}d\cos\theta\int_{0}^{2\pi}d\beta\sum|\mathcal{M}|^{2}\frac{4\pi}{(2\pi)^{7}}$
$\displaystyle\times\frac{\delta\left(E_{1}-\xi(E_{2},E_{3},\alpha,\beta,\theta)\right)}{2|E_{2}-E_{3}-|\mathbf{p}_{2}|\cos\alpha+|\mathbf{p}_{3}|\cos\theta|}f_{1}f_{2}(1+f_{3})(1+f_{4})\,,$
(17)

with $\beta$ the angle between the scattering planes defined by
$(\mathbf{p}_{1},\mathbf{p}_{2})$ and $(\mathbf{p}_{3},\mathbf{p}_{4})$, and
the function $\xi$ given by

$\xi(E_{2},E_{3},\alpha,\beta,\theta)=\frac{2\left(E_{2}E_{3}-|\mathbf{p}_{2}||\mathbf{p}_{3}|\left(\sin\alpha\sin\theta\cos\beta+\cos\alpha\cos\theta\right)\right)-m_{\pi}^{2}}{2(E_{2}-E_{3}-|\mathbf{p}_{2}|\cos\alpha+|\mathbf{p}_{3}|\cos\theta)}\,.$
(18)

Eq. (I) is then integrated numerically, leading to Eq. (Axion hot dark matter
bound, reliably).