BIG BANG NUCLEOSYNTHESIS WITH A NON-MAXWELLIAN DISTRIBUTION

The cosmological big bang model is in agreement with many observations relevant for our understanding of the universe. However, comparison of calculations based on the model with observations is not straightforward because the data are subject to poorly known evolutionary effects and systematic errors. Nonetheless, the model is believed to be the only probe of physics in the early universe during the interval from 3 to 20 minutes, after which the temperature and density of the universe fell below that which is required for nuclear fusion and prevented elements heavier than beryllium from forming. The model is inline with the cosmic microwave background (CMB) radiation temperature of 2.275 K (Noterdaeme et al. 2011) and provides guidance to other areas of science, such as nuclear and particle physics. Big bang model calculations are also consistent with the number of light neutrino families N_ν = 3. According to the numerous literature on the subject, the big bang model can accommodate values between N_ν = 1.8 and 3.9 (see, e.g., Keith 2002). From the measurement of the Z₀ width by LEP experiments at CERN one knows that N_ν = 2.9840 ± 0.0082 (LEP Collaboration 2006).

In the big bang model nearly all neutrons end up in ⁴He, so that the relative abundance of ⁴He depends on the number of neutrino families and also on the neutron lifetime τ_n. The sensitivity to the neutron lifetime affects big bang nucleosynthesis (BBN) in two ways. The neutron lifetime τ_n influences the weak reaction rates because of the relation between τ_n and the weak coupling constant. A shorter (longer) τ_n means that the reaction rates remain greater (smaller) than the Hubble expansion rate until a lower (larger) freezeout temperature, having a strong impact on the equilibrium neutron-to-proton ratio at freezeout. This n/p ratio is approximately given in thermal equilibrium by n/p = exp [ − Δm/k_BT] ∼ 1/6, where k_B is the Boltzmann constant, T is the temperature at weak freezeout, and Δm is the neutron–proton mass difference. The other influence of τ_n is due to their decay in the interval between weak freezeout (t ∼ 1 s) and when nucleosynthesis starts (t ∼ 200 s), reducing the n/p ratio to n/p ∼ 1/7. A shorter τ_n implies a lower predicted BBN helium abundance. In this work we will use the value of τ_n = 878.5 ± 0.7 ± 0.3 s, according to the most recent experiments (Serebrov et al. 2005; a recent review on the neutron lifetime is found in Wietfeldt & Greene 2011). Recently, the implications of a change in the neutron lifetime on BBN predictions have been assessed in Mathews et al. (2005).

The baryonic density of the universe deduced from the observations of the anisotropies of the CMB radiation constrains the value of the number of baryons per photon, η, which remains constant during the expansion of the universe. Big bang model calculations are compatible with the experimentally deduced value from Wilkinson Microwave Anisotropy Probe (WMAP) observations, η = (6.16 ± 0.15) × 10⁻¹⁰ (Komatsu et al. 2011).

Of our interest in this work are the abundances of light elements in BBN. At the very early stages (first 20 minutes) of the universe evolution, when it was dense and hot enough for nuclear reactions to take place, the temperature of the primordial plasma decreased from a few MeV down to about 10 keV, light nuclides as ²H, ³He, ⁴He and, to a smaller extent, ⁷Li were produced via a network of nuclear processes, resulting in abundances for these species which can be determined with several observational techniques and in different astrophysical environments. Apparent discrepancies for the Li abundances in metal-poor stars, as measured observationally and as inferred by WMAP, have promoted a wealth of new inquiries on BBN and on stellar mixing processes destroying Li, whose results are not yet final. Further studies of light-element abundances in low-metallicity stars and extragalactic H ii regions, as well as better estimates from BBN models are required to tackle this issue, integrating high-resolution spectroscopic studies of stellar and interstellar matter with nucleosynthesis models and nuclear physics experiments and theories (Fields 2011).

The Maxwell–Boltzmann distribution of the kinetic energy of the ions in a plasma is one of the basic inputs for the calculation of nuclear reaction rates during the BBN. The distribution is based on several assumptions inherent to the Boltzmann–Gibbs (BG) statistics: (1) the collision time is much smaller than the mean time between collisions, (2) the interaction is local, (3) the velocities of two particles at the same point are not correlated, and (4) that energy is locally conserved when using only the degrees of freedom of the colliding particles (no significant amount of energy is transferred to and from collective variables and fields). If (1) and (2) are not valid, the resulting effective two-body interaction is non-local and depends on the momentum and energy of the particles. Even when the one-particle energy distribution is Maxwellian, additional assumptions about correlations between particles are necessary to deduce that the relative-velocity distribution is also Maxwellian. Although the BG description of statistical mechanics is well established in a seemingly infinite number of situations, in recent years an increasing theoretical effort has concentrated on the development of alternative approaches to statistical mechanics which includes the BG statistics as a special limit of a more general theory (Tsallis 1988; see also Rényi 1960). Such theories aim to describe systems with long-range interactions and with memory effects (or non-ergodic systems). A very popular alternative to the BG statistics was proposed by C. Tsallis (Tsallis 1988; Gell-Mann & Tsallis 2004), herewith denoted as non-extensive statistics (for more details on this subject, see the extensive reviews Gell-Mann & Tsallis 2004; Tsallis et al. 2005; Tsallis 2009). Statistical mechanics assumes that energy is an "extensive" variable, meaning that the total energy of the system is proportional to the system size; similarly the entropy is also supposed to be extensive. This might be justified due to the short-range nature of the interactions which hold matter together. But if one deals with long-range interactions, most prominently gravity, one can then find that entropy is not extensive (Fa & Lenzi 2001; Lima et al. 2002; Taruya & Sakagami 2003, 2004; Chavanis & Sire 2005).

In classical statistics, to calculate the average values of some quantities, such as the energy of the system, the number of molecules, the volume it occupies, etc., one searches for the probability distribution which maximizes the entropy, subject to the constraint that it gives the right average values of those quantities. As mentioned above, Tsallis proposed to replace the usual (BG) entropy with a new, non-extensive quantity, now commonly called the Tsallis entropy, and maximize that, subject to the same usual constraints. There is actually a whole infinite family of Tsallis entropies, indexed by a real-valued parameter q, which quantifies the degree of departure from extensivity (one gets the usual entropy back again when q = 1). It was shown in many circumstances that the classical results of statistical mechanics can be translated into the new theory (Tsallis 2009). The importance of these families of entropies is that, when applied to ordinary statistical mechanics, they give rise to probabilities that follow power laws instead of the exponential laws of the standard case (for details on this see Tsallis 2009). In most cases that Tsallis formalism is adopted, e.g., Pessah et al. (2001), the non-extensive parameter q is taken to be constant and close to the value for which ordinary statistical mechanics is obtained (q = 1). Some works have also probed sizable deviations of the non-extensive parameter q from unity to explain a variety of phenomena in several areas of science (Gell-Mann & Tsallis 2004).

In the next sections, we show that the Maxwell–Boltzmann distribution, a cornerstone of the big bang and stellar evolution nucleosynthesis, is strongly modified by the non-extensive statistics if q strongly deviates from unity. As a consequence, it also affects strongly the predictions of the BBN. There is no "a priori" justification for a large deviation of q from the unity value during the BBN epoch. In particular, as radiation is assumed to be in equilibrium with matter during the BBN, a modification of the Maxwell distribution of velocities would also impact the Planck distribution of photons. Recent studies on the temperature fluctuations of the CMB radiation have shown that a modified Planck distribution based on Tsallis statistics adequately describes the CMB temperature fluctuations measured by WMAP with q = 1.045 ± 0.005, which is close to unity but not quite (Bernui et al. 2007). Perhaps more importantly, Gaussian temperature distributions based on the BG statistics, corresponding to the q → 1 limit, do not properly represent the CMB temperature fluctuations (Bernui et al. 2007). Such fluctuations, allowing for even larger variations of q, might occur during the BBN epoch, also leading to a change of the exponentially decaying tail of the Maxwell velocity distribution.

Based on the successes of the big bang model, it is fair to assume that it can set strong constraints on the limits of the parameter q used in a non-extensive statistics description of the Maxwell–Boltzmann velocity distribution. In the literature, attempts to solve the lithium problem have assumed all sorts of "new physics" (Fields 2011). This work adds to the list of new attempts, although our results imply a much wider impact on BBN as expected for the solution of the lithium problem. If the Tsallis statistics appropriately describes the deviations of tails of statistical distributions, then the BBN would effectively probe such tails. The Gamow window (see Figure 1) contains a small fraction of the total area under the velocity distribution. Thus, only a few particles in the tail of the distributions contribute to the fusion rates. In fact, the possibility of a deviation of the Maxwellian distribution and implications of the modification of the Maxwellian distribution tail for nuclear burning in stars have already been explored in the past (Lissia & Quarati 2005; Haubold & Kumar 2008; Degl'Innocenti et al. 1998; Coraddu et al. 1999). As we show in the next sections, a strong deviation from q = 1 is very unlikely for the BBN predictions, based on comparison with observations. Moreover, if q deviates from the unity value, the lithium problem gets even worse.

Zoom In Zoom Out Reset image size

Nuclear reaction rates in the BBN and in stellar evolution are strongly dependent on the particle velocity distributions. The fusion reaction rates for nuclear species 1 and 2 are given by 〈σv〉₁₂, i.e., an average of the fusion cross section of 1 + 2 with their relative velocity, described by a velocity distribution. It is thus worthwhile to study the modifications of the stellar reaction rates due to the modifications introduced by the non-extensive statistics.

2.1. Non-extensive Statistics

Statistical systems in equilibrium are described by the BG entropy,

$$\begin{equation} {\cal S}_{{\rm BG}}=-k_B \sum _i p_i \ln p_i, \end{equation} \tag{ 1 }$$

where k_B is the Boltzmann constant and p_i is the probability of the ith microstate. For two independent systems A and B, the probability of the system A+B being in a state i+j, with i a microstate of A and j a microstate of B, is

$$\begin{equation} p^{A+B}_{i+j} = p^A_i \cdot p^B _j. \end{equation} \tag{ 2 }$$

Therefore, the BG entropy satisfies the relation

$$\begin{equation} {\cal S}_{A+B}={\cal S}_A +{\cal S}_B. \end{equation} \tag{ 3 }$$

Thus, the entropy based on the BG statistic is an extensive quantity.

In the non-extensive statistics (Tsallis 1988), one replaces the traditional entropy by the following one:

$$\begin{equation} {\cal S}_q=k_B {1-\sum _i p^q_i \over q-1}, \end{equation} \tag{ 4 }$$

where q is a real number. For q = 1, ${\cal S}_q= {\cal S}_{{\rm BG}}$. Thus, the Tsallis statistics is a natural generalization of the BG entropy.

Now it follows that

$$\begin{equation} {\cal S}_q(A+B)={\cal S}_q(A) +{\cal S}_q(B)+{(1-q)\over k_B}{\cal S}_q(A){\cal S}_q(B). \end{equation} \tag{ 5 }$$

The variable q is a measure of the non-extensivity. Tsallis has shown that a formalism of statistical mechanics can be consistently developed in terms of this generalized entropy (Tsallis 2009).

A consequence of the non-extensive formalism is that the distribution function which maximizes S_q is non-Maxwellian (Silva et al. 1998; Lima et al. 2000; Muñoz 2006). For q = 1, the Maxwell distribution function is reproduced. But for q < 1, high-energy states are more probable than in the extensive case. On the other hand, for q > 1 high-energy states are less probable than in the extensive case, and there is a cutoff beyond which no states exist.

2.2. Maxwellian Distribution

In stars, the thermonuclear reaction rate with a Maxwellian distribution is given by (Fowler et al. 1967)

$$\begin{eqnarray} R_{ij}&=&\frac{N_{i}N_{j}}{1+\delta _{ij}} \langle \sigma v \rangle = \frac{N_{i}N_{j}}{1+\delta _{ij}}\left(\frac{8}{\pi \mu }\right)^{{1}/{2}}\left(\frac{1}{k_BT}\right)^{{3}/{2}}\nonumber \\ && \times \int _{0}^{\infty }dES(E) \exp \left[-\left(\frac{E}{k_BT}+2\pi \eta (E)\right)\right], \end{eqnarray} \tag{ 6 }$$

where σ is the fusion cross section, v is the relative velocity of the ij-pair, N_i is the number of nuclei of species i, μ is the reduced mass of i+j, T is the temperature, S(E) is the astrophysical S-factor, and η = Z_iZ_je²/ℏv is the Sommerfeld parameter, with Z_i the ith nuclide charge and E = μv²/2 is the relative energy of i+j.

The energy dependence of the reaction cross sections is usually expressed in terms of the equation

$$\begin{equation} \sigma (E) = {S(E)\over E}\exp \left[-{2\pi \eta (E)}\right]. \end{equation} \tag{ 7 }$$

We write $2\pi \eta =b/\sqrt{E}$, where

$$\begin{equation} b=0.9898 Z_iZ_j\sqrt{A} \ {\rm MeV}^{1/2}, \end{equation} \tag{ 8 }$$

where A is the reduced mass in amu. The factor 1 + δ_ij in the denominator of Equation (6) corrects for the double-counting when i = j. The S-factor has a relatively weak dependence on the energy E, except when it is close to a resonance, where it is strongly peaked.

2.3. Non-Maxwellian Distribution

The non-extensive description of the Maxwell–Boltzmann distribution corresponds to the substitution f(E) → f_q(E), where (Tsallis 2009)

$$\begin{eqnarray} &&f_q(E)=\left[1-\frac{q-1}{k_BT}E\right]^{{1}/{q-1}}\nonumber \\ &&\quad\quad\stackrel{q\rightarrow 1}{\longrightarrow }\exp \left(-\frac{E}{k_BT}\right), \quad 0<E<\infty. \end{eqnarray} \tag{ 9 }$$

If q − 1 < 0, Equation (9) is real for any value of E ⩾ 0. However, if q − 1 > 0, f(E) is real only if the quantity in square brackets is positive. This means that

$$\begin{eqnarray} && 0 \le E \le {k_BT\over q-1}, \ \ \ \ {\rm if} \ \ q\ge 1 \nonumber \\ &&\quad\quad\quad\quad 0\le E, \ \ \ \ {\rm if} \ \ q\le 1. \end{eqnarray} \tag{ 10 }$$

Thus, in the interval 0 < q < 1 one has 0 < E < ∞ and for 1 < q < ∞ one has 0 < E < E_max = k_BT/(q − 1).

With this new statistics, the reaction rate becomes

$$\begin{equation} R_{ij}= \frac{N_{i}N_{j}}{1+\delta _{ij}} I_q, \end{equation} \tag{ 11 }$$

and the rate integral, I_q, is given by

$$\begin{equation} I_{q}=\int _{0}^{E_{{\rm max}}}{\it dES}(E){\cal M}_{q}(E,T), \end{equation} \tag{ 12 }$$

where the "modified" Gamow energy distribution

$$\begin{eqnarray} {\cal M}_{q}(E,T)&=&{\cal A}(q,T)\left(1-\frac{q-1}{k_BT}E\right)^{{1}/{q-1}}e^{-b/\sqrt{{E}}} \nonumber \\ &=&{\cal A}(q,T) \left(1-\frac{q-1}{0.08617T_9}E\right)^{{1}/{q-1}}\nonumber \\ &&\times \exp \left[-0.9898 Z_iZ_j\sqrt{A\over E}\right] \end{eqnarray} \tag{ 13 }$$

is the non-extensive Maxwell velocity distribution, E_max = ∞ for 0 < q < 1 and E_max = k_BT/(1 − q) for 1 < q < ∞, and E is in MeV units. ${\cal A} (q,T)$ is a normalization constant which depends on the temperature and the non-extensive parameter q.

2.4. Non-Maxwellian Distribution for Relative Velocities

It is worthwhile to note that if the one-particle energy distribution is Maxwellian, it does not necessarily imply that the relative velocity distribution is also Maxwellian. Additional assumptions about correlations between particles are necessary to deduce that the relative-velocity distribution, which is the relevant quantity for rate calculations, is also Maxwellian. This has been discussed in detail in Lissia & Quarati (2005), Tsallis (1988), and Kaniadakis et al. (2005) where non-Maxwellian distributions, such as in Equation (13), were shown to arise from non-extensive statistics.

Here we show that if the non-Maxwellian particle velocity distribution is given by Equation (9), then a two-particle relative can be modified to account for the center of mass recoil. Calling the kinetic energy of a particle E_i, this distribution is given by

$$\begin{equation} f_q^{(i)} = \left(1 - \frac{q - 1}{kT}E_i\right)^{{1}/{q - 1}} \rightarrow \exp {\left[-\left(\frac{E_i}{kT}\right)\right]}. \end{equation} \tag{ 14 }$$

The two-particle energy distribution is $f_q^{(1)}f_q^{(2)}$. We now exponentiate the Tsallis distribution

$$\begin{equation} f_q^{(i)} = \exp \left\lbrace \frac{1}{q - 1} \left[\ln \left(1 - \frac{q - 1}{kT} E_i\right)\right]\right\rbrace \end{equation} \tag{ 15 }$$

and the product $f_q^{(12)}=f_q^{(1)}f_q^{(2)}$ reduces to

$$\begin{equation} f_q^{(12)} = \exp \left\lbrace \frac{1}{q - 1} \left[\ln \left(1 - \frac{q - 1}{kT} E_1\right)\left(1 - \frac{q - 1}{kT}E_2)\right)\right]\right\rbrace. \end{equation} \tag{ 16 }$$

Since $E_i = {m_i v_i^2}/{2}$, and thus, E₁ + E₂ = μv²/2 + MV²/2, where μ is the reduced mass of the two particles, M = m₁ + m₂, v is the relative velocity, and V is the center of mass velocity, the product inside the natural logarithm can be reduced to

$$\begin{eqnarray} && 1 - \frac{1 - q}{kT} \left(\frac{\mu v^2}{2} + \frac{M V^2}{2}\right) + \left(\frac{1 - q}{kT}\right)^2 \frac{\mu v^2}{2}\frac{M V^2}{2} \nonumber \\ && \quad =\left(1 - \frac{1 - q}{kT} \frac{\mu v^2}{2}\right)\left(1 - \frac{1 - q}{kT} \frac{M V^2}{2}\right). \end{eqnarray} \tag{ 17 }$$

Thus, the two-body distribution factorizes into a product of relative and center of mass parts:

$$\begin{eqnarray} &&f_q^{(12)}(v,V,T) = f_q^{{\rm (rel)}}(v,T) f_q^{{\rm (cm)}}(V,T), \end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray} f_q^{{\rm (rel)}}(v,T)&=&{\cal A}_{{\rm rel}}(q,T)\left(1 - \frac{1 - q}{kT} \frac{\mu v^2}{2}\right)^{{1}/{q - 1}}\nonumber \\ f_q^{{\rm (cm)}}(V,T)&=&{\cal A}_{{\rm cm}}(q,T)\left(1 - \frac{1 - q}{kT} \frac{M V^2}{2}\right)^{{1}/{q - 1}}, \end{eqnarray} \tag{ 19 }$$

with the normalization constants obtained from the condition

$$\begin{equation} \int d^3v d^3Vf_q^{(12)}(v,V,T) = 1. \end{equation} \tag{ 20 }$$

Because the distribution factorizes, the unit normalization can be achieved for the relative and center of mass distributions separately. The distribution needed in the reaction rate formula is, therefore,

$$\begin{equation} f_q(v,T) = \int d^3V f_q^{(12)}(v,V,T)= f_q^{{\rm (rel)}}(v,T), \end{equation} \tag{ 21 }$$

which attains the same form as for the absolute distribution.

In the limit q → 1 the two-particle distribution reduces to a Gaussian, with the last term in the left-hand side of Equation (17) dropping out,

$$\begin{equation} f_q^{(12)}= {\cal A}(q,T) \exp \left\lbrace \left[-\left(\frac{{\mu v^2}/{2} + {MV^2}/{2}}{kT}\right)\right]\right\rbrace , \end{equation} \tag{ 22 }$$

as expected.

2.5. Equilibrium with Electrons, Photons, and Neutrinos

One of the important questions regarding a plasma with particles (i.e., nuclei) described by the non-extensive statistics is how to generalize Fermi–Dirac (FD), Bose–Einstein (BE), and Tsallis statistics to become more unified statistics with the distribution for the particles. This has been studied in Buyukkiliç & Demirhan (2000), where it was shown that a similar non-extensive statistic for the distribution can be obtained for fermions and bosons and is given by

$$\begin{equation} n_q^{\pm } (E)= {1\over \left[1-(q-1){(E-\mu)\over kT}\right]^{{1}/{q-1}}\pm 1}, \end{equation} \tag{ 23 }$$

where μ is the chemical potential. This reproduces the Fermi distribution, n⁺, for q → 1 and the BE distribution for photons, n⁻, for μ = 0 and q → 1. Planck's law for the distribution of radiation is obtained by multiplying n⁻ in Equation (23) (with μ = 0) by ℏω²/(4π²c²), where E = ℏω. The number density of electrons can be obtained from n⁺ in Equation (23) with the proper phase factors depending if the electrons are non-relativistic or relativistic. Normalization factors ${\cal A}^\pm (q,T)$ also need to be introduced, as before.

The electron density during the early universe varies strongly with the temperature. At T₉ = 10, the electron density is about 10³² cm⁻³, much larger than the electron number density at the center of the Sun, $n_e^{{\rm sun}}\sim 10^{26}$ cm⁻³. The large electron density is due to the e⁺e⁻ production by the abundant photons during the BBN. However, the large electron densities do not influence the nuclear reactions during the BBN. In fact, the enhancement of the nuclear reaction rates due to electron screening was shown to be very small (Wang et al. 2011). The electron Fermi energy for these densities is also much smaller than kT for the energy relevant for BBN, so that one can also use μ = 0 in Equation (23) for n⁺.

In Figure 2, we plot the relative difference ratio $(n_q^{\pm }-n^\pm)/n^\pm$ between non-extensive, n_q, and extensive, n = n_{q → 1}, statistics. For both distributions, we use μ = 0. The solid curves are for FD statistics, n⁺, and dashed curves are for BE statistics, n⁻. We show results for q = 2 and q = 0.5, with T₉ = 10. One sees that the non-extensive distributions are enhanced for q > 1 and suppressed for q < 1, as compared to the respective FD and BE quantum distributions. The deviations for the FD and BE statistics also grow larger with the energy. For example, the non-extensive electron distribution is roughly a factor of two larger than the usual FD distribution at E_e = 1 MeV, at T₉ = 10.

Zoom In Zoom Out Reset image size

While we do not obtain numerical results with modified DG and BE distributions here, it can be expected that these generalizations will have a strong influence on the freezeout temperature and the neutron to proton, n/p, ratio. A numerical study of this problem may be presented in another paper. The freezeout temperature occurs when the rate, Γ ∼ <σv >, for weak reaction ν_e + n → p + e⁻ becomes slower than the expansion rate of the universe. Because during the BBN the densities of all particles, including neutrinos, are low compared to kT, the chemical potential μ can be set to zero in the calculation of all reaction rates. The adoption of non-extensive quantum distributions such as in Equation (23) will lead to the same powers of the temperature as those predicted by the FD distribution and the BE distribution. For example, Planck's law for the total blackbody is U∝T⁴, being form invariant with respect to non-extensive entropic index q which determines the degree of non-extensivity (Buyukkiliç et al. 2002). This result means that the weak decay reaction rates do not depend on the non-extensive parameter q. The freezeout temperature and n/p ratio remain the same as before.

More detailed studies have indicated that Planck's law of blackbody radiation and other thermodynamical quantities arising from non-extensive quantum statistics can yield different powers of temperature than for the non-extensive case (Aragao-Rego et al. 2003; Nauenberg 2003; Tsallis 2004; Chamati et al. 2004). If that is the case, then a study of the influence of non-extensive statistics on the weak decay rates and electromagnetic processes during BBN is worth pursuing.

2.6. Thermodynamical Equilibrium

Based on the abundant literature on non-extensive statistics (see, e.g., Lissia & Quarati 2005; Haubold & Kumar 2008; Degl'Innocenti et al. 1998; Coraddu et al. 1999; Tsallis 1988), we do not expect that the non-extensive parameter q differs appreciably from the unity value. However, in order to study the influence of a non-Maxwellian distribution on BBN we will explore values of q rather different than unity, namely, q = 0.5 and q = 2. This will allow us to pursue a better understanding of the nature of the physics involved in the departure from the BG statistics. In Figure 1, we plot the Gamow energy distributions of deuterons relevant for the reaction ²H(d, p)³H at T₉ = 1. The solid line, for q = 1, corresponds to the use of the Maxwell–Boltzmann distribution. Also shown are results for non-extensive distributions for q = 0.5 (dotted line) and q = 2 (dashed line). One observes that for q < 1, higher kinetic energies are more accessible than in the extensive case (q = 1). For q > 1, high energies are less probable than in the extensive case, and there is a cutoff beyond which no kinetic energy is reached. In the example shown in the figure for q = 2, the cutoff occurs at 0.086 MeV, or 86 keV.

We will explore the modifications of the BBN elemental abundances due to a variation of the non-extensive statistics parameter q. We will express our reaction rates in the form N_A〈σv〉 (in units of cm³ mol⁻¹ s⁻¹), where N_A is the Avogadro number and 〈σv〉 involves the integral in Equation (6) with the Maxwell distribution f(E) replaced by Equation (9). First, we show how the reaction rates are modified for q ≠ 1.

In Figure 3, we show the S-factor for the reaction ²H(d, p)³H as a function of the relative energy E. Also shown is the dependence on T₉ (temperature in units of 10⁹ K) for the effective Gamow energy:

$$\begin{equation} E=E_0=0.122\left(Z_i^2Z^2_jA\right)^{1/3}T_9^{2/3} \ {\rm MeV}, \end{equation} \tag{ 24 }$$

where A is the reduced mass in amu. The data are from Schulte et al. (1972), Krauss et al. (1987), Brown & Jarmie (1990), Bosch & Hale (1992), and Greife et al. (1995) and the solid curve is a chi-square polynomial function fit to the data.

Zoom In Zoom Out Reset image size

Using the chi-square polynomial fit obtained to fit the data presented in Figure 3, we show in Figure 4 the reaction rates for ²H(d, p)³H as a function of the temperature T₉ for two different values of the non-extensive parameter q. The integrals in Equation (12) are performed numerically. For charge particles, a good accuracy (within 0.1%) is reached using the integration limits between E₀ − 5ΔE and E₀ + 5ΔE, where ΔE is given by Equation (25) below. The rates are expressed in terms of the natural logarithm of N_A〈σv〉 (in units of cm³ mol⁻¹ s⁻¹). The solid curve corresponds to the usual Maxwell–Boltzmann distribution, i.e., q = 1. The dashed and dotted curves are obtained for q = 2 and q = 0.5, respectively. In both cases, we see deviations from the Maxwellian rate. For q > 1, the deviations are rather large and the tendency is an overall suppression of the reaction rates, especially at low temperatures. This effect arises from the non-Maxwellian energy cutoff which for this reaction occurs at 0.086T₉ MeV and which prevents a great number of reactions from occurring at higher energies.

Zoom In Zoom Out Reset image size

For q  <  1, the nearly similar result as with the Maxwell–Boltzmann distribution is due to a competition between suppression in reaction rates at low energies and their enhancement at high energies. The relevant range of energies is set by the Gamow energy which for a Maxwellian distribution is given by Equation (24) and by the energy window,

$$\begin{equation} \Delta E=0.2368\left(Z_i^2Z^2_jA\right)^{1/6}T_9^{5/6} \ {\rm MeV}, \end{equation} \tag{ 25 }$$

which for the reaction ²H(d, p)³H amounts to $0.2368T_9^{5/6}$ MeV. This explains why, at T₉ = 1, the range of relevant energies for the calculation of the reaction rate is shown by the solid curve in Figure 1. For q < 1, the Gamow window ΔE is larger and there is as much a contribution from the suppression of reaction rates at low energies compared to the Maxwell–Boltzmann distribution as there is a corresponding enhancement at higher energies. This explains the nearly equal results shown in Figure 4 for q = 1 and q < 1. This finding applies to all charged particle reaction rates, except for those when the S-factor has a strong dependence on energy at, and around, E = E₀. But no such behavior exists for the most important charged-induced reactions in the BBN (neutron-induced reactions will be discussed separately).

The findings described above for the reaction ²H(d, p)³H are not specific but apply to all charged particles of relevance to the BBN. We demonstrate this with one more example: the ⁷Li(p, α)⁴He reaction, responsible for ⁷Li destruction. In Figure 5, we show the S-factor for this reaction as a function of the relative energy E. One sees prominent resonances at higher energies. Also shown in the figure is the dependence of the reaction on T₉. The data are from Freeman et al. (1958), Cassagnou et al. (1962, 1963), Mani et al. (1964), Fiedler & Kunze (1967), Lerner & Marion (1969), Spinka et al. (1971), Rolfs & Kavanagh (1986), Werby (1973), and Engstler et al. (1992a, 1992b) and the solid curve is a chi-square function fit to the data using a sum of polynomials plus Breit–Wigner functions.

Zoom In Zoom Out Reset image size

Using the chi-square function fit obtained to fit the data presented in Figure 5, we show in Figure 6 the reaction rates for ⁷Li(p,α)⁴He as a function of the temperature T₉ for two different values of the non-extensive parameter q. The rates are given in terms of the natural logarithm of N_A〈σv〉 (in units of cm³ mol⁻¹ s⁻¹). The solid curve corresponds to the usual Maxwell–Boltzmann distribution, i.e., q = 1. The dashed and dotted curves are obtained for q = 2 and q = 0.5, respectively. As with the reaction presented in Figure 4, in both cases we see deviations from the Maxwellian rate. But, as before, for q = 2 the deviations are larger and the tendency is a strong suppression of the reaction rates as the temperature decreases. It is interesting to note that the non-Maxwellian rates for q = 0.5 are more sensitive to the resonances than for q > 1. This is because, as seen in Figure 1, for q < 1 the velocity distribution is spread to considerably larger values of energies, being therefore more sensitive to the location of high-energy resonances.

Zoom In Zoom Out Reset image size

We now turn to neutron-induced reactions, which are only a few cases of high relevance for the BBN, notably the p(n, γ)d, ³He(n, p)t, and ⁷Be(n, p)⁷Li reactions. For neutron-induced reactions, the cross section at low energies is usually proportional to 1/v, where $v=\sqrt{2mE}/\hbar$ is the neutron velocity. Thus, it is sometimes appropriate to rewrite Equation (7) as

$$\begin{equation} \sigma (E) = {S(E)\over E}= {R(E) \over \sqrt{E}}, \end{equation} \tag{ 26 }$$

where R(E) is a slowly varying function of energy similar to an S-factor. The distribution function within the reaction rate integral (12) is also rewritten as

$$\begin{eqnarray} {\cal M}_{q}(E,T)&={\cal A} (q,T)f_q(E)=&{\cal A} (q,T)\left(1-\frac{q-1}{k_BT}E\right)^{{1}/{q-1}}. \nonumber \\ \end{eqnarray} \tag{ 27 }$$

The absence of the tunneling factor $\exp (-b/\sqrt{E})$ in Equation (27) inhibits the dependence of the reaction rates on the non-extensive parameter q.

In Figure 7, we plot the kinetic energy distributions of nucleons relevant for the reaction p(n, γ)d at T₉ = 0.1 (upper panel) and T₉ = 10 (lower panel). The solid line, for q = 1, corresponds to the usual Boltzmann distribution. Also shown are results for the non-extensive distributions for q = 0.5 (dotted line) and q = 2 (dashed line). One observes that, as for the charged particles case, with q < 1 higher kinetic energies are more probable than in the extensive case (q = 1). With q > 1, high energies are less accessible than in the extensive case, and there is a cutoff beyond which no kinetic energy is reached. A noticeable difference from the case of charged particles is the absence of the Coulomb barrier and a corresponding lack of suppression of the reaction rates at low energies. As the temperature increases, the relative difference between the Maxwell–Boltzmann and the non-Maxwellian distributions decreases appreciably. This will lead to a rather distinctive pattern of the reaction rates for charged-induced reactions as compared to neutron-induced reactions.

Zoom In Zoom Out Reset image size

For neutron-induced reactions, a good accuracy (within 0.1%) for the numerical calculation of the reaction rates with Equation (12) is reached using the integration limits between E = 0 and E = 20 k_BT. As an example we will now consider the reaction ⁷Be(n, p)⁷Li. The energy dependence of $R(E)=S(E)\sqrt{E}$ for this reaction is shown in Figure 8. The experimental data were collected from Gibbons & Macklin (1959), Borchers & Poppe (1963), Sekharan et al. (1976), Poppe et al. (1976), and Koehler et al. (1988).

Zoom In Zoom Out Reset image size

Using the chi-square fit with a sum of polynomials and Breit–Wigner functions obtained to reproduce the data in Figure 8, we show in Figure 9 the reaction rates for ⁷Be(n, p)⁷Li as a function of the temperature T₉ for different values of the non-extensive parameter q. The rates are given in terms of the natural logarithm of N_A〈σv〉 (in units of cm³ mol⁻¹ s⁻¹). The solid curve corresponds to the usual Boltzmann distribution, i.e., q = 1. The dashed and dotted curves are obtained for q = 2 and q = 0.5, respectively. In contrast to reactions induced by charged particles, we now see strong deviations from the Maxwellian rate both for q > 1 and q < 1. For q < 1, the deviations are larger at small temperatures and decrease as the energy increases, tending asymptotically to the Maxwellian rate at large temperatures. This behavior can be understood from Figure 7 (for ⁷Be(n, p)⁷Li the results are nearly the same as in Figure 7). At small temperatures, e.g., T₉ = 0.1, the distribution for q = 0.5 is strongly enhanced at large energies and the tendency is that the reaction rates increase at low temperatures. This enhancement disappears as the temperature increases (lower panel of Figure 7). For q = 2 the reaction rate is suppressed, although not as much as for the charged-induced reactions, the reason being compensation by an increase because of normalization at low energies.

Zoom In Zoom Out Reset image size

Having discussed the dependence of the reaction rates on the non-extensive parameter q for a few standard reactions, we now consider the implications of the non-extensive statistics to the predictions of the BBN. It is clear from the results presented above that an appreciable impact on the abundances of light elements will arise.

The BBN is sensitive to certain parameters, including the baryon-to-photon ratio, number of neutrino families, and the neutron decay lifetime. We use the values η = 6.19 × 10⁻¹⁰, N_ν = 3, and τ_n = 878.5 s for the baryon–photon ratio, number of neutrino families, and neutron decay lifetime, respectively. Our BBN abundances were calculated with a modified version of the standard BBN code derived from Wagoner et al. (1967), L. Kawano (1988, unpublished),⁵ and Kawano (1992).

Although BBN nucleosynthesis can involve reactions up to the CNO cycle (Coc et al. 2012), the most important reactions which can significantly affect the predictions of the abundances of the light elements [⁴He, D, ³He, ⁷Li] are n-decay, p(n, γ)d, d(p,γ)³He, d(d, n)³He, d(d, p)t, ³He(n, p)t, t(d, n)⁴He, ³He(d, p)⁴He, ³He(α, γ)⁷Be, t(α, γ)⁷Li, ⁷Be(n, p)⁷Li, and ⁷Li(p, α)⁴He. Except for these reactions, we have used the reaction rates needed for the remaining reactions from a compilation by NACRE (Angulo et al. 1999) and that reported in Descouvemont et al. (2004). For the 11 reactions mentioned above, we have collected data from Smith et al. (1993), Angulo et al. (1999), and Descouvemont et al. (2004) and references mentioned therein (data for n(p, γ, d) reaction were taken from the online ENDF database⁶; see also Cyburt 2004; Ando et al. 2006), fitted the S-factors with a sum of polynomials and Breit–Wigner functions, and calculated the reaction rates for Maxwellian and non-Maxwellian distributions.

4.1. Elemental Abundances

In Figure 10, we show the calculated deuterium abundance. The solid curve is the result with the standard Maxwell distributions for the reaction rates. Using the non-extensive distributions yields the dotted line for q = 0.5 and the dashed line for q = 2. It is interesting to observe that the deuterium abundances are only moderately modified due to the use of the non-extensive statistics for q = 0.5. Up to temperatures of the order of T₉ = 1, the abundance for D/H tends to agree for the extensive and non-extensive statistics. This is due to the fact that any deuterium that is formed is immediately destroyed (a situation known as the deuterium bottleneck). But, as the temperature decreases, the reaction rates for the p(n, γ)d reaction are considerably enhanced for q = 2 (see Figure 7), and perhaps more importantly, they are strongly suppressed for all other reactions involving deuterium destruction, as clearly seen in Figure 4. This creates an unexpected overabundance of deuterons for the non-extensive statistics with q = 2. The deuterium, a very fragile isotope, is easily destroyed after the BBN and astrated. Its primordial abundance is determined from observations of interstellar clouds at high redshift, on the line of sight of distant quasars. These observations are scarce but allow to set an average value of ${\rm D/H} = 2.82^{ + 0.20}_{ - 0.19} \times 10^{-5}$ (Pettini et al. 2008). The predictions for the D/H ratio with the q = 2 statistics (D/H = 5.70 × 10⁻³) are about a factor of 200 larger than those from the standard BBN model, clearly in disagreement with the observation.

Zoom In Zoom Out Reset image size

A much more stringent constraint for elemental abundances is given by ⁴He, which observations set at about ⁴He/H ≡ Y_p = 0.2561 ± 0.0108 (Boesgaard & Steigman 1985; Aver et al. 2010; Izotov & Thuan 2010). The ⁴He abundance generated from our BBN calculation is plotted in Figure 11. The solid curve is the result obtained with the standard Maxwell distribution for the reaction rates. Using the non-extensive distributions yields the dotted line for q = 0.5 and the dashed line for q = 2. Again, the predicted abundances for q = 2 deviate substantially from standard BBN results. This time only about half of ⁴He is produced with the use of a non-extensive statistics with q = 2. The reason for this is the suppression of the reaction rates for the formation of ⁴He with q = 2 through the charged particle reactions, t(d, n)⁴He and ³He(d, p)⁴He.

Zoom In Zoom Out Reset image size

A strong impact of using non-extensive statistics for both q = 0.5 and q = 2 values of the non-extensive parameter is seen in Figure 12 for the ³He abundance. While for q = 2 there is an overshooting in the production of ³He, for q = 0.5 one finds a smaller value than the one predicted by the standard BBN. This is due to the distinct results for the destruction of ³He through the reaction ³He(n, p)t, which is enhanced for q = 0.5 and suppressed for q = 2, in the same way as it happens for the reaction ⁷Be(n, p)⁷Li, shown in Figure 8. ³He is both produced and destroyed in stars and its abundance is still subject to large uncertainties, ³He/H = (1.1 ± 0.2) × 10⁻⁵ (Bania et al. 2002; Vangioni-Flam et al. 2003).

Zoom In Zoom Out Reset image size

Non-extensive statistics for both q = 0.5 and q = 2 values also alter substantially the ⁷Li abundance, as shown in Figure 13. For both values of the non-extensive parameter, q = 2 and q = 0.5, there is an overshooting in the production of ⁷Li. The increase in the ⁷Li abundance is more accentuated for q = 2. The lithium problem is associated with a smaller value of the observed ⁷Li abundance as compared to the predictions of BBN. There have been many attempts to solve this problem by testing all kinds of modifications of the parameters of the BBN or the physics behind it (a sample of this literature is found in Fields 2011; Cyburt 2004; Nollett & Burles 2000; Cyburt et al. 2008; Cyburt & Pospelov 2009; Boyd et al. 2010; Wang et al. 2011; Kirsebom & Davids 2011; Iocco & Pato 2012 and references therein). In the present case, the use of a non-Maxwellian velocity distribution seems to worsen this scenario. A recent analysis yields the observational value of ⁷Li/H = (1.58 ± 0.31) × 10⁻¹⁰ (Sbordone et al. 2010).

Zoom In Zoom Out Reset image size

4.2. Sensitivity Study

We have calculated a window of opportunity for the non-extensive parameter q with which one can reproduce the observed abundance of light elements. We chose the data for the abundances, Y_i, of ⁴He/H, D/H, ³He/H, and ⁷Li as reference. We then applied the ordinary χ² statistics, defined by the minimization of

$$\begin{equation} \chi ^2 = \sum _i {[Y_i(q)-Y_i({\rm obs})]^2\over \sigma _i}, \end{equation} \tag{ 28 }$$

where Y_i(q) are the abundances obtained with the non-extensive statistics with parameter q, Y_i(obs) are the observed abundances, and σ_i are the errors for each datum, and the sum is over all data mentioned in Table 1. From this chi-square fit we conclude that $q=1^{+0.05}_{-0.12}$ is compatible with observations.

Table 1. Predictions of the BBN (with η_WMAP = 6.2 × 10⁻¹⁰) with Maxwellian and Non-Maxwellian Distributions Compared with Observations

	Maxwell	Non-Max.	Non-Max.	Observation
	BBN	q = 0.5	q = 2
4He/H	0.249	0.243	0.141	0.2561 ± 0.0108
D/H	2.62	3.31	570	$2.82^{+0.20}_{-0.19}$(× 10−5)
3He/H	0.98	0.91	69.1	(1.1 ± 0.2)(× 10−5)
7Li/H	4.39	6.89	356.	(1.58 ± 0.31)(× 10−10)

Notes. All numbers have the same power of 10 as in the last column.

Download table as:  ASCIITypeset image

No attempt has been made to determine which element dominates the constraint on q. This might be important for a detailed study of the elemental abundance influence from nuclear physics inputs, namely, the uncertainty of the reaction cross sections. A study along these lines might be carried in a similar fashion as described in Nollett & Burles (2000). Weights on the reliability of observational data should also be considered for a more detailed analysis. For example, constraints arising from ³He/H abundance may not be considered trustworthy because of uncertain galactic chemical evolution. On the other hand, a constraint from the observation of ³He/D is more robust, and so on. Based on our discussion in Section 2, it is more important to determine how the non-extensive statistics can modify more stringent conditions during the big bang, such as the modification of weak decay rates and its influence on the n/p ratio which strongly affects the ⁴He abundance.

In Table 1, we present results for the predictions of the BBN with Maxwellian and non-Maxwellian distributions. The predictions are compared with data from observations reported in the literature. It is evident that the results obtained with the non-extensive statistics strongly disagree with the data. The overabundance of ⁷Li compared to observation gets worse if q > 1. The three light elements D, ⁴He, and ⁷Li constrain the primordial abundances rather well. For all these abundances, a non-extensive statistics with q > 1 leads to a greater discrepancy with the experimental data.

Except for the case of ³He the use of non-extensive statistics with q < 0.5 does not rule out its validity when the non-Maxwellian BBN results are compared to observations. ³He is at present only accessible in our Galaxy's interstellar medium. This means that it cannot be measured at low metallicity, a requirement to make a fair comparison to the primordial generation of light elements. This also means that the primordial ³He abundance cannot be determined reliably. The result presented for the ³He abundance in Table 1 is quoted from Bania et al. (2002). Note that our analysis does not include the changes that the non-extensive statistics would bring to the n/p conversion rates. The electron distributions would also be expected to change accordingly. This would change the freezeout temperature and a corresponding influence on the ⁴He abundance.

We conclude that it does not seem possible to change the Maxwell–Boltzmann statistics to reproduce the observed abundance of light elements in the universe without destroying many other successful predictions of BBN. A chi-square fit of our calculations with the observations of elemental abundance concludes that the non-extensive parameter is constrained to $q=1^{+0.05}_{-0.12}$. This means that, should a non-Maxwellian distribution due to the use of the Tsallis non-extensive statistics be confirmed (with a sizable deviation from q = 1), our understanding of the cosmic evolution of the universe would have to be significantly changed.

We acknowledge beneficial discussions with Alain Coc and Richard Cyburt. This work was partially supported by the US-DOE grants DE-FG02-08ER41533 and DE-FG02-10ER41706, and by the Brazilian agencies, CNPq and FAPESP. One of the authors (C.B.) acknowledges the Helmholtz International Center for FAIR (HIC for FAIR) for supporting his visit to the GSI Helmholtzzentrum für Schwerionenforschung, where much of this work was done.