Relic Neutrino Freeze-out: Dependence on Natural Constants

Analysis of cosmic microwave background radiation fluctuations favors an effective number of neutrinos, $N_\nu>3$. This motivates a reinvestigation of the neutrino freeze-out process. Here we characterize the dependence of $N_\nu$ on the Standard Model (SM) parameters that govern neutrino freeze-out. We show that $N_\nu$ depends on a combination $\eta$ of several natural constants characterizing the relative strength of weak interaction processes in the early Universe and on the Weinberg angle $\sin^2\theta_W$. We determine numerically the dependence $N_\nu(\eta,\sin^2\theta_W)$ and discuss these results. The extensive numerical computations are made possible by two novel numerical procedures: a spectral method Boltzmann equation solver adapted to allow emerging chemical non-equilibrium, and a method to evaluate Boltzmann equation collision integrals that generates a smooth integrand.


Introduction
The relic neutrino background is believed to be a well preserved probe of a Universe only a second old. The properties of the neutrino background are influenced by the details of the freeze-out or decoupling process at a temperature T = O(1 MeV), which is in turn controlled by the standard model (SM) of particle physics parameters. In this paper, we study the influence of SM parameters on the neutrino distribution after freeze-out. This exercise is of interest because: • There is a tension between the known number of neutrinos (three flavors) and the effective number of neutrinos N ν ≃ 3.5 deduced from the study of the cosmic microwave background (CMB). Detailed analysis of the CMB by the Planck satellite collaboration (Planck) [1] tests both gravitational and standard model (SM) interactions in the early Universe . So far very little effort has been devoted to the understanding how these results characterize SM properties in the early Universe.
• The topic of time variation of natural constants is a very active field with a long history [2]. Consideration of neutrino freeze-out dependence on natural constants provides new insights on the time and/or temperature variation of several SM parameters considered at in early era O(1)s in the Universe's evolution.
The presence of relativistic particles, such as neutrinos, strongly impacts the dynamics of the expansion of the Universe, constrained by direct measurement and analysis of cosmic microwave background (CMB) temperature fluctuations [1]. The effect is described by N ν which quantifies the amount of radiation energy density, ρ r , in the Universe prior to photon freeze-out and after e ± annihilation and is defined by where ρ γ is the photon energy density. The factor 7/8 is the ratio of Fermi to Bose normalization in ρ and the neutrino to photon temperature ratio R ν is the result of the transfer of e + e − entropy into photons after Standard Model (SM) left handed neutrino freeze-out: This well known value R 0 ν arises in the limit where no entropy from the annihilating e ± pairs is transferred to neutrinos, i.e. all entropy feeds and reheats the photon background.
We expect to measure a value N ν = 3, i.e. the number of SM left handed neutrino flavors only if: 1. Photons and are the only other effectively massless particle species in the Universe between the freeze-out of the left handed neutrinos at T γ = O(1) MeV and photon freeze-out at T γ = 0.25 eV, and 2. No flow of entropy from e ± annihilation to neutrinos occurs. The current status is that computation of the neutrino freeze-out process employing SM two body scattering interactions and carried out using the Boltzmann equation gives N th ν = 3.046 [3], a value close to the number of flavors.
The value of N ν can be measured by fitting to observational data, such as the distribution of CMB temperature fluctuations. The Planck [1] analysis gives N ν = 3.36 ± 0.34 (CMB only), N ν = 3.30 ± 0.27 (CMB+BAO), and N ν = 3.62 ± 0.25 (CMB+H 0 ) (68% confidence levels). Combinations of of Planck results with other priors are also reported by the Planck collaboration with most resulting in central values N ν ∈ (3.3, 3.6). Efforts to reconcile Planck results with the recent BICEP2 polarization results also pushes N ν up [4,5] and values N ν > 4 are considered. With more dedicated CMB experiments on the drawing board it is believed that a significantly more precise value of N ν is forthcoming in the next decade.
The tension between the values inferred from observation and the SM prediction has inspired various theories, including the consideration of: modified neutrino interactions [6]; a model in which the temperature of decoupling was a model parameter [7]; a model of a new spontaneously broken symmetry associated with massless Goldstone bosons that freeze out prior to the disappearance of muons [8], and a similar consideration motivated by recognition that the quark-gluon plasma phase transition offers the required physics context [9].
In this paper we explore the dependence of N ν on the value of natural constants within realm of known interactions. We show that N ν depends only on the magnitude of the Weinberg angle in the form sin 2 θ W , and a dimensionless relative interaction strength parameter η, a combination of the electron mass m e , Newton constant G N , and the Fermi constant G F . The magnitude of sin 2 θ W is not fixed within the SM and could be subject to variation as function of time or temperature. We find that sin 2 θ W which differs substantially from the value measured today in vacuum is capable of significantly altering the value of N ν that is generated during neutrino freeze-out. We further show the combined effect of modified η and sin 2 θ W and argue that this can remove or at least reduce the tension of N ν with the Planck data.
In section 2 we introduce the two body scattering description of neutrino freeze-out. In particular we discuss the Boltzmann equation, which is used to model the neutrino freeze-out process, and present the matrix elements that control neutrino freeze-out in subsection 2.1. We then discuss in subsections 2.2 and 2.3 the dependence of the Boltzmann equation on SM parameters, the Weinberg angle sin 2 θ W , and the interaction strength parameter η respectively. In section 3 we introduce the technical methods we use to solve the Boltzmann equation. In subsection 3.1 we outline our solution method for the Boltzmann equation. In subsection 3.2 we detail a new method for analytically simplifying the collision integrals in order to reduce the numerical integration costs. In subsection 3.4 we compare with the results by previous authors, highlighting the improvements we have made.
In section 4 we show the impact of SM parameter values on neutrino freeze-out and the effective number of neutrinos. In particular, we show how the impact of the strength parameter η and sin 2 θ W on N ν . We give our concluding analysis in section 5. Appendix A contains some mathematical background that is useful for the collision integral calculations of subsection 3.2 and in Appendix B we apply the method to the processes involved in neutrino freeze-out. Finally, Appendix C contains additional plots and numerical fits that show the impact of SM parameters on various quantities characterizing the neutrino distributions after freeze-out.

Einstein-Boltzmann Equation
To model the flow of energy and entropy into the relic neutrino distribution, and hence obtain the value of N ν after freeze-out, we must solve the general relativistic Einstein-Boltzmann equation. Several references discuss this generalization of the Boltzmann equation in detail [10,11,12,13,14,15] and other works specialize this to the question of neutrino freeze-out [16,17,18,19,20,3]. Here we provide a quick overview of this literature. In the context of general relativity the Boltzmann equation is given by Γ j µν are the Christoffel symbols and so the left hand side expresses the fact that particles undergo geodesic motion in between point collisions.
The term C[f ] on the right hand side of the Boltzmann equation is called the collision operator and models the short range scattering processes that cause deviations from 3 geodesic motion. For 2 ↔ 2 reactions between fermions, such as neutrinos and e ± , the collision operator takes the form Here |M| 2 is the process amplitude or matrix element, S is a numerical factor that incorporates symmetries and prevents over-counting, f i are the fermi blocking factors, δ(∆p) enforces four-momentum conservation in the reactions, and the δ 0 (p 2 i −m 2 i ) restrict the four momenta to the future timelike mass shells.
We now restrict our attention to systems of fermions under the assumption of homogeneity and isotropy. We assume that the particles are effectively massless, i.e. the temperature is much greater than the mass scale. Homogeneity and isotropy imply that the distribution function of each particle species under consideration has the form f = f (t, p) where p is the magnitude of the spacial component of the four momentum. In a spatially flat FRW universe the Boltzmann equation reduces to This, combined with the Einstein equations and the matrix elements for the relevant processes, constitutes the dynamical equations governing neutrino freeze-out.
We also obtain formulas for the rate of change in the number density and energy density For free-streaming particles the vanishing of the collision operator implies conservation of 'comoving' particle number of species 1. From the associated powers of a in Eq. (7) and Eq. (8) we see that the energy of a free streaming particle scales as 1/a. This means that the distribution of a free streaming massive particle species will, once the mass scale becomes relevant, evolve into non-thermal shape [7,21]. The matrix elements for weak force scattering processes involving neutrinos and e ± are given in tables 1 and 2. They were obtained from Ref. [17] and are valid in the limit |p| ≪ M W , M Z , where in vacuum the gauge boson masses are M W = 80.4 GeV, M Z = 91.19 GeV.

Weinberg Angle
The SU (2) × U (1) gauge coupling constants g, g ′ are constrained by the two physical parameters, the Weinberg angle θ W and the electric charge e 4 Process S|M| 2 ν e +ν e → ν e +ν e 128G 2  [17] for electron neutrino processes where particle j = µ, τ ,  [17] for µ and τ neutrino processes where i = µ, τ , j = e, µ, τ , j = i, An alternative way to write these constraints is θ W enters the matrix elements presented in tables 1 and 2 by way of The Fermi constant G F fixes the vacuum expectation value of the Higgs field The mass of the W and Z gauge bosons can be written in terms of v We show the dependence of M W and M Z on Weinberg angle in figure 1, using ve/2 = 38.4 GeV at the Z-scale.
Fixing v and e, M W is minimized when sin θ W = 1. Thus we find a factor of 1/2 as a maximal possible reduction in M W , also seen in figure 1. This implies that M Z > 5 M W ≫ |p| for neutrino momentum p in the energy range of neutrino freeze-out, around 1 MeV, even as we vary sin θ W . Even if v is allowed to vary, for this approximation to cease to be valid it would have to be reduced by a factor of 10 5 , in our view an extreme amount. Therefore we can carry out the computation neutrino decoupling within the effective Fermi theory of weak interactions. The ratio M W /M Z = cos θ W = 0.8815 implies sin 2 θ W = 0.223. Considering that there is a rapid change with scale the actual values are sin 2 θ W MZ = 0.23116 ± 0.00012 and sin 2 θ W MW = 0.22296 ± 0.00028. We will present our results as a function of sin 2 θ W which we consider to be an unknown parameter in the hot Universe aged about one second. The other SM parameter of the electro-weak theory is the electric charge e. A variation in e is also possible, for example, due to time evolution of the grand unified scale [22].
The symmetry breaking parameter sin 2 θ W is at present a measured but theoretically unconstrained SM parameter. However, should a grand unified approach in which the strong interactions are merged into the electroweak interactions be discovered, then presumably sin 2 θ W could become fixed by the particular group structure. Such models are strongly constrained by proton decay limits [23], hence a fundamental constraint on sin 2 θ W is not (yet) in sight.

Interaction Strength Parameter η
In order to isolate the combination of natural constants which controls the neutrino freeze-out process, we cast the Einstein-Boltzmann model of neutrino freeze-out into dimensionless form. In the first step we look at the expansion of the Universe i.e. the 6 Hubble parameter H. The Einstein equations contain the Hubble equation where ρ = T 0 0 is the total gravitating energy density of the Universe and, as is usual in the context of general relativity, the Planck mass M p incorporates the factor 8π in the definition, Eq. (3).
The divergence freedom of the Einstein equations requires divergence freedom of the stress energy tensor T µν , a condition which reads for a homogeneous Universė Combining Eq. (14) with Eq. (15) shows that time change occurs at scale τ ∝ M p / √ ρ. In the domain of interest the energy density ρ is characterized by the electron mass ρ ∝ m 4 e . The scale m e is related both to the key energy component of the Universe at the time of neutrino freeze-out and the ambient temperature. We thus recognize the time scale to be characterized by τ ≡ M p /m 2 e = 6.12 s; the actual time scale is close to 1s considering the presence of many degrees of freedom.
Using the timescale τ , and scaling all momenta, energies, energy densities, pressures, and temperatures by the appropriate power of m e we can combine all scale dependent parameters in the Einstein-Boltzmann equation. We thus find where in the interaction strength η, Eq. (3), we include the G 2 F factor common to all of the neutrino interaction matrix elements.
Aside from the θ W dependence of the matrix elements seen in tables 1 and 2, the complete dependence on natural constants is now contained in a single dimensionless interaction strength parameter η with the vacuum present day value, If the dominant component of the electron mass originates in the Higgs mechanism, we find somewhat different scaling η ∝ g 3 Y e M p /v, where Yukawa electron coupling is introduced The discussion we presented is only focused on the normalization by natural constants of the collision term. The magnitude of the scattering integrals also depends on the magnitude of the scaled temperature T /m e . In particular, in the limit T /m e < 1 the scattering integrals involving e ± neutrino scattering are suppressed exponentially by a factor e −me/T or e −2me/T due to the diminished presence of e ± pairs. However, our objective in writing Eq. (16) was not to isolate the leading order behavior, but rather to separate out all dependence on dimensioned natural constants and isolated them in the interaction strength parameter η. This means that, as a dynamical system, the solutions of the dimensionless form Eq. (16) depends only on the parameters η and sin 2 θ W , and hence all quantities computed from solutions of the Boltzmann equation that are dimensionless, such as N ν , can also only depend on η and sin 2 θ W . Of course dimensioned quantities, for example the magnitude of freeze-out temperatures, still have to be scaled appropriately i.e. energies must be multiplied by m e and times must be multiplied by the timescale τ , and so dimensioned quantities will show an additional dependence on natural constants.
Our argument that there are only two dimensionless variables of interest, η, sin 2 θ W , relies on the fact that there is only one particle scale parameter that enters the energy density and collision integrals, namely m e . This is so since for T ∈ (0.1, 3) MeV muons are too heavy m µ = 105.66, and the baryon energy density controlled by M B − µ B is too small and all other energy components in Universe are completely negligible. Thus though in principle N ν (η, sin 2 θ W , m e /m µ , m e /(M B − µ B )), we can safely ignore all additional dimensionless quantities. Furthermore, given our hypothesis that a modification of SM parameters in the early Universe could contribute to N ν , there is also a contribution to the Universe dynamics from the rate of change of these parameters. We assume that any such rate of change is small enough to be insignificant and will not discuss it further.
The dependence of N ν (η, sin 2 θ W ) will be the key result of this work and is presented below in the section 4. Qualitatively, it is apparent that an increase in N ν requires increased coupling strength η. The dependence on sin 2 θ W is much less obvious in view of the gauge boson mass M W , M Z variation, see figure 1. The key question we aim to resolve in this work is how sensitive is N ν to a change in η and sin 2 θ W .

Emerging Chemical Nonequilibrium Method
We solve the Boltzmann equation Eq. (6) by the spectral method detailed in [24]. We give only a brief outline of the method here. Our approach is adapted to systems near kinetic equilibrium (i.e. equilibrium momentum distribution) but not necessarily chemical equilibrium (i.e. allowing for non-equilibrium particle number yield). In other words, the method performs best when the distribution is of the form where φ is small and T and Υ are the dynamical effective temperature and fugacity (i.e. phase space occupation parameter) respectively. Since we adapt both T and Υ as function of time, we employ a moving (in Hilbert space) frame, in which the orthogonal polynomial basis dynamically evolves to suit the problem. Our approach should be contrasted with the method used in [19,20], which we call the chemical equilibrium method, that studied neutrino freeze-out using a fixed orthogonal polynomial basis generated by the chemical equilibrium weight We note in the above that the temperature scaling is also assumed, that is T a(t) =Const. In our approach we allow for reheating of the effective temperature to occur and thus also T , like Υ, evolves in time independently. The deviation from chemical equilibrium is characterized in f Υ by the fugacity Υ. A non-equilibrium Υ = 1 builds up during neutrino freeze-out, specifically in the temperature range where the process e + e − → νν is too slow to equilibrate particle number 8 but e ± ν → e ± ν scattering is still able to equilibrate momentum. The introduction of chemical non-equilibrium through Υ = 1 contrasts to the spectral method employed in [19,20] which is adapted to chemical equilibrium. The chemical equilibrium method is appropriate for the physical regime studied in [19,20], wherein neutrinos are almost entirely decoupled by the time of e ± annihilation and therefore there is little time for reheating of neutrinos or the development of chemical nonequilibrium. However, for our purposes, namely the characterization of N ν (η, sin 2 θ W ), we must use a method that does not rely on small coupling for its effectiveness, hence we were motivated to develop the method described here. A detailed comparison with the method used by previous authors is presented in section 3.4. We refer to Ref. [24] for further discussion and detailed validation of the method we present.
After changing variables z = p/T , we will solve Eq. (6) by expanding φ in the basis of orthonormal polynomials,ψ i (z), generated by the parametrized weight function By convention, they are indexed so thatψ j has degree j. This choice of the weight is physically motivated by the phase space of practically massless neutrinos emerging into a chemical non-equilibrium distribution. The Boltzmann equation then results in an equation for the mode coefficients [24] b where the matrices A and B are For details on how to construct the inner products ∂ψi ∂Υ ,ψ k we refer to Appendix A of Ref. [24].
The dynamics of the effective temperature and fugacity are fixed by the requirement that f Υ captures the number density and energy density of the full distribution f , leaving φ to describe only the non-thermal distortions. In practice, this implies that b 0 (t) = b 1 (t) = 0 and a minimum of only two degrees of freedom (or modes), T and Υ, are required for our method. See Ref. [24] for details on the resulting evolution equations for T (t) and Υ(t). In contrast, we note that the minimum number modes required for the chemical equilibrium method is four.

Collision Integral Inner Products
In order to solve for the mode coefficients, the inner products of collision integrals with respect to the weight function Eq. (20), must be computed.
The matrix element for a 2 − 2 reaction can be written as a function of the Mandelstam variables s, t, u, of which only two are independent, defined by and we will consider that done for the analysis that follows. Note that R k only uses information about the distributions at a single spacetime point, and so we can work in a local orthonormal basis for the momentum. Among other things, this implies that p 2 = p α p β η αβ where η is the Minkowski metric From Eq. (26), we see that a crucial prerequisite of our spectral method is the capability to evaluate integrals of the type for some functions g i . Even after eliminating the delta functions in Eq. (36), we are still left with an 8-dimensional integral. To facilitate numerical computation, we analytically reduce this expression down to fewer dimensions. Fortunately, the systems we are interested in have a large amount of symmetry that we can utilize for this purpose.
The distribution functions we are concerned with are isotropic in some frame defined by a unit timelike vector U , i.e. they depend on the four-momentum p i only through the quantities p i · U and p 2 i = m 2 i . The same is true of the basis functionsψ k and hence we can assume the g i depend only on p i · U as well. In [16,17] approaches are outlined that reduce integrals of this type down to 3 dimensions. However, the integrand one obtains from these methods is only piecewise smooth or has an integration domain with a complicated geometry. This can present difficulties for numerical integration routines and so we take an alternative approach that, for the scattering kernels found in e ± , neutrino interactions, reduces the problem to three iterated integrals (but not quite to a three dimensional integral) and results in an integrand with better smoothness properties. Depending on the integration method used, this can significantly reduce the numerical cost of evaluating the collision integrals. The derivation presented expands on what is found in Ref. [25].

Simplifying the Collision Integral
Our strategy for simplifying the collision integrals is as follows. We first make a change of variables designed to put the 4-momentum conserving delta function in a particularly simple form, which allows us to analytically use that delta function to reduce the integral from 16 to 12 dimensions. The remaining four delta functions, which impose the mass shell constraints, are then seen to reduce to integration over a product of spheres. The simple form of the submanifold that these delta function restrict the integration to allows the method described in Appendix A to analytically evaluate all four of the remaining delta functions simultaneously. During this process, the isotropy of the system in the frame given by the 4-vector U allows us to reduce the dimensionality further, by analytically evaluating several of the angular integrals.
The change of variables that simplifies the 4-momentum conserving delta function is given by The Jacobian of this transformation is 1/2 8 . Therefore, changing variables in the delta functions we find where Θ(x) denotes the Heaviside function, b = 1/256(2π) 8 , and U is the four velocity characterizing the isotropic frame as discussed above.
Using the coarea formula, theorem 2 in Appendix A, we decompose this into an integral over s = p 2 , the center of mass energy, and also eliminate the integration over The lower bound on s comes from the fact that both p 1 and p 2 are future timelike and hence The other inequality is obtained using p = p ′ . Note that the integral in brackets in Eq. (40) is invariant under SO(3) rotations of p in the frame defined by U . Therefore we obtain where | p| denotes the norm of the spacial component of p and in the formula for K(s, p), p is any four vector whose spacial component has norm | p| and timelike component | p| 2 + s. Note that in integrating over δ(p 2 − s)dp 0 , only the positive root was taken, due to the Heaviside functions in the K(s, p).
We now simplify K(s, p) for fixed but arbitrary p and s that satisfy p 0 = | p| 2 + s and s > s 0 . These conditions imply p is future timelike, hence we can we can change variables in q, q ′ by an element of SO(1, 3) so that Note that the delta functions in the integrand imply p ± q is timelike (or null if the corresponding mass is zero). Therefore p 0 > ±q 0 iff p ∓ q is future timelike (or null). This condition is preserved by SO(1, 3), hence p 0 > |q 0 | in one frame iff it holds in every frame. Similar comments apply to p 0 > |(q ′ ) 0 | and so K(s, p) has the same formula in the transformed frame as well.
We now evaluate the measure that is induced by the delta functions, using the method given in Appendix A. We have the constraint function and must compute the solution set Φ(q, q ′ ) = 0. Adding and subtracting the first two components and the last two respectively, we have the equivalent conditions If we let (q 0 , q), ((q ′ ) 0 , q ′ ) denote the spacial components in the frame defined by p = ( √ s, 0, 0, 0) we have the equivalent conditions Note that the above formulas, together with s ≥ s 0 , imply and similarly for q ′ . Hence the Heaviside functions are identically equal to 1 under these conditions and we can drop them from the formula for K(s, p). The conditions Eq. (48) imply that our solution set is a product of spheres in q and q ′ , as long as the conditions are consistent i.e. so long as | q|, | q ′ | > 0. To see that this holds for almost every s, first note Therefore, for s > s 0 we have | q| > 0 and similarly for q ′ . Hence we have the result where B r denotes the radius r ball centered at 0. We will parametrize this by spherical angular coordinates in q and q ′ . We now compute the induced volume form. First consider the differential Evaluating this on the coordinate vector fields ∂ q 0 , ∂ r we obtain 13 Similar results hold for q ′ . Therefore we have the determinant By corollary 1 and Eq. (A.18) in Appendix A, this implies that the induced volume measure is and i is the interior product (i.e. contraction) operator as described in Appendix A. Consistent with our interest in the Boltzmann equation, we assume F factors as For now, we suppress the dependence on p, as it is not of immediate concern. In our chosen coordinates where U = (α, 0, 0, δ) we have and similarly for q ′ .
To compute Together, these imply that the integral in brackets in Eq. (61) equals Therefore This is as far as we can simplify things without more information about the form of the matrix elements. In Appendix B we apply this method and analytically simplify Eq. (66) for each of the processes in tables 1 and 2 as much as possible and in the process we show that M can be written in terms of three iterated integrals for each of these processes.

Validation
We solve the Boltzmann equation, Eq. (4), for both the electron neutrino distribution and the combined µ, τ neutrino distribution, including all of the processes from tables 1 and 2 in the scattering operator, together with the Hubble equation for a(t), Eq. (14). The total energy density appearing in the Hubble equation consists of the contributions from both neutrino distributions as well as chemical equilibrium e ± and photon distributions at some common temperature T γ . The dynamics of T γ are fixed by the divergence freedom of the total stress energy tensor, Eq. (15). In addition, we include the QED corrections to the e ± and photon equations of state as described in [20].
We compared the results of numerically evaluating the collision integrals using our method as given in sections 3.2 and Appendix B with the method used by Ref. [17] and validated that results agree, up to numerical integration tolerance. To compare our results from solving the Boltzmann equation with Ref. [3], where neutrino freeze-out was simulated using sin 2 θ W = 0.23 and η = η 0 , in table 3 we present N ν together with the following quantities  Table 3: Comparison of neutrino freeze-out results obtained in Ref. [3] (top line) with those obtained using the methods described above, which allow for a reduced number of expansion modes.
The quantities presented in Eq. (69) and table 3 were introduced in Ref. [3], but some additional discussion is in order.
1. The quantity z fin measures the deviation of the photon temperature T γ from the 'free streaming' temperature T ∝ 1/a, i.e the temperature of a (hypothetical) particle species that is completely decoupled throughout the domain of temperature considered. Therefore, z fin = T γ /T is the measure of the amount of reheating photons underwent due to the annihilation of e ± . For the case of already completely decoupled neutrinos, whose temperature is in this case just the free-streaming temperature, according to Eq. (2) For the case of some e ± annihilation occurring while neutrinos are still coupled, one expects this value to be slightly reduced, due to the transfer of some e ± entropy into neutrinos. This is reflected in the values seen in table 3. 2. ρ ν0 is the energy density of a single massless fermion with two degrees of freedom and temperature equal to the free-streaming temperature. In other words, it is the energy density of a single neutrino species, assuming it decoupled before reheating. Consequently, δρ ν is the fractional increase in the energy density of a coupled neutrino species, due to its partial participation in reheating.
The top entry in table 3 correspond to the reference values from Ref. [3]. The next two lines present our results that use the chemical non-equilibrium method which, as we show, allows for a smaller basis set. We show 2 and 3 modes, respectively which compare with the 4 modes case for the equilibrium method. The value of N ν we obtain agrees for the case of 2 modes with that found by [3], up to their cited error tolerance.
Considering both the improved smoothness properties of integrands developed in sections 3.2 and Appendix B and the reduction in the required number of modes for the chemical non-equilibrium method, our approach with the minimum number 2 of required modes is found to be more than 20× faster than the chemical equilibrium method with its minimum number 4 of required modes. This computational performance improvement makes it possible to explore the neutrino freeze-out process for many different circumstances and parameter sets employing a desktop PC.

Neutrino Freeze-out Temperature
SM parameters impact N ν by changing how long the neutrinos remain coupled to the annihilating e ± and thereby impacting the amount of energy and entropy transfer. In other words, the neutrino freeze-out temperature is modified. Before we present the dependence of N ν on θ W and η we first consider in detail the freeze-out temperatures of neutrinos with an initial focus on the conventional SM parameters.
In the literature one finds estimates of freeze-out temperatures based on a comparison of Hubble expansion with neutrino scattering length and considering only number changing (i.e. chemical) processes, see e.g. Ref. [14]. We employ a similar definition of freeze-out temperature in the context of the Boltzmann equation and refine the results by noting that there are three different freeze-out processes: 1. Neutrino chemical freeze-out: the neutrino pair number changing annihilation process l +l ⇔ ν l +ν l (71) which we will see decouples at the highest temperature. 2. Neutrino kinetic freeze-out: The sharing of energy between leptons and neutrinos by way of scattering l + ν ⇔ l + ν stops at a lower energy compared to neutrino number changing processes. 3. Collisions between neutrinos are capable of re-equilibrating energy within and between flavor families. These processes end at a yet lower temperature and the neutrinos will be truly free-streaming from that point on.
The attentive reader will notice that we have omitted here a discussion of flavor neutrino oscillations. If it weren't for the differences between the matrix elements for the interactions between e ± and ν e on one hand and e ± and ν µ , ν τ on the other, oscillations would have no effect on the flow of entropy into neutrinos and hence no effect on N ν . However, there are differences and they do lead to a modification of N ν . In Ref. [3] the impact of oscillations on neutrino freeze-out for the present day measured values of θ W and η was investigated. It was found that while oscillations redistributed energy amongst the neutrino flavors, the impact on N ν was negligible. We have therefore ignored neutrino oscillation effects in our study as we do not have a clear idea why for other values of η and θ W the redistribution of neutrino energy would create any larger effect than already determined. Once the relevant neutrino properties are fully understood, the precision of the results we present could possibly be improved by incorporating the effect of neutrino oscillations.

Scattering Length and Freeze-out Temperature
The notion of freeze-out temperature is conceptually useful, but within the Boltzmann approach there is no such precise temperature, as the freeze-out process is gradual, with low energy neutrinos freezing-out before high energy ones. Thus a procedure to determine the freeze-out condition can only be approximate. However, the differences that arise through investigating the freeze-out of the three different classes of processes while natural constants are varied help us to understand the results which will be presented below.
To define the freeze-out condition we follow the standard procedure [14]: we compare the distance L traveled by a particle between two scattering processes to the characteristic Universe expansion length L H = c/H. The crossing of the Hubble-length with the neutrino scattering length produces an estimate of the decoupling temperatures. To obtain the scattering length L we begin with Eq. (7) for the fractional rate of change of comoving particle number d dt (a 3 n) This expression includes both forward and back-reactions, producing the net change. However, we would rather count the number of interactions. For that reason, we consider only one direction for the process and define as the rate of interest where the forward-reaction operatorC[f ] is computed as in Eq. (5) except with F replaced byF If particle type 1 also participates in the reverse of the reaction 1 + 2 → 3 + 4 (i.e. it is the same as the final particle 3, 4) then a corresponding term for the reverse reaction must also be added. The key point is that we are counting reactions, and not the net particle number change which requires by detailed balance also a negative contribution. Using the average velocity, which for neutrinos isv = c = 1, we obtain according to this procedure the scattering length where the sum is over the one way scattering operators for the collection of processes of interest. Like in Ref. [14], L can be compared to the Hubble length L H = c/H and the temperature at which L = L H we call the freeze-out temperature for that reaction. Figure  2 shows L H together with the scattering length for the three types of neutrino reactions described above, on the left for ν e and the right for ν µ , ν τ . The flavor dependence is due to charge current W-mediated interactions being present only for ν e . The solid lines in ν − e ± Chemical Freeze-out ν − e ± Kinetic Freeze-out ν − ν Freeze-out ν − e ± Chemical Freeze-out ν − e ± Kinetic Freeze-out ν − ν Kinetic Freeze-out ν − e ± Chemical Freeze-out ν − e ± Kinetic Freeze-out ν − ν Kinetic Freeze-out  We see in figure 3 that as a function of θ W the behavior of T νe is opposite to that of T νµ and T ντ for sin 2 θ W < 0.25, and for sin 2 θ W > 0.25 all neutrino processes tend to decouple at lower temperature with increasing sin 2 θ W . Neutrino-neutrino scattering process remains constant, as their matrix elements are independent of Weinberg angle. An increased coupling strength η is equivalent to an increase in G F , resulting in the neutrinos interacting with the e ± plasma down to lower temperatures. Hence the monotonic decreasing behavior of the freeze-out temperature as a function of η/η 0 seen in figure 3 is expected. The SM values are seen at the left margin of the bottom panels. As discussed above, neutrino oscillations are not considered in these results. We expect that incorporating oscillations would lead to a smaller difference between the freeze-out temperatures of the different neutrino flavors, and would likely pull up the drop in ν µ , ν e freeze-out temperature at small sin 2 θ W , at least to some degree. We recall that for other observable quantities we discuss in the following, the effect of neutrino oscillations is expected to be negligible [3].

Dependence of N ν on Standard Model Parameters
The main result of this paper is the dependence of N ν on the SM parameters sin 2 θ W and η, Eq. (3). These results are shown in figure 4, presented as a function of Weinberg angle sin 2 θ W for η/η 0 = 1, 2, 5, 10. The effect of an increase in both parameters above the vacuum values superpose in the parameter range considered, amplifying the effect and generating a significant increase in N ν → 3.5. The present day vacuum value of Weinberg angle puts the ν µ , ν τ freeze-out temperature, seen in the right pane of figure 3, near its maximum value. This is why a comparatively large change in sin 2 θ W is needed to produce a change in N ν for sin 2 θ W ≈ 0.23.
We performed a least squares fit of N ν over the range 0 ≤ sin 2 θ W ≤ 1, 1 ≤ η/η 0 ≤ 10 shown in figure 4, obtaining a result with relative error less than 0.2%,  N ν is monotonically increasing in η/η 0 with dominant behavior scaling as η/η 0 . Monotonicity is to be expected, as increasing η decreases the freeze-out temperature and the longer neutrinos are able to remain coupled to e ± , the more energy and entropy from annihilation is transferred to neutrinos.
The bounds on N ν from the Planck analysis [1] can be used to constrain time or temperature variation of sin 2 θ W and η. In Figure 5 the dark (green) color shows the combined range of variation of natural constants compatible with CMB+BAO and the light (teal) color shows the extension in the range of variation of natural constants for CMB+H 0 , both at a 68% confidence level. The dot-dashed line within the dark (green) color delimits this latter domain. The dotted line shows the limit of a 5% change in N ν . Any increase in η/η 0 and/or sin 2 θ W moves the value of N ν into the domain favored by current experimental results.
Further parameter study is found in Appendix C. In the figure C.6 and the data fits Eqs.(C.1-C.4) we complement the N ν results by showing the variation of the parameters that characterize the neutrino distributions after freeze-out: the neutrino temperature, shown through the ratio of the reference photon to neutrino temperature T γ /T ν separately for ν e and ν µ , ν τ and well as the two fugacities Υ νe and Υ νµ = Υ ντ .

Summary, Discussion and Conclusions
We have employed a novel spectral method Boltzmann solver and a new procedure for evaluating the Boltzmann scattering integrals in order to characterize the impact of a potential time and/or temperature variation of SM parameters on the effective number of neutrinos. Specifically, we identified a dimensionless combination of m e , M p , and G F , called the interaction strength η, that, along with the Weinberg angle sin 2 θ W , control neutrino freeze-out and the resulting value of the effective number of neutrinos, N ν .

Novel Mathematical Tool
In order to carry this comprehensive study we have developed a novel approach to obtain Boltzmann Equation solutions. Our spectral method, which we call the emergent chemical non-equilibrium method, employs a moving (in Hilbert space) frame, in which the orthogonal polynomial basis dynamically evolves to suit the problem. Our approach as presented here makes several modifications that both improve its numerical speed and make it better suited to the regime we are investigating, namely the stronger coupling between neutrinos and e ± that is obtained when SM parameters are varied, and that lead to an increase in N ν . As detailed in the general presentation of the method [24], the improvements are 1. We allow a general time dependence of the effective temperature parameter T i.e.
we do not assume redshift temperature scaling T ∝ 1/a -this accommodates the effect of reheating. 2. We have introduced a chemical non-equilibrium distribution in the weight function i.e. we introduced an evolving, time dependent Υ which equals 1 at high temperature, corresponding to chemical equilibrium, and allows for the emergence of chemical non-equilibrium Υ = 1 during freeze-out. 3. We have introduced an additional factor of z 2 to the functional form of the weight as proposed in a different context in Refs. [26,27] which accounts in our approach for the effectively massless neutrino phase space.
The salient feature is that we are letting the fugacity, Υ and temperature T be time dependent and there is no requirement that Υ → 1. This should be contrasted with the method used in [19,20], which we call the chemical equilibrium method, that studied neutrino freeze-out using a fixed orthogonal polynomial basis generated by the chemical equilibrium weight and without the z 2 factor. The chemical equilibrium method also assumes a particular temperature scaling T a(t) =Const. In other words, the neutrino momenta are scaled by 1/a(t) instead of a dynamical effective temperature T (t) as in our method.
Due to the inclusion of the neutrino phase space z 2 factor in the weight, and facilitated by the near thermal shape of the distribution, only two modes corresponding to T and Υ are required to capture the energy density and number density of the neutrino distribution. In comparison, the chemical equilibrium method, because it lacks the z 2 factor, requires a minimum of four modes. We discussed how further important savings in computation time are arrived at by making the integrands of the collision integrals smooth functions. Overall, the speed up of solutions is at level 20 times or more.

Primordial Variation of Natural Constants
The question which we answer in this paper is: What neutrino decoupling in the early Universe can tell us about the values of natural constants when the Universe was about one second old and at an ambient temperature near to 1 MeV (11.6 billion degrees K). Our results were presented assuming that the Universe contains no other effectively massless particles but the three left handed neutrinos and corresponding, three right handed anti-neutrinos.
We found that near to the physical value of the Weinberg angle sin 2 θ W ≃ 0.23 the effect of changing sin 2 θ W on the decoupling of neutrinos is small. Thus as seen in Figure  4 the dominant variance is due to the change in the coupling strength η/η 0 , Eq. (3) and Eq. (17). The dotted line in Figure 5 shows that in order to achieve a change in N ν at the level of up to 5%, that is N ν 3.2, both sin 2 θ W and η/η 0 must change significantly, with e.g. η increasing by an order of magnitude.
Let us look closer at what an increase in the strength parameter η by factor 10 means, looking case by case on all the natural constant contributions as if each were responsible for the entire change: • Considering that η ∝ M p ∝ G −1/2 N this translates into a decrease in the strength of Gravity at neutrino freeze-out by a factor 100. This effect would need to become much smaller by the time the age of the Universe is 1000 times longer (1s compared to 10 min) for Big Bang nucleosynthesis to be unaffected. This presumably means that, conversely, as we go further back in time we would need gravity to continue to rapidly become very much weaker yet. In models of emergent gravity we can imagine a 'melting' of gravity in the hot primordial Universe. Whether such a model can be realized will be a topic for future consideration. The attractive aspect of Gravity weakening rapidly with increasing temperature is that for exponentially disappearing G N → 0 as t → 0 and/or T → ∞ the dynamics can be arranged to be similar to an inflationary Universe.
• Since η ∝ m 3 e , electron mass would need to go up 'only' by factor 2.15 . Compared to all other particles the electron mass has an anomalously low value. Appearance of a mechanism just when T ≃ m e that 'restores' the electron mass to where intuition would like it to be, a few MeV, arising from the systematics of other Yukawa Higgs coupling g Y e compared to the Yukawa coupling of other charged light particles, where m e = g Y e v seems to us also a possible scenario. Interestingly, laboratory limits for these conditions could be attainable in the foreseeable future.
• Since η ∝ G 2 F ∝ 1/v 4 we would need to find a mechanism that would decrease the vacuum value v 0 ≃ 246 GeV by factor 1.8 already at temperature T ≃ m e . Allowing three powers of v to cancel by using the Higgs minimal coupling formula for electron mass we need to change v by an order of magnitude near to T ≃ m e . This appears impossible.
While ideas justifying strong variation of η can be developed as two of the above three cases argue, a model for temperature or time dependence of sin 2 θ W seems at this time without a theoretical anchor point, mainly so since we do not have a valid grand unified theoretical framework in which the electro-weak mixing or equivalently the masses M W , M Z would be anchored.

Two Different Ways to Change N ν
There are additional challenges we have not at this time addressed. This is so since the immediate observable is the energy content of the invisible Universe as defined by the effective number of neutrinos N ν , see Eq. (1). Considering the value of N ν > 3 there could 23 be a contribution from, presently not discovered, weakly interacting massless particles that decoupled even before neutrinos, and which therefore could contribute fractionally to N ν , see our discussion in Ref. [9]. Of particular relevance could be a so called light sterile neutrino [28], possibly the right handed complement to the left handed neutrinos. If such particles exist and freezeout well before regular neutrinos, their contribution to N ν would be subject to dilution by reheating [9] and thus would depend on when precisely they begin free-streaming.
These unknown dark 'radiation' particles, as well as neutrinos, could have a mass that is at the scale of the temperature of photon decoupling T γ0 = 0.25 eV, for which an analysis of the Universe density fluctuations akin to Planck [1] would need to be adapted. We have discussed in Ref. [7] a consistent treatment of neutrino mass and N ν in the case of a particular type of delayed massive neutrino freeze-out. This approach is the same as for dark radiation: Near to T γ0 = 0.25 eV massive neutrinos are indistinguishable from massive dark radiation, which contributes as an additional particle with reduced contribution to N ν [9].
The alternative explanation of N ν > 3 in terms of variation of of natural constants that we have presented comprises speculative beyond the standard model ideas akin in this aspect to new dark 'radiation' particles. We believe that our present contribution provides a viable alternative mechanism capable of influencing N ν . In order to achieve an increase in N ν the change in natural constants must cause a delay in neutrino freeze-out and thus a greater participation of neutrinos in reheating during e ± annihilation. The changes in the natural constants which are required to make a large and visible change in N ν appear at first sight to reach beyond a variation that one could tolerate at the time of big bang nucleosynthesis, only a factor 1000 in time later. We have argued that a change in the electron mass m e by a factor larger than two, and/or Newton's constant G N even by several orders of magnitude can be considered.
where we omit the wedge product signs. In the following we will only be concerned with the results up to sign (i.e. with the density defined by a volume form). Iterated contractions with the vectors X j will be denoted by i (X1,...,Xm) . As we are only concerned with the result up to sign, the order in which we contract is irrelevant.
The Riemannian method of inducing volume measures extends to submanifolds of codimension greater than one as well as to semi-Riemannian manifolds, as long as the metric restricted to the submanifold is non-degenerate, by contracting with an orthonormal basis for the normal vectors. However, there are many situations where one would like to define a natural volume form on a submanifold that is induced by a volume form in the ambient space, but where the above method is inapplicable, such as defining a natural volume form on the light cone or other more complicated degenerate submanifolds in relativity. In this appendix, we will describe a method for inducing volume forms on regular level sets of a function that is applicable in cases where there is no metric structure and show its relation to more widely used semi-Riemannian case.
Let M , N be smooth manifolds, c be a regular value of a smooth function F : M → N , and Ω M and Ω N be volume forms on M and N respectively. Using this data, we will be able to induce a natural volume form on the level set F −1 (c). The absence of a metric on M is made up for by the additional information that the function F and volume form Ω N on N provide. The following theorem makes our definition precises and proves the existence and uniqueness of the induced volume form. Proof. F * is onto T F (x) N for any x ∈ F −1 (c). Hence there exists {v i } n 1 ⊂ T x M such that Ω N (F * v 1 , ..., F * v n ) = 1. (A.5) In particular, F * v i is a basis for T F (x) N . Define ω x = i (v1,...,vn) Ω x . This is obviously a nonzero m − n form on T x F −1 (c) for each x ∈ F −1 (c). We must show that this definition is independent of the choice of v i and the result is smooth.
since the terms involving ker F * will vanish on T x F −1 (c) = ker F * x . Therefore This proves that ω is independent of the choice of v i . If we can show ω is smooth then we are done. We will do better than this by proving that for any v i ∈ T x M the following holds To see this, take w i satisfying Eq. (A.4). Then F * v i = A j i F * w j . This determinant can be computed from Therefore, the same computation as Eq. (A.7) gives as desired. To prove that ω is smooth, take a smooth basis of vector fields {V i } m 1 in a neighborhood of x. After relabeling, we can assume {F * V i } n 1 are linearly independent at F (x) and hence, by continuity, they are linearly independent at F (y) for all y in some neighborhood of x. In that neighborhood, Ω N (F * V 1 , ..., F * V n ) is non-vanishing and therefore ω = (Ω N (F * V 1 , ..., F * V n )) −1 i (V1,...,Vn) Ω (A.14) which is smooth.
Corollary 2. If φ : M → R is smooth and c is a regular value then by equipping R with its canonical volume form we have where v ∈ T x M is any vector satisfying dφ(v) = 1.
A coarea formula can be proved for the induced volume forms.
Theorem 2 (Coarea formula). Let M be a smooth manifold with volume form Ω M , N a smooth manifold with volume form Ω N and F : M → N be a smooth map. If F * is surjective at a.e. x ∈ M then for f ∈ L 1 (Ω M ) L + (M ) where ω M z is the volume form induced on F −1 (z) as in theorem 1.
The induced measure defined above allows for a coordinate independent definition of a delta function supported on a regular level set. Such an object is of great use in performing calculations in relativistic phase space in a coordinate independent manner. Definition 1. Motivated by the coarea formula, we define the composition of the Dirac delta function supported on c ∈ N with a smooth map F : M → N such that c is a regular value of F by δ c (F (x))Ω M ≡ ω M (A.18) in place of It is useful to translate the induced volume element into a form that is more readily applicable to computations in coordinates. Choose arbitrary coordinates y i on N and write Ω N = h N (y)dy n . Choose coordinates and write Ω M = h M (x)dx m . The coordinate vector fields ∂ x i are transverse to F −1 (c) and so Therefore we obtain Using Eq. (A.24), along with the coordinates described there, we can (at least locally) write the integral with respect to the delta function in the more readily usable form The absolute value comes from the fact that we use δ c (F (x))Ω M to define the orientation on F −1 (c). where the coefficient C is given in table B.4. From here we obtain Therefore, as we claimed above, M νν→νν can be written in a form that requires the numerical evaluation of only three iterated integrals, but not quite as a three dimensional integral. If we want to emphasize the role of C then we write M νν→νν (C). Note that if one scales p and s by the appropriate powers of T in order to convert to dimensionless variables, one obtains a prefactor of T 8 .

Appendix B.2. νν → νν
Using Eq. (31), the matrix elements for neutrino anti-neutrino scattering can be simplified to where the coefficient C is given in table B.5. Using this we find Again, by converting to dimensionless variables we see that this scales with T 8 . If we want to emphasize the role of C then we write M νν→νν (C). Note that due to the polynomial form of the matrix element integral, the double integral in brackets breaks into a linear combination of products of one dimensional integrals, meaning that the nesting of integrals is only three deep.
Appendix B.3. νν → e + e − Using Eq. (31), the matrix elements for neutrino anti-neutrino annihilation into e ± can be simplified to where the coefficients A, B, C are given in table B.6.
Process The integral of each of these terms is By scaling s, p, and m e by the appropriate powers of T we again obtain a prefactor of T 8 . If we want to emphasize the role of A, B, C then we write M νν→e + e − (A, B, C). Note that this expression is linear in (A, B, C) ∈ R 3 . Also note that, under the assumptions that the distributions of e + and e − are the same (i.e. ignoring the small matter antimatter asymmetry), the G ij terms that contain the product of e ± distributions are even functions. Hence the term involving the integral of yz vanishes by antisymmetry.
Appendix B.4. νe ± → νe ± Using Eq. (31), the matrix elements for neutrino e ± scattering can be simplified to where the coefficients A, B, C are given in table B.7.
Process The integral of each of these terms is As above, after scaling s, p, and m e by the appropriate powers of T we obtain a prefactor of T 8 . If we want to emphasize the role of A, B, C then we write M νe→νe (A, B, C). Note that this expression is also linear in (A, B, C) ∈ R 3 .

Appendix B.5. Total Collision Integral
We now give the total collision integrals for neutrinos. In the following, we indicate which distributions are used in each of the four types of scattering integrals discussed above by using the appropriate subscripts. For example, to compute M νeνµ→νeνµ we set G 1,2 =ψ j f 1 f 2 , G 3,4 = f 3 f 4 , f 1 =ψ j f νe , f 3 = f νe , and f 2 = f 4 = fν µ in the expression for M νν→νν from section Appendix B.3 and then, to include the reverse direction of the process, we must subtract the analogous expression whose only difference is G 1,2 =ψ j f 1 f 2 , G 3,4 = f 3 f 4 . With this notation the collision integral for ν e is Symmetry among the interactions implies that the distributions of ν µ and ν τ are equal. We also neglect the extremely small matter anti-matter asymmetry and so we take the distribution of each particle to be equal to that of the corresponding antiparticle. Therefore there are only three independent distributions, f νe , f νµ , and f e and so we can combine some of the terms in Eq. (B.20) to obtain Introducing one more piece of notation, we use a subscript k to denote the orthogonal polynomial basis element that multiplies f 1 or f 1 in the inner product. The inner product of the kth basis element with the total scattering operator for electron neutrinos is therefore Under these same assumptions and conventions, the total collision integral for the combined ν µ , ν τ distribution (which we label ν µ ) is Appendix B.6. Conservation Laws and Scattering Integrals For some processes, some of the R k 's vanish exactly. As we now show, this is an expression of various conservation laws. First consider processes in which f 1 = f 3 and f 2 = f 4 , such as e ± ν → e ± ν. Since m 1 = m 3 and m 2 = m 4 we have r = r ′ , q 0 = (q ′ ) 0 . The scattering terms are all two dimensional integrals of some function of s and p multiplied by 2π 0 S|M| 2 (s, t(cos(ψ) 1 − y 2 1 − z 2 + yz))dψ f1(h1(y))f2(h2(y))dy × f 1 k (h1(z))f 2 (h2(z))dz (B.25) 2π 0 S|M| 2 (s, t(cos(ψ) 1 − y 2 1 − z 2 + yz))dψ f 1 (h1(y))f 2 (h2(y))dy × f 1,k (h1(z))f2(h2(z))dz (B.26) h1(y) =(p 0 + (q ′ ) 0 α − r ′ δy)/2, h2(y) = (p 0 − q 0 α + rδy)/2, f 1,k =ψ k f1, f 1 k =ψ k f 1 .

(B.27)
For k = 0,ψ 0 is constant. After factoring it out of I k , the result is obviously zero and so R 0 = 0. We further specialize to a distribution scattering from itself i.e. f 1 = f 2 = f 3 = f 4 . Since m 1 = m 2 and m 3 = m 4 we have q 0 = (q ′ ) 0 = 0 and h 1 (y) = (p 0 − r ′ δy)/2, h 2 (y) = (p 0 + rδy)/2. (B.28) By the above, we know that R 0 = 0.ψ 1 appears in I 1 in the formψ 1 (h 1 (z)), a degree one polynomial in z. Therefore R 1 is a sum of two terms, one which comes from the degree zero part and one from the degree one part. The former is zero, again by the above reasoning. Therefore, to show that R 1 = 0 we need only show I 1 = 0, except witĥ ψ 1 (h 1 (z)) replaced by z. Since h 1 (−y) = h 2 (y), changing variables y → −y and z → −z in the following shows that this term is equal to its own negative, and hence is zero S|M| 2 (s, t(cos(ψ) 1 − y 2 1 − z 2 + yz))dψ f 1 (h 1 (y))f 1 (h 2 (y))dy × zf 1 (h 1 (z))f 1 (h 2 (z))dz. (B.30) We note that the corresponding scattering integrals do not vanish for the chemical equilibrium spectral method employed in [19,20]. This is another advantage of the method outlined in section 3.1. Further differences are discussed in section 5.1. Finally, we point out how the vanishing of these inner products is a reflection of certain conservation laws. From Eq. (7), Eq. (26), and the fact thatψ 0 ,ψ 1 span the space of polynomials of degree ≤ 1, we have the following expressions for the change in number density and energy density of a massless particle for some c 0 , d 0 , d 1 . Therefore, the vanishing of R 0 is equivalent to conservation of comoving particle number. The vanishing of R 0 and R 1 implies ρ ∝ 1/a 4 i.e. that the reduction in energy density is due entirely to redshift; energy is not lost from the distribution due to scattering. These findings match the situations above where we found one or both of R 0 = 0, R 1 = 0. R 0 vanishes for scattering processes that exchange momentum but don't change particle number. Both R 0 and R 1 vanished for a distribution scattering from itself and in such a process one expects that no energy is lost from the distribution by scattering, it is only redistributed among the particles corresponding to that distribution.