Probing neutrino physics with a self-consistent treatment of the weak decoupling, nucleosynthesis, and photon decoupling epochs

We show that a self-consistent and coupled treatment of the weak decoupling, big bang nucleosynthesis, and photon decoupling epochs can be used to provide new insights and constraints on neutrino sector physics from high-precision measurements of light element abundances and cosmic microwave background observables. Implications of beyond-standard-model physics in cosmology, especially within the neutrino sector, are assessed by comparing predictions against five observables: the baryon energy density, helium abundance, deuterium abundance, effective number of neutrinos, and sum of the light neutrino mass eigenstates. We give examples for constraints on dark radiation, neutrino rest mass, lepton numbers, and scenarios for light and heavy sterile neutrinos.


Introduction
In this paper we demonstrate how a self-consistent treatment of the physics of the early universe (from temperature T ∼ 10 MeV down to T ∼ 0.1 eV) enables cosmic microwave background (CMB) observations, in concert with light-element abundances, to be used as new probes of beyond-standard-model (BSM) physics in the neutrino sector. In a sense, these probes are tantamount to probes of the CνB (relic Cosmic Neutrino Background). Recent observations [1][2][3][4] of the CMB and observationally-inferred primordial abundances of deuterium and helium [5][6][7][8] formed in big bang nucleosynthesis (BBN) already place tight constraints on both the cosmological standard model (CSM) and BSM physics. However, future observations will usher in a higher level of precision with even better prospects for BSM probes, see for example Ref. [9].
We anticipate that observations will bring about an overdetermined situation where BSM physics may manifest itself if it is present. The existing or anticipated observations and measurements of greatest utility for the present purposes are: (1) high-precision measurements of the baryon-to-photon ratio, or the equivalent baryon density (ω b ≡ Ω b h 2 ); (2) high-precision measurements of the "effective number of relativistic degrees of freedom" (N eff ); (3) high-precision measurements of the primordial deuterium abundance (D/H) from quasar absorption lines; (4) measurements of the primordial helium abundance (Y P ) directly from CMB polarization data; and (5) measurements of the sum of the light neutrino masses ( m ν ), i.e. the collisionless damping scale associated with neutrinos.
The physics that determines the relic neutrino energy spectra in weak decoupling and that of primordial nucleosynthesis affects observables of the CMB. These distinct, disparate epochs, however, depend not only on the values of parameters describing the cosmology (such as ω b , N eff , relevant cosmic constituents, etc.) but also on the parameters derived from these base, input parameters (such as primordial abundances, recombination history, etc.). This is the basis for what we will term "self-consistency." The requirement for self-consistency between the BBN and CMB epochs, for example, in current data analyses depends on the parametrization of the energy density in terms of N eff and upon the baryon-to-photon ratio. These two parameters, as measured at photon decoupling, are the sole determinants of the primordial helium abundance at the epoch of alpha particle formation (T ∼ 0.1 MeV) in "standard" (i.e. zero lepton numbers, no BSM physics) BBN calculations. Though this procedure is adequate for the CSM and standard model physics, here we argue that it can be insufficient to probe varieties of BSM physics when confronting model cosmologies with next-generation, high-precision CMB and light element abundance data. In the case of the helium and N eff self-consistency mentioned above, the relationship between the helium yield in BBN and N eff is a non-trivial function of the interplay of expansion rate and neutron-toproton ratio, as is well known [10]. The latter ratio is, in turn, a sensitive function of the ν e andν e energy distribution functions, and these can be affected by BSM issues like lepton numbers, flavor mixing, sterile neutrino states, heavy particle decay, etc.
Handling and analyzing the observed data while imposing self-consistency over multiple epochs in the early universe can require new procedures. The approach developed in Ref. [10] and employed in Planck XVI [11] has been highly successful. There, self-consistency is imposed approximately by the addition of a term in the log-likelihood function. A very small theoretical uncertainty ensures the posterior distributions are close to the corresponding theoretical values. However, when we impose self-consistency among the different epochs the need may arise to solve for the various derived parameters iteratively. Returning to the helium/N eff example discussed above, Y P and the recombination history of the universe both depend on the radiation energy density, usually parametrized by N eff . Subsequently, the recombination history is affected by the the primordial abundance of helium. The values of these two derived parameters, in turn, affect CMB observables, which recommends an iterative and therefore self-consistent approach.
Ideally, we would want a procedure to self-consistently treat neutrino and BSM physics from the post-QCD epoch to the onset of non-linearities in large scale structure (LSS). Recently, the quantum-kinetic equations (QKEs) governing neutrino flavor evolution in dense environments have been derived from first principles [12][13][14]. Such a program to treat neutrino physics in the early universe would incorporate a QKE treatment through weak decoupling at the very least. Weak freeze-out and BBN necessitate a neutrino-energy Boltzmannequation or QKE treatment fully coupled to a nuclear reaction network. The recombination, photon decoupling, and advent of LSS epochs require a Boltzmann-equation treatment of neutrino clustering. Verification of any models related to neutrinos and BSM physics would then rely on agreement with direct cosmological observables, including but not limited to: primordial abundances; CMB power spectra; and the total matter power spectrum.
As outlined, this is a challenging undertaking. We therefore employ a limited approach to show why self-consistency is required and how such a treatment can be efficacious in probing some issues in neutrino sector physics. To that end, we calculate the ratio of sound horizon to photon diffusion length, r s /r d , as a simple parametrization of the CMB, as described in detail in Sec. 3. It would, of course, be preferable to compute the full CMB power spectra in the self-consistent manner described above. There are two limiting factors that temper this ambitious proposal, however. The first is that the observable r s /r d is largely insensitive to the poorly constrained equation of state of the dark, vacuum energy component while the CMB power spectra are not. Additionally, current computations of the CMB power spectra [15] would need to be generalized to the BSM scenarios we contemplate here.
Our procedure for calculating derived cosmological quantities utilizes the neutrino occupation probabilities. Since we are developing the capacity to compute the effects of BSM physics that couples to the active neutrino sector, we do not restrict the form of the neutrino distribution functions to that of equilibrium Fermi-Dirac distributions. General forms are permitted and may be handled analytically or numerically. Neutrino occupation probabilities are taken in the present work in the flavor eigenbasis during weak decoupling, weak freeze-out and BBN. 1 During recombination and last scattering, we transform the occupation probabilities to the mass eigenbasis. When we consider the case of massive neutrinos, the statistic for the sum of the light neutrino masses, m ν , is simply the sum of the three lightest mass eigenstates.
We might further clarify the fact that this limited approach does not rely on the assumption that the radiation energy density can be described in terms of the single parameter N eff . Indeed, an original motivation for using r s /r d as a proxy for N eff (see Sec. 3) is the fact that the usual definition for N eff does not apply to non-equilibrium neutrino distributions. A corollary of this observation is the possibility that N eff may depend on the scale factor a(t); that is, we do not assume that N eff is independent of scale factor. When extended, our approach is able to handle BSM scenarios such as the decay of massive sterile neutrinos and other weakly interacting massive particles, which result in neutrino distributions that may differ substantially from a Fermi-Dirac distribution [16].
Our approach aims to be general; it does not rely on particular cosmological models but is appropriate to the class of Λ cold-dark matter (ΛCDM) models. It extends and explicates recent work [9,10,17,18] on the consistent incorporation of precision observations of the CMB [11] and observations of the light element abundances of helium and deuterium (mass fraction Y P and relative abundance D/H). The deuterium abundance [6,8] is now more precisely measured by a significant factor (4 or 5) in relative precision than Y P . Our results indicate that deuterium is sufficiently well determined to be incorporated into CSM analyses as a fixed prior.
The present work is a first result in an ongoing campaign to incorporate neutrino energy transport into a self-consistent treatment of BBN and recombination. Our approach is being implemented in Fortran90/95 as a suite of codes under the working title of "BBN Unitary Recombination Self-consistent Transport" (burst), which will be made publicly available for use on parallel computing platforms using openmpi. The current work is a prerequisite in an ongoing collaboration to develop codes that consistently handle BBN and neutrino energy transport from weak decoupling to the advent of LSS.
The paper is outlined as follows. In the following section, Sec. 2, we discuss our treatment of BBN and relevant CMB physics. Section 3 explains the self-consistent, iterative approach to determining N eff using r s /r d . Section 4 employs a dark-radiation model to test the accuracy of predictions made by burst against observations. Section 5 investigates a number of 'test problems' from BSM physics suited for burst. We conclude in Sec. 6 with future applications of the results given in this work. Throughout this paper, we use natural units where = c = k B = 1.

Big-bang nucleosynthesis
Recently there has been a dramatic increase in precision in the determination of D/H [8]. This improvement in observational precision drives the need to improve the standard tools used to calculate observables during the BBN and CMB epochs. This development allows us to consider improved descriptions of nuclear physics and non-standard particle and cosmological models.
We have developed a generalized BBN subroutine, part of the larger burst code. The current BBN routine handles general distribution functions for particles out of equilibrium. It numerically integrates the binned neutrino occupation probabilities to obtain thermodynamic variables and number and energy densities for ν e ,ν e , ν µ ,ν µ , ν τ ,ν τ . It is worthwhile, therefore, to summarize our approach to BBN even though it is substantially similar in some respects to that given in Refs. [19][20][21].
We assume that the cosmic fluid is homogeneous and isotropic and that the comoving number density of baryons is covariantly conserved. Neutrinos can experience, in addition to gravity, essentially arbitrary interactions within the Boltzmann transport approximation. 2 We allow for the possibility that neutrinos may decouple from the plasma with non-equilibrium distributions. This assumption implies that there may be deviations from the Fermi-Dirac momentum distribution [22][23][24]. In fact, it was this observation that provided initial motivation for developing the current approach.
BBN codes evolve light nuclide abundances Y i (t), defined as as a function of comoving time t in the background of the Friedmann-Lemaître-Robertson-Walker (FLRW) geometry. Here n i (t) is the proper number density of nuclide i = n, p, 2 H, 3 He, 4 He, . . . and n b (t) is the proper baryon number density. We will focus on two specific nuclides in this paper: the 4 He mass fraction (Y P ≡ 4Y4 He ), and the 2 H relative abundance (D/H ≡ Y D /Y H ). The evolution starts from a time when the plasma temperature T is near 30 MeV. Weak equilibrium obtains at this temperature. The time dependence of the metric is determined by the energy density ρ(a) as where H(a) is the Hubble parameter, m pl is the Planck mass, and the 'dot' indicates differentiation with respect to the FLRW-coordinate time t. The energy density is computed as the integral of the single-particle energy over the momentum distribution: where the sum over j reflects all species contributing to the energy density, including but not limited to: photons, baryons, dark matter, e ± , and neutrinos. Consistent with the above discussion, we do not need to assume a well defined temperature for any of the cosmic species. The time dependence of the distribution function f i (p; a) is indicated by the presence of the scale factor a(t) in its argument. Evolution in time of the distribution functions is accomplished by solving transport equations, such as the Boltzmann equation, in FLRW geometry. Given initial conditions for the temperature T and momentum distributions of the cosmological constituents (assumed to be in equilibrium through weak and strong interactions), the BBN code determines the evolution, with respect to the scale factor a(t) (related to time as dt = da/(aH(a))), of the temperature of the plasma T (a), the electron chemical potential µ e (a), and the nuclide abundances relative to baryon number Y i (a).
The interactions of the light nuclear species is governed by the nuclear reaction network. 3 The reaction network, which is determined by a chosen set of nuclides and the thermally averaged reactivities n α 1 n α 2 σ βα v α , is proportional to the rate of change of the number density of nuclides participating either in the initial state n α i or final state n β j where i, j indexes particles (up to 3) in the initial α or final β channel for the process α → β; that is, α and β are two-or three-body reaction channels and α i is the i th nuclide of the channel α. The nuclear reaction network includes nuclides with mass number A ≤ 9. Our code allows for the inclusion of additional nuclides with A > 9, but we maintain the smaller network as the larger network provides no new insights for this paper. The nuclides are taken to be in thermal equilibrium with the photon-electron plasma. We are currently employing an updated [25] version of the reaction network from Ref. [21]. We also couple in all relevant e ± and neutrino-induced weak interactions (charged and neutral current) [26].

Cosmic Microwave Background
The physics of BBN and photon decoupling are related through the epoch of recombination by several mechanisms including, importantly, He i and He ii recombination. The recombination rates depend sensitively on the amount of 4 He created in BBN [10,27,28]. The 4 He dependence also determines the photon diffusion damping wave number k d or, equivalently, the photon diffusion length, defined as r d ≡ π/k d . This quantity depends strongly on the recombination history through the free-electron fraction, X e (a), which, in turn, depends on the helium abundance (see Sec.

below).
We relate the diffusion length directly to an observed angular size θ d through the angular diameter distance D A at the epoch of last scattering (as described in detail below). Since D A depends on the vacuum (dark) energy equation of state, which is poorly understood, we employ the ratio θ s /θ d , equivalent to r s /r d (see Eqs. (2.6) and (2.13)), where θ s is the angle subtended by the sound horizon, which eliminates explicit dependence on the dark energy equation of state.
Modern computations of the recombination history of the early universe are available in accurate and well tested codes such as HyRec [29], CosmoRec [30], and the fast, approximate computations given in Ref. [31] by the code RecFast [32]. Our recombination history compares well with RecFast, agreeing at the level of < 3% over the recombination epoch. The present accuracy suffices for the purposes of the current exploratory work.

Sound horizon
The sound horizon is defined as where a γd is the scale factor at photon decoupling. 4 The ratio R(a) is given as where ρ b and ρ γ are the baryon rest mass and photon field energy densities, respectively. The angular size of the sound horizon θ s 5 , determined from the spacing of the acoustic peaks in the CMB temperature power spectrum, is related to the sound horizon by the angular diameter distance D A (a γd ) at photon decoupling (scale factor a γd ) as since the angle θ s is small. The angular diameter distance is The quantity D A (a γd ) depends on the vacuum (dark) energy equation of state, which is not very well understood. Our approach will eliminate dependence on this poorly constrained component of the energy density.

Free-electron fraction
We have written an independent code to calculate the free-electron fraction, X e . The results agree well with Ref. [32] (RecFast). In fact, the agreement is within 2% for most of the range 10 −4 a/a 0 10 −3 (10 4 z 10 3 ). The code allows significant deviations from the model parameters of ΛCDM; it should be used with caution, however, since the effective three-level treatments for helium and hydrogen recombination have been optimized for near-standard model parameters.
The number density of free and total electrons is where n b , n p , n He ii , n He iii are the proper number densities for baryons, protons, singly-and doubly-ionized helium, respectively. We write the free-electron fraction as so defined to take values in the range 0 ≤ X e ≤ 1.
We follow Refs. [33,34] and consider the simplification of the multi-level hydrogen and helium atoms to that of an effective three-level system which includes the ground n = 1 state, the first excited n = 2 states, and the continuum. All other excited states are assumed to be in equilibrium with the 2s state. We treat He ii recombination approximately [31] (via the Saha equation) since it is essentially complete at the advent of the epoch of He i recombination. The He iii contribution, therefore, in the Boltzmann equation for He ii is negligible.
Boltzmann equations for H ii and He ii contain a thermally-averaged cross section and relative velocity σv . We use Case B recombination coefficients for the σv of both neutral hydrogen [35] and helium [36]. Along with a Saha equation for He iii, the Boltzmann equations for H ii and He ii are a coupled set of ordinary differential equations constituting a recombination network to model the ionization history of the universe prior to photon decoupling.

Photon diffusion damping length
In the tightly coupled limit, photon diffusion damping is characterized through the damping wave number [37][38][39][40] where σ T is the Thomson cross section. Here, we have assumed that moments of the temperature fluctuation higher than the quadrupole make a negligible contribution in the linearized Boltzmann equation for the photon distribution. Equation (2.11) requires the free-electron fraction, n e (a), determined in the recombination history, as discussed in the previous section. The free-electron fraction, over the course of the recombination history, depends strongly on the primordial helium mass fraction Y P . The fraction Y P , in turn, depends on the relativistic energy density, which is often parametrized in terms of the parameter N eff : Equation (2.12) is only valid after the epoch of e ± annihilation. The point that BBN and recombination are related is well known [10] but has not been implemented self-consistently as a constraint for general, BSM physics model cosmologies. We return to it after describing the relation between r s and k d to directly observable quantities given by the CMB power spectrum.
The observed diffusion angle θ d (a γd ) [11] is related to the diffusion damping length, with the same stipulation regarding the smallness of the angle in Eq.(2.6).
3 r s /r d as a proxy for N eff In this section we describe our method for determining N eff in detail. We introduce two variants of N eff . When referring to the radiation energy density equation (2.12), which takes as an input the quantity N eff , we designate N eff as N (th) eff . When considering general cosmologies, perhaps with BSM physics, we deduce the value of N eff from the observable quantity r s /r d = θ s /θ d , described in this section, and designate it asÑ eff . The simplest cosmologies for which Eq. (2.12) obtains, having negligible neutrino mass, standard model constituents and no energy transfer between species have N (th) eff =Ñ eff . We consider a test input cosmology that is non-standard yet substantively similar to ΛCDM. We proceed by determining Y P at temperature T ∼ 0.1 MeV using the BBN network of burst [25]. The principal observational cosmological input at this time is ω b . Also input and incorporated into the BBN subroutine is any neutrino and BSM physics that constitute the test cosmology. Subsequently, we compute the recombination history of the universe from early times (a/a 0 ∼ 10 −7 ) to the current epoch. Specific observational inputs include ω b , ω c (where for cold dark matter ω c ≡ Ω c h 2 ), Y P , and H 0 (the Hubble constant, H 0 = H(a 0 )). For the purposes of the present discussion, other inputs of particular importance include neutrino occupation probabilities (as output from a neutrino transport calculation of weak decoupling that is fully coupled to BBN) and neutrino rest masses. The neutrino energy density of the neutrino seas is calculated by writing the occupation probabilities in the mass eigenbasis [41]. We emphasize the fact that N (th) eff is not input as a base parameter; this is of paramount import in the present approach. The recombination history, n e (a) determines the optical depth as a function of scale factor: We define the scale factor at photon decoupling a γd such that τ (a γd ) ≡ 1. In this definition, we do not include the effects of cosmic reionization when calculating n e (a) for use in Eq.
(3.1) [11]. We apply a γd and the input cosmology to equations (2.4) and (2.11) to compute the sound horizon and photon diffusion length, respectively. We arrive in this way at a ratio (r s /r d ) (inp) for our input cosmology. Our immediate objective is to determine a value of N eff (here termedÑ eff ) corresponding to this value for (r s /r d ) (inp) . We map out a range of values of r s /r d that correspond to the same input cosmology as that used in calculating (r s /r d ) (inp) , with one significant difference. We parametrize all of the neutrino and BSM physics into the single N (th) eff parameter. We then use N (th) eff to calculate the radiation energy density in Eq. (2.12) to determine the Hubble rate. We vary N (th) eff to compute the function shown in Fig. 1 The final step is to evaluate the previous function with our input cosmology ratio, i.e.Ñ eff = N , to obtain a value ofÑ eff . As an example, we take the best-fit values from Ref. [11] combined with WMAP Polarization data (100θ s = 1.04136 & 100θ d = 0.161375) to obtain (r s /r d ) (inp) = 100θ s /100θ d = 6.45304. This corresponds to a valueÑ eff = 3.31 on Fig. 1, in line with the best-fit value N eff = 3.25 (3.51 +0.80 −0.74 at 95% limits) [11]. We again note that for the simplest cosmologies the two generally distinct functions N   We can understand Fig. 1 qualitatively by a simple scaling argument. We expect that, as the radiation energy density increases with increasing N (th) eff , the sound horizon and the diffusion length will decrease with the increasing Hubble rate H(a). The sound horizon decreases due to the increased energy density driving a more rapid expansion and a decrease in the sound speed, c s = [3(1 + R(a))] −1/2 . The diffusion length increases, naively, due to a decrease in the scattering rate driven by the reduced Hubble time. Caution should be taken when using such naive scaling arguments. For example, the non-trivial dependence of the recombination history leads to counterintuitive effects in theÑ eff dependence on m ν [42]. Ref. [42] demonstrates that where a naive scaling argument would suggest an increase inÑ eff with increasing m ν , the non-trivial dependence of the recombination history on m ν impliesÑ eff decreases monotonically and rapidly with increasing m ν (See Sec. 5.1).

Verification of cosmological parameters with dark radiation
We adopt the point of view advocated in Refs. [10,27,28] that the constraint provided by predictions of BBN should be incorporated simultaneously with constraints due to recombination effects in the extraction of cosmological parameters. We also require, as previously discussed, that BBN and recombination be solved iteratively. In this section, we demonstrate the self-consistent extraction of these parameters by employing a simple model for the radiation energy density that avoids solving, for example, the Boltzmann equation, for the set of distribution functions of the cosmic constituents, in particular the neutrinos. The model we explore in this section is identical to ΛCDM except for one additional constituent: dark radiation. The dark radiation energy density is radiation at all epochs and does not interact with the other energy-density constituents through any force except gravitation.
The total energy density as given in terms of radiation, matter, and vacuum energy components is which depend on the scale factor a as a −4 , a −3 , a 0 , respectively, as long as there is no energy transfer between the species. We assume for the purposes of this section that the neutrinos are massless, always acting as radiation energy density. The matter energy density consists of contributions from baryons and cold dark matter. Modeling the matter as a pressureless gas and observing that the comoving matter energy density is conserved, we write the proper energy density as The vacuum energy is the least understood of the energy densities. Assuming the universe to be critically closed, the sum of the three energy densities must be equal to the critical energy density, specified only by the Hubble rate at the current epoch. Hence, ρ v = ρ c − ρ r,0 − ρ m,0 . The vacuum energy density is negligible at all epochs of interest in this paper but is included for completeness.
The radiation component is given as: where ρ γ is the photon energy density, ρ ν is the neutrino energy density, and ρ dr is the dark-radiation energy density. We parametrize ρ dr as ρ dr = 7 4 4 11 where δ dr is the dark-radiation parameter and is always assumed to be non-negative. In principle, we could entertain negative values of δ dr since it is an adjustable parameter of ρ r . Such a change requires a fundamental reworking of the ΛCDM model so that i Ω i = 1. These non-standard cosmologies obtain when considering, for example, neutrino oscillations. These models, however, are not continuously connected with our model at δ dr = 0, for any values of the parameters, thus motivating the maintenance of δ dr > 0. We write N eff for 'standard' cosmologies. It is, therefore, unnecessary for this simple dark-radiation model, to deduceÑ eff from r s /r d since N (th) eff =Ñ eff by construction. This model of dark radiation is the usual model applied, for example in Ref. [11] and we explore, in this section, the predictions of the present burst code to verify our results against those of prior results within the community. We will useÑ eff , for the remainder of this section, to denote the "effective number of relativistic degrees of freedom 6 ." Figures 2-5 show the results of computations in which the four parameters ω b , Y P , D/H, andÑ eff are varied. We begin by varying the two model inputs: ω b andÑ eff . The upper panel in Fig. 2 shows the dependence ofÑ eff on ω b for curves of constant Y P ; the vertical band is the Planck value of ω b = 0.02207 ± 0.00033. The figure is generated by first choosing a value for the baryon density in the range 0.004 ≤ ω b ≤ 0.029. Each selected value of the dark radiation parameter in the range 3 ≤Ñ eff ≤ 4.5 allows for the prediction of Y P and D/H by parametrizing the radiation energy density as in Eq. (2.12). The values so obtained are plotted as contours in theÑ eff -ω b plane in the upper and lower panels of Fig. 2. The solid curve is the preferred value Y P = 0.2465 ± 0.0097 of Ref. [7], which is a selection of observations of metal poor extragalactic H ii regions. The contours are spaced by roughly 0.0097/3 showing thatÑ eff is not strongly constrained by values of Y P alone; this is a manifestation of the degeneracy ofÑ eff and Y P . For example, at Y P = 0.2465 ± 0.0035, corresponding to the contours closest to Y P = 0.2465, the range allowedÑ eff is nearly consistent with both the standard, calculated value N eff = 3.046 and the Ref. [11] derived value of N eff = 3.30 ± 0.27.
Predictions of the primordial deuterium abundance are a much more sensitive constraint upon allowed values ofÑ eff . This can be seen in the lower panel of Fig. 2. The solid line contour with 10 5 × D/H = 2.530 ± 0.04 corresponds to the recent measurement of Ref. [8]. Contours in this figure are separated by the one standard deviation of Ref. [8]. There are two points of interest regarding the deuterium figure. First, as noted in Ref. [8], observation of the primordial component of deuterium is precise enough to begin to constrain the microscopic physics of the thermally averaged nuclear reaction rates and their cross sections. Additionally, given the precision of the current and forthcoming deuterium measurements and the strong dependence ofÑ eff on its value (at constant ω b ), we advocate using D/H as a prior, over Y P , for future base model parameter searches.
The degeneracy betweenÑ eff and Y P is again evident in Fig. 3. Each plot explores the D/H vs. ω b contour space, where the upper plot contains contours of constant Y P and the lower plot contains contours of constantÑ eff . The shaded bands in each figure indicate the one-sigma observations of ω b and D/H from Refs. [8,11], respectively. Deuterium is not an input parameter into our model. We compute it by choosing a baryon number and iteratively change the dark-radiation parameter, δ dr until matching the chosen deuterium target. The outputs from the process areÑ eff and Y P . Values of ω b , D/H andÑ eff are in satisfactory agreement with the standard cosmology at the precision of current observations. The quantities D/H and ω b are the tightest observationally-constrained parameters we are currently investigating. Figure 4 shows two plots in the Y P -N eff plane with contours of constant ω b (upper plot), and contours of constant 10 5 × D/H (lower plot). The horizontal band in each figure indicates the one-sigma observation of Y P from Ref. [7]. Like D/H, Y P is not an input into our model. Consequently, we adopt the identical iterative method for Y P in Fig. 4 as we do for D/H in Fig. 3. For the ω b (upper) plot of Fig. 4, the solid contour line is the best-fit value of Ref. [11] with nine-sigma spacing of the contours. The contours exist in a subspace of the Y P -N eff plane which is well within current observations, but nevertheless could span a range of radically different physics. Similarly, the bottom plot shows the   Y P plotted against ω b for contours of constantÑ eff . The contours are spaced by ∆Ñ eff = 0.33. Figure 5 shows howÑ eff changes in the Y P -ω b plane. The shaded bands in each figure indicate the one-sigma observations of ω b and Y P from Refs. [7,11], respectively. We may conclude from this plot that there could exist many different values ofÑ eff consistent with the observations of Y P and ω b . However, we caution against such a conclusion without considering the effect on D/H. In fact, we choose not to include a figure of contours of D/H in the Y P -ω b plane because the strong sensitivity of D/H to ω b produces too large a range of values for Y P to be useful as a constraint of the cosmological model.
All calculations show a consistency between ω b ,Ñ eff , D/H, and Y P to a conservative limit of two-sigma error range in each observation. We expect the uncertainties in each observation to improve in the coming years with large ground-based CMB experiments [43,44] and 30meter class telescopes [45][46][47]. Future high-precision measurements may result in tensions for the best-fit values of ω b ,Ñ eff , and D/H. These tensions could be indications of the need for more precise theoretical and numerical approaches or could signal the presence of physics beyond the standard model. As it stands here, the bottom panels of Figs. 2 and 3 shows that the tension between ω b and D/H cannot be resolved with the addition of extra radiation energy density. Uncertainties in nuclear reactions may produce disagreement between ω b and D/H, allowing D/H to become a probe of nuclear physics. It is also possible that the spectroscopic determination of D/H may be subject to small systematic errors only recognizable at such precision. An exciting prospect is the need to revise the CSM to resolve tensions with the observations of D/H and ω b , possibly leading to the conclusion of BSM physics active during BBN.

Examples for neutrino sector BSM physics
We describe three examples of unresolved issues in BSM/CSM physics, which require the fully self-consistent parameter determination described in previous sections. We consider, in turn, models incorporating neutrino rest mass, sterile neutrinos, and non-zero lepton numbers. We use, throughout this section, the observationally-inferred definition of N eff ,Ñ eff .

Neutrino rest mass
Section 4 details a self-consistent treatment of the BBN observables Y P , D/H,Ñ eff , and ω b . The sum of the light neutrino masses, denoted m ν , has no bearing on the determination of primordial abundances in BBN calculations due to the high temperatures relevant there. Here, however, we explore the epochs and energy scales in the history of the universe associated with the m ν energy scale in order to investigate the relationship between m ν and the other four observables of interest (ω b ,Ñ eff , Y P and D/H). Specific examples of such epochs that we might consider include the surface of last scattering (z ∼ 1100) and the advent of LSS (z 10). We focus on the surface of last scattering and implications for the CMB in this paper.
Reference [42] (hereafter GFKPI) investigates the effect of neutrino rest mass onÑ eff using the burst suite of codes. Conventional estimates based on the energy density added by non-zero neutrino rest masses suggest an increase in N (th) eff . However, using the method outlined in Sec. 3, GFKPI shows a decrease inÑ eff . The decrease is due to an effect on the recombination history stemming from an increase in the Hubble rate, which results in a larger freeelectron fraction. This counterintuitive result is termed the "neutrino-mass/recombination (νMR) effect." The νMR effect manifests itself only in a self-consistent treatment, such as that employed by GFKPI.
In addition, GFKPI investigates the dependence of the νMR effect on ω b . We revisit this physics here in preparation for a discussion on the effect of non-zero lepton number L ν later, in Sec. 5.3.2. Increasing ω b leads to an enhancement of the νMR effect when m ν is held constant, as is evident by the curvature of the contours in Fig. 6. The enhancement is a consequence of the effect that changing ω b has on the the recombination history. We might naively expect a larger change in the Hubble rate relative to the massless neutrino case resulting in an enhanced νMR effect for the smaller ω b case. This is opposite to that observed in Fig. 6. This result also is counterintuitive based on expectations from a simple scaling of the energy density and the resulting change in the recombination history [42].
The origin of the enhancement of the νMR effect can be understood by considering a simplification of the Boltzmann equation that determines the recombination history [Eq. (2.10)]. We take Y P = 0 for the purposes of this argument since the νMR enhancement is insensitive to Y P , as we have verified numerically for the ranges of parameters we are considering. In this simple scenario, X e = X p and we obtain an expression for the change in the free-electron fraction: where β ≡ α (2) (m e T /2π) 3/2 e −∆Q/T is the ionization coefficient and α (2) is the recombination coefficient with ∆Q = 13.6 eV. Equation (5.1) for the recombination history shows that the free-electron disappearance rate is proportional to the total electron number density, which in turn is proportional to ω b through Eq. (2.9). Equation (2.9) also shows how n (tot) e relates to Y P . Note that ω b and Y P affect n (tot) e differently. However, due to the relative insensitivity of Y P to ω b , the ω b dependence dominates in Eq. (2.9). The increase in energy density from m ν = 0 and the increase in the free-electron disappearance rate combine to alter the recombination history so as to enhance the νMR effect for increasing ω b .

Sterile neutrinos
We next consider the possibility that there exists either single or multiple sterile-neutrino species, which could have profound implications in cosmology. We entertain two possibilities of either light or heavy sterile neutrinos.

Light sterile neutrinos
Observations of neutrino events in large scintillating detectors may have revealed anomalies that could be interpreted as sterile neutrinos with rest masses m νs ∼ 1 eV [48][49][50]. We investigate the presence of a single sterile neutrino in the early universe by employing a model where the sterile state populates a thermal Fermi-Dirac shaped distribution with temperature parameter T s , possibly through flavor mixing. The sterile neutrino temperature, T s is taken to be less than or equal to the active neutrino temperature T ν . The ratio T s /T ν is assumed to be the same throughout weak decoupling, BBN, and recombination. For this analysis, we do not investigate smaller active-sterile neutrino mixing angles with resultant non Fermi-Dirac-shaped energy spectra [41,51]. Future work will consider such physics [52].  3. The deviation of the contours from the dotted lines is again due to an effect similar to the νMR effect but, in this instance, due to the sterile state. Fig. 7 takes the sum of the active neutrino masses to be 0.06 eV. This is inconsequential for large T s /T ν 1. For T s /T ν 0.1, ∆Ñ eff < 0 due to the νMR effect in the active neutrino sector. As a consequence, the contour forÑ eff = 3 is not coincident with the m νs axis. Since m νs is too small to be of any significant kinematic effect during BBN, we need only compute BBN once for a given value of T s /T ν . During recombination, m νs is kinematically important and affects the Hubble rate. Hence, for every point in the T s /T ν -m νs plane of Fig. 7 we calculate recombination. This figure clearly emphasizes the need for a self-consistent treatment between BBN and recombination when considering this BSM physics.

Heavy sterile neutrinos
Heavy sterile neutrinos that decay out of equilibrium in the early universe can affect weak decoupling and, as a consequence, primordial nucleosynthesis [16,53,54]. Sterile neutrinos in the rest mass range 0.1 GeV ≤ m νs ≤ 1.0 GeV, with lifetimes 1 s decaying during the weak decoupling, weak freeze-out, and/or BBN epochs can have constrainable, sometimes dramatic, cosmological effects.
Such sterile neutrinos have mass and vacuum mixings with ν e , ν µ , ν τ constrained by accelerator and other laboratory oscillation experiments/observations [55][56][57][58][59][60][61][62], beta-decay experiments [63], and cosmological considerations, including constraints on m ν and N eff [64][65][66][67][68][69]. In fact, stringent constraints can be obtained from N eff limits alone [16], as sterile neutrinos decaying out of equilibrium can lead to dilution (entropy production) which, in the weak decoupling epoch, can lead to distortions in the relic neutrino energy spectrum, affecting m ν , and have significant impact on the relativistic energy content and, hence, N eff . A sophisticated theoretical and computational treatment of dilution physics is a challenging endeavor. Even though these heavy sterile neutrinos may decay away before an epoch where T ∼ 10 keV, they can nevertheless alter the relationship between Y P , D/H, N eff , m ν , and ω b , necessitating the need for a self-consistent treatment between the weak decoupling, weak freeze-out, BBN, recombination, photon decoupling, and advent of LSS epochs [70].

Lepton numbers
We examine how lepton numbers affect the primordial abundances andÑ eff . We define the lepton number L ν for a neutrino species ν in a flavor eigenstate as where n ν is the number density of neutrino species ν, nν is the number density of antineutrino speciesν, and n γ is the number density of photons. Here, for illustrative purposes, we take L νe = L νµ = L ντ ≡ L ν . In fact, flavor oscillations will bring the various lepton numbers into near equality [71][72][73].
We use the comoving-invariant neutrino degeneracy parameter ξ ν ≡ µ ν /T ν , where µ ν is the chemical potential of neutrino ν, to compute L ν . The present model assumes that the lepton number evolves through the epoch of e ± annihilation only in response to the relative increase in n γ . We do not consider any BSM physics which could alter the difference n ν − nν; that is, we fix ξ ν throughout weak decoupling, BBN, and photon decoupling. We relate the degeneracy parameter to the lepton number using the following expression [51]:

Effect on nucleosynthesis
The helium mass fraction is sensitive to the neutron-to-proton ratio, n/p. We determine n/p by calculating the weak rates associated with neutrino-nucleon reactions, namely: though unlike the brief time/temperature range of α-formation, weak freeze-out occurs over several Hubble times at this epoch. The rates in reactions (5.4) through (5.6) are sensitive to the neutrino and e ± distributions. We follow Ref. [21] to evolve T and the electron chemical potential in order to maintain equilibrium between the electrons, positrons and photons. For the electron-flavor neutrinos, we use the comoving invariants aT ν and ξ νe to compute the neutrino distributions. We set m ν = 0 as neutrinos of sub-eV rest mass remain ultrarelativistic throughout weak freeze-out. Figures 8 and 9 show the helium mass fraction and the relative deuterium abundance, respectively. Each plot is in the L ν -ω b plane for contours of constant primordial abundance. The relationships between lepton number and nucleosynthesis are well known [74,75]. Increasing L νe leads to an overabundance of neutrinos compared to anti-neutrinos. The forward rate of reaction (5.4) freezes-out after the forward rate of reaction (5.5). The imbalance lowers n/p which lowers Y P as seen in Fig. 8. The decrease in n/p also leads to a decrease in D/H, although deuterium is not as sensitive to L ν as helium. However, D/H is known to much higher precision than is Y P .
Comparing with recent observations [7,8], the two light element abundances achieve consistency at 2σ. Y P prefers a value of L ν < 0 whereas D/H prefers a positive value of L ν . If future observations of the light-element abundances were to show a larger disagreement than 2σ, lepton numbers of identical value could not solely rectify the tension. Future analyses will consider scenarios with multiple facets of BSM physics including non-zero lepton numbers [76]. This analysis will useÑ eff as a discriminating factor.

Effect on N eff
We consider how two aspects of non-CSM/BSM physics (L ν = 0 and/or m ν = 0) modifỹ N eff . If we set the neutrino rest mass to m ν = 0, we can investigate whether the νMR effect still applies with non-zero L ν . Figure 10 shows the changes toÑ eff in the m νω b plane for three values of L ν = −0.05, 0, 0.05 corresponding to blue, magenta, and red contours, respectively. The non-zero lepton number increasesÑ eff for small values of m ν in accordance with Ref. [9]. However, for values of m ν ∼ 1.0 eV, the νMR effect overwhelms the extra energy density from more particles to lowerÑ eff below three. Figure 10 also shows that, despite the total energy density being insensitive to the sign of L ν , the contours for non-zero values of L ν with opposite sign do not overlap because Y P depends sensitively on its value. Y P is largest for the blue contours, so more helium suppresses the νMR effect.
Note that taking L ν = 0 conflates the interpretation that the effect of m ν is identical to perturbations in the matter power spectrum, which are used in calculating the suppression of power on small scales. This is borne out by the present model where the m ν statistic cannot be equated to the cosmological measurement. In our model, m ν is simply the sum of the active vacuum neutrino mass eigenvalues. The observationally determined value of m ν depends on quantities other than the sum of the active neutrino masses such as their energy distributions.

Conclusions and outlook
Cosmological considerations are a key route to exploring BSM physics. This is especially true for the neutrino sector, where there are many outstanding questions and where laboratory experiments are limited in what aspects of this physics can be addressed. In this paper we have argued that a self consistent treatment of BSM issues, across all epochs from weak decoupling to photon decoupling, is the best way to take advantage of the expected coming increase in precision of CMB measurements and observationally-inferred primordial abundances of the light elements. We employ a limited prescription to link the salient features of self consistency between early-time neutrino dynamics and the surface of last scattering. We couple the weak decoupling and nucleosynthesis of early times to CMB observables, including baryon-to-photon ratio (equivalently ω b ), sound horizon, and photon diffusion length.
We have shown that such a self consistent treatment is necessary, in part because new neutrino physics can alter the relationships between different cosmological epochs. For example, ω b and other CMB observables affect the calculated yields of deuterium and helium. In addition, the calculated relic neutrino energy spectra after weak decoupling affects the predicted value ofÑ eff at photon decoupling. The principal tool in our analysis is the suite of burst codes for nucleosynthesis and neutrino interactions and energy transport.
A case in point is our investigation of the relationship between neutrino rest masses, i.e. m ν , and the four potential observables ω b ,Ñ eff , Y P , and D/H. This analysis reveals the "neutrino-mass/recombination" (νMR) effect first described in Ref. [42]. The νMR effect is below the threshold of current CMB capabilities, but may not be in future observations [77].
There have been spectacular advances in the measurements/observations of neutrino properties. We know the neutrino mass-squared differences and three of the four parameters in the unitary transformation between the energy eigenstates (mass states) and the weak interaction eigenstates (flavor states) of the three active neutrinos (only the CP -violating phase remains unmeasured). As active neutrinos mix in vacuum and have non-zero rest masses, the question indubitably arises of whether there exist "sterile" neutrino states. If indeed sterile neutrinos do exist, we acknowledge that the parameter space of mass, vacuum mixing angle, and number is enormous. However, sterile neutrinos could have profound effects in all of the epochs under study in this paper. This possibility makes a self consistent treatment of these effects a powerful basis for constraining sterile neutrino states.
In this paper we have considered scenarios for both "light" (mass ∼ 1 eV) and "heavy" (mass ∼ 0.1 − 1 GeV) sterile neutrinos. In the former case we consider cases where the sterile neutrino relic energy spectra are Fermi-Dirac black-body shaped, though with a temperature parameter T s differing from that characterizing the relic energy spectra of active neutrino species. We show here that the νMR effect has interesting consequences and that this case demands a self consistent treatment of recombination and BBN. Additionally, heavy sterile neutrino decay out of equilibrium can lead to dilution and high energy relic active neutrinos, and both of these features potentially can have dramatic and constrainable effects on CMB-epoch observables. This implies that CMB observations can indirectly probe the CνB and explore active-sterile mass/vacuum-mixing parameter space unavailable to current accelerator-based experiments.
We have also studied the effects of non-zero lepton numbers on the relationship between CMB observables, nucleosynthesis, and neutrino physics. Our conclusion is that the primordial deuterium abundance is a potentially powerful probe of lepton number. However, an eventual CMB-only measurement of the primordial helium abundance Y P will be the most powerful probe of lepton number and many other issues. Determining Y P from CMB observables will require a sophisticated self-consistent approach to BBN, neutrino physics, and photon decoupling transport physics.
Finally, our study has revealed a potential tension betweenÑ eff , ω b , and the primordial deuterium abundance, D/H, inferred from high redshift QSO absorption systems. In fact, if the advent of 30-m class telescopes in the near future allows for a decrease in errors in observationally-inferred D/H to the ∼ 1 % level, while observed ω b andÑ eff maintain their respective current central values, then tension is unavoidable. This may signal BSM or non-CSM physics, likely in the neutrino sector, or it could point to not understanding systematics in the damped Lyman-α cloud measurements of the isotope-shifted hydrogen absorption lines. We advocate using future instruments to explore the rich physics of weak decoupling, nucleosynthesis, and photon decoupling to discover what role BSM neutrino physics has in these epochs.