Very different processes characterize the decoupling of neutrinos to form the cosmic neutrino background CB and the much later decoupling of photons from thermal equilibrium to form the cosmic microwave background CMB. The CB emerges from the fuzzy energydependent neutrinosphere and encodes the physics operating in the early universe in the temperature range T  10 MeV to T  10 keV. This is the epoch where beyond Standard Model BSM physics especially in the neutrino sector may be influential in setting the light element abundances the necessarily distorted fossil neutrino energy spectra and other light particle energy density contributions. Here we use techniques honed in extensive CMB studies to analyze the CB as calculated in detailed neutrino energy transport and nuclear reaction simulations of the protracted weak decoupling and primordial nucleosynthesis epochs. Our moment method relative entropy and differential visibility approach can leverage future high precision CMB and light element primordial abundance measurements to provide new insights into the CB and any BSM physics it encodes. We demonstrate that the evolution of the energy spectrum of the CB throughout the weak decoupling epoch is accurately captured in the Standard Model by only three parameters per species a nontrivial conclusion given the deviation from thermal equilibrium and the impact of the decrease of electronpositron pairs. Furthermore we can interpret each of the three parameters as physical characteristics of a nonequilibrium system. Though the treatment presented here makes some simplifying assumptions including ignoring neutrino flavor oscillations the success of our compact description within the Standard Model motivates its use also in BSM scenarios. We further demonstrate how observations of primordial light element abundances can be used to place constraints on the CB energy spectrum deriving response functions that can be applied for general deviations from a thermal spectrum. Combined with the description of those deviations that we develop here our methods provide a convenient and powerful framework to constrain the impact of BSM physics on the CB.
['big bang nucleosynthesis', 'cosmological neutrinos', 'particle physics - cosmology connection', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'High Energy Physics - Phenomenology', 'High Energy Physics - Theory']
Cosmic neutrino decoupling and its observable imprints insights from entropicdual transport
{'Cosmic Neutrino Decoupling and its Observable Imprints: Insights from Entropic-Dual Transport': '- J. Richard Bond, a George M. Fuller, b Evan Grohs, c Joel Meyers, d and Matthew James Wilson e\n- a Canadian Institute for Theoretical Astrophysics, Toronto, Ontario M5S 3H8, Canada\n- b Department of Physics, University of California, San Diego, La Jolla, California 92093, USA\n- c Department of Physics, North Carolina State University, Raleigh, NC 27695, USA\n- d Department of Physics, Southern Methodist University, Dallas, TX 75275, USA\n- e Institute for Astroparticle Physics, Karlsruhe Institute of Technology, Karlsruhe, 76131, Germany', 'Abstract': 'Very different processes characterize the decoupling of neutrinos to form the cosmic neutrino background (C ν B) and the much later decoupling of photons from thermal equilibrium to form the cosmic microwave background (CMB). The C ν B emerges from the fuzzy, energy-dependent neutrinosphere and encodes the physics operating in the early universe in the temperature range T ∼ 10 MeV to T ∼ 10 keV. This is the epoch where beyond Standard Model (BSM) physics, especially in the neutrino sector, may be influential in setting the light element abundances, the necessarily distorted fossil neutrino energy spectra, and other light particle energy density contributions. Here we use techniques honed in extensive CMB studies to analyze the C ν B as calculated in detailed neutrino energy transport and nuclear reaction simulations of the protracted weak decoupling and primordial nucleosynthesis epochs. Our moment method, relative entropy, and differential visibility approach can leverage future high precision CMB and light element primordial abundance measurements to provide new insights into the C ν B and any BSM physics it encodes. We demonstrate that the evolution of the energy spectrum of the C ν B throughout the weak decoupling epoch is accurately captured in the Standard Model by only three parameters per species, a non-trivial conclusion given the deviation from thermal equilibrium and the impact of the decrease of electron-positron pairs. Furthermore, we can interpret each of the three parameters as physical characteristics of a non-equilibrium system. Though the treatment presented here makes some simplifying assumptions including ignoring neutrino flavor oscillations, the success of our compact description within the Standard Model motivates its use also in BSM scenarios. We further demonstrate how observations of primordial light element abundances can be used to place constraints on the C ν B energy spectrum, deriving response functions that can be applied for general deviations from a thermal spectrum. Combined with the description of those deviations that we develop here, our methods provide a convenient and powerful framework to constrain the impact of BSM physics on the C ν B.', '1 Introduction': "Our focus in this paper is to provide a new way of analyzing the cosmic neutrino background (C ν B). The 'fingerprints' of the C ν B can be found on the Cosmic Microwave Background (CMB) and large scale structure as well as on the primordial light element abundances which emerge from the epoch in the early universe when the weak interactions involving neutrinos and the nuclear reactions freeze out. Revealing these fingerprints may leverage future experimental and observational results into new insights of the early universe and the physics operating therein. There is urgency about this effort, as observation and experiment are poised to provide higher precision data in cosmology and on physics in the neutrino sector [1-7], while the origin of neutrino mass is a frontier issue in elementary particle physics [8].  \nNext generation ground-based cosmic microwave background (CMB) experiments such as CMB-S4 [9, 10] promise to measure N eff (a measure of the energy density in particles with relativistic kinematics around the epoch of photon decoupling) with a precision of about one percent, place tight constraints on the primordial helium abundance, and provide a significant measurement of the 'sum of the light neutrino masses' (a measure of the suppression of matter clustering associated with massive cosmic neutrinos [11-16]). In parallel, laboratory- and accelerator-based experiments promise to elucidate key issues in neutrino physics [2]. Within five years, long baseline neutrino oscillation experiments will give us the neutrino mass hierarchy to fair confidence [17]. Likewise, we can expect important neutrino rest mass constraints from front line and future tritium endpoint experiments (e.g., KATRIN [18], Project-8 [19]), and tonne-scale neutrino-less double beta decay experiments (e.g., nEXO [20]) that may pin down the character (Majorana or Dirac) of neutrinos and provide neutrino mass and hierarchy data complementary to that derived from cosmology. A known neutrino mass hierarchy not only extends the utility of the CMB- and large scale structure-derived constraints on, or measurements of the sum of the light neutrino masses, but it also feeds back on the results of the laboratory neutrino mass probes [2].  \nThe existence of the C ν B is a central prediction of hot big bang cosmology. Within the Standard Model of particle physics, assuming a standard thermal history of our universe, we can make very definite predictions about the properties of cosmic neutrinos. Up to some small corrections, the Standard Model predicts that cosmic neutrinos and antineutrinos of all flavors were well described by a relativistic, nearly zero degeneracy parameter Fermi-Dirac (FD) distribution. In this simple and approximate description, these distributions share a common temperature with photons when the temperature is high. At lower temperatures when the electron-positron pairs have disappeared, co-moving entropy-conservation yields a neutrino-to-photon-temperature ratio of (4 / 11) 1 / 3 . Subsequent to decoupling, neutrinos free-fall and their momenta redshift with the expanding universe, implying that some may have non-relativistic kinematics at the current epoch [11]. In what follows, we will use the results of detailed neutrino energy transport calculations to characterize deviations from this simplistic picture.  \nCosmic neutrinos made up a sizable fraction of the energy budget of the early universe. The gravitational influence of cosmic neutrinos had a significant impact on the expansion of the universe and the evolution of structure during the radiation-dominated era. Observations of primordial light element abundances, the angular power spectra of the CMB, and the power  \nspectrum of large-scale structure are therefore sensitive to the presence of cosmic neutrinos. These observations are most sensitive to the total energy density of the cosmic neutrino background, parameterized through N eff , defined such that the total radiation energy density of the universe is  \nρ r = ρ γ ( 1 + 7 8 ( 4 11 ) 4 / 3 N eff ) , (1.1)  \nwhere ρ γ is the energy density of photons. The Standard Model prediction is N eff = 3 . 044(1) [2123]. Current observational constraints on N eff are consistent with the predictions of the Standard Model; data from the Planck satellite provides a measurement N eff = 2 . 99 ± 0 . 17 (68% CL) [24] in ΛCDM+ N eff cosmology. In the same model, measurements of the primordial helium and deuterium abundances along with BBN calculations suggest a constraint N eff = 2 . 89 ± 0 . 23 [25], and CMB combined with helium and deuterium measurements gives N eff = 2 . 90 ± 0 . 14 [25].  \nUpcoming CMB experiments will greatly improve the precision with which N eff can be measured [3, 9, 10, 26-28]. On the other hand, future CMB observations by themselves are unlikely to give much insight into more detailed properties of cosmic neutrinos such as neutrino flavor content or the energy spectrum. However, when coupled with Big Bang Nucleosynthesis (BBN) and primordial light-element abundance determinations [29], we can gain insights into these issues. BBN depends not just on the gravitational influence of cosmic neutrinos but also on their weak interactions - with exquisite sensitivity [30]. Primordial abundances therefore in principle provide a window into measuring properties of cosmic neutrinos which are inaccessible through other observations.  \nPhysics Beyond the Standard Model (BSM) can cause the cosmic neutrinos to have an energy distribution function which differs significantly from a relativistic FD distribution [31, 32]. For example, the decay or annihilation of some new particles could inject additional neutrinos compared to the standard thermal history, or oscillation into new sterile states could alter the neutrino distribution [33]. If these processes occur at sufficiently early times, the neutrino distribution will tend to relax into thermal equilibrium. However, late changes to the neutrino distribution may become frozen-in leading to an altered spectrum at late times, perhaps accompanied by alterations in N eff and light element abundances relative to the Standard Model. New light degrees of freedom may also contribute to N eff and impact associated observables [4].  \nLarge-scale numerical simulations have been employed to follow this evolution of the neutrino component and the associated Standard Model baryonic physics like nuclear freeze-out and BBN [21-23, 34-45]. However, these coupled neutrino transport and baryonic physics problems are not simple to model with Boltzmann calculations, and will be even more difficult to model with general treatments of neutrino flavor quantum coherence. The heart of this paper is the introduction of a technique for capturing and physically interpreting the salient features in these large-scale numerical simulations. Our goal for the present work is not to directly replace these simulations, but rather to introduce a framework that can be used to describe the relevant physics and which may be further developed to compactly treat more general BSM scenarios. As such, we leave out some effects, such as neutrino flavor oscillations, that are included in other treatments of C ν B transport that are aimed at precision calculation of N eff in the Standard Model, for example.  \nWe demonstrate that we can capture the essential aspects of neutrino transport in the early universe through use of generalized statistical techniques. The generalized entropy approach that we develop allows us to faithfully represent the evolution of the cosmic neutrino background, as derived from detailed Boltzmann transport, using just three parameters per species. We show how the distribution function of each neutrino species can be described at high precision as n ( t ) = [exp( α ( t )) + 1] -1 , with the function α = α 0 + α 1 ϵ + α 2 ϵ 2 , where ϵ is related to the energy of the neutrinos, to be defined more precisely below. The two parameters α 0 and α 1 are familiar in the context of the Grand Canonical Ensemble (GCE) from statistical mechanics, and relate to average number and average energy. On the other hand, the third parameter is a new aspect of the neutrino distributions. It is conjugate to the average of the energy-squared of the distribution, and by extension is related to the variance in energy.  \nIn general the total energy-dependent quantity α is conjugate to the Wigner distribution function matrix of Boltzmann transport theory, and fully encodes the transport information, as the response in entropy due to a change in the distribution of number. This duality is general for all transport networks. For example, in our nuclear reaction network, each nuclide J is described by a chemical potential µ J , with a common temperature T . The abundance is conjugate to a chemical potential, which would give a Maxwell-Boltzmann distribution of energies, and hence entropy, when used in tandem with the common temperature. The nuclear reaction network we use in this work follows abundances, but it is useful to consider how it maps into the { α J } dual description: the α J give the response of the entropy to changes in the number of J nuclei. If the nuclear reactions which connect a subset of the nuclide network are rapid, equilibrium is locally attained and the total chemical potentials of the products and reactants in the subset are in detailed balance, with entropy differences related to the latent heat, i.e., the nuclear binding energies. At temperatures well above 1 MeV, the system is in chemical and weak equilibrium, so all nuclide abundances are describable in terms of a single chemical potential related to the baryon density. As the temperature decreases, more chemical potentials ( i.e. , α J ) are needed, first in quasi-nuclear statistical equilibrium among a subset of nuclides. The partition becomes finer with decreasing reaction rates until eventually every nuclide requires its own α J in an out-of-equilibrium configuration. As in the neutrino network case of distribution function numerical solutions, our nucleosynthesis network is solved numerically in abundances, but the dual description gives insight into such chemical equilibrium relations. The total response in the nuclear entropy for a change in baryon number is ∑ J ⟨ α J ⟩ Y J , where ⟨ α J ⟩ is the energy-averaged value of α J and Y J is the abundance. Similarly, the total response in the neutrino entropy is a weighted sum over the individual responses for each energy and flavor. If oscillations are ignored, the ν µ -ν τ neutrino flavor blocks are related by the same microscopic rates and have the same responses (likewise for the ν µ -ν τ blocks). The distribution function progression we describe below as the temperature drops is akin to a full neutrino equilibrium needing one α 1 , i.e. , one neutrino temperature, then a α 1 for temperature and α 0 for number, then an additional α 2 . To complete the network description in α terms, in the photon-electron-positron plasma, the equilibrium relation µ e -+ µ e + = 2 µ γ = 0 is obeyed, implying µ e + = -µ e -: the γ -e --e + plasma reduces to a full statistical equilibrium description, with the asymmetry in α e ± encoding the residual electron number per baryon once the e + e -→ γγ phase transition is complete.  \nThis technique applied to the neutrino network is a simple-to-use parameterization which  \ncaptures the most relevant physics and can offer cogent physical interpretations of key features emerging from the numerical simulations. It allows us to understand the origin of the deviations from the thermal spectrum in the Standard Model and BSM scenarios. That, in turn, helps us to quantify the degree to which observations could be sensitive to the energy distribution of cosmic neutrinos. Using our method, we can describe a 'last scattering surface' for cosmic neutrinos in close analogy with the CMB, finding a broad peak with energy-dependent width around z ∼ 10 10 (corresponding to a temperature T ∼ 1 MeV). We show how primordial abundances respond to changes in the cosmic neutrino energy spectrum, and thereby show how abundance measurements and CMB parameters can be used to place constraints on physics beyond the Standard Model.  \nIn Section 2, we give an overview of the physics of weak decoupling and weak freeze-out. Section 3 describes the general framework of neutrino transport in the language of Wigner functions and their conjugates. We compare aspects of neutrino decoupling to CMB transport in Section 4. Section 5 discusses the weak interactions relevant to neutrino decoupling. We describe the numerical implementation of neutrino transport and primordial nucleosynthesis in Section 6. Section 7 describes neutrino transport in the Standard Model and develops a compact treatment of how the occupation probability of neutrinos deviate from equilibrium. In Section 8, we introduce the neutrino differential visibility and productivity as a novel means to visualize the neutrino decoupling process. We discuss how weak nucleon interactions affect primordial abundance yields in Section 9. Section 10 introduces the linear response of light element abundances to distortions of the cosmic neutrino spectrum. We conclude in Section 11. In what follows, we use natural units where c = ℏ = k B = 1.", '2 Overview of Weak Decoupling and Weak Freeze-Out': 'In this section we review the process of weak decoupling and weak freeze-out in the Standard Model, assuming a standard thermal history, following Ref. [44]. We will describe how out-ofequilibrium processes cause small distortions in the cosmic neutrino spectrum even in the standard case. In the next section we introduce a full parameterization of these distortions in terms of a thermodynamic potential conjugate to the general neutrino distribution functions. We will show how this distribution function can be represented in an expansion in terms of the coefficients { α ( f ) } , with f denoting a flavor label. Surprisingly, with very few of these coefficients, our procedure can reproduce the results of detailed transport calculations for C ν B decoupling. This conjugate description can illuminate when and how the C ν B in the Cosmological Standard Model, the CSM ν , is altered by BSM physics. Figure 1 shows how the transport-derived distribution functions evolve away from Fermi-Dirac character as the universe expands and cools. We will show the subsequent sections how we can compactly describe the evolution of the CSM ν , accurately accounting for this deviation.', '2.1 Neutrino Last Scattering': 'In this section we focus on the transition between the times of frequent neutrino scattering and free expansion, defining an energy-dependent region of last scattering for cosmic neutrinos. As described in the previous section, weak decoupling is a gradual process, and we will therefore show that the neutrino last scattering surface is much broader than its CMB analog.  \nFigure 1 . Difference between neutrino distribution functions and a relativistic Fermi-Dirac distribution with vanishing chemical potential plotted against ϵ = E/T cm at epochs labeled by T cm . E is the neutrino energy and T cm is the comoving temperature quantity. Here, the neutrino distribution function for flavor f is parameterized as n ν f = (exp[ α ( f ) ( ϵ, T cm )] + 1) -1 . For this figure, f = e . A relativistic FermiDirac distribution with vanishing chemical potential is described by α eq = ϵ , and represented by the dashed black line. Here we show how the evolution of the ν e distribution function resulting from a numerical transport calculation, described by α burst , differs from the Fermi-Dirac distribution as neutrinos go out of equilibrium. Note that in this and similar figures below, we label only every other line with its corresponding value of T cm ; the unlabeled lines should be understood to have T cm given by the geometric mean of that of the neighboring lines. The Density Fraction color contours indicate the neutrino occupation fraction as a function of ϵ .  \n<!-- image -->  \n/epsilon1  \nIn the early universe during the weak decoupling and BBN epochs, SM neutrinos interact only via weak interactions. At tree-level, the scattering cross sections are proportional to the square of the Fermi constant multiplied by the square of the energy  \nσ weak ∝ G 2 F E 2 = G 2 F ϵ 2 T 2 cm , (2.1)  \nwhere G F ≃ 1 . 166 × 10 -11 MeV -2 and E is an energy-scale typical of the weakly-interacting particles. We have used Eq. (2.1) to define two quantities which we will return to frequently throughout this work. T cm is the comoving temperature quantity and is an energy scale which redshifts with increasing scale factor. It replaces the conventional neutrino temperature scale which is no longer applicable in the presence of spectral distortions. Using T cm , we define ϵ ≡ E/T cm as an alternative to E . For massless particles, ϵ is a comoving invariant as E and T cm redshift in the same manner. The cross section in Eq. (2.1) increases with increasing energy, implying that more energetic interactions will maintain equilibrium kinetics to lower temperatures. Furthermore, we can approximate the scattering rates of the neutrinos as the weak cross section multiplied by the neutrino number density  \nΓ weak ∼ σ weak n ν ∼ G 2 F T 5 , (2.2)  \nwhere we have used the ultra-relativistic kinematics of n ν ∝ T 3 and taken the typical energy scale also to be proportional to T . When compared to the Hubble expansion rate for radiation-  \ninated conditions ( H ∝ T 2 ), we find  \nΓ weak /H ∝ T 3 . (2.3)  \nNote that we will always approximate neutrinos as ultra-relativistic during weak decoupling, an approximation that is well-justified by current upper bounds on the sum of neutrino masses from beta decay endpoint measurements by KATRIN [46] and also from cosmological constraints [24]. Alternatively, if neutrinos were massive, and these masses were comparable to or smaller than the temperature, their number densities would include a Boltzmann factor e -m/T - a much more sensitive function of temperature than the relationship in Eq. (2.3). Combined with the fact that the cross section is energy dependent, neutrinos with different energies decouple at different times. Neutrino decoupling is a protracted process and the surface of last scattering is broad; see Fig. 2.  \nThe process of CMB photon decoupling differs significantly from cosmic neutrino decoupling. Photons remain in thermal equilibrium with the plasma through Thomson scattering with free electrons, since the typical photon energy during the relevant period ( ≲ 1 eV) is small compared to the electron rest mass, m e . The number density of free electrons drops precipitously when the formation of a neutral gas becomes energetically favorable (which occurs when the number density of photons with energy exceeding the atomic binding energy becomes small compared to the number density of nuclei). The ionization fraction, x e , is governed during the early stages of recombination by the Saha equation  \n1 -x e x 2 e ∝ η ( T m e ) 3 / 2 exp( B H /T ) , (2.4)  \nwhere η is the baryon-to-photon ratio and B H = 13 . 6 eV is the binding energy of hydrogen. A more careful treatment of recombination shows a somewhat slower decline in the later stages of recombination and a larger residual ionization fraction than predicted by the Saha equation [47-52]. Small deviations from a purely thermal blackbody spectrum are expected due to dissipation of acoustic waves and photon interactions with adiabatically cooling electrons and baryons, though the fractional change to the intensity is expected to be O (10 -8 ) in standard cosmology [53, 54]. Additionally, the frequency-dependent Rayleigh scattering of CMB photons on neutral species leads to spectral distortion of CMB anisotropies, but does not lead to a change in the average frequency spectrum of the CMB [55-59]. The end result is well-described by a sharp and energy-independent surface of last scattering and a nearly perfect blackbody spectrum for cosmic photons [60]. This is in stark contrast to the protracted and energy-dependent neutrino decoupling process.', '2.2 Neutron-to-proton interconversion': 'In our use of the term neutrino/weak decoupling, we describe the freeze-out of the weak interactions between neutrinos and charged leptons, and also between neutrinos and themselves. Very importantly, we do not group the baryon-isospin-changing weak-interaction processes into neutrino decoupling. Although the baryons and neutrinos interact via both charged-current and neutral-current reactions, the baryon density is so small that neutrino-baryon reactions have a  \nFigure 2 . The broad energy-dependent last scattering surface of the cosmic neutrino background can be seen from this plot of the differential visibility function for electron-type neutrinos, namely, V ν e = d ( e -τ νe ) /d ln a for scale factor a and optical depth τ ν e (see Sec. 8 for more details). Contours of V ν e are shown in the ϵ versus T cm space. The black dot-dash line gives the peak value of V ν e , and the dashed pink (white) give the peak values of the energy (number) density for massless fermions. This should be contrasted with the last scattering surface of the cosmic microwave background, which is more localized in temperature and is essentially independent of the photon frequency.  \n<!-- image -->  \nnegligible effect on the dynamics of neutrino decoupling. Conversely, the neutrino fluences are so large that they can maintain weak-chemical equilibrium among the baryons until late times. For the specific process of baryons going out of weak-equilibrium, we use the term weak freeze-out to describe the breakdown of chemical equilibrium. Although weak freeze-out is not simultaneous with the process of neutrino decoupling, they are both protracted processes that occur concurrently.  \nThe charged-current, isospin-changing, weak interactions that we consider are  \nν e + n ↔ p + e -, (2.5)  \nn ↔ p + e -+ ¯ ν e . (2.7)  \ne + + n ↔ p + ¯ ν e , (2.6)  \nThe three forward reactions yield neutron destruction, or conversely proton production, and viceversa for the reverse reactions. For brevity, we call this set of 6 isospin-changing reactions the neutron-to-proton interconversion reactions.  \nSimilar to the scattering reactions between neutrinos and charged leptons, the interconversion rates are fast at high temperature and are able to maintain a steady-state equilibrium between the number densities of neutrons and protons. The cross section for these interconversion processes scales roughly like that of Eq. (2.1), but with an important caveat. Neutrons are heavier than protons. Consequently, forward reactions of Eqs. (2.5) - (2.7), the neutron destruction channels, have no energy threshold. Correspondingly, there is an energy threshold for the reverse reactions  \nthat create neutrons from protons. The upshot is that as the universe expands and cools (absent electron neutrino degeneracy), the competition between the forward and reverse of Eqs. (2.5) (2.7) leads to an excess of protons over neutrons. Phase space considerations imply that the cross sections for the reactions in Eqs. (2.5) - (2.6) are larger than those of neutrino-electron scattering for the energy scales characteristic of weak decoupling and BBN. The energy scale for neutronto-proton interconversion is a function of the neutron-proton mass difference, and so the cross section is much larger than the neutrino scattering rates at BBN energy scales. The implications are two-fold. First, the larger cross section implies larger rates and later epochs for freeze-out. Second, the scaling relation between Γ and H is an even more gradual function in temperature than the analog in Eq. (2.3). The end result is another protracted period of decoupling: the weak freeze-out epoch.', '2.3 Weak Decoupling, Weak Freeze-Out, and Precision Observables': 'The processes of weak decoupling and weak freeze-out have implications for cosmological observables. Weak decoupling determines the number density, energy density, and energy spectrum of cosmic neutrinos. Gravitational influence of the cosmic neutrinos impacts the expansion history and the evolution of density perturbations (especially during the radiation-dominated era when neutrinos make up a significant fraction of the total energy density) and also the formation of cosmological structure (particularly at late times when neutrinos are non-relativistic). Weak freeze-out is intimately related to the process of BBN and thereby determines the primordial abundance of light elements. Precise prediction of these observables thus requires a precise treatment of these events. Furthermore, weak decoupling and weak freeze-out overlap in time, and so it is not possible to treat them as separate, isolated processes. Instead, one requires self-consistent modeling of neutrino transport throughout the whole epoch.  \nDetailed calculation of neutrino transport in the early universe is a non-trivial task even within the Standard Model where all relevant interactions are well-understood in principle [2123, 44, 61, 62]. Neutrinos interact not only among themselves but can also exchange energy with electrons and positrons in the plasma. Interaction with the plasma is complicated by the fact that there exist charged-current channels for electron-type neutrinos that are not present for the ν µ and ν τ . Neutrino decoupling occurs at temperatures around m e , such that electron-positron annihilation takes place while neutrinos are going out of equilibrium. Flavor oscillations among the neutrinos adds yet more complexity.  \nEven in the Standard Model, the neutrino decoupling process with flavor evolution included is a vexing problem. Mean-field quantum kinetic approaches to neutrino flavor evolution, as well as many-body alterations to the mean field picture, are frontier problems in computational physics. There has been ongoing rapid development of numerical techniques for solving quantum kinetic equations (QKEs) [63-65] and accompanying physical understanding of this problem [2123, 43, 62, 66]. We will neglect flavor oscillations in the analysis presented here, though the framework we develop can be extended to accommodate this additional physics.  \nThe problem becomes even more challenging when one introduces physics beyond the Standard Model. New ingredients can alter the process in a myriad of complicated ways. For example, cosmological asymmetry in lepton number [67-69], flavor oscillations among active and sterile states [33, 70-81], out-of-equilibrium decays of massive particles [31, 32, 82-86], neutrino self-  \ninteractions [75, 76, 80, 87-91], or scattering of neutrinos with new species [92-95] could alter early universe neutrino transport. The various observational handles on the cosmic neutrino background provided by primordial light element abundances, the cosmic microwave background, and large-scale structure can therefore be used to explore and constrain this new physics, but doing so requires precise calculation of neutrino transport in the early universe.  \nIn the next section, we introduce a formalism that is well-suited to handle the complexities of neutrino transport in the early universe. This formalism, a generalized entropy approach, is built on Wigner operator conjugate-control variables. As we will demonstrate, this approach can capture the salient features of the phase space evolution of the cosmic neutrinos. In later sections, we show how this generalized entropy approach allows for a compact yet highly precise description of the rich physics of neutrino transport in the Standard Model. The techniques we employ are flexible enough to accommodate a wide range of phenomena expected in models of physics beyond the Standard Model.', '3.1 Motivation and Grand Canonical Ensemble': "This section describes a generalized entropy approach to neutrino transport in the early universe. We use the term 'neutrino transport' to denote a time-evolution calculation of the six neutrino energy-momentum distribution functions in the Standard Model ( ν e , ¯ ν e , ν µ , ¯ ν µ , ν τ , ¯ ν τ ). In a full mean field quantum kinetic equation (QKE) treatment this is a daunting problem. Our goal here is to develop a complementary approach, one that is based on generalized statistical techniques and that has been eminently successful in modeling the physics at the photon decoupling epoch. The generalized entropy approach we develop, and hone on Boltzmann neutrino transport calculations, gives a description of the dynamically evolving system in the early universe based on physical constraints. The success of this treatment and the physical insights into neutrino decoupling it yields, all of which we present in subsequent sections, suggest that it could be applied to more sophisticated fully quantum mechanical neutrino transport calculations in future work. In addition, we anticipate an extension of this treatment to beyond the Standard Model (BSM) scenarios. To begin, we outline a generalized entropy approach that is inspired by full quantum mechanical many-body system calculations, but is applied here in a classical statistical guise.  \nAs described above in Sec. 2, the neutrino component at very high temperature in the early universe will be in a thermal, chemically equilibrated distribution of neutrino phase space (momentum and position) occupation probabilities: a Fermi-Dirac (FD) distribution. As the universe expands and cools and the energy-dependent weak interactions of neutrinos slow down relative to the expansion, the phase space occupation probability distributions will deviate from the FD form. Our objective is to follow these deviations with a procedure that can capture the important aspects of detailed Boltzmann transport calculations, but does so in an efficient manner with as few parameters as possible.  \nWe have found a surprising result: the essential physical properties and effects as derived from detailed Boltzmann transport can be captured with just three parameters per species. Refs. [39, 96] discuss a mathematical fitting procedure for the evolution of neutrino distribution functions. By contrast, our work here presents a statistical approach to this problem, one that is rooted in  \nquantum transport and that allows for ready physical interpretation of fitting parameters. We start with an exposition of this procedure directed at a quite general system, later specializing to the neutrino decoupling case.  \nIn describing this general procedure, it will prove useful to think first about how one uses the grand canonical ensemble in equilibrium statistical mechanics to derive the occupation numbers per unit volume of phase-space. After establishing that as a basis for comparison, we will generalize to non-equilibrium systems. There exist multiple derivations of the occupation numbers per unit volume of phase-space - called the distribution function n - in equilibrium statistical mechanics. For our purposes of exposition, the general n is a function of kinematic variables (position and momentum, i.e., x , p ) and quantum numbers (e.g., spin, flavor, lepton number, etc.).  \nOne approach to derive n employs the use of Lagrange multipliers with generalized operatorcharges assuming a form for the entropy distribution S in terms of n . For concreteness, we label the operator-charges as Q A with index A and associate to each Q A a Lagrange multiplier α A . In the example of a fermionic system in equilibrium subject to no external constraints, only two operator-charges are needed, namely, Q 0 = ∑ n for total particle number and Q 1 = ∑ En for total energy, where the summation is over all of the independent kinematic and quantum variables which describe n . The form of the entropy for this system is S = -∑ [ n ln n + (1 -n ) ln(1 -n )] derived from independent particle combinatorics. In this specific fermionic system without external constraints, the only kinematic variable is the energy E = √ | p | 2 + m 2 for particle mass m .  \nThe equilibrium condition is equivalent to the extremization of the entropy. By proceeding to extremize the entropy with the Lagrange multipliers on Q 0 and Q 1 , one can deduce the form of n as the standard Fermi-Dirac distribution  \nn { ℓ } ( E ) = 1 1 + e ( E -µ ) /T , (3.1)  \nwhere the subscript { ℓ } denotes the set of quantum numbers { ℓ } . Identifying the values of the Lagrange multipliers as  \nα 0 = -µ T , α 1 = 1 T , (3.2)  \nwhere µ and T are the chemical potential and temperature, respectively, makes the connection to classical thermodynamics and the grand potential.  \nIn writing Eq. (3.1), we have assumed no dependence of the quantum numbers { ℓ } on the distribution function n { ℓ } - the quantum numbers simply supply a label for the entropy summation. This implicit assumption that the quantum numbers only act to give the degeneracy in a given phase-space volume need not hold for a fermionic system under external constraints. Indeed, if the kinematics depend on quantum numbers, we would need to adopt different temperature and chemical potential variables for each set of quantum numbers. In fact, this is the case for the fermionic system of neutrinos in the early universe with an asymmetry in lepton number. Momentum-changing interactions between neutrinos tend to bring all species to the same temperature T . On the other hand, momentum-preserving-interactions would tend to bring all neutrinos to the same chemical potential, µ , while simultaneously bringing all anti-neutrinos to  \nthe same chemical potential opposite in sign to the neutrinos, -µ . In this case, the temperature is the same, but the chemical potential depends on the quantum lepton number.  \nAs emphasized in Sec. 2, the parameters describing the distribution function for each neutrino species evolve slowly with time until the rates of neutrino scattering begin to fall below the expansion rate of the universe and the neutrino component as a whole can no longer be described adequately by an equilibrium approximation. The out-of-equilibrium scattering of neutrinos even in the SM, and similarly in BSM physics, can lead to spectral distortions away from the thermal FD blackbody 2-parameter solution of Eq. (3.1). Capturing the effects of these spectral distortions will require more than the two parameters α 0 and α 1 .", '3.2 Entropy-content operator': "Consequently, we pursue an alternative procedure that we term the generalized entropy approach. We stress that this approach is general and can be used for other systems, although we will tailor it to capture moments of the early universe neutrino energy distributions. Inspired by a quantummechanical approach, we define a density operator ρ  \nIn writing Eq. (3.3), the density operator refers to an ensemble of multi-particle states. This ensemble may contain many different kinds of multi-particle states, i.e., fully-uncorrelated states, fully-entangled states, or a combination of the two. Each state has a well defined number of particles, N , and total energy, E , but not all states have the same N and E . As a result, we need to consider elements of a Fock space to describe the ensemble. In quantum mechanics, each N -particle state would have an associated wavefunction, which we denote as the ket | ψ ⟩ , and can be written using a basis of wavefunctions in a 3 N Hilbert space. The elements of the Fock space which we need to describe the ensemble are denoted as | k ⟩ and consist of the union of bases from each 3 N Hilbert space. We stay general and assume that | k ⟩ is not coincident with the eigenbasis of ̂ ρ . Hence, this is why we need the double summation in Eq. (3.3). Therefore, w k,k ' gives the density matrix element ⟨ k | ̂ ρ | k ' ⟩ . Very importantly: ̂ ρ contains off-diagonal elements. A member of the ensemble of ̂ ρ need not be limited to a single multi-particle wavefunction, or even a linear superposition of multi-particle wavefunctions in a single 3 N Hilbert space, but instead could be a superposition of states in the general Fock space. In this work, we consider the former possibility but neglect the later. In other words, ̂ ρ is a block-diagonal matrix in the Fock space, where each block corresponds to a 3 N Hilbert sub-space. We borrow the notion of superposition in a Hilbert space to write the quantity | Ψ ⟩ as a superposition of ensemble members in the Fock space. We will call this quantity a 'configuration' to discriminate it from the quantum-mechanical concept of wavefunction, although we utilize the bra-ket notation of quantum mechanics. For some configuration | Ψ ⟩ , we can write the weights as w k,k ' = ⟨ k | Ψ ⟩⟨ Ψ | k ' ⟩ and normalize them such that the trace of ρ is unity, i.e., Tr[ ρ ] ≡ ∑ k ⟨ k | ρ | k ⟩ = 1.  \n̂ ̂ ρ = ∑ k,k ' w k,k ' | k ⟩⟨ k ' | . (3.3)  \n̂ ̂ ̂ Now that we have oriented ourselves in the statistical Fock space, the next step of the generalized entropy approach is to promote the operator-charges Q A to operators ̂ Q A in the Fock space, and associate to each operator a control parameter α A , analogous to the Lagrange multipliers in the equilibrium treatment. To provide another example, ̂ Q N may be the number operator:  \na summation of tensor products of appropriate creation and annihilation operators ̂ a † ℓ ̂ a ℓ with identity operators  \n̸  \nwhich would give the number of particles in a many-body state | k ⟩ . In Eq. (3.4), we have used the shorthand notation ℓ = { ℓ } , x to denote a single particle state in the position basis (although we could use p instead to write the state in the momentum basis). A member of the Fock space which is also a member of a 3 N Hilbert space is an eigenvector of ̂ Q N , namely, ̂ Q N | ψ ⟩ = N| ψ ⟩ . A configuration in the Fock space is in general not an eigenvector. However, we can write  \n̂ Q N = ∑ ℓ   ̂ a † ℓ ̂ a ℓ ⊗ ∏ ℓ ' = ℓ ̂ I ℓ '   , (3.4)  \n⟨ Ψ | ̂ Q N | Ψ ⟩ ≡ ⟨ ̂ Q N ⟩ Ψ = N, (3.5) where we have defined the expectation value of ̂ Q N . In this case, for the configuration | Ψ ⟩ , we would expect N number of particles when selecting a multi-particle state from | Ψ ⟩ . In other words, N is the average number of particles in the configuration.  \nIn addition, note that the control parameters α A are not Lagrange multipliers. They will be functions of the various conditions of the environment under study, including time.  \nUsing the operators ̂ Q A , we define the entropy content operator which drives a system from an ( i )nitial configuration, e.g., | Ψ i ⟩ to a ( f )inal one ( | Ψ f ⟩ )  \nunder a specific realization of the control parameters { α ( fi ) A } . These control parameters are specific to the initial and final configurations, | Ψ i ⟩ and | Ψ f ⟩ , or equivalently ̂ ρ i and ̂ ρ f . F is a complex number and a functional of the control parameters. We scale it with the Fock space identity operator and explain its role below after we have introduced other key concepts.  \n̂ s ( { α ( fi ) A } ) = [ ∑ A α ( fi ) A ̂ Q A ] -F ( { α ( fi ) A } ) ̂ I, (3.6)  \nThe transformation under { α ( fi ) A } will take an initial operator ̂ ρ i to the final operator ̂ ρ f in the following manner  \n1 = ⟨ Ψ i | e -̂ s | Ψ i ⟩ ≡ ⟨ e -̂ s ⟩ i . (3.9)  \n̂ ρ f = exp[ -̂ s ( { α ( fi ) A } ) / 2] ̂ ρ i exp[ -̂ s ( { α ( fi ) A } ) / 2] . (3.7) We consider a set of transformations above such that the trace of the density matrix is preserved, i.e., Tr[ ̂ ρ f ] = Tr[ ̂ ρ i ] = 1 for our normalization. Note we use the same operator in the argument of the exponentials, implying the transformation is not necessarily unitary for an arbitrary choice of ̂ s . For the trace to be preserved, we must have the following condition Tr[ e -̂ s ̂ ρ i ] = 1 , (3.8) where we have used the cyclic property of the trace, the commutativity of identical operators, and abbreviated the exponential operator for brevity in notation. If we associate the configuration | Ψ i ⟩ to the density operator ̂ ρ i , we can show  \nBy Jensen's inequality  \n⟨ e -̂ s ⟩ i ≥ e -⟨ ̂ s ⟩ i , (3.10)  \nimplying a thermodynamic-second-law analog for the entropy-content operator  \n⟨ ̂ s ⟩ i ≥ 0 . (3.11) We return to the definition of ̂ s in Eq. (3.6) to discuss the functional F . Our goal is to use control parameters { α ( fi ) A } to follow the deformation of the density operator. However, the first term on the right-hand side of Eq. (3.6) does not preserve unitarity for general deformations. As a result, we introduce a normalization factor, i.e., F ̂ I , to ensure the unitarity of the transformation in Eq. (3.7). The identity operator commutes with any/all Fock-space operators ̂ Q A . Using Eq. (3.9), we can show  \nTo explain the significance of Eq. (3.13), we return to equilibrium statistical mechanics. One can write the Grand Canonical Partition function as the following  \ne F = 〈 exp [ -∑ α ( fi ) A ̂ Q A ]〉 -1 i , (3.12) = ⇒ F = -ln {〈 exp [ -∑ α ( fi ) A ̂ Q A ]〉 i } = -ln    ∑ k,k ' w k,k ' 〈 k ∣ ∣ ∣ ∣ exp [ -∑ α ( fi ) A ̂ Q A ] ∣ ∣ ∣ ∣ k ' 〉    . (3.13)  \nZ = ∑ i 〈 ψ i ∣ ∣ ∣ ∣ exp ( -˜ H i -µ ˜ N i T ) ∣ ∣ ∣ ∣ ψ i 〉 (3.14)  \nfor the ensemble of states indexed by i , with associated multi-particle wavefunctions | ψ i ⟩ , and quantum-mechanical Hamiltonian/Number operators ˜ H i / ˜ N i for the i th Hilbert space. Upon examination of Eq. (3.14) and the argument in the natural logarithm of (3.13), we see the connection between the generalized-entropy approach and the grand canonical ensemble. In equilibrium statistical mechanics, the grand potential, Φ G , is a generalization of the Helmholtz free energy from the canonical ensemble. Φ G is related to the grand canonical partition function by the following expression  \nΦ G = -T ln Z . (3.15)  \nIdentifying the argument of the natural logarithm in Eq. (3.13) as a non-equilibrium analog of Z from Eq. (3.14), we will call F the free-energy generating functional.  \nOne may be tempted to continue the analogy and extend the control parameters α ( fi ) A to the thermodynamic variables 1 /T and -µ/T if/once the system reaches an equilibrium state. This would be incorrect in our formalism. In the grand canonical ensemble, all states of the combined system and reservoir are equally likely, and physical processes act to exchange momentum, particle number, quantum numbers, etc., to change the system from one microstate to another. In our formalism, the control parameters act to drive the initial configuration into a final one. If the final density matrix is diagonal with coefficients which scale with the number of states the reservoir can obtain, then the system has reached a state of equilibrium. If ̂ ρ f is in equilibrium, the  \nsystem stays in this configuration in the absence of external influences 1 . As a result, the control parameters are all zero and trivially F = 0. For the rest of this subsection, we will drop the superscript label ( fi ) from the control parameters for brevity in notation, although stressing that the α A always take an initial configuration into a final one. Following the analogy of the GCE, we will take derivatives of Eq. (3.13) with respect to the control parameters to determine various thermodynamic/statistical analogs. However, the α A are functions of time, and also position perhaps. Therefore, F is a functional and we need a functional derivative. As no covariant derivatives, i.e., ∂ ν α ν A , appear in Eq. (3.13), the functional derivative mimics a partial derivative.  \nFunctional derivatives of the free-energy generating functional with respect to the α A parameters provide the connected correlation functions of the operators Q A  \n∣ where the subscript 'cc' denotes connected component, and the parentheses in the subscripts indicates symmetrization over the indices. The free-energy generating functional therefore plays the same role as the natural log of the partition function in quantum field theory and statistical mechanics, or the cumulant generating functional in statistics. In Appendix A, we provide justification for this statement and show explicit derivation of the low-order connected correlation functions from the free-energy generating functional.  \n̂ δ M F ( { α ( fi ) A } ) δα m 1 ...δα m M ∣ ∣ ∣ ∣ α A =0 = ( -1) M +1 〈 ̂ Q ( m 1 ... ̂ Q m M ) 〉 i, cc , (3.16)", '3.3 Connection to Wigner functions and the Mean field': "Our formalism for calculating the deformation of a systems' distribution function is completely general to this point. In principle, it could encompass the dynamics of systems in Fock space with inherent quantum coherence and entanglement. This generality may be necessary for following the energy distribution function evolution in the neutrino decoupling epoch if, for example, we include neutrino flavor transformation or various BSM scenarios involving macroscopic coherence. In a general approach the next step would be to re-configure the problem by using quantum correlations and translation symmetries to reduce two-point positions into an average position and momentum representation. This could be done via a generalized N -particle Wigner transformation. Although the equivalence of this N -particle approach to full quantum mechanics may seem obvious, it was controversial in its day, suggested by Weyl and by J. Groenewold [97]. Moyal and collaborators further developed this phase space formulation of quantum mechanics [98], and more recently C. Zachos, D. Fairlie, and T. Curtright [99] explored it in more detail.  \nHowever, for our example case of neutrino decoupling with no assumed quantum flavor coherence, our task is much simpler. We seek to find phase-space representations for the neutrino wavefunctions via a generalized N -particle Wigner transformation as alluded to above. Our configurations in Fock space are for multiple particles and anti-symmeterized to comply with FD statistics. As our goal at this point is only to find single-particle correlation functions, we will employ a mean field statistical approach to reduce the N -particle quasi-distribution to a single-particle correlation function. These objects will be necessary for evaluating the functional derivatives and operators in (3.16).  \n̸  \nWe give a sketch of how we construct the single-particle correlation function. First, we start with a generic multi-particle state | ψ ⟩ . If we write the multi-particle wavefunction as a tensor product of single-particle states, we can arbitrarily pick an entry in the tensor product which we label as the j th particle. We will take a generalized N -particle Wigner transform on the outer product | ψ ⟩⟨ ψ | . To calculate the matrix element of | ψ ⟩⟨ ψ | in the integrand of the Wigner transform, we write a tensor-product of position-quantum-number eigen-kets and bras, where the j th entry would have | X j + δ x j , s j , f j ⟩ , and ⟨ X j -δ x j , s ' j , f ' j | , for spins s and flavors f . If we take the Fourier conjugate of δ x j to be the momentum p j , we can calculate the generalized N -particle Wigner transform as a function of 6 N kinematic variables and 4 N quantum variables. We marginalize, i.e., integrate/sum over all of the other phase-space/quantum numbers, to obtain the single-particle Wigner function, which we denote as W s ' j f ' j ; s j f j ( X j , p j , t ). The Wigner transform and associated function yield a phase-space distribution from input single-particle quantum states. With the phase-space distribution and the free-energy generating functional F , we can proceed to derive a form for the distribution in terms of the α A control parameters.", '3.4 Free-energy generating Functional for Fermions': "If we start in an initial configuration i we will have a set of occupation numbers and corresponding constrained expectation values ⟨ ̂ Q A ⟩ of the generalized charges. The entropy-content operator that will generate the change in the distribution functions in going from the initial configuration to a final 'deformed'configuration f is  \n̂ s fi | α ≡ ̂ s ( q i ; { α ( fi ) A } ) = ∑ A α ( fi ) A ̂ Q A ( q i ) -F ( { α ( fi ) A } ) ̂ I , (3.17)  \nwhere q i denote the set of parameters characterizing the initial configuration of the system (e.g., a set of single particle occupation probabilities for each neutrino type) and { α ( fi ) A } are a corresponding set of deformation parameters. Hereafter, we again drop the superscript ( fi ) on the α parameters for brevity. For a given realization of a set of α parameters, { α } , the mean of the entropy content operator (hereafter the 'mean entropy') is  \nand the entropy fluctuation from the mean is the difference  \nδ ̂ s fi ( q ) | α = ̂ s fi ( q ) | α -⟨ ̂ s fi ⟩| α ̂ I = ∑ α A ( ̂ Q A ( q ) -⟨ ̂ Q A ⟩| α ̂ I ) , (3.19)  \n⟨ ̂ s fi ⟩| α = ∑ α A ⟨ ̂ Q A ⟩| α -F ( { α A } ) , (3.18)  \nwhere again ⟨ ̂ Q A ⟩| α denotes an expectation value for ̂ Q A for a given set { α } . The (non-equilibrium) evolution equation for the density operator ̂ ρ i is the familiar N -point quantum Liouville equation, or master equation when written in terms of matrix elements. Here, for our problem of neutrino energy distribution evolution in the early universe, we will make the usual approximation of mean-field evolution of neutrinos so that neutrinos can be treated as quasiparticles. For example, this is the usual approach in neutrino flavor transformation in the early universe and other astrophysical environments [100-102].  \nThe initial density matrix, ̂ ρ i , is in an equilibrium configuration, where the many-body wave functions are appropriately weighted in the GCE (i.e., both thermal and chemical equilibrium),  \nour starting point for the transport and evolution of the neutrino distribution functions. Full transport here is captured with a time-sequence of non-equilibrium density operators. We have used the transformation in Eq. (3.7) and the symbol ̂ ρ f for the deformed density operator. At this point, we will denote the time-sequence of deformed density operators as ̂ ρ ( t ), with the initial condition ̂ ρ ( t = t i ) = ̂ ρ i . The code burst has been used to evolve the neutrino distribution functions using appropriate collision kernels [44]. Here, our intent is to capture that result in terms of the evolution of an appropriate expansion of mathematical objects derived from the density operator ̂ ρ ( t ), namely Eq. (3.16). We cast the density operator in terms of a singleparticle representation. This representaion follows from an underlying independent quasiparticle state basis, where each pure state in this basis has been Wigner-transformed in terms of a phase space, i.e., X , p . We can summarize this picture with the following form of the density operator  \n̂ ρ ( t ) = e F ⊗ Xp s ' f ' sf e -α ̂ n , (3.20)  \nwhere the spin and flavor labels of the control parameters α s ' f ' ; sf ( t ; X , p ) and the phase-space number operator ̂ n s ' f ' ; sf ( X , p ) as well as their dependence on time, X and p are implicit. In writing Eq. (3.20), we have put in two conventions. Firstly, we have placed time dependence in the control parameters. Recall that as defined, the control parameters are used for deformations of the density operator, and are zero when ̂ ρ f remains in equilibrium. In Eq. (3.20), the initial α ( t = t i ) are used to specify the initial density operator (e.g., a FD equilibrium configuration) and α ( t > t i ) give deviations from that initial configuration. Secondly, we no longer have the subscript A used to associate a control parameter with a Fock space operator ̂ Q A . Now, we have used the quantum numbers and phase-space variables to specify a control parameter for each single-particle quantum state.  \nUsing the unitarity of the trace, we can derive an expression for the mean value of ⟨ n s ' f ' ; sf ( X , p ) ⟩ .  \n̂ 1 = Tr[ ̂ ρ ( t )] = Tr   e F ⊗ Xp s ' f ' sf e -α ̂ n   = e F Tr   ⊗ Xp s ' f ' sf e -α ̂ n   = e F ∏ Xp s ' f ' sf Tr[ e -α ̂ n ] (3.21)  \nFor particles obeying Fermi-Dirac statistics, the number operator for each cell in phase space has two distinct eigenvalues: 0 and 1. Thus, the trace over the exponentiated number operator is  \nTr[ e -α ̂ n ] = 1 + e -α , (3.22)  \nwhich implies the trace expression becomes  \n1 = e F ∏ Xp s ' f ' sf [1 + e -α ̂ n ] (3.23)  \nSolving for F gives  \nF ( t ) = -∑ Xp s ' f ' sf [ln(1 + e -α )] , (3.24)  \nwhich we can substitute back into Eq. (3.20) to yield  \n̂ ρ ( t ) = exp   -∑ Xp s ' f ' sf ln(1 + e -α )   ⊗ Xp s ' f ' sf e -α ̂ n = ∏ Xp s ' f ' sf exp [ -ln(1 + e -α ) ] ⊗ Xp s ' f ' sf e -α ̂ n = ∏ Xp s ' f ' sf (1 + e -α ) -1 ⊗ Xp s ' f ' sf e -α ̂ n = ⊗ Xp s ' f ' sf (1 + e -α ) -1 e -α ̂ n . (3.25)  \nWe can find the expectation value for a given phase-space cell ( X 1 p 1 s ' 1 f ' 1 s 1 f 1 ) by the following  \n̸  \n⟨ ̂ n X 1 p 1 s ' 1 f ' 1 s 1 f 1 ⟩ = Tr[ ̂ n ̂ ρ ] = Tr   ( ̂ ne -α ̂ n 1 + e -α ) X 1 p 1 s ' 1 f ' 1 s 1 f 1 ⊗ ⊗ Xp s ' f ' sf = X 1 p 1 s ' 1 f ' 1 s 1 f 1 e -α ̂ n 1 + e -α   = ( e -α 1 + e -α ) X 1 p 1 s ' 1 f ' 1 s 1 f 1 = ( 1 1 + e α ) X 1 p 1 s ' 1 f ' 1 s 1 f 1 . (3.26)  \nThis result can be rearranged to give α X 1 p 1 s ' 1 f ' 1 s 1 f 1 in terms of ⟨ n X 1 p 1 s ' 1 f ' 1 s 1 f 1 ⟩ as  \nThe density matrix then can be written as  \nρ ( t ) = e F ⊗ Xp s ' f ' sf exp([ ̂ n ln ⟨ ̂ n ⟩ -̂ n ln(1 -⟨ ̂ n ⟩ )]) , (3.28)  \n̂ α X 1 p 1 s ' 1 f ' 1 s 1 f 1 = [ln(1 -⟨ ̂ n ⟩ ) -ln ⟨ ̂ n ⟩ ] X 1 p 1 s ' 1 f ' 1 s 1 f 1 . (3.27)  \nand the free energy can then be expressed as  \nsuch that the density matrix is  \nF ( t ) = ∑ Xp s ' f ' sf [ln(1 -⟨ ̂ n ⟩ )] , (3.29)  \nρ ( t ) = ⊗ Xp s ' f ' sf exp[ ̂ n ln ⟨ ̂ n ⟩ +(1 -̂ n ) ln(1 -⟨ ̂ n ⟩ )] . (3.30)  \nThe mean Wigner matrix can be obtained as the derivative of the free energy with respect to α  \n⟨ ̂ n ⟩| α = ∂ F ∂α = ( e α +1) -1 , (3.31)  \nand the mean entropy is given by the usual expression  \n⟨ ̂ s fi ⟩| α ( t ) = ∑ Xp s ' f ' sf α ⟨ ̂ n ⟩| α -F = -∑ Xp s ' f ' sf [ ⟨ ̂ n ⟩ ln ⟨ ̂ n ⟩ +(1 -⟨ ̂ n ⟩ ) ln(1 -⟨ ̂ n ⟩ )] . (3.32)  \nWithin the independent particle approximation, we can express the full density matrix in terms of a set of density matrices deformed by α s ' f ' ; sf ( t ; X , p ). That is, the transport is fully encoded by these control parameters, but there are as many control parameters as there are dimensions of phase space. This change of perspective therefore retains all of the information, but it also retains all of the complexity of describing the evolution of each of the individual particles in the system.  \nIn many situations, however, it is possible to describe the relevant features of the system with a much smaller set of parameters. When such a compression is possible, we can express the evolution of the system in terms of a smaller set of deformation parameters { α A } . Therefore, at a given time t we can explore ρ ( t ; q, q ' ) = ρ ( q, q ' ) | α A ( t ) via its deformations by imposing various applied α A , where the α A would be a truncated set relative to the full α phase-space control-matrix. This set of control parameters α A specifies the constraints to which the system is subjected at any fixed time. The interesting question is how close we can get to an accurate description of the full evolution of the system ρ ( t + δt ) /ρ ( t ) by using a truncated set of evolving control parameters α A ( t + δt ). If the number of variables A = 0 , 1 , . . . is small then we have reduced the transport problem considerably. Of course this is familiar if we have the densities described by a sequence of equilibrium states, characterized at each time by a chemical potential α 0 = -βµ and an inverse temperature α 1 = βE . This equilibrium case provides an example of how the salient features of the evolution can be captured with just a few parameters even when the system contains many degrees of freedom. We will now apply this formalism to the case of neutrino decoupling.  \nThe energy dependence of the neutrino cross sections implies that low energy neutrinos are not as tightly coupled as high energy ones. This might be thought to motivate a simple nearequilibrium distribution at high energy with a modification at low energy. Instead we have found that the expansion of α Xp s ' f ' sf ( t ) in powers of momenta to one order higher than the equilibrium case captures the essence of the Standard Model C ν B. Allowing for control parameters which describe the number density current and energy-momentum current, along with one additional control parameter related to the current of the energy-momentum variance, the expansion takes the form  \nα s ' f ' ; sf ( Xp ) = α (0) µs ' f ' ; sf p µ E + α (1) µν s ' f ' ; sf p µ p ν T cm E + α (2) µνλs ' f ' ; sf p µ p ν p λ T 2 cm E + . . . . (3.33)  \nThis expansion remains quite general and could, in principle, be applied to a wide array of problems, including neutrino transport in neutron stars and supernovae (in which case the factors of T cm appearing in this expression would be replaced by a characteristic energy scale relevant to the problem). However, for describing the cosmic neutrino background, the homogeneous, isotropic, and hot conditions of the early universe lead to further simplifications of this expression. Cosmic neutrinos are ultra-relativistic in the early universe such that E ≃ | p | , and we further focus on the case diagonal in spin and flavor, leaving the problem of describing neutrino flavor  \noscillations within this framework to future work. As a result, this expression simplifies to the form  \nα = α (0) + α (1) ϵ + α (2) ϵ 2 + · · · , (3.34)  \nfor each neutrino species. A new ingredient is α (2) , which is, in the Lagrange multiplier framework, indicating that fluctuations in the total energy ⟨ (∆ U ) 2 ⟩ are constrained as well as the number and energy densities. As our figures show, a description of the system using just α (0) and α (1) is in fact quite accurate over much of the neutrino decoupling epoch, but adding α (2) considerably sharpens the agreement at low T cm .  \nHow do we measure the accuracy of such a radical compression of full transport into just three variables per neutrino species? A very general approach is based on the the relative entropy (also called the negative of the Kullback-Liebler divergence) between the two distributions, which for fermionic systems takes the form  \nD ( a ∥ b ) = -∑ Xp s ' f ' sf [ ⟨ ̂ n ⟩| α a ln ⟨ ̂ n ⟩| α a ⟨ ̂ n ⟩| α b +(1 -⟨ ̂ n ⟩| α a ) ln (1 -⟨ ̂ n ⟩| α a ) (1 -⟨ ̂ n ⟩| α b ) ] . (3.35)  \nThis treatment can be used to assess the similarity of the deformed state that results from the full transport calculation α a (in our case provided by the burst code) to the state that is approximated by our three-parameter fit α b . This prescription works even if the distributions being compared are far from equilibrium. In practice, for the near-equilibrium evolution of the Standard Model thermal history, we use a goodness-of-fit parameter that measures the weighted mean square deviation of the distribution at each sampled point throughout the decoupling epoch; see Sec. 7.  \nBefore applying this formalism to the C ν B, we will briefly discuss an application of the same formalism to the (perhaps more familiar) case of CMB transport.", '4 Relation to CMB transport': "The closest analogue of the calculations we are doing here is to the development of spectral distortions of the CMB [53, 103-107], and it is instructive to compare CMB transport with the C ν B case. This section will summarize the treatment of CMB transport described in Ref. [108].  \nFor the decoupling of photons from the plasma in the early universe, there are three identifiable epochs:  \n- 1. The blackbody epoch at z > z Pl ≈ 10 6 . 9 ( Ω b h 2 0 . 01 ) -0 . 39 , above which Bremsstrahlung ep → ep + γ and Double Compton scattering γe → γe + γ are sufficiently fast that a Planck equilibrium distribution forms. The redshift z Pl therefore defines the cosmic photosphere. The photon distribution in this epoch is characterized by just α (1) = βE .\n- 2. The chemical potential µ γ epoch, at z Pl > z > z BE ≈ 10 5 . 6 ( Ω b h 2 0 . 01 ) -1 / 2 when Compton scattering dominates, photon number is conserved with injected energy redistributed fast enough to give a Bose-Einstein shape, and α is well characterized by α = α (0) + βE , with α (0) = -βµ γ .  \n- 3. The Compton scattering-dominated epoch of thermal Sunyaev-Zeldovich y -distortions, at z < z y ≈ 10 5 ( Ω b h 2 0 . 01 ) -1 / 2 , when the energy redistribution is slow, resulting in the characteristic perturbative thermal Sunyaev-Zeldovich y -distortion shape: ∆ α = -yxψ K ( x ) with y ≡ τ ( T e -T ) /m e , τ is the optical depth, x ≡ ω/T , and ψ K ( x ) ≡ 2 -x 2 e x +1 e x -1 . Below z y , energy redistribution by Compton scattering is achieved through small positive energy kicks, ⟨ δE/E ⟩ ∝ T e /m e , and root mean square diffusion in energy ⟨ ( δE/E ) 2 ⟩ ∝ T e /m e , encoded in the Kompaneets diffusion equation for the photon distribution function.  \nThere are intermediate transport regimes between the three epochs; for example, in between z Pl and z BE some photons lock into a Planck distribution and others do not [103, 104, 108, 109].  \nTo understand how these epochs arise, and to estimate magnitudes of the three characteristic redshifts, an analytic treatment based on Ref. [103] is quite adequate. The transport equation for the photon distribution function can be written as  \nd ⟨ ̂ n ⟩ dt = C K + C bremss + C DC , (4.1)  \nwhere the right hand side gives the sum of the Kompaneets, bremsstrahlung, and double Compton source terms. We can write the photon distribution function as  \n⟨ ̂ n ⟩ = ( e α -1) -1 . (4.2) We linearize the transport equation in ∆ α = α -ϵ . In the tight coupling regime, C K + C bremss + C DC approximately vanishes; this condition is satisfied for small ϵ e = E γ /T e , where T e is the electron temperature, if α = α 0 ( t ) exp( -ϵ e /ϵ 0 ). Thus for low frequencies, ϵ e < ϵ 0 , bremsstrahlung and the double Compton process dump photons in fast enough to yield a Planck form, but for ϵ e > ϵ 0 the Bose-Einstein form prevails. Here ϵ 0 = (4 ϵ 3 (Γ bremss +Γ DC ) / Γ K ) 1 / 2 , where the 'Kompaneets' rate is Γ K ≡ 4 n e σ T T e /m e , in terms of electron density and Thomson scattering cross section. The approximate constancy of ϵ 3 (Γ bremss +Γ DC ) is exploited to obtain this result. If we assume an injection rate of ˙ Y γ photons per baryon with average energy ¯ E γ at time t , adding to the Y γ 0 photons per baryon already there, then the scaling parameter α 0 evolves according to  \ndα 0 dt = -α 0 τ D + ( ¯ E γ 3 . 6 T γ -1 ) 1 . 87 ˙ Y γ Y γ 0 , τ D = 1 . 29 x 0 [(Γ bremss +Γ DC ) x 3 ] -1 = 1 . 29(Γ K / 4) -1 x -1 0 . (4.3)  \nThus there is a damping term with timescale τ D which pushes α 0 towards zero, i.e. toward a Planck distribution. The injection term works in the opposite direction, driving the distortion. When the damping time is shorter than the expansion rate of the universe, any injected energy input would be re-thermalized into a Planckian distribution in equilibrium with the electrons within one Hubble time. The transition between these epochs defines z Pl . When the Kompaneets rate is a few times the expansion rate, x 0 will be low but α 0 will not be zero, and the Bose-Einstein form is appropriate. This defines z BE . However, it is not until the Kompaneets rate is a few times below the expansion rate that the perturbative y -distortion solution prevails. This defines z y ≈ z BE / 4. Naturally z BE and z y scale in the way defined by Γ K /H .  \nWe can see from this example that following the evolution of α provides a convenient method for calculating spectral distortions during CMB decoupling. We would like to see if there is any relatively simple form for α for the C ν B as there is for the CMB. A key ingredient in the CMB case is that the photons that are shaken off in ep → ep and γe → γe are soft photons which must be redistributed over all energies by Compton scattering. For neutrino production, the energies are not soft, with highly energy-dependent reaction rates ∝ E . Indeed one of the main highlights in this paper is identifying the extreme energy-dependence of the C ν B neutrino-sphere above which equilibrium holds. By contrast, the CMB photosphere is essentially energy-independent.  \nThe inhomogeneous transport of photons through CMB photon decoupling has a largely fixed spectral shape, and has usually been done in ∆ T s ' s variables (see e.g. [108]), but for accuracy we now recommend using α s ' s ( X , p , t ) or, ignoring the spectral distortions given the dominance of the unperturbed Planck spectrum, the comoving inverse-temperature and its fluctuations a -1 β γ ; s ' s ( X , p , t ), since it does not have the nonlinear corrections that ∆ T has.  \nThe inhomogeneous C ν B transport through the time of CMB photon decoupling and to later times as the effects of neutrino masses become important is also best done with the comoving inverse temperature, a -1 β ν ; s ' f ' sf ( X , p , t ). Unlike for photons, the dependence of the distribution on the magnitude of the comoving momenta, ϵ ≡ p/T cm , is important because neutrinos of different momenta experience the effects of neutrino masses at different cosmological epochs. While the formalism presented here allows for treatment of the full inhomogeneous neutrino transport, without a means to directly detect the C ν B, cosmological observations are only sensitive to the homogeneous evolution. Perturbations in the C ν B density do, however, leave observational imprints on the CMB and large scale structure power spectra since they freely stream after C ν B decoupling [110-112]. The shift in phase of the acoustic peaks due to the effects of neutrino free streaming has been observed in the CMB [113] and in baryon acoustic oscillation data [114].", '5 The Collision Rate Operator': "In the previous sections we described the general transport framework that we will utilize and described how the framework can be applied to the (perhaps more familiar) case of CMB transport. We now turn to a description of the interactions relevant to neutrino transport. As mentioned above, we will neglect neutrino flavor oscillations in the treatment presented here.  \nThe Boltzmann equation for neutrino transport can often be thought of as consisting of a coarse-grain change in n as a result of quasiparticle free-streaming evolution plus a fine-grain flow of information to the coarse-grain from the net result of collisions: D ⟨ ̂ n ν f ⟩| α /dt = C ν f [ n ν ], where C ν f [ n ν ] dt encodes the net change in the mean Wigner function due to short scale collisions. In the absence of collisions, ⟨ ̂ n ν ⟩| α is conserved along geodesic trajectories. If we force α to be tightly-coupled to n through α = ln( ⟨ ̂ n ν ⟩| -1 α -1), then Dα/dt = ( d ⟨ ̂ n ν ⟩| α /dα ) -1 C = ⟨ δ ̂ nδ ̂ n ⟩ -1 C . Here ⟨ δ ̂ nδ ̂ n ⟩ = ⟨ ̂ n ν ⟩| α (1 -⟨ ̂ n ν ⟩| α ). Generally the collision operator evaluated in the independent particle approximation is full of products of ̂ n up to the fourth power, and even if α is transported, the collision rate operator will still use this product expansion in ̂ n . The collision term is derivable directly from the density matrix formalism as long as we go beyond the independent particle approximation to include the interaction Hamiltonian in ρ i and ρ f , and expand to second order in the Fermi 4-point coupling G F . This is equivalent to the usual derivation of C where  \nfactors of ̂ n are assumed to be correlated only through the independent particle approximation in which correlators of higher order can be expressed as polynomial products of the first order correlators ⟨ ̂ n ⟩| α , ⟨ ∏ ̂ n f ⟩| α,cc = ( -1) N ∂ N F / ( ∂α 1 · · · ∂α N ). Given that, there is a clean separation of Fermi Golden rule squared matrix elements and products of the statistical factors. For 2-body interactions, the only factors entering are ⟨ n ⟩ and Fermi-blockers 1 -⟨ n ⟩ .  \nT cm = T in a in [ 1 a ( t ) ] , (5.1)  \n̂ ̂ Neutrinos ( ν f ) and anti-neutrinos (¯ ν f ) of flavors f = e, µ, τ collide with each other and with electrons ( e -) and positrons ( e + ) through various two-body interaction channels, the latter resulting in entropy flow between the e -/ e + /photon plasma and the neutrino sea, resulting in distortions of the mean neutrino occupation probabilities. Here and below, where the context avoids potential confusion, we will use the notation n ν f to refer to the expectation value of the occupation probability ⟨ ̂ n ⟩ f . The n ν f distributions are time-dependent functions of coarse-grained positions X and relative 3-momenta p , and time, t . We assume that the neutrino sea is spatially homogeneous, i.e. , X -independent, and independent of the 3-momentum direction ̂ p , so only a function of its magnitude p and time t . The coarse-grained evolution is Dn/dt = [ ∂ ∂t -H ( a ) p ∂ ∂p ] n , if we assume the neutrino masses are small compared with the momenta of interest. In that case it is advantageous to transform from p ≈ E p to a comoving energy parameter ϵ = E ν T cm , where T cm ∝ a -1 is a measure of the scale factor expressed in comoving temperature units:  \nwhere T in and a in are the initial values of temperature and scale factor of the plasma before the weak decoupling epoch.  \nIn the flavor eigenbasis, in the Boltzmann equation for neutrinos of flavor f as a function of time and comoving energy ϵ ,  \nd dt n ν f | ϵ = C ν f [ n f ' ] = ∑ r C ( r ) ν f [ n f ' ] , (5.2)  \nthe collision rate operator is decomposed into interaction channels r in which flavor f ' stimulates changes in flavor f .  \nEquation (5.2) gives the equations of motion for the neutrino distribution functions indexed by ϵ . Section 7 presents results from the code burst which integrate Eq. (5.2) forward in time along with fits using the dual-description α from Sec. 3. The description of n ν f | ϵ is indeed dual to α ν f as can be seen by transforming Eq. (5.2) into the new dependent variable  \nd dt α ν f = C ν f [ n f ' ] e -α ν f ( e α ν f +1) 2 , (5.3)  \nwhere the collision term would need to be rewritten by substituting 1 / [ e α f ' ( ϵ ) + 1] for n f ' | ϵ . Although we show fits using α ν f later in Sec. 7, one could in principle discretize Eq. (5.3) and integrate the set of ODEs forward in time. We leave that program to future work (see Appendix B for details).  \nReturning to the discussion of the collision term using occupation numbers n ν f | ϵ , the general form for a 2-body interaction 1 + 2 ↔ 3 + 4 is given by Fermi's Golden Rule,  \nC ( r ) ν 1 [ n f ] = 1 2 E 1 ∫ d 3 p 2 (2 π ) 3 2 E 2 d 3 p 3 (2 π ) 3 2 E 3 d 3 p 4 (2 π ) 3 2 E 4 (2 π ) 4 δ (4) ( p 1 + p 2 -p 3 -p 4 ) × S r ⟨|M r | 2 ⟩ F r ( p 1 , p 2 , p 3 , p 4 ) . (5.4)  \n⟨|M r | 2 ⟩ is the squared amplitude of the weak interaction r , averaged over initial spin states and summed over final spin states. S r is the symmetrization factor for identical particles, where S r = 1 2 for identical particle interactions, and S r = 1 for all other interactions. δ (4) ( p 1 + p 2 -p 3 -p 4 ) is the energy and momentum conserving delta-function. The function F r ( p 1 , p 2 , p 3 , p 4 ) includes the independent particle Fermi-blocking factors 1 -⟨ ̂ n ⟩| α in the outgoing channels as well as the occupation number product in the incoming channels:  \nF r ( p 1 , p 2 , p 3 , p 4 ) = (1 -n 1 )(1 -n 2 ) n 3 n 4 -n 1 n 2 (1 -n 3 )(1 -n 4 ) = F (in) r -F (out) r . (5.5)  \nHere F (in) r is the statistical factor associated with the probability of scattering into a state, and F (out) r is the statistical factor associated with the probability of scattering out of a state. It is assumed that electrons and positrons maintain a thermal Fermi-Dirac distribution during this time, with the electrons at the plasma temperature T e ≈ T γ . Thus we need only solve the Boltzmann equations for neutrinos and anti-neutrinos. Thus p 1 always refers to either a neutrino or anti-neutrino, while p 2 , p 3 , p 4 refer to a neutrino, anti-neutrino, or charged lepton. Table (1) outlines all the weak interaction processes, along with their scattering matrix term G -2 F S r ⟨|M r | 2 ⟩ , where G F = 1 . 166 × 10 -11 MeV -2 is the Fermi 4-point interaction constant. For completeness we list the terms r as derived for use in the burst code, as given by an analogous table in [44].  \nThere is an analogous set of weak interactions for anti-neutrinos, except for processes r = 6 , . . . , 9, where the squared matrix term is the parity-conjugate of the respective process in Table 1. We are operating at temperatures far below the muon and τ mass, yet far above the neutrino masses, so in terms of reaction rates there is no difference between the behavior of ν µ and ν τ . 2 Of course we are just in the regime in which the temperature is comparable to the electron mass, hence the amplitudes of ν e interactions includes charged current interactions, modifying the (2 sin 2 θ W -1) to a (2 sin 2 θ W +1) factor, where θ W is the weak mixing angle (or Weinberg angle), with sin 2 θ W ≈ 0 . 23.  \nThe computation of neutrino transport is iterative: collision rates are calculated for all neutrino flavor distribution functions at a given time step t which kick the distributions into a new state at t + dt . Those newly calculated distributions are then used to calculate the collision rates at the next time step, and so on. The 9-dimensional integral in Eq. (5.4) can be reduced to a 2-dimensional integral over ϵ variables. The reduction and subsequent formulation of the rates is spelled out in detail in Appendix B and C of Ref. [44].  \nTable 1 . This table displays the weak interaction processes for neutrinos. r = 5 , 10 , 11 represent annihilation channels, while the other processes are scattering channels. P i represents the momenta of neutrinos or anti-neutrinos, while Q i represents momenta of charged leptons. Anti-neutrinos undergo the same set of interactions, except the matrix amplitudes are the parity-conjugate for processes r = 6 , . . . , 9.  \n-  \n-  \n·", '6 The Neutrino BURST Numerical Computations': 'One application of the maximum-entropy distributions discussed above is to evaluate the response of the light-element abundances to the cosmic neutrino distributions. To accomplish this task, we need to run numerical experiments on BBN using a computation with sufficient energy and time resolution. We use the code burst [44] for these experiments. burst calculates neutrino-energy transport without resorting to any approximations of equilibrium. Using a suitable numerical implementation of the collision integrals in Sec. 5, we can time-evolve the neutrino distributions through the weak decoupling epoch and couple those distributions to the neutron-to-proton interconversion rates. The free nucleon abundances are then input into the nuclear reaction network to calculate the light-element abundances. The progression of neutrino scattering, isospin-changing reactions, nucleosynthesis is all within the context of the standard cosmology. By extending the standard-model burst computation, we have the ability to conduct numerical experiments with perturbed neutrino distributions.  \nWe refer the reader to Ref. [44] for details of the numerical implementations in burst . We give a few pertinent details here for the specific experiments in this work. We use 6 distributions (3 flavors of neutrinos and anti-neutrinos) with 101 energy bins, linearly spaced from ϵ = 0 to 25. All 101 energy bins are used when calculating energy densities and the weak-isospin-changing rates. For the neutrino scattering rates, we do not evolve the distribution for ϵ = 0. We do not  \nset the cosmological parameter N eff as an input for our nucleosynthesis calculations. N eff is an output from the neutrino-energy-transport and other electromagnetic processes. Finally, we only consider neutrino kinetics for the energy scales of BBN - scales much larger than the neutrino rest masses. We always treat neutrinos as ultra-relativistic and ignore any rest-mass contribution to the total energy density of the universe at the time of BBN. Our results are agnostic in regard to the cosmological parameter ∑ m ν .  \nWe neglect the role of neutrino flavor oscillations in this work. The computational problem of including oscillations with collisions has been done within the standard cosmology with a mean field, quantum kinetic equation (QKE) approach [21-23] and in a non-zero lepton-number extension [66]. For the case of the standard cosmology, the occupation numbers for the neutrino energy distributions have the same qualitative features between the Boltzmann and QKE calculations throughout the weak decoupling epoch. Importantly, the electron neutrino (and anti-neutrino) distributions have smaller deviations from FD equilibrium for the QKE case compared to the Boltzmann case. This result derives from the fact that with oscillations present, an electron neutrino has a non-zero probability to oscillate to a µ or τ flavor neutrino. As an example, our results for the linear response of the nuclide abundances to spectral distortions using Boltzmann transport could be applied to the distributions with oscillations included if we were to rescale those responses to smaller values. In theory, neutrinos could follow a much more complicated evolution during and after weak decoupling if collisions were to damp any coherence in the neutrino wave functions at early times and the energy-dependent vacuum Hamiltonian introduced a non-zero flavor coherence at later times. However, the results of Refs. [21-23] indicate the neutrino density matrices and the vacuum Hamiltonian are coincident with one another, i.e., that neutrino wave functions are parallel to the mass eigenvectors instead of the flavor eigenvectors and are therefore immutable from the vacuum potential at late times. This finding holds over the entirety of weak decoupling implying that the distribution functions follow the same trajectories in both the Boltzmann and QKE treatments, but with different amplitudes. As a result, we expect that our results here are transferable to the more general QKE problem in the standard cosmology, but the abundance-responses would need to be scaled down. For a BSM model, the results could be qualitatively different between Boltzmann and QKE treatments. Moreover, even in the Standard Model case, the mean field QKE approach may not be adequate to capture the neutrino flavor evolution. Quantum many body effects might be important [115], though the status of studies of this issue are mixed [116]. Some show the development of entanglement entropy with time and with the number of neutrinos, a marker for the break down of the mean field approach. However, these studies [117-119] focus on lower entropy, compact object-like conditions and it is not clear how that will scale to the early universe. In any case, a potential target for our analysis will be to look for signatures of BSM or many-body effects in altering the light element and N eff observables.  \nThe code begins evolving cosmological variables at a temperature T = 30MeV, where all of the neutrino species are in thermal equilibrium with the electrons and positrons. During this time, the comoving temperature is equal to the plasma temperature. Neutrino transport begins at a temperature of T = 10MeV, after which the comoving temperature begins to deviate from the plasma temperature. Neutrino transport terminates at a temperature of T cm = 0 . 015 MeV, which is after the annihilation of the thermal electrons and positrons. burst uses an adaptive  \n〉  \n̂  \n〈  \n〉  \n̂  \n〈  \nFigure 3 . The mean relativistic Fermi-Dirac occupation probability ⟨ ̂ n ( ϵ ) ⟩ at temperature T = T cm in terms of ϵ = p/T cm and the normalized density fraction ϵ 2 ⟨ ̂ n ( ϵ ) ⟩ / ( ϵ 2 ⟨ ̂ n ( ϵ ) ⟩ ) max of a species described by this distribution. The equilibrium distribution function accurately describes cosmic neutrinos at high temperatures, before weak decoupling. The normalized density fraction is used below to weight parametric fits of the non-equilibrium distribution function.  \n<!-- image -->  \ntime step based on convergence criteria from the previous time step. The code takes on order 10 4 time steps to progress through weak decoupling and BBN.', '7 The Evolution of Neutrino Distribution Functions and their Conjugates in Weak Decoupling': 'At high temperatures, cosmic neutrinos of each flavor f were in thermal equilibrium with the plasma and were well described by a relativistic Fermi-Dirac distribution  \n⟨ ̂ n ν f ⟩ = 1 e ϵ +1 , (7.1)  \nwhich is shown in Figure 3. As neutrinos decouple from the plasma, their distribution begins to deviate from the equilibrium Fermi-Dirac distribution. We will characterize the evolution of the neutrino distributions by  \n⟨ ̂ n ν f ⟩| α = 1 e α ( f ) ( ϵ,T cm ) +1 , (7.2)  \nwhere the quantity α ( ϵ, T cm ), introduced in Sec. 3, describes the out-of-equilibrium neutrino transport. In this notation, the equilibrium occupation probability is given by α ( ϵ, T cm ) = α eq = ϵ .  \nAs described in Sec. 3, the definition of the occupation probability in terms of α is motivated by its connection to the entropy. The entropy per unit phase space is given by  \n⟨ ̂ s f ⟩ = -⟨ ̂ n ν f ⟩ ln ⟨ ̂ n ν f ⟩ -(1 -⟨ ̂ n ν f ⟩ ) ln(1 -⟨ ̂ n ν f ⟩ ) . (7.3)  \nWith the occupation probability given by Eq. (7.2), the derivative of the entropy with respect to the distribution function is given by  \nd ⟨ ̂ s f ⟩ d ⟨ n ν f ⟩ = α ( f ) . (7.4)  \n̂ This shows the direct link between the response in entropy to the stimulus of distribution function change through α , which acts here as a stimulus-response transport coefficient.  \nAs described below, we can capture the evolution of the non-equilibrium distribution function through the epoch of weak decoupling in the Standard Model very precisely by modeling the α ( f ) ( ϵ, T cm ) as polynomials in ϵ . Furthermore, deviations from the equilibrium distribution are small in the Standard Model, and can be treated perturbatively in this case.', '7.1 Linear Model for α Evolution': 'First, we model the departure from equilibrium with a simple two-parameter model for the neutrino occupation probability, defined as  \nα ( f ) fit = α ( f ) 0 ( T cm ) + α ( f ) 1 ( T cm ) · ϵ . (7.5)  \nIn this model, the parameter α 0 ( T cm ) plays a role similar to the chemical potential divided by the temperature, but here the chemical potentials of the neutrinos and antineutrinos have the same sign. A positive value of α 0 ( T cm ) represents a deficit of both neutrinos and antineutrinos compared to an equilibrium distribution with vanishing chemical potential. Meanwhile the parameter α 1 ( T cm ) is related to the mean energy, and plays a role similar to the effective temperature of the neutrinos.  \nIn order to determine the values of the parameters α ( f ) 0 ( T cm ) and α ( f ) 1 ( T cm ), we performed a density-weighted fit of the model shown in Eq. (7.5) to the output of the burst code for each value of T cm and for each neutrino flavor. The weighting function is chosen to be the normalized density fraction for the equilibrium Fermi-Dirac distribution ϵ 2 ⟨ ̂ n ( ϵ ) ⟩ / ( ϵ 2 ⟨ ̂ n ( ϵ ) ⟩ ) max , plotted in Figure 3. Since we have not included flavor oscillations, the µ and τ neutrino distributions are identical to one another, though they are distinct from the electron neutrino distribution.  \nAs shown in Figure 4, this model provides a reasonable fit to the full transport calculation early in the neutrino decoupling epoch, but it fails to accurately capture the occupation probability at the later stages of neutrino decoupling. The evolution of the best-fit parameters for this linear model is shown in Figure 5. We also show the goodness of fit for this model, defined as  \nρ 2 ( T cm ) ≡ 1 -∫ dϵ W ( ϵ ) ( α burst ( ϵ, T cm ) -α fit ( ϵ, T cm )) 2 ∫ dϵ W ( ϵ ) ( α burst ( ϵ, T cm ) -α eq ( ϵ, T cm )) 2 , (7.6)  \nwhere W ( ϵ ) = ϵ 2 ⟨ ̂ n ( ϵ ) ⟩ / ( ϵ 2 ⟨ ̂ n ( ϵ ) ⟩ ) max is the weighting function, taken to be the normalized density fraction plotted in Fig. 3, and the range of integration is taken to be 0 < ϵ < 10. A perfect fit would have α fit = α burst and thus ρ 2 ( T cm ) = 1.  \nOne can see that this linear model provides a reasonable approximation to the evolution of the neutrino occupation probability at high temperature, but it begins to deviate from the numerical results as the temperature drops below about 1 MeV. Next, we will describe a model which provides a better fit to the full transport calculation throughout the entire process of neutrino decoupling.  \nα  \nα  \nFigure 4 . We plot here the evolution of α ( ϵ, T cm ) for electron neutrinos through the weak decoupling epoch. The quantity α is related to the occupation probability through 〈 ̂ n ν f 〉 = [exp( α ) + 1] -1 . We show the difference of α ( ϵ, T cm ) calculated with burst , labeled α burst (solid lines), from the value that it would have in thermal equilibrium, α eq = ϵ . We also show α ( ϵ, T cm ) obtained from weighted linear fit to the results of burst , α fit ( ϵ, T cm ) = α 0 ( T cm ) + α 1 ( T cm ) ϵ (dotted lines). Shaded regions denote the fraction of neutrino density in a given range of ϵ = p/T cm . Note that the 2-parameter model provides a good fit to α ( ϵ ) at high temperatures during the initial phases of weak decoupling, but it begins to deviate at low temperatures during the later stages of weak decoupling.  \n<!-- image -->  \n/epsilon1', '7.2 Quadratic Model for α Evolution': 'Now we define a 3-parameter polynomial model of the form  \nα ( f ) fit = α ( f ) 0 ( T cm ) + α ( f ) 1 ( T cm ) · ϵ + α ( f ) 2 ( T cm ) · ϵ 2 . (7.7)  \nAs in the 2-parameter model described above, the parameter α 0 ( T cm ) is conjugate to the mean number of neutrinos and α 1 ( T cm ) is conjugate to the mean energy. The parameter α 2 ( T cm ) is conjugate to the average square of neutrino energy.  \nJust as was done for the 2-parameter model, we performed a density-weighted fit of the model shown in Eq. (7.7) to the output of the burst code, using the same weighting function. Figure 6 shows that this 3-parameter model provides an excellent fit to the neutrino occupation probability calculated from full transport calculation throughout the epoch of neutrino decoupling. While the goodness of fit for the 2-parameter model was 1 -ρ 2 ≈ 0 . 07, the 3-parameter model fit achieves 1 -ρ 2 < 0 . 0025 throughout the whole epoch of interest, and describes the late-time asymptotic spectrum with 1 -ρ 2 < 0 . 001. The best-fit parameters and goodness of fit as defined in Eq. (7.6) for the quadratic model are shown as a function of T cm in Figure 7. To further  \nFigure 5 . The shape of α ( ϵ, T cm ) can be reasonably approximated with a 2-parameter fit early in the weak decoupling epoch. Here we show the best fit values for the coefficients in a 2-parameter model of α ( ϵ, T cm ) with the form α fit ( ϵ, T cm ) = α 0 ( T cm ) + α 1 ( T cm ) · ϵ . These coefficients are related to the extensive variables of the neutrino ensemble: α 0 is conjugate to the expected number of neutrinos in the system, α 1 is conjugate to the expected energy. The goodness of fit, ρ 2 defined in Eq. (7.6), for this 2-parameter model is also shown as a function of T cm , showing that this model is an excellent fit early in the decoupling epoch, though less so as decoupling progresses.  \n<!-- image -->  \nemphasize how accurately this 3-parameter model describes the neutrino distribution, we show the fractional difference between the best fit model for α and the output of burst in Figure 8 and a comparison of the distribution function, the energy density, and the entropy density in Figure 9.  \nThe quadratic model is particularly successful at describing the neutrino occupation probability after neutrinos have fully decoupled. At late times, each of α ( f ) 0 , 1 -α ( f ) 1 , and α ( f ) 2 have magnitude ≲ 10 -3 . We note, however, that even for small values for α ( f ) 2 , at large ϵ the total α ( f ) will become large in magnitude and negative in sign, giving unphysical occupation probabilities. As a result, we stress that our expansion in three parameters is a good fit for the ϵ range shown in Fig. 6. This motivates a perturbative treatment of the deviation from an equilibrium Fermi-Dirac distribution that we will explore further in the next subsection.  \nOne could also consider a higher-order polynomial fit to the evolution of α . We tested this possibility and found that the α 3 parameter in a cubic fit was significantly smaller than the other  \nα  \nα  \nFigure 6 . Same as Figure 4, but for a 3-parameter fit of the form α ( ϵ, T cm ) = α 0 ( T cm ) + α 1 ( T cm ) ϵ + α 2 ( T cm ) ϵ 2 . Note that the 3-parameter model provides a good fit throughout weak decoupling, especially in the region with the highest number density of neutrinos, and matches very well the asymptotic occupation probability at low temperatures.  \n<!-- image -->  \n/epsilon1  \nα coefficients, and the goodness of fit only marginally improved compared to the quadratic model. For this reason, we focus on the quadratic model in what follows.', '7.3 Perturbed Densities': 'The main cosmological observables resulting from the weak decoupling epoch are the asymptotic BBN products (to be discussed below) and the energy density in neutrinos which is often parameterized in terms of a number N eff , defined in Eq. (1.1). The radiation density as parameterized by N eff affects the expansion rate, and hence affects BBN and the evolution of cosmological perturbations. For the Standard Model, N eff = 3 . 044(1) [21-23, 43] with 3 neutrino species, the excess over 3 a reflection of extra entropy transferred from leptonic weak interactions in the late stages of weak decoupling. Its strongest constraints come from analysis of the CMB power spectrum, and current bounds from Planck give N eff = 2 . 99 ± 0 . 17 (68% CL) [24]. In general, the value of N eff (as well as quantities such as number density, pressure, and entropy density) depends on α after weak decoupling and may differ from the value predicted in the standard scenario.  \nAt a given order of the power series expansion of α in energy, α ν ≈ ∑ M m =0 α m E m , we can relate the α m ( t ) to the mean m -th powers of the energy. Since we have shown M = 2 is all that is needed for an excellent description of α ν ( E,t ), the mean number density, the mean energy density, and the density of the fluctuations in the energy ⟨ δE 2 ⟩ are the quantities of interest. Significant insight can be gained by calculating the perturbations to these quantities that are  \nα  \nFigure 7 . Same as Figure 5 but for a 3-parameter model of the form α fit ( ϵ, T cm ) = α 0 ( T cm )+ α 1 ( T cm ) · ϵ + α 2 ( T cm ) · ϵ 2 . In this model, α 0 , α 1 , α 2 are conjugate to the expected number, energy, and squared-energy, respectively, of the neutrino distributions. The goodness of fit ρ 2 for this model can be compared directly to that of the 2-parameter model shown in Figure 5, demonstrating that the 3-parameter model provides an excellent approximation to the true evolution of α throughout the decoupling epoch.  \n<!-- image -->  \nFigure 8 . Fractional difference between the best fit 3-parameter model describing the occupation probability of electron neutrinos α fit ( ϵ, T cm ) = α 0 ( T cm ) + α 1 ( T cm ) · ϵ + α 2 ( T cm ) · ϵ 2 and the value of α ( ϵ, T cm ) calculated with burst , labeled α burst . While the 3-parameter model is a remarkably good fit for all values of T cm and ϵ , the largest fractional deviations occur for low ϵ at low T cm and high ϵ at intermediate T cm .  \n<!-- image -->  \nlinear in ∆ α νf = α νf -ϵ .  \n〈 ̂ n ν f 〉 = 1 e α ( ϵ,T cm ) +1 ≃ 1 e ϵ +1 -( α ( f ) 0 ( T cm ) + ( α ( f ) 1 ( T cm ) -1) ϵ + α ( f ) 2 ( T cm ) ϵ 2 ) e ϵ ( e ϵ +1) 2 + . . . (7.8) The number density for such a distribution is  \nFor small values of the parameters defining the model for α ( ϵ, T cm ), it is straightforward to calculate the perturbations to quantities like the neutrino energy density, pressure, and entropy. To first order in α 0 , α 1 -1, and α 2 , the distribution function of neutrinos is given by  \nN ν f = g f 2 π 2 T cm 3 ∫ dϵ ϵ 2 ∂ F ∂α 0 = g f 2 π 2 T cm 3 ∫ dϵ ϵ 2 〈 ̂ n ν f 〉 ≃ g f T cm 3 2 π 2 ( 3 ζ (3) 2 -π 2 6 α ( f ) 0 ( T cm ) -9 ζ (3) 2 ( α ( f ) 1 ( T cm ) -1) -7 π 4 30 α ( f ) 2 ( T cm ) ) + . . . , (7.9)  \nthe energy density is  \nρ ν f = g f 2 π 2 T cm 3 ∫ dϵ ϵ 2 ∂ F ∂α 1 = g f 2 π 2 T cm 4 ∫ dϵ ϵ 3 〈 ̂ n ν f 〉 ≃ g f T cm 4 2 π 2 ( 7 π 4 120 -9 ζ (3) 2 α ( f ) 0 ( T cm ) -7 π 4 30 ( α ( f ) 1 ( T cm ) -1) -225 ζ (5) 2 α ( f ) 2 ( T cm ) ) + . . . , (7.10)  \n2  \nFigure 9 . Comparison of the distribution function at three values of ϵ , the energy density, and the entropy density derived from burst and from the quadratic α model. It is evident from these comparisons that the quadratic model provides an excellent fit to the output of burst .  \n<!-- image -->  \nwhile the pressure is P ν f = ρ ν f / 3, and the entropy density is  \ns ν f = -g f 2 π 2 T cm 3 ∫ dϵ ϵ 2 [〈 ̂ n ν f 〉 ln 〈 ̂ n ν f 〉 + ( 1 -〈 ̂ n ν f 〉) ln ( 1 -〈 ̂ n ν f 〉)] ≃ g f T cm 3 2 π 2 ( 7 π 4 90 -9 ζ (3) 2 α ( f ) 0 ( T cm ) -7 π 4 30 ( α ( f ) 1 ( T cm ) -1) -225 ζ (5) 2 α ( f ) 2 ( T cm ) ) + . . . . (7.11)  \nThe perturbative expansions of the above quantities only rely on a single energy scale, namely, T cm . This energy scale is inherent in our dualα description of the distribution functions, as the  \nneutrino energy E ν = ϵ T cm and α is a function of ϵ . We can also determine perturbative expansions of other neutrino-related quantities, although doing so may require additional energy scales indirectly related to T cm . For example, we adopt the standard cosmological parameter N eff ≡ 8 7 ( 11 4 ) 4 / 3 ρ ν ρ γ to characterize the radiation energy density at epochs after electron-positron annihilation. The photon energy density is explicit in the definition of N eff , and therefore the photon temperature T γ is a required ingredient to determine N eff perturbatively with α . More accurately, we need the ratio T cm /T γ to determine N eff . This ratio depends on the transfer of energy from the electromagnetic plasma to the neutrino seas. As we have stressed and explicitly shown in this work, that transfer of energy distorts the neutrino distributions and precipitates non-zero values for α 0 , α 1 -1 , and α 2 . In practice, T cm /T γ cannot be extracted from examining the non-zero α coefficients in isolation of the other components in the system. To obtain a perturbative expansion of N eff with the α coefficients, we must adopt the output ratio T cm /T γ from the burst Boltzmann-transport simulations and use that value to calculate N eff perturbatively. We find T cm /T γ = 0 . 71477 in the burst output for times well after neutrino decoupling, which is slightly larger than the value (4 / 11) 1 / 3 predicted in the simplified case where neutrinos instantaneously decouple prior to electron-positron annihilation. Treating the change to the fourth power of the temperature ratio and the perturbations to the distribution function as being of the same order and working to first order in these quantities, we find that the change to N eff is  \n∆ N eff ≃ ∑ f g f 2 [ ( 11 4 ) 4 / 3 ( T cm T γ ) 4 -1 -540 ζ (3) 7 π 4 α ( f ) 0 ( T cm ) -4( α ( f ) 1 ( T cm ) -1) -13500 ζ (5) 7 π 4 α ( f ) 2 ( T cm ) ] , (7.12)  \nwhich for the best-fit values of the 3-parameter model for α ( ϵ, T cm ) defined in Sec. 7.2 gives N eff = 3 . 034 for T cm < 10 keV. This agrees very well with the value of N eff computed directly from burst as described in detail in Ref. [44], giving further support to the notion that our compact description provides an accurate model of the out-of-equilibrium distributions. If we were to do a full calculation of neutrino decoupling using the α coefficients - and coupling those distributions to the other components of the system - we would follow the ratio T cm /T γ to late times and subsequently deduce the value of N eff , essentially coupling the α evolution to the thermodynamics of the early universe. We leave this calculation to future work.  \nWe can also consider a more general class of quantities that can be computed from our framework. For example the perturbed heat capacity of the cosmic neutrino background is given by  \nC ν f = g f 2 π 2 T cm 3 ∫ dϵ ϵ 2 ( -α 2 1 ∂ 2 F ∂α 2 1 ) = g f 2 π 2 T cm 3 ∫ dϵ ϵ 2 α 2 1 ϵ 2 e α (1 + e α ) 2 ≃ g f T cm 3 2 π 2 ( 7 π 4 30 -18 ζ (3) α ( f ) 0 ( T cm ) -7 π 4 10 ( α ( f ) 1 ( T cm ) -1) -675 ζ (5) α ( f ) 2 ( T cm ) ) + . . . . (7.13)  \nWe have shown in this section that we can achieve a very precise description of the C ν B for a standard thermal history using just three parameters. In the subsequent sections, we further apply our formalism to gain some intuition for how physics beyond the Standard Model may alter the properties of C ν B and its impact on observables.', '8.1 Determining the Differential Visibility': "Processes beyond the Standard Model, such as dark matter decay or flavor oscillation involving sterile neutrino states, may alter the standard neutrino decoupling process and the resulting neutrino occupation probabilities described above. These processes may add or remove neutrinos from the system, ultimately leading to changes in late-time cosmological observables. Disturbances to the neutrino distribution introduced at sufficiently high temperatures would be expected to be smoothed away due to the high rate of weak interactions compared to the expansion rate which would quickly relax the neutrinos to their equilibrium distribution. If, however, a disturbance was introduced at low temperatures, the scattering rates may be small compared to the Hubble rate such that the disturbance is not able to relax, resulting in a deviation from the neutrino occupation probability that may be imprinted in late-time observables. We can achieve a more quantitative description of the conditions on disturbances that leave a lasting imprint on the neutrino occupation probability by drawing an analogy with CMB physics. In photon decoupling, it is convenient to define the optical depth and differential visibility to describe the probability that CMB photons experienced their last scattering at a given redshift. Here we will use a similar treatment to describe the last scattering of neutrinos.  \nFirst, we can express the Boltzmann equation for the neutrinos [Eq. (5.2)] in a slightly more convenient form  \nd dt n ν f = -Γ ν f n ν f + P ν f (1 -n ν f ) , (8.1)  \nwhere Γ ν f is the rate of neutrinos scattering out of a state, and P ν f is the production rate of neutrinos flowing into a state. As in Sec. 5, we will refer to the expectation value of the distribution function 〈 ̂ n ν f 〉 as simply n ν f in this section in order to streamline our notation. Fermi blocking is responsible for the factor (1 -n ν f ) multiplying the production rate, which suppresses the flow of neutrinos into a given state. We introduce the variable Γ ' ν f = Γ ν f + P ν f , which is the effective scattering rate out of states, taking account of Fermi-blocking. The Boltzmann equation can then be rewritten as  \nd dt n ν f = -Γ ' ν f n ν f + P ν f . (8.2)  \nIn this notation, we can relate the probability of neutrino last scattering to the combination Γ ' i H , which gives the effective scattering rate relative to the expansion rate.  \nThe optical depth for a given neutrino species ν f is a summation of the scattering rates along the line of sight, and is given as  \nτ ν f = ∫ a 0 a Γ ' ν f H d ln a ' = ∫ T cm T 0 Γ ' ν f H 1 T cm ' dT cm ' , (8.3)  \nwhere a is the scale factor, and a subscript 0 refers to a quantity evaluated today. The neutrino differential visibility for neutrino species ν f is given by  \nV ν f ≡ d d ln a e -τ ν f = Γ ' ν f H e -τ ν f . (8.4)  \nFigure 2 shows the differential visibility of electron-type neutrinos as a function of T cm and ϵ . The magnitude of the differential visibility provides a measure of probability that a neutrino of comoving energy ϵ last scattered at a particular value of T cm .  \nX  \nFigure 10 . ( Top: ) Differential ν e visibility V ν e as a function of T cm and ϵ . Peak visibility is shown in a black dash-dot line along with contours denoting 50%, 75%, and 95% of peak visibility in solid colored lines. ( Middle: ) Total (free and bound) neutron to proton ratio as a function of T cm for the full BBN calculation compared to the equilibrium expectation given by ( n/p ) eq = exp( -Q/T ). ( Bottom: ) Evolution of nuclide abundances during nucleosynthesis. Solid lines show the evolution of the mass fractions by integrating the nuclear reaction network. Dashed lines show NSE mass fractions, as given by Eq. (9.5).  \n<!-- image -->  \nFigure 2 showcases a novel way of visualizing the weak decoupling of neutrinos from the plasma and provides useful insight about the neutrino freeze-out process. The differential visibility peak defines a border that can be referred to as the neutrino-sphere. Inside the neutrino-sphere (to  \nFigure 11 . Differential visibility for µ - and τ -type neutrinos. As can be seen from the difference in the peak visibility curves, the neutrino-sphere for ν µ ( τ ) appears at a slightly higher T cm compared to ν e , as expected from the larger interaction rate of electron-type neutrinos.  \n<!-- image -->  \nthe left of the peak), scattering rates are large compared to the Hubble expansion rate, so energy injections/non-equilibrium perturbations are easily dispersed across the distribution functions. Outside the neutrino-sphere, scattering rates are small compared to the Hubble rate, so any perturbations or injection of energy cannot be easily dissipated. The latter situation is what causes the late-time neutrino distribution function to deviate from equilibrium, and which can have impacts on late-time observables.  \nThe process of primordial nucleosynthesis provides one example of how neutrino spectral distortions may impact late-time observables. The neutrino differential visibility is shown compared to the evolution of the neutron-to-proton ratio and the light nuclide abundances in Figure 10, to be discussed further in Section 9. One can see from Figure 10 that the broad peak of neutrino last scattering overlaps with the period when the neutron-to-proton ratio begins to deviate from the equilibrium expectation; this deviation plays an essential role in determining the primordial light element abundance yields.  \nDue to the fact that there are charged-current interaction channels for electron-type neutrinos with the plasma that are not present for the ν µ and ν τ , we expect that the differential visibility will differ for the µ and τ neutrinos compared to the electron neutrinos. We compare the differential visibility and neutrino-sphere of ν e and ν µ ( τ ) in Figure 11. Recall that since we neglect flavor oscillations and ν µ undergo the same interactions as ν τ , the scattering rates are identical, and thus their differential visibilities are the same as well. The neutrino-sphere for ν µ ( τ ) appears at slightly larger values of T cm than that of ν e . This is due to the higher scattering rate of the latter as can be seen in Table 1. Only electron-neutrinos have a direct effect on BBN production processes via charged-current interactions, while all neutrino species have an indirect effect through their impact on the Hubble expansion rate.", '8.2 Differential productivity from neutrino decoupling': "Out of equilibrium processes involving neutrinos including, for example, e -+ e + ↔ ν + ¯ ν , cause a transfer of entropy from the photone ± -baryon plasma into the decoupling neutrino seas. This redistribution of energy and degrees of freedom results in a change in the phasing of the temperature, scale factor, and time that, in turn, results in altered light-element abundance yields. We define the differential productivity by rewriting Eq. (8.2) in the form  \nd d ln a ( e -τ ν f n ν f ) = Γ ' ν f H e -τ ν f ∑ r P r Γ ' ν f = Γ ' ν f H e -τ ν f ∑ r n eq P r P tot . (8.5)  \nHere, P r is the production from process r , P tot is the total production rate from all neutrinoscattering processes, and n eq = P tot / Γ ' ν f is the equilibrium distribution function found by setting Eq. (8.2) equal to zero. Using this formalism, we are able to observe how a particular production process, or set of production processes, will impact the distribution functions over time. In other words, we can visualize when a given process will have the most impact on the occupation probabilities.  \nFigure 12 shows the logarithm of the visible production of ν e from the process e -+ e + → ν e +¯ ν e , given by the product of the ν e differential visibility V ν e with P e -e + Γ ' νe . There is also a weighting of ϵ 2 · n ν e ( ϵ ) in order to highlight the region with the highest neutrino number density. The change to the total entropy as a function of T cm is also shown.  \nFigure 12 provides a powerful example of a novel way to visualize processes during the weak decoupling epoch. We can clearly observe the temperatures and energies where the process e -+ e + → ν e + ¯ ν e has the largest impact on the neutrino occupation probability. The differential productivity provides a new tool that can be applied to processes both within the Standard Model and beyond, allowing new insights into how various mechanisms may leave imprints on the late-time neutrino distribution function.", '9 Charged-current weak interactions on neutrons and protons and BBN': "Next, we turn our focus to how neutrinos impact the formation of light nuclei in the early universe. As discussed in Sec. 2, cosmic neutrinos play a significant role in the charged-current, isospin-changing weak interactions that interconvert protons and neutrons. This in turn impacts the abundances of all nuclides in the nuclear network that is ultimately responsible for the light element abundance yields produced through primordial nucleosynthesis.  \nWe begin by giving a qualitative description of the nuclear reaction network based upon an α formalism dual to nuclide abundances; the later method being the standard treatment which the burst code employs. The light element abundance of species J is given by Y J = n J /n b , and similarly Y e for electrons, Y p for protons, and Y n for neutrons. In all cases the equilibrating reactions are fast enough to ensure the respective conjugates to the distribution functions are linear in energy. For species J with momentum p , spin s , and energy E psJ we have α X , p ,t = α 0 J + α 1 J E psJ . The conjugate to the abundance of species J is α 0 = -β J µ J and the conjugate to the energy is of the form α 1 J ≡ β J = T -1 J where µ J is the chemical potential and T J is the temperature of species J . The interactions among all species apart from neutrinos are rapid  \nFigure 12 . ( Top: ) Logarithm of the density-weighted visible production of ν e from e -e + annihilation ( P e -e + ), found by multiplying the differential visibility V ν e by P e -e + Γ ' νe (Γ ' is the effective scattering rate for ν e ) with a weighting of ϵ 2 · n ν e ( ϵ ). We also show contours for 50%, 75%, and 95% of peak visibility for ν e for the sake of comparison. ( Bottom: ) Logarithmic change to the neutrino, plasma, and total entropy as a function of T cm . Note that the derivative of the total entropy is scaled up by a factor of 10 2 relative to the neutrino and plasma entropy derivatives. The purple dashed line shows the logarithmic change to the neutrino entropy as computed from Eq. (7.11) with the 3-parameter fit to α ( ϵ, T cm ) shown in Figure 6, while the solid purple line shows the same quantity calculated directly from burst . The excellent agreement between the dashed line and the solid line further demonstrates that our 3-parameter model provides a precise description of the cosmic neutrino background in the Standard Model.  \n<!-- image -->  \nenough during the BBN epoch to ensure that the temperatures of all nuclides, electrons, and photons are the same. This is the first large reduction in control parameters needed for the reaction networks for abundances and energy.  \nThe neutrino energy-distribution occupation probabilities are almost completely determined by purely leptonic processes, since the number density of nucleons relative to the number density of neutrinos is only ∼ 10 -10 . The leptonically-computed neutrino distributions determine the neutron-to-proton ratio, which in turn influences primordial abundances. There are small perturbative changes in α ν e and α ¯ ν e that accompany these processes and therefore could give insight into the weak interaction history of weak decoupling and BBN.  \nThe neutron-proton interconversion channels relevant to BBN are given in Eqs. (2.5) - (2.7). In burst , the neutron-to-proton ratio along with the nuclear reaction network for the nuclear abundances and the neutrino Boltzmann transport network are all solved together and selfconsistently at each timestep throughout the weak decoupling and BBN regime. The details of the BBN nuclear reaction network are explored in [120], [121], and [122]. The neutron-to-proton  \nratio at any time t , follows from the integration of  \nd ( n/p ) dt = Y p ( λ pe -+ λ p ¯ ν e + λ pe -¯ ν e ) -Y n ( λ nν e + λ ne + + λ n decay ) , (9.1)  \nwhere the reaction rates λ ν e n , λ e + n , λ n decay (forward), and λ pe -, λ p ¯ ν e , λ pe -¯ ν e (reverse) [120] correspond to the reactions in Eqs. (2.5) - (2.7). When the neutron-to-proton interconversion rates are fast compared to the Hubble expansion rate, the instantaneous n/p ratio at time t is near the steady-state solution, namely  \nn p ≃ λ pe -+ λ p ¯ ν e + λ pe -¯ ν e λ ν e n + λ e + n + λ n decay . (9.2)  \nAt high temperatures, where neutrino and charged-lepton distributions have comparable temperatures, Eq. (9.2) becomes coincident with the result from chemical equilibrium [123]  \n( n p ) eq = e -( m n -m p ) /T , (9.3)  \nwhere the neutron-proton mass difference is m n -m p ≃ 1 . 293 MeV. The correction in the neutrino-energy-distribution functions stemming from Boltzmann neutrino transport and out-ofequilibrium neutrino scattering tends to increase the n/p ratio [120].  \nburst actually includes the possibility of any nuclear element of atomic number A J = Z J + N J being present, where Z J is the proton number and N J is the neutron number. At high temperature the nuclear system is in Nuclear Statistical Equilibrium (NSE), since the time scale of the strong nuclear force is much shorter than the characteristic expansion time, so that each reaction is governed by chemical potential balance. In the NSE regime, the many α 0 J describing the abundances of the nuclides may all be related to just the n/p ratio and the total baryon density n b in terms of the binding energy of each nuclide B X , thus representing another major network compression.  \nburst yields the light-element abundances as a function of time. The abundance relative to baryon of nuclear species J (for example, n , p , 2 H, 3 He, etc.) is  \nY J ( t ) = n J ( t ) n b ( t ) , (9.4)  \nand the corresponding mass fraction of element J with atomic number A J is X J = A J Y J . The main asymptotic BBN observables are the abundances of 4 He, deuterium (D), 3 He, and 7 Li. (By convention, the mass fraction of helium, X α , is denoted by Y P , not to be confused with the abundance of protons denoted Y p .) The NSE abundance for nuclide X is  \nY X = Y Z p Y A -Z n 2 ( A -3) / 2 π 3( A -1) / 2 g X A 3 / 2 [ n b ( Tm b ) 3 / 2 ] A -1 e B X /T , (9.5)  \nwhere Z is the atomic number, A is the atomic mass number, Y p is the abundance of free protons, Y n is the abundance of free neutrons, g X is the nuclear partition function for nuclide X , n b is the baryon number density, T is the plasma temperature, m b ≃ 931 . 5 MeV is the baryon rest mass, and B X is the binding energy of nuclide X .  \nEquipped with the knowledge of weak equilibrium and NSE in Eqs. (9.3) and (9.5), we return to discussing Fig. 10. The middle panel of Fig. 10 shows curves for two tracks of n/p . The solid curve gives n/p from an explicit calculation of the n ↔ p rates, encapsulated in Eq. (9.1), whereas the dashed line gives the weak-equilibrium expression from Eq. (9.3). In the bottom panel of Fig. 10, we show the evolution of the mass fractions X J = A J Y J of the nuclides. Similar to the middle panel, solid curves are output from the reaction-network solver, and dashed lines are NSE fractions from Eq. (9.5). We observe that at high temperatures, n/p follows the equilibrium trajectory, and each of the X J for nuclides with A J ≤ 4 follow NSE trajectories ( A J = 6 , 7 are in NSE at higher temperatures T ≃ 10 MeV). In this case, only one chemical potential - and therefore only one energy-averaged ⟨ α b ⟩ for baryons - is needed to characterize the entire set of nuclides. Equation (9.5) simplifies as Y n = Y p e -( m n -m p ) /T . The temperature drops, and the weak-interaction rates no longer keep n/p in weak equilibrium. The departure from weak equilibrium occurs somewhat earlier than what one may guess [30], with deviations between dashed and solid starting as early as T cm = 2MeV, and becoming starkly apparent by 1 MeV. However, strong and electromagnetic reactions keep d , t , 3 He, and 4 He (not visible on the scale of the vertical axis) in NSE with one another. This system at T cm ∼ 1 MeV requires two, independent values {⟨ α n ⟩ , ⟨ α p ⟩} to describe NSE. In other words, Eq. (9.5) requires unique values of Y n and Y p in this temperature range.  \nAs the temperature drops further the network evolution is characterized by fast scrambling among sub-blocks of elements with the interaction between coarser blocks slow enough that chemical potentials are required for each. This is called quasi-nuclear statistical equilibrium (qNSE), familiar from stellar nucleosynthesis when interaction between the iron and silicon blocks is slow. As opposed to a stellar interior, the entropies are so high in the early universe that the abundances in NSE are highest for n and p to epochs well below T cm ∼ 100 keV. As BBN develops, there is an extended freeze-out where progressively more chemical potentials are needed to characterize the network blocks. Nuclides with larger A J depart from NSE earlier, and hence require independent chemical potentials α 0 ,J conjugate to their mass fraction X J . 4 He departs from NSE at T cm ≃ 700 keV and the A J = 3 nuclides at T cm ≃ 200 keV. Deuterium breaks the equilibrium relation last, at a plasma temperature around 80 keV.  \nOut-of-equilibrium neutrino-transport effects result in a decrease in Y p and an increase in deuterium, 3 He, and 7 Li abundances [44]. There are other small corrections that affect late time observables, such as Coulomb corrections and zero-temperature radiative corrections which are fairly independent of each other and of neutrino transport [44]. Other effects, such as finitetemperature radiative corrections, and in-medium renormalization of electron and positron rest masses, may induce nonlinear changes to the abundances. However, all of these conclusions are based on Standard Model physics. BSM physics may alter a number of features of weak decoupling and BBN which can be modeled in high-fidelity with, for example, an appropriately modified burst code. The corresponding alterations in abundances and weak decoupling history can be captured in changes to the generalized entropy parameters, i.e., the α f ( p, t ).", '10 Abundance responses to perturbations': 'By observing light element abundances, we can place constraints on the energy spectrum of cosmic neutrinos. In order to achieve this, we first need to know how these abundances respond to changes to the distribution function of neutrinos. A target of this paper is the linear response of the asymptotic BBN abundances δY J to distortions of the C ν B from the CSM ν . We have only a few light element abundances that we observe (deuterium, helium, and possibly lithium), yet an entire energy and time-dependent set of distorted neutrino distribution functions for the Y J to respond to. At the linear response level, we could determine the neutrino spectrum that is most likely to produce the observed abundances, but about this mean there are very large fluctuations. That is, from observations of light element abundances we can learn only a few modes among the very many modes possible in the full range of time and energy-dependent perturbations to the neutrino spectrum, leaving the spectrum almost completely unconstrained. Using observed Y J in order to make C ν B inferences is a bottom-up approach.  \nThis limitation suggests that a top-down approach is better: we specify a BSM ν model, and compute simultaneously the distorted neutrino distribution functions and the BBN abundances. This is the experimental approach, using numerical calculations like those in burst . What we present here is a framework to understand what is needed in a BSM ν in order to have impact on BBN, but rely on the burst code for the ν -BBN computations which provide the experimental results for our musings. An important aspect of this work is it shows clearly where in the energytime plane perturbations to the CSM ν would be needed to get a measurable BBN response. For example, early energy injection just equilibrates to a Fermi-Dirac distribution and BBN proceeds as in the CSM ν .  \nWe therefore use the burst code [44] to calculate the linear response of light element abundances to small perturbations in the neutrino distribution function. In order to isolate the effect of changes to the spectrum, we take a set of perturbations for which the change to the distribution function of the electron-type neutrinos is equal to the electron-type antineutrinos in order to conserve flavor and lepton numbers. Furthermore, we perturb the distribution function of muontype neutrinos and anti-neutrinos with the same amplitude and opposite sign in order to conserve the total radiation energy density. With these restrictions, we then compute the perturbations to the final primordial abundances (well after BBN) numerically using burst [44].  \nThe perturbations to the occupation numbers, δn , are proportional to the deviations from the zero-chemical-potential FD equilibrium occupation numbers  \nδn ( ϵ, T cm ) ≡ ⟨ ̂ n ⟩ ( ϵ, T cm ) -n (eq) ( ϵ ) n (eq) ( ϵ ) , n (eq) ( ϵ ) = 1 e ϵ +1 . (10.1)  \nPerturbing the ν e , ν e , ν µ , and ν µ occupation numbers in the manner outlined above in Standard Model evolution with out-of-equilibrium neutrino scattering, induces a response in the abundances. Introduction of BSM physics can result in different perturbations in neutrino occupation numbers and, correspondingly, a different response in the abundances. We characterize the abundance responses as a relative difference from a baseline  \nδX ≡ X -X (base) X (base) , (10.2)  \nwhere the baseline abundance X (base) is the abundance calculated when δn = 0 for all ϵ .  \nWe demonstrate the impact of two types of neutrino spectral distortions on the primordial abundances. First, we calculate the impact of perturbations to the neutrino spectra that are persistent and independent of T cm . In this case, the neutrino occupation numbers are perturbed by  \nδn ( ϵ, T cm ) = { ∆ for ϵ ∈ ϵ perturb 0 otherwise (10.3)  \nfor all T cm . The effects of this sort of perturbation for bins in ϵ of width 0 . 15 from 0 to 15 to the primordial helium-4 and deuterium yields are shown in Figures 13 and 14, respectively. Next, we calculate the impact of spectral perturbations that are isolated in both T cm and ϵ , with the form  \nδn ( ϵ, T cm ) = { ∆ for ϵ ∈ ϵ perturb , T cm ∈ T cm , perturb 0 otherwise (10.4)  \nfor bins in ϵ of width 0 . 15 and bins in log 10 ( T cm / MeV) of width 0 . 03. From this latter type of perturbation, we calculate the response function of the primordial abundances to neutrino spectral distortions  \nK X ( ϵ, T cm ) = d 2 δX δn ( ϵ, T cm ) dϵ d log 10 ( T cm / MeV) , (10.5)  \ndefined such that the fractional change to each abundance can be calculated by integrating over the perturbation to the neutrino occupation number  \nδX = ∫ K X ( ϵ, T cm ) δn ( ϵ, T cm ) dϵ d log 10 ( T cm / MeV) . (10.6)  \nTo an excellent approximation, the response functions for each of the primordial abundances can be calculated as a constant scaling of the response function for the neutron-to-proton ratio K X ( ϵ, T cm ) = A X K ( n/p ) ( ϵ, T cm ), shown in Figure 15. The scaling factor A X for each nuclide is given in Table 2.  \nIt should be noted that neither of the types of perturbations to the neutrino spectra discussed in the previous paragraph represent realistic histories of the decoupling epoch, since we did not allow for evolution of the neutrino occupation numbers in the calculations of this section. In reality, if a physical process led to an increased occupation number of neutrinos over some range of ϵ and T cm , weak interactions would cause the perturbation to be redistributed. For example, a localized increase to the neutrino occupation number inserted at high T cm and high ϵ would quickly be smoothed away as the system rapidly approached a thermal distribution. Perturbations at lower T cm , particularly those introduced after the peak of the neutrino visibility function, would persist with less modification due to the low rate of weak interactions in that regime. Despite the artificial nature of the perturbations, their effect on the primordial abundances provides insight into the regions of n ( ϵ, T cm ) to which the primordial abundances are most sensitive.  \nWe observe two regions different from zero in Fig. 15. In the blue region centered at T cm ≃ 0 . 8 MeV and ϵ ≃ 3, K n/p is less than zero. Conversely, in the light red region at T cm ≃ 1 . 0 MeV and ϵ ≃ 7, K n/p is greater than zero. The main contributors to changing n/p in these regions are  \nFigure 13 . Relative deviations in Y P as a function of ϵ for differing changes in occupation number n ( ϵ ). Each color represents a different change to n ( ϵ ). Solid lines correspond to enhancements in n ν,ν e , while dash-dot lines correspond to suppressions in n ν,ν e . The changes to n ( ϵ ) persist for the duration of the calculation.  \n<!-- image -->  \nFigure 14 . Same as Fig. 13 except for relative deviations in D / H  \n<!-- image -->  \nFigure 15 . Response of the total (free and bound) neutron-to-proton ratio after the conclusion of BBN to perturbations to the ν e energy spectrum at specified values of ϵ and T cm in terms of the response function K ( n/p ) defined in Eq. (10.5). The perturbations involve an increase in the occupation number of electron-type neutrinos and antineutrinos, along with an equal and opposite change to the occupation number of muon-type neutrinos and antineutrinos. The changes to the neutron-to-proton ratio shown here are therefore calculated at fixed N eff , and reflect how changes to the energy spectrum of electron neutrinos impact the primordial abundance yields from BBN. The impact on each primordial abundance can be obtained from Eq. (10.6) using the scaling functions shown in Table 2. Also shown are the contours of 50%, 75%, and 95% of peak visibility for electron neutrinos (as shown in Figure 2). The grid on which the response was calculated is overlaid.  \n<!-- image -->  \nthe four processes in Eqs. (2.5) and (2.6) which we repeat below  \nν e + n ↔ p + e -, (10.7)  \ne + + n ↔ p + ¯ ν e , (10.8)  \nwhere the forward rates imply neutron destruction, d ( n/p ) /dt < 0, and the reverse rates imply neutron production, d ( n/p ) /dt > 0. For the first process, abbreviated n ( ν e , e -) p , there are no neutrino energy thresholds. When we increase the distributions by adding δn ( ϵ, T cm ), we enhance the neutron destruction rate, and likewise suppress the neutron production rate via the Pauli blocking factor. The net result is a decrease in n/p as evidenced by the blue region in Fig. 15. For n ( e + , ¯ ν e ) p , there is a neutrino energy threshold in the neutron creation rate, namely m n -m p + m e ∼ 1 . 8 MeV. An increase δn ( ϵ, T cm ) only affects the rates if ϵ > ( m n -m p + m e ) /T cm . If δn ( ϵ, T cm ) occurs above the ϵ threshold, then we enhance the neutron creation rate and suppress the neutron destruction rate, opposite of the effect of δn ( ϵ, T cm ) on n ( ν e , e -) p . However, the enhancement of the neutron production/destruction rates is not equal and does not offset one another. Recall from Eq. (9.1), that d ( n/p ) /dt is a sum of the n ↔ p rates weighted by the abundances of neutrons and proton, Y n and Y p , respectively. For the standard cosmology and  \nTable 2 . Factor relating the response function for the primordial abundance of each nuclide K X ( ϵ, T cm ) defined in Eq. (10.5) to the response function of the neutron-to-proton ratio plotted in Figure 15, such that K X ( ϵ, T cm ) = A X K ( n/p ) ( ϵ, T cm ).  \nperturbations we consider here, Y p ≥ Y n for all times. Therefore, the neutron production rate ( ∼ Y p λ p ¯ ν e ) has a larger effect than the neutron destruction rate ( ∼ Y n λ nν e ). The net result is that the change in n/p is positive, as evidenced by the light-red region in Fig. 15.', '11 Conclusions': 'We introduced a framework which allows for compact and efficient description of the cosmic neutrino background and its deviation from a relativistic Fermi-Dirac distribution. We demonstrated that the spectrum of the cosmic neutrino background in the Standard Model is accurately captured by a three-parameter description of the thermodynamic potential throughout the epoch of weak decoupling and neutrino last scattering. While we neglected neutrino flavor oscillations in the analysis presented here, our framework could be extended to include this and other additional physics.  \nAccurate description of photon decoupling has been essential in extracting detailed information from observations of the cosmic microwave background. The precision with which we can observe the properties of the cosmic neutrino background is unfortunately quite limited compared to the rich data set provided by observations of the cosmic microwave background. Nonetheless, indirect probes of the cosmic neutrino background can be quite incisive, and are becoming ever more so as cosmological observations improve. It is therefore essential that the subtle aspects of neutrino decoupling are treated with sufficient precision to make the most out of the various observable imprints that are left by the influence of the cosmic neutrino background.  \nThe energy density of the cosmic neutrino background affects the expansion rate of the Universe, leaving observable imprints in the CMB anisotropies and the baryon acoustic oscillations observed with galaxy surveys. This arises primarily through the impact of the radiation density, N eff , on the diffusion damping scale [4]. Additionally, fluctuations in the density of neutrinos, which propagate at the speed of light after decoupling, impart a characteristic phase shift of the acoustic peaks appearing in CMB and matter power spectra [110-112]. Both of these effects of the cosmic neutrino background depend only on the integrated energy density, and not on its energy spectrum.  \nThe non-zero mass of neutrinos implies that at least some cosmic neutrinos are non-relativistic and contribute to the matter density of the Universe at late times. Unlike cold dark matter and baryons, cosmic neutrinos retain a large thermal velocity through much of cosmic history, implying that cosmic neutrinos act like hot dark matter. Their presence suppresses the growth of cosmic structure on scales smaller than the neutrino free streaming length compared to a universe containing only massless neutrinos [11-16]. Searches for this suppression of clustering are typically characterized in terms of constraints on the sum of neutrino masses ∑ m ν . However,  \nthis parameterization assumes that cosmic neutrinos are well described by a relativistic FermiDirac distribution. A modified energy spectrum of cosmic neutrinos could manifest itself in suppression of matter clustering that differs from the standard picture, or turning things around, what should be inferred from the measurement of suppressed matter clustering is affected by the energy spectrum of cosmic neutrinos [124].  \nImpacts on the expansion rate and the suppression of matter clustering, typically framed in terms of N eff and ∑ m ν , depend only on the energy spectrum of cosmic neutrinos well after neutrino decoupling. However, cosmic neutrinos play a significant role in BBN in a way that is more sensitive to the evolution of the energy spectrum. Neutrinos participate in the interconversion of neutrons and protons through the weak interaction during the time period of weak decoupling [29]. As a result, the production of light elements during BBN is sensitive to the time evolution of the cosmic neutrino energy spectrum, and not just to the form that it approaches asymptotically at late times.  \nWithin the Standard Model, the cosmic neutrino background deviates only slightly from a relativistic Fermi-Dirac spectrum. However, even in that case, the evolution of the effective temperature of the neutrinos relative to that of the photons evolves in a non-trivial way through weak decoupling and weak freeze-out, due in large part to the fact that electron-positron annihilation overlaps with the period of weak decoupling. Figures 4 and 5 show that this non-trivial evolution cannot be precisely captured through two-parameter fits consistent with extending the GCE to out-of-equilibrium conditions. Conversely, the generalized entropy approach described here provides a compact and precise description of the deviations of the Standard Model C ν B from a Fermi-Dirac distribution using only one additional parameter per species. Furthermore, our formalism provides a physical interpretation that each parameter is conjugate to macroscopic properties of the system, namely: total number of particles, total energy, and fluctuations in energy (i.e., energy-squared).  \nWe demonstrated that the formalism we developed to describe the evolution of the cosmic neutrino background spectrum is sufficiently precise to capture the small time-dependent deviations from an equilibrium distribution that are present in the Standard Model (noting that we neglect flavor oscillations in the present work). An important aspect of the framework we developed is that it provides a valuable set of tools to describe and constrain physics beyond the Standard Model. There are many ways that new physics can impact the neutrino sector, the thermal history of the universe, and the weak decoupling epoch. A very broad class of models of new physics can be usefully mapped onto the description presented here in terms of the evolution of the thermodynamic potential for neutrinos. Our formalism thereby provides a bridge between the diverse phenomenology that may operate in the early universe and the constraints that can be obtained from indirect observations of the cosmic neutrino background. For example, constraints derived from the requirement that BBN not be disrupted by new physics is often treated as a blunt instrument, requiring that the thermal history and neutrino sector are unmodified around temperatures of 0.1 to 1 MeV. The formalism presented here provides a mechanism by which BBN constraints can be made more precise by systematically treating the observational impacts of deviations from the standard scenario.  \nAs a concrete example, consider the out-of-equilibrium decay of a massive particle species around the time of weak decoupling. It is relatively straightforward to compute how this type  \nof physics impacts integrated quantities like N eff , but in general accounting for more detailed evolution must be treated on a case-by-case basis. However, our formalism provides a means by which more general conclusions can be drawn. Any physics which results in a change to the number density, energy density, or energy fluctuations of the cosmic neutrino background can be described in terms of changes to the α parameters we introduced. We can derive observational constraints on these generalized modifications from measurements of primordial light element abundances using the response functions described in Section 10. In this way, very general conclusions can be drawn about how BBN can be used to constrain diverse classes of new physics, without the requirement of carrying out the full calculation of neutrino transport coupled to BBN for each scenario. All that is required is to map the effects of any particular model onto its impact on the α parameters.  \nFurthermore, our description of the temperature and energy dependent cosmic neutrino visibility function gives us a good sense of how new physics operating in various regimes is likely to impact the thermodynamic potential. Injection of high energy neutrinos at early times, well before weak decoupling are likely to become completely thermalized, resulting in little change to the neutrino and plasma distributions described by the standard scenario. Production of neutrinos at late times will simply be added to the standard spectrum, since weak interactions at late times are inefficient and will not redistribute neutrinos in phase space. Energy injection around the time of weak decoupling will result in the new contributions being partially thermalized, contributing mostly to the α 0 , α 1 , and α 2 parameters. Even though we demonstrated this explicitly only in the Standard Model, electron-positron annihilation provides a non-trivial example of this behavior. One may have naively expected that electron-positron annihilation results in a bump in the neutrino spectrum at the value of ϵ associated with the electron mass at the time of annihilation. However, we showed that partial thermalization of the resulting neutrinos allows for a precise description of the resulting neutrino spectrum in terms of just α 0 , α 1 , and α 2 .  \nThere is a close analogy of this description with spectral distortions of the cosmic microwave background. Energy injection during different eras in the early universe tends to result in characteristic µ - and y -distortions of the CMB spectrum [53, 54, 108]. We have shown how Standard Model processes lead to distortions of the C ν B spectrum in terms of α 0 , α 1 , and α 2 . One can use observational constraints on these generalized spectral distortions to gain insight into a broad class of new physics, without the requirement of carrying out a detailed transport calculation for each individual model.  \nThe types of new physics that can be constrained through limits on CMB spectral distortions can also be constrained with primordial abundance measurements through the formalism we developed here, though BBN constraints tend to apply to physics operating at an earlier time ( z ≥ 10 10 for BBN compared to z ≥ 10 3 for CMB spectral distortions). Examples of new physics that can be constrained by both CMB and C ν B spectral distortions include annihilating or decaying particles, evaporating primordial black holes, and dissipation of small-scale acoustic oscillations [53, 54, 125]. Spectral distortions of the cosmic neutrino background may also arise from neutrino self-interactions, right-handed neutrinos, sterile neutrino oscillations, dark matterneutrino interactions, light dark matter freeze-in scenarios, among many others.  \nThe rich landscape of well-motivated new physics scenarios that may impact the cosmic neutrino spectrum calls out for a treatment of the constraints imposed by primordial abundance  \nobservations with more nuance than a simple rigid exclusion over some range of temperatures. The formalism presented here enables this more precise treatment, and it does so in a convenient and compact way that allows for broad and valuable constraining power on models of new physics.', 'Acknowledgments': "The authors would like to thank Daniel Green, Marilena Loverde, Mark Paris, and Nashwan Sabti for helpful discussions. JRB's work was supported by Canada's National Science and Engineering Research Council's Discovery Grant and his Fellowship in the Canadian Institute for Advanced Research. GMF is supported in part by National Science Foundation grant PHY-2209578 at UCSD. EG is supported by the Department of Energy Office of Nuclear Physics award DE-FG0202ER41216. JM is supported by the US Department of Energy under Grant DE-SC0010129. This work was partially enabled by the National Science Foundation under grant No. PHY-2020275: Network for Neutrinos, Nuclear Astrophysics, and Symmetries (N3AS). Additional support was provided by the Heising-Simons foundation under grant No. 2017-228. Computational resources for this research were provided by SMU's Center for Research Computing. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001. GMF, EG, and JRB thank the Institute for Nuclear Theory (INT) at University of Washington, and all authors thank the Canadian Institute for Theoretical Astrophysics (CITA) at University of Toronto for their hospitality while parts of this research were being conducted.", 'A Connected Correlation Functions': "In this appendix, we provide some further details on how the free energy generating functional is related to the connected correlation functions in our generalized entropy formalism.  \nWe start by noting that e -F acts as the generating functional for the correlation functions of the operators Q A . This can be seen as follows  \n̂ e -F = 〈 exp [ -∑ A α ( fi ) A ̂ Q A ]〉 i = 〈 ∞ ∑ j =0 ( -1) j j ! [ ∑ A α ( fi ) A ̂ Q A ] j 〉 i = ∞ ∑ j =0 ( -1) j j ! 〈[ -∑ A α ( fi ) A ̂ Q A ] j 〉 i = 1 -∑ A α ( fi ) A 〈 ̂ Q A 〉 i + ∑ AB α ( fi ) A α ( fi ) B 〈 ̂ Q ( A ̂ Q B ) 〉 i -∑ ABC α ( fi ) A α ( fi ) B α ( fi ) C 〈 ̂ Q ( A ̂ Q B ̂ Q C ) 〉 i + . . . , (A.1)  \nwhere parentheses indicate symmetrization on the indices  \n̂ Q (1 . . . ̂ Q m ) = 1 m ! ∑ σ ̂ Q σ (1) . . . ̂ Q σ ( m ) , (A.2)  \nand σ labels the permutations of the indices. It is then straightforward to see that correlation functions can be obtained from functional derivatives of e -F  \nδ M e -F ( { α ( fi ) A } ) δα m 1 ...δα m M ∣ ∣ ∣ ∣ ∣ α A =0 = ( -1) M 〈 ̂ Q ( m 1 ... ̂ Q m M ) 〉 i . (A.3)  \nWe can then construct the generating functional for the connected correlation functions by taking the logarithm of the correlation function generating functional, which in this case is just the negative of the free energy generating functional  \nF = -ln ( e -F ) = -∞ ∑ j =0 ( -1) j j ! 〈[ -∑ A α ( fi ) A ̂ Q A ] j 〉 i, cc = -1 + ∑ A α ( fi ) A 〈 ̂ Q A 〉 i + ∑ AB α ( fi ) A α ( fi ) B 〈 ̂ Q ( A ̂ Q B ) 〉 i, cc -∑ ABC α ( fi ) A α ( fi ) B α ( fi ) C 〈 ̂ Q ( A ̂ Q B ̂ Q C ) 〉 i, cc + . . . . (A.4)  \nConnected correlation functions can then be computed as functional derivatives of the free energy generating functional  \nδ M F ( { α ( fi ) A } ) δα m 1 ...δα m M ∣ ∣ ∣ ∣ ∣ α A =0 = ( -1) M +1 〈 ̂ Q ( m 1 ... ̂ Q m M ) 〉 i, cc . (A.5)  \nWe can verify that the connected correlation functions denoted here have the expected relationship to the ordinary correlation functions by explicit computation of the first few M -point functions.  \n〈 ̂ Q A 〉 i = -δe -F δα A ∣ ∣ ∣ ∣ α m =0 = e -F δ F δα A ∣ ∣ ∣ ∣ α m =0 = 〈 ̂ Q A 〉 i, cc (A.6)  \n〈 ̂ Q ( A ̂ Q B ) 〉 i = δ 2 e -F δα A δα B ∣ ∣ ∣ ∣ α m =0 = -e -F δ 2 F δα A δα B ∣ ∣ ∣ ∣ α m =0 + e -F δ F δα A δ F δα B ∣ ∣ ∣ ∣ α m =0 = 〈 ̂ Q ( A ̂ Q B ) 〉 i, cc + 〈 ̂ Q A 〉 i, cc 〈 ̂ Q B 〉 i, cc (A.7)  \n〈 ̂ Q ( A ̂ Q B ̂ Q C ) 〉 i = -δ 3 e -F δα A δα B δα C ∣ ∣ ∣ ∣ α m =0 = e -F δ 3 F δα A δα B δα C ∣ ∣ ∣ ∣ α m =0 -e -F δ F δα A δ 2 F δα B δα C ∣ ∣ ∣ ∣ α m =0 -e -F δ F δα B δ 2 F δα C δα A ∣ ∣ ∣ ∣ α m =0 -e -F δ F δα C δ 2 F δα A δα B ∣ ∣ ∣ ∣ α m =0 + e -F δ F δα A δ F δα B δ F δα C ∣ ∣ ∣ ∣ α m =0 = 〈 ̂ Q ( A ̂ Q B ̂ Q C ) 〉 i, cc + 〈 ̂ Q A 〉 i, cc 〈 ̂ Q ( B ̂ Q C ) 〉 i, cc + 〈 ̂ Q B 〉 i, cc 〈 ̂ Q ( C ̂ Q A ) 〉 i, cc + 〈 ̂ Q C 〉 i, cc 〈 ̂ Q ( A ̂ Q B ) 〉 i, cc + 〈 ̂ Q A 〉 i, cc 〈 ̂ Q B 〉 i, cc 〈 ̂ Q C 〉 i, cc (A.8) We can then see that  \n〈 ̂ Q A 〉 i, cc = 〈 ̂ Q A 〉 i 〈 ̂ Q ( A ̂ Q B ) 〉 i, cc = 〈 ̂ Q ( A ̂ Q B ) 〉 i -〈 ̂ Q A 〉 i 〈 ̂ Q B 〉 i 〈 ̂ Q ( A ̂ Q B ̂ Q C ) 〉 i, cc = 〈 ̂ Q ( A ̂ Q B ̂ Q C ) 〉 i -〈 ̂ Q A 〉 i 〈 ̂ Q ( B ̂ Q C ) 〉 i, cc -〈 ̂ Q B 〉 i 〈 ̂ Q ( C ̂ Q A ) 〉 i, cc -〈 ̂ Q C 〉 i 〈 ̂ Q ( A ̂ Q B ) 〉 i, cc -〈 ̂ Q A 〉 i 〈 ̂ Q B 〉 i 〈 ̂ Q C 〉 i = 〈 ̂ Q ( A ̂ Q B ̂ Q C ) 〉 i -〈 ̂ Q A 〉 i 〈 ̂ Q ( B ̂ Q C ) 〉 i -〈 ̂ Q B 〉 i 〈 ̂ Q ( C ̂ Q A ) 〉 i -〈 ̂ Q C 〉 i 〈 ̂ Q ( A ̂ Q B ) 〉 i +2 〈 ̂ Q A 〉 i 〈 ̂ Q B 〉 i 〈 ̂ Q C 〉 i . (A.9)  \nIt is easy to verify for A = B = C , these expressions give the usual relations between moments and cumulants. We note in passing that one could construct a generating functional for the fluctuations in the operators as  \n〈 exp [ -∑ A α ( fi ) A δ ̂ Q A ]〉 i = 〈 exp [ -∑ A α ( fi ) A ( ̂ Q A -〈 ̂ Q A 〉) ]〉 i = ⟨ exp[ -δ ̂ s fi ] ⟩ i , (A.10)  \nand the corresponding generating functional for the connected correlation functions of the fluctuations is then -ln ⟨ exp[ -δ ̂ s fi ] ⟩ i . B Evolution of α coefficients  \nWe start with our covariant form of the Boltzmann equation  \nDn dt = [ ∂ ∂t -Hp ∂ ∂p ] n ( t, p ) = C [ n ( t, p )] , (B.1)  \nwhere the symbol D/dt is a covariant derivative along a worldline defined by the time coordinate. Our first step is to define ϵ ≡ p/T cm , where p is an independent variable, and T cm is a function of time as given in Eq. (5.1)  \nT cm = T in a in a ( t ) , (B.2)  \nfor constants T in and a in . We will change coordinates to the following t ' and ϵ by the transformation  \nt ' = t ; (B.3)  \nϵ = p T cm (B.4)  \nTherefore, the time derivative transforms to  \n∂n ∂t = ∂n ∂t ' ∂t ' ∂t + ∂n ∂ϵ ∂ϵ ∂t = ∂n ∂t ' + ∂n ∂ϵ [ -ϵ T cm ( -HT cm ) ] = ∂n ∂t ' + ϵH ∂n ∂ϵ , (B.5)  \nand the momentum derivative transforms to  \n∂n ∂p = ∂n ∂t ' ∂t ' ∂p + ∂n ∂ϵ ∂ϵ ∂p = 1 T cm ∂n ∂ϵ . (B.6)  \nIn this change of coordinates, the lhs of the Boltzmann equation becomes  \n[ ∂ ∂t -Hp ∂ ∂p ] n ( t, p ) = ∂n ∂t ' + ϵH ∂n ∂ϵ -HϵT cm 1 T cm ∂n ∂ϵ = ∂n ∂t ' (B.7)  \nFor ease in notation, assign t = t ' so that the Boltzmann equation becomes  \nDn dt = ∂n ∂t = C [ n ( t, ϵ )] . (B.8)  \nWith the simplified form of the Boltzmann equation, we will make a transformation into a new dependent variable α = α ( t, ϵ )  \nn = 1 e α +1 = ⇒ α = ln ( 1 -n n ) . (B.9)  \nThe covariant derivative of α with respect to time is  \nDα dt = Dn dt dα dn = -Dn dt 1 n (1 -n ) = -Dn dt e -α ( e α +1) 2 . (B.10)  \nTo solve Eq. (B.10), we write α ( t, ϵ ) in terms of orthogonal polynomials  \nα ( t, ϵ ) = ∑ i a i ( t ) P i ( ϵ ) , (B.11)  \nwhere the coefficients a i ( t ) are time dependent [and bear no relation to scale factor a ( t )], and the polynomials P i ( ϵ ) are time-independent. We have a choice for how we pick the P i . We could use Laguerre polynomials which form an orthogonal basis over the range 0 < ϵ < ∞ . This would be a natural choice given that the range of ϵ is positive real numbers. However, there is very little density at large values of ϵ so capturing the tail of the distribution would not be important for the precision at which we use the quadratic expansion in α . An alternative choice is to use Legendre polynomials for the orthogonal basis. This requires a finite range, i.e., 0 < ϵ < ϵ max and is the protocol burst uses where ϵ max ∼ 20. Depending on the problem at hand, it may behoove us to use one basis over the other. For the problem of weak decoupling in the standard cosmology, we surmise that the Legendre polynomials are a better choice.  \nOnce we pick a set of orthogonal polynomials, we can integrate Eq. (B.10) over ϵ with a weight of a polynomial to find the time evolution of the coefficients  \n̸  \n-∫ dϵ P j ( ϵ ) Dn dt e -α ( e α +1) 2 = ∫ dϵ P j ( ϵ ) Dα dt = ∫ dϵ P j ( ϵ ) D dt [ ∑ i a i ( t ) P i ( ϵ ) ] = ∑ i Da i dt ∫ dϵ P j ( ϵ ) P i ( ϵ ) = ∑ i Da i dt c j δ ij , (B.12)  \nwhere c j = 0 is the coefficient for the particular orthogonal polynomial algebra. For example: c j = 2 / (2 j + 1) for Legendre polynomials (modulo a normalization if the domain has been rescaled). We can solve for the time derivative of the coefficient in the expansion of α  \nDa i dt = -1 c i ∫ dϵ P i ( ϵ ) Dn dt e -α ( e α +1) 2 . (B.13)  \nWehave used the symbol Dn/dt for the collision term implying that we still use the n -distributions to write that term. There is no need to transform Dn/dt into an expression for α : it is the same expression, but wherever we see the symbol n ( ϵ j ), we use 1 / (exp[ α ( ϵ j )] + 1). As a result, we will set A [ α ] ≡ Dn/dt to delineate the difference in notation. A [ α ] is a functional of α and a function of t and ϵ . For clarity, we write the form of A [ α ] explicitly for 2 × 2 neutrino scattering ( r = 1 in Table 1) using the nomenclature of Sec. 5  \nA ( r =1) ν 1 [ α f ] = 1 2 E 1 ∫ d 3 p 2 (2 π ) 3 2 E 2 d 3 p 3 (2 π ) 3 2 E 3 d 3 p 4 (2 π ) 3 2 E 4 (2 π ) 4 δ (4) ( p 1 + p 2 -p 3 -p 4 ) × S r =1 ⟨|M r =1 | 2 ⟩ F r =1 ( p 1 , p 2 , p 3 , p 4 ) , (B.14)  \nwhere in isotropic conditions, the statistical factor becomes  \nF r =1 ( ϵ 1 , ϵ 2 , ϵ 3 , ϵ 4 ) = e α f ( ϵ 1 ) e α f ( ϵ 2 ) -e α f ( ϵ 3 ) e α f ( ϵ 4 ) [ e α f ( ϵ 1 ) +1][ e α f ( ϵ 2 ) +1][ e α f ( ϵ 3 ) +1][ e α f ( ϵ 4 ) +1] , (B.15)  \nwhere the f subscript on α gives the flavor index, and all of the neutrino lines for the r = 1 process in Table 1 have the same flavor.  \nOur final results are time derivatives of the coefficients of the orthogonal polynomials  \nDa i dt = -1 c i ∫ dϵ P i ( ϵ ) A [ α ] e -α ( e α +1) 2 . (B.16)  \nIf we expand α ( t, ϵ ) into a power series, namely  \nα ( t, ϵ ) = ∑ i α i ( t ) ϵ i , (B.17)  \nwe can solve for the power-series coefficients in terms of the orthogonal polynomial coefficients to find the total time derivative for α i  \n∑ i Dα i dt ϵ i = ∑ j Da j dt P j ( ϵ ) , (B.18) = ⇒ ∑ i,k Dα i dt ∫ dϵ P k ( ϵ ) ϵ i = ∑ j,k Da j dt ∫ dϵ P k ( ϵ ) P j ( ϵ ) , = ∑ j,k Da j dt c k δ jk = ∑ k Da k dt c k (B.19) = ⇒ M · D α dt = D a dt (B.20) = ⇒ D α dt = M -1 · D a dt , (B.21)  \nwhere the matrix M is given as  \nM ki = 1 c k ∫ dϵ ϵ i P k ( ϵ ) . (B.22)", 'References': "- [1] M. A. Amin et al., 'Snowmass2021 Theory Frontier White Paper: Data-Driven Cosmology,' in 2022 Snowmass Summer Study. 3, 2022. arXiv:2203.07946 [astro-ph.CO] .\n- [2] M. Gerbino et al., 'Synergy between cosmological and laboratory searches in neutrino physics,' Phys. Dark Univ. 42 (2023) 101333, arXiv:2203.07377 [hep-ph] .\n- [3] C. L. Chang et al., 'Snowmass2021 Cosmic Frontier: Cosmic Microwave Background Measurements White Paper,' in 2022 Snowmass Summer Study. 3, 2022. arXiv:2203.07638 [astro-ph.CO] .\n- [4] C. Dvorkin et al., 'The Physics of Light Relics,' in Snowmass 2021. 3, 2022. arXiv:2203.07943 [hep-ph] .\n- [5] DUNE Collaboration, A. Abed Abud et al., 'Snowmass Neutrino Frontier: DUNE Physics Summary,' arXiv:2203.06100 [hep-ex] .\n- [6] Hyper-Kamiokande Collaboration, J. Bian et al., 'Hyper-Kamiokande Experiment: A Snowmass White Paper,' in Snowmass 2021. 3, 2022. arXiv:2203.02029 [hep-ex] .\n- [7] JUNO Collaboration, F. An et al., 'Neutrino Physics with JUNO,' J. Phys. G 43 no. 3, (2016) 030401, arXiv:1507.05613 [physics.ins-det] .\n- [8] A. de Gouvˆea, 'Neutrino Mass Models,' Annual Review of Nuclear and Particle Science 66 no. 1, (Oct., 2016) 197-217.\n- [9] CMB-S4 Collaboration, K. N. Abazajian et al., 'CMB-S4 Science Book, First Edition,' arXiv:1610.02743 [astro-ph.CO] .\n- [10] K. Abazajian et al., 'CMB-S4 Science Case, Reference Design, and Project Plan,' arXiv:1907.04473 [astro-ph.IM] .\n- [11] A. D. Dolgov, 'Neutrinos in cosmology,' Phys. Rept. 370 (2002) 333-535, arXiv:hep-ph/0202122 .\n- [12] J. Lesgourgues and S. Pastor, 'Massive neutrinos and cosmology,' Phys. Rept. 429 (2006) 307-379, arXiv:astro-ph/0603494 .\n- [13] Y. Y. Y. Wong, 'Neutrino mass in cosmology: status and prospects,' Ann. Rev. Nucl. Part. Sci. 61 (2011) 69-98, arXiv:1111.1436 [astro-ph.CO] .\n- [14] J. Lesgourgues and S. Pastor, 'Neutrino mass from Cosmology,' Adv. High Energy Phys. 2012 (2012) 608515, arXiv:1212.6154 [hep-ph] .\n- [15] C. Dvorkin et al., 'Neutrino Mass from Cosmology: Probing Physics Beyond the Standard Model,' arXiv:1903.03689 [astro-ph.CO] .\n- [16] D. Green and J. Meyers, 'Cosmological Implications of a Neutrino Mass Detection,' arXiv:2111.01096 [astro-ph.CO] .\n- [17] R. B. Patterson, 'Prospects for Measurement of the Neutrino Mass Hierarchy,' Ann. Rev. Nucl. Part. Sci. 65 (2015) 177-192, arXiv:1506.07917 [hep-ex] .\n- [18] KATRIN Collaboration, M. Aker et al., 'KATRIN: Status and Prospects for the Neutrino Mass and Beyond,' arXiv:2203.08059 [nucl-ex] .\n- [19] Project 8 Collaboration, A. A. Esfahani et al., 'The Project 8 Neutrino Mass Experiment,' in 2022 Snowmass Summer Study. 3, 2022. arXiv:2203.07349 [nucl-ex] .  \n- [20] nEXO Collaboration, G. Adhikari et al., 'nEXO: neutrinoless double beta decay search beyond 10 28 year half-life sensitivity,' J. Phys. G 49 no. 1, (2022) 015104, arXiv:2106.16243 [nucl-ex] .\n- [21] K. Akita and M. Yamaguchi, 'A precision calculation of relic neutrino decoupling,' JCAP 08 (2020) 012, arXiv:2005.07047 [hep-ph] .\n- [22] J. Froustey, C. Pitrou, and M. C. Volpe, 'Neutrino decoupling including flavour oscillations and primordial nucleosynthesis,' JCAP 12 (2020) 015, arXiv:2008.01074 [hep-ph] .\n- [23] J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor, and Y. Y. Y. Wong, 'Towards a precision calculation of N eff in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED,' JCAP 04 (2021) 073, arXiv:2012.02726 [hep-ph] .\n- [24] Planck Collaboration, N. Aghanim et al., 'Planck 2018 results. VI. Cosmological parameters,' Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO] . [Erratum: Astron.Astrophys. 652, C4 (2021)].\n- [25] T.-H. Yeh, J. Shelton, K. A. Olive, and B. D. Fields, 'Probing physics beyond the standard model: limits from BBN and the CMB independently and combined,' JCAP 10 (2022) 046, arXiv:2207.13133 [astro-ph.CO] .\n- [26] Simons Observatory Collaboration, P. Ade et al., 'The Simons Observatory: Science goals and forecasts,' JCAP 02 (2019) 056, arXiv:1808.07445 [astro-ph.CO] .\n- [27] NASA PICO Collaboration, S. Hanany et al., 'PICO: Probe of Inflation and Cosmic Origins,' arXiv:1902.10541 [astro-ph.IM] .\n- [28] N. Sehgal et al., 'CMB-HD: An Ultra-Deep, High-Resolution Millimeter-Wave Survey Over Half the Sky,' arXiv:1906.10134 [astro-ph.CO] .\n- [29] R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, 'Big Bang Nucleosynthesis: 2015,' Rev. Mod. Phys. 88 (2016) 015004, arXiv:1505.01076 [astro-ph.CO] .\n- [30] E. Grohs and G. M. Fuller, 'The surprising influence of late charged current weak interactions on Big Bang Nucleosynthesis,' Nuclear Physics B 911 (Oct., 2016) 955-973, arXiv:1607.02797 [astro-ph.CO] .\n- [31] G. M. Fuller, C. T. Kishimoto, and A. Kusenko, 'Heavy sterile neutrinos, entropy and relativistic energy production, and the relic neutrino background,' arXiv e-prints (Oct., 2011) arXiv:1110.6479, arXiv:1110.6479 [astro-ph.CO] .\n- [32] H. Rasmussen, A. McNichol, G. M. Fuller, and C. T. Kishimoto, 'Effects of an intermediate mass sterile neutrino population on the early Universe,' Phys. Rev. D 105 no. 8, (2022) 083513, arXiv:2109.11176 [hep-ph] .\n- [33] K. Abazajian, N. F. Bell, G. M. Fuller, and Y. Y. Y. Wong, 'Cosmological lepton asymmetry, primordial nucleosynthesis and sterile neutrinos,' Phys. Rev. D 72 no. 6, (Sept., 2005) 063004, arXiv:astro-ph/0410175 [astro-ph] .\n- [34] D. A. Dicus, E. W. Kolb, A. M. Gleeson, E. C. G. Sudarshan, V. L. Teplitz, and M. S. Turner, 'Primordial Nucleosynthesis Including Radiative, Coulomb, and Finite Temperature Corrections to Weak Rates,' Phys. Rev. D 26 (1982) 2694.\n- [35] J.-L. Cambier, J. R. Primack, and M. Sher, 'Finite Temperature Radiative Corrections to Neutron Decay and Related Processes,' Nucl. Phys. B 209 (1982) 372. [Erratum: Nucl.Phys.B 222, 517-517 (1983)].  \n- [36] N. C. Rana and B. M. Seifert, 'Effect of neutrino heating in the early universe on neutrino decoupling temperatures and nucleosynthesis,' Phys. Rev. D 44 (1991) 393-397.\n- [37] S. Dodelson and M. S. Turner, 'Nonequilibrium neutrino statistical mechanics in the expanding universe,' Phys. Rev. D 46 (1992) 3372-3387.\n- [38] S. Hannestad and J. Madsen, 'Neutrino decoupling in the early universe,' Phys. Rev. D 52 (1995) 1764-1769, arXiv:astro-ph/9506015 .\n- [39] A. D. Dolgov, S. H. Hansen, and D. V. Semikoz, 'Nonequilibrium corrections to the spectra of massless neutrinos in the early universe,' Nucl. Phys. B 503 (1997) 426-444, arXiv:hep-ph/9703315 .\n- [40] N. Y. Gnedin and O. Y. Gnedin, 'Cosmological neutrino background revisited,' Astrophys. J. 509 (1998) 11-15, arXiv:astro-ph/9712199 .\n- [41] A. D. Dolgov, S. H. Hansen, and D. V. Semikoz, 'Nonequilibrium corrections to the spectra of massless neutrinos in the early universe: Addendum,' Nucl. Phys. B 543 (1999) 269-274, arXiv:hep-ph/9805467 .\n- [42] J. Birrell, C.-T. Yang, P. Chen, and J. Rafelski, 'Relic neutrinos: Physically consistent treatment of effective number of neutrinos and neutrino mass,' Phys. Rev. D 89 no. 2, (Jan., 2014) 023008, arXiv:1212.6943 [astro-ph.CO] .\n- [43] M. Escudero Abenza, 'Precision early universe thermodynamics made simple: N eff and neutrino decoupling in the Standard Model and beyond,' JCAP 05 (2020) 048, arXiv:2001.04466 [hep-ph] .\n- [44] E. Grohs, G. M. Fuller, C. T. Kishimoto, M. W. Paris, and A. Vlasenko, 'Neutrino energy transport in weak decoupling and big bang nucleosynthesis,' Phys. Rev. D93 no. 8, (2016) 083522, arXiv:1512.02205 [astro-ph.CO] .\n- [45] A.-K. Burns, T. M. P. Tait, and M. Valli, 'PRyMordial: the first three minutes, within and beyond the standard model,' European Physical Journal C 84 no. 1, (Jan., 2024) 86, arXiv:2307.07061 [hep-ph] .\n- [46] KATRIN Collaboration, M. Aker et al., 'Direct neutrino-mass measurement with sub-electronvolt sensitivity,' Nature Phys. 18 no. 2, (2022) 160-166, arXiv:2105.08533 [hep-ex] .\n- [47] P. J. E. Peebles, 'Recombination of the Primeval Plasma,' Astrophys. J. 153 (1968) 1.\n- [48] Y. B. Zeldovich, V. G. Kurt, and R. A. Sunyaev, 'Recombination of hydrogen in the hot model of the universe,' Zh. Eksp. Teor. Fiz. 55 (1968) 278-286.\n- [49] S. Seager, D. D. Sasselov, and D. Scott, 'A new calculation of the recombination epoch,' Astrophys. J. Lett. 523 (1999) L1-L5, arXiv:astro-ph/9909275 .\n- [50] S. Seager, D. D. Sasselov, and D. Scott, 'How exactly did the universe become neutral?,' Astrophys. J. Suppl. 128 (2000) 407-430, arXiv:astro-ph/9912182 .\n- [51] Y. Ali-Haimoud and C. M. Hirata, 'HyRec: A fast and highly accurate primordial hydrogen and helium recombination code,' Phys. Rev. D 83 (2011) 043513, arXiv:1011.3758 [astro-ph.CO] .\n- [52] J. Chluba and R. M. Thomas, 'Towards a complete treatment of the cosmological recombination problem,' Mon. Not. Roy. Astron. Soc. 412 (2011) 748, arXiv:1010.3631 [astro-ph.CO] .\n- [53] J. Chluba and R. A. Sunyaev, 'The evolution of CMB spectral distortions in the early Universe,' Mon. Not. Roy. Astron. Soc. 419 (2012) 1294-1314, arXiv:1109.6552 [astro-ph.CO] .  \n- [54] J. Chluba, 'Which spectral distortions does ΛCDM actually predict?,' Mon. Not. Roy. Astron. Soc. 460 no. 1, (2016) 227-239, arXiv:1603.02496 [astro-ph.CO] .\n- [55] F. Takahara and S. Sasaki, 'Effect of Rayleigh Scattering on Cosmic Microwave Background Anisotropies,' Progress of Theoretical Physics 86 no. 5, (Nov., 1991) 1021-1030.\n- [56] Q.-J. Yu, D. N. Spergel, and J. P. Ostriker, 'Rayleigh scattering and microwave background fluctuations,' Astrophys. J. 558 (2001) 23-28, arXiv:astro-ph/0103149 .\n- [57] A. Lewis, 'Rayleigh scattering: blue sky thinking for future CMB observations,' JCAP 08 (2013) 053, arXiv:1307.8148 [astro-ph.CO] .\n- [58] E. Alipour, K. Sigurdson, and C. M. Hirata, 'Effects of Rayleigh scattering on the CMB and cosmic structure,' Phys. Rev. D 91 no. 8, (2015) 083520, arXiv:1410.6484 [astro-ph.CO] .\n- [59] B. Beringue, P. D. Meerburg, J. Meyers, and N. Battaglia, 'Cosmology with Rayleigh Scattering of the Cosmic Microwave Background,' JCAP 01 (2021) 060, arXiv:2008.11688 [astro-ph.CO] .\n- [60] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, 'The Cosmic Microwave Background spectrum from the full COBE FIRAS data set,' Astrophys. J. 473 (1996) 576, arXiv:astro-ph/9605054 .\n- [61] G. Mangano, G. Miele, S. Pastor, and M. Peloso, 'A Precision calculation of the effective number of cosmological neutrinos,' Phys. Lett. B 534 (2002) 8-16, arXiv:astro-ph/0111408 .\n- [62] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico, 'Relic neutrino decoupling including flavor oscillations,' Nucl. Phys. B 729 (2005) 221-234, arXiv:hep-ph/0506164 .\n- [63] C. Volpe, D. Vaananen, and C. Espinoza, 'Extended evolution equations for neutrino propagation in astrophysical and cosmological environments,' Phys. Rev. D 87 no. 11, (June, 2013) 113010, arXiv:1302.2374 [hep-ph] .\n- [64] Y. Zhang and A. Burrows, 'Transport equations for oscillating neutrinos,' Phys. Rev. D 88 no. 10, (Nov., 2013) 105009, arXiv:1310.2164 [hep-ph] .\n- [65] A. Vlasenko, G. M. Fuller, and V. Cirigliano, 'Neutrino quantum kinetics,' Phys. Rev. D 89 no. 10, (May, 2014) 105004, arXiv:1309.2628 [hep-ph] .\n- [66] J. Froustey and C. Pitrou, 'Primordial neutrino asymmetry evolution with full mean-field effects and collisions,' arXiv e-prints (Oct., 2021) arXiv:2110.11889, arXiv:2110.11889 [hep-ph] .\n- [67] M. J. Savage, R. A. Malaney, and G. M. Fuller, 'Neutrino Oscillations and the Leptonic Charge of the Universe,' ApJ 368 (Feb., 1991) 1.\n- [68] K. N. Abazajian, J. F. Beacom, and N. F. Bell, 'Stringent constraints on cosmological neutrino-antineutrino asymmetries from synchronized flavor transformation,' Phys. Rev. D 66 no. 1, (July, 2002) 013008, arXiv:astro-ph/0203442 [astro-ph] .\n- [69] A. D. Dolgov, S. H. Hansen, S. Pastor, S. T. Petcov, G. G. Raffelt, and D. V. Semikoz, 'Cosmological bounds on neutrino degeneracy improved by flavor oscillations,' Nuclear Physics B 632 no. 1-3, (June, 2002) 363-382, arXiv:hep-ph/0201287 [hep-ph] .\n- [70] K. Abazajian, G. M. Fuller, and M. Patel, 'Sterile neutrino hot, warm, and cold dark matter,' Phys. Rev. D 64 no. 2, (July, 2001) 023501, arXiv:astro-ph/0101524 [astro-ph] .\n- [71] A. D. Dolgov and S. H. Hansen, 'Massive sterile neutrinos as warm dark matter,' Astroparticle Physics 16 no. 3, (Jan., 2002) 339-344, arXiv:hep-ph/0009083 [hep-ph] .  \n- [72] C. T. Kishimoto, G. M. Fuller, and C. J. Smith, 'Coherent Active-Sterile Neutrino Flavor Transformation in the Early Universe,' Phys. Rev. Lett. 97 no. 14, (Oct., 2006) 141301, arXiv:astro-ph/0607403 [astro-ph] .\n- [73] M. Drewes et al., 'A White Paper on keV Sterile Neutrino Dark Matter,' JCAP 01 (2017) 025, arXiv:1602.04816 [hep-ph] .\n- [74] A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens, and O. Ruchayskiy, 'Sterile neutrino Dark Matter,' Progress in Particle and Nuclear Physics 104 (Jan., 2019) 1-45, arXiv:1807.07938 [hep-ph] .\n- [75] L. Johns and G. M. Fuller, 'Self-interacting sterile neutrino dark matter: The heavy-mediator case,' Phys. Rev. D 100 no. 2, (July, 2019) 023533, arXiv:1903.08296 [hep-ph] .\n- [76] A. de Gouvˆea, M. Sen, W. Tangarife, and Y. Zhang, 'Dodelson-Widrow Mechanism in the Presence of Self-Interacting Neutrinos,' Phys. Rev. Lett. 124 no. 8, (Feb., 2020) 081802, arXiv:1910.04901 [hep-ph] .\n- [77] G. B. Gelmini, P. Lu, and V. Takhistov, 'Cosmological dependence of resonantly produced sterile neutrinos,' J. Cosmology Astropart. Phys. 2020 no. 6, (June, 2020) 008, arXiv:1911.03398 [hep-ph] .\n- [78] G. B. Gelmini, M. Kawasaki, A. Kusenko, K. Murai, and V. Takhistov, 'Big bang nucleosynthesis constraints on sterile neutrino and lepton asymmetry of the Universe,' J. Cosmology Astropart. Phys. 2020 no. 9, (Sept., 2020) 051, arXiv:2005.06721 [hep-ph] .\n- [79] G. B. Gelmini, P. Lu, and V. Takhistov, 'Visible sterile neutrinos as the earliest relic probes of cosmology,' Physics Letters B 800 (Jan., 2020) 135113, arXiv:1909.04168 [hep-ph] .\n- [80] G. M. Fuller, L. Gr'af, A. V. Patwardhan, and J. Spisak, 'Dark population transfer mechanism for sterile neutrino dark matter,' arXiv e-prints (Feb., 2024) arXiv:2402.13878, arXiv:2402.13878 [hep-ph] .\n- [81] S. Gariazzo, P. F. de Salas, and S. Pastor, 'Thermalisation of sterile neutrinos in the early Universe in the 3+1 scheme with full mixing matrix,' JCAP 07 (2019) 014, arXiv:1905.11290 [astro-ph.CO] .\n- [82] R. J. Scherrer and M. S. Turner, 'Primordial Nucleosynthesis with Decaying Particles. 1. Entropy Producing Decays. 2. Inert Decays,' Astrophys. J. 331 (1988) 19-32.\n- [83] R. H. Cyburt, J. R. Ellis, B. D. Fields, and K. A. Olive, 'Updated nucleosynthesis constraints on unstable relic particles,' Phys. Rev. D 67 (2003) 103521, arXiv:astro-ph/0211258 .\n- [84] A. V. Patwardhan, G. M. Fuller, C. T. Kishimoto, and A. Kusenko, 'Diluted equilibrium sterile neutrino dark matter,' Phys. Rev. D 92 no. 10, (Nov., 2015) 103509, arXiv:1507.01977 [astro-ph.CO] .\n- [85] M. Kawasaki, K. Kohri, T. Moroi, and Y. Takaesu, 'Revisiting Big-Bang Nucleosynthesis Constraints on Long-Lived Decaying Particles,' Phys. Rev. D 97 no. 2, (2018) 023502, arXiv:1709.01211 [hep-ph] .\n- [86] K. N. Abazajian and H. Garc'ıa Escudero, 'Visible in the laboratory and invisible in cosmology: Decaying sterile neutrinos,' Phys. Rev. D 108 no. 12, (Dec., 2023) 123036, arXiv:2309.11492 [hep-ph] .\n- [87] N. F. Bell, E. Pierpaoli, and K. Sigurdson, 'Cosmological signatures of interacting neutrinos,' Phys. Rev. D 73 (2006) 063523, arXiv:astro-ph/0511410 .  \n- [88] F.-Y. Cyr-Racine and K. Sigurdson, 'Limits on Neutrino-Neutrino Scattering in the Early Universe,' Phys. Rev. D 90 no. 12, (2014) 123533, arXiv:1306.1536 [astro-ph.CO] .\n- [89] Neutrino Non-Standard Interactions: A Status Report, vol. 2. 2019. arXiv:1907.00991 [hep-ph] .\n- [90] E. Grohs, G. M. Fuller, and M. Sen, 'Consequences of neutrino self-interactions for weak decoupling and big bang nucleosynthesis,' J. Cosmology Astropart. Phys. 2020 no. 7, (July, 2020) 001, arXiv:2002.08557 [astro-ph.CO] .\n- [91] C. Chichiri, G. B. Gelmini, P. Lu, and V. Takhistov, 'Cosmological dependence of sterile neutrino dark matter with self-interacting neutrinos,' J. Cosmology Astropart. Phys. 2022 no. 9, (Sept., 2022) 036, arXiv:2111.04087 [hep-ph] .\n- [92] G. Mangano, A. Melchiorri, P. Serra, A. Cooray, and M. Kamionkowski, 'Cosmological bounds on dark matter-neutrino interactions,' Phys. Rev. D 74 (2006) 043517, arXiv:astro-ph/0606190 .\n- [93] R. Laha, B. Dasgupta, and J. F. Beacom, 'Constraints on New Neutrino Interactions via Light Abelian Vector Bosons,' Phys. Rev. D 89 no. 9, (2014) 093025, arXiv:1304.3460 [hep-ph] .\n- [94] A. Olivares-Del Campo, C. Bœhm, S. Palomares-Ruiz, and S. Pascoli, 'Dark matter-neutrino interactions through the lens of their cosmological implications,' Phys. Rev. D 97 no. 7, (2018) 075039, arXiv:1711.05283 [hep-ph] .\n- [95] J.-T. Li, G. M. Fuller, and E. Grohs, 'Probing dark photons in the early universe with big bang nucleosynthesis,' J. Cosmology Astropart. Phys. 2020 no. 12, (Dec., 2020) 049, arXiv:2009.14325 [astro-ph.CO] .\n- [96] S. Esposito, G. Miele, S. Pastor, M. Peloso, and O. Pisanti, 'Nonequilibrium spectra of degenerate relic neutrinos,' Nucl. Phys. B 590 (2000) 539-561, arXiv:astro-ph/0005573 .\n- [97] H. Groenewold, 'On the Principles of elementary quantum mechanics,' Physica 12 (1946) 405-460.\n- [98] J. Moyal, 'Quantum mechanics as a statistical theory,' Proc. Cambridge Phil. Soc. 45 (1949) 99-124.\n- [99] C. K. Zachos, D. B. Fairlie, and T. L. Curtright, eds., Quantum Mechanics in Phase Space. World Scientific Publishing Co Pte Ltd, 12, 2005.\n- [100] S. P. Mikheyev and A. Y. Smirnov, 'Resonance enhancement of oscillations in matter and solar neutrino spectroscopy,' Yadernaya Fizika 42 (Jan., 1985) 1441-1448.\n- [101] L. Wolfenstein, 'Neutrino oscillations in matter,' Phys. Rev. D 17 no. 9, (May, 1978) 2369-2374.\n- [102] L. Wolfenstein, 'Neutrino oscillations and stellar collapse,' Phys. Rev. D 20 no. 10, (Nov., 1979) 2634-2635.\n- [103] Y. B. Zeldovich and R. A. Sunyaev, 'The Interaction of Matter and Radiation in a Hot-Model Universe,' Astrophys. Space Sci. 4 (1969) 301-316.\n- [104] R. A. Sunyaev and Y. B. Zeldovich, 'The Interaction of matter and radiation in the hot model of the universe,' Astrophys. Space Sci. 7 (1970) 20-30.\n- [105] C. Burigana, L. Danese, and G. De Zotti, 'Formation and evolution of early distortions of the microwave background spectrum - A numerical study,' Astron. Astrophys. 246 no. 1, (1991) 49-58.\n- [106] W. Hu and J. Silk, 'Thermalization and spectral distortions of the cosmic background radiation,' Phys. Rev. D 48 (1993) 485-502.  \n- [124] J. Alvey, M. Escudero, and N. Sabti, 'What can CMB observations tell us about the neutrino distribution function?,' JCAP 02 no. 02, (2022) 037, arXiv:2111.12726 [astro-ph.CO] .\n- [125] D. Jeong, J. Pradler, J. Chluba, and M. Kamionkowski, 'Silk damping at a redshift of a billion: a new limit on small-scale adiabatic perturbations,' Phys. Rev. Lett. 113 (2014) 061301, arXiv:1403.3697 [astro-ph.CO] ."}