Breakdown of chiral perturbation theory for the axion hot dark matter bound
AAxion hot dark matter bound, reliably
Luca Di Luzio, ∗ Guido Martinelli, † and Gioacchino Piazza ‡ DESY, Notkestraße 85, D-22607 Hamburg, Germany Physics Department and INFN Sezione di Roma La Sapienza,Piazzale Aldo Moro 5, 00185 Roma, Italy IJCLab, Pˆole Th´eorie (Bˆat. 210), CNRS/IN2P3 et Universit´e Paris-Saclay, 91405 Orsay, France
We show that the commonly adopted hot dark matter (HDM) bound on the axion mass m a (cid:46) Introduction.
The axion originally emerged as alow-energy remnant of the Peccei Quinn solution tothe strong CP problem [1–4], but it also unavoid-ably contributes to the energy density of the Uni-verse. There are two qualitatively different pop-ulations of relic axions, a non-thermal one com-prising cold dark matter (DM) [5–8], and a ther-mal axion population [9] which, while still relativis-tic, would behave as extra dark radiation. Suchhot dark matter (HDM) component contributes tothe effective number of extra relativistic degrees offreedom [10] ∆ N eff (cid:39) / / [4 g S ( T D )]) / , with g S ( T D ) the number of entropy degrees of freedomat the axion decoupling temperature, T D . Thevalue of ∆ N eff is constrained by cosmic microwavebackground (CMB) experiments, such as the Plancksatellite [11, 12], while planned CMB Stage 4 (CMB-S4) experiments [13] will provide an observable win-dow on the axion interactions.There are several processes that can keep the ax-ion in thermal equilibrium with the Standard Model(SM) thermal bath. From the standpoint of the ax-ion solution to the strong CP problem, an unavoid-able process arises from the model-independent cou-pling to gluons, α s π af a G ˜ G . For T D (cid:38) T D (cid:46) f a or,equivalently, by m a (cid:39) . × (10 GeV /f a ) eV, andit encompasses the range m a ∈ [0 . , .
1] eV (withheavier axions leading to lower decoupling temper-atures). Although the transition region cannot beprecisely determined due to the complications of thequark-hadron phase transition, for heavier axionsapproaching the eV scale the main thermalization Other thermalization channels arise from model-dependentaxion couplings to photons [9], SM quarks [14–17] and lep-tons [18]. channel is aπ ↔ ππ [22, 23], with T D (cid:46)
200 MeV.In this regime, scatterings off nucleons are subdomi-nant because of the exponential suppression in theirnumber density.The highest attainable axion mass from cosmo-logical constraints on extra relativistic degrees offreedom, also known as HDM bound, translate into m a (cid:46) m a (cid:46) It is the purpose of this Letter to revisit the ax-ion HDM bound in the context of the next-to-LO(NLO) axion-pion chiral EFT. This is motivated bythe simple observation that the mean energy of pi-ons (axions) in a heat bath of T (cid:39)
100 MeV is (cid:104) E (cid:105) ≡ ρ/n (cid:39)
350 MeV (270 MeV), thus questioningthe validity of the chiral expansion for the scatter-ing process aπ ↔ ππ . The latter is expected to failfor √ s ∼ (cid:104) E π (cid:105) + (cid:104) E a (cid:105) (cid:38)
500 MeV, correspondingto temperatures well below that of QCD deconfine-ment, T c = 154 ± Tree-level axion couplings to electrons are absent in KSVZmodels [34, 35], thus relaxing the constraints from RedGiants and White Dwarfs. The axion coupling to photons,constrained by Horizontal Branch stars evolution, can beaccidentally suppressed in certain KSVZ-like models [36–38]. Finally, the SN1987A bound on the axion couplingto nucleons can be considered less robust both from theastrophysical and experimental point of view [39–42]. a r X i v : . [ h e p - ph ] J a n of the aπ ↔ ππ thermalization rate (that can be castas an expansion in T /f π , with f π (cid:39)
92 MeV) andshow that the NLO correction saturates half of theLO contribution for T χ (cid:39)
62 MeV. The latter can beconsidered as the maximal temperature above whichthe chiral description breaks down for the processunder consideration. On the other hand, the regionfrom T χ up to T c , where chiral perturbation theorycannot be applied, turns out to be crucial for theextraction of the HDM bound and for assessing thesensitivity of future CMB experiments.We conclude with a proposal for extracting theaxion-pion thermalization rate via a direct LatticeQCD calculation, in analogy to the well-studied caseof π - π scattering. Axion-pion scattering at LO.
The constructionof the LO axion-pion Lagrangian was discussed longago in Refs. [46, 47]. We recall here its basic in-gredients (see also [22, 48]), in view of the exten-sion at NLO. Defining the pion Goldstone matrix U = e iπ A σ A /f π , with f π (cid:39)
92 MeV, π A and σ A ( A = 1 , ,
3) denoting respectively the real pionfields and the Pauli matrices, the LO axion-pion in-teractions stem from L LO a - π = f π (cid:2) U χ † a + χ a U † (cid:3) + ∂ µ a f a Tr (cid:2) c q σ A (cid:3) J Aµ , (1)where χ a = 2 B M a , in terms of the quark conden-sate B and the ‘axion-dressed’ quark mass matrix M a = e i a fa Q a M q e i a fa Q a , with M q = diag ( m u , m d )and Tr Q a = 1. The latter condition ensuresthat the axion field is transferred from the opera-tor α s π af a G ˜ G to the phase of the quark mass ma-trix, via the quark axial field redefinition q → exp( iγ a f a Q a ) q . In the following, we set Q a = M − q / Tr M − q , so that terms linear in a (including a - π mass mixing) drop out from the first term inEq. (1). Hence, in this basis, the only linear axioninteraction is the derivative one with the conservedSU(2) A pion current. The latter reads at LO J Aµ | LO = i f π Tr (cid:2) σ A (cid:0) U ∂ µ U † − U † ∂ µ U (cid:1)(cid:3) , (2)while the derivative axion coupling in Eq. (1) isTr (cid:2) c q σ A (cid:3) = ( m u − m d m u + m d + c u − c d ) δ A , where thefirst term arises from the axial quark rotationthat removed the aG ˜ G operator and the secondone originates from the model-dependent coefficient c q = diag( c u , c d ), defined via the Lagrangian term ∂ µ a f a qc q γ µ γ q . For instance, c u,d = 0 in the KSVZmodel [34, 35], while c u = cos β and c d = sin β in the DFSZ model [49, 50], with tan β the ratio be-tween the vacuum expectation values of two Higgsdoublets. Expanding the pion matrix in Eq. (1) oneobtains L LO a - π ⊃ (cid:15) ∂ µ a∂ µ π + C aπ f a f π ∂ µ a [ ∂πππ ] µ , (3) with the definitions [ ∂πππ ] µ = 2 ∂ µ π π + π − − π ∂ µ π + π − − π π + ∂ µ π − , (cid:15) = − f π C aπ f a and C aπ = 13 (cid:18) m d − m u m u + m d + c d − c u (cid:19) . (4)At the LO in (cid:15) the diagonalization of the a - π termis obtained by shifting a → a − (cid:15)π and π → π + O ( (cid:15) ) a , where we used the fact that m a /m π = O ( (cid:15) ).Hence, as long as we are interested in effects that arelinear in a and neglect O ( (cid:15) ) corrections, the axion-pion interactions in Eq. (3) are already in the basiswith canonical propagators.For temperatures below the QCD phase transi-tion, the main processes relevant for the axion ther-malization rate are a ( p ) π ( p ) → π + ( p ) π − ( p ),whose amplitude at LO reads M LO aπ → π + π − = C aπ f π f a (cid:2) m π − s (cid:3) , (5)with s = ( p + p ) , together with the crossedchannels aπ − → π π − and aπ + → π + π . Theamplitudes of the latter are obtained by replacing s → t = ( p − p ) and s → u = ( p − p ) , re-spectively. Taking equal masses for the neutral andcharged pions, one finds the squared matrix element(summed over the three channels above) [23] (cid:88) |M| = (cid:18) C aπ f a f π (cid:19) (cid:2) s + t + u − m π (cid:3) . (6) Axion-pion scattering at NLO.
To compute theaxion thermalization process beyond LO we needto consider the one-loop amplitudes from the LOLagrangian in Eq. (1) as well as the tree-level am-plitudes stemming from the NLO axion-pion La-grangian, both contributing to O ( p ) in the chi-ral expansion. The NLO interactions include thederivative coupling of the axion to the NLO axialcurrent, which has been computed here for the firsttime.We stick to the expression of the NLO chiral La-grangian given in Ref. [51] (see for example Ap-pendix D in [52] for trace notation), which, consid-ering only two flavours, depends on 10 low-energyconstants (LECs) (cid:96) , (cid:96) , . . . , (cid:96) , h , h , h . The ax-ion field has been included in the phase of the quarkmass matrix, as described after Eq. (1). Note thatsince we are interested in 2 → O ( p ) Wess-Zumino-Witten term[53, 54] since it contains operators with an odd num-ber of bosons.To compute the axial current J Aµ at NLO, we pro-mote the ordinary derivative to a covariant one, de-fined as D µ U = ∂ µ U − ir µ U + iU l µ , with r µ = r Aµ σ A / l µ = l Aµ σ A / l Aµ and r Aµ , respectively. Taking the R − L combina-tion and switching off the external fields, the NLOaxial current reads J Aµ | NLO = i (cid:96) Tr (cid:2) σ A (cid:8) ∂ µ U † , U (cid:9)(cid:3) Tr (cid:2) ∂ ν U ∂ ν U † (cid:3) + i (cid:96) Tr (cid:2) σ A (cid:8) ∂ ν U † , U (cid:9)(cid:3) Tr (cid:2) ∂ µ U ∂ ν U † + ∂ ν U ∂ µ U † (cid:3) − i (cid:96) Tr (cid:2) σ A (cid:8) ∂ µ U, χ † a (cid:9) − σ A (cid:8) U, ∂ µ χ † a (cid:9) + σ A (cid:8) ∂ µ χ a , U † (cid:9) − σ A (cid:8) χ a , ∂ µ U † (cid:9) (cid:3) , (7)where curly brackets indicate anti-commutators.New a - π mixings arise at NLO, both at treelevel from the NLO Lagrangian and at one loopfrom L LO a - π . These mixings are explicitly taken intoaccount in the Lehmann-Symanzik-Zimmermann(LSZ) reduction formula [55] (focussing e.g. on the aπ → π + π − channel) M aπ → π + π − = 1 (cid:112) Z a Z π Π i =1 lim p i → m i (cid:0) p i − m i (cid:1) × G aπ π + π − ( p , p , p , p ) , (8)where the index i runs over the external particles, Z a ( Z π ) is the wave-function renormalization of theaxion (pion) field and the full 4-point Green’s func-tion is given by G aπ π + π − = (cid:88) i,j = a,π G ijπ + π − (9) × G π + π + ( m π ) G π − π − ( m π ) G ai ( m a = 0) G π j ( m π ) . The first term is the amputated 4-point function,multiplied by the 2-point functions of the externallegs with the axion mass to zero. Working with LOdiagonal propagators, the 2-point amplitude for the a - π system reads P ij = diag ( p , p − m π ) − Σ ij ,where Σ ij encodes NLO corrections including mix-ings. The 2-point Green’s function G ij = ( − i P ) − ij is hence G ij = i (cid:32) p Σ aπ p ( p − m π − Σ ππ )Σ aπ p ( p − m π − Σ ππ ) 1 p − m π − Σ ππ (cid:33) . (10)Plugging Eq. (9) and (10) into the LSZ formula forthe scattering amplitude and neglecting O (1 /f a ) terms, one finds (with Z a = 1, Z π = 1 + Σ (cid:48) ππ ( m π )and primes indicating derivatives with respect to p ) M aπ → π + π − = (cid:18) (cid:48) ππ ( m π ) (cid:19) G LO aπ π + π − − Σ aπ ( m a = 0) m π G LO π π π + π − + G NLO aπ π + π − , (11)where the G ’s are evaluated at the physical massesof the external particles. The one-loop amplitudeshave been computed in dimensional regularization.To carry out the renormalization procedure in the (modified) MS scheme, we define the scale indepen-dent parameters (cid:96) i as [51] (cid:96) i = γ i π (cid:20) (cid:96) i + R + ln (cid:18) m π µ (cid:19)(cid:21) , (12)with R = d − − log(4 π ) + γ E −
1, in order to can-cel the divergent terms (in the limit d = 4) with asuitable choice of the γ i . Eventually, only the termsproportional to (cid:96) , , contribute to the NLO ampli-tude, which is renormalized for γ = 1 / γ = 2 / γ = 0. The latter coincide with the valuesobtained in Ref. [51] for the standard chiral theorywithout the axion.The renormalized NLO amplitude for the aπ → π + π − process (and its crossed channels) is givenin Supplemetary Material. We have also checkedthat the same analytical result is obtained via a di-rect NLO diagonalization of the a and π propaga-tors, without employing the LSZ formalism with off-diagonal propagators. For consistency, we will onlyconsider the interference between the LO and NLOterms in the squared matrix elements, (cid:80) |M| (cid:39) (cid:80) |M| + (cid:80) M LO M ∗ NLO ], since the NLOsquared correction is of the same order of the NNLO-LO interference, which we neglect.
Breakdown of the chiral expansion at finitetemperature.
The crucial quantity that is neededto extract the HDM bound is the axion decouplingtemperature, T D , obtained via the freeze-out condi-tion (following the same criterium as in [23])Γ a ( T D ) = H ( T D ) . (13)Here, H ( T ) = (cid:112) π g (cid:63) ( T ) / T /m pl is the Hubblerate (assuming a radiation dominated Universe) interms of the Planck mass m pl = 1 . × GeVand the effective number of relativistic degrees offreedom, g (cid:63) ( T ), while Γ a is the axion thermalizationrate entering the Boltzmann equationΓ a = 1 n eq a (cid:90) d p (2 π ) E d p (2 π ) E d p (2 π ) E d p (2 π ) E × (cid:88) |M| (2 π ) δ ( p + p − p − p ) × f f (1 + f )(1 + f ) , (14)where n eq a = ( ζ /π ) T and f i = 1 / ( e E i /T − c u, d = 0 (cf. Eq. (4)), to complywith the standard setup considered in the literature[22, 23, 25–33] (see [56] for an exception). Moreover,we will neglect thermal corrections to the scatteringmatrix element, since those are small for T (cid:46) m π [57–59]. By integrating numerically the phase spacein Eq. (14) and neglecting third order terms in theisospin breaking, we find (see Supplementary Mate-rial for a useful intermediate analytical step, or [60] h LO h NLO - - m π / T h - f un c t i on s FIG. 1. Numerical profile of the h LO and h NLO functionsentering the axion-pion thermalization rate in Eq. (15). for a slightly different approach)Γ a ( T ) = (cid:18) C aπ f a f π (cid:19) . T (cid:104) h LO ( m π /T ) − . T f π h NLO ( m π /T ) (cid:105) , (15)where for the numerical evaluation we used thecentral values of the LECs (cid:96) = − . (cid:96) = 4 . (cid:96) = 3 . (cid:96) = 4 . (cid:96) = 7(4) × − [45], m u /m d = 0 . f π = 92 . m π = 137 MeV(corresponding to the average neutral/charged pionmass). The h -functions are normalized to h LO (0) = h NLO (0) = 1 and plotted in Fig. 1.We have checked that h LO reproduces the resultof Ref. [23] within percent accuracy. It should benoted that Eq. (15) is meaningful only for m π /T (cid:38)
1, since at higher temperatures above T c pions aredeconfined and the axion thermalization rate shouldbe computed from the interactions with a quark-gluon plasma. Nevertheless, we are interested inextrapolating the behaviour of Eq. (15) from thelow-temperature regime, where the chiral approachis reliable.In Fig. 2 we compare the LO and NLO rates con-tributing to Γ a = Γ LO a + Γ NLO a . In particular, the | Γ NLO a / Γ LO a | ratio does not depend on f a . Requir-ing as a loose criterium that the NLO correction isless than 50% of the LO one, yields T χ (cid:39)
62 MeVas the maximal temperature at which the chiral de-scription of the thermalization rate can be reliablyextended.Fig. 3 shows instead the extraction of the decou-pling temperature (defined via Eq. (13)) for two ref-erence values of the axion mass (setting the strengthof the axion coupling via f a ), namely m a = 1 eV and0.1 eV. Assuming a standard analysis employing theLO axion thermalization rate [23], the former bench-mark (1 eV) corresponds to the most conservativeHDM bound [33], while the latter (0.1 eV) saturatesthe most stringent one [33] and also represents thetypical reach of future CMB-S4 experiments [13].However, from Fig. 3 we see that T LO D (cid:39)
59 MeV f π T χ T c % % → chiral expansionbreaks down
20 40 60 80 100 120 140 1600.00.51.01.52.02.53.03.5 T [ MeV ] | Γ N L O / Γ L O | FIG. 2. Ratio between the NLO and the LO axion-pionthermalization rate. T χ (cid:39)
62 MeV corresponds to aNLO correction of 50%. f π T χ T c H ( T ) m = m = Γ Γ LO | Γ NLO |
20 40 60 80 100 120 140 1600.0010.0100.100110100 T [ MeV ] R a t e s i nun i t s o f [ T / m P l ] FIG. 3. Axion-pion thermalization rate vs. Hubble ratefor two reference values of the axion mass, m a = 1 eVand 0.1 eV. The full Γ a has been stopped at T (cid:39) | Γ NLO a / Γ LO a | = 90%. for m a = 1 eV and T LO D ∼
200 MeV for m a = 0 . T χ (cid:39)
62 MeV, in the lat-ter is well above it. In particular, the region wherethe chiral expansion fails, T D (cid:38) T χ , corresponds to m a (cid:46) . m a ∈ [0 . ,
1] eV, the decoupling temperatureand consequently the axion HDM bound cannot bereliably extracted within the chiral approach. Note,also, that in the presence of model-dependent axioncouplings c u,d (cid:29) c u,d = 0case is obtained for larger f a , thus shifting down themass window relevant for the axion HDM bound. Towards a reliable axion HDM bound.
Thefailure of the chiral approach in the calculation ofthe axion-pion thermalization rate can be tracedback to the fact that in a thermal bath with temper-atures of the order of T (cid:39)
100 MeV the mean energyof pions is (cid:104) E π (cid:105) (cid:39)
350 MeV, so that π - π scatteringshappen at center of mass energies above the valid-ity of the 2-flavour chiral EFT. The latter can berelated to the scale of tree-level unitarity violationof π - π scattering resulting in √ s (cid:46) √ πf π (cid:39) aπ → ππ amplitudes using lat-tice QCD simulations. To this end one may em-ploy the standard techniques used to compute weaknon-leptonic matrix elements [66, 67] and π - π scat-tering amplitudes as a function of the energy at fi-nite volume [68–71]. Although this approach haslimitations with respect to the maximum attain-able center of mass energy, we believe that it canbe used to compute the amplitudes up to values of √ s ∼ −
900 MeV or higher [72].We conclude by stressing the importance of ob-taining a reliable determination of the axion-pionthermalization rate, not only in view of the extrac-tion of a notable bound in axion physics, but also inorder to set definite targets for future CMB probesof the axion-pion coupling, which could represent a‘discovery channel’ for the axion.
Acknowledgments.
We thank Enrico Nardi andMaurizio Giannotti for helpful discussions. Thework of LDL is supported by the Marie Sk(cid:32)lodowska-Curie Individual Fellowship grant AXIONRUSH(GA 840791) and the Deutsche Forschungsgemein-schaft under Germany’s Excellence Strategy - EXC2121 Quantum Universe - 390833306. The work ofGP has received funding from the European Union’sHorizon 2020 research and innovation programmeunder the Marie Sk(cid:32)lodowska-Curie grant agreementN ◦ ∗ [email protected] † [email protected] ‡ [email protected][1] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. ,1440 (1977).[2] R. D. Peccei and H. R. Quinn, Phys. Rev. D16 ,1791 (1977).[3] F. Wilczek, Phys. Rev. Lett. , 279 (1978).[4] S. Weinberg, Phys. Rev. Lett. , 223 (1978).[5] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. , 127 (1983).[6] L. F. Abbott and P. Sikivie, Phys. Lett. , 133(1983).[7] M. Dine and W. Fischler, Phys. Lett. , 137(1983).[8] R. L. Davis, Phys. Lett. B , 225 (1986).[9] M. S. Turner, Phys. Rev. Lett. , 2489 (1987), [Er-ratum: Phys.Rev.Lett. 60, 1101 (1988)].[10] E. W. Kolb and M. S. Turner, The Early Universe ,Vol. 69 (1990).[11] N. Aghanim et al. (Planck), Astron. Astrophys. , A1 (2020), arXiv:1807.06205 [astro-ph.CO].[12] N. Aghanim et al. (Planck), Astron. Astrophys. , A6 (2020), arXiv:1807.06209 [astro-ph.CO].[13] K. N. Abazajian et al. (CMB-S4), (2016),arXiv:1610.02743 [astro-ph.CO]. [14] A. Salvio, A. Strumia, and W. Xue, JCAP , 011(2014), arXiv:1310.6982 [hep-ph].[15] D. Baumann, D. Green, and B. Wallisch, Phys.Rev. Lett. , 171301 (2016), arXiv:1604.08614[astro-ph.CO].[16] R. Z. Ferreira and A. Notari, Phys. Rev. Lett. ,191301 (2018), arXiv:1801.06090 [hep-ph].[17] F. Arias-Aragon, F. D’Eramo, R. Z. Ferreira,L. Merlo, and A. Notari, (2020), arXiv:2012.04736[hep-ph].[18] F. D’Eramo, R. Z. Ferreira, A. Notari, and J. L.Bernal, JCAP , 014 (2018), arXiv:1808.07430[hep-ph].[19] E. Masso, F. Rota, and G. Zsembinszki, Phys. Rev.D , 023004 (2002), arXiv:hep-ph/0203221.[20] P. Graf and F. D. Steffen, Phys. Rev. D , 075011(2011), arXiv:1008.4528 [hep-ph].[21] Z. Berezhiani, A. Sakharov, and M. Khlopov, Sov.J. Nucl. Phys. , 1063 (1992).[22] S. Chang and K. Choi, Phys. Lett. B , 51 (1993),arXiv:hep-ph/9306216.[23] S. Hannestad, A. Mirizzi, and G. Raffelt, JCAP , 002 (2005), arXiv:hep-ph/0504059.[24] P. Zyla et al. (Particle Data Group), PTEP ,083C01 (2020).[25] A. Melchiorri, O. Mena, and A. Slosar, Phys. Rev.D , 041303 (2007), arXiv:0705.2695 [astro-ph].[26] S. Hannestad, A. Mirizzi, G. G. Raffelt, andY. Y. Wong, JCAP , 019 (2008), arXiv:0803.1585[astro-ph].[27] S. Hannestad, A. Mirizzi, G. G. Raffelt, andY. Y. Wong, JCAP , 001 (2010), arXiv:1004.0695[astro-ph.CO].[28] M. Archidiacono, S. Hannestad, A. Mirizzi, G. Raf-felt, and Y. Y. Wong, JCAP , 020 (2013),arXiv:1307.0615 [astro-ph.CO].[29] E. Giusarma, E. Di Valentino, M. Lattanzi, A. Mel-chiorri, and O. Mena, Phys. Rev. D , 043507(2014), arXiv:1403.4852 [astro-ph.CO].[30] E. Di Valentino, S. Gariazzo, E. Giusarma, andO. Mena, Phys. Rev. D , 123505 (2015),arXiv:1503.00911 [astro-ph.CO].[31] E. Di Valentino, E. Giusarma, M. Lattanzi,O. Mena, A. Melchiorri, and J. Silk, Phys. Lett. B , 182 (2016), arXiv:1507.08665 [astro-ph.CO].[32] M. Archidiacono, T. Basse, J. Hamann, S. Hannes-tad, G. Raffelt, and Y. Y. Wong, JCAP , 050(2015), arXiv:1502.03325 [astro-ph.CO].[33] W. Giar`e, E. Di Valentino, A. Melchiorri, andO. Mena, (2020), arXiv:2011.14704 [astro-ph.CO].[34] J. E. Kim, Phys. Rev. Lett. , 103 (1979).[35] M. A. Shifman, A. I. Vainshtein, and V. I. Za-kharov, Nucl. Phys. B166 , 493 (1980).[36] D. B. Kaplan, Nucl. Phys. B , 215 (1985).[37] L. Di Luzio, F. Mescia, and E. Nardi, Phys. Rev.Lett. , 031801 (2017), arXiv:1610.07593 [hep-ph].[38] L. Di Luzio, F. Mescia, and E. Nardi, Phys. Rev.D , 075003 (2017), arXiv:1705.05370 [hep-ph].[39] G. G. Raffelt, Phys. Rept. , 1 (1990).[40] J. H. Chang, R. Essig, and S. D. McDermott, JHEP , 051 (2018), arXiv:1803.00993 [hep-ph].[41] P. Carenza, T. Fischer, M. Giannotti, G. Guo,G. Mart´ınez-Pinedo, and A. Mirizzi, JCAP ,016 (2019), [Erratum: JCAP 05, E01 (2020)],arXiv:1906.11844 [hep-ph].[42] N. Bar, K. Blum, and G. D’Amico, Phys. Rev. D , 123025 (2020), arXiv:1907.05020 [hep-ph]. [43] A. Bazavov et al. , Phys. Rev. D , 054503 (2012),arXiv:1111.1710 [hep-lat].[44] M. Spalinski, Z. Phys. C , 87 (1988).[45] G. Grilli di Cortona, E. Hardy, J. Pardo Vega,and G. Villadoro, JHEP , 034 (2016),arXiv:1511.02867 [hep-ph].[46] P. Di Vecchia and G. Veneziano, Nucl. Phys. B ,253 (1980).[47] H. Georgi, D. B. Kaplan, and L. Randall, Phys.Lett. B , 73 (1986).[48] L. Di Luzio, M. Giannotti, E. Nardi, andL. Visinelli, Phys. Rept. , 1 (2020),arXiv:2003.01100 [hep-ph].[49] A. R. Zhitnitsky, Sov. J. Nucl. Phys. , 260 (1980),[Yad. Fiz.31,497(1980)].[50] M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B104 , 199 (1981).[51] J. Gasser and H. Leutwyler, Annals Phys. , 142(1984).[52] S. Scherer, Adv. Nucl. Phys. , 277 (2003),arXiv:hep-ph/0210398.[53] J. Wess and B. Zumino, Phys. Lett. B , 95 (1971).[54] E. Witten, Nucl. Phys. B , 422 (1983).[55] H. Lehmann, K. Symanzik, and W. Zimmermann,Nuovo Cim. , 205 (1955).[56] R. Z. Ferreira, A. Notari, and F. Rompineve,(2020), arXiv:2012.06566 [hep-ph].[57] J. Gasser and H. Leutwyler, Phys. Lett. B , 83(1987).[58] J. Gasser and H. Leutwyler, Phys. Lett. B , 477(1987).[59] P. Gerber and H. Leutwyler, Nucl. Phys. B ,387 (1989).[60] S. Hannestad and J. Madsen, Phys. Rev. D , 1764(1995), arXiv:astro-ph/9506015.[61] G. Colangelo, J. Gasser, and H. Leutwyler, Nucl.Phys. B , 125 (2001), arXiv:hep-ph/0103088.[62] S. Aoki et al. (Flavour Lattice Averaging Group),Eur. Phys. J. C , 113 (2020), arXiv:1902.08191[hep-lat]. [63] L. Darm´e, L. Di Luzio, M. Giannotti, and E. Nardi,(2020), arXiv:2010.15846 [hep-ph].[64] S. Weinberg, Phys. Rev. Lett. , 616 (1966).[65] U. Aydemir, M. M. Anber, and J. F. Donoghue,Phys. Rev. D , 014025 (2012), arXiv:1203.5153[hep-ph].[66] T. Blum et al. , Phys. Rev. D , 074502 (2015),arXiv:1502.00263 [hep-lat].[67] R. Abbott et al. (RBC, UKQCD), Phys. Rev. D , 054509 (2020), arXiv:2004.09440 [hep-lat].[68] M. Luscher, Nucl. Phys. B , 531 (1991).[69] K. Rummukainen and S. A. Gottlieb, Nucl. Phys.B , 397 (1995), arXiv:hep-lat/9503028.[70] C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl.Phys. B , 218 (2005), arXiv:hep-lat/0507006.[71] M. T. Hansen and S. R. Sharpe, Phys. Rev. D ,016007 (2012), arXiv:1204.0826 [hep-lat].[72] R. A. Briceno, J. J. Dudek, R. G. Edwards, andD. J. Wilson, Phys. Rev. Lett. , 022002 (2017),arXiv:1607.05900 [hep-ph].[73] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr,and B. Fuks, Comput. Phys. Commun. , 2250(2014), arXiv:1310.1921 [hep-ph].[74] N. D. Christensen and C. Duhr, Comput. Phys.Commun. , 1614 (2009), arXiv:0806.4194 [hep-ph].[75] T. Hahn, Comput. Phys. Commun. , 418 (2001),arXiv:hep-ph/0012260.[76] V. Shtabovenko, R. Mertig, and F. Orellana,Comput. Phys. Commun. , 107478 (2020),arXiv:2001.04407 [hep-ph].[77] V. Shtabovenko, R. Mertig, and F. Orel-lana, Comput. Phys. Commun. , 432 (2016),arXiv:1601.01167 [hep-ph].[78] R. Mertig, M. Bohm, and A. Denner, Comput.Phys. Commun. , 345 (1991).[79] H. H. Patel, Comput. Phys. Commun. , 276(2015), arXiv:1503.01469 [hep-ph]. SUPPLEMENTARY MATERIAL
The calculation of the amplitudes was carried out using the computational tools
FeynRules [73, 74],
FeynArts [75],
FeynCalc [76–78] and
Package-X [79]. The full analytical expression of the renormalizedNLO amplitude for the aπ → π + π − process reads M NLO aπ → π + π − = C aπ π f π f a (cid:40) m π ( u + t ) − u − ut − t − (cid:96) (cid:0) m π − s (cid:1) (cid:0) m π − s (cid:1) − (cid:96) (cid:0) − m π ( u + t ) + 4 m π + u + t (cid:1) + 9 m π (cid:96) + 18 (cid:96) m π ( m π − s ) + 576 π (cid:96) m π (cid:18) m d − m u m d + m u (cid:19) + 3 (cid:34) (cid:114) − m π s s (cid:0) m π − s (cid:1) ln (cid:32) (cid:112) s ( s − m π ) + 2 m π − s m π (cid:33) + (cid:114) − m π t (cid:0) m π ( t − u ) + 3 m π + t ( u − t ) (cid:1) ln (cid:32) (cid:112) t ( t − m π ) + 2 m π − t m π (cid:33) + (cid:114) − m π u (cid:0) m π ( u − t ) + 3 m π + u ( t − u ) (cid:1) ln (cid:32) (cid:112) u ( u − m π ) + 2 m π − u m π (cid:33)(cid:35) (cid:41) + 4 (cid:96) m π m d (cid:0) s − m π (cid:1) m u ( m d − m u ) f π f a ( m d + m u ) , (16)where the terms proportional to (cid:96) , (cid:96) and (cid:96) in the second row arise from the LO amplitude, via theNLO corrections to m π and f π (see Ref. [51]). The latter include the charged-neutral pion mass differencearising at second order in the isospin breaking parameter m d − m u , which has been neglected in the numericalintegration of the rate. The amplitudes for the crossed channels aπ − → π π − and aπ + → π + π are obtainedby cross symmetry through the replacements s ↔ t and s ↔ u , respectively.Next, we describe our procedure to analytically reduce the 12-dimensional phase space integral of Eq. (14)down to a 5-dimensional one. We first integrate the fourth-particle phase space in Eq. (14) using therelation d p / (2 E ) = d p δ (cid:0) p − m (cid:1) θ ( p ). Therefore, defining the angles α and θ via cos α = p · p | p || p | andcos θ = p · p | p || p | , the thermalization rate becomesΓ a = 1 n eq a (cid:90) dp | p | E dp | p | E dp | p | E (cid:90) − d cos α (cid:90) − d cos θ (cid:90) π dβ (cid:88) |M| π (2 π ) × δ ( E − ξ ( E , E , α, β, θ ))2 | E − E − | p | cos α + | p | cos θ | f f (1 + f )(1 + f ) , (17)with β the angle between the scattering planes defined by ( p , p ) and ( p , p ), and the function ξ given by ξ ( E , E , α, β, θ ) = 2 ( E E − | p || p | (sin α sin θ cos β + cos α cos θ )) − m π E − E − | p | cos α + | p | cos θ ) ..