Axion-pion scattering at finite temperature in chiral perturbation theory and its influence in axion thermalization

Axion-pion scattering amplitudes at finite temperatures are calculated within chiral perturbation theory up to the one loop level. Unitarization procedure is implemented to these amplitudes in order to extend the applicable range of energy and temperature. The influence of the thermal axion-pion scattering amplitudes on the $a\pi\to\pi\pi$ cross sections and the axion thermalization rate is investigated, with the emphasis on the comparison with the zero-temperature-amplitude case. A brief discussion on the cosmological implication of the axion thermalization rate, that is calculated by using the $a\pi\to\pi\pi$ amplitudes at finite temperatures, is also given. The thermal corrections to the axion-pion scattering amplitudes can cause around a $10\%$ shift of the determination of the axion decay constant $f_a$ and its mass $m_a$, comparing with the results by using the $a\pi\to\pi\pi$ amplitudes at zero temperature.


Introduction
The conjectural particle axion constitutes a dynamical research subject in many areas of physics, including particle physics, cosmology, astronomy, etc [1][2][3][4][5][6].The primary aim of the QCD axion is to solve the strong CP problem, which can be elegantly addressed within the Peccei-Quinn (PQ) mechanism [7][8][9].The axion is conjectured to be a Nambu-Goldstone boson that arises from the spontaneous breaking of a global U (1) P Q symmetry at some high energy scale f a , also known as the axion decay constant.The key idea of the PQ mechanism is to promote the anomalous term θG G in QCD, being G and G the gluon filed tensor and its dual, by the dynamical axion field a as a fa G G, which softly breaks the U (1) P Q symmetry and makes the axion actually a pesudo-Nambu-Goldstone boson (pNGB).In fact the anomalous a fa G G term is commonly regarded as the model independent part in the construction of various axion models.Being the essential QCD nature, axion inevitably interacts with the hadrons and intensive investigations are carried out to find evidence for such interactions [4].
One attractive way to trace the axion imprint is to study the axion production in the thermal medium of the early Universe.Below the critical temperature (T c ∼ 150 MeV) of QCD phase transition, the quarks and gluons are confined inside hadrons and consequently one needs the axion-hadron interactions as inputs to reliably calculate the axion production in the thermal bath.As the lightest QCD hadron, the pions are expected to provide the most important sources for the axion production below T c .Indeed, continuous efforts have been made to refine the axion-pion interactions in recent years, in order to improve the cosmological constraints on the axion properties.In a pioneer work [10], the axion-pion scattering amplitude is calculated by taking the leading-order (LO) chiral perturbation theory (χPT).The LO χPT axion-pion amplitude is adopted in many sequential works to constrain the axion decay constant f a or its mass m a [11][12][13][14][15]. Recently, the calculation of the axion-pion scattering amplitude is pursued up to the next-to-leading order (NLO) in χPT, and an important finding is that the perturbative χPT results become unreliable for the temperature T > 62 MeV, due to the loss of convergence of the chiral expansion for the axion thermalization rate [16].Later on, different strategies are proposed to extend the applicable region of the χPT amplitudes in Refs.[17,18].The inverse amplitude method (IAM) is employed in Ref. [18] to unitarize the NLO axion-pion amplitude.While, in Ref. [17] it is assumed that the axion-pion amplitudes can be approximated by scaling the pion-pion amplitudes with the axion-pion mixing strength factor.It is pointed out that in all the previous works the axion thermalization rates are calculated by taking the axion-pion scattering amplitudes at zero temperature, instead of the thermal ones, as an approximation.In this work, we fill this gap by performing the aπ → ππ scattering amplitudes at finite temperatures in χPT up to the one-loop level.Then we briefly study the influence of the thermal aπ → ππ amplitudes on the axion thermalization rates and the determinations of the axion parameters by confronting the extra effective number of relativistic thermal degrees of freedom ∆N eff .
This article is organized as follows.In Sec. 2, we elaborate the relevant axion chiral Lagrangian and the calculation of aπ → ππ scattering amplitudes at finite temperatures up to the one-loop level.A survey of thermal corrections to the scattering amplitudes will be also given in this section, with the emphasis on comparing with the zero-temperature amplitudes used in literature.In Sec. 3, a brief discussion about the cosmological constraints on the axion parameters relying on the averaged axion thermalization rates calculated from the thermal aπ → ππ scattering amplitudes will be presented.We give a short summary and conclusions in Sec. 4. Technical details about the calculations of thermal loops and phase space integrals are given in Appendices.
2 Axion-pion scattering amplitudes at finite temperatures

Relevant axion χPT Lagrangians
In this work we stick to the minimal QCD axion Lagrangian with the dual of the gluon tensor G µν c ≡ ε µνρσ G c ρσ /2.Various extensions of the axion interactions with quarks, leptons and photons can be found in many recent reviews [1][2][3][4][5][6].By performing the quark field transformation q → e i a 2fa Qaγ 5 q, being Q a a 2 × 2 hermite matrix in quark flavor space with Tr(Q a ) = 1, one can eliminate the anomalous term aG G and in the meantime also introduce additional axial-vector coupling and a modification to the quark mass term as where M q corresponds to the two-flavor quark mass matrix M q = diag(m u , m d ) .
Based on these ingredients, the axion χPT Lagrangian can be constructed order by order as the usual χPT [19,20] The LO axion χPT Lagrangian at O(p 2 ) in the SU (2) reads where F corresponds to the pion decay constant at LO, χ a = 2B 0 M q (a), and the axion dressed quark mass matrix M q (a) is given by The pion fields are collected in J µ A LO stands for the axial-vector current −qγ µ γ 5 Q a q in Eq. ( 2) at the hadron level from the LO χPT In literature sometimes different assignments for the chiral transformation behaviors of the chiral building blocks, especially the axion dressed quark mass terms, are assumed in the calculation.Here we follow the conventions of the seminal χPT works in Refs.[19,20] for the chiral transformation behaviors of the various building blocks: (so that the QCD quark mass term −q R (s + ip)q L + h.c. is chiral invariant), and , the latter of which implies that the covariant derivative takes the form D µ U = ∂ µ U − ir µ U + iU l µ , with the relations between the right/left-hand (r µ /l µ ) and vector/axial-vector (v µ /a µ ) external sources: r µ = v µ + a µ and l µ = v µ − a µ .The sign conventions introduced in Eq. ( 6) for the axion dressed quark mass and Eq. ( 8) for the axial-vector current are consistent with the just mentioned definitions of chiral transformation behaviors used in Refs.[19,20].It is mentioned that there is a minus sign difference for the definition of hadronic axial-vector currents between ours in Eq. ( 8) and the one in Ref. [18].
Regarding the choice of Q a , it is noticed in Ref. [21] that by taking the mass mixing between axion and pion will be automatically eliminated at LO.Nevertheless, the last term in Eq. ( 5) will lead to the kinematic mixing for axion and pion.To be definite in our calculation, we will take the form of Q a in Eq. ( 9) throughout, although the physical quantities should remain the same with different choices of Q a .The axion-pion scattering amplitudes receive the NLO contributions at O(p 4 ) both from the loops by taking the LO vertices in Eq. ( 5) and the O(p 4 ) local chiral operators [19,22] where the NLO piece of the axial-vector current is Notice that only the terms relevant to this work are kept in Eqs.(10) and (11).
We use the imaginary time approach to include the finite-temperature effects [23][24][25].Within the framework of χPT [26,27], the thermal corrections only enter via the chiral loops, while the low energy constants (LECs) l i and h i accompanying the local operators in Eqs.(10) and (11) do not depend on the temperatures.This in turn implies that once the values of the unknown LECs are fixed at zero temperature one can make pure predictions at finite temperatures.In Appendix A, we give relevant formulas for the one-loop integrals at finite temperatures.

Calculation of axion-pion scattering amplitudes at finite temperatures up to one loop
For the calculation of the amplitudes, one needs to address the LO a-π 0 mixing first.After taking the specific form of Q a in Eq. ( 9), the a-π 0 mixing at LO is exclusively caused by the ∂ µ aJ µ A LO term in Eq. ( 5), which leads to δ aπ ∂ µ a∂ µ π 0 with By taking the field redefinition: , one can eliminate the mixing term of axion and π 0 at LO and in the meantime render the coefficients of their kinetic terms as 1/2 at the level of O(1/f a ) (see Appendix B).Nevertheless, the a-π 0 mixing will appear again at NLO in the chiral expansion, which leads to a non-diagonal twopoint functions G ij with i, j = a, π 0 .In Appendix B, we give a detailed discussion about the two-point functions up to NLO, along with the temperature dependence of the axion and pion masses.The terms of O(1/f 2 a ) will be neglected in the calculation of the axion-pion scattering amplitudes throughout.
The axion-pion scattering amplitudes can be extracted from the four-point Green functions by using the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula, as done in Ref. [16].
For the a(p 1 )π(p 2 ) → π(p 3 )π(p 4 ) scattering amplitude at finite temperatures, it can be decomposed into two parts: the zero temperature part and the thermal correction one, The amplitude M (T =0) aπ;ππ at zero temperature up to one loop has been computed and then plugged into the phase space integrals weighted by the Bose-Einstein (BE) distribution factors to estimate the axion thermal rates in Refs.[16,18].As an improvement, we calculate the thermal correction part ∆M (T ) aπ;ππ up to one loop in this work, and its full expressions are given in Appendix C. It is noted that at finite temperatures the commonly adopted Mandelstam variables are not enough to describe the thermal amplitudes ∆M (T ) aπ;ππ , due to the loss of the Lorentz invariance at finite temperature T .As shown later in Appendix C it is convenient to introduce the four momenta p s , p t , p u , corresponding to s, t, u, respectively, to describe the thermal amplitudes [28] p In general the thermal amplitude can depend on the temporal and spatial components of the four momenta in a separate manner, as can be seen in Appendix C. Due to the rotation invariance, the number of independent kinematic variables in the finite-temperature amplitudes is five, and this number is also the effective dimension of the phase space integral in the evaluation of axion rate.In Appendix D we will illustrate in detail about the five specific kinematic variables to calculate the phase space integral.Due to the momentum expansion feature of χPT, its amplitudes are not expected to be convergent well at relatively high energy regions, especially where the resonances can appear.In Ref. [16], it is demonstrated that the axion thermalization rates calculated from the NLO χPT amplitudes are questioned to be valid above T > 62 MeV.One efficient way to extend the applicable regions of energy and temperature for the perturbative χPT amplitudes is to restore the unitarity relations of the partial-wave (PW) amplitudes.We will use the inverse amplitude Figure 1: Feynman diagrams of the amputated four-point Green functions in Eqs. ( 13) and ( 14).The dashed line stands for the axion field and the solid line corresponds to the pion field.
The number n for each vertex denotes its chiral order of O(p n ).Diagrams 1(a)-1(f) denote A (2,4) aπ;ππ , and diagram 1(g) corresponds to A ππ;ππ which will be multiplied by the NLO a-π 0 mixing to give one type of contributions to the aπ → ππ scattering amplitudes at O(p 4 ), i.e. the last terms in Eqs. ( 13) and (14).
At finite temperature T , we will work in the center of mass (CM) frame to perform the PW expansion in order to proceed the unitarization procedure, with P J (cos θ) the Legendre polynomials.In the CM frame all the kinematic variables appearing in the thermal amplitudes can be related to the energy E cm and the scattering angle θ defined in that frame.The PW amplitudes with definite angular momentum J can be obtained via Eq.( 22) by performing the PW integrals We follow the IAM procedure to construct the unitarized PW aπ → ππ amplitude at finite temperature with and the BE distribution factor In the zero-temperature limit, the BE distribution factor n B (E) tends to vanish and the thermal unitarity relation in Eq. ( 25) recovers the standard one at zero temperature given in Ref. [18].A subtle difference comes up in the finite-temperature case that additional thermal Landau cuts generally appear [33][34][35].E.g., in the Kπ → Kπ scattering process, the thermal Landau cuts generated in the u channel show up in the physical energy region above the Kπ threshold [35].
In the current study of the aπ → ππ process, apart from the thermal unitarity cuts shown in Eq. ( 25), there are also thermal Landau cuts in the energy region E cm > 2m π contributed from both the t and u channels, which magnitudes turn out to be generally much smaller than the unitarity cuts.
Figure 3: Cross sections in the charged bases.For the notaions of σ LO,NLO,NLO',IAM , see the text for details.Notice that the LO contributions to the aπ 0 → π 0 π 0 amplitude and cross section vanish.

Surveys of the thermal corrections to the axion-pion scattering amplitudes and cross sections
In this part, we examine to what extent the finite-temperature corrections can affect the axion-pion scattering amplitudes and cross sections, comparing with the zero-temperature results.To make close comparison with the zero-temperature study of Ref. [18], we will take the same inputs for the O(p 4 ) LECs and other parameters as those used in the former reference, i.e., in perturbative calculation we take l1 = −0.36(59),l2 = 4.31 (11), l3 = 3.53 (26), l4 = 4.73 (10), l 7 = 0.007(4), together with m u /m d = 0.50(2), F π = 92.1 (8) MeV and m π = 137 MeV; while for the IAM results we take the combinations l1 − l2 = −5.95(2),l1 + l2 = 4.9 (6), and others are same as the perturbative case.
In Fig. 2, we give the magnitudes of different PW amplitudes at T = 0, 100 and 155 MeV.The first lesson we can learn from this figure is that the thermal corrections to the perturbative NLO amplitudes are insignificant up to T = 155 MeV.The largest thermal effects appear in the IAM amplitudes for the IJ = 11 channel, while the thermal corrections to the IAM amplitudes in other channels look small up to T = 155 MeV.
We then calculate the cross sections for the aπ → ππ in charged bases, both at zero and finite temperatures, where |⃗ p i | and |⃗ p f | are the magnitudes of three-momenta of the initial and final particles in the CM frame, and N = 1/2 for identical particles and otherwise N = 1.The comparisons are shown in Fig. 3.For the evaluation of the cross sections by using the perturbative amplitudes, one can expand the amplitude squared as where the property that the LO M (2) is real has been used.In later discussions, we will denote σ LO as the LO cross section when taking |M (2) | 2 in Eq. ( 29), and designate σ NLO when taking the first two terms of Eq. ( 29), and σ NLO ′ when taking all the three terms of Eq. ( 29) to evaluate the cross sections in Eq. (28).When using the unitarized IAM amplitude M IAM in Eq. ( 24) to calculate the cross sections, the result is simply denoted as σ IAM .One can gain a rough estimate about how the next-to-next-to-leading order perturbative corrections may affect the cross sections by comparing σ NLO with σ NLO ′ at low energies.The curves corresponding to the IAM amplitudes in Figs. 2 and 3 are considered to be the preferred results in this work.By taking a closer look at different curves, one can gain a rough conclusion that the perturbative NLO amplitudes begin to be unreliable around E ≃ 0.5 GeV, above which the next-to-next-to-leading order and IAM results begin to deviate noticeably.

Brief discussions on the axion thermalization rate and its cosmology implication
We will closely follow Refs.[16,18] to proceed a brief discussion of the cosmological constraint on the axion parameters based on our updated calculation of the thermal aπ → ππ amplitudes.The axion decoupling temperature T D plays the key role, and according to the proposal in Ref. [11] it can be estimated as the temperature at which the axion interaction rate turns to be lower than the expansion rate of the universe.The explicit freeze out condition advocated in Ref. [11] to determine the axion decoupling temperature is where H(T ) = T 2 4π 3 g * (T )/45/m Pl denotes the Hubble expansion parameter, with g * (T ) the effective number of relativistic thermal degrees of freedom (d.o.f) and m Pl the Plank mass.Γ a (T ) in Eq. ( 30) stands for the averaged axion thermalization rate in the pion thermal bath, and it can be calculated via where the sum runs over all the possible aπ reaction channels, the axion number density in equilibrium is n eq a = ζ(3) T 3 /π 2 with ζ the Riemann Zeta function, and the phase space integral is The recipe given in Ref [36] will be used to calculate the above phase space integral.For the sake of completeness, we give some details about the evaluation of this integral in Appendix D.
A more sophisticated way to solve the momentum dependent Boltzmann equations for the axion-pion scatterings, rather than taking the criteria in Eq. ( 30), is proposed in Refs.[17,37].However, the primary aim of this work is to examine to what extent the thermal corrections to the aπ → ππ scattering amplitudes can influence the determinations of the axion parameters, comparing with the situation by using the amplitudes at zero temperature, and therefore we will stick to the criteria in Eq. ( 30) to proceed.
For the perturbative calculation up to NLO, we expand the squared amplitude as |M| 2 ≃ |M (2) | 2 + 2M (2) Re(M (4) ), and then the perturbative axion rate will be casted into the following form [11] Γ pert a (T ) = where the dimensionless functions h lo and h nlo denote the contributions from |M (2) | 2 and 2Re(M (2) M (4) ), respectively.The numeric factors hlo and hnlo are introduced so that h lo (m π /T ) and h nlo (m π /T ) are separately normalized at unity at the temperature T c , which will be set as 155 MeV in order to make close comparisons with the results in Refs.[16,18].The values of these factors are: hlo ≃ 0.164, hnlo ≃ −0.0550 (when using zero-temperature amplitudes) and −0.0755 (when using thermal amplitudes).When using IAM amplitudes as inputs, the rate is simply expressed as For comparison, h IAM (m π /T ) is also normalized to unity at T c = 155 MeV by the numeric factor hIAM ≃ 0.117 (when using zero-temperature amplitudes) and 0.0817 (when using thermal amplitudes).The results of h lo,nlo,IAM (m π /T ) are shown in Fig. 4(a).
For the axion thermalization rates estimated from the perturbative χPT amplitudes, the corresponding curves are given in Fig. 4(b).After considering the NLO corrections of the amplitudes, the axion rate even drops to a negative value at high temperatures, due to the extreme breakdown of the chiral expansion in these regions.In Refs.[16,18], it is proposed to inspect the valid region of χPT by using the ratios of the NLO and LO parts of the thermalization rates as a criteria.Our results for such ratios as a function of T are given in Fig. 4(c).The breakdown temperature of χPT was extracted from the maximun of the ratio for the aπ + channel [18].By approximating the aπ scattering amplitude from its zero-temperature expression, we confirm that the χPT breakdown temperature deduced from the aπ + channel is around 70 MeV, and at this point the ratio reaches the maximum value that is around 60%, see the blue solid line in Fig. 4(c).For the result after taking the thermal correction to the aπ scattering amplitude, i.e. the blue dashed line, the maximum point of the ratio for aπ + channel increases to around 85 MeV, with the maximum ratio around 65%.The dashed lines in Fig. 4(c) are all above the solid lines, thus the convergence of the chiral expansion for axion thermalization rate gets worse after taking the thermal corrections into the aπ scattering amplitudes.Therefore after including the thermal corrections to the aπ → ππ amplitudes it is still unlikely to extract reliable bounds from the axion thermalization rates by using the perturbative χPT amplitudes, validating the previous conclusion in Refs.[16,18] that was obtained by taking the zero-temperature aπ → ππ amplitudes.
The axion averaged thermalization rates as a function of T from the IAM amplitudes are illustrated in Fig. 5(a), where two groups of lines corresponding to two specific values By ZT M IAM , f a =2.510 7 GeV By FT M IAM , f a =2.510 7 GeV By ZT M IAM , f a =610 6 GeV of f a are given.We then compare the axion rates by scanning the values of f a with the Hubble parameters by taking the cosmological inputs from Ref. [38].According to the criteria in Eq. ( 30), the interception point of the two curves of Γ a (T ) and H(T ) correspond to the axion decoupling temperature T D , which varies with the value of f a .In Fig. 5(b), we show the relations between T D and f a by taking the scattering amplitudes both at zero and finite temperatures to evaluate Γ a (T ).The effects from the thermal corrections to the scattering amplitudes become noticeably important at higher values of f a .
To take a comprehensive study of the cosmological analysis goes beyond the scope of this work, here we will concentrate on exploring the constraint of the extra effective number of relativistic d.o.f: on the axion parameters as done in Refs.[16,18], where g ⋆s (T ) corresponds to the effective number of entropy d.o.f at temperature T , and the cosmological determinations in Ref. [38] will be used in our analysis.The constraint of ∆N eff on the axion decay constant f a is given in Fig. 5(c), where the gray area in the bottom corresponds to the region for T > 155 MeV and the IAM χPT is considered to be untrustworthy in this region.Our study indicates that the inclusion of the thermal corrections to the scattering amplitudes can slightly lower the bound f a ≳ 2.1 × 10 7 GeV, comparing with the constraint f a ≳ 2.3 × 10 7 GeV obtained by taking the aπ → ππ scatering amplitudes at zero temperature.That is to say the thermal corrections to the scattering amplitudes can cause around a 10% shift to the constraint of the axion decay constant f a .
In turn one could translate the bounds of f a into the constraints of the axion mass m a , relying on the LO/NLO χPT predictions in Appendix B. The constraints on the QCD axion mass m a are shown in Fig. 5(d), where we have distinguished the constraints by including/excluding the thermal corrections to the aπ scattering amplitudes and also using the LO/NLO expressions for the axion masses.When excluding the thermal corrections of aπ scattering amplitudes, the bounds for the axion masses are m a ≲ 0.24 eV by m LO a and m a ≲ 0.25 eV by m NLO a .When including the thermal corrections of aπ scattering amplitudes, the bounds for the axion mass are m a ≲ 0.27 eV by m LO a and m a ≲ 0.28 eV by m NLO a .Therefore the thermal corrections of the amplitudes introduce a 10% variation to the determination of the upper limits of m a .

Summary and conclusions
In this work we have performed the complete calculation of the aπ → ππ scattering amplitudes up to the one-loop level at finite temperatures within the SU (2) QCD axion chiral perturbation theory.The inverse amplitude method is further used to unitarize the partialwave aπ → ππ scattering amplitudes, which extends the applicable energy and temperature regions of the amplitudes in chiral perturbation theory.The largest effect from the thermal correction to the aπ → ππ scattering amplitudes appears in the aπ ± → π ± π 0 channels with IJ = 11 for the ππ system, while the thermal effects in other scattering amplitudes are small.
The axion averaged thermalization rates are then calculated by using the updated aπ → ππ scattering amplitudes at finite temperatures.The axion decoupling temperatures are obtained with the averaged axion thermalization rates calculated from the thermal aπ → ππ scattering amplitudes within the inverse amplitude method.The extra effective number of relativistic degrees of freedom ∆N eff is then used to constrain the axion decay constant f a , which gives f a > 2.1 × 10 7 GeV, comparing with the bound of f a > 2.3 × 10 7 GeV when neglecting the thermal corrections in the aπ scattering amplitudes.This 10% correction effect for the lower limit of f a also delivers a similar magnitude of shift to the upper limit of the axion mass.Therefore our study shows that the inclusion of the thermal effects in the aπ → ππ scattering amplitudes can cause a mild shift (around 10%) of the axion parameters when confronting with the cosmological constraint of ∆N eff .It is necessary to include the thermal corrections in the aπ scattering amplitudes in the future attempts that aim to improve the cosmological determination of the axion parameters at this level.

A Relevant one-loop integrals at finite temperatures
In the imaginary time formalism, the finite-temperature effects enter via the chiral loops by replacing the temporal integrals with the Matsubara sums, i.e., where the replacement of q 0 → iω n = i2πnT with n ∈ Z is taken and the symbol β is introduced to highlight the integrals with finite-temperature corrections.Meanwhile, the zeroth components of all external momenta also need to be replaced by the Matsubara frequencies, which can be extrapolated to real energies after completing the Matsubara sums.When performing such sums, one can always separate the results into the zero-temperature part and the finite-temperature correction part: where the first term corresponds to the zero-temperature integral and the second one stands for the finite-temperature correction.One can use the standard dimensional regularization to calculate the zero-temperature integral, in which the UV divergence will appear.We use the conventional M S − 1 renormalization scheme that is widely employed in χPT [19] for the zero-temperature integrals G(T = 0).Next we give relevant formulas for the one-loop integrals at finite temperatures, i.e. the ∆G(T ) part.The finite-temperature integrals in this work are the same as those in thermal ππ scattering [28], which include the tadpole integrals, like the diagram 1(c), and the one-loop two-point integrals for diagram 1(d)-1(f).

A.1 Tadpole loop integral
According to Eq. (A.2), the basic tadpole integral reads where the zero-temperature integral in the M S − 1 scheme is and the finite-temperature correction is For the integrals with p • q or q 2 in the numerators, the results are A.2 Basic two-point one-loop integrals at finite temperatures Four types of the two-point one-loop integrals will appear in the calculation of χPT whose zero-temperature parts can be given in terms of the due to the Lorentz invariance.In the following, we refrain from discussing the zero-temperature parts of the integrals in Eq. (A.8) and focus on the finite-temperature correction (FT C).While for the latter part, the Lorentz invariance is lost and one needs to separately treat the temporal and spatial components of the integration variable q µ = (q 0 , q i ) in the evaluation.It turns out that the FT C parts of the four types of integrals in Eq. (A.8) can be expressed in terms of three independent thermal loop functions [28] ∆J l (T, k 0 , | ⃗ k|) ≡ −i where the subscript FT C is introduced to emphasize that ∆J l (T, k 0 , | ⃗ k|) only include the finite-temperature parts of the integrals.After completing the Matsubara sums, the integrals in Eq. (A.10) can be written as depend on k 0 and | ⃗ k| in a separate manner.Before giving the technical details about the calculation of the thermal loop functions in Eq. (A.12), we first elaborate the reduction formulas for the temperature correction parts of the two-point one-loop integrals in Eq. (A.8).

A.3 Reduction of the two-point one-loop integrals at finite temperatures
The numerators with spatial components q i in the integrals of Eq. (A.8) can be written in terms of the basic ones in Eqs.(A.3) and (A.10).For the sake of completeness, we give a brief discussion about the finite-temperature parts of the various integrals in Eq. (A.8).
For the three types of numerators with (p 1 •q), q 2 , (p 1 •q)(p 2 •q) in the integrals of Eq. (A.8), their FT C parts can be reduced to ∆F β (T ) (A.3) and ∆J l (T, k 0 , | ⃗ k|) (A.10).The FT C reduction formulas take the form It is reiterated that the above reduction formulas with the subscript of FT C include only the finite-temperature correction parts.

A.4 Technical details about the evaluation of basic thermal loop functions
We give some technical details below about how to evaluate the integrals in Eq. (A.12).It is noted that f (l) λ 1 λ 2 (T, k 0 , | ⃗ k|) satisfy the following reflection relations with respect to k 0 [34]  Next we use a similar procedure in Ref. [34] to simplify the expressions of f (l) given in Eq. (A.12).To proceed the evaluation of the integrals in Eq. (A.12), we take E q = ⃗ q 2 + m 2 and E k−q = ⃗ k 2 + ⃗ q 2 − 2| ⃗ k||⃗ q | cos θ + m 2 as the two integral variables, i.e., The real and imaginary parts of the f (l) 28) where P denotes taking principle value of the integrals.The δ functions in Imf E max k−q , and to combine with k 2 > 4m 2 one can get Therefore, we have The remaining integration of E q in Imf (l) ) is straightforward and can be even analytically evaluated , nevertheless we do not explicitly show the results here, due to the rather lengthy expressions.

B Two-point functions and masses of axion and pion at finite temperatures
The bilinear terms of a and π 0 from the LO Lagrangian (5) read where δ aπ is given in Eq. ( 12) and the QCD axion mass at LO is The barred quantities correspond to their LO expressions.Up to O(1/f 2 a ) , the a-π 0 mixing in Eq. (B.1) can be eliminated by the field redefinition and the Lagrangian in Eq. (B.1) becomes The correction term δ 2 aπ in the pion mass behaves O(1/f 4 a ) and can be safely ignored.Then the LO propagators of a and π read With contributions from the various two-point 1PI amplitudes Σ ij , the two-point Green functions of a and π take the form Up to NLO in chiral expansion, the explicit expressions of the matrix elements of G and Σ read  where Σ (4)  aa and the superscript (4) denotes the O(p 4 ) contribution.In this work, the isospin breaking correction in m π is neglected, and therefore we have In the calculation of the axion-pion scattering amplitudes one will need G aπ 0 (m 2 a ) to account for the a-π 0 mixing at NLO. Ignoring the terms of O(1/f 2 a ) and pulling out the axion pole in Eq. (B.11), by keeping the terms up to NLO one has Up to NLO in χPT, the ratios of the pion and axion masses at finite and zero temperatures can be written as where the masses at zero temperature read In Eq. (B.19) we have replaced m2 π with the "physical" pion mass m 2 π up to NLO.The results are consistent with Ref. [39] by ignoring δ I term in m 2 π .When calculating the boundary of the axion mass, the renormalization scale µ is set at 770 MeV, and h 1 = h 3 = 0 is taken.Our result confirms the temperature dempendence of the QCD axion mass in the former reference as well.The curves of m a (T )/m a (T = 0) and m π (T )/m π (T = 0) are given in Fig. 6.As shown in this figure, the thermal correction decreases the axion mass around 10% up to T = 155 MeV, while it slightly increases the pion mass at the level around 3%.

C Explicit expressions of the aπ → ππ scattering amplitudes at finite temperatures
Here we give the explicit one-loop expressions of the aπ 0 → π + π − and aπ 0 → π 0 π 0 scattering amplitudes at finite temperatures.For completeness we also give the zero-temperature parts of the scattering amplitudes, although they have been given in Ref. [18].The finitetemperature corrections to the aπ → ππ scattering amplitudes represent the new ingredients of this work.In this work, we will use in the scattering amplitudes the LO pion decay constant F , which is independent of the temperature T .Up to NLO, F is related to the physical pion decay constant F π via F = F π 1 − m 2 π (4πFπ) 2 l4 [19], which leads to F = 86.0MeV by taking the same value of l4 = 4.73 from Ref. [18].
The NLO amplitude is decomposed into the zero-temperature part and finite-temperature correction part.The zero-temperature part is [18] M (4), (T =0) The integrals over cos θ 1,3 can be performed by the remaining two δ functions.Using δ (g(x)) = i δ(x − xi )/|g ′ (x i )| with xi the root of g(x) , we have aπ are the two-point one-particle-irreducible (1PI) amplitudes whose expressions can be found in Appendix B, and Σ (4) ′ ππ stands for the derivative of Σ (4) ππ .The superscripts (2) and (4) denote the chiral orders.The corresponding Feynman diagrams for the amputated amplitudes A (2,4) aπ;ππ and A (2)

Figure 2 :
Figure 2: Magnitudes of the various aπ → ππ partial-wave amplitudes in the isospin bases.The curves labeled as NLO include both the contributions at O(p 2 ) and O(p 4 ).

Figure 4 :
Figure 4: The symbols of ZT and FT indicate that the results are calculated by taking the zero-temperature and finite-temperature scattering amplitudes, respectively.(a) The curves for the h lo, nlo, IAM (m π /T ) functions are defined in Eqs.(33) and (34).(b) The axion rates in this figure are calculated by taking the perturbative aπ → ππ scattering amplitudes upto NLO, with f a = 1.5 × 10 7 GeV for illustration.Γ LO a only includes the contribution from h lo , and Γ NLO a contains the contribution both from h lo and h nlo .(c) Ratio between the NLO parts and LO parts of the axion thermal averaged rates.Γ lo a denotes the LO part of the perturbative rate and Γ nlo a denotes the part of perturbative rate contributed only by h nlo .

Figure 5 :
Figure 5: The symbols of ZT and FT indicate that the results are calculated by taking the zero-temperature and finite-temperature scattering amplitudes, respectively.(a) Axion thermalization rates calculated with the IAM amplitudes as a function of T .The solid lines denote the results by taking the amplitudes at zero temperature as an approximation and the dashed lines correspond to the updated results by including the thermal corrections to the aπ → ππ amplitudes.(b) Axion decoupling temperature as a function of the axion decay constant f a .(c) Extra effective thermal relativistic d.o.f. as a function of f a .(d) The bounds on the axion masses m a .The notations of m LO a and m NLO a refer to the situations by using the LO/NLO expressions for the axion masses.

Figure 6 :
Figure 6: The ratios of the axion and pion masses at finite and zero temperatures predicted by the NLO χPT.