The Effects of QCD Equation of State on the Relic Density of WIMP Dark Matter

Weakly Interactive Massive Particles (WIMPs) are the most widely studied candidate particles forming the cold dark matter (CDM) whose existence can be inferred from a wealth of astrophysical and cosmological observations. In the framework of the minimal cosmological model detailed measurements on the cosmic microwave background by the PLANCK collaboration fix the scaled CDM relic density to $\Omega_{c}h^2=0.1193\pm0.0014$, with an error of less than 1.5%. In order to fully exploit this observational precision, theoretical calculations should have a comparable or smaller error. In this paper we use recent lattice QCD calculations to improve the description of the thermal plasma. This affects the predicted relic density of"thermal WIMPs", which once were in chemical equilibrium with Standard Model particles. For WIMP masses between 3 and 15 GeV, where QCD effects are most important, our predictions differ from earlier results by up to 9% (12%) for pure S-wave (P-wave) annihilation. We use these results to compute the thermally averaged WIMP annihilation cross section that reproduces the correct CDM relic density, for WIMP masses between 0.1 GeV and 10 TeV.


INTRODUCTION
Assuming Newtonian gravity, and its extension into General Relativity, describes gravitational forces correctly, astronomical and cosmological observations show that most of the matter in our Universe is a non-luminous, neutral form of matter, called dark matter (DM). Quantitatively, it accounts for ∼ 85% of all matter [1][2][3]. There are many possible DM candidates [4]. Among those, weakly interactive massive particles (WIMPs) have been studied in most detail. There are several reasons for the popularity of WIMPs. If they once were in full chemical equilibrium with the particles of the Standard Model (SM), which is true if the largest temperature after the most recent period of entropy production exceeded about 5% of the WIMP mass, the WIMP relic density can be calculated unambiguously for a given cosmological model, independent of initial conditions, using only particle physics quantities (masses and couplings, or cross sections) as input. This calculation yields approximately the observed relic density for roughly weak-strength WIMP annihilation cross sections. This not only hints at a deep connection between DM and extensions of the SM addressing some of its shortcomings [5][6][7][8][9], it also allows various ways in which the existence of WIMP DM can be probed. Unfortunately these probes have so far not yielded an unambiguous signal. Successful WIMP candidates must therefore not only satisfy the relic abundance constraint, but also constraints coming from direct [10][11][12] and indirect [13][14][15][16][17] detection experiments.
In standard cosmology it is assumed that the comoving entropy density remained constant since the epoch when WIMPs were in full thermal equilibrium with SM particles. The quantity of interest is then the ratio Y χ of the WIMP number density n χ and the entropy density s. In order to compute the temperature dependence of n χ , one thus needs to know the precise temperature dependence of s. In addition, the expansion rate, described by the Hubble parameter H, is proportional to the square root of the total energy density ρ; hence the temperature dependence of ρ also has to be known. In early analyses [5,18] s(T ) and ρ(T ) were calculated ignoring all interactions between SM particles, i.e. treating the thermal plasma as a relativistic free gas. This is quite a good approximation for most temperatures. However, it was realized early on [18,19] that this approach fails for temperatures near the deconfinement transition, i.e. the transition from hadronic (pions, kaons, . . . ) to partonic (quarks, gluons) degrees of freedom. More recently it was realized that also for some range of temperatures above this transition the interactions between quarks and gluons should not be ignored [20,21].
In this paper we carefully model the effect of strong interactions on the temperature dependence of the energy and entropy densities. We exploit recent calculations performed in the framework of the Hadron Resonance Gas (HRG) model well below the deconfinement transition [22], smoothly matching to results from lattice QCD (LQCD) at T = 100 MeV [23]. We find that this can change the predicted relic density by more than 5% compared to earlier treatments [19][20][21] of strong interaction effects on the thermal plasma.
We then use our improved treatment of strong interaction effects to determine the value of σv , the thermally averaged product of WIMP annihilation cross section times relative velocity of the annihilating WIMPs, that is required to reproduce the correct relic density, under the assumption that σv is independent of temperature. As already noted in ref. [24], depending on the WIMP mass, this product can deviate by up to a factor of ∼ 1.5 in either direction from the "canonical" value of 3 · 10 −26 cm 3 s −1 . This is important since experiments are now beginning to be sensitive to the canonical cross section. For example, analyses of FERMI-LAT gamma ray observations of nearby dwarf galaxies and of the diffuse gamma ray emission in the galaxy have, for some combinations of WIMP masses and WIMP annihilation final states, constrained the WIMP annihilation cross section to be lower than the canonical one [13][14][15].
The remainder of this paper is organized as follows. In Sec. 2 we review the calculation of the relic abundance of thermal WIMPs. We outline earlier estimates of s(T ) and ρ(T ), and explain in detail how we compute these quantities. In Sec. 3 we present our result for σv , and compare with previous results. In Sec. 4 we compare our σv with bounds on this quantity from indirect WIMP detection experiments and CMB anisotropies. Finally, we summarize the main results of our work in Sec. 5.

Basic Framework
The starting point of our calculation is the Boltzmann equation, which describes the time evolution of the number density n χ of a stable particle species χ that can annihilate into a pair of SM particles [18]: Here t is the cosmological time, H is the Hubble parameter, σv is the thermal average of the product of WIMP annihilation cross section and Møller velocity 1 , and n χ,eq is the density of χ particles in thermal equilibrium. We have assumed that χ particles are self-conjugate (Majorana) particles. With minor modifications our results can also be applied to particles that are not self-conjugate, such as Dirac fermions or complex scalar particle. In this case one will in general have to track the densities of particles and antiparticles separately. Since the total entropy is assumed to be conserved, S = sa 3 =const. where a is the scale factor in the Friedman-Robertson-Walker metric and s is the entropy density, it is convenient to define Y χ ≡ n χ /s. Using the definition H = (1/a)da/dt it is easy to see that this definition absorbs the dilution term 3Hn χ on the right-hand side of the Boltzmann eq.(1).
Moreover, since s as well as n χ,eq depend explicitly on temperature rather than on time, one replaces time t by x ≡ m χ /T . To that end, the entropy density is written as Here h(T ) counts the effective number of relativistic degrees of freedom (d.o.f.) that contribute to s, i.e. if all relevant particles have a common temperature T and masses T , h(T ) → i g i , where g i is the number of effective internal degrees of freedom for particle species i. For example, a massless photon contributes g γ = 2, a massless Dirac fermion (with both helicities being in thermal equilibrium) contributes g e = 4 · 7/8 = 3.5, and so on. Using d dt h(T )T 3 a 3 = 0 one then finds We thus also need to know the temperature dependence of the Hubble parameter. To that end we use the Friedmann equation H 2 = 8πG N ρ/3, where ρ and G N are the energy density and Newton's gravitational constant, respectively.
During the radiation dominated epoch, in which the decoupling of WIMPs falls, the energy density can be written as Here g(T ) counts the effective number of relativistic d.o.f. that contribute to ρ. In the limit where all particle masses can be ignored, g(T ) = h(T ), i.e. the numbers of d.o.f. defined via the entropy density and via the energy density are the same, but in general this is not the case.
Putting everything together, we have [19,25] dY where we have introduced and m χ being the WIMP mass. Numerically, λ ≡ 2.76 × 10 9 for a WIMP mass of 1 GeV and σv = 10 −26 cm 3 s −1 . The quantities g and h appearing in eq.(6) are the effective number of relativistic d.o.f. defined via the energy density, eq.(4), and entropy density, eq.(2), respectively. Note that we need to know both g and h: both appear in the definition (6) of g 1/2 * , and h also appears in the explicit expression for the scaled equilibrium density Y χ,eq . WIMPs decouple when they are essentially non-relativistic, so we can write where g χ = 2 for a neutral Majorana fermion.

The Functions g(T ) and h(T )
Most papers on the calculation of the WIMP relic density focus on the annihilation cross section. Our concern is instead an accurate calculation of the effective numbers of degrees of freedom encoded in the functions g(T ) and h(T ). To that end, we have to compute the energy density ρ(T ) and the pressure p(T ); the entropy density is then given by s(T ) = [ρ(T ) + p(T )]/T . If the single particle distribution functions f i ( k, T ) for a particle species i is known, the contribution of this species to the energy density and pressure are given by [18] ; here g i is again the number of degrees of freedom of species i, and k is the three-momentum.
The simplest case is that of a free (non-interacting) particle with mass m i . The distribution function is then the Bose-Einstein or Fermi-Dirac distribution function, which depends only on the ratio of the energy E i = m 2 i + k 2 and the temperature. Introducing the dimensionless quantities y i = E i /m i and x i = m i /T , one has [18,19]: In the denominators of the integrals +1 applies to fermions and −1 to bosons.
In the presence of strong interactions, Eqs. (9) and (10) no longer provide good approximations. Before describing our own calculation of these functions, we briefly review the state of the art, as encoded in the widely used program packages for the calculation of the WIMP relic density DarkSUSY [26], micrOMEGAs [27] and SuperIso [28].
It was realized quite early on that care has to be taken to describe the thermodynamics of the early universe around the deconfinement transition. In [29] the interactions between hadrons and between partons were approximated by simple non-relativistic potentials. Ref. [30] instead used free particles, and defined the transition temperature from the hadronic to the partonic phase as that temperature where the two calculations give the same entropy density; note that the hadronic phase gives a quickly rising h(T ), since the number of contributing hadrons quickly increases with temperature. This "hadron resonance gas model" was used in all subsequent calculations at sufficiently low temperatures, including our own.
One problem of the simple definition of the transition temperature used in ref. [30] is that it leads to a discontinuity in g(T ). Ref. [25] therefore used smooth functions interpolating between the hadronic and partonic phases. While these functions ensure that not only g(T ) and h(T ), but also their derivatives are smooth, they were not based on dynamical considerations. The authors advocated estimating the uncertainty by using two quite different values, 150 and 400 MeV, for the transition temperature. The same functions were used in Ref. [19]; the functions for a transition temperature of 150 MeV are still used by default in the computer packages mentioned above.
The first attempt to include the results of lattice QCD calculations was due to Hindmarsh and Philipsen [20]. At that time the most accurate lattice QCD calculations did not include dynamical quarks. There was some evidence that the ratio of the true pressure to the corresponding value for non-interacting particles shows little dependence on the number of quark flavors [31]. Ref. [20] therefore scaled the contribution of all strongly interacting partons by the same correction function, determined from pure glue lattice calculations [31]; at T = 1.2 GeV, these were matched to perturbative calculations [32].
The treatment by Laine and Schroeder [21] is rather similar. However, their results are based on a different set of pure glue lattice QCD calculations [33]. Moreover, they match to perturbative calculations at the much lower temperature of 350 MeV. Finally, they include the quark mass dependence up to next-to-leading order, O(g 2 ), in the perturbative expansion. In particular, they point out that charm quarks make non-negligible contributions already at temperatures of a few hundred MeV.
We now describe own treatment. At temperatures well above the electron mass, we treat all SM particles without strong interactions as free particles, i.e. we use Eqs. (9) and (10) for these particles. This includes the leptons, the electroweak gauge bosons as well as the single physical Higgs boson of the SM. Note that for the physical Higgs mass, m H 125 GeV, in the Standard Model electroweak symmetry breaking leads to a smooth cross-over, not a phase transition [34][35][36]; hence the comoving entropy density remains constant, as assumed in the derivation of eq.(5). For simplicity we compute g(T ) and h(T ) using free, massive W ± and Z bosons (with three d.o.f. each) and a single physical Higgs boson even for temperatures above the electroweak cross-over, where a more accurate treatment would use massless gauge bosons (with two d.o.f. each) and a massive complex Higgs doublet (with four d.o.f.). The difference between these treatments is very small, and will only affect the relic density of WIMPs with masses above 2 TeV.
The main focus of our work is on the effect of QCD interactions. These are most important around the deconfinement transition, which was also a smooth crossover [37] with conserved total entropy. We use the results of a recent lattice calculation with N f = 2 + 1 active flavors (meaning equal masses are used for u and d quarks, but the larger mass of the strange quark has been taken into account) [23]. This calculation covers temperatures between 100 and 400 MeV. Ref. [23] provides a parameterization of the pressure due to u, d, s quarks and gluons in this temperature range: Heret = T /T c , T c = 154 MeV being the QCD transition temperature. In this parameterization, p id = 19π 2 /36 is the ideal gas value of p/T 4 for QCD with three massless quarks. The values of the numerical coefficients appearing in eq.(11) are listed in Table 2.2. Recall that tanh(z) approaches unity for large argument z. Therefore eq.(11) automatically approaches the ideal gas value for T T c , i.e.t 1. Moreover, ref. [23] shows that eq.(11) matches quite well to available perturbative calculations at higher temperature. We therefore use this parameterization to describe the contribution from u, d, s quarks and gluons for all temperatures above 100 MeV. Once the pressure p is known, the energy density ρ can be computed from the relation between the trace of the energy-momentum tensor, also called the trace anomaly, and the pressure [23]: Since we have the analytical expression (11) for p/T 4 , we can easily obtain its derivative, and thus I(T ); the first equation in (12) then allows to compute ρ(T ) from I(T ) and p(T ). Finally, once ρ(T ) and p(T ) are known, s(T ) = [ρ(T ) + p(T )]/T can also be computed. As noted above, the effect of the charm quark is not negligible at temperatures near T c [21]. We included its contribution to the functions g and h using lattice QCD results from Table 6 of [38], using the physical ratio of charm and strange quark masses, m c /m s = 11.85 [39]. This gives us p c /T 4 , from which the charm contributions to ρ and s can be obtained as outlined in the previous paragraph. This description is valid for T ≤ 1 GeV. For larger temperatures we smoothly match to the ideal gas results (9) and (10), using a fit function like (11) with p id = 7π 2 /60 and different values for coefficients and powers to interpolate between the two regimes. This ensures that not only the functions g and h but also their first derivatives are smooth everywhere.
Bottom and top quarks contribute significantly only at high temperatures, where even QCD interactions have become relatively small. We therefore treat these quarks as free particles, with on-shell masses given by the Particle Data Collaboration [40].
At temperatures lower than T c , the thermodynamic behavior of QCD can be described by the hadron resonance gas model, in which all the hadrons and hadron resonances are considered to contribute to the thermodynamics quantities as non-interacting particles. As noted above, this has been used already in the early treatments [29,30]. Ref. [23] shows that for temperatures between 100 MeV and T c it matches well to the QCD results parameterized in eq. (11). A convenient parameterization of the trace anomaly in this model can be found in ref. [22]: with We use it to describe the contribution from strongly interacting particles for all temperatures T < 100 MeV, using cubic splines to interpolate smoothly to QCD results at T > 100 MeV. The contribution from charmed particles is negligible at these low temperatures. 2 By inverting eq. (12), the pressure can be calculated: using the numerical result p(T 0 )/T 4 0 = 0.1661 at T 0 = 70 MeV [22]. The integral in eq. (14) can easily be evaluated analytically if I(T ) is given by eq.(13). Eq. (14) thus again provides us with an analytical parameterization of the pressure, from which we can compute the energy and entropy densities as described above.
At temperatures below 1 MeV, the effect of neutrino decoupling should be included. The rate of reactions changing the ν µ and ν τ number densities actually becomes smaller than the Hubble parameter, indicating decoupling, at a temperature of several MeV. However, at first the expansion of the universe affects photons and neutrinos in the same way even after neutrino decoupling, i.e. the photon and neutrino temperatures remain the same. This changes only once e + e − pairs begin to annihilate, at T m e . Since neutrinos are already (almost) decoupled by this time, the entropy that was stored in electrons and positrons gets transferred (almost) entirely to photons, not to neutrinos. In the limit where neutrino decoupling was complete when electron decoupling began, this argument shows that for T m e the ratio of relic photon and neutrino temperatures is T γ /T ν = (11/4) 1/3 . Actually (electron) neutrinos were not completely decoupled at T m e . This can be described by writing with N eff 3.046. Note that these expressions include the contribution from the photon, with g γ = 2. Eqs. (15) and (16) are applicable for T m e , in practice for T ≤ 50 keV. As mentioned above, for T > 1 MeV we have T ν = T γ . For 50 keV < T < 1 MeV we use numerical results from Fig. 1 of [41] to determine the evolution of T ν with respect to T γ . This can then be plugged into eqs. (15) and (16) instead of T ν /T γ = (4/11) 1/3 to compute the photon and neutrino contribution to g(T ) and h(T ). Note that the temperature T is defined to be that of the photons, T = T γ .  2) and (6). The original calculation by Gondolo and Gelmini [19], based on results from ref. [25], are shown by the green dashed curves. The red dot-dashed and blue dotted curves show results from refs. [20] and [21], which are based on pure glue lattice QCD calculations. The black solid curves depict our results, which are based on lattice calculations with N f = 2 + 1 dynamical quark flavors.

RESULTS AND COMPARISON WITH PREVIOUS STUDIES
We are now ready to present some numerical results. Figure 1 shows the functions h(T ) and g 1/2 * (T ) that parameterize thermodynamic effects in the Boltzmann equation (5). The black solid curves show our results, while the green dashed, red dot-dashed and blue dotted curves show results from refs. [19], [20] and [21], respectively.
Since we include the effect of e + e − decoupling on the neutrino background, described by N eff 3.046 in eq.(15), we obtained h (T γ,0 ) = 3.9387 in our calculation; here the present photon temperature T γ,0 = 2.7255 ± 0.0006 K [42]. Our current value of h is thus slightly higher than h (T γ,0 ) = 3.9138 in ref. [20] and h (T γ,0 ) = 3.9139 in ref. [26]. Note that for a given value of Y χ (T γ,0 ) the final relic density Ω χ h 2 is directly proportional to h(T γ,0 )T 3 γ,0 . Because both deconfinement of quarks and gluons and the restoration of the electroweak gauge symmetry are associated with smooth cross-overs rather than true phase transitions, the functions g(T ) and h(T ) are smooth everywhere. As noted in the previous Section, we ensured smoothness of these functions by using cubic splines to interpolate between different temperature regions.
There clearly are some differences between the four calculations. These are most visible near the QCD deconfinement transition. Moreover, the differences are more visible in g 1/2 * , largely due to the derivative term in eq.(6), which accentuates the differences between the various treatments. We see that the older calculation [25] used in [19] overestimates g 1/2 * somewhat for T 0.1 GeV, compared to all three calculations using results from lattice QCD. The treatment of ref. [20] gives a discontinuity at T = T c , and hence a divergent derivative, yielding a formally infinite spike in g 1/2 * . The continuity has been smoothed out in micrOMEGAs [27], from which we took the numerical results for g and h. We see that this treatment still gives a prominent spike in g 1/2 * . Apart from this spike, ref. [20] predicts a smaller value of g 1/2 * in this temperature range than we do. Finally, the prediction for g 1/2 * from ref. [21] is quite close to our own result, except for some oscillatory behavior just above the QCD transition temperature.
In order to explore the effect of changes in h and g 1/2 * on the WIMP relic density, we solve the Boltzmann equation (5) numerically. If we choose the initial value x i of x to be < ∼ 10, the final result is independent of the input value Y χ (x i ). We numerically track the behavior of Y χ (x) up to x = 1, 000, at which point Y χ has become practically constant. The present scaled relic density times squared scaled Hubble constant can then be computed from where Y χ,0 = Y χ (x = 1000) and s 0 = 2891.2 cm −3 ; the subscript 0 denotes quantities evaluated at the present time [42]. Note that H 0 cancels on the right-hand side of eq.(17); this is why the relic density is usually quoted in the combination Ω χ h 2 . Numerically, Ω χ h 2 = 2.7889 · 10 8 Y χ,0 m χ /(1 GeV).
The change of the predicted WIMP relic density due to our more refined treatment of the functions h and g 1/2 * is illustrated in fig. 2. The upper frame shows results for a temperature independent σv , while in the lower frame we have assumed σv ∝ 1/x. These behaviors describe the thermally averaged cross section at small velocity away from poles (i.e. if the WIMPs cannot annihilate into any particle φ with m φ 2m χ ) and thresholds (i.e. if the WIMPs are significantly heavier than all relevant final-state particles), if the annihilation occurs from a pure S−wave and pure P −wave initial state, respectively. The latter occurs, for example, for a Majorana WIMP annihilating into light SM fermions, or for a complex scalar annihilating through s−channel exchange of a gauge boson. Not unexpectedly, we observe the largest differences for WIMP masses of a few GeV, which decouple just above the QCD transition temperature. The differences amount to up to 9% for the pure S−wave, and up to 12% for the pure P −wave.
The results of fig. 2 can be understood in more detail using the approximate analytical solution of the Boltzmann equation developed in ref. [18]: Very roughly, x F ∼ 20 for WIMP masses and annihilation cross sections of interest. This equation shows that g 1/2 * affects the final result more strongly than h does, which appears only logarithmically. The derivation of eq.(18) assumes that g 1/2 * and h are constant around the WIMP decoupling temperature T F = m χ /x F . This is not a very good approximation near the QCD deconfinement transition, where these functions change rapidly, as we saw in fig. 1. However, we can see directly from the Boltzmann equation that the most relevant temperature range is that around the decoupling temperature. At higher temperatures, Y χ is in any case close to its equilibrium value, which does not depend on g 1/2 * . At temperatures well below the decoupling temperature, i.e. for x x F , the right-hand side of the Boltzmann equation (5) Fig. 1, as function of the WIMP mass. The upper frame is for a constant σv , chosen such that our prediction for Ωχh 2 = 0.1193, while the lower frame is for a pure P −wave annihilation, with σv = 1.2 · 10 −24 cm 3 s −1 · T /mχ. These results are almost independent of the numerical size of the annihilation cross section.
suppression at x > x F is even stronger. Sharp features in g 1/2 * therefore give sharper features, with larger amplitudes, for pure P −wave annihilation than for S−wave annihilation.
We noticed earlier that the older treatment of ref. [19] overestimates g 1/2 * for some range of temperatures above T c . Eq. (18) indicates that this should lead to a smaller predicted relic density, which is confirmed by fig. 2. Since g 1/2 * is over-estimated for an extended range of temperatures, the effect on the relic density is about the same for S− and P −wave annihilation, amounting to about 5% near the peak of the ratio shown in fig. 2. The second, much lower peak near m χ = 1 TeV is probably due to lack of knowledge of the top mass at the time when ref. [19] was written.
The spike in g 1/2 * predicted by the micrOMEGAs treatment of the results of ref. [20] gives prominent spikes in the ratios shown in fig.2. These spikes are numerical artefacts that result from the smoothing procedure used in micrOMEGAs. As argued in ref. [20], a true δ function spike in g 1/2 * should not affect the numerical solution of the Boltzmann equation, which necessarily entails some discretization. The probability that the program then has to evaluate the right-hand side of the Boltzmann equation at the precise value of x where the δ function diverges is zero. We nevertheless show results including this spike since it results from the "standard treatment" encoded in micrOMEGAs. Outside the mass range affected by this spike, the results of ref. [20] predict a slightly too large relic density, consistent with our observation that it predicts smaller values of g 1/2 * and h than our treatment does. Note that for fixed g 1/2 * , reducing h will (slightly) increase x F , leading to an increase of the predicted relic density. A decrease of h therefore goes into the same direction as a decrease of g 1/2 * . However, the fact that the relative difference between our calculation and the prediction based on ref. [20] is almost the same in both frames of fig. 2 at large WIMP masses shows that the main effect still comes from the change of g 1/2 * . We saw in fig. 1 that the prediction for g 1/2 * from ref. [21] lies below our prediction, except for a very narrow range of temperatures around T c . As a result, for pure S−wave annihilation the prediction for the relic density based on the treatment of ref. [21] lies above our prediction for all WIMP masses larger than 2 GeV. We argued above that the relevant range of temperatures is (even) smaller for pure P −wave annihilation. This explains why the blue curve in the lower frame of fig. 2 goes slightly above 1 for m χ 25 GeV. Note also that the predictions using our treatment agrees with the prediction using ref. [21] to better than 1% for all WIMP masses, except in the range between 3 and 15 GeV where the difference reaches 9 (12)% for pure S − (P −)wave annihilation.
In order to put these results into perspective, it should be noted that the lattice QCD predictions for the energy and entropy densities listed in Table 1 of [23], on which our treatment is based, still have significant uncertainties, which decrease from about 14% at T = 130 MeV to about 3% at T = 400 MeV. The corresponding uncertainty in the relic density is up to 2.5% for 2 GeV ≤ m χ ≤ 20 GeV. Finally, we note that treating the charm quark as a free particle would increase the predicted relic density by about 2.2% for m χ 30 GeV.
The results shown in fig. 2 are almost independent of the assumed value of the WIMP annihilation cross section as long as the relic density comes out at least roughly correctly, although they evidently are sensitive to the functional dependence of σv on the temperature. Eq. (18) shows that thermodynamic effects enter primarily through g 1/2 * (x F ), and x F depends on the annihilation cross section only logarithmically.
On the other hand, the exact value of the annihilation cross section that reproduces the correct relic density, now (within standard ΛCDM cosmology) constrained to be Ω χ h 2 = 0.1193±0.0014 [43], does depend on g 1/2 * (T ∼ T F ) and, to a lesser extent, on h(T ∼ T F ). Precise knowledge of this cross section is important to constrain the free parameters of models of thermal WIMPs. Moreover, as we will see in more detail below, indirect DM searches now begin to probe annihilation cross sections close to the required value of σv , if the latter is (approximately) independent of the temperature.
In fig. 3 we show the required value of σv , assumed to be independent of the temperature, for a Majorana fermion, obtained from our refined calculation of g 1/2 * and h. This updates the results of ref. [24], which assumed Ω χ h 2 = 0.11 and used [21] to compute g 1/2 * and h. We see that for 10 TeV > m χ > 10 GeV the required value of σv is in fact closer to 2 · 10 −26 cm 3 s −1 than to the "canonical", often cited value of 3 · 10 −26 cm 3 s −1 . The near constancy of the required value over such a large range of WIMP masses results from an "accidental" cancellation of two effects. This can again be understood from the approximate analytical solution (18) of the Boltzmann equation. On the one hand, increasing m χ increases x F , which increases the relic density. Since x F depends only logarithmically on m χ , the freeze-out temperature T F = m χ /x F still increases as m χ is increased. As shown in fig. 1, this increases g 1/2 * (T F ), which in turn reduces the relic density. For m χ > 10 TeV, all SM particles are essentially fully relativistic at T F , i.e. g 1/2 * becomes independent of T , reaching its asymptotic value of 106.75. For these very large WIMP masses the required value of σv would thus increase logarithmically with m χ , in order to cancel the effect of the increase of x F . However, since by dimensional analysis and unitarity arguments [44] σv ∝ 1/m 2 χ , it is very difficult to find scenarios with sufficiently large WIMP mass for m χ > 10 TeV.
On the other hand, for WIMP masses below 10 GeV the rapid decrease of g 1/2 * (T F ) with decreasing T F shown in fig. 1 requires a rather rapid increase of σv , to a peak value of about 4.5 · 10 −26 cm 3 s −1 . Finally, for m χ < 0.35 GeV, g 1/2 * (T F ) becomes approximately constant again, with electrons, positrons, neutrinos, and photons contributing so that g 1/2 * 3.29. Since x F keeps decreasing with decreasing m χ , keeping the relic density constant requires that σv also decreases logarithmically with decreasing WIMP mass for these very light WIMPs.

EXPERIMENTAL CONSTRAINTS ON σv
In this section, we will compare experimental constraints from indirect WIMP searches and from analyses of the cosmic microwave background (CMB) with our prediction for σv shown in fig. 3. Since the CMB decoupled much later than WIMPs did, and hence also at a much lower temperature (∼ 0.3 eV rather than ∼ m χ /20), while WIMPs in galaxies now have an average kinetic energy of ∼ 10 −6 m χ , such a comparison is meaningful only if σv is largely independent of the temperature. If σv ∝ T , as in pure P −wave annihilation, or for even stronger T −dependence, the bounds on σv from the CMB and from indirect WIMP searches are still several orders of magnitude above the value required to obtain the correct relic density.
Currently the strongest and most robust upper bounds on σv from indirect WIMP searches come from searches for hard γ rays by the FermiLAT collaboration. Photons travel in straight lines through our galaxy, whereas charged particles get deflected by the galactic magnetic field. This not only introduces a sizable uncertainty, since our knowledge of this magnetic field is far from perfect; it also isotropizes the arrival directions of the WIMP annihilation products, making it impossible to focus on regions of space where the WIMP signal should be particularly strong. Another advantage of γ rays is that they are present in nearly all possible final states: hard photons can be emitted directly off final or intermediate state particles, can originate from the decay of neutral pions and other hadrons that result from the hadronization of qq final states or the decay of τ leptons, and can be produced from energetic electrons or positrons through "inverse Compton" upscattering of ambient photons.
The strongest WIMP signal is expected from near the center of our own galaxy. Unfortunately this region also hosts several backgrounds, both in form of point sources and in form of extended emission. It has been claimed that there is evidence for an additional component in the GeV γ flux from near the galactic center which can be explained through WIMP annihilation [45], but other interpretations of this additional component exist [46]. We also note that the FermiLAT collaboration itself has not published any analysis of their data on the galactic center.
In this paper we therefore focus on FermiLAT observations of nearby dwarf galaxies [13][14][15][16]. In contrast to big galaxies like our own, the mass density of dwarf galaxies should be dominated by dark matter even in the central region, yielding a much better signal-to-background ratio for indirect WIMP signals. No such signal has been seen. Our analysis is based on the very recent 6-year "Pass 8" analysis [16].
The results are shown in fig. 4. We see that the upper bound on σv is strongest if WIMPs predominantly annihilate into uū final states, but the bound for WIMP annihilation into bb is only slightly weaker. For the τ + τ − final state the upper bound on the cross section is similar for WIMP masses below 40 GeV, but is somewhat weaker for heavier WIMPs; hadronic final states have higher multiplicity, and hence higher γ flux per WIMP annihilation, for larger WIMP masses, whereas for the τ + τ − final state the photon multiplicity is essentially independent of the WIMP mass. These constraints exclude WIMPs with mass m χ ≤ 70 to 100 GeV annihilating into hadrons or τ leptons with temperature independent σv .  fig. 3 is compared with several observational upper bounds on σv , which is assumed to be independent of temperature. The cyan, green and magenta curves follow from the FermiLAT upper bound [16] on the γ flux from dwarf galaxies, for different dominant WIMP annihilation channel (uū, bb or τ + τ − ), whereas the red curve results from an upper bound on spectral distortions of the CMB, assuming WIMP annihilation into e + e − pairs.
As mentioned above, WIMP annihilation into e + e − can produce hard photons through upscattering ambient, e.g. visible, photons. The corresponding upper bound on σv (not shown) is worse than that for WIMP annihilation into τ + τ − by a factor of about two to three [16], excluding WIMPs with mass m χ ≤ 15 GeV annihilating into e + e − pairs for the value of σv shown in fig. 3.
WIMP annihilation can also affect the CMB. The strongest limits on σv originate from the fact that WIMP annihilation heats up the plasma in the "recombination" epoch when neutral atoms first formed [47], thereby delaying the decoupling of the CMB photons and distorting the pattern of CMB anisotropies. In fig. 4 we show the bound on the WIMP annihilation cross section into e + e − pairs that results from an analysis [17] of data from the WMAP and ACT collaborations. It excludes a thermal WIMP with m χ < 8 GeV.
PLANCK data will lead to considerable stronger constraints [43,48]. Unfortunately these papers only cite upper bounds on the product of the WIMP annihilation cross section and an efficiency factor f eff with which the energy of the WIMP annihilation products is absorbed in the thermal plasma. Using results from ref. [49], we estimate that the latest PLANCK data exclude WIMPs with m χ < ∼ 40 GeV annihilating into e + e − pairs with temperature independent cross section; see also the recent analysis [50].
Since the efficiency factor should be similar for the other final states considered in fig. 4, the CMB constraint should also be similar for all channels. The current CMB constraint is thus weaker than the bounds derived from the most recent FermiLAT data if WIMPs mostly annihilate into qq or τ + τ − final states, but is stronger for WIMPs annihilating predominantly into e + e − . However, one should keep in mind that the CMB constraint is less direct. It is conceivable that additional non-standard ingredients to the CMB fit -e.g., the presence of sterile neutrinos, a significant running of the spectral index of inflation, and/or a large contribution from tensor modes -can (partly) compensate the distortions caused by early WIMP annihilation, thereby weakening the constraint on σv . On the other hand, the constraint derived from the observation of dwarf galaxies depends on the assumed dark matter distribution [13]. In any case, it is encouraging that recent astrophysical and cosmological observations begin to probe relatively light thermal WIMPs with temperature independent annihilation cross section.

SUMMARY AND CONCLUSIONS
Using recent lattice QCD results with dynamical quarks for the equation of state, we have computed the energy and entropy densities of the SM with emphasis on temperatures around the deconfinement transition at T c = 154 MeV. These results are described by the functions g(T ) and h(T ). Of particular relevance for the calulcation of the relic density of thermal WIMPs is the quantity g 1/2 * defined in eq. (6), which also depends on the derivative of h with respect to the temperature. Our results for these functions can readily be embedded in the public codes computing the relic density [26][27][28].
Our predictions for the WIMP relic density differ from earlier treatments that relied on phenomenological models or pure glue lattice QCD calculations. These differences are most pronounced for WIMP masses between 3 and 15 GeV; they can reach about 9% for a constant σv , and up to 12% if σv ∝ T as in pure P −wave annihilation away from poles and thresholds. These differences are partly due to a spurious divergence in h(T ) leading to a divergence in g 1/2 * , which results in a sharp spike in the standard implementation. This illustrates the importance of ensuring that h(T ) is continuous everywhere, not just during the deconfinement transition but also during electroweak symmetry breaking at T 100 GeV, which we essentially ignored in our treatment.
It should be noted that the uncertainties in the recent lattice calculation we used still translates into an uncertainty of the predicted relic density of up to 2.5%, considerably larger than the observational uncertainty (at least in the framework of the standard ΛCDM cosmology). On the other hand, for m χ > 20 GeV our result for the relic density agrees to better than 1% with the best previous calculation [21].
We used our improved treatment of the thermodynamics of the very early universe to update the calculation of the WIMP annihilation cross section required to reproduce the observed dark matter relic density, assuming the thermal average σv to be independent of temperature. We found that this cross section is indeed nearly constant for 10 GeV < m χ < 10 TeV, thanks to a fortuitous cancellation between two competing effects; however, as also pointed out in ref. [24], the required value is closer to 2 · 10 −26 cm 3 s −1 than to the often-cited "canonical" value of 3 · 10 −26 cm 3 s −1 . On the other hand, for m χ < ∼ 3 GeV the required value of σv exceeds 4 · 10 −26 cm 3 s −1 . Producing a sufficiently large annihilation cross section for such light WIMPs typically requires the introduction of new light "mediators" between the WIMPs and SM particles. If the mediator mass is < ∼ m χ /10, the contribution of the mediator particles to the entropy and energy densities should be included, slightly modifying our results for the required annihilation cross section shown in fig. 3.
We also compared the required value of σv with upper bounds on this quantity that come from searches for energetic γ rays produced in WIMP annihilation as well as from CMB constraints. The most recent PLANCK data most likely give the strongest upper bounds, although a dedicated analysis for specific final states has not yet been performed. Moreover, the CMB bound assumes standard ΛCDM cosmology, and can thus more easily be evaded than the bounds from γ ray searches. All of these bounds can be evaded if the WIMPs predominantly annihilate into neutrinos [51].
It should be noted that both our calculation of the required WIMP annihilation cross section and our analysis of observational upper bounds on this quantity assumed that the thermal average σv is independent of temperature, or, equivalently, that the annihilation cross section is independent of the invariant center-of-mass energy. Theoretically this assumption is not well motivated. For non-relativistic WIMPs the annihilation cross section can usually (but not always [52]) be expanded in terms of the relative velocity v, σv = a + bv 2 + · · · . Even if the constant term a is not suppressed, i.e. if annihilation from an S−wave is allowed, one would generically expect b to be of the same order as a. This would reduce the required value of σv in today's universe by about 10%. If the constant term in σv is suppressed, the upper bounds on σv that follow from observations in today's universe do not constrain thermal WIMP models significantly.
However, even in this case it is important to calculate the relic density as accurately as possible. This leads to a constraint on the free parameters of the underlying WIMP model, which can hopefully one day be compared to direct measurements at colliders. Only then will we be able to say with some confidence that this WIMP model is indeed correct. This paper makes a small contribution to this ambitious long-term program.