Neutrino non-standard interactions meet precision measurements of $N_{\rm eff}$

The number of relativistic species, $N_{\rm eff}$, has been precisely calculated in the standard model, and would be measured to the percent level by CMB-S4 in future. Neutral-current non-standard interactions would affect neutrino decoupling in the early Universe, thus modifying $N_{\rm eff}$. We parameterize those operators up to dimension-7 in the effective field theory framework, and then provide a complete, generic and analytical dictionary for the collision term integrals. From precision measurements of $N_{\rm eff}$, the most stringent constraint is obtained for the dimension-6 vector-type neutrino-electron operator, whose scale is constrained to be above about 195 (331) GeV from Planck (CMB-S4). We find our results complementary to other experiments like neutrino coherent scattering, neutrino oscillation, collider, and neutrino deep inelastic scattering experiments.


Introduction
The great triumph of the Standard Model (SM) of particle physics was the discovery of the Higgs particle in 2012 [1,2]. However, SM can not be the complete theory as there are still several unsolved puzzles, such as neutrino masses, dark matter and baryon asymmetry of the Universe, which require new physics beyond the SM. Tremendous new physics models have been invented and studied to address these issues, yet no definite signals of any of these models have been observed at colliders or from low-energy precision measurements. This has in turn motivated physicists to search for new physics in a model-independent and systematic way.
Effective Field Theories (EFTs) provide such a systematic and model-independent framework for the study of new physics, especially if its characteristic scale is above the weak scale. The EFT, obtained by integrating out the newly introduced heavy particles to vestigated [38], while the dimension-7 neutrino-electron operators are not yet considered to be constrained.
In the early Universe where only neutrinos, electrons, positrons, and photons are present, these neutrino-neutrino, neutrino-electron/positron and neutrino-photon NC NSIs would affect neutrino decoupling, thus modifying the effective number of relativistic degrees of freedom, viz., N eff . In light of the precision measurements of N eff from LEP [113] and Planck [114], the upcoming SPT-3G [115] and the Simon Observatory [116], as well as the proposal from Cosmic Microwave Background-Stage 4 (CMB-S4) [117], CORE [118], PICO [119] and CMB-HD [120], one naturally expects constraints on these NC NSIs from N eff .
In this work, we investigate all kinds of neutrino-neutrino, neutrino-electron/positron and neutrino-photon NC NSI operators up to dimension-7, as well as their impact on N eff . Since the light mediator directly serves as one additional degree of freedom and thus resulting in large N eff , these NC NSI operators we consider in this work are assumed to be induced by integrating out some heavy new physics above ∼ O(100 MeV) that is about the muon mass or heavier.
To obtain new physics corrections to N eff , the SM prediction to N eff has to be known precisely in the first place. However, it has been known for a long time that the precision calculation of N eff is very challenging. Within the SM, the precision calculation of N eff has been carried out through the density matrix formalism. Due to its complexity, however, the density matrix formalism is very difficult to generalize to other scenarios, for example, when new physics is present. For recent development of precision calculation of N eff , see Refs. [121][122][123][124][125].
In this work, we adopt the strategy developed in Refs. [124,125] that reproduces the SM prediction of N eff , works fast, and can be easily generalized to include effects from various new physics. To find corrections to N eff from some new physics, this strategy has already been applied in Refs. [126,127] with the introduction of right-handed partners of neutrinos, Ref. [128,129] with dark matter and/or sterile neutrinos, Ref. [130] for dark photon, Ref. [131] with the introduction of neutrino-scalar interactions, Ref. [132] with the inclusion of neutrino flavor oscillation and primordial nucleosynthesis, and Ref. [133] with a light Z to explain the recent XENON1T excess [134]. Applying this strategy to the calculation of N eff with the inclusion of NC NSIs up to dimension-7, in this work, we • provide a complete, generic and analytical dictionary for the collision term integrals in section 4. This dictionary can be used directly for computing corrections to N eff from some new physics, either in the EFT framework up to dimension-7, or in some UV models as long as the new physics is above ∼ O(100 MeV); • present our constraints on the NC neutrino NSI operators up to dimension-7 in section 5, and also compare our results with previous ones.
The rest of this work is organized as follows. We briefly review neutrino decoupling in the early Universe and the definition of N eff in section 2. In section 3, we discuss the strategy developed in Refs. [124,125], and then summarize our strategy for calculating the collision terms integrals. Since these collision term integrals are essential to boost the calculation of N eff , we provide a complete generic and analytical dictionary of the collision term integrals, as well as the NSI operators we study in this work in section 4. Constraints on these NC NSI operators are presented in section 5. We conclude in section 6.
2 Brief review of neutrino decoupling and N eff In the early Universe when the temperature is above O(10) MeV and below the muon mass, electrons, positrons, neutrinos and photons are in thermal equilibrium from electroweak interactions. As the Universe expands and the temperature cools down, neutrinos decouple from the rest of the plasma at around T dec = 2 MeV. The neutrinos then undergo simple dilution from the expansion of the Universe, while e ± and photons are still in thermal equilibrium. However, when the photon temperature cools further down below the electron mass m e , γγ → e + e − becomes suppressed while the inverse process is still permitted, heating up the photons.
The number of relativistic degrees of freedom during this period can be parameterized by N eff [135][136][137]: with ρ γ the photon energy density, and ρ R the total energy density from all relativistic species during this epoch. Equivalently, with ρ 0 ν the energy density of a single massless neutrino, and ρ 0 γ the energy density of photons in the instantaneous decoupling limit. Obviously, in the instantaneous limit, ρ 0 γ = ρ γ and ρ R = 3ρ 0 ν + ρ γ , resulting in the well-known N eff = 3. On the other hand, due to the tininess of neutrino masses, the three flavor neutrinos can be effectively taken as massless, permitting to express N eff in eq. (2.2) also in terms of the photon temperature T γ and the neutrino temperature T ν as, upon assuming T γ = T e which is valid since photons and electrons are tightly coupled during neutrino decoupling, Similarly, in the instantaneous decoupling limit, T ν /T γ = (4/11) 1/3 [138] and once again N eff = 3. However, it has been known for decades that the instantaneous decoupling picture is not accurate. Indeed, neutrinos are still slightly interacting with the electromagnetic plasma, and neutrino oscillations are also active during neutrino decoupling [139][140][141][142]. Furthermore, the electromagnetic plasma also receives corrections from finite temperature QED corrections. Taking all these effects into account, one finds N eff = 3.044 [123,132]. Corrections from these effects will be discussed further in detail in sections 3.1.1 and 4.3.

Setup of the Boltzmann equation
Evolution of phase space distribution (PSD) of any particle in the early Universe is governed by the Boltzmann equation, which we briefly review in this subsection. As mentioned in the introduction, we follow the discussion in Ref. [125], which simplifies the calculation of N eff significantly and reproduces the prediction for N eff by using the density matrix formalism.

The Boltzmann equation
The Boltzmann equation reads with f i (p, t) the PSD for particle i, H the Planck constant that accounts for the dilution effect from the expansion of the Universe, and C the collision term defined as 2 and "+ (−)" is for bosonic (fermionic) particles. Note that the difference in the last line above correctly accounts for the production and annihilation of particle i. Upon integrating over the phase space of particle i on both sides of eq. (3.1), one finds 3 where g is the intrinsic degree of freedom of particle i, E is its energy, and n and ρ are the number and the energy densities of particle i respectively. Note that after the phase space integration on the right hand side of eqs.(3.3-3.4), δn/δt and δρ/δt are functions of the temperature T , the chemical potential µ and the model parameters only. 4 Thus, in terms of the Hubble parameter, one can readily obtain the following equations through the application of the chain rule: These two equations effectively describe the evolution of T and µ for any particles in the early Universe, and can thus be used to solve the decoupling of neutrinos from the rest of the plasma as we will see later in this section.
3.1.1 Evolution of T γ , T ν and µ ν in the SM At the time of neutrino decoupling, since photons and electrons are still tightly coupled, one can safely set µ γ = µ e = 0 and T γ = T e . By applying eqs. (3.5-3.6), one obtains [125] dT γ dt = − 4Hρ γ + 3H (ρ e + p e ) + δρνe δt + δρνµ δt + δρντ δt ∂ργ Note that the above equation for T γ is derived assuming the finite temperature corrections are negligible, while it has been known for a long time that this is not the case especially given the precision measurements of N eff from future experiments. To be clearer, in the future, N eff will be measured to the percent level, while finite temperature corrections to N eff is also at the percent level [117][118][119][120][143][144][145][146]. Therefore, to correctly interpret the results from future experiments and/or to disentangle contributions to N eff from any potential new physics from the SM, the QED corrections have to be included.
The leading-order QED corrections were obtained decades ago [147,148], and higherorder corrections up to O(e 4 ) were recently calculated in Ref. [121], where the authors found corrections to N eff are about −0.0009 and 10 −6 at O(e 3 ) and O(e 4 ) respectively. Since both corrections at O(e 3 ) and O(e 4 ) exceed the proposed precision target of the future experiments, we neglect those in our setup and only keep the finite temperature corrections up to O(e 2 ). On the other hand, neutrino oscillations also lead to a correction to N eff , which is about 0.0007 as reported in Ref. [139,140,149]. Note that, as was pointed out in Ref. [121], since contributions to N eff from neutrino oscillations and the finite temperature corrections at O(e 3 ) are comparable, they shall both be included for a consistent precision calculation of N eff . However, as stated above, due to their smallness, we also neglect contributions from neutrino oscillations in this work.
To conclude this subsection, we show the result for T γ with the inclusion of the aforementioned finite temperature corrections following the notations of Refs. [121,125]: where P int and ρ int ≡ −P int + dP int /d ln T γ are finite temperature corrections to the electromagnetic pressure and the electromagnetic energy density respectively, whose analytical expressions can be found in Ref. [121].

Brief review of the collision term integrals
From eqs.(3.9) and (3.8), one can then solve T γ (t) and T ν (t), and thus N eff at the time of neutrino decoupling. From eqs.(3.3-3.4), we conclude that to solve T γ (t) and T ν (t), the remaining task is to first finish these phase space integrals, which are in general very challenging with no analytical expressions. This in turn slows down numerical calculation of N eff , especially in the presence of new physics. However, as pointed out in Ref. [125], analytical results for those collision term integrals exist in the Maxwell-Boltzmann limit, as a result, numerical calculation of N eff can be boosted significantly. For specific processes in the SM, the author of Refs. [125] presented analytical results for the collision terms integrals in Ref. [124], which, however, can not be generalized to processes in the presence of new physics. In light of this, and since the analytical forms of the collision term integrals are essential to boost the numerical calculation of N eff , we present a full generic and analytical dictionary for the collision term integrals in section 4, while present the method we use to obtain these collision term integrals in this subsection.
To start, we consider the collision term integrals for 2 → 2 processes as these are the only interaction types relevant for N eff calculation in the SM and with the inclusion of NSI operators considered in this work. 5 Upon leaving out the irrelevant factor g, for a generic process 1 + 2 → 3 + 4, one can write the collision term integrals on the right hand of eqs.(3.3-3.4) generically as with p i and E i the four-momentum and the energy of the i-th particle. To simplify the collision term integral, we reproduce the results presented in Appendix D of Ref. [161] and cite the result here: 6 where α (θ) is the angle between p 1 and p 2 ( p 3 ), and we define For decay or inverse decay, the collision term integral, defined is eq. (3.10) is a nine-fold one that can be reduced to a two-fold integral. There are many literatures in the past discussing the collision term integral, see, for example, Refs. [121,137,139,140,[149][150][151][152][153][154][155][156][157][158][159][160]. 6 We assume CP conservation that allows the factorization of M 2 1+2→3+4, which is a well-justified approximation for our purpose here.
Note that a ≤ 0 from above definition, and the integrating regions for α and θ are altered, where resulting from the requirement of the existence of physical solutions to the collision term integral. Note also that M 2 1+2→3+4 is in general a function of p ij , defined as which also depends on the angles α and θ, making the integral in eq. (3.11) too cumbersome to be completed. Surprisingly, one huge simplification that eventually allows the completion of the integral in eq. (3.11) can be realized when (1) all the particles involved are massless, i.e., m i = 0 (i = 1, . . . , 4), leading to Q = 0, min(1, cos θ + ) = 1 = cos θ + , (3.21) and (2) all the particles obey the Maxwell-Boltzmann distribution, permitting analytical expressions for almost all possible forms of M 2 1+2→3+4 . The only exception is when there exists a light mediator in the t and/or the u channels, where it has been well-known that IR divergence emerges when all external particles become massless, the Compton scattering for example. 7 However, this IR divergence has to cancel out for any sufficiently inclusive quantities, as is guaranteed by the Kinoshita-Lee-Nauenberg (KLN) theorem [162,163]. Recently, it is also shown in Ref. [164] that to have IR finiteness, one does not necessarily need to sum over both the initial and the final states as stated by the KLN theorem, rather, one only needs to sum over all possible final (initial) state for a given fixed initial (final) state, as long as the forward scattering is also included.
In this work, since we are interested in constraints on new physics from N eff in a model independent manner within the EFT framework, we will mainly focus on scenarios with heavy mediators such that the IR divergence issue mentioned above never shows up. The only exception is the dimension-5 neutrino magnetic dipole operator, which we will discuss in section 4. Furthermore, in order to obtain analytical results, as discussed in last paragraph, we assume all particles (1) are massless, and (2) obey the Maxwell-Boltzmann distribution only when calculating the collision term integrals from eq. (3.11). We then present a complete dictionary for all possible M 2 1+2→3+4 up to products with three p ij in M 2 1+2→3+4 in section 4, which are also provided in auxiliary Mathematica notebook files. Corrections from non-vanishing masses and Fermi-Dirac/Bose-Einstein distribution are also discussed in section 4.

EFT operators and the collision term integrals
As no new particles have been observed after the discovery of the Higgs particle in 2012 [1,165], EFTs have become the natural framework for the study of any new heavy physics. In this work, we are interested in corrections to N eff from higher-dimensional operators in the early Universe. The active degrees of freedom at that time are neutrinos, photons, electrons and positrons. Since the neutrinos decouple from the rest of the plasma at around 2 MeV, for the EFTs to be valid, the potential new physics could be as light as ∼ O(100 MeV) that is about the muon mass. Note that the lower bound of the new physics scale would also be constrained from, for example, Big Bang Nucleosynthesis. Given that future experiments like CMB-S4 could constrain ∆N eff < 0.06 at 95% CL [117,143,144,146] with ∆N eff the corrections to SM prediction of N eff , one naturally expects the EFT operators could also be constrained from the precision measurements of N eff .
In this section, we will first enumerate the relevant EFT operators up to dimension-7 in section 4.1. The resulting invariant amplitudes from these operators turn out to be functions of p ij defined in eq. (3.20), model parameters and the Wilson coefficients. Depending on how explicitly the M 2 depends on p ij , the collision term integral in eq. (3.11) needs to be calculated case by case. From momentum-energy conservation, the redundancy in collision term integral computation can be reduced to a set of limited number of bases as presented in section 4.2. Starting from these bases, we then present a complete dictionary of the collision term integrals in section 4.4.

List of relevant EFT operators
We start from the SMEFT, obtained by integrating out the heavy new degrees of freedom introduced to the SM, where the Lagrangian can be expressed as the SM Lagrangian, plus a tower of higher-dimension operators O (j) : where C j 's are the Wilson coefficients and Λ is the characteristic scale of new physics. In this setup, the neutrino masses can be naturally generated through the dimension-5 Weinberg operator [3]. However, for the rest of this work, we neglect the masses of neutrinos due to their tininess compared with the other scales involved in our calculation. In the early Universe where the active fields are neutrinos, photons, electron and positrons, the Universe can be described by the LEFT, where the top quark, the SU(2) gauge bosons, and the Higgs boson have also been integrated out within the SMEFT, inducing both CC and NC neutrino NSIs. See, for example, Ref. [50] for the discussion. However, we point out that CC and NC NSIs are not necessarily generated by heavy particles above the weak scale, instead, it can also be generated by some light particles above the O(100 MeV) scale. To illustrate this point, we briefly discuss a toy Z and a toy pseudoscalar models here: • The toy U(1) model we consider, without restricting ourselves to any other constraints such as anomaly cancellation, collider and cosmological constraints etc., is the Z

Dimensions Operators
Wilson coefficients Table 1. Effective operators relevant for N eff up to dimension-7 with α, β, α , β = e, µ, τ , the neutrino flavor indices, and f = e. Operators with ♣'s are the extra operators we consider in this work and the symbol "c" along with related operators means charge conjugation. The last column shows our convention for the Wilson coefficients.
model that can be written as where L and R are the left-handed lepton doublet and the right-handed lepton singlet under SU(2) L respectively, and Z is the new vector boson charged under the U(1) group. For our purpose, Z needs not to be above the weak scale, and as long as Z is above ∼ O(100 MeV) or equivalently the muon mass, Z can be integrated out: leading to the NC neutrino-electron and electron-electron contact interactions as seen above when m 2 Z is larger than the momentum transfer p 2 . During neutrino decoupling, since p 2 is of O(10 MeV) 2 , thus as long as m Z is above O(100 MeV), the EFT after integrating out Z serves as a good framework for the study of this new physics.
• The toy pseudo-scalar model we consider, without considering any theoretical and/or experimental constraints, can be expressed as where φ is the pseudo-scalar with mass m φ . Similarly, as long as m φ is above ∼ O(100 MeV), one can integrate out the particle φ, and obtain the contact neutrino self-interacting operators.
The CC NSIs have been recently studied in Refs. [64,65]. In Ref. [65], the authors took the running and the matching effects at different EFT scales into account, and the resulting constraint on the UV scale Λ was found to be as large as about 20 TeV from neutrino oscillation data. The NC operators are the relevant ones for our study in this work, part of which has been previously studied in Refs. [51,75,[107][108][109][110][111][112], and a comprehensive study was recently presented in Ref. [38]. Note that, since the authors in Ref. [38] were interested in constraints on these NSIs from neutrino experiments, they did not consider any dimension-6 neutrino self-interacting operators as neutrinos only feebly interact with our matter world. However, since these operators are closely related to N eff by modifying the neutrino number and the energy densities directly through neutrino self-interactions, we include these operators in this work and study constraints on these neutrino self-interacting operators from precision measurements of N eff . 8 We summarize all the EFT operators up to dimension-7 in table 1 following the notations in Ref. [38].
One immediate observation from table 1 is that, the dimension-5 operator in the first row, i.e., the neutrino magnetic dipole operator, corresponds to the light-mediator scenario we discussed at the end of section 3.2. When the intermediate photon shows up in the s-channel, there is no IR divergence, and the collision term integral can be calculated analytically. However, since the intermediate photon can also appear in the t-and/or u-channels for the ν α ν β → γ * → ν α ν α 9 process for example, the collision term integrals would exhibit the IR divergence discussed earlier. In principle, one could remove this divergence by applying the KLN theorem or following the procedure discussed in Ref. [164] for any inclusive observables. However, since the scenario with a light mediator resides in a different regime compared with all the other operators listed in table 1, we leave the light mediator scenario for a future project. We also point out that this operator is very stringently constrained from the magnetic moment of ν e using Borexino Phase-II solar neutrino data [38,171], which justifies our ignorance of the O (5) 1 operator for the calculation of N eff . We will discuss more on this in section 5. Now, including corrections from the dimension-6 and dimension-7 operators in table 1, the invariant amplitude M 2 1+2→3+4 in eq. (3.11) can generically be written as: where M SM and M EFT are the amplitudes from L SM and the EFT operators in table 1 respectively. Clearly, when the potential new physics scale Λ is about or above O(Λ W ) with Λ W the weak scale, the interference term in eq. (4.5) would be of the same order as the SM contributions, and the M 2 EFT term can then be safely neglected. However, we point out that in the case where Λ Λ W , though contributions from the EFT operators dominate, one can not simply discard the M 2 SM term in eq. (4.5) since for some of the operators in table 1, for example, the O 6 3,4,5 operators, the M 2 SM is the only part that tells how T γ and T να evolve with time as we will see below. In light of this, we always keep all the three terms in eq. (4.5) during our calculation, while using the large Λ limit to cross check our results.
By plugging in M SM and M EFT in eq. (4.5), one obtains M 2 1+2→3+4 as a function of p ij , the SM model parameters, the scale of new physics Λ, and the Wilson coefficients C j in the last column of table 1. One can then calculate the collision term integrals through finishing the integral shown in eq. (3.11). However, due to momentum-energy conservation, redundancy exists in the calculation of collision term integrals. This redundancy can be sufficiently removed by first choosing a set of basis, which we discuss in the next subsection.

Choices of the independent bases
To start, we realize that for both the SM and the EFT contributions, the relevant processes are either 2 → 2 scattering or 2 → 2 annihilating processes. Denoting these processes generically as 1 + 2 → 3 + 4 with momentum p i for the i-th particle, the invariant amplitude M 2 1+2→3+4 , defined in eq. (4.5), can be expressed as a tower of the momentum scalar product p ij defined in eq. (3.20): k's up to k = 3 in eq. (4.6), which then gives, by leaving out the arguments of c ij···mn···st 's, Note that for k = 1, and similarly for the k = 2, 3 cases, we have included contributions from i = j in the c 0 term from the on-shell conditions. On the other hand, terms in the second and the third lines of eq. (4.7) are not all independent from momentum-energy conservation, resulting in redundancy when one calculates the collision term integrals from eqs. (4.7) and (3.11). However, this redundancy can be sufficiently removed by first choosing an independent basis in terms of p ij , and then rewrite M 2 1+2→3+4 as a linear combination of these bases. Depending on k, the independent bases we choose are presented in table 2. One can then readily express the momentum tower as linear combinations of these bases. For example, → p 12 · p 13 · p 14 in the massless limit. (4.9) At this stage, the collision term integral in eq. (3.11) boils down to the collision term integral with M 2 1+2→3+4 being replaced by the independent bases shown in table 2. Particularly, when all the masses involved are vanishing, M 2 1+2→3+4 generically simplifies significantly as seen from the example above. Therefore, as also discussed at the end of section 3.2, to obtain analytical results for the collision term integrals, we assume m i = 0 (i = 1, . . . , 4) 10 and all particles obey the Maxwell-Boltzmann distribution. We then present the complete generic and analytical dictionary of the collision term integrals in section 4.4. Corrections from finite electron mass m e , spin statistics and neutrino chemical potentials are discussed in section 4.3.

Corrections from m e , spin statistics and chemical potentials
As already noticed in Ref. [125], corrections from finite electron mass m e and spin-statistics have to be included to reproduce N eff obtained from the density matrix formalism. Furthermore, as one can see from Table 1 of Ref. [125], these corrections are of the same order as the finite temperature corrections discussed at the beginning of this section. Thus, to be consistent, these corrections have to be included.
In Ref. [125], finite m e corrections are obtained by finding the ratios of the collision term integrals in eq. (3.11) by switching on and off m e . Similarly, corrections from spin statistics are computed by finding the ratios of the collision term integrals with Fermi-Dirac/Bose-Einstein and Maxwell-Boltzmann distributions respectively. For a detailed discussion, see Refs. [125,126]. Though the methods used for the collision term integrals are different, we reproduce the numbers in Table 6 of Ref. [125] and/or Table III of Ref. [126]. These corrections are then included in the Boltzmann equations and used to solve N eff . The results are presented and discussed in more detail in section 5.
We also comment on that neutrino chemical potentials are highly suppressed due to the rapidνν ↔ e + e − ↔ γγ conversion and that the electron chemical potential is negligibly small compared to the plasma temperature in the early Universe. In our setup, in order to be generic, we keep neutrino chemical potentials throughout our analytical calculations, but stress that they would have no visible impact on the current/planned precision measurement of N eff . For our numerical calculations, we choose µ ν = µν and |µ ν /T γ | = 10 −4 with T γ = T ν = 10 MeV as our initial conditions, and verify that this setup is numerically equivalent to vanishing neutrino chemical potentials as expected.

A complete generic and analytical dictionary of the collision term integrals
In last subsection, we list in table 2 the independent bases by which the invariant amplitudes M 2 1+2→3+4 can be expressed, and conclude that the redundancy of collision term integrals from momentum-energy conservation can be removed by working with these bases directly. In this subsection, we provide the complete analytical dictionary of the collision term integrals for particle "1" and up to k = 3, with k the number of p ij 's in the invariant amplitude. We note that a subset of this complete dictionary was presented in the appendices of Ref. [124,126], which agrees with our results presented in this subsection as long as one specifies T i and µ i accordingly.
For the collision term integrals, we follow the procedure briefly summarized in section 3.2 and stick to our notation in eq. (3.11), where j = 0 represents the collision term integral for the number density and j = 1 that for the energy density. The dependence on p ij of M 2 1+2→3+4 is reflected by the argument of C (j) . 11 For example, C (0) (p 12 ) means the collision term integral for the number density with M 2 1+2→3+4 = p 12 12 in eq. (3.11).
From eq. (4.10), one notes that in the j = 0 case, the collision term integral, corresponding to the number density, vanishes when T 1 = T 3 and T 2 = T 4 . This is expected since it actually stands for a scattering process where the number density for each species remains the same before and after the interaction. This conclusion holds generically and is independent of the form of M 2 1+2→3+4 , as one can also see clearly from the results below. On the other hand, if T 1 = T 2 = T 3 = T 4 and µ 1 = µ 2 = µ 3 = µ 4 , then all the C (j) 's vanish, which is also as expected since particle self interactions do not modify the number and the energy densities as long as thermal equilibrium is maintained. This observation also acts as a cross-check of our analytical results presented in these subsections.

4.4.2
13) 11 We stress that the argument of C (j) here is only used to reflect the dependence on pij of M 2 1+2→3+4, and this argument does not mean that C (j) depends on pij. Instead, C (j) only depends on the model parameters, the new physics scale, the Wilson coefficients, the temperatures Tγ,ν α and the chemical potentials µν α . 12 We leave out any overall factors in M 2 1+2→3+4 that are independent of pij here and in the following.
The complete results of M 2 from SM and the EFT operators listed in table 1 are given in an auxiliary Mathematica notebook file for all relevant processes, together with all the replacement rules to rewrite M 2 in terms of the bases listed in table 2. 13 We point out that when all the Wilson coefficients vanish therein, we reproduce the SM results as presented in, for example, Ref. [173]. Using the complete dictionary summarized in this section and building our code upon nudec_BSM from Ref. [125], we study corrections to N eff from the NSI operators in table 1, and discuss the results in section 5.

Constraints on EFT operators from N eff
With the complete dictionary presented in section 4, one can readily solve the Boltzmann equations for T γ and T να , and thus obtain corrections to N eff . In what follows, we define these corrections as where N SM+EFT eff is the theoretical prediction of N eff with the inclusion of the NC NSI operators, and N SM eff = 3.044 [123,132] that from the pure SM. For Planck, we use the current result N eff = 2.99 +0. 34 −0.33 [114] at the 95% CL to obtain the constraints, and ∆N eff < 0.06 at 95% CL for CMB-S4 [117,143,144,146].
Our code is built upon nudec_BSM from Ref. [125], and is then used to solve eqs. (3.8-3.9) numerically by Mathematica. During our numerical solutions, we keep terms proportional to m e in M 2 and assume T νµ = T ντ = T νe and µ νµ = µ ντ = µ νe . In the very large Λ limit, we reproduce the results in Table 1 of Ref. [125] for both T νe = T νµ,τ and T νe = T νµ,τ . We then show our results for varying Wilson coefficients or the new physics scale Λ in the following subsections.

Constraints on Λ with fixed Wilson coefficients
Following the notations clarified in table 1 and fixing the Wilson coefficients at unity, we present our results in figure 1. The constraints shown in figure 1 are obtained by assuming only one nonvanishing NSI operator at a time, and the results are presented from considering the latest results from Planck [114] in orange and the proposed precision goal of CMB-S4 [117,143,144,146] in purple. Several points from this plot merit emphasizing: [123,132] is the SM prediction of N eff , and N SM+EFT eff is that from SM and new physics. Note that the plot is obtained by fixing all Wilson coefficients at unity and considering only one non-vanishing operator only at a time. See the main text for more discussion.
• Constraints on dimension-6 EFT operators are generically stronger than those on the dimension-7 ones. The reason is that dimension-7 operators are more suppressed by one more power of Λ. Moreover, among the dimension-6 operators, currently, the Planck data leads to the most stringent constraint on the O 1,e operator, whose lower bound is presently constrained to be about 195 GeV. In the future, CMB-S4 would improve this lower bound to about 240 GeV, as one can see from the first purple histogram in figure 1. Quantitively, we summarize the lower bounds on Λ's for all the operators shown in figure 1 in table 3.
• As one can see from table 3, for the O 2,e results from the destructive interference between the SM and the NSI operators, while that for O (7) 5,f results from a negative shift to N eff when Λ is small. We show this point in the second row of figure 2.
• Constraints on dimension-6 operators O  3.47 6.17 Table 3. Constraints on EFT operators from current Planck data and future CMB-S4 proposal. All lower bounds are obtained by assuming one non-vanishing EFT operator at a time and fixing the Wilson coefficients at unity. Note that for O 2,e and O 5,f , there are exception intervals for CMB-S4 as a result of destructive interference or a negative shift in N eff from the EFT operators. See the main text for more discussion.
one might naïvely think the M 2 SM term in eq. (4.5) can be safely discarded and very large N eff would be predicted from O 3,4,5 . However, as we already point out right below eq. (4.5), the SM part can not be ignored since in this case, it is the only part that governs the evolution of T γ . Furthermore, when Λ Λ W , neutrino self interactions are rapid enough to eliminate any difference between T νe and T νµ,τ , and neutrinos of all flavors have exactly the same temperature. 14 This in turn results in vanishing corrections to the collision term integrals for O (6) 3,4,5 as discussed right after eq. (4.10). Thus, when Λ is very small, corrections from O to N eff vanish and the SM prediction is restored.
• While it is a very good approximation to neglect the temperature and the chemical differences among neutrinos when calculating N eff within the SM framework, see Table 1 of Ref. [125] for example, this approximation does not stay valid any more in the case where new physics introduces only neutrino self interactions as the O 3,4,5 operators. In this scenario, if one takes the equal neutrino temperature and the equal chemical potential approximation for all the three-flavor neutrinos, then the collision term integrals simply vanish such that the effects of this new physics can never be observed.
• Constraints on O are also missing in figure 1, due to the suppression factors α/(12π) and α/(8π), respectively: At the amplitude level, both these two factors lead to suppression of O(10 −4 ), thus the invariant amplitude M 2 is suppressed by a factor of 14 With our choice of the Wilson coefficients and the small Λ, neutrino self-interacting rates are always larger than the Hubble rate such that the three flavor neutrinos always stay in thermal equilibrium. However, neutrino decoupling is not affected since photon-electron-positron sector is governed by weak interactions. We emphasize that the equal temperature of neutrinos is the direct result of neutrino self-interactions introduced by O 3,4,5 , the moderate Wilson coefficients, and the small Λ.
5,e only. In all these subfigures, the black curve stands for corrections to N eff from the related NSI operator, the horizontal red dashed line is the constraint on ∆N eff from Planck, and the horizontal red dashed line is that from CMB-S4.
O(10 −8 ). Note also that there is no interference between the SM and O 1 , it was concluded in Ref. [38] that the most stringent lower bound on Λ was 2.7 × 10 6 GeV from the magnetic moment of ν e using Borexino Phase-II solar neutrino data [171]. On the other hand, translating this constraint on the magnetic moment of ν e from O 1,2 , they found Λ > 328 GeV and Λ > 1081 GeV for O 1 and O 2 respectively. Furthermore, O 1,e was also constrained to have a lower bound of 1005 GeV from a global fitting of neutrino oscillating data [38,106]. As one can see from table 3, the constraint on O 1,e from the global fitting is stronger than that from N eff . However, all the other operators in figure 1 are not constrained in Ref. [38], making our work complementary to theirs as well as that in Ref. [65].
We emphasize that conclusions above are obtained by fixing the Wilson coefficients at one and considering only one non-vanishing NSI operator at a time. In Ref. [65], we find that if multiple operators exist at the same scale, then the correlation among them may change the constraints by orders of magnitude. However, due to the computation challenge, this correlation effect is in general ignored except for some UV models where the number of operators at certain dimension is limited. In this work, we find when Λ ∼ Λ W or smaller where NSI contributions to M 2 are comparable to or dominate over those from the SM, numerical computation of the Boltzmann equations is extremely slow or even impossible even with high-performance clusters. For this reason, the correlation mentioned above will not be discussed.

Constraints on Wilson coefficients with fixed Λ
Alternatively, we present the constraints for the Wilson coefficients with Λ = 1 TeV and 100 GeV in this subsection, and the results are shown in the first and the second rows of figure 3, respectively. Constraints are shown for O (6) 1,e and O For C 2,e , the two exception intervals for CMB-S4 in the last line of the two bullets above result from destructive interference as already discussed in last subsection -For C 1,e G F sin 2 θ W T 9 γ π 5 Λ 2 , (5.6) where θ W is the weak mixing angle and we only show the interfering terms here by omitting any sub-leading effects in them in each case. Note that a larger (smaller) neutrino energy density would be equivalent to a higher (lower) neutrino temperature. Thus, as can be understood from eq.(2.3), the constructive (destructive) interference also explains the positive (negative) shift feature of N eff from O 1,e (O 2,e ) in figures 2 and 3 when Λ ∼ Λ W . For the other NSI operators not shown in figure 3, since they are at least suppressed by one more power of Λ, Planck and CMB-S4 are not able to constrain those Wilson coefficients when Λ is fixed at 1 TeV. Similar observation is obtained for Λ = 100 GeV, with stronger constraints on C 1,e and C (6) 2,e , whose magnitudes are 100 times smaller compared with the Λ = 1 TeV case as expected.

Comparison with current constraints on NC NSIs
To compare with constraints on the dimension-6 operators O (6) 1 (2),e from other experiments, we first review the parameterization commonly used in the literatures to describe neutrino NSIs: with f = e, u, d the charged fermioins, α, β = e, µ, τ the flavor of neutrinos, and P = L, R the chiral projector operators where L = (1 − γ 5 )/2 and R = (1 + γ 5 )/2. Note that the nine f,P αβ 's are all real, and Hermiticity of the Lagrangian guarantees that only six of them are independent. The relevant operators for our study in this work are e,P αβ . One readily finds, in terms of the Wilson coefficients used in this work, the parameters can be expressed as 1,e + C Fixing Λ 174.10 GeV from the Λ 2 · 2 √ 2G F = 1 condition such that the LEFT in our notation mimics that in eq. (5.8), we present our results for L,R αβ in figure 4 by including all NC NSIs in eq. (5.8) while ignoring all neutrino flavor dependence of C (6) (1,2),e . The colors in each subgraph of figure 4 have exactly the same meaning as those in figure 2, and to obtain the constraints, we once again ignore the neutrino flavor dependence of the LEFT Wilson coefficients and consider only one non-vanishing at a time in our analysis. However, we stress that, as one can see directly from eq. (5.9), one non-vanishing in general includes contributions from both O In comparison, we list constraints on these parameters from previous studies and ours obtained in this work in table 4 by ignoring bounds from loops [24,94]. Note that constraints from Ref. [106] in the second column are originally presented in terms of e,L+R αβ ≡ e,L αβ + e,R αβ . We translate them on individual e,L(R) αβ by assuming only one of them non-vanishing. Constraints from TEXONO are obtained at the Kuo-Sheng Nuclear Power Station in Ref. [100], the LEP, LSND and CHARM-II experiments in Ref. [24], a global analysis of ν e e andν e e scattering data from LSND, Irvine, Rovno and MUNU experiments in Ref. [86], a combination of OPAL, ALEPH, L3, DELPHI, LSND, CHARM-II, Irvine, Rovno and MUNU experiments in Ref. [87], solar and reactor neutrino experiments in Ref. [88], low-energy solar neutrinos at source and detector from the Borexino experiment in Ref. [93], a global analysis of short baseline νe andνe data from LSND, LAMPF, Irvine, Rovno, MUNU, TEXONO and KRANOYARSK in Ref. [101], and the DUNE experiment in Ref. [39]. For constraints from Ref. [86], we cite their results in the one-parameter case since it leads to the most stringent constraints on these NC NSIs, and similarly for results in Refs. [93,101]. For constraints from Ref. [39], the upper and the lower intervals are obtained using an exposure of 300 and 850 kt.MW.yr for DUNE respectively. For constraints from Ref. [93], the upper number is obtained from a detector-only study using low-energy solar neutrinos at Borexino, while the lower is the future prospect from a combined analysis of the detector and the source. For all the other cases 15 These results can be readily obtained by using our complete dictionary in section 4.4 or the analytical expressions in the auxiliary Mathematica notebook file. 16 Since we ignore the neutrino flavor dependence, we thus leave out the neutrino flavor indices here and in the following.  Table 4. Summary of constraints on dimension-6 neutrino-electron NC NSIs from previous studies and this work. Constraints from a global fitting of all kinds of neutrino oscillation data plus the COHERENT result are obtained in Ref. [106], the TEXONO collaboration in Ref. [100], the LEP, LSND and CHARM-II experiments in Ref. [24], a global analysis of ν e e andν e e scattering data from LSND, Irvine, Rovno and MUNU experiments in Ref. [86], OPAL, ALEPH, L3, DELPHI, LSND, CHARM-II, Irvine, Rovno and MUNU experiments in Ref. [87], solar and reactor neutrino experiments in Ref. [88], low-energy solar neutrinos at source and detector from the Borexino experiment in Ref. [93], a global analysis of short baseline νe andνe data from LSND, LAMPF, Irvine, Rovno, MUNU, TEXONO and KRANOYARSK in Ref. [101], and DUNE in Ref. [39].
in table 4, whenever two intervals appear, it means two "disjoint" ranges that are simultaneously allowed from their analyses. We refer the reader to the original references for more details.
As one can see from table 4, in general, constraints from other experiments are stronger than those we obtain from Planck. However, from the last column of table 4, the results from CMB-S4 would be improved by a factor of ∼ 3 (5) for e,L(R) . As a result, all the 's would be bounded at the 10% level. On the other hand, in table 4, one notes that seven of these 's are constrained at the 10% level from previous experiments, except the following five's: e,L ee [87,88], e,(L,R) µµ [24,87], e,L µτ [106], e,R ee [100], and e,R µτ [106] that are stringently constrained at the 1% level. Therefore, constraints on most of these 's from CMB-S4 are basically comparable to the existing ones. For example, On the other hand, taking both e,L and e,R into account, we obtain simultaneous constraints on e,L and e,R as shown in figure 5, where the orange and the purple regions are still allowed by Planck and CMB-S4 respectively. The permitted regions are along the diagonal region on the e,Le,R plane since it is where contributions from O from CMB-S4 near the origin. Since we assume neutrino flavor independence and N eff is more sensitive to light degrees of freedom, our constraints are slightly less stringent than, but again very comparable to, those discussed in last paragraph. Our results presented in this work complement those from previous studies on NC neutrino NSIs from collider, neutrino coherent scattering and neutrino oscillation experiments.

Conclusions
Null observation of any new resonances after the discovery of the Higgs particle at the LHC has gradually changed our strategy in searching for new physics from specific UV models to modelindependent studies. EFTs provide a systematic and model independent approach to heavy new physics. In the early Universe where the active fields are neutrinos, electrons, positrons and photons, the system can be described by the LEFT, even with the introduction of some new physics above the ∼ O(100 MeV) scale. NC NSIs induced by the new physics would affect neutrino decoupling in the early Universe, thus would also modify the prediction of N eff . In light of the very precision measurements of N eff from current Planck data and the precision target from CMB-S4, we present constraints on NC NSIs from N eff up to dimension-7 in this work by assuming that all NC NSIs are induced by heavy mediators above ∼ O(100 MeV).
To that end, we adopt the strategy developed in Refs. [124,125], which permits a fast and precision calculation of N eff , and can also be easily generalized to include various new physics. The fast calculation of N eff largely seeds in the pre-calculated collision term integrals, which are only obtained for several specific processes in the SM. In this work, we provide a complete, generic and analytical dictionary for these collision term integrals in section 4. With our results, as long as the invariant amplitudes are known, one can refer to this dictionary to write down the Boltzmann equations, and then solve the prediction of N eff from the SM or some new physics with few efforts. We also show an example for the application of this dictionary at the end of section 4.
Including the NC NSIs and using the dictionary described above, we study constraints on these operators from precision measurements of N eff . Our results are presented in figure 1 and summarized in table 3, where the lower bounds on the scale of new physics Λ is obtained by fixing the Wilson coefficients at unity and considering only one non-vanishing NSI operator at a time. We find that, the dimension-6 NSI operators O (6) 1,e and O (6) 2,e are constrained to be above ∼ 331 GeV and ∼ 240 GeV respectively from CMB-S4. On the other hand, due to suppression from the new physics scale, the couplings and m e , dimension-7 operators O (7) (5,6,7,8,9,10,11),e only have visible corrections to N eff when the new physics is relatively light, thus the current lower bounds on these operators are about 6 GeV and 3 GeV for O (7) (7,8,9,10,11),e and O (7) (5,6),e , respectively. Operators O (6) 3,4,5 are not constrained from N eff due to (1) negligible corrections to N eff when Λ Λ W and (2) realization of thermal equilibrium among the three flavor neutrinos that results in vanishing contributions to N eff . Operators O (1,2),e are also not constrained from N eff due to suppression of tiny couplings. On the other hand, we also study constraints on the Wilson coefficients with Λ fixed at 1 TeV and 100 GeV. The results are shown in figure 3 by taking only one non-vanishing NSI operator into account at a time. We find that only C 1,e and C (6) 2,e are constrained by N eff since the dimension-7 operators are all suppressed by one more power of Λ, as well as m e and the small couplings. At Λ = 100 GeV, we find the magnitude of C 1,e is constrained to be around 0.3 (0.1) from Planck (CMB-S4), while it is about 1.4 (0.2) for C (6) 2,e from Planck (CMB-S4). The results are summarized in eqs. (5.2-5.5).
Constraints on the dimension-6 neutrino-electron NC NSI operators O (1,2),e from precision measurements of N eff are also compared with previous results from, for example, a global fitting of neutrino oscillation experiments and collider experiments. To that end, we first obtain constraints on the NC NISs using the parameterization, and then present the results in figure 4 and table 4 for one non-vanishing NC NSI operator at a time, and figure 5 for the inclusion of both operators. We find that constraints from precision measurements of N eff from Planck are in general weaker than those from other experiments mentioned above. However, the improved results from CMB-S4 in future would become comparable for certain operators. Our work complements previous studies on NC NSIs from other experiments. In the future, if the cosmic neutrino background (CνB) could be directly measured, N eff would be determined with a much better precision, and one could then expect also much stronger constraints on these NC neutrino NSIs from CνB.