Lepton asymmetry, neutrino spectral distortions, and big bang nucleosynthesis

Lepton asymmetry, neutrino spectral distortions, and big bang nucleosynthesis

E. Grohs    George M. Fuller    C. T. Kishimoto    Mark W. Paris Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA Department of Physics, University of California, San Diego, La Jolla, California 92093, 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
July 5, 2019

We calculate Boltzmann neutrino energy transport with self-consistently coupled nuclear reactions through the weak-decoupling-nucleosynthesis epoch in an early universe with significant lepton numbers. We find that the presence of lepton asymmetry enhances processes which give rise to nonthermal neutrino spectral distortions. Our results reveal how asymmetries in energy and entropy density uniquely evolve for different transport processes and neutrino flavors. The enhanced distortions in the neutrino spectra alter the expected big bang nucleosynthesis light element abundance yields relative to those in the standard Fermi-Dirac neutrino distribution cases. These yields, sensitive to the shapes of the neutrino energy spectra, are also sensitive to the phasing of the growth of distortions and entropy flow with time/scale factor. We analyze these issues and speculate on new sensitivity limits of deuterium and helium to lepton number.

cosmological parameters from CMBR, big bang nucleosynthesis, cosmological neutrinos, neutrino lepton numbers, neutrino theory
preprint: LA-UR-16-29176

I Introduction

In this paper we use the burst neutrino-transport code Grohs et al. (2016) to calculate the baseline effects of out-of-equilibrium neutrino scattering on nucleosynthesis in an early universe with a nonzero lepton number, i.e. an asymmetry in the numbers of neutrinos and antineutrinos. Our baseline includes: a strong, electromagnetic, and weak nuclear reaction network; modifications to the equation of state for the primeval plasma; and a Boltzmann neutrino energy transport network. We do not include neutrino flavor oscillations in this work. Our intent is to provide a coupled Boltzmann transport and nuclear reaction calculation to which future oscillation calculations can be compared. In fact, the outstanding issues in achieving ultimate precision in big bang nucleosynthesis (BBN) simulations will revolve around oscillations and plasma physics effects. These issues exist in both the zero and nonzero lepton-number cases, but are more acute in the presence of an asymmetry.

We self-consistently follow the evolution of the neutrino phase-space occupation numbers through the weak-decoupling-nucleosynthesis epoch. There are many studies of the effects of lepton numbers on light element, BBN abundance yields. Early work Wagoner et al. (1967); Schramm and Wagoner (1977) briefly explored the changes in the helium-4 () abundance in the presence of large neutrino degeneracies. Later work considered how lepton numbers could influence the yield Shi (1996); Kirilova and Chizhov (1998) through neutrino oscillations. In addition, other works employed lepton numbers to constrain the cosmic microwave background (CMB) radiation energy density Hansen et al. (2001); Simha and Steigman (2008) or the sum of the light neutrino masses Shiraishi et al. (2009). Refs. Kneller et al. (2001); Mangano et al. (2011) simultaneously investigated BBN abundances and CMB quantities using lepton numbers. The most recent work has used the primordial abundances to constrain lepton numbers which have been invoked to produce sterile neutrinos through matter-enhanced Mikheyev-Smirnow-Wolfenstein (MSW) resonances Abazajian et al. (2005); Smith et al. (2006); Chu and Cirelli (2006). Currently, our best constraints on these lepton numbers come from comparing the observationally-inferred primordial abundances of either or deuterium (D) with the predicted yields of and D calculated in these models.

Previous BBN calculations with neutrino asymmetry have made the assumption that the neutrino energy distribution functions have thermal, Fermi-Dirac (FD) shaped forms. In fact, we know that neutrino scattering with electrons, positrons and other neutrinos and electron-positron annihilation produce nonthermal distortions in these energy distributions, with concomitant effects on BBN abundance yields Grohs et al. (2016). Though the nucleosynthesis changes induced with self-consistent transport are small, they nevertheless may be important in the context of high precision cosmology. Anticipated Stage-IV CMB measurements Carlstrom and et al. (2016); Bond et al. (2017) of primordial helium and the relativistic energy density fraction at photon decoupling, coupled with the expected high precision deuterium measurements made with future 30-meter class telescopes Kirkman et al. (2003); Pettini and Cooke (2012); Cooke et al. (2014); Cooke and Pettini (2016); Cooke et al. (2016) will provide new probes of the relic neutrino history.

In the standard cosmology with zero lepton numbers, neutrino oscillations act to interchange the populations of electron neutrinos and antineutrinos (, ) with those of muon and tau species (, , , ) Kostelecký and Samuel (1995). Once we posit that there are asymmetries in the numbers of neutrinos and antineutrinos in one or more neutrino flavors, then neutrino oscillations will largely determine the time and temperature evolution of the neutrino energy and flavor spectra Savage et al. (1991); McKellar and Thomson (1994); Casas et al. (1999); Wong (2002); Abazajian et al. (2002); Dolgov et al. (2002); Gava and Volpe (2010); Johns et al. (2016); Barenboim et al. (2016). In this paper we ignore neutrino oscillations and provide a baseline study of the relationship between neutrino spectral distortions arising from the lengthy ( Hubble times) neutrino decoupling process and primordial nucleosynthesis. This is an extension of the comprehensive study of this physics in the zero lepton-number case with the burst code Grohs et al. (2016), and in other works Dolgov and Fukugita (1992); Dolgov et al. (1997, 1999); Esposito et al. (2000); Serpico and Raffelt (2005); Mangano et al. (2005); Smith et al. (2008); Saviano et al. (2013); de Salas and Pastor (2016). We will introduce alternative descriptions of the neutrino asymmetry to study the individual processes occurring during weak decoupling. Our studies in this paper, together with the methods in other works, will be important in precision calculations for gauging the effects of flavor oscillations in the early universe.

As we develop below, a key conclusion of a comparison of neutrino-transport effects with and without neutrino asymmetries is nonlinear enhancements of spectral distortion effects on BBN in the former case. This suggests that phenomena like collective oscillations may have interesting BBN effects in full quantum kinetic treatments of neutrino flavor evolution through the weak decoupling epoch.

The outline of this paper is as follows. Section II gives the background analytical treatment of neutrino asymmetry, focusing on the equations germane to the early universe. Sec. III presents the rationale in picking the neutrino-occupation-number binning scheme and other computational parameters. We use the same binning scheme throughout this paper as we investigate how the occupation numbers diverge from FD equilibrium, starting in Sec. IV. In Sec. V, we present a new way of characterizing degenerate neutrinos in the early universe. Sec. VI details the changes to the primordial abundances from the out-of-equilibrium spectra. We give our conclusions in Sec. VII. Throughout this paper we use natural units, , and assume neutrinos are massless at the temperature scales of interest.

Ii Analytical Treatment

To characterize the lepton asymmetry residing in the neutrino seas in the early universe, we use the following expression in terms of neutrino, , antineutrino, , and photon, , number densities to define the lepton number for a given neutrino flavor


where . The photons are assumed to be in a Planck distribution at plasma temperature , with number density


where . The neutrino spectra have general nonthermal distributions and their number densities are given by the integration


Here, is the comoving temperature parameter and scales inversely with scale factor


where the subscripts reflect a choice of an initial epoch to begin the scaling. In this paper, we will choose such that is coincident with the plasma temperature when . For , the plasma temperature and comoving temperature parameter are nearly equal as the neutrinos are in thermal equilibrium with the photon/electron/positron plasma. and diverge from one another once electrons and positrons begin annihilating into photon and neutrino/antineutrino pairs below a temperature scale of . The dummy variable in Eq. (3) is the comoving energy and related to , the neutrino energy, by . The sets of are the phase-space occupation numbers (also referred to as occupation probabilities) for species indexed by . In equilibrium the occupation numbers for neutrinos behave as FD


where is the neutrino degeneracy parameter related to the chemical potential as . Unlike the lepton number for flavor in Eq. (1), the corresponding degeneracy parameter is a comoving invariant. If we consider the equilibrium occupation numbers in the expression for number density, Eq. (3), we find


where is the relativistic Fermi integral given by the general expression


We can define the following normalized number distribution


Figure 1 shows plotted against for three different values of .

Figure 1: Normalized number density plotted against for three choices of degeneracy parameter: nondegenerate (, solid blue), degenerate with an excess (, dashed red), and degenerate with a deficit (, dash-dotted green).

The expressions for the number densities in Eqs. (2) and (3) have different temperature/energy scales. As the temperature decreases, electrons and positrons will annihilate to produce photons primarily, thereby changing with respect to . As a result, Eq. (1) decreases from the addition of extra photons. To alleviate this complication, we calculate the lepton number at a high enough temperature such that the neutrinos are in thermal and chemical equilibrium with the plasma. We can take at high enough temperature and write Eq. (1) as


where we call the comoving lepton number. Eq. (9) simplifies further if we use the FD expression in Eq. (5) and recognize that in chemical equilibrium the degeneracy parameters for neutrinos are equal in magnitude and opposite in sign to those of antineutrinos


where is the degeneracy parameter for neutrinos of flavor . Eq. (10) provides an algebraic expression for relating the lepton number to the degeneracy parameter with no explicit dependence on temperature. We will give our results in terms of comoving lepton number and use Eq. (10) to calculate the degeneracy parameter for input into the computations. In this paper, we will only consider scenarios where all three neutrino flavors have identical comoving lepton numbers. Unless otherwise stated, we will drop the subscript and replace it with the neutrino symbol, i.e. , to refer to all three flavors.

Degeneracy in the neutrino sector increases the total energy density in radiation. The parameter is defined in terms of the plasma temperature and the radiation energy density


Eq. (11) can be used at any epoch, even one in which there exists seas of electrons and positrons, e.g. Eq. (31) in Ref. Grohs et al. (2016). We will consider and at the epoch , after the relic seas of positrons and electrons annihilate. Assuming equilibrium spectra for all neutrino species, the deviation of , , from exactly 3 would be


where the summation assumes the possibility of different neutrino degeneracy parameters for each flavor Esposito et al. (2000); Shimon et al. (2010).

We begin by presenting the case of instantaneous neutrino decoupling with pure equilibrium FD distributions. Table 1 shows the deviations in energy densities for neutrinos and antineutrinos with respect to nondegenerate FD equilibrium, the asymptotic ratio of to , and the change to , for various comoving lepton numbers. In this paper, we will colloquially refer to the asymptote of any quantity as the “freeze-out” value. For the values of in Table 1, a decade decrease in produces comparable decreases in and . is related to the energy densities through the degeneracy parameter derived from Eq. (10), which is approximately linear in for small . The change in is quadratic in which is discernible for and at the level of precision presented in Table 1. The freeze-out value of is not identically , the canonical value deduced from covariant entropy conservation Kolb and Turner (1990); Weinberg (2008). Although the neutrino-transport processes are inactive for Table 1 and therefore the covariant entropy is conserved, finite-temperature quantum electrodynamic (QED) effects act to perturb away from the canonical value Heckler (1994); Fornengo et al. (1997).

Table 1: Observables and related quantities of interest for zero and nonzero comoving lepton numbers without neutrino transport. Column 1 is the comoving lepton number. Columns 2 and 3 give the relative changes of the and and energy densities compared to a FD energy distribution with zero degeneracy parameter at freeze-out. Column 4 shows the ratio of comoving temperature parameter to plasma temperature also at freeze-out. For comparison, . Column 5 gives . does not converge to precisely in the nondegenerate case due to the presence of finite-temperature-QED corrections to the equations of state for photons and electrons/positrons.

Iii Numerical Approach

For this work, changes to the quantities of interest will be as small as a few parts in . To ensure our results are not obfuscated by lack of numerical precision, we need an error floor smaller than the numerical significance of a given result. In burst, we bin the neutrino spectra in linear intervals in -space. The binning scheme has two constraints: the maximum value of to set the range; and the number of bins over that range. We denote the two quantities as and , respectively, and examine how they influence the errors in our procedure.

The mathematical expressions for the neutrino spectra have no finite upper limit in . We need to ensure is large enough to encompass the probability in the tails of the curves in Fig. 1. As an example, consider the normalized number density in Eq. (8). We would numerically evaluate the normalization condition as


For large , , and so we exclude a contribution to the above integral on the order of if we take . If we are using double precision arithmetic, the contribution becomes numerically insignificant for , which corresponds to . This value of would seem like the natural value to take without loss of a numerically significant contribution to the integral in Eq. (13). However, if we fix the number of abscissa in the partition used when integrating Eq. (13) (i.e. fixing in the binned neutrino spectra), we lose precision in the evaluation of the contribution to the integral from each bin as we increase . Clearly, there is a trade off between and .

Figure 2 examines the versus parameter space by looking at the calculation of the equilibrium comoving lepton number, in a scenario where . We take to be exactly and solve the cubic equation in Eq. (10) for . Next, we calculate neutrino and antineutrino spectra with the equal and opposite degeneracy parameters. We proceed to integrate Eq. (9) with the two spectra for different pairs of values. The integration is carried out using Boole’s rule, a fifth-order integration method for linearly spaced abscissas. Figure 2 shows the filled contours of values for the error in


for a given pair . We immediately see the loss of precision in the upper-left corner of the parameter space, corresponding to small and large . Furthermore, for , the error value flat lines with increasing , implying that the error is a result of a too small choice for . The black curve superimposed on the heat map gives the value of with the lowest error as a function of . It monotonically increases for , at which point it reaches and begins to fluctuate. The fluctuations are a result of reaching the double-precision floor, implying that increasing adds no more numerical significance.

Figure 2: versus for contours of constant error in . The exact value of is . The black line on the contour space gives the value of with the smallest error as a function of .

The computation time required to run burst scales as . In this paper, we attempt to be as comprehensive as possible when exploring neutrino energy transport with nonzero lepton numbers. Therefore, we will choose bins for the sake of expediency. Fig. 2 guides us in picking , and dictates a floor of for our best possible precision. It would appear that the choice would give the absolute best precision for calculating . This is valid if using the linearly spaced abscissas as a binning scheme. We highlight both the precision and timing needs for a comprehensive numerical study on binning schemes. Such a study would be germane for the more general problem which includes neutrino oscillations and disparate lepton numbers in the three active species de Salas and Pastor (2016); Barenboim et al. (2016); Johns et al. (2016).

For more details on the numerics of burst, we refer the reader to Ref. Grohs et al. (2016). We have essentially preserved the computational parameters except for a quantity related to the determination of nonzero scattering rates. Ref. Grohs et al. (2016) used , and in this work, we use .

Iv Neutrino Spectra

In this section we give a detailed accounting of how the neutrino energy spectra evolve through weak decoupling in the presence of zero and nonzero lepton numbers. In the first subsection we integrate the complete transport network, including all the neutrino scattering processes in Table I of Ref. Grohs et al. (2016), from a comoving temperature parameter down to . In the second subsection we investigate how the different interactions between neutrinos and charged leptons affect the spectra.

We compare our results to that of FD equilibrium. For the neutrino occupation numbers, we use the following notation to characterize the deviations from FD equilibrium


Here, is the FD equilibrium occupation number for degeneracy parameter given in Eq. (5). When it is obvious, we will drop the argument , i.e. . As an example, gives the relative difference of the occupation number from the nondegenerate, zero chemical potential FD equilibrium value.

We also examine the absolute changes for the number and energy distributions


When using the absolute change expressions, we normalize with respect to an equilibrium number or energy density in order to compare to dimensionless expressions. For the energy density, we use the appropriate degeneracy factor


For the number density, we will exclusively use zero for the degeneracy factor


The out-of-equilibrium evolution of the neutrino occupation numbers driven by scattering and annihilation processes with charged leptons does not proceed in a unitary fashion. Consequently, the total comoving neutrino number density increases. The increase in number results in an increase in energy density, and so we use to normalize the absolute changes in differential energy density distribution to compare with the initial distribution at high temperature. However, the difference in number density between neutrinos and antineutrinos, characterized by the comoving lepton number in Eq. (9), does not change with kinematic neutrino transport. In practice, burst follows the evolution of neutrino and antineutrino occupation numbers separately, precipitating the possibility of numerical error. We will use the same normalization for neutrino and antineutrino differential number density distributions to study the relative error in . We will take the normalization quantity to be that of the nondegenerate number density in Eq. (19).

iv.1 All processes

Table 2 shows how neutrino transport alters neutrino energy densities, , the ratio of comoving temperature parameter to plasma temperature, and entropy per baryon in the plasma, . These quantities are computed for a range of values and refer to the results at the end of the transport calculation, , well after weak decoupling. In this table, we focus on the energy-derived quantities. The relative changes in energy density are with respect to a nondegenerate FD distribution at the same comoving temperature, i.e.


Columns 2 - 5 of Table 2 show the relative changes in energy density at , once the neutrino spectra have converged to their out-of-equilibrium shapes. We see a monotonic decrease in for the neutrinos, and a monotonic increase in for the antineutrinos with decreasing . Column 6 gives the ratio of at the end of the simulation. increases with decreasing lepton number. However, the decrease is less than one part in between and . The larger lepton number implies a larger total energy density which increases the Hubble expansion rate. The faster expansion implies a smaller time window for the entropy flow out of the plasma and into the neutrino seas. As a result, the evolution of the plasma temperature is such that larger lepton numbers will maintain at higher values, and the ratio at freeze-out will decrease, albeit by an amount which is numerically insignificant. With the changes in energy densities and temperature ratios, we can calculate


The coefficient in front of the second parenthetical expression, equal in value to , results from the approximation in taking the and flavors to behave identically. The approximation employed here is valid as there are no and charged leptons in the plasma and is the same in all flavors. Both Refs. Mangano et al. (2005); de Salas and Pastor (2016) calculate weak decoupling with a network featuring neutrino flavor oscillations, which are absent in our calculation in Table 2. However, Ref. de Salas and Pastor (2016) states that oscillations have no affect on the value of at the level of precision which they use. The difference in our value of versus the standard calculation of Ref. Mangano et al. (2005) is most likely due to a different implementation of the finite-temperature QED effects detailed in Refs. Heckler (1994); Fornengo et al. (1997). Ref. Mangano et al. (2005) uses the perturbative approach outlined in Ref. Mangano et al. (2002) compared to our nonperturbative approach. We leave a detailed study of the finite-temperature-QED-effect numerics to future work.

The final column of Table 2 shows the change in the entropy per baryon in the plasma. The relative changes in entropy for varying lepton numbers are large enough to see a difference at the level of precision Table 2 uses, unlike . With the faster expansion, neutrinos have less time to interact with the plasma, yielding a smaller entropy flow.

Table 2: Observables and related quantities of interest in zero and nonzero lepton-number scenarios with neutrino transport. Column 1 is the lepton number. Columns 2, 3, 4, and 5 give the relative changes of the , , and energy densities compared to a FD energy distribution with zero degeneracy parameter. Comparisons are given for , after the conclusion of weak decoupling. Column 6 shows the ratio of comoving temperature parameter to plasma temperature also at the conclusion of weak decoupling. For comparison, . Column 7 gives as calculated by Eq. (21). Column 8 gives the fractional change in the entropy per baryon in the plasma, .

An increase in lepton number implies a larger energy density for the neutrinos over the antineutrinos. Fig. 3 shows four neutrino spectra after the conclusion of weak decoupling in a scenario where . Plotted against is the relative difference in the neutrino occupation number with respect to a nondegenerate spectrum. As seen in the first data row of Table 2, obtains the largest difference from equilibrium. The thick red line in Fig. 3 shows the final out-of-equilibrium spectrum for . The spectrum has the largest deviation from equilibrium, congruent with Table 2. The black dashed lines show equilibrium spectra for nondegenerate (flat, horizontal line) and degenerate cases. As increases, the neutrino curves diverge from the positive spectrum in much the same manner as the antineutrino curves diverge from the negative . The primary difference in the out-of-equilibrium spectra is due to the initial condition that the neutrinos have larger occupation numbers over the antineutrinos for a positive lepton number.

Figure 3: Relative differences in neutrino/antineutrino occupation numbers plotted against at . The relative differences are with respect to FD with zero degeneracy parameter. The solid lines show the evolution for a scenario where . and curves are colored red, and and curves are colored green. Neutrinos have thick line widths and antineutrinos have thin line widths. Plotted for comparison are black dashed curves representing the equilibrium relative differences. The top dashed curve corresponds to a spectrum with , the middle horizontal curve corresponds to a spectrum with , and the bottom dashed curve corresponds to a spectrum with .

We would like to compare the out-of-equilibrium spectra to their respective equilibrium spectra. Such a comparison allows us to examine how the initial asymmetry propagates through the Boltzmann network. Fig. 4 shows the evolution of for (thick solid lines) and (thin solid lines) for a scenario where . We only show the relative differences for three unique values of , namely . The and spectra follow similar shapes, but are suppressed relative to the electron flavors. For comparison, we also plot the out-of-equilibrium spectrum for in the case of no initial asymmetry, i.e.  identically zero. It is unnecessary to show the spectrum for when because it is exceedingly near the spectrum (see Fig. [3] of Ref. Grohs et al. (2016)). For the degenerate spectra, the show a larger divergence from equilibrium than the at these three specific values. This is consistent with Ref. Esposito et al. (2000) (see Figs. 8 and 9 therein) and is the case for all after the neutrino spectra have frozen out. Fig. 5 shows the final freeze-out values of the relative changes in the neutrino occupation numbers as a function of . Fig. 5 is similar to Fig. 3 except for the use of the general instead of . We have also included the transport-induced out-of-equilibrium spectra for and in the nondegenerate scenario. For a given flavor, the relative changes in the nondegenerate spectrum are nearly averages of those in the and spectra. We also note that for , all of the relative differences are negative, although this is obscured in Fig. 5 due to the clustering of lines. For small , the antineutrino occupation numbers are larger than those of the neutrino, i.e. the are not as negative.

Figure 4: Relative differences in electron neutrino/antineutrino occupation numbers plotted against . The relative differences are with respect to FD with the same degeneracy parameters as Fig. 3. The solid lines show the evolution for a scenario where . The (thin red curves) has a larger relative change than the (thick red curves). Plotted for comparison is the relative difference for in a scenario (blue dash-dot curves). The relative differences are plotted for three values of , from bottom to top: .
Figure 5: Relative differences in neutrino/antineutrino occupation numbers plotted against at . The relative differences are with respect to FD with the same corresponding degeneracy parameters as Fig. 3. The solid lines are for a scenario where . Plotted for comparison is the relative difference for (blue dash-dot curve) and for (blue dotted curve) in a scenario.

The weak interaction cross sections scale as , where is the Fermi constant () and is the total lepton energy. We would expect a larger difference from equilibrium for increasing . Except for the range , Figures 4 and 5 clearly show an increase. The change in the energy distribution does not follow from a scaling relation. Fig. 6 shows the normalized absolute difference in the energy distribution plotted against at the conclusion of weak decoupling. The nomenclature for the six lines in Fig. 6 is identical to that of Fig. 5. The energy distributions all show a maximum at . Similar to Fig. 5, the nondegenerate curves of Fig. 6 appear to be averages of the and curves in the degenerate scenario.

Figure 6: Absolute change in the neutrino/antineutrino energy distributions plotted against at . The changes are with respect to the same degeneracy parameters as those in Fig. 5. Furthermore, the line colors and styles correspond to the same species and scenarios as Fig. 5.

In the positive lepton-number scenarios, the always have larger occupation numbers than the , when compared against the equilibrium degenerate spectrum/distribution. This is not surprising as the occupation numbers for antineutrinos are suppressed, implying less blocking. When compared against its equilibrium distribution, the have larger rates, leading to a larger distortion. In Fig. 7, we compare the out-of-equilibrium number density distributions with those of the nondegenerate case solely. In other words, the normalizing factor is the same for each of the six curves in Fig. 7. We have adopted this nomenclature for the comparison of number density distributions to study the change in the comoving lepton number. None of the weak decoupling processes modify the lepton number in our model. The total change in number density for should be identical to the total change in number density for . Fig. 7 shows this indirectly. We can see a difference; the curves are skewed to higher and have a larger maximum than the . The negative change in the distributions for the range is much more noticeable in Fig. 7 than in Fig. 5. It is clear that the changes in become positive for smaller than those of , implying there are more than for . Overall, when integrating the curves in Fig. 7, the total changes in number density for should be the same as for . We have calculated this quantity and expressed it as a relative change in the , taken to be exactly


Eq. (22) gives the relative error in our calculation. We conserve the comoving lepton number for both electron and muon flavor at approximately . Also plotted in Fig. 7 are the absolute changes for and in the nondegenerate scenario. We do not directly compare the lepton-number relative errors as the quantity is not defined for the symmetric case. We do note that the nondegenerate curves are close to the average of the and distributions, similar to that of Figs. 5 and 6.

Figure 7: Absolute change in the neutrino/antineutrino number distributions plotted against at . The changes are with respect to the same degeneracy parameters, and the nomenclature of line colors and line styles is the same as those in Fig. 5. The numbers given on the plot ( for ) show the relative error accumulated over the course of a simulation.

In Figs. 3 through 7, we have only presented the scenario. Fig. 8 shows the relative differences in occupation number for plotted against for other values of . The behavior of each curve is in line with those of Fig. 5. Not plotted are the curves for . They also behave in a similar manner, where becomes larger than for increasing . The result is that with transport, acts to increase the asymmetries in the occupation numbers, which manifest in differences in the absolute changes of the differential energy density.

Figure 8: Relative differences in neutrino occupation numbers plotted against at for various values of . The relative differences are with respect to FD with the corresponding degeneracy parameter, namely (), (), (), and (). Shown are two sets of curves: the set with the larger relative differences correspond to the spectral distortions, and the set with the smaller differences are . Within each set, the increase in leads to a decrease in . For the antineutrinos, the relative differences behave in the opposite manner: increase in leads to an increase in .

iv.2 Individual processes

Figs. 5 and 6 demonstrate that the initial asymmetry in the neutrino energy density is maintained and even amplified by scattering processes. We can dissect the relative contribution of various scattering processes to this amplification.

Figs. 9 and 10 show the absolute changes in the number density distribution versus when we include only certain transport processes. Fig. 9 contains three annihilation processes, schematically shown as:


In this scenario, we have included only the annihilation channel into electron/positron pairs when computing transport. The changes are with respect to the same degeneracy parameters as those in Fig. 5. The line colors in Fig. 9 correspond to the same species as Fig. 5. Because of the close proximity of the neutrino and antineutrino curves, we depart from the previous nomenclature of emphasizing the curves with a thicker line width so as not to obscure the curves. For this plot, the absolute differences are normalized with respect to the equilibrium number density at temperature with degeneracy parameter . For a given neutrino species, the total change in number density should be equal to the change in number density for the corresponding antineutrino.

Figure 9: Absolute change in the neutrino/antineutrino number distributions plotted against at . The numbers given on the plot ( for ) show the relative error accumulated over the course of a simulation.

Fig. 10 shows the effect of including 12 elastic scattering processes:


and the opposite- reactions, for neutrino flavors . In this scenario, we have included only the elastic scattering channel with electrons/positrons (while neglecting the neutrino-antineutrino only channels) when computing transport. The changes are with respect to the same degeneracy parameters as those in Fig. 5. Furthermore, the line colors and styles in Fig. 10 correspond to the same species and scenarios as Fig. 5. For this plot, the absolute differences are normalized with respect to the equilibrium number density at temperature with degeneracy parameter . In an identical manner to the processes in Fig. 9, the total change in number density should be equal to the change in number density for the corresponding in Fig. 10.

Figure 10: Absolute change in the neutrino/antineutrino number distributions plotted against at . The numbers given on the plot ( for ) show the relative error accumulated over the course of a simulation.

The elastic scattering processes of Eqs. (24) and (25) (and the opposite- reactions) preserve the total number of neutrinos and antineutrinos. The plasma of charged leptons acts to upscatter low energy neutrinos and antineutrinos to higher energies, precipitating an entropy flow. Fig. 10 vividly shows a deficit of neutrinos in the range , and the corresponding excess for . The deficit is more pronounced in Fig. 10 but also appeared in Figs. 5, 6, and 7 when computing the entire neutrino-transport network. The annihilation processes, shown in Fig. 9, do not preserve the total numbers of neutrinos and antineutrinos and can fill the phase space vacated by the upscattered neutrinos. The complete transport network, which includes annihilation, elastic scattering on charged leptons, and elastic scattering among only neutrinos/antineutrinos, is able to redistribute the added energy by filling the occupation numbers for lower epsilon.

V Integrated asymmetry measures

In our presentation to this point, we have used the comoving lepton number to describe the asymmetry in the early universe. does not evolve with temperature in our model, except for errors in precision encountered by our code. Therefore, we introduce two integrated quantities to examine how the initial asymmetry propagates to later times. The quantities provide new means to analyze the out-of-equilibrium spectra.

The first integrated quantity we define is the lepton energy density asymmetry


where is the flavor index. Like the comoving lepton number in Eq. (9), we divide Eq. (26) by so that is comoving and dimensionless. This will allow us to follow the evolution of to later times. At large , all flavors have identical equilibrium FD spectra and lepton numbers/degeneracy parameters. For degeneracy parameter , we calculate the equilibrium value of


where is the sign function with real-number argument , and is the Lerch function (see Sec. 9.55 of Ref. Gradshteyn and Ryzhik (2007))

Figure 11: The relative changes in plotted against for and with different processes included in the transport calculation. Red lines correspond to and green lines correspond to . The process scheme is all processes (solid curves), annihilation only (dash-dot curves), or elastic scattering only (dotted curves).

Figure 11 shows the relative changes in from the baseline (), plotted against for different combinations of transport processes. Solid lines (All) are for the complete calculation, whereas dash-dot curves only include the annihilation channels (Annih.) of the reaction shown in (23), and dotted curves only include the elastic scattering channels (Scatt.) of the reactions shown in (24), (25), and the opposite- reactions. Red lines correspond to and green lines to . increases for all six combinations of flavor and transport process, until an eventual freeze-out. Indirectly, Figs. 7, 9, and 10 all show that the neutrinos have larger changes in the energy density distributions, increasing the asymmetry. Because of the charged-current process, experiences a greater enhancement. What is important to note is that the total , for either flavor, is not an incoherent sum of the two transport processes taken individually. There are two reasons for this.

First, there are other transport processes in the full calculation. Neutrinos scattering on other neutrinos and antineutrinos will redistribute energy density. Second, the transport processes with the charged leptons are dependent on one another. Positron-electron annihilation into neutrino-antineutrino pairs populates the lower energy levels. Those particles upscatter on charged leptons through elastic scattering. Positron-electron annihilation is then suppressed by the Pauli blocking of the upscattered particles. Both reasons change the evolution of the total , but do so in a flavor-dependent manner. For , the incoherent sum of annihilation and elastic scattering is smaller than that of the total asymmetry. For , the total asymmetry is dominated by the contribution from elastic scattering.

In analogy with the lepton energy density asymmetry, we define the lepton entropy asymmetry as


where is the entropic density for particle , given by


and we have suppressed the arguments of for brevity in notation. Under the equilibrium assumptions, we find


Fig. 12 shows the evolution of the relative change in away from when divided into processes. The nomenclature for the line styles and colors is identical to that in Fig. 11. The evolution of the lepton entropy asymmetry shows more features than that of the lepton energy density asymmetry.

Figure 12: The relative changes in plotted against . Line colors and styles correspond to the same transport processes and neutrino flavors in Fig. 11.

To understand the dynamics of in Fig. 12, we begin by considering how the entropy depends on perturbations to the occupation numbers. We write the occupation numbers as differences from FD equilibrium


We can calculate the change in the entropy produced by the out-of-equilibrium occupation numbers by substituting Eq. (32) into Eq. (30). After dropping the argument, argument, and species index for notational brevity, we find for small


where and are the changes in number and energy density, respectively, from equilibrium. The expression for the lepton entropy asymmetry is


Lepton number is conserved in our scenarios, implying . As a result, we can write the lepton entropy asymmetry as


Eq. (38) shows how the lepton entropy asymmetry changes for small perturbations to the occupation numbers. Two trends are evident from this equation. First, adding particles () decreases the asymmetry. Second, increasing the asymmetry in energy density (), leads to an increase in the lepton entropy asymmetry. For the annihilation processes, the changes in the number density distribution for neutrinos and antineutrinos vary in the same way across space for all flavors (see Fig. 9). Therefore, the corresponding changes in the energy density will also be the same, and there will be no contribution to the change in from the energy density terms. The dash-dot curves in Fig. 12 shows the relative change in for a run with only the annihilation channels active. Both the and flavors show a suppression in with decreasing . Figure 10 shows that for elastic scattering of neutrinos and charged leptons, the neutrino and antineutrino number density distributions are not coincident. Overall, each neutrino species has zero net change in number density, as elastic scattering can only redistribute the number. Therefore, there will be no contribution to the change in from the number density term. As there are more neutrinos over antineutrinos for , elastic scattering enhances the neutrino spectra over the antineutrino spectra. The result is a net positive change in the energy density differences. Fig. 12 shows an increase in the relative change in for the elastic-scattering-only runs for both flavors. When we add the elastic-scattering and annihilation channels together, along with the other transport processes which do not involve charged leptons, we see that the two processes essentially cancel, leaving only a modest change in as shown by the solid lines in Fig. 12.

Figure 13: Same as Fig. 12 except zoomed in on the solid curves.

The interesting thing to note in Fig. 12 is the asymmetry between flavors. Fig. 13 is a zoomed-in version of the solid lines in Fig. 12. We see that is monotonically increasing for decreasing . The incoherent sum of the relative changes from the annihilation and elastic-scattering processes in Fig. 12 nearly gives the relative change in that we obtain when all transport processes are active. The same cannot be said for . The sum of the two transport processes is not incoherent, the evolution of is not monotonic, and the final freeze-out value of is of opposite sign from . Although the elastic scattering would appear to produce a larger enhancement of over the suppression of annihilation, the two processes do not have equal weight. We observe this by looking at the maxima in the number density distributions in Figs. 9 and 10. The ratio of maxima in Fig. 9 for annihilation is


The ratio of maxima in Fig. 10 for elastic scattering is


This shows that annihilation is more dominant in the electron neutrino/antineutrino sector than it is in the sector. In Figs. 9 and 10, we have only showed the final distributions at freeze-out. Electron-positron annihilation into neutrinos is not always so dominant, as evidenced by the positive values of for .

The analysis of the lepton entropy asymmetry focused on the transport processes which involve the charged leptons. The other scattering processes redistribute occupation number and therefore change . However, we have verified that the contributions from the transport processes which involve only neutrinos or antineutrinos do not alter enough to explain the full evolution shown in Fig. 13. The transport processes which involve the charged leptons play the dominant roles.

We have considered the evolution of the integrated asymmetry measures for only. Table 3 gives the relative changes in and at freeze-out for various values of . Note that the positive relative changes for imply an absolute decrease in either quantity. We see that the differences between the various values of are beneath the error floor.

Table 3: Relative changes in integrated asymmetry measures at freeze-out for select values of comoving lepton number .

Vi Abundances

Our calculations show potentially significant changes in lepton-asymmetric BBN abundance yields with neutrino transport relative to those without. With the inclusion of transport we find that the general trends of the yields of and D with increasing or decreasing lepton number are preserved: positive decreasing the yields of both, while negative lepton numbers increase both. In broad brush, Boltzmann transport makes little difference for helium, but gives a reduction in the offset from the FD, zero lepton-number case with transport. This change in the reduction is comparable to uncertainties in BBN calculations arising from nuclear cross sections and from plasma physics and QED issues. For all BBN calculations, the baryon to photon ratio is fixed to be (equivalent to the baryon density given by Ref. Planck Collaboration et al. (2014)). In addition, the mean neutron lifetime is taken to be .

Table 4 contains relative differences in the primordial abundances with and without transport. Columns with the label “FD Eq.” are the calculations without any active transport processes. The spectra freeze-out at high temperatures where they are in FD equilibrium with a degeneracy parameter corresponding to . Columns with the label “Boltz.” are the calculations in the full Boltzmann neutrino-transport calculation. Relative differences are with respect to the appropriate abundance in the zero-degeneracy Boltz. calculation. The relative changes in the abundances for the two different calculations are quite close: differs by 2 - 3 parts in ; and differs by 3 - 4 parts in . Both differences are consistent across . We caution against any interpretation that links the two calculations together, as the No Trans. calculations ignore important physics related to non-FD spectra, entropy flow, and the Hubble expansion rate.

We have examined the detailed evolution of the spectra and integrated asymmetry measures in the Boltz. calculations. The electron neutrinos and antineutrinos behave differently compared to muon and tau flavored neutrinos. This behavior will have ramifications for the neutron-to-proton ratio and nucleosynthesis. To facilitate the analysis of the effects of neutrino transport on BBN, we will introduce a model which uses additional radiation energy density. We will try to determine whether this simplistic “dark radiation” model Grohs et al. (2015); Mukohyama (2000) – which includes radiation energy density distinct from photons and active neutrinos, but does not include transport – can mock up the effects of the extra energy density which arise from neutrino scattering and the associated spectral distortions. We will compare this dark-radiation model to the full neutrino-transport case. For ease in notation when comparing the two scenarios, we will abbreviate the dark-radiation model as “DR” and the full Boltzmann neutrino-transport calculation as Boltz.

(FD Eq.) (Boltz.) (FD Eq.) (Boltz.)
Table 4: Relative changes in primordial abundances of and D in two calculations of BBN with nonzero comoving lepton numbers . FD Eq. signifies the calculation without transport. Boltz. signifies the full Boltzmann neutrino-transport network calculation. The abundances are given as relative changes from the zero degeneracy, full Boltzmann calculation. Column 1 is the comoving lepton number. Column 2 gives the relative change of at freeze-out in the no-transport model. Column 3 gives the relative change of at freeze-out in the Boltz. calculation. The relative changes for D/H are given in columns 4 and 5. The four rows where is not a power of 10 are projected sensitivity limits for changes in the primordial abundances.

In the DR model, we introduce extra radiation energy density, , described at early times by the dark-radiation parameter


The FD Eq. calculation in Table 4 used . We mandate that the dark radiation be composed of relativistic particles which are not active neutrinos. We have chosen the specific form of Eq. (41) for conformity with , namely . The relation is not a strict equality due to the presence of finite-temperature-QED corrections to the electron rest mass Grohs et al. (2016); Heckler (1994); Fornengo et al. (1997); Cambier et al. (1982); Lopez and Turner (1999). The DR model differs from the Boltz. calculation in multiple respects. First, the DR model fixes the neutrino spectra to be in degenerate FD equilibrium. Second, neutrino transport induces an entropy flow from the plasma into the neutrino seas, absent in the DR model. Third, the entropy flow changes the phasing of the plasma temperature with the comoving temperature parameter as compared to the case of instantaneous neutrino decoupling in the DR model. The phasing is dependent on the Hubble expansion rate and the flow of entropy. Although the expansion rates are identical in the two scenarios, the entropy flows are not.

For all calculations, we will fix