Neutrino energy transport in weak decoupling and big bang nucleosynthesis

Neutrino energy transport in weak decoupling and big bang nucleosynthesis

E. Grohs    G. M. Fuller    C. T. Kishimoto    M. W. Paris    A. Vlasenko Department of Physics, University of California, San Diego, La Jolla, California 92093, USA Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA Department of Physics and Biophysics, University of San Diego, San Diego, California 92110, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA Department of Physics, North Carolina State University, Raleigh, North Carolina 27695, USA
July 9, 2019

We calculate the evolution of the early universe through the epochs of weak decoupling, weak freeze-out and big bang nucleosynthesis (BBN) by simultaneously coupling a full strong, electromagnetic, and weak nuclear reaction network with a multi-energy group Boltzmann neutrino energy transport scheme. The modular structure of our code provides the ability to dissect the relative contributions of each process responsible for evolving the dynamics of the early universe in the absence of neutrino flavor oscillations. Such an approach allows a detailed accounting of the evolution of the , , , , , energy distribution functions alongside and self-consistently with the nuclear reactions and entropy/heat generation and flow between the neutrino and photon/electron/positron/baryon plasma components. This calculation reveals nonlinear feedback in the time evolution of neutrino distribution functions and plasma thermodynamic conditions (e.g., electron-positron pair densities), with implications for: the phasing between scale factor and plasma temperature; the neutron-to-proton ratio; light-element abundance histories; and the cosmological parameter . We find that our approach of following the time development of neutrino spectral distortions and concomitant entropy production and extraction from the plasma results in changes in the computed value of the BBN deuterium yield. For example, for particular implementations of quantum corrections in plasma thermodynamics, our calculations show a increase in deuterium. These changes are potentially significant in the context of anticipated improvements in observational and nuclear physics uncertainties.

preprint: LA-UR-15-29163

I Introduction

In this paper we concurrently solve for the evolution of the neutrino and matter/radiation components in the early universe. A key result of this work is that there is, in fact, nonlinear feedback between these components during the time when the neutrinos go from thermally and chemically coupled with the plasma of photons/electrons/positrons/baryons, to completely decoupled and free streaming. This feedback can be important for high precision calculations of the primordial light element abundances emerging from big bang nucleosynthesis (BBN). The work we describe here builds on the many previous studies of the evolution of the neutrino energy distribution functions in the early universe (see Refs. Dicus et al. (1982); Cambier et al. (1982); Herrera and Hacyan (1989); Rana and Mitra (1991); Dodelson and Turner (1992); Hannestad and Madsen (1995); Dolgov et al. (1997); Gnedin and Gnedin (1998); Dolgov et al. (1999); Lopez and Turner (1999); Esposito et al. (2000); Mangano et al. (2002, 2005); Birrell et al. (2015) and Appendix A). Higher precision in theoretical calculations of neutrino transport and nucleosynthesis in the early universe is warranted by recent and anticipated improvement in the precision of cosmological observations.

The advent of high precision cosmological observations will demand a deeper understanding and higher precision in modeling the known microphysics of the standard model relevant during early universe through the neutrino weak decoupling (or simply “weak decoupling”) and BBN epochs. For example, future cosmic microwave background (CMB) polarization experiments promise increased sensitivity to issues closely associated with relic neutrino energy distribution functions such as the “sum of the light neutrino masses” and measures of the radiation energy density Abazajian et al. (2015). Additionally, the advent of extremely large optical telescopes, with adaptive optics, can improve the precision in primordial abundance determinations Silva et al. (2007); McCarthy and Bernstein (2014); Hook (ed.) (2005).

Leveraging the increased observational precision to achieve better probes of and constraints on beyond standard model (BSM) physics will demand higher precision in simulation of standard model physics. Many possible BSM scenarios (e.g. sterile neutrinos, light scalars, out-of-equilibrium particle decay, etc.) could affect weak decoupling and, hence, nucleosynthesis in subtle but potentially measurable ways. Accurate and self-consistent treatments of the standard nuclear and particle physics furthers the objective of a clear interpretation of potential BSM issues (see Refs. Fuller et al. (2011); Nollett and Steigman (2015); Vassh et al. (2015)).

Neutrino kinetics affect the neutrino distributions and primordial nuclide abundances in the early universe in three principal respects. First, the transfer of entropy from the photon/electron/positron plasma to the neutrino seas cools the plasma temperature relative to the case of no transport. The cooler temperature alters the ratio of comoving to plasma energy scales from the canonical value Kolb and Turner (1990); Dodelson (2003); Weinberg (2008).

The second out-of-equilibrium effect is the distortion of the thermal Fermi-Dirac (FD) spectrum of high energy neutrinos. Upscattering of low energy neutrinos and the production of neutrino-antineutrino pairs contribute to this distortion through a variety of mechanisms. An important consequence of this mechanism is the effect that the high-energy distortion has on the neutron-to-proton ratio (). A running theme throughout the present study is that such changes induced by the distortion of the neutrino distributions away from equilibrium have effects that must be calculated concurrently with the evolution of the nuclide abundances. In this way, we reveal nonlinearities in feedback mechanisms between the neutrino transport and the thermodynamics of the plasma. These changes to the temperature evolution have an effect on relative changes in the nuclide abundances through the reaction rates and the sensitive dependence of, for example, Coulomb barriers on them.

The third out-of-equilibrium effect is entropy production. The Boltzmann H theorem implies that the entropy of a closed system is a non-decreasing function of time. In this paper, we investigate the conventional assumption Kolb and Turner (1990); Weinberg (2008) of comoving entropy conservation. We find that there is a small change in the total entropy of the universe due to the non-equilibrium kinetics of the neutrinos, which generates entropy. In essence, out-of-equilibrium neutrino energy transport and associated entropy flow changes the phasing between scale factor and plasma temperature evolution.

A common feature of past works is that the effect of these transport/entropy issues on the primordial abundances is small, typically on the order of 0.05% for helium-4 and lithium (in particular, see Ref. Dolgov et al. (1997), hereafter DHS). Our work shows that the magnitude of these effects can be significantly larger, depending on assumed microphysics.

The present work employs a non-perturbative method to calculate the evolution of active neutrino occupation probabilities for flavor . Homogeneity and isotropy has been assumed to restrict the dependence of the to only the magnitude of the three-momentum and the comoving time . The evolution is computed in the presence of two-body to two-body () collisions, the rates of which are given by the collision integrals , where refers to the occupation probabilities of neutrinos, antineutrinos, and charged leptons. These are functionals of the set of neutrino and antineutrino occupation probabilities and evolve, within the Boltzmann equation approach, as


where is the Hubble expansion rate at scale factor . We define the independent variable using the neutrino energy and the comoving temperature parameter . The comoving temperature parameter is not a physical temperature. It is simply an energy scale that redshifts like the energy of a massless particle in free fall with the expansion of the universe and is, in essence, a proxy for inverse scale factor. Therefore, we can write as a function of scale factor, where and are the plasma temperature and scale factor at an initial epoch of our choosing. For neutrinos in the range of plasma temperatures 3 MeV 10 keV, is equivalent to the commonly used quantity . Equation (1) can be cast in terms of as


The independent variable is chosen so that energy conservation takes the simple form for the scattering process .

The evaluation of the collision integral in Eq. (1) or (2) for the weak-interaction processes of interest is numerically intensive. However, the required integrations (described in detail in Sec.II.2.3 and Appendices B and C) are performed in parallel with the code burst (BBN/Unitary/Recombination/Self-consistent/Transport) in Fortran 90/95 under openmpi. We have developed a routine to evaluate the collision term for the Boltzmann equation in burst (using methods detailed in the Appendices) which reduce the number of required integrations to two. Numerical integration, effected under a combination of quadrature techniques (detailed in Sec.II), has been tested by ensuring conservation of lepton number; it is satisfied at the level of (see Sec.II.3.2).

The code has been developed to address the problem of weak-decoupling collision terms and for self-consistent coupling to nuclear reactions assuming that a Boltzmann equation treatment is sensible. The “embarrassingly parallel” structure of the problem allows for the simultaneous evaluation of the occupation probabilities for each energy, implying a nearly linear scaling of code performance with the number of cores. The present calculational approach is readily generalizable to treat the full neutrino quantum kinetic equations (QKEs) developed in Ref. Vlasenko et al. (2014) and therefore neutrino flavor oscillations (see Refs. Barbieri and Dolgov (1991); Akhmedov and Berezhiani (1992); Raffelt and Sigl (1993); Strack and Burrows (2005); Balantekin and Pehlivan (2007); Volpe et al. (2013); Balantekin and Fuller (2013); de Gouvêa and Shalgar (2013); Serreau and Volpe (2014); Cirigliano et al. (2015) for discussion on the QKEs). As mentioned, the present work neglects neutrino flavor oscillations. A detailed calculation that concurrently solves the neutrino QKE equations, incorporating both effects of flavor oscillations and energy transport, and the primoridal nucleosynthesis is required and currently underway. An example of the need for such a calculation is indicated by the high sensitivity of the ratio at weak freeze-out to the electron neutrino energy distrubtion (see Sec.III). One of the primary effects of flavor oscillation, whose subsequent effect on primordial nucleosynthesis is difficult to estimate in a self-consistent approach in the dynamic environment of the BBN-epoch of the early universe, is the suppression of the rate. This suppression occurs when an electron neutrino oscillates to either a or state, which do not convert . A detailed, self-consistent calculation will account for the phasings of various such mechanisms, which may be important at the level of precision anticipated for in the next generation of cosmological observations.

We emphasize that we couple neutrino-energy transport self-consistently and concurrently to evaluation of the neutron-to-proton rates and nucleosynthesis reaction network. At each time step in burst, the weak interaction neutron-proton conversion rates ( rates),


are determined using the evolved, non-equilibrium and spectra. The thermodynamics of the electromagnetic plasma is coupled to the neutrino seas to account for heat flow between the plasma and the neutrinos. Non-equilibrium effects generate entropy, increasing the total entropy of the plasma and the neutrinos, through a timelike entropy-current flux. Finally, we integrate the neutrino occupation probabilities to determine the energy density for calculating the Hubble expansion rate. In this way, self-consistency within the neutrino sector is maintained over approximately Hubble times. The overall architecture employed in burst differs from the approaches used in previous treatments (see Appendix A).

The nuclear reaction network employed in the current code is based on those of Refs. Wagoner (1969); Smith et al. (1993) as augmented in Ref. Boyd et al. (2010); details are discussed in Ref. Grohs et al. (2015). Ongoing work is focused on incorporating into the present approach a nuclear reaction network based on a reaction formalism that respects unitarity.

The outline of this work is as follows. In Sec. II, we present details of the transport code and weak-decoupling calculations. We investigate, in Sec. III, the contributions of the scattering processes to the out-of-equilibrium neutrino spectra. Section IV describes the evolution of the entropy during the weak-decoupling process. Section V discusses primordial nucleosynthesis resulting from the self-consistent coupling to the transport code. We conclude in Sec. VI. Appendix A contains a summary of the calculations of different groups. Appendices B and C describe the analytical derivations of the collision terms. We should emphasize that the current manuscript represents a preliminary step toward the objective of coupling neutrino kinetics to the nucleosynthesis reaction network. The proper treatment of neutrino flavor oscillations and possible coherent effects requires a quantum kinetic approach Vlasenko et al. (2014). Flavor oscillations have been estimatedMangano et al. (2002) to change the production of at the 20% level. The self-consistent approach that we consider here might be expected to enhance this change; a detailed calculation is required to estimate the actual effect. We detail further, ongoing efforts in this work in the conclusion, Sec. VI. Throughout this paper we use natural units where .

In this manuscript we have provided a pedagogical presentation of some familiar topics. This is done in the interest of giving a clear presentation of our work and in the hopes of making our analytical and numerical computations reproducible.

Table 1: Weak interaction processes relevant for neutrino weak decoupling. The left column labels the scattering, production, and annihilation processes in the middle column by an index . The right column gives the spin-averaged and summed square of the matrix element for process with the Fermi constant and symmetry factor divided out. Indices and in the middle column for processes , which describe neutrino and antineutrino scattering, are distinct. Processes with an antineutrino scattering on a charged lepton, correspond to the parity-conjugate reactions of . Since they have identical matrix elements to these they are not shown in the table, although their effect is explicitly accounted for in antineutrino energy transport. is unity for all processes except , where .

Ii Neutrino weak decoupling calculations

Neutrinos decouple from the plasma, roughly speaking, when typical rates of the weak processes, given in Table 1, fall below the Hubble rate:


where MeV is the Fermi constant and MeV. By numerically evolving the neutrino distributions for , and we find, however, that the neutrinos exchange entropy with the plasma until a temperature of nearly 100 keV, for many Hubble times beyond the estimate in Eq. (6) (see Fig. 9). This is in part explained by the fact that given the large entropy of the early universe, which is carried by both photons/electrons/positrons and neutrinos, a significant fraction of the neutrinos have energies larger than the temperature. This effect is enhanced by plasma particles scattering from the neutrinos, which preferentially up-scatter the neutrinos and distort the high-momentum tails of the neutrino distributions. In this section, we present the details of the numerical evaluation of the collision integrals, the solution of the Boltzmann equation, and performance statistics of the code, followed by details of the weak decoupling calculations.

ii.1 Weak interaction processes

We discuss the weak interactions relevant for neutrino weak decoupling here and their implementation in the collision integral in the Boltzmann equation, Eq. (1).

Expressions for the neutral and charged current weak interaction processes involving neutrinos, antineutrinos and the charged leptons of the plasma are given in Table 1. The table gives the squared amplitudes , where labels two-body processes that are important during neutrino weak decoupling Tubbs and Schramm (1975); Flowers and Sutherland (1976), averaged over initial spin states and summed over final spins. The initial state particle four-momenta in Table 1 are given particle numbers 1 and 2; final states are 3 and 4. That is:


where particle 1 is always a neutrino (or antineutrino). We label neutrino four-momenta as and charged lepton four-momenta as .

The are different for electron-flavor neutrinos compared to or -flavor neutrinos due to the charged-current interaction, which alters the factor to .111We note some typographical differences between Table 1 and Tables I and II in DHS. Row 10 here corresponds to Row 6 of Table I in DHS. While the expression is the same as that of DHS, the particle indexed 3 of our Row 10 is an electron, and particle 3 of Row in Table I of DHS is a positron, which should result in a different expression. This discrepancy also occurs between our Row 11 and Row 6 of Table (2) in DHS. Our expression for , however, agrees with that of Row 7 of Table I in Ref. Hannestad and Madsen (1995). The Weinberg angle is taken as . At the energy scales of interest here the and neutrino species have the same interactions.

ii.1.1 Collision integrals

Given the amplitudes of Table 1, we may calculate the collision integral of Eq. (1):


where is the symmetrization factor for identical particles, and


Here we have suppressed time dependence and written the occupation probability functions in abbreviated form. For example, for would read . The quantities , corresponding to the first and second lines of Eq. (II.1.1), give the probability for scattering into or out of the phase space volume for particle “1”; they include Pauli blocking factors . The phase space measure for particles , and , and the arguments of the four-momentum conserving delta function and of are written schematically with the dependence of on , which can either be four-momentum or , suppressed. The factor ensures that an integral over of the collision integral for vanishes in number-conserving processes; this is discussed in more detail in Sec.II.3. All amplitudes in Table 1 are proportional to , the Fermi coupling constant. The square of the Fermi coupling and a factor of may be taken outside of the collision integral [Eq. (8)] to give a dimensionless expression with integration variable , the binning parameter for the occupation probabilities. The product has dimensions of energy or inverse time, appropriate to that for a rate. The expression for the collision integral appearing in Eq. (1) is


for processes that include .

In general, for processes, Eq. (8) is a nine-dimensional integral over the phase space of particles 2, 3, and 4. The four-momentum conserving function reduces the collision integral to five dimensions. Homogeneity and isotropy further reduce Eq. (8) to a two-dimensional expression in terms of single-particle energies of either species 2 and 3, or 2 and 4, or 3 and 4. The method of the reduction to two-dimensions is distinct from and independent of that of DHS. The reduction is tuned for the specific process in Table 1 to ensure speed and accuracy in a parallel computation. Appendix B details our straightforward but lengthy method to obtain the two-dimensional expression for the process in the first row of Table 1. Appendix C gives the reduction algorithm for collision integrals for the remaining processes of Table 1. In both appendices, we relabel the indices of the active particle species in Table 1 to simplify the presentation.

ii.2 Numerical evaluation

In the interest of providing a complete description of the numerical evaluation of the collision integrals of Eq. (8) we describe here our choices for the energy () binning, numerical quadrature, interpolation and extrapolation, and convergence criteria.

ii.2.1 Binning

We employ a linear binning scheme for the occupation probabilities in terms of the comoving invariant quantity . The interval from to is partitioned into equal-width bins. For a linear binning scheme, we use abscissas with the lowest abscissa at . The must be chosen large enough to support the and . We compare the numerically integrated equilibrium energy spectrum to the analytical FD calculation at high temperature and find agreement to a few parts in .

We have performed test calculations with values of from 100 to 1000. The computing time has been verified empirically to scale as . Computation of the nuclear reaction network and thermodynamic quantities associated with charged leptons and photons incurs minimal computational overhead. Parallel code implementation of the calculations results in reasonable wall-clock times days on processors even with fine binning. Typically, we find convergence for =100, as discussed later in this section.

ii.2.2 Charged lepton quantities

For the processes relevant to weak decoupling, the occupation probabilities for the charged leptons are required. We assume these are given by the FD equilibrium spectra with chemical potential and temperature :


where is the electron degeneracy parameter. Here is determined by the requirement of charge neutrality in the electron/positron/baryon plasma. We assume zero lepton number residing in the neutrino seas222Our approach allows non-zero lepton asymmetry. The assumption of zero neutrino lepton number is stipulated for the present work and is in accord with the standard model. and neglect neutrino-nucleon charged-current transfer of electron lepton number between the electrons/positrons and electron neutrino/antineutrinos. This is plausible since the baryon-to-photon ratio is small. Finite electron mass is taken into account in Eq. (12) where . We define, for future use, the scaled mass as


We employ the comoving temperature, as its evolution is simple.

ii.2.3 Numerical quadrature

Upon reduction of the three-body, nine-dimensional momentum integrals as detailed in Appendix B, we may effect the remaining momentum integrations, which are transformed to integrals over , via numerical quadrature. We refer to the integration performed first (second) as “inner” (“outer”). We neglect the neutrino rest mass and divide the energy (or, equivalently, momentum) variable by to obtain the variables , where refers to either inner or outer integrations. If the integral is over a charged-lepton kinematic variable, we use its energy. The squared-amplitude expressions require both energies and three-momenta. We determine dimensionless momenta as , where is given in Eq. (13).

Depending on the specific process in Table 1, the integral may be over a neutrino or a charged lepton. For the inner integral, irrespective of the species, the integration method is a Gaussian quadrature method Elhay and Kautsky (1987). When the limits of the inner integral are finite, we use Gauss-Legendre. For finite intervals over a range of larger than 200 and semi-infinite intervals, we use Gauss-Laguerre.

When the outer integral is over an -value of a charged lepton, we use either a Gauss-Legendre or Gauss-Laguerre method, depending on the integration limits. In the case that the outer integral is over a neutrino energy, we use a five-point (Boole’s) rule Press et al. (1993) with abscissas aligned with the bin points. This affords a slight improvement in performance by avoiding interpolation for this integration of the occupation probabilities for the neutrino energy of the outer integral.

ii.2.4 Interpolation and Extrapolation

As detailed in Appendix C, we have the freedom to choose which single-particle -values to use in calculating the collision integral. The processes in Table 1 have at least two neutrinos in the combined initial and final states. We use three of the four energy- or momentum-conserving functions to eliminate an integral over the phase space of one of the neutrino species. This procedure requires an interpolation over the -value of that species to determine the occupation probability. Processes that involve four neutrinos or antineutrinos require an additional interpolation over for the occupation probability of the inner integration variable species. The outer integration is either a Gaussian quadrature method over a charged lepton, or a Boole’s rule method over the bin points. In either case, no interpolation is required. There is no situation in which we need to interpolate the occupation probabilities for the charged leptons since they are given by the equilibrium FD expressions for electrons and positrons.

We use a fifth-order polynomial interpolator Press et al. (1993) for the neutrino occupation probabilities if the energy of the third or fourth neutrino does not fall on an abscissa. The accuracy of this interpolation is better if we interpolate on the logarithms of the occupation probabilities, as opposed to the occupation probabilities themselves. The domain of integration is extended beyond , to , by extrapolation. Beyond this point the occupation probability is taken to be zero. None of our results are sensitive to these extrapolations.

ii.2.5 Acceptance tolerance for rates

When the occupation probabilities in Eq. (II.1.1) are all equilibrium-distribution values, the collision integral is zero, independent of the value of the squared matrix elements. Numerical quadrature and interpolation, however, incur errors at the precision limitations of these methods and the collision integrals attain small values when calculated under equilibrium conditions. During the computation the need arises to set the tolerance to accept a collision integral value as non-zero or, conversely, to reject a value as the result of imprecision. To accomplish this task, we use the net rate and forward-reverse-summed (FRS) rate. The net rate is the value given by the collision integral in Eq. (8). The FRS rate corresponds to the sum of contributions to the collision integral by substituting for [Eq. (II.1.1)] into Eq. (8).

We calculate the net and FRS rates for each neutrino and antineutrino species in each bin for all processes (and the antineutrino versions of interactions ) in Table 1 assuming thermal and chemical equilibrium between the three flavors of neutrinos, antineutrinos, positrons and electrons. We sum over all of the processes to obtain the collision integral for the net rate, and a modified collision integral for the FRS rate. For each neutrino species and each bin, we calculate the precision ratio, defined as:


where is the equilibrium FD occupation probability for species at a given bin. The FRS rate is numerically strictly positive. The absolute value of the net rate is required to obtain a strictly positive precision ratio since negative values can arise in and near equilibrium due to finite numerical precision. For diagnostic purposes only (i.e., not in our transport calculations) we take , meaning no temperature dependence in Eq. (14).

In production runs of burst, during weak decoupling, we calculate the collision integrals for both the net and FRS rates at each time step. We compare the ratio of values of the net and FRS rates for the evolved, in general non-equilibrium distributions , to those of the equilibrium distributions [Eq. (14)]. If the ratio in (14) is larger than the tolerance threshold


the collision integral is accepted as non-zero and used in the evaluation of the time derivative of the occupation probability . If the left-hand side of Eq. (15) is smaller than the threshold, we set the collision integral to zero. The precision ratio never gets larger than a few parts in .

ii.3 Conservation sum rules

We have tested the convergence of the numerical quadrature of the collision integrals by studying number and energy sum rules. Accurate evaluation of the collision integral is necessary to maintain the conservation of energy-momentum, particle number (for species with conserved charges), and neutrino lepton number. These are discussed in the following two sections.

ii.3.1 Number and energy sum rules

We define the total scaled errors in the number and energy densities as:


respectively. The summation over is for the three flavors of neutrinos and antineutrinos and the denominators in these expressions are strictly positive. We evaluate the sum rules including contributions only from processes isolated within the neutrino seas, i.e., in Table 1 to gauge the effectiveness of the numerical evaluation in respecting number and energy conservation. The spectra of the charged leptons are assumed to be described by equilibrium distributions so scattering processes involving electrons and positrons will not preserve the sum rules as written in Eqs. (16) and (17).

The neutrinos are assumed, in our computational approach, to be in thermal equilibrium with the electrons and positrons until a temperature MeV. The comoving temperature and plasma temperature are equal for all temperatures greater than the input temperature: . At , we commence evaluation of the collision integrals and evolve the neutrino occupation probabilities until a comoving temperature . The computation approach adopted in burst utilizes an adaptive Cash-Karp Press et al. (1993) time step. It evolves observables at steps on the interval defined by and with a fifth-order Runge-Kutta (RK5) algorithm. All simulations in this paper have , , , and . The terminal temperature is , corresponding to a plasma temperature of . In thermal equilibrium, the total scaled errors are small but non-zero and evaluate to for both the number and energy sum rules for 100 bins.

All 0.9282 0.3771
10, 11 0.9383 0.2867
1, 2, 10, 11 0.9268 0.2963
1, 2, 3, 4, 5, 10, 11 0.8557 0.3465
6, 7, 8, 9 0.1853 0.0639
1, 2, 6, 7, 8, 9 0.1724 0.0778
1, 2, 3, 4, 5, 6, 7, 8, 9 0.1559 0.0886
Table 2: Process-dependent changes in neutrino energy density properties. For all runs , , , , . The first column gives the processes used for a given run. The second column is the ratio of comoving to plasma temperature. For column two reference, . Columns three and four are the relative changes of the and energy densities. The quantity is given by Eq. (24). Round-off error of the neglected fifth significant digit in columns 2, 3, and 4 accounts for the one part in discrepancy with column 5.

We monitor the total scaled errors of Eqs. (16) and (17) at each time step during our weak decoupling calculations. On average, we maintain accuracy to better than one part in over the entire run.

ii.3.2 Neutrino lepton number conservation

Elastic processes satisfy


since the processes and (and their antineutrino counterparts) conserve neutrino (antineutrino) number. The annihilation processes, and 11 satisfy, for example:


Analogous relations hold for other annihilation processes that fall under the reaction classes and 11.

We have confirmed that the neutrino lepton numbers are conserved at the level of for all values of the scale factor .

Iii Results in the neutrino sector

Our treatment of the Boltzmann-equation evolution of the neutrino energy transport reveals novel features of the transport characteristics of the active neutrino sector. We focus first on these results, which are largely independent of the coupling to BBN through the nuclear reaction network. The present calculations reveal, in particular, that the history of annihilation to photons displays a rich set of behaviors that has not been discussed before. We also look into the role of QED radiative corrections. These results are in line with previous work but they indicate that a more comprehensive treatment of the plasma physics during the epochs we consider is warranted.

iii.1 Neutrino interactions and energy transport

Table 2 summarizes the neutrino energy transport properties in the present calculations, which as mentioned are carried out for computational parameters , and . In this section, we focus on the first row of the table, when all of the weak interactions of neutrinos (and the antineutrino reactions corresponding to the parity conjugates of the reactions ) are computed. We discuss the results for selective process evaluations corresponding to the remaining rows of this table in the next section, Sec. III.2. We briefly describe this table to orient the subsequent discussion.

The first column of Table 2 lists the processes from Table 1 used for a given run. The second column gives the ratio of the comoving to plasma temperatures. Columns three and four give the relative changes of the and energy densities, respectively, with respect to the equilibrium energy density:


The last column is the change in . is defined through the energy density in ultra-relativistic particles – the radiation energy density – after the epoch of photon decoupling as


where is the scale factor at universal comoving time , is the plasma or photon temperature and is the photon temperature at the conclusion of the epoch of photon decoupling. We make the assumption that and do not change significantly for . Therefore, if we set the radiation energy density equal to the sum of the photon and neutrino densities in Eq. (22), we can determine from and :


In writing Eq. (III.1), we have assumed that antineutrinos have the same relative change in energy density as neutrinos. The change in is given as:


where is given by Eq. (III.1). It is clear from this table that the dominant contribution to the parameter is due to annihilation processes & 11. Additionally, for , the value typically quoted in the literature Mangano et al. (2005), the effect of charged lepton scattering is small but not negligible.

Another feature apparent in Table 2 is non-linearity in the combination of processes. Adding, for example, the values of for Table 2 rows 4 and 5, which sums all of the processes (and the implied charged-lepton antineutrino scattering processes) with is not equivalent to the Table 2 first row with . Finally, the ratio of the comoving temperature to the plasma temperature is largely set by the annihilation processes. We note, however, that this does not uniquely determine as Eq. (22) implies.

Figures 4, 2, 3 and 1 show relative changes in the neutrino spectra for a calculation with transport versus a no-transport calculation. All processes are active in the transport calculation, i.e. row 1 of Table 2. The no transport calculation maintains FD-like distributions at temperature parameter . We compare our present results, in detail, to DHS and Ref. Esposito et al. (2000). To this end, we first define several quantities to facilitate this comparison and then turn to a detailed discussion of each of these figures.

We define at a given time and to be the relative change in the occupation probabilities with respect to the FD occupation probability:




We note that does not depend explicitly on time or temperature. Figures 1 and 2 show as a function, respectively, of for and 7 and as a function of at a comoving temperature . Figure 3 displays the difference in the relative change for neutrinos and antineutrinos:


Figure 4 shows the normalized change in the differential energy density:


The antineutrino behavior is nearly identical to the neutrino behavior for all flavors.

For each -value in Fig. 1 the relative change in the electron-flavor () is larger than the relative change in the muon-flavor () neutrino sea. The annihilation and scattering rates with electrons and positrons are faster due to the contribution of the charged-current diagrams for , which are absent for , as noted in Ref. Esposito et al. (2000). In addition to the larger affect on the spectra, the charged-current processes keep the in thermal contact with the charged leptons longer than . This is apparent from Fig. 1 where freeze-out corresponds to the point where the derivative of the curves goes to zero. The freeze-out occurs at an earlier epoch than the freeze-out. Additionally, freeze-out occurs later for larger values, as noted in DHS. These results are generally consistent with DHS and Ref. Esposito et al. (2000). For example, the , curve rises at a more rapid rate than, and crosses, the , curve. Figure (4) of Ref. Esposito et al. (2000) also exhibits this crossing between the , curve and the , curve. Comparing Fig. 1 with Figs.(3a) and (3b) of DHS, confirms the similar behavior of burst and DHS.

Figure 1: (Color online) The relative change, as in Eq. (25), in the occupation probability as a function of the comoving temperature . Three values of are evaluated at and . The solid lines are for electron-flavor neutrinos, and the dashed lines are for muon-flavor neutrinos. The larger correspond to larger values.
Figure 2: (Color online) The relative change, as in Eq. (25), in the occupation probability as a function of for keV. The larger change is the electron-flavor neutrinos, over the muon-flavor neutrinos. The antineutrino evolution is nearly identical to the neutrino evolution for all flavors.

Figure 2, plotted at a temperature of well after weak decoupling, shows that the have a larger distortion than the and that this effect is enhanced at large . An interesting feature of Fig. 2 is the negative relative change for . It appears to occur in Fig. (5) of both DHS and Ref. Esposito et al. (2000) but is not explicitly mentioned in either reference. We investigate this phenomenon in more detail in the subsections below. In addition, our relative changes are in good agreement with those of DHS for both and .

In Fig. 3, we exhibit the difference in the relative change of and , and also and . The electron-flavor shows an enhanced effect over the muon-flavor for all -values. For both flavors, at -values and in Fig. 3, the relative changes are positive. The negative differences for indicate there is an abundance of antineutrinos over neutrinos, independent of flavor. The differences between neutrino and antineutrino distributions for both and are small for all epsilon values considered here. This raises, however, the important issue of how neutrino flavor evolves under the full quantum kinetic evolution Vlasenko et al. (2014).

Figure 3: (Color online) The difference in relative changes in the occupation probabilities of and [Eq. (27)] as a function of comoving temperature . Three values of are plotted at and . The solid lines are for electron-flavor neutrinos, and the dashed lines are for muon-flavor neutrinos. The experience a larger change than the .

Figure 4 shows where the largest change in the energy-density spectrum occurs. Figure 4 is approximately equivalent to Fig. 2 multiplied by . The peak of the normalized change in the differential energy density is located at , for both and . Figure (6) of Ref. Esposito et al. (2000) also shows a peak at an . Although Fig. 2 shows that the deviation from equilibrium of the occupation probabilities increases for increasing -values, the probability is small enough in the high- bins that the large changes from equilibrium have little effect on the total energy density.

Figure 4: (Color online) The normalized change in the differential energy density [Eq. (29)] as a function of . The electron neutrinos exhibit a larger change compared to the muon neutrinos. The antineutrino evolution is nearly identical to the neutrino evolution for all flavors.

Integrating the neutrino energy distributions in Fig. 4, we find relative changes in the energy density of: , and . The temperature ratio is given in the first row of Table 2 as . Using Eq. (24), we find . The quantities , , , and all agree closely with both Ref. Esposito et al. (2000) and DHS.

Figure 5 shows how the energy densities, , and evolve with until they reach their asymptotic values. The and are computed from Eq. (21) and the relative change in is computed by comparing the evolution of the temperature with transport and without :


Finally, to calculate the time evolution change in , we use


where the subscript denotes time dependence, in contrast to the asymptotic limit of Eq. (24). As may be seen in Fig. 5, does not converge to 3.034, the value consistent with DHS. The reason its asymptotic value is instead 3.033 is due to the fact that the run with no transport has in the asymptotic limit. If we assume the neutrinos are in thermal equilibrium for , the temperature ratio incurs a modification from the finite electron rest mass as:


to second order in . Setting , we find an altered gives .

Figure 5: (Color online) Quantities related to energy density and temperature are plotted against the comoving temperature parameter. The blue solid curve shows the change in using Eq. (31). The red dashed curve shows the relative change in the energy density of . The green dash-dot curve shows the relative change in the energy density of . The magenta dotted curve shows the relative change in using Eq. (30). At a given , is equivalent to a lower plasma temperature in the transport case compared to no transport.

The evolution of in Fig. 5 displays interesting features that are driven by the specifics of the loss of entropy in the plasma from the annihilation of electrons and positrons to neutrinos and the transfer of entropy from electrons/positrons to photons through annihilation (see Sec. IV for a detailed discussion of entropy). The annihilation of electrons and positrons into neutrinos can be seen in the rise of the curves in Fig. 5. For , entropy is lost from the plasma into the neutrino seas resulting in a lower plasma temperature for the transport case (where entropy is lost) versus the no-transport case (where entropy is not lost). The increase in for is caused by this entropy loss.

To analyze the entropy transfer from the electron/positron components to the photons, we need the total number densities of electrons and positrons


which, in local thermodynamic equilibrium, are solely functions of the plasma temperature and the electron chemical potential. Fig. 5 shows a different phasing of scale factor and temperature for the two cases. At a given , the plasma temperature is always lower in the transport case versus the no-transport case.

Using a notation similar to that of Eq. (30), we define the absolute change in the number density of charged leptons as


and the relative change in the number as


The quantity is proportional to the total number of electrons and positrons in a comoving volume. The dot-dashed curve in Fig. 6 shows this quantity, while the solid curve in Fig. 6 shows the relative change in the number of charged leptons. The absolute change in the total number of charged leptons in a comoving volume is negative. This implies that there are fewer charged leptons and hence a lower plasma temperature in the transport case than in the no-transport case at a given . The slope of the absolute change represents the different annihilation rates. The negative slope (at ) indicates that the annihilation rate is greater in the transport case than in no-transport, while the opposite is true for .

The rate at which entropy is transferred from the electrons and positrons to the photons is proportional to the annihilation rate divided by the plasma temperature. For , the competition between a larger annihilation rate in the no-transport case and the lower plasma temperature in the transport case results in a slight decrease in . Finally for , the greater annihilation rate in no-transport results in an increasing until virtually no electrons and positrons remain and reaches its asymptotic value.

We turn now to the study of the individual and joint contributions of the annihilation and elastic processes.

Figure 6: (Color online) Quantities related to charged lepton number density are plotted against the comoving temperature parameter. The blue solid curve is the relative change in the sum of positron and electron number densities as calculated in Eq. (35). The red dash-dot curve is the absolute change in the sum of positron and electron number densities, divided by .

iii.2 Neutrino energy transport analysis

We return to the discussion of the weak interaction processes in Table 1 and their effect on neutrino observables that probe the changes to their energy density – the comoving-plasma temperature ratio, energy density changes, and . The component contributions are collected in the following subsections in terms of the annihilation and elastic channels.

iii.2.1 Annihilation Channel

We focus here on annihilation-channel effects, which for the present purposes are defined according to Table 1 as processes and 11. Figure 7 shows the relative changes in the occupation probabilities for the annihilation channels as a function of . The solid curves are when all weak interaction transport processes are neglected, except for annihilation of a neutrino and antineutrino into an electron-positron pair. The dashed curves are the same as the solid curves with the addition of processes and 2. The dotted curves further include processes and 5. These processes do not exchange population among individual -values for a given flavor. Instead, there is an equilibration between the and flavors. This is clearly seen in the figure since the difference between the dashed and dotted curves decreases relative to the solid, annihilation curves.

Figure 7: (Color online) The change in the neutrino occupation probabilities relative to the equilibrium distributions as a function of . Electron neutrinos exhibit a larger change compared to the muon neutrinos. Solid lines correspond to processes and 11. Dashed lines correspond to processes and 11. Dotted lines correspond to processes and .

iii.2.2 Elastic Scattering Channel

Figure 8 shows changes relative to equilibrium for the combinations of processes involving elastic scattering of neutrinos on the charged leptons. It is distinguished by its behavior at low , where the change relative to equilibrium goes negative for , a feature not found in Fig. 7. Neutrinos, whose number are conserved in processes , upscatter and populate the larger epsilon bins. increases, and conversely, the energy in the plasma decreases. Processes behave much the same way as they do for the annihilation combinations. These processes act to equilibrate the occupation probabilities among flavors for a given bin. Perhaps surprisingly, the five neutrino-only processes do not appear to equilibrate the bins for a given flavor; the addition of the processes in rows do not change the intersection with the horizontal axis at . The further additions of the processes in rows also preserve the intersection point with the horizontal axis. Note that in both Figs.7 and 8 there is a larger divergence between the solid, dashed, and dotted lines with increasing -value for the as compared to the . This is because the population is transferred from into both and .

Figure 8: (Color online) The change in the neutrino occupation probabilities relative to the equilibrium distributions as a function of . The electron neutrinos again exhibit larger changes relative to the muon neutrinos. Solid lines correspond to processes . Dashed lines correspond to processes . Dotted lines correspond to processes .

When adding the annihilation and elastic scattering channels together, as in Fig. 2, the annihilation channels are able to repopulate the low -values states. Annihilation erases much of the deficit caused by elastic scattering, although Fig. 2 shows that annihilation cannot entirely erase the deficit for .

iii.3 Finite temperature QED radiative corrections

Our calculations show a new sensitivity to finite temperature QED radiative corrections. The feedback between neutrino energy transport and plasma conditions is especially sensitive to electron-positron pair number density as these are key targets for neutrino scattering. Previous works Dicus et al. (1982); Cambier et al. (1982); Heckler (1994); Lopez and Turner (1999) include the effects of finite temperature QED radiative corrections. The corrections have been calculated for the electron mass and wave function renormalization, the electron-photon vertex and infrared photon emission and absorption. These corrections have an effect on a variety of quantities including the dispersion relation of the electron mass, weak interaction rates, and the equation of state.

In the present study, we include radiative corrections to the self-interaction energy for electrons, positrons, and photon dispersion relations. We follow the approach employed in Ref. Mangano et al. (2002). This formulation is adopted primarily for comparison with previous results. We note, however, that the feedback between neutrino energy distribution evolution and plasma conditions can depend sensitively on the electron-positron pair density and that, in turn, can depend on these corrections. This highlights the need for a more complete treatment of finite temperature plasma effects. In any case, finite temperature QED radiative corrections by themselves have a small effect on, say, the relative abundance change, which is at the level of , smaller than the neutrino transport effects that we are primarily concerned with in this work.

We apply finite temperature QED radiative corrections to contributions to the thermodynamic quantities , , , , , and . For the collision integrals of Sec. II.1, we take the vacuum value of electron rest mass. In a preview of Sec.V, the weak interaction terms involving neutrinos and free nucleons utilize the electron rest mass at its vacuum value and do not include the higher-order effects detailed in Refs. Dicus et al. (1982); Cambier et al. (1982); Lopez and Turner (1999). Our rates do not take into account the renormalization of the electron rest mass in any BBN computations. We have computed the effects of non-zero electron degeneracy and found them to be negligible, so we take in these radiative corrections. We maintain in calculations of neutrino transport, weak rates, and overall charge neutrality with baryons.

Following Refs. Heckler (1994); Mangano et al. (2002) we take the shift in the electron rest mass to be