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∗=mχ/T∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_*=m_\chi /T_*$$\end{document}. The point x∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_*$$\end{document}, which coincides with the stationary point of the equation for the quantity Δ=Y-Y0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta = Y-Y_0$$\end{document}, is where the maximum departure of the WIMPs abundance Y\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y$$\end{document} from the thermal value Y0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_0$$\end{document} is reached. For each mass mχ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_\chi $$\end{document} and total annihilation cross section ⟨σannvr⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \sigma _\text {ann}v_\text {r}\rangle $$\end{document}, the temperature x∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_*$$\end{document} and the actual WIMPs abundance Y(x∗)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y(x_*)$$\end{document} are exactly known. This value provides the true initial condition for the usual differential equation that have to be integrated in the interval x≥x∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x\ge x_*$$\end{document}. The matching of the two abundances at x∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_*$$\end{document} 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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S$$\end{document}-wave and P\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P$$\end{document}-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.


Introduction and motivations
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. Earlier reviews of the connection between cosmology and particle physics are found in Refs. [3,4]. a e-mail: mirco.cannoni@dfa.uhu.es Today we know that most of the mass of the Universe is constituted by some form of non-luminous dark matter as evidenced by the determination of the cosmological parameters from the data provided by the WMAP and Planck satellites [5,6]. A massive, neutral, stable particle that interacts with the standard model particles with a strength similar to that of weak interaction (weakly interacting massive particles or WIMPs) is the paradigmatic relic and could explain the properties of dark matter.
Using natural units where = c = k B = 1, the wellknown [1,2] rate equation for the number density n of relics is where a is the Friedmann-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 number density at zero chemical potential determined by the equilibrium statistics obeyed by WIMPs. It is useful to write Eq. (1) in a different form [7][8][9][10]. If the expansion proceeds adiabatically, the entropy per comoving volume S = sa 3 is conserved, with entropy density s = (2π 2 /45)g s T 3 . In the radiation dominated epoch 19 GeV the Planck mass and ρ = (π 2 /30)g ρ T 4 the energy density. Using the abundance Y = N /S = n/s as a function of x = m χ /T instead of the time, with the change of variable d/dt = H xd/dx, the equation becomes [7][8][9][10] where √ g * = g s / √ g ρ (1 + T /3 d(ln g s )/dT ) accounts for the temperature dependence of the relativistic degrees of freedom g ρ and g s [9,[11][12][13].
The present relic abundance is given by Y (x 0 ), which is the constant asymptotic value of the solution of Eq. (2), with T 0 = 2.725 K, the temperature of the microwave background today. The dark matter relic mass density is ρ DM = m χ Y (x 0 )s 0 , with s 0 = 2890 cm −3 the entropy density today. The ratio over the critical density ρ c = 3H 2 0 /(8π G) is Ω DM h 2 = 2.75 × 10 8 m χ /GeV Y (x 0 ) 0.11, where h = H 0 /(100 km s −1 Mpc −1 ) is the reduced Hubble constant. The numerical value of the ratio is the present experimental central value.
To a very good approximation, quantum statistical effects can be neglected [9], thus the thermal number density for non degenerate relativistic gas is given by the Maxwell-Boltzmann-Juttner statistics [14] with zero chemical potential μ i = 0, Here g are the spin degrees of freedom of the particle and K i are modified Bessel functions of the second kind. In the relativistic limit x = 0 become n m=0 An approximate estimation of the relic density can be obtained with the approximate analytical solution of Eq. (1), which is known as freeze-out approximation, already formulated by Zeldovich et al. [1]. The evolution can be divided into two stages: in the first stage, at earlier times during the expansion, WIMPs go slowly out of chemical equilibrium, thus n ∼ n 0 . The relative variation of the density because of the expansion has characteristic time t −1 1 = 1/(n 0 a 3 )d(n 0 a 3 )/dt, while the characteristic time of the interactions in equilibrium is τ −1 = σ ann v r n 0 . The instant t f of "freezing" [2], or "quenching" [1], separating the two stages is estimated by requiring t −1 1 = τ −1 . In modern notation, we can write We have inserted a minus sign because the derivative of Y 0 is always negative. In the second stage, at later times, the creation term proportional to Y 2 0 in (2) can be neglected and the Equation (5) is used in Refs. [1,4,15] as the approximate initial condition for the integration. The origin of the factor 2 is not explained. In this respect, we note that the freezeout condition used by Lee and Weinberg [2] and successive earlier works [16][17][18][19][20][21] is equivalent to Eq. (4), but in the initial condition for the integration there is no factor 2 as in (5). We will discuss this point in Sect. 4.
The criterion (4) takes a more familiar form in the nonrelativistic limit. From Eq. (3) it is easy to see that Using C/x 2 = s/(H x) and Γ = σ ann v r Y 0 s, Eq. (4) becomes which is the more precise version (see the Erratum of [8]) of the popular criterion "the annihilation rate per particle Γ equals the expansion rate H ". The exact solution of the non-linear Riccati equation (2), on the other hand, in general can only be obtained with numerical methods. In order to do that, an initial condition is necessary. The differential equation is of first order, thus its solution is determined completely by the initial condition. The true initial condition is Y (0) = Y 0 (0), which anyway is only an asymptotic condition because the equation has a singularity in x = 0. The numerical integration hence is done by setting Y = Y 0 at some small x, for example x = 2 [10], where the solution is expected to be very close to Y 0 . It must be stressed that this is a false initial condition because Y > Y 0 at any x > 0. The difference Y − Y 0 at small x is anyway so small that the numerical result is correct for all practical purposes.
If the differential equation gives the actual abundance in all the range (0, ∞), the freeze-out temperature x f does not exist, is an artificial arbitrary point introduced by the freezeout approximation [22]. Also the "exact" solution is in reality an approximation because uses a false initial condition.
The problem of the initial condition is important for the relic abundance calculation because it determines the precision with which the asymptotic value is obtained. For example, public codes use very different approaches. In Dark-SUSY [23] the numerical integration is started with Y (x = 2) = Y 0 (x = 2) with a step-size changing adaptive implicit method to avoid stiffness problems. In MicrOMEGAs [24] the initial point is ; then the analytical approximation neglecting Y 0 is employed. SuperIso [25] uses the freeze-out approximation as given of Ref. [9]. In MadDM [26] the asymptotic value is calculated integrating the equation from x = 30 up to a large vale, x = 1000. The process is repeated backward varying the initial is less than the required precision. It is clear that these codes, leaving aside the uncertainties intrinsic in the numerical methods, for a given WIMP mass and total annihilation cross section inevitably will differ in the value of the relic abundance, even if within a few percent in normal situations.
The problem that we pose and answer in this paper hence is the following. Suppose that for each WIMP mass and total annihilation cross section the actual abundance is exactly known in one intermediate point x * . Then we can use it as true initial condition for the differential equation, and the asymptotic value at large x calculated in this way would be the true value. But a conceptual problem arises: it would tell us that the relic abundance really depends on an intermediate temperature, in a way similar to the freeze-out approximation, and not on the asymptotic value at x = 0 as it should be if the first order differential equation described the evolution over the whole range. In other words, in that case the usual differential equation would have a boundary value at x * , the evolution before x * would be given by something else.
In the next section we show that such a point exists.

A special temperature
In present literature it is typically adopted the freeze-out approximation in the form given by Scherrer and Turner [8], which is similar to the earlier treatment of Steigman [3,13].
Calling λ(x) = C σ ann v r /x 2 , the method consists in studying the differential equation for the difference Δ = Y − Y 0 , the "distance" from equilibrium, The matching point of the two approximate solutions in the two regimes is found with the ansatz Δ(x f ) = cY 0 (x f ), being c being an arbitrary constant numerical factor O(1) determined by fitting the numerical solution of Eq. (2) for various masses and cross sections. Values for c such as √ 2 − 1 [8,27], 1/2 [3,28], ( √ 5 − 1)/2 [13], 1.5 [9,29], has been found to provide the relic abundance accurate at the level of a few percent. The value of x f is determined neglecting dΔ/dx in (8) and substituting Δ(x f ) = cY 0 (x f ). It depends logarithmically on c. The equation is numerically stiff, thus is not surprising that different authors find different values for numerical coefficients c.
In reality Eq. (8) tells something much more precise. Equation (8) possesses non-trivial stationary points where dΔ/dx = 0. The difference Y − Y 0 is initially zero, then grows (always positive) and tends to Y asymptotically because Y 0 is rapidly decreasing. Hence the extremal point is a maximum. We indicate the stationary point with x * to distinguish it from the arbitrary point x f of the freeze-out approximation. At x = x * Eq. (8) is a quadratic algebraic equation with constant coefficients, Disregarding the non-physical negative solution, and indicating with a subscript * all the quantities evaluated at x * , the physical solution is and the abundance at x * then takes the form If we now define then we have At x = x * , the ansatz is an exact result with c = c * . Now c is not an arbitrary numerical value, but it is determined in each case by the WIMP mass, the total annihilation cross section and Y 0 (x). 1 At x * the actual abundance Y and c have a well-defined analytical form as seen in Eqs. (11) and (12). We call the functions defined by those formulas Y 1 (x) and δ(x), The point x * where dΔ/dx = 0 is evidently given by dY 1 /dx| x=x * = dY 0 /dx| x=x * , or, using Eqs. (15), (16), The root of Eq. (17) gives x * . Using Eq. (11) also the actual abundance at the point x * is known. The differential equation As far as the calculation of the relic density is concerned, the problem where to start the integration is solved: for each mass and total annihilation cross section, Eqs. (17) and (15) furnish the true initial condition, and what happens before x * does not matter. There is no reason to start the integration before (or after) x * with a false initial condition. We remark that there is no approximation in this procedure. But if the evolution before x * does not affect the abundance after x * , it means that the differential evolution in temperature really starts at x * , while at higher temperatures it should be fixed by some constraint.
Let us look at the properties of the function Y 1 (x). The function Y 1 (x) is not solution of the differential equation, as it is easily seen by substitution, and it has the correct behavior to describe the evolution before x * . In fact, it tends to Y 0 for x → 0 asymptotically. At small x we may expand the square root and find which is the usual freeze-out approximation [8,30]. More important, given that is continuous and differentiable in x * as it must be for the exact evolution. In the freeze-out approximation, on the contrary, x f is an elbow point. Then the natural questions: is Y 1 (x) the exact actual abundance at x ≤ x * ? Where does it come from?

Exact theory
Usually Eq. (1) is derived from the Boltzmann equation [9,[30][31][32][33][34]. The integrated Boltzmann equation (1) has the form of a kinetic chemical equation, which suggests to approach the problem from the point of view rate equations that describe non-equilibrium chemical reactions. This thermodynamic approach allows one to arrive at the results and answer the previous questions more easily.

Generalities on rate equations
Let ν 1 X 1 +ν 2 X 2 ν 3 X 3 +ν 4 X 4 be an elementary stoichiometric 2 → 2 reaction between species X i with stoichiometric coefficients ν i that can proceed at the same time in both directions. We follow the convention that the stoichiometric coefficients ν i are negative for reactants X 1,2 and positive for products X 3,4 . It is convenient to consider the concentration given by the number density n i = N i /V , with N i the numbers of particles instead of the numbers of moles in a volume V .
The forward (left to right) and backward (right to left) rates are where k f , k b are the related temperature dependent rate constants. In chemical equilibrium, from the detailed balance principle, the rates are equal, R f = R b , and the concentrations stay constant at their equilibrium value that defines the If concentration changes because both the number of particles, due to reactions, and the volume, due to expansion, change in time, then the rate equation for the specie X i is given by Consider χ , a self-conjugated neutral particle that annihilates into particle-antiparticles pairs of the thermal bath constituted by N b species ψ i , i = 1, . . . , N b of the standard model. The ψ particles are all the particles that interact with the WIMP in a given model, including massless particles like the photons. In this case clearly ψ =ψ. We further assume that the only reactions that change the number of WIMPs is pair annihilation and pair creation, In the case of 2 → 2 collisions in a gas of elementary particles, the temperature dependent rate is given by The factor 1/(1 + δ 12 ) accounts for the fact that if 1 and 2 are the same species and can react, then each reaction is counted twice. Here v r is the relative velocity of the colliding particles and P(v r ) its probability distribution. In the case of relativistic particles obeying Maxwell-Boltzmann-Juttner statistics, the probability distribution of the relativistic relative velocity was discussed in Ref. [35].

Standard assumption for thermal relics
In the early Universe the temperature was so high that all fields behaved as thermal radiation with zero chemical potential, i.e. the number density of particles depends only on the temperature as given by Eq. (3). The reactions 2 χ ψ iψi at the initial time were in chemical equilibrium, While the number density of the bath particles is much larger than that of the WIMPs and maintains its thermal form during the freeze-out process in reason of the fast interactions between them, WIMPs have to go out of thermal equilibrium, n χ = n 0,χ . Altogether also chemical equilibrium would be maintained, with densities changing only because of the expansion/temperature variation and there would be no relic abundance. Note that by assumption (25), equilibrium statistics with chemical potential maintained to zero means that the bath particles behave like thermal photons of the black body radiation, i.e. their density is completely fixed by the temperature and does not depend on the amount of material of the cavity, which in our case is the amount of WIMPs in the comoving volume. As a consequence, in the process of thermal WIMP production and decoupling the total number of particles is not conserved.

Rate equation for WIMPs
Using the above formalism with V = a 3 , and the notation σ χχ v r i and σ ψ iψi v r for the annihilation and creation rates, respectively, we can write the rate equation for WIMPs, The assumption (25), together with Eqs. (24), implies that the creation rates are determined by thermal densities: Doing these substitutions in Eqs. (26) We can now define the total annihilation cross section, and sum over the channels. Equation (28) finally reads which is the usual differential equation for WIMPs 2 derived under the usual assumptions. 2 Note the cancelation between the stoichiometric and the statistical factors [20]. Had we started with a particle-antiparticle WIMP system with equal abundances, there would be another equation forχ that is equal to the one for χ, with stoichiometric and rate statistical factors equal to 1. The doubling of all the terms when summing the two equations cancels, hence Eqs. (36) and (30) are obtained again.

Equation for bath particles
Let us pass now to the equations for the bath particles, which read The rate equation forψ i coincides with Eq. (31) settingψ i in the time derivative, thus we do not write it explicitly. Given that the equations for n ψ i and nψ i are similar we sum Eq. (31) over the N b channels and multiply by 2. Employing the assumption (25), Using now (27) and (29) as we did for the WIMPs, Eq. (32) results in Let us pause and look at this equation. Already at this stage it is clear that this is an algebraic equation, not a differential equation. In fact, the only unknown is the actual density n χ of the WIMPs, which appears only on the right-hand side. The derivative on the left-hand side is the sum of the time/temperature derivative n 0,ψ i for each of the N b species. The left-hand side of (33) is hence a known function as much as n 0,χ is.
If the bath particles had a chemical potential, hence their density n ψ,i would evolve in time as dictated by the rate equation in a similar way to n χ . Conservation of the total number of particles would hold and we would have 2 N b i=1 d(n ψ,i a 3 )/dt = −d(n χ a 3 )/dt. Equation (31) would give the same equation (26). We would have only one equation, but the creation term would be σ ψ iψi v r n ψ i nψ i and could not be expressed in terms of the annihilation cross sections and of the thermal equilibrium density of the WIMPs.
Anyway, by hypothesis, the bath particles are locked to the thermal state with zero chemical potential, i.e. the number density of the bath particles depends only on the temperature and is not affected by the reactions with the WIMPs. The total number of particles is not conserved: the time/temperature derivative 2 N b i=1 d(n 0,ψ i a 3 )/dt cannot be equal to −d(n χ a 3 )/dt.
Can we further simplify the left-hand side of Eq. (33)? Note that the expansion changes the density n 0 but the product n 0 a 3 is not affected. At the same time the thermal density n 0 depends only on temperature and is not changed by the reactions. The quantity then remains constant during the freeze-out process. Hence which substituted in Eq. (33) finally gives In terms of the abundance Y = n/s and dropping the subscript χ , Eqs. (36) become which solved for Y directly gives the function Y 1 of Eq. (15). Given that the algebraic equation (37) and the differential equation (2) must hold at the same time and the right-hand side is the same, besides the trivial solution Y 1 , the only other possibility is that the function Y 1 gives the actual abundance in the interval (0, x * ], while the function Y 2 gives the actual abundance in the interval [x * , ∞). The solutions can coincide only at point x * and the matching is continuous and differentiable as we have already seen.
The functional form of the abundance discussed in Sect. 2, Eq. (11), hence, is the true non-equilibrium WIMPs abundance at x ≤ x * .

Numerical example
In order to show how this works in practice, we take m χ = 100 GeV, g χ = 2, σ ann v r = 10 −10 GeV −2 , g s = g ρ = g = 100, and √ g * = √ g, neglecting the temperature dependence of the degrees of freedom. From Eq. (17) we find that the freeze-out temperature is x * = 20.32. The inlay in the top panel of Fig. 1 shows the behavior of Δ and the position of its maximum. In the top panel of Fig. 1 we show the piecewise function built with Y 1 (x ≤ x * ), Eq. (15), solid turquoise curve, and Y 2 (x ≥ x * ), dashed-red curve, obtained by numerical integration of Eq. (18) with the initial condition (19). The standard picture is obtained.
The actual abundance Y 1 in the early stage at small x closely tracks the equilibrium function Y 0 because the term −x 2 /(C σ ann v r )dY 0 /dx is very small. Only with x approaching x * it becomes relevant; numerically we find that δ(x * ) = 0.24 but, for example, δ(2) ∼ 10 −10 . The smallness of δ(2) explains why if one integrates the differential equation with the false initial condition Y (2) = Y 0 (2), once having obviated the stiffness problems, one obtains a curve that is practically indistinguishable from that in Fig. 1. The numerical advantages of the exact theory should be clear: integrating the differential equation from x * , the solution varies smoothly less that two orders of magnitudes, instead of eight, before reaching the asymptotic constant value: there is no stiffness problem. The asymptotic value obtained in this way is the true value.

Discussion
If one writes down the Boltzmann equations for the bath particles and simplify them under the same assumptions used in deriving the equation for the WIMPs, one would arrive at the same results. The crucial point is not in the method, but in using the hypothesis of the thermal bath in all its consequences.
In a normal chemical process ν 1 X 1 + ν 2 X 2 ν 3 X 3 + ν 4 X 4 with atoms, molecules, nuclei or elementary particles where the total number of particles is conserved, the stoichiometry of the reaction implies that The four differential equations for the concentration of each specie are thus not independent. Equation (38) suggests that actually only one quantity is necessary to describe the changes in the number of moles/particles. This state variable in modern thermodynamics is called the extent of reaction [36]; it is defined by At the beginning of the reaction ξ(0) = 0. If the initial abundances at t = 0 are N 0,i , then at each time during the reaction progress the abundances of all species are given by where ξ(t) determined by the rate equation [36] 1 V with the initial condition ξ(0) = 0. Note that in Eq. (40) the rates are expressed in terms of ξ and the initial abundances appear explicitly in the equation.
If we remove the assumption of the thermal bath n ψ = nψ = n 0,ψ , the equations for WIMPs and for the bath particles would conform to the usual scheme. As we have already discussed above, we would have only one equation, but that equation would not have the standard form. In the creation rates the densities n ψ would appear instead of n 0,ψ , and we should know the initial abundance of both WIMPs and ψ i particles.
The presence of dY 0 /dx in Eq. (37) should not hide the fact that (37) is just another form of the algebraic equation (33), which accounts for the properties of the thermal bath.
Assumption (25) is basic for the theory of thermally produced WIMPs. It not only constrains the creation rate, but also the extent of reaction. Comparing Eqs. (26) and (31) with Eqs. (39) and (41), we see that what appears in the left hand sides of the former is the time derivative of the extent of reaction, while Eq. (40) tells us that Δ, up to the stoichiometric coefficient, is the same as the integrated extent of reaction. Imagine we can stop the expansion, that is, we fix the temperature, and observe the evolution in time of the creation reaction ψψ → 2χ as if it were an ordinary reaction. Let the initial amount of both bath particles and WIMPs be N 0 . The maximum extent to which this reaction can proceed corresponds to the total consumption of the ψ, ξ max = (0 − N 0 )/(−1) = N 0 where we used our convention of negative stoichiometric coefficients for the reactants. The amount of WIMPs at the maximum hence would be N 0 + N 0 = 2N 0 . This heuristically explains the Zeldovich approximate ansatz at the freeze-out point.
The resulting physical picture is clear. During the expansion WIMPs are produced by the thermal bath up to the temperature T * where the extent of reaction for production reaches the maximum value. Both T * and the maximum depends on the WIMP mass and on the strength of interactions between the WIMPs and the bath in a known constrained way. Once having reached this maximum, the production becomes negligible, the annihilation rate reduces the abundance that relaxes asymptotically to a constant value because the same annihilation rate diminishes with further expansion/cooling.
From the mathematical point of view the whole evolution in time/temperature of the WIMPs abundance is given by a differential-algebraic system. The differential equation (2) has the form dY/dx = f (Y, Y 0 , x), while the algebraic equation (37), which can be written as has the form g(Y, Y 0 , dY 0 /dx, x) = 0 and acts as a constraint in the interval x ≤ x * . The initial value problem Y (x * ) = Y 1 (x * ) is a consistent initial value problem for the differential equation given the physical assumptions and constraints.
In fact, we have seen that the production and decoupling of WIMPs depend both on dY 0 /dx and on Y 0 . This is taken into account by the algebraic equation g(Y, Y 0 , dY 0 /dx, x) = 0. The fact that the evolution of the abundance is differential only starting at x * is due to the fact that the function f (Y, Y 0 , x) does not depend on dY 0 /dx.

Approximation for the relic abundance: revising the Zeldovich criterion
The fact that the abundance at x ≤ x * is an already fixed function Y 1 (x) allows one to find a very accurate approximate formula for the relic abundance.
Let us see what happens if we apply the freeze-out approximation with x f = x * . Neglecting the term Y 2 0 in (18) and integrating with the initial condition (19) we have As can be seen in the bottom panel of Fig. 1, referring to Eq. (43), the green dotted-dashed line underestimates the exact solution Y 2 (x), red dashed line, by 20 % in the freeze-out zone. In fact, near x * , Y 0 is not much smaller than the actual abundance Y 2 (x) and the approximation is not so good. As far as concerns the analytical approximation, one is not obliged to use x * , but one can choose a better point. The function Δ posses another characteristic point, where it is intersected by the equilibrium function Y 0 , Δ(x 2 ) = Y 0 (x 2 ); see the inlay in the top panel and the bottom panel of Fig. 1. At x = x 2 , hence Y 1 (x 2 ) = 2Y 0 (x 2 ) and δ(x 2 ) = 1. With the same numbers as before, requiring δ(x 2 ) = 1 we find the numerical value x 2 = 21.97.
As it evident in the bottom panel, in the interval between x * and x 2 , Y 1 and Y 2 are still very close, while only after x 2 the true abundance Y 2 starts to deviate significantly from Y 1 . This behavior suggests one to use x 2 and as initial conditions in Eq. (43). In the bottom panel we clearly see the improvement: now the black solid line underestimates the red dashed line by 4 % around x 2 . At larger x Y 2 is practically superimposed, thus asymptotically the approximation will be even better. In fact, the relic abundance is hence given by where we have taken a power law dependence on the temperature, σ ann v r = σ ann v r 0 x −n [8], and sent x to infinity in (43). We find, for example, Y 2 (10 8 )/Y fo 2 (∞) = 1.01. In the case of P-wave scattering, σ ann v r = σ ann v r 0 /x, 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 asymptotic ratio using x 2,P is Y 2 (10 8 )/Y fo 2 (∞) = 1.02. Equation (45) gives the relic abundance with an accuracy at the level of 1-2 %, which is extremely good. We emphasize that the function Y 1 (x) extended up to x 2 , together with the function Y fo 2 (x) starting at x 2 , furnishes a semi-exact analytical description of the whole evolution of the actual WIMPs abundance.
Note that from Eq. (16), the condition δ(x 2 ) = 1 is equivalent to In the non-relativistic limit as in Eqs. (4), (6), (7) it becomes or, in terms of characteristic times, we have t −1 1 = 3τ −1 . We can establish now the connection with the Zeldovich criterion [1,4,15]. In Ref. [15] the freeze-out temperature was derived with the criterion The origin of this criterion can be understood reasoning as in the Introduction. The characteristic time of the expansion t −1 1 = (1/n 0 )dn 0 /dt is evaluated using the nonrelativistic density n 0 = cost.(mT ) 3/2 exp(−m/T ) and the relations between time and temperature in the radiation era, T ∝ t −1/2 and dT /dt = (1/2)T /t. Changing variable from t to T , and taking the derivative of n 0 neglecting the variation of the term T 3/2 , we have t 1 2t T /m.
The characteristic time of variation of the abundance due to interactions can be estimated from the kinetic equation written as dY/dt = − σ ann v r s(Y 2 − Y 2 0 ). Near equilibrium, in the right-hand side we can approximate s( , giving a characteristic time of variation τ = 1/(2 σ ann v r n 0 ). By requiring t 1 = τ we obtain Eq. (48).
The criterion thus gives very good approximation for the asymptotic value when compared with the numerical integration, as stated in Ref. [15], just a little bit worse than our proposed criterion. This is clear from the bottom panel of Fig. 1. Selecting a point between x * and x 2 , gives a freezeout approximation curve that lays between the dot-dashed and the full-black line: the point x V DZ is near x 2 , thus the freeze approximation curve starting from there will approximate well the exact numerical solution given by the dashed line.

Chemical potential of WIMPs and entropy production
We have seen that at the initial instant all chemical potentials are zero and chemical equilibrium holds. With the expansion/cooling of the universe and the bath particles with μ ψ i = 0, the WIMPs start to develop a nonzero chemical potential.
A non-equilibrium chemical transformation is characterized by the state variable called affinity, which is defined as [36] where μ i are the chemical potentials of the species. The affinity is zero in chemical equilibrium, A eq = 0, and is related to the rates by [36] A nonzero affinity drives the process out of chemical equilibrium. From Eq. (51), with the convention that stoichiometric coefficients are negative for reactants, we have Since the explicit form of the rates is known, using Eq. (52) we find and from Eqs. (53) and (55) the chemical potential is Here Y is the piecewise abundance (21). At small temperatures, large x, when the rates are virtually zero and the WIMPs abundance becomes constant, μ χ cannot grow arbitrarily and must tend to an asymptotic value. The chemical potential being the free energy per particle that can be released or gained when particles are produced or annihilated, the asymptotic value is the mass of the WIMP: In fact, at large x, ln(Y 2 /Y 0 ) ∼ ln(e x ) ∼ x, thus from Eq. (56) we have μ χ ∼ m. Instead, ln(Y 1 /Y 0 ) = ln(1 + δ(x)) ∼ ln(e x/2 ) ∼ x/2 thus the chemical potential would tend to m χ /2. This is consistent with the fact Y 1 gives the abundance only up to x f while for x > x f the evolution is given by Y 2 . The behavior is clearly seen in Fig. 2. The rate equation for the extent of reaction can also be written in terms of the affinity as [36] 1 Note that (1/V )dξ/dt = −1/2(1/a 3 )d(na 3 )/dt. Then, from Eq. (55), it follows that 1 − exp(−A/T ) = 1 − Y 2 0 /Y 2 = 1 − n 2 0 /n 2 and R f = ( σ ann v r /2)n 2 , thus Eq. (58) gives the standard equation.
The freeze-out is a nonequilibrium irreversible process that increases the entropy of the Universe. We can easily verify that anyway the produced entropy is negligible compared with the entropy of the radiation. The entropy production in  Fig. 1. The abundance Y is the piecewise function (21). The solid line corresponds to Y 1 and the dashed line to Y 2 as in Fig. 1 chemical reactions is determined by affinity and the extent of reaction through [36] 1 V This quantity is always positive as required by the second law of thermodynamics. By changing variable from time to x we find where S r is the entropy produced in the reaction and S = sa 3 is cosmic entropy density. It is easily seen that the scale of S r relative to S is set by the magnitude of σ ann v r . In our numerical example hence (dS r /dx)/S 10 −10 : the entropy produced in the production and decoupling of WIMPs is thus completely negligible compared with S, which thus remains constant during the process.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited. Funded by SCOAP 3 / License Version CC BY 4.0.