Thermal production of axino Dark Matter

We reconsider thermal production of axinos in the early universe, adding: a) missed terms in the axino interaction; b) production via gluon decays kinematically allowed by thermal masses; c) a precise modeling of reheating. We find an axino abunance a few times larger than previous computations.


Introduction
The strong CP problem can be solved by the Peccei-Quinn symmetry [1], that manifests at low energy as a light axion a with a decay constant f > ∼ 5 10 9 GeV [2]. In supersymmetric models the axion a gets extended into an axion supermultiplet which also contains the scalar saxion s and the fermionic axinoã [3]. Depending on the model of supersymmetry breaking, the axino can easily be lighter than all other sparticles, becoming the stable lightest supersymmetric particle (LSP) and consequently a Dark Matter candidate [4,5,6]. It is thereby interesting to compute its cosmological abundance. The axino can be produced: i) from decays of the next to lightest sparticle (NLSP), such that Ωã = mãΩ NLSP /m NLSP [5]; plus ii) thermally in the early universe when the temperature T was just below the reheating temperature T RH . We here reconsider the thermal axino abundance, improving on previous computations [7] in the following ways: a) in section 2 we show that, beyond the well known axino/gluino/gluon interaction, there is a new axino/gluino/squark/squark interaction, unavoidably demanded by superymmetry, that contributes at the same order to the usual 2 → 2 scatterings that produce axinos; b) in section 3 we show that axinos are also thermally produced by 1 → 2 decays kinematically allowed by the gluon thermal mass. The gluon → gluino + axino process gives the dominant contribution in view of the large value g 3 ∼ 1 of the strong coupling constant; c) in section 4 we precisely model the reheating process.
As a result the axino production rate is significantly enhanced.

Axino couplings
The effective coupling between the axion supermultiplet A = (s + ia)/ √ 2 + √ 2θã + · · · and the strong gauge vectors, described by the gluon G a and gluinog a in the super-multiplet where θ is the super-space coordinate, f is the axion decay constant, related in a well known way to the parameters of each axion model at hand. The RGE evolution of α 3 = g 2 3 /4π encodes the RGE renormalization of the full operator.
By expanding in components and converting to Dirac 4-component notation we get: The first term is the usual axion coupling to the gluon; the second term is the corresponding saxion coupling and both are accompanied by couplings to gluinos and to squarks. The terms in the lower row are the axino couplings; the latter term was not considered in the literature [4,7,8]. It contributes to the axino production rate at the same order as the well known third term, as can be seen by inserting the explicit value of the strong D-term, where the sum runs over all squarksq.
In addition to eq. (1), the axion supermultiplet can have extra non-minimal couplings to the electroweak vectors [4,8]. We only consider the presumably dominant strong interaction contribution to the thermal axino production rate.

Thermal axino production rate
The axino production thermal rate has been computed in [7] at leading order in the strong gauge coupling g 3 . This roughly amounts to compute the 2 → 2 scatterings listed in table 1, with thermal effects ignored everywhere but in the propagator of virtual intermediate gluons: indeed a massless gluon exchanged in the t-channel gives an infinite cross-section because it mediates a long-range Coulomb-like force; the resulting infra-red logarithmic divergence is cut off by the thermal mass of the gluon, m ∼ g 3 T , leaving a ln T /m. The Hard Thermal Loop (HTL [9,10]) approximation (m T i.e. g 3 1) gives the following result for the space-time density of scatterings into axinos [7]: This production rate unphysically decreases for g 3 > ∼ 0.7 becoming negative for g 3 > ∼ 1.2 [7]. Fig. 1 illustrates that the physical value, g 3 ≈ 0.85 at T ∼ 10 10 GeV, lies in the region where the leading-order rate function F HTL (g 3 ) (dashed line) is unreliable. Fig. 1 also illustrates our final result: F HTL is replaced by the function F (continuous lines); they agree at g 3 1 and differ at g 3 ∼ 1.
process Table 1: Squared matrix elements for axino production in units of g 6 3 /128π 4 f 2 summed over all polarizations and gauge indices. g,g, q,q denote gluons, gluinos, quarks, squarks. The gauge factors equal C = |f abc | 2 = 24 and C = q |T a ij | 2 = 48 after summing over all quarks.
To improve the computation going beyond the leading-order HTL approximation, we notice that (analogously to what discussed in [11] for gravitino thermal production) the new decay process gluon → gluino + axino (4) first contributes to the axino production rate γ at next order in g 3 . Indeed this process is made kinematically allowed by the gluon thermal mass. Despite being higher order in g 3 , the decay rate is enhanced by a phase space factor π 2 , because a 1 ↔ 2 process has a phase space larger by this amount than 2 ↔ 2 scatterings. Thereby such decay gives a correction of relative order (πg 3 ) 2 to the axino production rate γ. Subsequent higher order corrections should be suppressed by the usual g 3 /π factors. Our goal is including the enhanced higher order terms, and this finite-temperature computation is practically feasible because a decay is a simple enough process.
Proceeding along the lines of [11], the axino production rate is precisely defined in terms of the imaginary part of the thermal axino propagator at one loop [10], that we compute using the resummed finite-temperature propagators for gluons and gluinos in the loop. The resulting expression can be interpreted as the thermal average of the decay process (4), taking into account that the gluon and gluino thermal masses break Lorentz invariance, and that actually a continuum of 'masses' is present, as precisely described by the gluon and gluino thermal spectral densities [10].
Thermal field theory is just a tool to describe the collective effect of scatterings, and the decay diagram indeed resums an infinite subset of scatterings, and in particular some lowestorder 2 → 2 scatterings. Thereby, in order to avoid overcounting, 2 → 2 scatterings must be added subtracting the contributions already included in resummed decay [11], which includes the modulus squared of single Feynman 2 → 2 diagrams. What remains are interferences between different Feynman 2 → 2 diagrams [11]. Table 1 gives explicit values for the total and subtracted axino scattering rates. Unlike the total rate, the subtracted rate is infra-red convergent: no 1/t factors appear because all divergent Coloumb-like scatterings are included in the decay diagram. Subtracted rates for processes C, D, E, G, I vanish: a single diagram contributes, such that no interference terms exist. This is not the case for scatterings H and J,  [7], which asymptotically agrees with our result in the limit of small g 3 , and behaves unphysically for g 3 ∼ 1.
where the second axino interaction in eq. (5) contributes, changing the coefficient of the first term among parenthesis with respect to the corresponding table 1 of [7], as well as giving a non-vanishing subtracted contribution.
Actually there is a stronger connection between the axino and the gravitino thermal production rate. The latter contains the contribution of the Goldstino χ production rate, induced by its coupling to the divergence of the MSSM super-current S µ : having including only the contribution from the gluino mass M 3 . By comparing with eq. (1) we notice that the Goldstino and the axino couple to the same combination of operators. Even in the Goldstino case, the second operator in ∂ µ S µ was initially missed and finally noticed by [13]. While the goldstino mass is predicted in terms of the supersymmetry-breaking F term as m 3/2 = F/ √ 3M Pl , the axino mass mã is a model-dependent free parameter. In conclusion, after converting M 3 /F ↔ √ 2α 3 /4πf the result of [11] for the Goldstino thermal production rate is translated into the axino production rate: Our result for F , plotted in fig. 1, behaves physically for physical values of g 3 ∼ 1.

Thermal axino abundance
The axino abundance can now be computed by integrating the relevant cosmological Boltzmann equation for the axino number density nã [14,7,11,12] dnã dt + 3Hnã = γ as well as the equations that describe reheating. The reheating temperature T RH is approximately defined as the maximal temperature of the universe ('instantaneous reheating'), or more precisely as the temperature at which the inflaton decay rate becomes smaller than the Hubble rate H R ≡ 8πρ R /3/M Pl (computed including only the energy density ρ R of the thermal bath), such that whatever happened earlier at higher temperatures got diluted by energy injection from inflaton decay [12]. The two different definitions give the following solutions [11] for the axino number abundance normalized to the entropy density s: Using the numerical factor appropriate for realistic reheating, the present axino mass density is mã GeV T RH 10 4 GeV with g 3 renormalized around T RH and F (g 3 ) plotted in fig. 1 (the analytic approximation is appropriate around g 3 ≈ 1). The axino energy density must be equal or smaller to the DM density Ω DM h 2 = 0.110 ± 0.006 [15], and the resulting phenomenology was studied in [16]. In fig. 2 we plot the values of f such that thermally produced axinos have an abundance equal to the observed DM abundance. We recall that f < ∼ 10 12 GeV in order to avoid a too large axion DM abundance, unless the initial axion vev is close to the minimum of the axion potential [17]. In eq. (8) and (9) we assumed that nã n eq a ≈ 1.8 10 −3 s, otherwise the axino reaches thermal equilibrium, giving Ωã Ω DM for mã 10 eV.

Conclusions
In section 2 we derived the axino coupling to the strong sector, finding that it is described by two terms: one is the well knownãgG coupling, the otherãgq * q term was missed in previous studies, although they contribute at the same order to the thermal axino production rate. This makes the axino interaction fully analogous to the Goldstino interaction, such that we could infer the axino thermal production rate from the Goldstino rate, computed in [11]. Such result goes beyond the leading order computations [7] based on the Hard Thermal Loop (HTL) approximation g 3 1, which gives unphysical results at the physical value of g 3 ∼ 1. As a consequence we find an enhancement, plotted in fig. 1, by a factor of 6 (3) at T RH = 10 4 (10 7 ) GeV. The function F determines the axino rate as described by eq.s (6) and (9).