Gravitino thermal production revisited

We calculate the gravitino production rate, computing its one-loop thermal self-energy. Gravitino production processes that do not result through thermal cuts of its self-energy, have been identified and taken into account. Correcting analytical errors and numerical approximations in the previous calculations, we present our result. This deviates from the latest estimation by almost 10\%. More importantly, we provide a convenient formula, for calculating the gravitino production rate and its thermal abundance, as a function of the reheating temperature of the Universe.

Introduction.-Extensions of the Standard Model (SM) in the context of supergravity provide us with a dark matter (DM) candidate particle the gravitino, the superpartner of graviton. It interacts purely gravitationally with other particles and thus naturally escapes direct or indirect detection, as the current experimental and observational data on DM searches suggest. Therefore, the precise knowledge of its cosmological abundance is essential to apply cosmological constraints on these models. After inflation, gravitinos may be produced in various ways: (i) non-thermally, from the inflaton decays, (ii) much later around the Big Bang Nucleosynthesis epoch, through the decays of unstable particles and (iii) last but not least, thermally, when the Universe cools down from the reheating temperature (T reh ) until now.
The effort in calculating the thermal gravitino abundance, using various techniques, methods and approximations, spans over almost the last four decades. Since the gravitinos are mainly thermally produced at very high temperatures, the effective theory of light gravitinos, the so-called non-derivative approach, involving only the spin 1/2 goldstino components, was initially used. In this context, as some of the production amplitudes exhibit infrared (IR) divergences, they were regularized by introducing either a finite thermal gluon mass or an angular cutoff. Thus, in [1] the basic 2 → 2 gravitino production processes had been tabulated for the first time and calculated, see Table I. This calculation was further improved in [2,3].
As the Braaten, Pisarski, Yuan (BPY) method [4,5] succeeded in calculating the axion thermal abundance, in [6] it was further applied to the gravitino, motivated by the fact that the gravitino like axion, interacts extremely weakly with the rest of the spectrum. Although in [7] the previous IR regularization technique was used, in [8,9] the BPY method was employed, taking in addition into account the contribution of the spin 3/2 pure gravitino components.
Eventually, in [10] the calculation method improved significantly. There it was argued that the basic requirement to apply the BPY prescription, that is g 1, where g is the gauge coupling constant, is not satisfied in the whole temperature range of the calculation, especially if g is the strong coupling constant g 3 . Therefore, the authors calculated the one-loop thermal gravitino self-energy numerically beyond the hard thermal loop approximation, with the benefit that this incorporates also the 1 → 2 processes besides the 2 → 2 ones. More importantly, it was noticed that the so-called subtracted part, i.e. pieces of the 2 → 2 squared amplitudes for which the self-energy may not account for, are IR finite. The main numerical result in [10] on the gravitino production rate differs significantly, almost by a factor of 2, with respect to the previous works [9,11]. Unfortunately, in [10] the main analytical results appear to be inadequate. In particular, in Section IV(A) the equations on the self-energy contribution for the gravitino production rate, seem to be inconsistent even dimensionally. In addition, the numerical estimation of this self-energy, was computed only inside the light-cone, due to the limited computation resources of that time. Furthermore, two out of the four non-zero subtracted parts in the corresponding Table I in [10], turn out to be zero.
Motivated by these, in this Letter we recalculate the thermally corrected gravitino self-energy and we compute it without numerical approximations, at the one-loop level. Eventually, since our final result for the gravitino production rate, like in [10] is numerical, following [12] we present an updated handy parameterization of this. Our final result differs from that shown in [10] approximately by 10%. Moreover, we calculate the gravitino thermal abundance and discuss possible phenomenological consequences.
The setup.-As the gravitino is the superpartner of the graviton, its interactions are suppressed by the inverse of the reduced Planck mass M P = (8π G) −1/2 . Hence, the dominant contributions to its production, in leading arXiv:2010.14621v1 [hep-ph] 27 Oct 2020 TABLE I. Squared matrix elements for gravitino production in SU (3)c in terms of g 2 3 Y3/M 2 P assuming massless particles, Y3 = 1 + m 2 g /(3m 2 G ), C3 = 24 and C 3 = 48.
order of the gauge group couplings, are processes of the form a b → c G, where G stands for gravitino and a, b, c can be three superpartners or one superpartner and two SM particles. The possible processes and the corresponding squared amplitudes in SU (3) c are given in Table I, where for their denotation by the letters A to J we follow the "historical" notation of [1]. In SU (3) c the particles a, b and c could be gluons g, gluinosg, quarks q or/and squarksq. Analogous processes happen in SU (2) L or the unity is related to the 3/2 gravitino components and the rest to the 1/2 goldstino part. For the calculation of the spin 3/2 part in the amplitudes, following [10], we have employed the gravitino polarization sum where Ψ µ is the gravitino spinor and P its momentum. As in [10], for the goldstino 1/2 part the non-derivative approach is used [11,13]. The result for the full squared amplitude has been proved to be the same, either in the derivative or the non-derivative approach [14]. The Casimir operators in Table I a,i,j denotes the sum over all involved chiral multiplets and group indices. f abc and T a are the group structure constants and generators, respectively. Processes A, B and F are not present in U (1) Y because C 1 = 0. The masses for the particles a, b and c are assumed to be zero. In the third column of Table I we present for each process the square of the full amplitude, which is the sum of individual amplitudes, where the indices s, t, u indicate the diagrams which are generated by the exchange of a particle in the corresponding channel and the index x stands for the diagram involving a quartic vertex. The so-called D−graph is illustrated in Fig. 1 for the case of the gluino-gluon loop. We call it D−graph following the terminology of [10]. Its contribution is the sum of the squared amplitudes for the s, t and u-channel graphs, The subtracted part of the squared amplitudes is the difference For the processes B, F, G and H the corresponding amplitudes are IR divergent. For this reason we follow the more elegant method comprising the separation of the total scattering rate into two parts, the subtracted and the D−graph part. It is worth to mention that for the processes with incoming or/and outgoing gauge bosons, we have checked explicitly the gauge invariance for |M X,full | 2 . On the other hand, we note that |M X,sub | 2 is gauge dependent [15].
To sum up, the gravitino production rate γ 3/2 consists of three parts: (i) the subtracted rate γ sub (ii) the D−graph contribution γ D and (iii) the top Yukawa rate γ top , Below, these three contributions are discussed in detail.
The subtracted rate.-In the fourth column of Table I we present the so-called subtracted part (4), which is the sum of the interference terms among the four types of diagrams (s, t, u, x), plus the x-diagram squared, for each process. The subtracted part is non-zero only for the processes A and B. Noting that in [10] the subtracted part for the processes H and J is also non-zero, we assume that the authors had used the squark-squark-gluino-goldstino Feynman rule as given in [16], where a factor γ 5 is indeed missing. In contrast we are using the correct Feynman rule as given in [11].
To calculate the subtracted rate for the processes a b → c G, we use the general form In the temperature range of interest all particles but the gravitino are in thermal equilibrium. For the gravitino the statistical factor f G is negligible. Thus 1 − f G ∼ 1, as it is already used in (6). Furthermore, backward reactions are neglected. In addition, the simplification 1 ± f c ∼ 1 is usually applied, making the analytic calculation of (6) possible. In our case there is no such reason. We keep the 1 ± f c and consequently we proceed calculating the subtracted rate numerically [17].
The contribution of the processes A and B, for each gauge group, can be read from Table I as (8) In (8) a factor 1/2 is already included for the process A due to the 2 identical incoming particles. Substituting (8) in (6), the subtracted rate is obtained as (9) The numerical factors, calculated by using the Cuba library [18], are C s BBF = 0.25957 × 10 −3 and C t BFB = −0.13286 × 10 −3 . The subscripts B and F specify if the particles are bosons or fermions respectively and the superscripts determine if the squared amplitude is proportional to s or t. It is easy to see that our result for the subtracted part unlike in [10] is negative. This is not unphysical, since the total rate and not the subtracted one is bound to be positive.
The D−graph contribution.-As it has been discussed above, Eq. (3) describes the relation between the D−graph and the sum of the squared amplitudes for the s, t, and u channels. In the D−graph contribution we will implement the resummed thermal corrections to the gauge boson and gaugino propagators [19]. Although in Fig. 1 the gluino-gluon thermal loop is displayed, the contributions of all the gauge groups have been included in our analysis. The momentum flow used to calculate the D−graph can be depicted in Fig. 1. That is G(P ) → g(K) +g(Q), with P = (p, p, 0, 0) , K = (k 0 , k cos θ k , k sin θ k , 0) and Q = (q 0 , q cos θ q , q sin θ q , 0), where θ k,q are the polar angles of the corresponding 3momenta k, q in spherical coordinates.
The spectral functions ρ L,T and ρ ± can be found in Eqs. (3.7) in [10]. The thermally corrected one-loop self-energy for gauge bosons, scalars and fermions that we have used in calculating these spectral functions, can be found in [23][24][25][26][27][28]. Comparing (13) with the corresponding analytical result given in Eqs. (4.6) and (4.7) in [10], one can notice that they differ on the overall factor and on the number of independent phase-space integrations. Our analytical result has been checked using various frames for the momenta flow into the loop.
The top Yukawa rate.-The production rate resulting from the top-quark Yukawa coupling λ t is given by [10] γ top = where A t is the trilinear stop supersymmetry breaking soft parameter and C s BBF = 0.25957 × 10 −3 . Since this contribution stems from the process squark-squark → higgsino-gravitino, only the collision term C s BBF is involved.  The parameterization of the result.-Following [12] we parametrize the results (9) and (13) using the gauge couplings g 1 , g 2 and g 3 . Thus where the constants c N and k N depend on the gauge group and their values are given in Table II. In Fig. 2 we summarize our numerical results for the gravitino production rates divided by Y N T 6 /M 2 P . Especially, for the case of the top Yukawa contribution, in Y N the m 2 λ N has to be replaced by A 2 t . The colored solid curves represent the SU (3) c (red), SU (2) L (blue) and U (1) Y (green) The gravitino production rates divided by YN T 6 /M 2 P . The solid curves represent in order, the total rate (black) given by (5), the SU (3)c (red), SU (2)L (blue) and U (1)Y (green) rates given by (15) and the top Yukawa rate (purple) given by (14). The upper dashed curve is the total production rate obtained in [10]. The top Yukawa coupling λt has been taken equal to 0.7 so that our result can be directly compared with that in [10]. rates given by (15) and the top Yukawa rate (purple) given by (14), while the black solid curve is the total result given by (5). The dashed black curve corresponds to the total result from [10]. For the sake of comparison, we have also chosen λ t = 0.7.
Despite the analytical and numerical discrepancies with [10], it is interesting that our result for the total gravitino production rate is only 5 − 11% smaller than that in [10]. Being unable to explain this quantitively in details, we assume that the aforementioned differences have opposing effects on the total result. For convenience, in Fig. 2 universal gauge coupling unification is assumed at the grand unification (GUT) scale ∼ 2 × 10 16 GeV, but certainly the result in (15) can be used independently of this assumption. Eq. (15) along with the numbers in Table II is the main result of this Letter.
The Gravitino abundance.-The Boltzmann equation for the gravitino number density n 3/2 iṡ where H is the Hubble constant and the dot denotes time differentiation. The gravitino abundance is defined as with n rad = ζ(3)T 3 /π 2 . Substituting (17) into (16) we obtain that the gravitino abundance for T T reh is given by [12] where g * s are the effective entropy degrees of freedom in the Minimal Supersymmetric SM. In [12] it was assumed that the inflaton decay is instantaneous as is the thermalization of the Universe. Following the latest data from the Planck Satellite, the cosmological accepted value for the DM density in the Universe is Ω DM h 2 = 0.1198 ± 0.0012 [29]. Assuming that the thermal gravitino abundance amounts to the observed DM, we obtain that where ρ cr = 3 H 2 0 M 2 P is the critical energy density. H 0 = 100 h km/(s Mpc) is the Hubble constant and T 0 = 2.725 K the cosmic microwave background temperature today. The entropy degrees of freedom at the associated temperatures are g * s (T 0 ) = 43/11 and g * s (T reh ) = 915/4. The last number equals to the effective energy degrees of freedom for H(T reh ) in the Minimal Supersymmetric SM too. Fig. 3 illustrates the 3σ regions resulting from (19), for various values of m 1/2 . In this figure the trilinear coupling A t has been ignored and the top Yukawa coupling is λ t = 0.7, as previously. As before, gauge coupling unification it is assumed, as well as a universal gaugino mass m 1/2 at the GUT scale.
For large gravitino mass the reheating temperature is m 1/2 independent, as the characteristic factor m 2 λ N /(3m 2 G ) becomes negligible for m 1/2 m G . Assuming that m 1/2 750 GeV, as it is suggested by the recent LHC data [30,31] on gluino searches, from Fig. 3 we infer that for maximum T reh 10 9 GeV the corresponding gravitino mass is m 3/2 550 GeV. Allowing for a reheating temperature an order of magnitude smaller, T reh 10 8 GeV, for the same gravitino mass, m 1/2 can go up to 3 − 4 TeV.
Conclusions.-In this Letter we have calculated the gravitino thermal abundance, using the full one-loop thermally corrected gravitino self-energy. Having rectified the main analytical formulae for the gravitino production rate, we have computed it numerically without approximation. We offer a simple and useful parametrization of our final result. In the context of minimal supergravity models, assuming gaugino mass unification, we have updated the bounds on the reheating temperature for certain gravitino masses. In particular, saturating the current LHC gluino mass limit mg 2100 GeV, we find that a maximum reheating temperature T reh 10 9 GeV is compatible to a gravitino mass m 3/2 500−600 GeV.
It should be noted that, trying to constrain the reheating temperature by applying the cosmological data on gravitino DM scenarios, implies that thermal leptogenesis is assumed as a mechanism for generating baryon asymmetry. In any case, there are many alternative models for baryogenesis. In addition, as it has been pointed out before, the thermal gravitino abundance is in general a part of the whole DM density and the inclusion of other components will affect the phenomenological analysis. Moreover, our results for the gravitino production rate and its thermal abundance can be used in a general supersymmetric framework, independently of the details of the supersymmetry breaking scenario.