Thermal axion production

We reconsider thermal production of axions in the early universe, including axion couplings to all Standard Model (SM) particles. Concerning the axion coupling to gluons, we find that thermal effects enhance the axion production rate by a factor of few with respect to previous computations performed in the limit of small strong gauge coupling. Furthermore, we find that the top Yukawa coupling induces a much larger axion production rate, unless the axion couples to SM particles only via anomalies.


Introduction
The strong CP problem can be solved by a Peccei-Quinn (PQ) symmetry [1], that manifests at low energy as a light axion a.The axion is a good dark matter (DM) candidate, if cold axions are produced non-thermally via the initial misalignment mechanism [2].The cosmological DM abundance is reproduced for an order one initial misalignment angle provided that f a ≈ 10 11 GeV, which is compatible with the experimental bound f a > ∼ 5×10 9 GeV [3,4].The ADMX experiment can probe such scenario in the next years [5].
Furthermore, thermal scatterings in the early universe unavoidably produce a population of hot axions.The goal of this paper is performing an improved computation of hot axion thermal production.The thermal axion production rate [6] was previously computed in [7] making use of the Hard Thermal Loop (HTL) [8,9] approximation (g 3 1), and considering only the axion coupling to gluons.The resulting space-time density of thermal axion production was: 3)g 4 3 T 6 64π 7 f 2 1 We show that scattering involving many particles can be neglected.This is trivially true in quantum field theory: for example the cross section for 2 → 3 scatterings is g 2 /(4π) 2 times smaller than the cross section of the dominant 2 → 2 scatterings.However, this is not generically true in thermal field theory, where the expansion parameter is g (rather than g 2 /(4π) 2 ) and where collinear kinematical configurations can enhance higher order scatterings by powers of 1/g, such that a resummation of 2 → n scattering becomes needed.This subtlety was noticed in the context of computations of photon emission from a quark-gluon plasma [10]: for example, the 2 → 3 process constructed adding to a qg → qg scattering a a q → qγ vertex, where q and γ are almost collinear (the directions of their moment differ by a small angle θ) allows a kinematical configuration where the propagator of the virtual gluon that mediates the scattering is enhanced by 1/θ 2 , while the gauge vertex q → qγ is only suppressed by θ.This results into a 1/θ enhancement of the scattering amplitude, cut-off by the thermal mass m ∼ gT , and thereby to a 1/g 2 enhancement of the 2 → 3 scattering rate.
We verified that no such collinear enhancement is present for axion production, because the axion vertices (such as the axion/gluon/gluon vertex aG µν Gµν ) are suppressed as θ 2 in the collinear limit.We also verified that the similar vertices relevant for graviton, gravitino [14], axino [15] production similarly lead to no collinear enhancement, being θ 2 suppressed.Thereby there is no need of adding higher order scatterings and previous computations remain valid.The functions F 1,2,3 (m/T ) are defined in eq. ( 5.1) and the thermal masses of the vectors within the SM are m/T = g 3 for gluons, m/T = 11/12g 2 for the W, Z and m/T = 11/12g Y for hypercharge.For comparison, the lower dashed curve is the result of [7] (eq.(1.1)) computed within the HTL approximation, valid in the limit g 3 1.
2 Outline of the computation

Effective axion Lagrangian
The effective action that describes axion couplings to SM particles at first order in the axion field a is written in the basis where the SM Lagrangian L SM does not contain the axion as [1] Here Gµν = 1 2 ε µναβ G αβ are the field strength duals; ε 0123 = 1; H is the Higgs doublet; the Weyl spinors ψ = {Q, U, D, L, E} are the SM fermions; f a is the effective axion decay constant in the convention where c 3 = 1; c 1 and c 2 are the axion couplings to electro-weak vectors; c H is the axion coupling to the Higgs; c ψ are the axion derivative couplings to the SM fermions.All c coefficients are real and dimensionless.In the full axion theory c H and c ψ are the PQ charges of the SM fields (they vanish in KSVZ axion models [3]), while c 1,2,3 also receive contributions from extra heavy fermions (not present in DFSZ axion models [4]). 2hile in previous computations only the axion/gluon coupling was considered, we want to consider all axion couplings.
For this purpose, it is convenient to perform a phase redefinition of the SM matter fields such that, at the first order in a, the c ψ and c H couplings are removed, at the price of shifting the axion coupling to vectors as follows where the sum runs over the 3 fermion generations.We used the fact that all SM matter field lie in fundamental representations with generators T a normalized as Tr(T a T b ) = δ ab /2.The transformation (2.2) also introduces an axion phase in the SM Yukawa couplings.For our purposes, all SM Yukawa couplings are negligibly small except the top Yukawa, for which the transformation induces the following axion phase: So, at first order in a, the transformation generates the following Lagrangian interaction: The thermal axion production rate will be computed in terms of c 3 (strong interactions), c 2 (weak interactions), c 1 (hypercharge) and c t (top Yukawa coupling).

Thermal production rate
According to the general formalism of thermal field theory [9], the thermal production rate of a weakly interacting scalar a is equivalently computed from the imaginary part of its propagator Π a as Here Π < a is the non time-ordered axion propagator and P = (E, p).Thermal field theory cutting rules allow to see that, at leading order in the SM couplings, eq.(2.6) is equivalent to the usual summing of all rates for the various tree-level processes that lead to axion production.
We illustrate the general discussion with the concrete example of the axion coupled to a simplified SM consisting only of gluons.In such a case the thermal axion production rate γ at leading order in g 3 can be obtained by computing the gg → ga scattering rate and thermally averaging it.
axion production rate at tree level, that would dominate with respect to one-loop contribution that we consider, encoded in the anomaly coefficients c 1,2,3 .The computation of such extra contribution would be analogous to the top Yukawa contribution discussed below.The leading-order gg → ga scattering rate in the thermal plasma is equivalently obtained by: a) summing the Feynman diagrams in the upper row, squaring the total amplitude, performing the thermal average; b) summing the imaginary parts of the two loop thermal diagrams in the lower row.In both cases the result is infra-red divergent, such that proper inclusion of higher order effects is needed.For simplicity, we here plotted the diagrams relative to a simplified world without quarks.
• The Feynman diagrams for gg → ga scatterings are plotted in the upper row of figure 2 and are named S (s-channel gluon exchange), T (t-channel), U (u-channel) and X (quartic vertex).When computing the rate in terms of scatterings, the rate is proportional to the modulus squared of the total amplitude, |S + T + U + X| 2 .
• The equivalent thermal diagrams at leading order in g 3 arise at two-loop level and are plotted in the lower row of figure 2, where they are named A, B, C, D. The rate is proportional to the their sum, that contains the various tree-level scatterings in the following way: We explicitly see that the thermal axion production rate γ However, both computations give an infra-red divergent result, because of the massless gluon in the T and U diagrams, or equivalently in the thermal diagram D. We employ the thermal field theory formalism because it is more appropriate for dealing with such issues.The infra-red divergence is regulated by the thermal gluon mass.We re-sum the thermal effects that modify the gluon dispersion relation by substituting the two-loop thermal diagram D with the one-loop 'Decay' diagram of fig.3, where the tree-level gluon propagator is replaced by the full thermal gluon propagator at leading order in the strong coupling.In a diagrammatic expansion, the 'Decay' diagram corresponds to diagram D, plus all higher order diagrams with any number of corrections to the gluon propagator, as illustrated in the right-handed side of figure 3. We give the name 'decay' to such resummed diagram because physically it describes the decay process g T → g T a of the thermal gluon g T , opened by the non-relativistic thermal corrections to the gluon propagator.The rationale for re-summing this class of higher-order effects is that they are enhanced by the 2 → 1 phase space factor, which is ∼ (4π) 2 bigger than the phase space relative to 2 → 2 scatterings.
In conclusion, the resummed total axion production rate is computed as (2.8) The computation of γ sub (subtracted scattering rates) is presented in section 3, and the computation of γ Decay is presented in section 4. Unlike in the HTL approximation, our technique does not need the introduction of an arbitrary splitting scale k * that satisfies the problematic conditions g 3 T k * T in order to control infra-red divergences.The total rate will be positive for any g 3 .
While we omitted quarks and other axion couplings to simplify the above discussion, of course we take them into account in the full computation.

Subtracted scattering rates
Table 1 lists the full scattering rates and the subtracted scattering contributions to the various axion production processes due to the axion/gluon/gluon interaction.It is important to notice that, unlike the total rate, the subtracted rates are infra-red convergent as expected: the infrared divergent factors 1/t and 1/u present in the full rate disappear from the subtracted rate.Actually, by performing computations in the Feynman gauge, we find that γ sub = 0.The same holds for the other SM vectors.
In order to double-check our result, we computed the subtracted scattering rates also as thermal diagrams that contribute to the non time-ordered axion propagator Π < a , see eq. (2.6).This computation is presented in the next part of this section.At the end of this section we also evaluate the top Yukawa contribution, which emerges from the axion interaction term in (2.5).

Diagram A
We first consider the contribution of the thermal diagram A in fig. 2. Making use of the Kobes-Semenoff rules (see e.g.[11]) we obtain the following contribution of the diagram A to Π < : where , P is the axion momentum and ∆ and ∆ < are respectively the tree level scalar propagator and non time-ordered propagator at finite temperature: where n B (x) ≡ (exp(|x|/T )−1) −1 .These emerge from the scalar part of the gluon propagators, while the contraction of Lorentz and color indices leads to the function α defined by where we have introduced C N ≡ N (N 2 − 1) with N = 3, 2, 1 for SU(3) c , SU(2) L and U(1) Y respectively.Looking at the thermal diagram we find that the kinematical configurations with non vanishing phase space in the integrand of eq.(3.9) are As anticipated in eq.(2.7), these contributions correspond to the interferences between the s-channel, t-channel and u-channel scattering diagrams.Although these contributions are separately non-trivial, we find that their sum is identically zero.So Π < A = 0. process Table 1: Axion production rate from the axion/gluon/gluon interaction.The total rate is gauge independent; the subtracted rate is gauge dependent and computed in the Feynman gauge.

Diagram B
We now turn to the thermal diagram B in fig. 2. We find: The function β, which results from the contraction of Lorentz and color indices, is The phase space of the integrand in (3.12) can be divided in three regions, which correspond to the interferences between the x-channel and the s, t and u-channels respectively: The sum of the three contributions in the phase space gives zero like for Diagram A: Π < B = 0.

Diagram C
Finally, we evaluate the thermal diagram C in fig.2, the contribution of the x-channel alone.It contributes to Π < an amount Like for thermal diagrams A and B we can divide the phase space in three parts, which this times correspond to the possibility to choose one out of three gluons and put it in the final state.The sum of these three contributions vanishes like for diagram A and B: Π < C = 0.

The top Yukawa contribution
The axion interaction in eq.(2.5) produces the following contribution to Π < a at the first nontrivial order in the perturbative expansion where ∆ < F is the tree level non time-ordered propagator at finite temperature for a fermion: and n F (x) ≡ (exp(|x|/T ) + 1) −1 .Like in the previous computation, the integral receives contributions from three distinct integration regions, which correspond to the effects that can be equivalently computed as Q 3 H * → a Ū3 , U 3 H * → a Q3 and Q 3 U 3 → aH scatterings (as well as their CP-conjugated processes): The first possibility, for example, leads to the following contribution to the production rate where the integral is performed on the intersection between the domains and These conditions emerge because k and q are the lengths of the three dimensional parts ( k and q) of the on shell momenta K and Q, once the delta functions in (3.16) are used, and z k and z q are the cosines of the angles between k and p and q and p respectively.The other two contributions lead to similar expressions.The total result due to the interaction in (2.5) is where the numerical factor is the Bose-Einstein and Fermi-Dirac correction with respect to the analytic result computed in Boltzmann approximation, n B,F (E) ≈ e −E/T .

Thermal vector decays
We start summarizing some well known results from quantum field theory at finite temperature that are relevant for our computations.

Thermal corrections to vector propagators
We list the full one-loop expressions for thermal corrections to a vector [12,13,9] with fourmomentum K = (ω, k) (K 2 = ω 2 −k 2 ) with respect to the rest frame of the thermal plasma.We denote by U µ the four-velocity U µ of the plasma.We use the Feynman gauge where all effects are condensed in two form factors even in the non-abelian case [12].Polarizations are conveniently decomposed in T ransverse (i.e.orthogonal to K and to k), Longitudinal (i.e.orthogonal to K and parallel to k) and parallel to K. The corresponding projectors (Π where the corrections are contained in the scalar functions π 0 (k, ω) (quantum corrections at T = 0) and π T (k, ω) and π L (k, ω) (thermal corrections), explicitly given in [14] for a general theory.The corresponding non-time ordered propagator is Here, ρ T , ρ L are the spectral densities for the transverse vectors and longitudinal vectors respectively Furthermore, k 0 > 0 describes a vector in the final state, and k 0 < 0 describes a vector in the initial state: this convention allows to compactly describe all possible processes.Indeed the factors gives the usual statistical factors: −n (number of particles in the initial state) or 1 ± n (stimulated emission or Pauli-blocking in the final state), where n B (E) ≡ 1/(e |E|/T − 1) is the Bose-Einstein distribution.
Fig. 4 shows the spectral densities in the HTL limit (g 1, left panels) and in a realistic situation (g ∼ 1, right panels).In the realistic case, the poles get smeared acquiring a finite width.A more significant difference arises below the light cone (ω < k), where both the HTL and the one-loop spectral densities do not vanish.This describes 'Landau damping' i.e. the fact that particles exchange energy with the thermal plasma.However, the HTL approximation cannot be applied at k ∼ T (a region relevant for us, since g ∼ 1), and indeed it misses that at k T spectral densities get suppressed by an exponential Boltzmann factor.
Figure 4: One loop thermal densities ρ T (ω, k) (upper row) and ρ L (ω, k) (lower row) for a vector in the Hard Thermal Loop limit g 1 (left) and for g ∼ 1 (right).In the HTL limit there is a pole above the light cone and a continuum below the light cone.In the full case the pole above the light cone acquires a finite width and becomes a continuum, and the continuum below the light cone gets Boltzman suppressed at k m.

Axion production via vector thermal decays
In section 3 we have computed the subtracted axion production rate; we here compute the resummed 'Decay' diagram of fig.3, which reduces at leading order to diagram D in fig. 2. As already stated, the rationale for re-summing only this class of higher-order effects is the phase space enhancement of the 2 → 1 processes (relative to the 2 → 2 scatterings).Therefore, the residual gauge dependence in our result is expected to be of relative order g 2 /π 2 .The computation applies to all SM vectors V = {g, W, Y } with gauge couplings α i = {α 3 , α 2 , α Y } and dimension of the gauge group d i = {8, 3, 1}.
The resummed contribution to the axion propagator Π < a is The vector quadri-momenta are K µ and Q µ = P µ −K µ , where P µ is the axion quadri-momentum.
Inserting the parametrization P = (p 0 , p, 0, 0) , K = (k 0 , k cos θ k .ksin θ k , 0) , Q = (q 0 , q cos θ q , q sin θ q , 0), ( Gauge group Table 2: Numerical coefficients for vector thermal mass m 2 i = 1 6 g 2 i T 2 (N + N S + N F /2) in the SM in terms of the SU(N ) factor, of the number of fermions N F and of scalars N S .
we obtain where we used the decomposition of the resummed propagator given in (4.3).In order to compute the integral, it is convenient to multiply by 1 = d 4 Q δ (K + Q − P ).After performing the angular integrations over θ k and θ q , and using the equations, we obtain Note that the integration range is restricted to |p−k| ≤ q ≤ p+k.Finally, for each factor of the SM gauge group we computed this integral numerically using the spectral densities described in the previous section.We followed the method provided in [14].

Result
The total axion production rate due to all axion couplings c 3 , c 2 , c 1 , c t and taking into account all large SM couplings, g 3 , g 2 , g Y , y t is where the thermal masses of SM gauge bosons of SU(3) c , SU(2) L and U(1) Y are (see table 2): Fig. 1 shows our results for the functions F 1,2,3 that parameterize the axion production rate due to gauge interactions, while our result for the top Yukawa part is given analytically.Fig. 5 shows the four contributions to the axion production rate as function of the temperature.As long as c t ∼ c 3 the top Yukawa axion production rate gives the dominant contribution because it arises at tree level, while the anomalous axion couplings arise at loop level.Previous works ignored the top Yukawa effect and computed only the function F 3 within the HTL approximation i.e. in the limit of small strong gauge coupling, g 3 1.We see that in this limit our improved computation reproduces to the HTL limit.However, when g 3 is set to its physical value, g 3 ≈ 1, the results differ: the HTL approximation breaks down and the HTL rate function F HTL 3 becomes unphysically negative for a large enough g 3 , while our result grows with increasing g 3 .

Cosmological axion yield
In the usual scenario of reheating after inflation, the inflaton φ with energy density ρ φ decays with width Γ φ into SM particles (excluding the axion).The reheating temperature T RH is defined as the temperature at which Γ φ equals H R , the expansion rate due to the radiation density only: Here we neglect possible non-equilibrium effects at T T RH [17].T RH effectively is the maximal temperature of the universe.Indeed, while higher temperatures exist, particles produced at T > T RH are diluted by the entropy released by inflaton decays, as described by the Z − 1 = −Γ φ ρ φ /4Hρ R term in the Boltzmann equations (5.4) Here H = 8πρ/3/M Pl is the Hubble parameter, z = T RH /T , s = 2T ( The latter factor in eq.(5.6) is the order-one term among square brackets in eq.(5.1).The approximated analytical expression of eq.(5.5) valid for intermediate values of r is obtained by fitting the numerical solution.
Even for the lowest possible value of f a > ∼ 5 × 10 9 GeV and taking into account the new top effect in γ a , eq. (5.6) implies that axions decoupled at T > ∼ M Z , when the number of relativistic SM degrees of freedom g SM still included all SM particles.Thereby, when SM particles later become non-relativistic, they annihilated heating photons and neutrinos, but not axions.This means that today, and at the epoch of CMB decoupling, thermal axions constitute a small fraction of the total relativistic energy fraction, conveniently parameterised by the usual "effective number of neutrinos" as 3  ∆N eff ν = 0.0264 Y a Y eq a . (5.7) The phenomenological manifestations in cosmology of a thermal axion component of the universe are analogous to having an extra freely-streaming neutrino component.Such effects can be parameterised by the axion contribution to effective number of neutrinos ∆N eff ν and by the axion mass m a ≈ 0.6 meV(10 10 GeV/f a ). 3 Today photons have temperature T γ and neutrinos have temperature T ν = T 0 (4/11) 1/3 , for a total of g * s = 43/11 effective entropy degrees of freedom.Axions went out of thermal equilibrium at T M Z would have a present temperature T a = T γ (g * s /g SM ) 1/3 = 0.903 K, which corresponds to ∆N eff ν = 4(T a /T ν ) 4 /7 ≈ 0.0264.We recall that the SM alone predicts N eff ν ≈ 3.046 where the small deviation from 3 is due to imperfect neutrino decoupling when electrons become non-relativistic [16].
Full cosmological bounds in the plane (∆N eff ν , m) were computed in fig.6a of [18] and, more recently, in fig.28 of [19].In practice, present global fits of cosmological data find ∆N eff ν = 0.48 ± 0.48 [19]: the uncertainty is more than one order of magnitude above the maximal thermal axion effect.Future experiments which are being discussed, such as CMBpol, can reduce the uncertainty on N eft ν to ±0.044 [20].

Conclusions
In conclusion, we improved over previous computations of the thermal axion density in two ways: 1.By including higher-order effect enhanced by the thermal decay kinematics gluon → gluon + axion.Unlike the leading order result, which becomes negative for g 3 > 1.5, our result behaves physically for all relevant values of the strong gauge coupling.
2. By considering all axion couplings to SM particles; not only to gluons.The strong interaction contribution receives new contributions from the axion couplings to quarks, as encoded in the difference between c 3 and c 3 , eq. (2.3).Electroweak effects are small.More importantly, as long as the axion couples to the top quark, there is a new effect related to the top Yukawa coupling, which dominates by 3 order of magnitude over the effect related to the strong gauge coupling (see fig. 5).
Our result for the thermal axion production rate is given in eq.(5.1) in terms of the axion couplings c 3 (strong interactions), c 2 (weak interactions), c 1 (hypercharge) and c t (top Yukawa coupling) defined in eq.(2.3) and (2.5).The thermal functions F 3 , F 2 , F 1 are numerically plotted in fig. 1.Furthermore, we have shown that there are no collinear enhancements (in analogous computations such effects are present and require resummation of extra classes of thermal diagrams).
The thermal axion abundance is then computed adopting the usual simplified model of reheating (rather than the instantaneous reheating approximation adopted in previous works) and allowing for the possibility of a thermalised axion.We provide in eq.(5.5) a simple numerical approximation for the final axion abundance.

Figure 1 :
Figure 1: Our result for the axion production rate as function of the thermal mass of vectors.The functions F 1,2,3 (m/T ) are defined in eq.(5.1) and the thermal masses of the vectors within the SM are m/T = g 3 for gluons, m/T = 11/12g 2 for the W, Z and m/T = 11/12g Y for hypercharge.For comparison, the lower dashed curve is the result of[7] (eq.(1.1)) computed within the HTL approximation, valid in the limit g 3 1.

Figure 2 :
Figure2: The leading-order gg → ga scattering rate in the thermal plasma is equivalently obtained by: a) summing the Feynman diagrams in the upper row, squaring the total amplitude, performing the thermal average; b) summing the imaginary parts of the two loop thermal diagrams in the lower row.In both cases the result is infra-red divergent, such that proper inclusion of higher order effects is needed.For simplicity, we here plotted the diagrams relative to a simplified world without quarks.

Figure 3 :
Figure3: The thermal diagram 'Decay', where the gluon propagator includes one-loop thermal corrections, is equivalent to the thermal diagram D (thick lines denote the propagator of the thermal gluon g T ) plus the resummation of higher order diagrams with corrections to the gluon propagator.

Figure 5 :
Figure 5: The four contributions to the thermal axion production rate γ a induced by the SM couplings y t (upper black curve), g 3 (red curve), g 2 (blue), g Y (green) for unity values of the axion couplings c t = c 3 = c 2 = c 1 = 1 in eq.(5.1).The red dashed line is the previous result for the strong coupling contribution computed in Hard Thermal Loop approximation.
3g SM π 2 /45 is the entropy density of SM particles (g SM = 427/4), n a is the axion number density, Y a = n a /s, and Y eq a = n eq a /s ≈ 0.00258 with n eq a = ζ(3)T 3 /π 2 is the thermal equilibrium value of Y a .The solution to the Boltzmann equations for the axion abundance at T T RH is