Axino Cold Dark Matter Revisited

Axino arises in supersymmetric versions of axion models and is a natural candidate for cold or warm dark matter. Here we revisit axino dark matter produced thermally and non-thermally in light of recent developments. First we discuss the definition of axino relative to low energy axion one for several KSVZ and DFSZ models of the axion. Then we review and refine the computation of the dominant QCD production in order to avoid unphysical cross-sections and, depending on the model, to include production via SU(2) and U(1) interactions and Yukawa couplings.


Introduction
The lightest supersymmetric particle (LSP) with R-parity conservation is absolutely stable and can contribute to the present Universe energy density as a dominant component of cold dark matter (CDM). If the strong CP problem is solved by introducing a light pseudoscalar, an axion, its fermionic SUSY partner, an axino, can be a natural candidate for CDM if it is the LSP. The relic axino CDM can be produced either non-thermally, through out-of-equilibrium decays, or thermally, through scatterings and decays in the hot plasma, as originally shown in Refs. [1,2]. The scenario was subsequently extensively studied during the last decade [3][4][5][6][7][8][9][10], with either a neutralino or a stau as the next-tolightest supersymmetric particle (NLSP). Ways of testing the axino CDM scenario at the LHC were explored in Refs. [8,[11][12][13] and implications for Affleck-Dine Baryogenesis in Ref. [14].
The strong CP problem is naturally solved by introducing a very light axion field a. The axion appears when the Peccei-Quinn (PQ) symmetry is broken at some scale f a . Below the PQ scale, after integrating out heavy quarks carrying PQ charges [15], an effective axion-gluon interaction is given by where α s is the strong coupling constant, G is the field strength of the gluon field and G µν ≡ 1 2 ǫ µνρσ G µν is its dual with ǫ 0123 = −1. Different types of axion models have been proposed, with distinctively different couplings to Standard Model (SM) fields, depending on their PQ charge assignment.
Very light axion models contain a complex SM singlet scalar field carrying a PQ charge. In the Kim-Shifman-Vainstein-Zakharov (KSVZ) class of models [15] the PQ charges are assigned to new heavy quarks, while in the Dine-Fischler-Srednicki-Zhitnitskii (DFSZ) approach [16] the PQ charges are assigned to the SM quarks. This difference is the origin of different phenomenological properties [17] of the KSVZ and DFSZ classes of models since in the low energy effective theory at the electroweak (EW) scale (after integrating out heavy fields) the gluon anomaly term is the source of all interactions in the KSVZ models while the Yukawa couplings are the source of all interactions in the DFSZ models. For solving the strong CP problem, one needs a coupling of the axion to the gluon anomaly and this is generated by a heavy quark loop in the KVSZ models, appearing as a non-renormalizable effective coupling when those heavy fields are integrated out. In the DFSZ models instead the coupling is generated by SM quark loops.
The couplings arising in these two popular classes of models will be discussed in the next section. For the KSVZ axion, constraints on the PQ scale f a have been obtained from astrophysical and cosmological considerations and the scale is limited to a rather narrow window 10 10 GeV ∼ < f a ∼ < 10 12 GeV [17], while for the DFSZ case precise constraints on f a have not been derived yet.
The axino as a candidate for CDM was originally studied mostly for the SUSY version of the KSVZ axion model, in an important production mode corresponding to the interaction term given by Eq. (1.1) [2]. The supersymmetrization of axion models introduces a full axion supermultiplet A which contains the pseudoscalar axion a, its scalar partner saxion s, and their fermionic partner axinõ a, where F A stands for an auxiliary field and ϑ for a Grassmann coordinate. The effective axion interaction of Eq. (1.1) can easily be supersymmetrized using the superpotential of the axion and the vector multiplet W α containing the gluon field Effective axion multiplet interactions with the other gauge bosons can be obtained in a similar way. Axino production from QCD scatterings due to interaction Eq. (1.3) has been considered in the literature in different approximations, which have led to somewhat different numerical results. The main technical difficulty and source of uncertainty is the question of how to regulate the infrared divergences due to the exchange of massless gluons. In the original study [2] a simple insertion of a thermal gluon mass to regulate infrared (IR) divergences was used and the leading logarithmic term was obtained, without much control over the subleading finite piece, since the thermal mass was introduced by hand. Subsequently, a hard thermal loop (HTL) resummation method was applied in Ref. [5], allowing for a self-consistent determination of the gluon thermal mass and more control over the constant term in the high energy region of axino production. Recently, in Ref. [6] a new calculation was presented which, although not gauge invariant, includes more terms of the perturbative series, in particular the decay of the gluon whose thermal mass can be larger than the gluino and axino mass taken together. The two latter methods have their own advantages and limitations, and they coincide in the high energy region where the convergence of the perturbation series is stable. In Ref. [6] a previously neglected dimension-5 term in the Lagrangian, containing purely interactions between supersymmetric particles was also included. This however changes the axino production rate by less than 1%. In our paper, we adopt the original way of effective mass approximation but we improve it to make the result positive definite. Although this method is not gauge invariant either, it is only known viable method at relatively low reheating temperatures, which is the regime important for the axino as cold dark matter candidate. We shall come back to this discussion in more detail below.
In the calculations published so far the couplings of the axino to the gauge multiplets other than the one of the gluon was also neglected. In fact, in Ref. [2] a chiral transformation of the left-handed lepton doublets was performed to remove the axion SU (2) anomaly interaction and to leave only the axion U (1) contribution. Then the axion SU (2) anomaly coupling re-appears in principle from the leptonic loops, which are independent of the fermion masses. The corresponding axino loops, on the other hand, are suppressed by the ratio m 2 lepton /M 2 where M is the largest mass in the loop. Therefore the error in neglecting the axion SU (2) anomaly is estimated to be at most of order m 2 τ /m 2 τ compared to the SU (3) term. However, in supersymmetric extensions of axionic models, the chiral rotation on the lepton fields involves also their scalar counterparts, the sleptons and the saxion, which generates additional couplings for those fields. It is also not clear if it is justified to perform a redefinition in a fully supersymmetric way when supersymmetry is in any case broken. Rather than going into the swamp of such interactions with unknown slepton and saxion masses as parameters, we will keep here the general axion SU (2) anomaly interaction and its SUSY counterpart, which is more tractable, and examine its effect on the axino abundance.
Moreover, in the present work we will also consider the role of Yukawa-type interactions between the axino and the matter multiplets, which arise either at one loop or at tree level, depending on the model. We will investigate what effect these terms may have and in particular how model dependent our results for the axino abundance are. In the previous study, only the KSVZ model was considered for the axino production. However in the DFSZ model axino production is dominated by the Yukawa coupling and the dependence on the reheating temperature is quite different from that in the KSVZ model. Only recently Refs. [33] and [47] considered axinos in the DFSZ model. In particular, Ref. [33] gave a simple approximate formulae for the relic density of light axino as dark matter. In our paper, we study axino production in the DFSZ model and the suitability of the DFSZ axino as dark matter in a more complete way.
In view of the recent developments in estimating the QCD contribution, in this paper we also update and re-examine relic CDM axino production, with an emphasis on an estimate of the uncertainties as well as on model dependence. In particular, our updated analysis of the axino CDM scenario improves on the previous works in the following aspects: (i) an inclusion of some previously neglected terms in the axino production and of subleading terms in the mass of the gluon beyond the logarithmic and constant pieces -these last parts ensure that the cross-section remains positive even for the invariant mass s smaller than the gluon thermal mass; (ii) an explicit calculation of axino production via SU (2) and U (1) interactions; (iii) a derivation of the axino abundance in specific implementations of the KSVZ and DFSZ models. (iv) an update on the constraints on the reheating temperature T R 1 for the both the neutralino and the stau as the NLSP using the current WMAP-7 result on the DM relic density and relevant structure formation data.
In Sect. 2, we define the axino by its relation to the shift symmetry of the KSVZ and DFSZ axions. In Sect. 3 we calculate axino production rate for the KSVZ and DFSZ axion models, and in Sect. 4 1 We assume the instant reheating approximation and define the reheating temperature as the maximum temperature at which standard Big Bang expansion with a thermalised bath of SM particles starts. We can easily translate the axino abundance given with this conventional definition of the reheating temperature to more specific reheating scenarios. E.g., Ref. [6] considers a reheating process from the decay of the inflaton and obtains a slightly smaller abundance of axinos, reduced by a factor of 0.745. In comparing our results we account for this difference.
we discuss several constraints on the scenario arising from cosmology. Finally, in Sect. 5 we present our conclusions.

The framework
In a recent review [17] low energy axion interactions were given in terms of the effective couplings with the SM fields c 1 , c 2 , and c 3 , which arise after integrating out all heavy PQ charged fields. The effective axion Lagrangian terms are where c 3 can be defined to be 1 (if it is non-zero) by rescaling f a . 2 The c 1 term is the PQ symmetry preserving derivative interaction of the axion field that can be reabsorbed into the c 2 term by a partial integration over on-shell quarks. For the c 2 terms, we have only kept the lowest order terms (proportional to 1/f a ), while in principle an infinite series of terms in a/f a arises. In the following we will consider the two popular scenarios mentioned earlier: the KSVZ and the DFSZ classes of axion models. The KSVZ class of axion models corresponds to the choice c 1 = c 2 = 0, c 3 = 1, and the DFSZ one to c 1 = c 3 = 0 and c 2 = 0, after integrating out the heavy field sector responsible for PQ symmetry breaking. In the latter model, if the Higgs doublets H u,d carry respective PQ charges Q u,d , the SM fields also carry PQ charges, 3 see Table I, and the anomaly interaction proportional to c 3 = 0 arises from SM quark loops. The axion mass is given by the strong interaction and is proportional to |c 2 + c 3 |. Hence the sum c 2 + c 3 , if it is non-zero, defines the QCD axion. Only this combination of the two couplings is physical, since a chiral axion-dependent PQ rotation of the quark fields can shift the values of c 2 and c 3 , while keeping c 2 + c 3 constant. This is connected to the well-known fact that, if one of the quark masses is zero then the anomaly becomes unphysical and can be reabsorbed in the rotation of the massless field.
In the supersymmetric version of axionic models, the interactions of the saxion and the axino with matter are related by supersymmetry to those of the axion. Hence the definition of the axion at low energy must be connected to the definition of the axion multiplet, and therefore of the saxion and the axino, at energies above the EW scale. 4 Below we will examine more closely the KSVZ and DFSZ models of the axion. In both models one imposes the PQ symmetry at a high energy scale and the axion emerges from its spontaneous symmetry breaking. As a specific example, let us consider the PQ sector at high energy with the PQ 2 The strong coupling g 3 is usually omitted by absorbing it to the field strengths. Here, fa denotes the so-called axion decay constant which is typically related to the vacuum expectation value V of the PQ symmetry breaking scalar field by fa = V /N where N is the domain wall number. 3 In variant axion models, Q u,d may have family dependence but here we suppress family indices. They can be inserted if needed. 4 The supersymmetrization of axion models was first discussed in Refs. [18][19][20]]. An explicit model was first constructed in [21], and the first cosmological study was performed in Ref. [22]. breaking implemented as in Ref. [21] W where Z, S 1 , and S 2 are gauge singlet chiral superfields, f Z is a Yukawa coupling and V a is a parameter in the Lagrangian which determines non-zero VEVs S 1 = S 2 = V a in the minimization process. The superfields transform under the U (1) PQ symmetry as The potential in the global supersymmetric limit is where g is the gauge coupling of the gauge groups and φ denotes collectively all the scalar fields. With the superpotential given by Eq. (2.2) one has At the tree-level both S 1 and S 2 develop VEVs and break the PQ symmetry. The fermionic partners of Z and S ′ = (S 1 + S 2 )/ √ 2 combine to become a Dirac fermion with mass m Z = √ 2f Z V a , while S = (S 1 − S 2 )/ √ 2 contains the axion field and can be identified with the axion multiplet. From Eq. (2.6) we note that, if Z = 0 and SUSY is not broken then both the axino and the saxion are mass degenerate with the axion.
However, with soft SUSY breaking terms included, Z can develop a non-zero value, in which case both the saxion and the axino become massive, independently of the axion. Therefore, a full specification of the SUSY breaking mechanism is needed in order to determine the mass of the axino exactly. This was first studied in Ref. [21] for the superpotential above, while another superpotential and SUSY breaking with a very small axino mass was given in Ref. [23]. Recently the case of a direct coupling between the axion and SUSY breaking sector was discussed in Ref. [24].

The KSVZ model
In the KSVZ approach in order to obtain the anomalous interaction of the axion and the gluon fields, one introduces the heavy quark fields Q L and Q R in the superpotential as Ref. [15] (2.7) After the U (1) PQ symmetry breaking takes place, as discussed above, Q L and Q R combine to become a heavy Dirac fermion with mass m Q = f Q V a , Assuming that V a ≫ m 3/2 , the scalar partner has practically the same mass and the whole supermultiplet can be integrated out in a supersymmetric way. The low-energy Lagrangian can then be obtained by integrating out the heavy quark multiplet or by using anomaly matching condition. In this way one finds that the axion anomaly term becomes where the axion a is the pseudoscalar component of the superfield S ≡ (S 1 − S 2 )/ √ 2, and the axion decay constant is f a = 2V a . N Q is the number of the heavy quarks and we consider N Q = 1 in our case. Then we have c 1 = c 2 = 0 and c 3 = 1 [17].
The low-energy interactions of the saxion and the axino fields can be obtained in the same way by integrating out the heavy (s)quark fields. However, in the limit of unbroken SUSY the low-energy effective Lagrangian should be given in a SUSY invariant form, including Eq. (2.8), as The axion superfield A and W a α are given by [25] (2.10) The effective Lagrangian in terms of the Bjorken-Drell gamma matrices then reads 5 .
., (2.11) where the gluino and axino 4-spinors are given bỹ (2.12) The effective Lagrangian (2.11), including the axino interaction with the D-term has recently been used in Ref. [6] to calculate the thermal production rate of axino dark matter and will be the basis also for the present analysis.

The DFSZ model
In the DFSZ framework, the SU (2) L × U (1) Y Higgs doublets carry PQ charges and thus the light quarks are also charged under U (1) PQ . The charge assignment is shown in Table 1. The anomaly coupling aG G can be obtained after electroweak symmetry breaking (EWSB) through the coupling of the axion to the Higgs doublets which couple to the light quarks. To this end, one adds a nonrenormalizable term to the PQ breaking superpotential [26] where W PQ is given in Eq. (2.2) and Note that here f s ∼ µM P /V 2 a ∼ 1 generates a phenomenologically acceptable supersymmetric µ-term.
The superpotential including light quarks is given by (2.14) Before the EW symmetry is broken but after the U (1) PQ symmetry is broken, the massless axion Goldstone boson is identified as Here we defined the component fields in the same way as in Eq. (2.10), and similarly for S 2 . Instead, the axino mass eigenstate can be obtained from the mass matrix in the basis of (ã 1 ,ã 2 ,Z) and we have used z 0 = Z . The lightest state, which we will identify with the axino, has mass f Z z 0 and is given byã which coincides with Eq. (2.15). Since EW symmetry is not broken, the axion (and thus the axino) do not mix with the Higgs (and higgsino) and therefore the axion cannot have the SU (3) c anomalous interaction with gluons generated at one loop via the quark triangle diagrams. 6 However, SU (2) L and U (1) Y anomalous interaction can be generated via a higgsino triangle loop through the Yukawa coupling derived from Eq. (2.13). They are given by These anomaly interactions appear also for the axino, but may obtain corrections of order of one from SUSY breaking, compared to the axion couplings, and such uncertainties will be later encoded in the coefficients C aW W , C aY Y .
After the EW and the U (1) PQ symmetries are both broken, the axion is identified with the Goldstone boson of the broken U (1) PQ symmetry, given by where a s = (a 1 − a 2 )/ √ 2, and we expanded the Higgs field as Equating (2.20) with the last term of Eq. (2.22), we obtain up to a common normalization constant. The interactions of the axion with the matter fields are obtained through the axion part of the phase of the Higgs fields, (2.24) Considering the Yukawa interaction from the superpotential, Eq. (2.14), the Lagrangian terms for the up-type quark axion couplings are given by 26) and similarly for the down-type quarks. We can now compare the above to the coefficients in the definition of the axion [17] and obtain for following values for c u,d 2 ; compare Eq. (2.1), After integrating out all the Higgs fields except the axion supermultiplet, which remains light, all the axion couplings arise from the c 2 terms given in Eq. (2.26) and Eq. (2.27) at low energies above the quark masses. At one loop in the SM fermions one obtains the axion-anomaly interaction term. It is then given by 28) where N = 6 and again f a = 2V a /N . In the supersymmetric limit, below EW symmetry breaking, the axino mass eigenstate can be read off from Eq. (2.20) and is given bỹ Hereh d,u denote the fermionic components of the electrically neutral parts of H d,u . However since supersymmetry is broken, in general in the DFSZ models the axino mixes with the higgsinos differently than the axion with the Higgs and the mass eigenstate is not exactly the state given in Eq. (2.29). Such a field is a good approximation to the physical axino only if the mixing generated from the superpotential (2.13) and SUSY breaking, which is of order v u,d /f a and z 0 /f a , is negligible. Otherwise the whole axino-neutralino mixing matrix has to be considered in detail. We will not discuss this case further, but point instead to the related studies in the case of the Next-to-Minimal Supersymmetric Standard Model with a singlino LSP [27].
In the DFSZ models, there are also axino tree-level Yukawa interaction terms to the quark and squark with a coupling of the order of m q /f a below the EWSB scale, with the Higgs and higgsino with coupling µ/f a even above the EWSB scale. These tree-level interactions are not suppressed by a gauge coupling or a loop factor, as the QCD anomaly term has, and thus they give the dominant contribution to axino production through the decay and/or scattering processes at low reheating temperature T R . At high reheating temperatures above EWSB instead the SU (3) c anomaly coupling is absent and the SU (2) L anomalous interaction dominates the production.

The production of axinos
There are two efficient and robust ways of populating the early Universe with axinos. Firstly, they can be produced through scatterings and decays of particles in thermal equilibrium. This mechanism, which we call thermal production (TP) depends on the reheating temperature after inflation. The other mechanism, which is independent of T R , involves non-thermal production (NTP) of axinos from the decay of the NLSP after it has frozen out from the plasma. Note also that, even though squarks are normally not the NLSP and remain in thermal equilibrium, for T R ∼ < m q and large gluino mass, axino yield from decay processes q → qã can dominate the abundance [3]. Additionally the decay of inflaton or moduli can produce axinos but such contributions are very model dependent and won't be considered here.
Thermally produced axinos in the keV mass range were considered as warm dark matter (WDM) in Ref. [28] and much lighter ones as hot DM in Ref. [22]. In Ref. [1] it was shown that axinos from neutralino NLSP decays can be a natural candidate for CDM and this was extended in Ref. [2] to include TP. If the axino mass is between around an MeV to several GeV, the correct axino CDM density is obtained when T R is less than about 5 × 10 4 GeV [2]. At higher T R and lower (∼ keV) mass, axinos could constitute WDM. As a digression, we note that, when axinos are very heavy, and it is the neutralinos that play the role of the LSP, their population from heavy axino decays could constitute CDM [29], leading in particular to the possibility of TeV-scale cosmic ray positrons, as pointed out in Ref. [30]. The possibility of either WDM or very heavy axinos is not discussed in this paper.
The interactions leading to CDM axinos were extensively studied in terms of cosmological implications in Refs. [2][3][4] and collider signatures in Refs. [8,[11][12][13]. If the LHC does not confirm the decay of heavy squarks or gluino to a lighter neutralino, the axino CDM idea with the R-parity conservation can not be saved unless some other mechanisms such as an effective SUSY is introduced [31].
In general the couplings of the axino to gauge and matter fields are analogous to those give in Eq. (2.9) Here we normalize the PQ scale by 2V a /N → f a which sets c 3 = 1 for the QCD anomaly coupling and therefore defines the coefficients C aW W and C aY Y in the axion models considered above.

Thermal production
At sufficiently high temperatures ( ∼ > 10 9 GeV) axinos can reach thermal equilibrium with SM particles and their superpartners. However, assuming that a subsequent period of cosmic inflation dilutes the population of such primordial axinos (and that they are not produced directly in inflaton decay), a postinflationary axino population comes firstly from the hot thermal bath. If the reheating temperature is very high, above the decoupling temperature of axinos (∼ 10 9 GeV), their relic number density reaches again thermal equilibrium and is the same as that of photons. In that case axinos must be so light ( ∼ < 1 keV) that they become warm or hot DM [28]. On the other hand, when the reheating temperature is below the decoupling temperature, axino number density is much smaller than that of photons and its time evolution is well described by the Boltzmann equation without backreaction [2].
In the KSVZ model, at high temperatures, the most important contributions come from two-body scatterings of colored particles into an axino and other particles. At lower reheating temperatures, on the other hand, the decay of squarks or gluinos can dominate the production of axinos [3,4]. In Ref. [2], an effective gluon thermal mass was introduced to regulate the infrared divergence in the scattering cross-section. Subsequently, a more consistent calculation using the hard thermal loop (HTL) approximation was applied to the axino production in Ref. [5]. However, the HTL approximation is valid only for small gauge coupling, g ≪ 1, which corresponds to the reheating temperature T R ≫ 10 6 GeV. Below 10 6 GeV the HTL approximation becomes less reliable [5]. In fact, the production rate becomes even negative at g 3 1.2, and therefore the result becomes unphysical. Strumia tried to improve this behaviour by using the full resummed finite-temperature propagators for gluons and gluinos in the loop [6]. This procedure includes axino production via gluon decays, which is kinematically allowed by thermal mass, and results in an enhancement compared to the HTL approximation. However, his method is gauge-dependent in the next-to-leading order [32], indicating that not all the contributions of that order are included, and therefore also does not give a fully satisfactory result in the large gauge coupling regime. For these reasons, we believe that it is worth pursuing an alternative method of computing the rate at large couplings and we apply the effective mass approximation as the only way to evaluate the axino thermal production at large coupling g and low reheating temperature T R 10 4 GeV. Even though this method is also not gauge invariant, it captures relevant physical effects like plasma screening and gives positive and physical cross-sections for each single scattering process.
As stated earlier, in the previous studies of TP of relic axinos the contributions from SU (2) L and U (1) Y gauge interactions were neglected, but for completeness, we will discuss here all SM gauge groups explicitly in order to examine their possible effects.
To start with, in evaluating the contributions from the strong interactions we will follow the method used previously in Ref. [2] to obtain the (dominant) axino TP cross-section. We will further update and correct the tables presented there by following Ref. [6] to include the previously ignored dimension-5 axino-gluino-squark-squark interaction term. This term changes only the contribution from the processes H and J in Table 2, while it does not affect the other terms. The processes B, F, G and H, with gluon t-or u-channel exchange, are infrared-divergent and thus, following Ref. [2], are regularized with the inclusion of an effective gluon thermal mass in the gluon propagator. This method gives always positive definite values for the single cross sections. The full results and the method how they are obtained is described in the Appendix. In Table 2, we give for simplicity only the first two leading terms in the expansion for s/m 2 eff >> 1. However the logarithm in the approximate formulae of Table 2 gives unphysical negative value for s < m 2 eff . Therefore in the numerical calculation we have used the full formulae which are positive definite for all values of s listed in the Appendix.
The total cross-sections σ n , where n = A, . . . , J, can be written as whereσ n (s) denotes the cross-section averaged over spins in the initial state and are given in Table 2.
The relevant multiplication factors are also listed: n s (the number of initial spin states), n F (the number of chiral multiplets) and η i (the number density factor, which is 1 for bosons and 3/4 for fermions). We assume particles in thermal equilibrium to have a (nearly) Maxwell-Boltzmann distribution, proportional to η i , and neglect Fermi blocking or Bose-Einstein enhancement factors, which are close to one at these temperatures. We restrict ourselves to temperatures above the freeze-out temperature of SM superpartners involved, so that the approximation of thermal equilibrium is always satisfied. 7 Finally, in Table 2 the group theory factors f abc and T a jk of the gauge group SU (N ) satisfy the relations a,b,c |f abc | 2 = N (N 2 − 1) and a jk |T a jk | 2 = (N 2 − 1)/2. Next we move to include contributions from the SU (2) L and U (1) Y gauge interactions. The relevant axino-gaugino-gauge boson and axino-gaugino-sfermion-sfermion interaction terms, in view n Process σ n ns n F Hg a +q k →ã +q j 1 Table 2. The cross-section for each axino thermal-production channel involving strong interactions. The particle masses are neglected except for the plasmon mass m eff . The H and J entries with an asterisk in the third column are changed due to including the missing term and cross-sections or nF in the others processes (A,B,D,E,F, and I) are corrected from those of Ref. [2]. The logarithm in the approximate formulae in this of Eq. (2.11), are given by where the terms proportional to α 2 come from the SU (2) L and the ones proportional to α Y from the U (1) Y gauge groups. C aW W and C aY Y are the model-dependent couplings for the SU (2) L and U (1) Y gauge group axino-gaugino-gauge boson anomaly interactions, which is defined after the standard normalization of f a , as in the first line for SU (3), as stated below Eq. (3.1). Here α 2 ,W , W µν and α Y ,Ỹ , Y µν are the gauge coupling, the gaugino field and the field strength of SU (2) L and U (1) Y gauge groups, respectively.f D represents the sfermions of the SU (2)-doublet andf are the sfermions carrying the U (1) Y charge. We start from Table 2 and replace quark triplets with SU (2) L -doublets with corresponding group factors. For the abelian U (1) Y factor, the processes A, B, and F vanish and we can replace |T a jk | 2 with the square of the corresponding hypercharges. Finally, we use N F = (12,14,11) to count the matter multiplets charged under the MSSM gauge groups (SU (3), SU (2) L and U (1) Y , respectively.
We include above the SUSY breaking scale the full 1-loop MSSM running of the gauge couplings and gaugino masses.
The second term for each gauge group in Eq. (3.3) also generates three-body gaugino decays into an axino and two sfermions, assuming that the gauginos are heavy enough. The three-body decay rate of the gluino is given by where mg denotes the gluino mass, 5) and the mass of the axino has been neglected. However, the three-body decay is suppressed by an additional power of the gauge coupling constant and is kinematically allowed only when the gluino mass is larger than the sum of the two final-state squark masses. Therefore gluino three-body decay through the second term in Eq. (3.3) is subdominant to the two-body decay.
As stated in Ref. [3], an effective dimension-4 axino-quark-squark coupling can be generated at a loop level also in the KSVZ model. Here we take into account this effective Yukawa interaction, which appears at a two-loop level in the KSVZ models and a tree level (with a tiny mixing) in the DFSZ models [3], where ψ j andψ j denote the SM fermions and their superpartners.
In the KSVZ class of models, the effective coupling comes predominantly from the logarithmically divergent part of the gluon-gluino-quark loop term and is proportional to mg [3], Subleading terms have not yet been computed and may give a correction of the order of 20 − 30%, in analogy with what has recently been obtained in Ref. [8] for the effective tau-stau-axino coupling.
In the DFSZ models there exists also a tree-level axino-quark-squark coupling which is proportional to the mass of the quark [33], coming from the c 2 interaction term in Eq. (2.27), where the upper row relates to the up-type quarks and the lower row to the down-type quarks. These tree-level couplings are always smaller than the one-loop ones for the light generations, but not for the third one. In fact for the top quark, the tree-level coupling dominates if tan β ∼ < 4 for the gluino mass of 700 GeV, while the bottom quark tree-level coupling dominates if tan β ∼ > 1 for the same choice of the gluino mass. Note that, at low reheating temperatures, only the top-stop-axino coupling is important for axino thermal production. This is because the lighter stop is usually the lightest colored superpartner and remains in equilibrium to rather low temperatures, even below the EWSB scale. Similarly, there exists an effective tau-stau-axino vertex, which was first obtained in Ref. [4] and more recently re-derived in Ref. [8] in a full two-loop computation. This coupling is smaller and not important for thermal production, but it is important for the non-thermal production when the stau is the NLSP.
In the DFSZ models, there is also a tree-level axino-Higgs-higgsino coupling [47]; compare the second term in Eq. (2.13), where µ ∼ f s V 2 a /M P and f a = 2V a . It contributes to axino TP through higgsino decays in thermal equilibrium, and can be comparable to that from squark decays through Eq. (3.7) and Eq. (3.8), or even become much larger if µ is larger than the top quark mass. The axino production due to this coupling in DFSZ models has been recently investigated also in Ref. [33]. Note that the Yukawa coupling and the axino-higgsino mixing contained in the neutralino mass matrix also give rise to additional scattering channels contributing to the axino production, but we will neglect them here since such dimension-4 scatterings are usually less important than the decays [3].
We have evaluated the thermal production of axinos numerically and present the results in Figures 1 and 2 for representative values of f a = 10 11 GeV and m q = mg = 1 TeV. We do not consider here the dependence on masses, however see Ref. [3]. For different values of f a the curves move up or down proportional to f 2 a . In Fig. 1, we show the axino yield Y (where Y ≡ n/s is the ratio of the number density to the entropy density) from strong interaction in the KSVZ model. Our result obtained with the effective mass approximation is shown with the solid black line. Compared to the previous plot in Ref. [2], the inclusion of the squark decay changes the plot at low reheating temperature, while the other new squark interactions do not have any noticeable effect. There is a factor 3 difference in the abundance at high reheating temperature compared to that in Figure 2 of Ref. [2], which was a numerical error at that time and was corrected later. For comparison, the axino yield from scatterings using the HTL approximation [5] is plotted with the blue (dashed) line and Strumia's result [6] is shown with the green line.
For T R 10 4 GeV the axino abundance using the effective mass approximation increases consistently with that of Strumia. We found that the difference between the two prescriptions is of order a factor of three. In principle we could reabsorb this difference in the definition of the effective gluon mass at high temperature, which in our scheme cannot be determined self-consistently. With this tuning we could match the perturbative result at high temperature. On the other hand, doing this would require a gluon thermal mass smaller than the expression ∼ gT used in our calculations. Hence we prefer instead to consider this factor as an estimate of the theoretical error of using the effective mass approximation at high reheating temperatures. We assume that this error does not increase for reheating temperatures less than 10 4 GeV and we apply the effective mass approximation to the DFSZ model.
For lower temperatures the contributions from the decays of squarks and gluinos in thermal plasma, which were not included in Ref. [6], start playing some role. We mark those in Fig. 1 with a green solid and red dashed curves, respectively. It is known that, at reheating temperatures above superpartner masses scattering diagrams involving dimension-4 operators are usually subdominant relative to those coming from dimension-5 operators and to decay terms, and are negligible. Also the decays do not give significant contribution to the TP of axinos, apart from very low T R [3], and this is confirmed in Fig. 1. While all the above contributions are generated by strong interactions only, for comparison, we show also (as a black dashed line) the relative contribution from an out-of-equilibrium bino-like  as a function of the reheating temperature TR from strong interactions using the effective mass approximation (black). We use the representative values of fa = 10 11 GeV and m q = mg = 1 TeV. For comparison, we also show the HTL approximation (dotted blue/dark grey) and that of Strumia (green/light grey). We also denote the yield from squark (solid green/light grey) and gluino decay (dotted red), as well as out-of-equilibrium bino-like neutralino decay (dashed black). Here we used the interactions in Eq. (3.3) and Eq. (3.7) for the KSVZ model. We use the same definition of reheating temperature in the instantaneous reheating approximation for the three methods. neutralino decay to an axino and a photon originally considered in Ref. [1]. It is clear that NTP is only important at very low T R , well below squark or gluino masses.
In Fig. 2, we show a contribution to the axino yield in thermal production from each SM gauge group interaction. Here we set the coefficients C aW W = 1 and C aY Y = 1 as a normalization. As shown in the figure, the contributions from scatterings due to SU (2) L and U (1) Y couplings (blue dotted and green solid lines, respectively) are significantly suppressed compared to that of SU (3) c (red dashed), by a factor of 10 or more. This is because the interaction between axinos and gauge bosons are proportional to a gauge coupling-squared so that the cross-section is σ ∝ α 3 . Thus it would be only for very large (and perhaps unnatural) values of the effective couplings C aW W , C aY Y that these channels could become comparable to the QCD contribution. To give an order of magnitude estimate of these effects, we included the SU(2) and U(1) contributions with a normalized value in figure 2. For different values of C aW W and C aY Y the curves move up and down. We note that in general SUSY breaking effects in the leptonic sector may bring a modification of the couplings here considered. The situation here is different from the case of gravitino production since the interactions of the gravitino to the three gauge groups are of the same order: the spin-3/2 gravitino component couples in fact universally, while the goldstino component proportionally to the gaugino masses. We therefore conclude that the QCD contribution is strongly dominant in the KVSZ models and so the axino production at high T R is practically model independent as long as the number of heavy PQ charged states can be absorbed into the definition of f a . However, at low T R the thermal production from scatterings becomes strongly suppressed by the Boltzmann factor. In the region where T R ∼ < 100 − 1000 GeV, axino production due to the decays of gaugino, squarks or neutralinos become important. Actually, the lightest neutralino decay via U (1) Y couplings becomes dominant in the very low reheating temperature regime since the number density of the heavier colored particles becomes strongly suppressed by the Boltzmann factor there. On the other hand the neutralinos, depending on their composition and the supersymmetric spectrum, can freeze-out with a still substantial number and then give rise also to non-thermal axino production, as we have seen in Fig. 1. This contribution to the axino yield is usually more important than the one due to neutralino decays in equilibrium, which is proportional to C aY Y and typically below 10 −12 .
For the case of the DFSZ instead, the role of the QCD interaction is played by the SU (2) interaction and the dominant decay term above the EW symmetry breaking is the higgsino one instead of the squark one.
Our results for the total thermal production yield for both KVSZ and DFSZ type of models can be  seen in Figure 3. 8 There we show the KSVZ model with different values of the C aY Y coupling. In the case of non-zero C aY Y the contribution from neutralino decay in equilibrium can be clearly seen for T R ∼ 10 GeV and it is very suppressed. Moreover we give also the yield for the DFSZ model in solid green, for µ = 200 GeV and the higgsino mass mh = 200 GeV. We can see that even for relatively small µ, axino production from higgsino decay dominates over the one from the anomaly terms for reheating temperature T R 10 6 GeV. The importance of axino-Higgsino-Higgss coupling in the DFSZ model was recently discussed in [33] and our result is consistent with that analysis. The abundance is so large that the CDM density can be reached with an axino mass as small as 100 keV, independently of the reheating temperature. 9 In this range of reheating temperature, the axino production from decay dominates that from scatterings. Therefore the use of the effective mass approximation or another IR screening prescription in the scattering process is irrelevant to the axino production in the DFSZ model in the range of reheating temperature where the decay dominates. For higher T R , the SU (2) L anomaly term starts dominating and the abundance is proportional to T R as in the KVSZ case, but with a smaller coefficient. In the same figure, we also mark horizontal lines corresponding to the axino  mass giving the correct DM relic density for the given relic abundance of Yã.

Non-thermal production
As stated above, axinos can be produced non-thermally in NLSP decays after they have frozen out of equilibrium. This NTP mechanism is dominant for reheating temperatures below the mass of the gluino and squarks [1,2]. In this case, the axino abundance is independent of the reheating temperature as long as the temperature is high enough for the NLSP to thermalize before freeze-out. Axino relic density from NTP is simply given by Clearly, in order to produce a substantial NTP population of axinos, the NLSP must itself have an energy density larger than the present DM density.
To see if such production is sufficient to give a dominant DM component, we need to know the yield of NLSPs after they have frozen out of the thermal plasma. For the neutralino NLSP yield, relevant processes include pair-annihilation and co-annihilation with the charginos, next-to-lightest neutralinos and sleptons. For the stau NLSP, the yield is determined by the stau-stau annihilation and stau-neutralino co-annihilation processes. A typical relic abundance is for a bino-dominated neutralino, and Yτ ≃ 0.01 × 10 −12 mτ 100 GeV , (3.12) for the stau. Note that in the latter case the NTP can produce sufficient axino abundance to explain the whole DM density only for stau masses above 1.9 TeV, which may thermalize only at accordingly higher temperatures. These two choices for the NLSP were considered in the Constrained Minimal Supersymmetric Standard Model (CMSSM) in Ref. [4], for f a < 10 12 GeV, for which even the stau lifetime is of order 1s, or less. Recently, the case of stau NLSP, including four-body hadronic decays, was discussed in Ref. [8] also for larger values of f a .
In conclusion, for neutralino NLSP, which decays mainly into an axino and a photon or a Z-boson, the NTP production is usually more efficient. For the stau NLSP, which can decay to an axino and a tau-lepton through a coupling of the type given in Eq. (3.6), the contribution is smaller, but can still be substantial. too warm (NTP from neutralino NLSP) Figure 6. The same as Fig. 4 but with DFSZ model used in Fig. 3.
Regarding other NLSPs, colored relics (or even a wino-like neutralino if it is lighter than 1.8 TeV) usually remain in thermal equilibrium so long that their number density after freeze-out becomes negligible and therefore cannot produce any substantial axino population after freeze-out [37,38].

The relic density of dark matter
For the total axino DM relic density, we apply the 3σ range derived from WMAP-7 data [39] 0.109 < Ωãh 2 < 0.113. (4.1) This produces a stripe in the parameter space and also plays a role of the upper bound on the relic density when there are additional DM components, e.g. the axion. This can be seen in Figs. 4 and 5, where we present the reheating temperature versus the axino mass plane for f a = 10 11 GeV and 5 × 10 9 GeV, respectively. We have included both the thermal and the non-thermal production contributions of axinos, the latter assuming Y NLSP = 0 (black solid) 10 −10 (green solid) and 10 −8 (red dash), denoted also as (I), (II), and (III), respectively. A typical stau and neutralino yield after freezout will lie between (I) and (III). A correct relic density of axinos, in the range given by Eq. (4.1), corresponds to the thin bands between like curves. The parameter space above the curves is excluded as giving too much relic abundance. Similar figures in the DFSZ model are presented in Figs. 6 and 7, which can be compared to the figures in Ref. [47].

Nucleosynthesis
An injection of high energy electromagnetic and hadronic particles during or after Big Bang Nucleosynthesis (BBN) epoch may disrupt the abundances of light elements. For axino DM, the lifetime of the NLSP, such as the neutralino or the stau, is typically around 1 sec, or less, and therefore constraints from the BBN are weak. However, for longer lifetimes such constraints become important [40,41]. This leads to an upper bound on the decay products of the NLSP for a given NLSP lifetime. For the stau NLSP, a bound state with 4 He severely constrains its lifetime to be less than roughly 5 × 10 3 sec [42] (although in specific cases with gravitino as DM, this can be up to an order of magnitude larger [43]). However for the parameters considered in our study, i.e. f a < 10 12 GeV, the BBN constraint can be avoided due to the small lifetime of the NLSP. For larger values of f a , on the other hand, non-trivial constraints arise, especially for the stau NLSP, as recently discussed in Ref. [8].

Structure formation
The density perturbation due to axino population is suppressed at scales below their free-streaming length. When the axino mass is larger than O( keV), thermally produced axinos become cold [2]. However, the non-thermal population of axinos from NLSP decays can still have a large velocity dispersion and can be too warm. Lyman-α [44] and reionization data [45] give a bound on the velocity of the WDM component and its fraction in the DM density. A recent analysis using the SDSS Lyman-α data [44] leads to an upper limit on the average velocity, v WDM < 0.013 km/s for pure WDM, or otherwise in the case of mixed cold/warm DM the WDM fraction is constrained to be Ω WDM /Ω DM < 0.35 in the larger velocity region.
The present velocity of thermally produced axinos is given by [46] v 0 ≃ 0.065 km/s 1 keV mã . (4.2) Therefore the above Lyman-α data implies mã 5 keV for TP axinos. This lower bound is marked in Figs. 4 and 5 with a vertical blue dashed line and an arrow.
The free-streaming velocity of axinos produced non-thermally can be obtained from the lifetime of the NLSP and the mass relations [2,4,45]. For the bino-like neutralino NLSP with C aYY = 8/3, we find  For NTP axinos and for the neutralino NLSP we find a lower limit on the axino mass mã 30 MeV and 1.5 MeV when f a = 10 11 GeV and 5×10 9 GeV, respectively. These limits are marked in Figs. 4 and 5 with a vertical blue dashed line and an arrow. For the stau NLSP the analogous lower bounds are mã 150 MeV and 7.5 MeV, respectively. We stress that these bounds apply solely if the population of axinos produced through NTP is substantial ( ∼ > 20 − 30%). In Figs. 8 and 9, we show contours of the reheating temperature in the plane spanned by the NLSP and the axino mass. 10 In Fig. 8 we have fixed f a = 10 11 GeV and assumed Y NLSP = 10 −12 (m NLSP /100 GeV), which is a typical value for bino-like neutralino NLSP. The cyan wedge in the upper right-hand corner is excluded by the overdensity of DM, while in the red wedge below it the axino is not the LSP. In Fig. 9 instead Y NLSP = 10 −14 (m NLSP /100 GeV) has been assumed, typical of stau NLSP. So long as the TP dominates, the curves of constant T R remain vertical and practically independent of mã, while as soon as the NTP becomes important, mã dependence arises leading to non-vertical curves. For example, with a bino-like neutralino NLSP of 100 GeV, as in Fig. 8, the NTP contributes only up to 10 % of the axino LSP DM density. For stau NLSP with small abundance Y ∼ < 10 −13 , as given in Eq. (3.12), NTP is always subdominant so that the contours remain vertical. In this case the bounds coming from the free-streaming velocity of the NTP axinos are absent.

Conclusion
We have performed an updated analysis of the relic axino production, taking into account some new calculations that have appeared after the initial study [1,2], especially [5,6] for the thermal part, and compared them with our own results to explore the question of uncertainty and model-dependence.
We have found that the uncertainty has not really decreased after the latest calculation [6]. This is probably not surprising since the QCD coupling is large and the convergence of the perturbative 10 Similar figures are shown in Ref. [8]. series is quite slow. Comparing the different results, we estimate the uncertainty in the relic density of axinos produced thermally to be still of order a factor of 10 or so at T R ∼ 10 4 GeV, and to that one has to add also some possible (unknown) contributions due to non-perturbative effects. Our result lies above both estimates given in Refs. [5,6], and this seems natural since we included subleading terms in m 2 eff /s, which do indeed increase the cross-sections for single channels ensuring their positivity in the whole range of integration, even in the limit of very large gauge coupling.
Regarding the model dependence, our conclusions are more optimistic: in the KSVZ-type the QCD anomaly term strongly dominates the axino thermal production mechanism, apart from the case of small reheating temperatures where sparticle decay contributions start playing the dominant role. The inclusion of the additional anomaly couplings is completely negligible, apart from unnaturally large values of the coefficients C aW W and C aY Y . Instead for DFSZ-type models, the Yukawa interaction dominates around the weak scale and can give the main production mechanism of axinos making it independent of the reheating temperature. Therefore the axino abundance is free from the uncertainty in the method used for the IR-divergence. At large temperatures, it is the EW anomaly term that dominates, giving a lower abundance than in the KSVZ case.
For both models, also the non-thermal production via NLSP decay can produce the required axino DM density, if the NLSP decouples with a sufficiently large abundance. But for this mechanism to dominate, the reheating temperature has to be very low and the axino and NLSP masses not too hierarchical. We find that, interestingly enough, a light bino NLSP decaying out of equilibrium can still produce the whole DM density at the cost of a very low reheating temperature (but sufficiently high for NLSP thermalization). For the stau case instead, it is quite unlikely that the NTP can dominate, unless the stau NLSP yield after freeze-out is unusually large.
During the final stages of completion of this work, the analysis Ref. [47] appeared which discusses in details axino couplings and finds a non-trivial momentum dependence in the one-particle irreducible one-loop axino-gluon-gluino couplings. The coefficient C 1P I of these interactions is suppressed when the external particle momentum is much larger than the mass M of the PQ-charged fermions in the loop. Due to this effect, the authors claim a suppression of order M 2 /T 2 of the axino production from the dimension-5 operators for the DFSZ case and for extremely small KSVZ quark masses (with Yukawa couplings less than 10 −5 , i.e. for the heavy quark mass less than 10 6 GeV for f a ≃ 10 11 GeV).
We investigated such suppression by inserting their C 1P I coupling in the relevant diagrams and we obtain different suppressions depending on the type of Feynman graph and a strong dependence on the IR regulator contained there, in most cases the gluon thermal mass. In particular, we find no suppression at all for the t-channel gluon exchange for vanishing gluon mass. Since graphs with the one loop C 1P I coupling and a gluon thermal mass insertion arise at lowest order at two loops, probably a full investigation of the two-loop diagrams in thermal field theory is needed to resolve this issue.
Note in any case that even without suppression, the DFSZ axino production is dominated by the decay term up to temperatures of the order of 10 6−7 GeV. For the KVSZ case, we show in Fig. 10 as Neutralino NLSP) Figure 10. The same as Fig. 4 but showing as well in magenta the yield suppression found in Ref. [47] for different values of the heavy quark masses.
violet lines how the yield changes according to Ref. [47] for small heavy quark masses, m Q = 10 6 GeV and m Q = 10 5 GeV. The Cold DM axino, on which our present study is based, is practically not affected. For large fermion masses in the loop, M > T , our anomalous couplings coincide fully with the C 1P I in Ref. [47] and our results are in perfect agreement.

A Appendix
In this Appendix, we present the computation of the axino relic density in detail. The time evolution of the axino number density, nã, is described by the Boltzmann equation: dnã dt + 3Hnã = i,j σ(i + j →ã + · · · )v rel n i n j + i Γ(i →ã + · · · ) n i .
Here H is the Hubble parameter, H(T ) = (π 2 g * )/(90M 2 P )T 2 , where g * is the number of effective massless degrees of freedom. The first term in the r.h.s. of Eq. (A.1) is the axino production from the scattering process of particles i and j in the thermal bath with cross section σ(i + j →ã + · · · ) and relative velocity v rel . n i is the i-th particle's number density in the thermal bath. The second term in the r.h.s. of Eq. (A.1) is the axino production from the decay of i particle which are in thermal bath with the decay rate Γ(i →ã + · · · ). · · · denotes the thermal averaging including averaging over initial spins and summing over the final spins. Here we neglect the inverse processes since the axinos are decoupled from the thermal bath and their number density is very small.
Using the axino yield defined as where s = (2π 2 /45)g s * T 3 is the entropy density, and normally g s * = g * in the early Universe, the solution of the Boltzmann equation is where Y scat i,j = TR 0 dT σ(i + j →ã + · · · )v rel n i n j sHT , dT Γ(i →ã + · · · ) n i sHT . (A.4) The explicit formulae for the thermal average are given in Ref. [48]. The relevant 2-body scattering cross sections are summarized in table 2, keeping in mind that the physical cross-sections are Eq. Among the processes, B,F,G and H have an infrared (IR) divergence due to the massless gluon exchange in the t-or u-channel. To regularize this IR divergence, we introduce, only in the t or u-channel gluon propagators, an effective thermal gluon mass m 2 eff = g 2 T 2 , which is generated by plasma effects [2,49]. Then we obtain finite and always positive cross sections. For example, the squared matrix element for process B with using massless gluon propagator is |M B | 2 = −4(2s 2 + 2ts + t 2 ) t |f abc | 2 . (A.7) However introducing effective thermal gluon mass in the t-or u-channel, it changes to |M B | 2 = −2t(3m 4 eff − 6m 2 eff t + 4s 2 + 4ts + 2t 2 ) (t − m 2 eff ) 2 |f abc | 2 , (A. 8) which is positive definite by definition. Thenσ B (s) is obtained after integration over t as given in the Eq. (A.6). Here we list the full cross-sectionsσ n for these processes. These formulas give back the expressions in Table 2, in the limit of large s/m eff . We mention here that those expressions in Table 2 provide unphysical result for small s/m eff and give negative contribution to the thermal averages. In fact substituting s ≃ T 2 , m 2 eff ∼ g 2 T 2 , one finds e.g. for the B processσ B ∝ log[1/g 2 (T )] − 7 4 < 0 for g 2 (T ) > 0. 17 (A.13) and this holds for all the interesting region up to the GUT scale where g 2 ∼ 0.5. So even if the total thermal cross-section may result positive, it is lower than the real one due to the negative contribution of these IR-divergent channels. On the other hand the full expressions above from Eq. (A.9) to Eq. (A.12) are positive for any value of s. We can see that easily in the limit of large m 2 eff /s because those reduce to: so they have the correct asymptotic behaviour corresponding to screening and physical decoupling of the intermediate gluon channels for large gluon mass. 11 So with these formulas the single cross-sections are always positive and the thermal average larger than in the previous estimates.
In our numerical computations we use the full formulae above for the IR divergent processes, while for the other processes, A, C, D, E, I, J, we use the expressions are given in Table 2.