{'fraction_non_alphanumeric': 0.044558860920416984, 'fraction_numerical': 0.01288244766505636, 'mean_word_length': 4.053211991434689, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical mechanics and the description of the early universe I.: Foundations for a slightly non-extensive cosmology

We analyze how the thermal history of the universe is influenced by the statistical description, assuming a deviation from the usual Bose-Einstein, Fermi-Dirac and Boltzmann-Gibbs distribution functions. These deviations represent the possible appearance of non-extensive effects related with the existence of long range forces, memory effects, or evolution in fractal or multi-fractal space. In the early universe, it is usually assumed that the distribution functions are the standard ones. Then, considering the evolution in a larger theoretical framework will allow to test this assumption and to place limits to the range of its validity. The corrections obtained will change with temperature, and consequently, the bounds on the possible amount of non-extensivity will also change with time. We generalize results which can be used in other contexts as well, as the Boltzmann equation and the Saha law, and provide an estimate on how known cosmological bounds on the masses of neutrinos are modified by a change in the statistics. We particularly analyze here the recombination epoch, making explicit use of the chemical potentials involved in order to attain the necessary corrections. All these results constitute the basic tools needed for placing bounds on the amount of non-extensivity that could be present at different eras and will be later used to study primordial nucleosynthesis.

1 Introduction 2, then, does not pretend to show new results. Section 3 discusses the values of the plasma parameter along cosmic history, and assess the a priori possibility that the universe as a whole can be considered as a non-extensive system. Section 4 and 5 are the basis for all further development; they state the thermodynamical picture and the cosmologically conserved quantities in this new setting. Then, some new applications follow: Section 6 is concerned with the decoupling of relativistic and non-relativistic species; Section 7, with the process of recombination; and Section 8, with the new Boltzmann equation and the study of the process of freezing. Then, a brief Section 9 states the predictions for the current values of some observables of cosmological importance; and Section 10 studies hot and cold relics in the new framework and generalizes some of the well known bounds on neutrino masses and on other species. Finally, in Section 11, it is briefly accounted for the correction that non-extensivity would introduce to the value of matter-radiation time and temperature. We end by summarizing and obtaining some conclusions of a general nature. Paper II in this series will be an analytical analysis of the nucleosynthesis process, going all the way up to the formation of Deuterium and Helium 4. We shall be specially concerned with the analysis of the nuclear principle of detailed balance and the influence it has in this new setting.

Basics of non-extensive statistics
Among physicist and astronomers, there is a non-stated consensus on that the Boltzmann-Gibbs (BG) statistical mechanics is always applicable. However, in analogy with Newtonian dynamics and Special Relativity, it could a priori justified to consider it just as part of a bigger framework, an enlarged statistical description where extensive as well as non-extensive phenomena could be taken into account [14]. An interesting generalization of the BG entropy form has been recently proposed by Tsallis [15] (for recent reviews see [16,17], for a full bibliography see [18]). This new entropy, that we discuss below, possesses the usual properties of positivity, concavity and irreversibility, and generalizes the additivity in a non-extensive way. Examples in which BG statistics seems to present serious problems are systems for which there are long-range forces, and/or which present memory effects, and/or are subject to evolution in a non-Euclidean space [19].
The generalized form for the entropy is where k a positive constant, W the total number of physical states accessible to the system, and the set of probabilities p i satisfies Eq. (1) recovers the usual (BG) form for the entropy in the limit q → 1, i.e.
The entropic index q, intimately related to, and determined by, microscopic dynamics, characterizes the system under consideration. This reflects itself in the pseudo-additivity law for the entropy S q , where A y B are two independent systems in the sense that the joint probabilities (those corresponding to the system A + B) are such that p ij (A + B) = p i (A)p j (B).

Expectation values
We now introduce the following non-normalized expectation value, such that A 1 corresponds to the standard mean value for the observable A. If the system is a quantum one, its description is given in terms of the density operator, ρ, with eigenvalues {p i }. Then, the generalized entropy is and the expectation value becomes A q ≡ T rρ q A.
Eq. (6) can be recast as S q = −k ln q ρ q .
If the system is classic, and the relevant variables are continuum ones, we can describe it by a probability distribution p( r), where r is a dimensionless variable, say, in a many-body phase space. Then, the generalized entropy will be and the expectation value, We shall present the formalism in the case in which the system is described by a set of microscopic probabilities, W .

Canonical Ensemble
The first non-trivial physical situation is that of a system in thermal contact with a thermostat at a temperature T . We must then extremize the entropy taking into account several constraints. The first one is just the definition of probability, For the rest, a detailed discussion is in order.

Internal energy
There is some freedom concerning the choice of the constraint imposing the relationship among the different energy levels with the total internal energy. The first choice was introduced by Tsallis [15]: to Eq. (11), it is added the following constraint W i=1 p i ǫ i = U (1) .
Here, an index (1) refers to the first choice, and {ǫ i } are the Hamiltonian eigenvalues of the system. Using standard techniques, it can be seen that the {p i }-values that extremize S q with the imposed constraints are It must be said that β ⋆ is not the Lagrange multiplier associated with the constraint upon the internal energy. This expression recovers the usual one (p i ∝ e −βǫ i ) in the limit q → 1 and it depends on the microscopic energies as a power law instead of the familiar exponential function. The second choice [15] postulates that The where we have defined the generalized partition function This result differs from the previous in that the role played by (1 − q) is now equivalent to what was played by (q − 1), and that in this case, β is the Lagrange multiplier associated with the constraint on the internal energy. The probability distribution can be conveniently recast as what allows to prove the following series of equalities [20], with F q , U q and C q standing for the corresponding generalizations of Helmholtz free energy, internal energy and specific heat. Then, the formal thermodynamical structure (Legendre transformations) remains valid. Three unwanted consequences, however, should be noted [21] 1. The distribution given by Eqs. (15) and (16) is not invariant against an uniform shift of the energy zero.
2. The mean value of a constant differs from the constant itself, 1 q =1.
3. Finally, if two systems A and B are such that p A+B what differs from U (2) q (A)/k + U (2) q (B)/k. This means that the energy is not an additive quantity.
A third choice for the energy constraint was proposed by Plastino, Mendes and Tsallis [21], and its key aspect is normalization. The constraint is now, i.e. it weights the eigenvalues with a set of escort probabilities p q i / W i=1 p q i . The set {p i } which now extremizes S q are (see the Appendix for the definition of the function involved) wherē It can be shown that, if T ≡ 1/kβ [21], and In addition, use of Eq. (27), together with β∂U q )/∂β leads to Note thatZ (3) q refers to the energy levels {ǫ i }, with respect to U q . We can choose 0 as the energy reference defining Z and recast (27) and (29) as, and Finally, it can be shown that the Legendre structure is preserved too. The appealing features of this third choice are based in that it avoids all three problems that were noted for the second set of possible constraints [21]. It is only necessary to define the new expectation values as where O is any observable. Now, 1 q = 1 ∀q and It is important to mention that the probabilities associated with the third choice coincide with those obtained with the second if we use a normalized temperature [21]. This is the reason by which all theorems that do not use an explicit temperature dependence of a given phenomenon will continue to be valid.

Distribution functions
Despite that the third option is conceptually simpler than the second, actual computations of thermal dependencies are much harder. This stems from the fact that in the second option, the equations for {p i } are explicit (see (15)) whereas they are implicit in the third case (see (24)). The objective of this work is to analyze the evolution of the universe when the distribution functions slightly differ from the standard ones. To this end, working in the second option will provide the simplest expressions possible for all observables, what will allow analytical computations thorough. For particular cases, in addition, it has been proven that the difference between the exact results (using the third choice) and the approximation we shall make here are in very good agreement [22].
Let us now consider a gas composed by N non-interacting particles. If the system is in thermal contact with a heat and a particle reservoir, the energy and the number of particles will be conserved. The stationary state of the system will be given by the solution of where H is the Hamiltonian and ψ R is the wave function of the system. The accessible states will be represented by R. In second quantization formalism, R is given by the set of occupation numbers {n 1 , n 2 , . . . , n k , . . .}, where n k denotes the number of particles in the state k. Taking n k = 0, 1 for fermions ensures Pauli's exclusion principle. We must now extremize the entropy with the imposed constraints, To obtain the solution we use the usual method, extremizing the expression (with units such that k = 1), where α, β are γ Lagrange multipliers. Making, ∂Q/∂P R = 0, we obtain Since it should be zero for all R, each term in the sum must vanish, Identifying β and γ by γ = −βµ, where µ is the chemical potential and T the temperature, the probability that the ensamble is in the state R is, where The occupation numbers automatically determine the quantum state of the system, where ǫ k is the energy of a particle in the state k. Using Eqs. (45) and (46) in the expressions (44) and (43) we obtain, with Z q = n 1 ,...,n k ,...
The partition function can be factorized as It can be shown that (see Ref. [13] for details) the generalized occupation numbers are with ξ = 0, +1 or −1 in the case of a Maxwell-Boltzmann (MB), a Bose-Einstein(BE), or a Fermi-Dirac (FD) gas. In the limit, (q − 1) → 0 standard distribution functions are recovered. The expressions for the distribution functions are not exact ones, since from Eq. (48) to (49) we have made and approximation, a.k.a. factorization approach [13]. In general, indeed, Since the equality is valid for q = 1 we can certainly wait that for values |q − 1| ≪ 1 the approximation will remain a good one. We are interested in the first order in (q − 1), then the expressions (50) becomes, In Fig. 1 we show the change in the distribution functions for different (q − 1)-values.  Eq. (52) is the expression we were looking for, it is simple enough, but yet accurate for small deviations, to give an interesting framework where to analyze (as a parameterization) the effect of a small change on the extensivity properties of the universe on its own evolution. Looking at Fig. 1 it can be seen that despite the particular analytical form for the new distribution functions, it is clear that they will do an excellent job in showing an slight deviation from the standard case, this being the ultimate reason why we use them.
Up to recent days, quantum distribution functions within the factorization approximation were regarded as a rather rough technique, see Pennini et al. [23]. These latter authors considered fermion and boson systems with very small occupation numbers. However, Wang and Lé Méhauté [24] analyzed the problem in detail and showed that there exist a temperature interval, which they called 'forbidden zone', where the deviation from the exact result maybe significant. Otherwise, outside this zone, the factorization approach results can be used with confidence. In addition, they verified that the magnitude of the forbidden zone remained constant with the increase of the number of particles. This fact motivated new efforts in the study of macroscopic systems (where the number of particles is large), since the generalized distribution functions of the factorization approach could be used at temperatures up to 10 20 K for such systems [24]. All these results encourage us further to use this approximation for the analytical study that follows.
3 The plasma parameter, the plasma era, and before One of the earliest applications of a non-Maxwellian distribution to plasmas was referred to the solar core. Clayton considered an ionic quantum distribution function that slightly deviated from a BG one, being the correction factor proportional to exp[−δ(E/kT ) 2 ] [25]. δ can be related with a non-extensive parameter q by means of the formula δ = (1 − q)/2, see [3]. There are several reasons why to expect a slight deviation from a BG behavior in the sun core. Particularly, we know that long range forces are present (Coulomb and gravitational) and that there are also effects related with long range microscopic memory. Then, the leading idea is try to see whether the solar core is a non-ideal plasma, or a quasi-plasma, or equivalently, to what extent it can be considered as ideal. While we forward the reader looking for a detailed account of these effects to the references quoted in [3], let us just mention here the computation of the plasma parameter.
The plasma parameter is the ratio between typical electrostatic energy of the nearest neighbors of a given particle and the thermal energy kT . For a plasma with number density n, the mean inter-particle distance will be n 1/3 , and it can be proven that the number of particles inside a Debye sphere is given by [26] The plasma parameter will be then defined as Γ = N −2/3 d . When Γ ≪ 1, the plasma behaves (just considering the long range forces) as an ideal gas, the electrostatic energy is small compared with kT and there are a large number of particles that shield the long range force in a Debye-radius range. If, instead, Γ ≫ 1, the plasma is strongly coupled and certainly non-ideal. For values of Γ near 1 we expect slight corrections. For the sun core, for instance, the plasma parameter is of the order of 0.1 Is it small enough to produce no appreciable deviation? The answer seems to be no. It have been recently proven, see last reference in [3], that the effects of random electric microfields are of crucial importance. These fields have in general long-time and long-range interactions that can generate anomalous diffusion that alters the distribution. Most importantly, the amount of the deviation has been quantified to be proportional to Γ 2 .
The plasma parameter in the early universe is used to define what is known as the plasma era. At the end of the plasma era is where the process of recombination occurs. Here, the plasma parameter is very small and an ideal gas approach is, a priori, a very good one. If the universe is well described within the standard model of cosmology, we do not expect the distribution functions to deviate from the usual BG ones within the plasma era and in later epochs. Then, if we force the distributions to deviate, the corrections should turn out to be very big, in such a way that experiments could then impose restrictive bounds on, for instance, the non-extensive parameter. This is indeed what happens, as we shall see below particularly when analyzing the recombination process. Do the same happens for all eras of cosmic evolution?
Near the temperature of freezing at about 1 MeV, something very important happens in the universe: the electron-positron annihilation. At about 0.5 MeV, the cosmological plasma changes its composition, before the annihilation being formed mainly by electrons, positrons, photons, neutrinos, and traces of nonrelativistic particles. Which is the value of the plasma parameter in this era? This plasma is relativistic, and then we cannot expect that the previously introduced parameter (although indeed it grows as the temperature of the universe increases) to be a complete description. We can, however, consider it as a first approximation. Its value is ∼ 0.07 at the time of freezing (t ∼ 1s, T ∼1 MeV). This value of Γ is close to the sun core case. Can we a priori expect a similar correction to the quantum distribution functions? As we shall later show, when relativistic particles dominate the problem, corrections introduced by a slight deviation in the distribution functions are not so big as in the non-relativistic case. Restrictions on the parameter space are then not so strong, and, a posteriori, it is not possible to straightforwardly discard such slight deviations. In addition, if Γ 2 gives the amount of the deviation also for this system, we would expect corrections of the order of 10 −2 − 10 −3 , what we do obtain (independently of these comments) when we analyze the epoch of primordial nucleosynthesis, see paper II in this series [27].
We do not claim here that these paragraphs provide evidence of such an strong nature as to consider a priori that deviations to BG distributions functions are unavoidable, as the analysis for the sun core case seem to show [3]. We do believe, however, that they make worth looking at the standard model predictions in an statistically modified setting, to see the allowed range of the deviations along the different cosmic eras. To that aim we devote this and the follow up paper in the series. On the other hand, we recall that several authors, with various motivations, have considered other effects that could slightly change the distribution functions as well: Non-equilibrium effects originating in residual interactions between neutrinos and electrons (see for instance [28]) and fluctuations in magnetic fields possibly present at that time (see for instance [29]) are some examples.
4 Thermodynamics in the expanding universe 4

.1 Energy density, pressure, number of particles
The number density, the energy density, and the pressure for a dilute gas are given in terms of their corresponding distribution function f (p) by (units are such thath = c = 1), where E(|p|) = |p| 2 + m 2 , and g is the degeneracy factor of the particle. From now on, E(|p|) = E(p) and |p| = p. Let us first take the non-degenerate case, in which kT ≫ µ. As before, T will stand for the temperature and µ for the chemical potential. Use of the form of f (p) in the non-extensive approach we adopted gives the correction to the standard result at first order in (q − 1). Several cases are worth mentioning.
When the particles are relativistic, E ≫ m, the energy density is with ρ st representing the usual result, and g = g bosons + 7/8g f ermions . For the number of particles, we obtain, For the pressure, starting from the definition, we can immediately prove that P = ρ/3 disregarding the degree of non-extensivity.
Remark 1: Note that, for relativistic particles, the temperature dependence of ρ, P and n in the nonextensive formalism is the same as the standard. This is because it comes only as a product of the change of variables in the integrals, which is done exactly in the same way for both contributions.
It is interesting now to compute the energy per particle, which results to be Note that to equal q, the fermions get more affected than the bosons by the change of the distribution. In the non-relativistic limit m ≫ T , we can neglect the ±1 term in the denominator of all integrals. For instance, the number density becomes, Here, the first term in the right hand side stands for the standard result and the second three for the correction at (q − 1) order. The energy density is then given by ρ q = mn q and the pressure is P q = n q T .
Remark 2: For non-relativistic particles, the change to a non-extensive setting does change the temperature dependence of the observables. This is not a minor effect and will show its impact in the recombination process.

Particle anti-particle excess
We focus now on the computation of the particle anti-particle excess. Consider the reaction r + +r − ↔ γ+γ, where r + and r − represent a generic particle and its anti-particle. If this reaction is in equilibrium, then the chemical potential of r + and r − are related by µ + = −µ − . Then, the net number of fermions is given by where the integration variable is, as before, x = E/T . In the relativistic limit, m/T ≪ 1, and the standard integrals can be solved as To compute the (q − 1) order corrections we need to calculate Since we would like to have an analytical solution, we think in a new application of the free parameter trick introduced in this context by Tırnaklı and Torres [22]. We know that Here, m is considered as a free parameter. Then, the integral we need can be obtained by making the following steps. Consider −d 2 I m /dm 2 , evaluated at m = 1. This last integral will give the one we want plus two additional terms. These terms are We can evaluate the relative strength of these extra terms with respect to the integrals we are looking for. The first thing one can analytically do is to consider the case in which ψ ∼ 0. This, we know, is not far from the real situation, since after electron positron annihilation, only a little excess of electrons will survive, in order to give the universe, together with protons, electrical neutrality. Because the ratio of the number density of protons to the number density of photons, n p st /n st γ ∼ 10 −9 , use of the standard expression for n + st − n − st and n st γ will show that and clearly the condition µ/T ≪ 1 is sustained. Taking this into account we can weight each of these two extra terms with respect to our integrals. These terms are, and If we call f A (u) and f B (u) to the integrands of the previous expressions, we note that Around 0, f B (u) ≃ f A (u); but the integrand B falls exponentially. Finally, a numerical integration of the two terms, that we show in Fig. 1, yields a difference bigger than one order of magnitude: A = 22.63 and B = 1.037. For the sake of analytical computation, we shall adopt the following approximation for the correction to the particle anti-particle excess n c , which finally yields the result, In the case of non-relativistic particles, a similar computation gives with the standard term being and the correction being

Effective number of degrees of freedom
The total energy density of all species in equilibrium can be expressed in terms of the photon temperature.
As most of the energy density is contributed by relativistic species, we shall firstly consider only them.
The total density will be the sum over all particular species, Note that we are using a sub-index i for boson particles, whereas j is used for fermions. From previous expressions we have, taking x r = E r /T , recall here that f r q (x r ) is the non-extensive distribution function, and we are in the limit in which E r ≫ m r and x r ≫ m r /T r are valid. If we as well consider T r ≫ µ r , then f r q (x r ) ≡ f r q (x). This gives the separate corrections, The standard contributions are the usual, Adding up both contributions, we can write where g q * is given by the addition of the standard, and the new correction, Remark 3: Note that, due to the same temperature dependence, we can define an effective number of degrees of freedom exactly in the same way as it is done in the standard case. Then, we can consider the same evolution equations, and hide all non-extensivity effects in the new g q * .
As an example of the previous remark, we can consider the Friedmann equation. It reads, where R(t) is the scale factor in the FRW metric, and k is the curvature of the space-time. In the early universe, it is a very good approximation to take k = 0. During the radiation dominated phase P = ρ/3, R(t) = t 1/2 , andṘ/R = H(t) = 1/t. Using the expression for the energy density, i.e. Eq. (83), we obtain Then, although the functional form is the same as in the standard case, the actual time-temperature relationship is different because of the change in the quantum distributions.

Conserved quantities
From the full set of Einstein field equations, or from the fact that T µν ;ν = 0, it may be established that This can be written as And using the previous equation we finally get d dT as can be checked by direct differentiation.
Remark 4: Disregarding the degree of non-extensivity, s q ≡ (ρ q + P q )/T is a conserved quantity in a comoving volume.
If the relativistic contribution is dominant, This yields where we have defined g q * ,s as g q * ,s = g st * ,s + g c * ,s , From the conserved comoving quantity, we may define the number S q = s q R 3 = constant, and obtain the temperature-scale factor relationship as Then, as in the standard case, T ∝ R −1 , but now only in the periods where g q * ,s (and not g * ,s ) is constant.

Explicit form for the conserved numbers
For reasons that clarify themselves when actually making computations, it is convenient to defineĝ c ⋆s such that g c ⋆s = (q − 1)ĝ ⋆s . As in the standard case, we have proven that for any particle which is not created nor destroyed during the evolution of the universe, is a conserved quantity. For future use, we quote here the explicit form of the conserved number for relativistic particles. For instance, for photons, the number density is n γ q = g/(2π 2 ) [2ζ(3) + 12ζ(4)(q − 1)] T 3 , with g = 2. Then, This gives, Similar computations, but now for fermions, would yield, with

Decoupling
Once a given species of particles has decoupled completely, each particle will travel along a geodesic of space-time. This will ensure that the form of the distribution function will remain the same during the subsequent evolution of the universe. Indeed, the only change will be given by a redshift correction. If the decoupling is produced at t = t d , then, due to the redshift, all particles that at the instant t have momentum p, should have had momentum pR(t)/R(t d ) at decoupling. Then, Remark 5: The previous result is independent of the particular f eq (p, t) we have adopted.

Relativistic decoupling
Let us write the rhs of Eq. (104) taken into account that E ≃ p and in the case in which µ/T ≃ 0, We see that the distribution function has the same functional form than f eq q but with a different "temperature" given by, Despite this species is no longer in thermodynamical equilibrium, the "temperature" in the distribution function falls as R −1 ; making S r q = s r q R 3 to conserve separately. Note that for the species still in equilibrium, the temperature falls as T ∝ (g q ⋆s (T )) −1/3 R −1 . The number density for these particles is then given by, and using that for any time t, after t d , R(t) , the previous equations really are, These number densities are comparable to that of photons at any given time. Such population will then continue to exist as a relativistic relic.

Non-relativistic decoupling
In the case of non-relativistic decoupling, the distribution function is given by and the number density is given by If we now consider that for a particle species which has already decoupled, the number density must be such that n q ∝ R −3 , we see that as in the standard case (see for instance Section 3.4 of Ref. [30]), we need to identify the temperature with T d R 2 (t d )/R 2 (t), and ask for the chemical potential to vary as µ = m − (m − µ d )T /T d . Then, the distribution function is similar to an equilibrium distribution but with a new temperature given by decreasing with the square of the scale factor. In the limit m ≫ T , the energy density will be given by ρ = nm.
Also in a non-extensive framework, the correct way to handle with decoupling process is through the Boltzmann equation, which we analyze below.

The relic neutrino background
The conservation of s q , applied to the particles which are in equilibrium with radiation shows that the quantity g q * ,s T 3 γ R 3 is constant throughout the expansion. During pair annihilation, g q * ,s will decrease, since we go towards a state with less number of degrees of freedom. Then, T 3 γ R 3 after the process will be bigger than its previous value. We say that photons are heated up by the annihilation process, in such a way that is sustained. This directly yields a relationship between both temperatures, (1 + 0.109(q − 1)) Note that the first term is the standard result, from the usual hot big bang model. The second term is the non-extensive correction.
Since before e + e − annihilation, neutrinos were in thermal equilibrium with photons, T bef ore γ will be the same as T bef ore ν . Now notice that neutrinos do not participate in the process e + + e − ↔ 2γ, and so, neutrino temperature is constant. Then, T bef ore ν = T af ter ν , and we can write After annihilation, g q * ,s does not change anymore, and both temperatures T af ter γ and T af ter ν fall as R −1 , conserving their ratio. Since now T γ ∼ 2.728, then T ν = 1.947(1 − 0.109(q − 1))K.
Remark 6: Even if the evolution of the universe diminish the degree of non-extensivity, the relic backgrounds, produced at early times, are sensitive to the possible non-extensivity present at that epoch.
Remark 7: Note that if (q − 1) could have any value, one can have a photon temperature unchanged after e + e − annihilation: T af ter γ /T bef ore γ = 1. This effect particularly show how important is the statistical description in the evolution of the early universe. Note, however, that a posteriori constraints on the value of (q − 1) do apply, and this situation will turn out to be impossible.
The difference between T ν and T γ can be bounded using information coming from the CMBR, see for instance the work by Torres [11].

Recombination
We focus now on the generalization of the Saha law. Let n H , n p and n e the number density of hydrogen atoms, protons, and electrons respectively. Electrical neutrality implies that n e = n p . Baryon total number is given by n B = n H + n p . In thermal equilibrium, a temperatures such that T < m i , we have seen that n i is given by Eq. (63). The process p + e → H + p guarantees that µ p + µ e = µ H . We can as well consider that m p ≃ m H . Then, we would like to find an expression for n H q /n e q n p q . Immediate algebra yields, Here, we have made use of the definition, We would like to introduce now the number of baryons, since this is a conserved number (as far as precision measurements can tell), i.e. N B = n B /s is constant. Note however that η q = n B q /n γ q does not remain unchanged, since g q * ,s may change in time. Using g p = g e = 2 and g H = 4 we can write, We now define X e q = n p q /n B q and 1 − X e q = n H q /n B q . Using them, we have where B = m p + m e − m H and n q γ = n γ,st + (q − 1)n γ,c . Note that we have shortened the notation on the last bracket, but we are actually implying the definition of Eq. (118). Also, Finally, we rewrite with and we use that Recalling that n γ st = (2ζ(3)/π 2 )T 3 , we can write, with g H =, and g p = g e = 2, And taking into account the standard result, we see that we have arrived to the following relationship: This is the generalized Saha law. Looking at Eq. (127) we may note that we are going to confront a problem that does not appear in the standard case: The chemical potentials appear in a separate way, in such a way that not only the relationship µ H − µ e + µ p = 0 enters in the computation. Since it is not common to find the µ-values, or even the method to find them, in standard cosmology books, we shall briefly outline here how they can be naively obtained. We shall then obtain the standard chemical potentials µ = µ st , and consider that generalized expressions for them will be written as µ q = µ st + (q − 1)C. Since the term including the factor C will provide a second order correction in Eq. (127) we shall be concerned only with the standard computation. We have to consider a system of three equations: 1. chemical equilibrium µ H − µ e + µ p = 0 2. electrical neutrality n e = n p 3. the standard result for 1 − X e st /(X e st ) 2 . The fraction X e st is shown in Fig. 3.
These three equations form a complicated system but with three unknowns, regarded here as (µ i − m i )/T , where i runs over the three species involved. Although this system cannot be solved analytically, it can be solved in a numerical way. The exact values of the solution will ultimately depend on the value of η, the baryon to photon ratio, which in turns depends on the parameter Ωh 2 . Results for the chemical potentials are shown in Table 1 and in Fig. 3.
With these values, the correction terms in the Saha law, using Eq. (127), can be immediately computed. The results are shown in Fig. 4, where two different cases for the baryon content of the universe are studied. As the corrections in the non-relativistic expression for the number of particles is proportional to the factors ((µ i − m i )/T ) 2 , we note that leading order corrections in the period we are interested in are ∼ (q − 1)50 2 = (q − 1)2500. These are large factors, and are the main reason by which all bounds imposed at the time of recombination, where matter and radiation cease to be in thermal contact, i.e. using CMBR data, are much stronger than those imposed at earlier epochs.

The Boltzmann equation
To analyze in an appropriate way the decoupling of particles we must consider the microscopic evolution of the distribution functions f (p µ , x µ ), via the Boltzmann equation. An interesting recent study on this issue  Figure 4: Final corrections in the Saha law. Note that the relative strength of the correction differs appreciably for different baryon densities. Note that the interesting range for the correction is limited by the equilibrium in the reaction D ←→ p + n, represented by µ d = µ n + µ p .
was presented in Ref. [12], but their emphasis was put on other aspects. We shall start here by writing the Boltzmann equation as,L [ where C is the collisional operator andL is the Liouville operator. The expression for the latter, in the non-relativistic case, for the phase space density f ( v, x) of a species with mass m under a force F = d p/dt, The relativistic generalization becomes,L Gravitational forces are present through the connection Γ α βγ . For a Friedmann-Robertson-Walker universe, f = f (| p|, t) (or equivalently f = f (E, t)) and the only non-null components of the connection Γ α βγ are the usual with h ij = −g ij and Latin indices running from 1 to 3. Taken this into account and recalling that p µ = (E, p) and x µ = (t, x), Eq. (131) becomeŝ Substituting this expression in Eq. (128) and taking into account the definition of the number density, we can integrate by parts to obtainṅ We can also rewrite the collisional term defininġ in order to have the Boltzmann equation as, The collisional term is usually dominated by annihilations between particles and anti-particles (pp) and we assume now that the number of particles is identical to the number of anti-particles. In this way,ḟ c is given by (see Section 9.2 of Ref. [31] for details): where < σv > is the product of the cross section times the velocity, averaged in velocities space. We can extract < σv > out of the integral, evaluating it in an averaged energy. Let us focus in ψ-particles and their anti-particlesψ. Taking into account Eq. (139) we integrate over momentum space in Eq. (138) to arrive aṫ Recall that the index q refers to the non-extensive distribution functions and generalized quantities. We can now account for the thermal production of particles. Naming Υ the term giving this contribution, we can write the Boltzmann equation aṡ This term can be spelled out as in the standard case [31], invoking thermodynamical equilibrium: in an static universe (Ṙ = H = 0), n ψ q would be a constant and equal to n ψ q,eq ; then, We can rewrite the previous equation in terms of the variable Y q ≡ n ψ q /s q , where Y q,eq ≡ n ψ q,eq /s q , resulting, All this resembles the usual derivation, and indeed we can go still further. We can writeṅ ψ q + 3Hn ψ q in terms of Y q . We shall take into account that s q R 3 = constant, and we derive with respect to time both the definition for Y q and the equation s q R 3 = constant. Finally recalling that H =Ṙ/R, we geṫ n ψ q + 3Hn ψ q = s qẎq .
Combining the last two expression forṅ ψ q + 3Hn ψ q , we obtain: Introducing explicitly the temperature, via x = m/T , with m being any appropriate mass scale, and using the radiation-dominated-era relation between x and T , Eq. (88), we get and we can writeẎ Using the Jacobian dt/dx, we can rephrase forẎ q where we have defined Note that H(m) is related with H by the simple expression Finally, using Eqs. (145) and (148) we can write the Boltzmann equation as where < σv > is given by Considering other ψψ annihilation channels, say ψψ into a final state F , an additional term would appear that would be similar to that in the rhs of the previous equation, but with < σ ψψ→xx v > replaced by < σ ψψ→F v >. Sum over all annihilation channels yields the final result in terms of the effective cross Defining Γ A = n ψ q,eq < σ A v >, we get Remark 8: Note that this derivation is independent of the particular form for the distribution functions, but the final product has the same aspect than that obtained in the standard scheme, see for instance Section (5.2) of Ref. [30].

Freezing
Following the standard procedure, we would like to study how much the freezing times and temperatures are modified by a change in the statistical description. Using Eq. (154), we see that the comoving number density of ψ particles is controlled by the ratio Γ/H, and the amount of the deviation of the distribution functions from their equilibrium values. When Γ/H is less than unity, the change in the number density is small, −∆Y q /Y q ∼ −(xdY /dx)/Y q,eq ∼ Γ/H < 1, annihilations stop and the number density freezes.
The annihilation rate Γ A varies as n q,eq times the averaged cross section < σ A v >. In the relativistic regime, n q,eq ∝ T 3 . In the non-relativistic regime, n q,eq ∝ e −m/T , and Γ A exponentially decreases. In both regimes, Γ A decreases with temperature, and eventually goes below the necessary rate to maintain equilibrium, which by definition occurs when x = x f (the "freeze out"). Then we expect that for We can compute the value of Y q,eq ≡ n ψ q,eq /s q in each of the previous cases recalling the definition of s q and Eqs. (59) and (60) for relativistic (x ≪ 3) particles and Eq. (63) for non-relativistic (x ≫ 3) particles. We obtain: for relativistic bosons, Y f q,eq = 0.278 for relativistic fermions, and Y nr q,eq = 0.145gx 3/2 e −x (g q ⋆s ) −1 1 + q − 1 2 with x ≫ 3, for non-relativistic particles.
Since we can write g q ⋆s as g q ⋆s = g st ⋆s + (q − 1)ĝ c ⋆s , we can get a first order result for the above mentioned quantities, Y f q,eq = 0.278 and Y nr q,eq = 0.145gx 3/2 e −x (g st ⋆s ) −1 1 + q − 1 2 In Fig. 5, we show the equilibrium abundance for non-relativistic particles together with the nonextensive correction. The correction diminishes exponentially with decreasing temperature, but this is an effect of the overall factor e −x also present in the standard result, and not of the non-extensivity introduced.

Current values
We are interested in the value of the corrected ratio where G = m −2 pl with m pl = 1.2211 × 10 19 GeV and H(t) is H 0 = 2.1332h × 10 −42 GeV (0.4 < h < 1). From the current value of ρ R q , we can get the contribution of relativistic particles To get ρ R q we must know the the value of g q ⋆ and g q ⋆s . Carrying out a first order computation using Eqs. (95), and (96), and Eqs. (84) and (85), and assuming a universe populated by photons and three types of neutrinos, we get (g q ⋆ ) today = 3.36 [1 + 9.43(q − 1)] (g q ⋆s ) today = 3.91 [1 + 7.11(q − 1)] Taking the current value of the photon temperature, T = 2.728 ± 0.004 K, we can obtain: ρ R q , s q , n γ q and Ω q,R , (Ω q,R ) today = 4.17 × 10 −5 [1 + 9.43(q − 1)] h −2 .
10 Wimps in non-extensive statistics We now consider weakly interacting massive particles, for instance neutrinos with very small mass, in non-extensive statistics.

Relativistic relics, x f < 3
In this case, Y q,eq is not changing with time, (see Eqs. (158) and (159)). The asymptotic value of Y q , (Y q (x → ∞)) ≡ Y q,∞ , is just the equilibrium value at freezing: The current abundance of the ψ species today is n ψ q,today = The mass density contribution of this kind of hot relics is given by and the fraction of the critical mass will be We know that Ω 0 h 2 < 1; then we can apply this bound to the contribution to the species ψ to Ω 0 h 2 and obtain a cosmological upper limit to the mass of this kind of particles m ψ < 13.08 Then, for an species with two components, g f = g ν = 2, and to first order in (q − 1), it results This is a non-extensive generalization of the Cowsik and McClelland bound, see for instance Ref. [30]. Even when is improbable that an experiment with enough sensitivity could differentiate among values of q using Eq. (177), further recalling that it is an upper bound what that equation is imposing, it is worth noticing that a change of statistics can have a direct influence upon the mass spectrum of the particles.

Cold relics, x f > 3
In this case, Now, Y q,eq strongly depends on m (since x f = m/T d ). To make a numerical estimation we have to determine T d , using the condition Γ ≃ H. Reactions able to change the number of wimps A, are like AĀ ↔ XX, where X is a generic particle, assumed in thermal equilibrium. The average value of σv can be expressed as The value of the power law index k will depend on the details of the annihilation process, usually it is of order 1. The value of σ 0 depends on m and has a simple form for the extreme cases m ≪ m Z and m ≫ m Z , where m Z ≃ 10 2 Gev is the mass of the Z-boson.
For wimps with m < m Z , the effective cross section σ 0 can be written as [26], where G F is the Fermi constant, G −2 F = 292.8 GeV. The value of c depends on the fermion. We shall consider Dirac (spin 1/2) particles for which c ≃ 5.
The reaction rate is given by or, equivalently, using Eq. (63), The expansion rate H as a function of temperature is Since g q ⋆ is a slowly varying function we can solve this equation in an iterative form. Consider the case of a wimp with m > 1 GeV. From Eq. (188) we see that m/T d ≃ 17.74, assuming for simplicity that k = 0, the term ln(m/T d ) corrects the value m/T d to m/T d ≃ 19.18, giving T d ≃ 52 MeV (m/GeV). At this temperature, all species are relativistic and all of them contribute to g q ⋆ , making it of the order of 100. Taking this into account, the term ln((g q ⋆ ) 1/2 /g A ) ≃ ln 5 ≃ 1.61, and the value of m/T d gets corrected to m/T d ≃ 17.57. The non-extensive term is ln 1 + (q − 1) 1 2 The fermion A and its anti-particleĀ, will provide twice this value of Ω q,today h 2 , i.e. the constraint Ω q,today h 2 < 1 is then a corrected lower bound for the mass of these particles m > 1.64 [1 − 37.6(q − 1)] GeV. (194) 11 Modification of the matter-radiation equality Let us characterize the matter-radiation equality by the time of its occurrence t = t eq , its scale-factor R = R eq , and its redshift z = z eq . From the equalities ρ matter (t eq ) = ρ matter (t 0 ) ρ radiation (t eq ) = ρ radiation (t 0 ) it follows 1 + z eq = Ω matter (t 0 ) Ω radiation (t 0 ) ≃ Ω Ω radiation (t 0 ) .

Conclusions
In this paper we have studied how the statistical description that we suppose is valid in the early moments of the universe affects different processes along the thermal history of our basic cosmological model. In particular, starting from the modifications on the energy and number densities, we explored all consequences introduced by non-extensivity in the processes of decoupling, particle-anti-particle excess, recombination, photon-neutrino temperature relationship, matter-radiation equality, and others. In order to avoid re-stating our obtained results, let us just mention a general conclusion. Suppose that all cosmological observables are measured with a precision of 1%. Suppose, too, that the standard values for these observables agree well with observations, i.e. the standard values are always within the error bar of the experiments. (This latter assumption is not always true, particularly when more than one experiment is involved at the same time, like in the case of the determination of the baryon density using nucleosynthesis and recent CMBR measurements.) When relativistic particles are involved, we can summarize our results as follows. Given any standard prediction, the correction introduced by non-extensivity is of the order of ten times the standard value, where the factor F is order 10X st . Then, given the previous assumptions, when relativistic particles are involved, |q − 1| is bound to be of the order 10 −3 . On the other hand, when non-relativistic particles are involved, the corrections are much bigger, and F can be at least of the order of 1000X st . In these cases, and also within the previous assumptions, |q − 1| is bound to be less than 10 −4 . It is worth noticing that this simple conclusion is sustained by all currently existing bounds, and it is also worth noticing that the ultimate reasons by which the constraints on |q − 1| are more restrictive in the epoch of recombination than those obtained in nucleosynthesis (see Paper II) can be tracked from this formal point of view. In our second paper in this series we shall study the nucleosynthesis epoch, and will provide an analytical assessment of the Helium 4 and Deuterium primordial production in the framework of non-extensive statistics, including corrections coming from the free neutron decay process. We shall particularly focus on how non-extensivity affects the principle of detailed balance.

Appendix: Useful mathematical formulae
It is convenient to introduce the following functions: supplemented by the definition, for q < 1, e x q = 0 if 1 + (1 − q)x ≤ 0, and for q > 1, e x q diverges in x = 1/(q − 1), and It is immediate to show that lim q→1 e x q = e x and that lim q→1 ln q x = ln x. Also the following properties are valid e lnq x q = ln q e x q = x, ∀(x, q).
With these definitions,