Exact theory of freeze out

We show that the standard theory of thermal production and chemical decoupling of WIMPs is incomplete. The hypothesis that WIMPs are produced and decouple from a thermal bath implies that the rate equation the bath particles interacting with the WIMPs is an algebraic equation that constraints the actual WIMPs abundance to have a precise analytical form down to the temperature $x_\ast=m_\chi /T_\ast$. The point $x_\ast$, which coincides with the stationary point of the equation for the quantity $\Delta= Y-Y_0$, is where the maximum departure of the WIMPs abundance $Y$ from the thermal value $Y_0$ is reached. For each mass $m_\chi$ and total annihilation cross section $\langle \sigma_\text{ann}v_\text{r}\rangle$, the temperature $x_\ast$ and the actual WIMPs abundance $Y(x_\ast)$ are exactly known. This value provides the true initial condition for the usual differential equation that have to be integrated in the interval $x\ge x_\ast$. The matching of the two abundances at $x_\ast$ is continuous and differentiable. The dependence of the present relic abundance on the abundance at an intermediate temperature is an exact result. The exact theory suggests a new analytical approximation that furnishes the relic abundance accurate at the level of $1\%-2\%$ in the case of $S$-wave and $P$-wave scattering cross sections. We conclude the paper studying the evolution of the WIMPs chemical potential and the entropy production using methods of non equilibrium thermodynamics.


I. INTRODUCTION
Freeze out [1,2] refers to the idea that stable particles, once in thermal and chemical equilibrium in the hot and dense early universe, left a relic abundance because at a certain stage of the evolution, the expansion and cooling made their density so small that the annihilation reaction rates became frozen, see Refs. [3,4] for earlier reviews on the connection between cosmology and particle physics.
Today we know that most of the mass of visible Universe is constituted by some form of nonluminous dark matter as evidenced by the determination of the cosmological parameters by the WMAP [5] and Planck [6] satellites. A massive, neutral, stable particle that interact with the standard model particles with a strength similar to that of weak interaction (weakly interacting massive particles or WIMP), is the paradigmatic frozen relic and could explain the properties of dark matter.
Using natural units where ℏ = c = k B = 1, the well known [1,2] rate equation for the number density n of relics is where a is the Friedman-Robertson-Walker scale factor, a 3 the comoving volume, H = 1 a ( da dt ) the Hubble parameter, σ ann v r the total thermally averaged annihilation unitary rate and n 0 the equilibrium number density that is a known function once the statistics obeyed by WIMPs is specified.
The approximate analytical solution of (1), which is known as freeze out approximation, was given by Zeldovich, Okun and Pilkener already in Ref. [1]: at earlier times during the expansion n ∼ n 0 while at later times the creation term proportional to n 2 0 in (1) can be neglected and the equation is easily integrated. The instant of "freezing" [2], or "quenching" [1], separating the two stages was estimated by requiring t −1 where τ −1 = σ ann v r n 0 ≡ Γ is the characteristic time of interactions in equilibrium, and t −1 is the characteristic time of variation of N 0 = n 0 a 3 in equilibrium. As initial condition for the second stage they set approximately N (t f ) ∼ 2N 0 (t f ). The freeze out condition used by Lee and Weinberg, after rewriting the equation in terms of f = n/T 3 and x = T /m, is equivalent to the previous one [2]. Earlier calculations of the abundance for various relics and dark matter candidates, solving numerically the equation or using the approximate solution, are found in Refs. [1,2,[7][8][9][10][11][12][13][14][15][16].
In present literature it is typically adopted the freeze out approximation in the form given by Scherrer and Turner [17], which is similar to the earlier treatment of Steigman [3,22]. In order to discuss this method it is useful to write Eq. (1) in a different form [16][17][18][19]. If the expansion proceeds adiabatically, the entropy per comoving volume S = sa 3 , with entropy density s = (2π 2 /45)g s T 3 , is conserved. In the radiation dominated epoch H = (8/3)πGρ with M P = 1/ √ G = 1.22 × 10 19 GeV the Plank mass and ρ = (π 2 /30)T 4 g ρ the energy density. Using the ratio Y = N/S = n/s as a function of the x = m χ /T instead of the time, the equation becomes [16][17][18][19] dT ) accounts for the temperature dependence of the relativistic degrees of freedom g ρ and g s [18,[20][21][22]. The present relic abundance is given by Y (x 0 ) with T 0 = 2.725 K the temperature of the microwave background today. Knowing Y (x 0 ) and the entropy density today s 0 = 2890 cm −3 , the dark matter relic mass density is ρ DM = m χ Y (x 0 )s 0 and the ratio over the critical density ρ c = 3H 2 0 /(8πG) is where h = H 0 /(100 km s −1 Mpc −1 ) is the reduced Hubble constant and the numerical value is the experimental central value.
The method of Ref. [17] consists in studying the equation for the difference ∆ = Y − Y 0 and neglecting the derivative d∆/dx in the region x < x f where Y ∼ Y 0 . One finds ∆ ∼ − 1 2 x 2 C σannvr Y0 dY0 dx , while for x > x f one can neglect the term proportional to Y 2 0 and integrate easily. The freeze out temperature is found with the condition ∆(x f ) = δY 0 (x f ) with δ a constant numerical factor O(1) determined by fitting the numerical solution of (2) for various masses and cross sections. Values for δ such as √ 2 − 1 [17,23], 1.5 [18,24], 1/2 [3,25], ( √ 5 − 1)/2 [22] has been found to provide the relic abundance accurate at level of few percent.
On the other hand, in order to solve numerically Eq. (2) and find the asymptotic value that gives the relic abundance, it is necessary an initial condition, the initial abundance that is unknown. The numerical integration hence is done by setting Y = Y 0 at some x, for example x = 2, where the solution is very close to Y 0 . Strictly speaking this is not correct because at any x > 0 it is Y > Y 0 but the difference Y − Y 0 anyway is so small that the numerical result is correct for all practical purposes. From this point of view the existence of the intermediate temperature that determines the asymptotic value of the solution is an artefact of the analytical approximation [26].
In this work we show that, under the same assumptions of validity Eq. (1), the problem is actually defined by a system of two equations: the solution Y 1 of the new equation gives the abundance up to an intermediate temperature, that we identify with the freeze out temperature x f ; the differential equation (2) is valid only for We also suggest that for the analytical approximation of the relic abundance at level of 1% − 2%, it is better to consider as initial condition the temperature

II. KINETIC EQUATIONS AND SOLUTIONS
The derivation of Eq. (1) from the Boltzmann equation has been widely discussed in literature [14,[27][28][29]. Here we note that the integrated Boltzmann equation (1) has also the form of a kinetic chemical equation, that suggests to approach the problem from the point of view of chemical reactions. Let be a a simple 2 → 2 reaction between species X i with stoichiometric coefficients ν i . The rate of reaction, for example, of the specie X 1 , is given by where c i is the molar concentration if the volume is constant or the quantity matter if the volume changes during the reaction, and K (1,2→3,4) is the temperature dependent reaction constant. We now consider χ, a self-conjugated neutral particle that annihilate into particle-antiparticles pairs of the thermal bath constituted by N b species ψ i ,ψ i , i = 1, ..., N b of the standard model. We further assume that the only reactions that change the number of WIMPs is pair annihilation, 2 χ → ψ iψi , and pair-creation ψ iψi → 2χ. The analogous of the reaction constant is determined by the microscopic expression for the averaged rate per volume R a 3 , with The factor 1/(1 + δ 12 ) account for the fact that if 1 and 2 are the same specie and can react, then each reaction is counted two times; v r is the relative velocity of the colliding particles and P (v r ) its probability distribution [30]. We thus identify c i with na 3 and K with R a 3 in (4). Accounting for both creation and annihilation reactions for all species, we are led to the system of rate equations for bath particles, and for WIMPs. We used the notation σ χχ v r i and σ ψiψi v r for the annihilation and creation rates, respectively. In all generality this is a system of 2N b + 1 coupled differential equations.
In order to simplify the problem we make the standard assumptions: (i) the momentum distribution, and consequently the distribution of the relative velocity, is determined by the same equilibrium statistics of n 0 and is the same for all the particles; (ii) the number density of the bath particles is equal to the equilibrium number density during all the evolution n ψi = nψ i = n 0,ψi ; (iii) the chemical potentials and quantum statistical effects are negligible.
In this way we can replace each pair-creation rate σ ψiψi v r n ψi nψ i with σ ψiψi v r n 2 0,ψi . In chemical equilibrium, the total annihilation rate is balanced by the total creation rate, and the total number of particles, WIMPS and bath, (n 0,χ a 3 ) + 2 Defining the total annihilation unit rate σ ann v r = N b i=1 σ χχ v r i , after substitution of (9) and (10) in Eqs. (6), (7), the sum of the 2N b equations for the bath gives the equation while Eq. (8) becomes Eq. (1), Note the cancellation between the stoichiometric and the statistical factors [13]. Had we started with a particleantiparticle WIMP system with equal abundances, there would be another equation forχ that is equal to the one for χ, with stoichiometric and rate statistical factors equals to one. The doubling of all the terms when summing the two equations cancels, hence Eqs. (11) and (12) are obtained again. In the following we drop the subscript χ in the densities and work with Eqs. (11), (12) in the form (2) for the abundance Y . The problem is thus defined by a system of two equations. Given that the right-hand side of (11) and (12) is the same, then or Y = Y 0 all over the range, that is not possible for a massive particle, or the two equations are valid in different intervals and their solutions are matched at some intermediate temperature x f , the freeze out point. Equation (11) is easily solved for, once specified the statistics obeyed by the WIMP, Y 0 and its derivative are known functions. Defining the solution Y 1 of (11) and the difference ∆ = Y 1 − Y 0 can be cast in the form When the second term in the square root of Eq. (13) is small, one can expand and ∆ reduces to the approximate formula reported in the Introduction. The Steigman-Scherrer-Turner ansatz is in reality exact for Y 1 for all x with δ(x) a known function once σ ann v r and Y 0 are given. For x ≥ x f the abundance Y 2 is determined by the differential equation (2) with the initial condition At the matching point Y 2 (x f ) = Y 1 (x f ) and Eq. (16) gives To put the freeze out condition in a different form, we note that the difference ∆ grows with increasing x, reaches a maximum and then rapidly tends to Y 1 . At the maximum, d∆/dx = 0 that implies dY 1 /dx = dY 0 /dx, or, In other words, the relative variation of the equilibrium function Y 0 equals the relative variation of the difference In all generality, Eqs. (13), (14) together with Eqs. (16), (17) and the matching condition (19) furnish the exact theory of freeze out and relic density.
To see how this work in practice, we use the relativistic Maxwell-Boltzmann-Juttner statistics [18,19,30,31] to specify the equilibrium function, where K i are modified Bessel functions of the second kind and g χ are the spin degrees of freedom of the WIMP. For a numerical example we take m χ = 100 GeV, g χ = 2, σ ann v r = 10 −10 GeV −2 and g s = g ρ = g = 100, √ g * = √ g, neglecting temperature dependence of the degrees of freedom. Using the formulas (20) in Eq. (19) the freeze-out temperature is x f = 20.32 and δ(x f ) = 0.24. In the top panel of Figure 1, the dashed curve is the numerical solution of Eq. (2) with the false initial condition Y (2) = Y 0 (2). Superimposed is the piecewise function Y 1 (x ≤ x f ), blue curve, and Y 2 (x ≥ x f ), red curve, given by the numerical solution of Eq. (16). Note that δ(2) ∼ 10 −10 , which explains the perfect matching of the curves. From this point of view solving numerically the equation with the false initial condition is the best possible approximation to the exact theory. Now the usual analytical approximation for Y 2 at x > x f can be made more precise. Neglecting the term Y 2 0 in (2) and integrating with x f and Y 1 (x f ) as initial condition we have We can send to infinity x in (21) and x 0 in (3) because T 0 ≪ m χ . The relic abundance is then given by where we have taken a power law dependence on the temperature, σ ann v r = σ ann v r 0 x −n [17]. As can be seen in the bottom panel of Fig. 1, a zoom on the freeze out zone evidenced with a box in the top panel, Eq. (21) underestimates the exact solution by 20% in the freeze out zone but we find that (22) underestimates by 4% in the asymptotic region at large x.
As far as concerns the analytical approximation, one is not obliged to use x f but can choose a better point. The function ∆ posses another characteristic point, where and δ(x 2 ) = 1, Y 1 (x 2 ) = 2Y 0 (x 2 ). In the above numerical example, x 2 = 21.97, see the inlay in the top panel of Fig. 1 In the bottom panel we show the improvement obtained by setting the as initial conditions x 2 and Y 1 (x 2 ) in Eq. (21). Now we find, for example, Y 2 (10 8 )/Y x2 ∞ = 1.009. We have also considered the case of σ ann v r = σ ann v r 0 /x, P -wave scattering, with σ ann v r 0 = 10 −10 GeV −2 . The new points are x f,P = 18.19 and x 2,P = 19.11. The accuracy of the asymptotic value (22) using x f is 5%, but Y 2 (10 8 )/Y x2,P ∞ = 1.018 using x 2,P . We conclude that Eq. (22), with x 2 and Y 1 (x 2 ) in the place of x f and Y 1 (x f ), gives the relic abundance with an accuracy at the level of 1% − 2% Finally, it is worth to note that using Eq. (13), the freeze out condition δ(x 2 ) = 1 becomes Going back to the original condition of Ref. [1] discussed in the Introduction, we see that to have the freeze out point when N (t f ) = 2N (t f ) the condition is t −1 1 = 3τ −1 .

III. SUMMARY AND FINAL COMMENTS
In summary, we have shown that if bath particles maintain their statistical equilibrium density during the expansion and freeze out process, then the departure from the equilibrium of the WIMPs is fixed through Eq. (11). The solution is known and is valid only up to x f given by Eq. (19). The validity of the usual equation is restricted to x ≥ x f . The fact that the asymptotic value of the latter, the relic abundance, depends on an intermediate point and not on an initial condition at small x is thus clarified: the complete evolution is actually determined by two equations and the solutions Y 1 and Y 2 are matched at x f . The existence of x f is not an artefact of the freeze out approximation.
Another computational advantage is in the fact that performing the numerical integration starting at x f , the equation is less stiff: there are no singular points and the solution varies at most by 2 orders of magnitudes. We think that our findings may also help to improve numerical algorithms of public codes for relic abundance calculation [32][33][34][35].