Neutrino Magnetic Moments Meet Precision $N_{\rm eff}$ Measurements

In the early universe, Dirac neutrino magnetic moments due to their chirality-flipping nature could lead to thermal production of right-handed neutrinos, which would make a significant contribution to the effective neutrino number, $N_{\rm eff}$. We present in this paper a dedicated computation of the neutrino chirality-flipping rate in the thermal plasma. With a careful and consistent treatment of soft scattering and the plasmon effect in finite temperature field theories, we find that neutrino magnetic moments above $2.7\times 10^{-12}\mu_B$ have been excluded by current CMB and BBN measurements of $N_{\rm eff}$, assuming flavor-universal and diagonal magnetic moments for all three generation of neutrinos. This limit is stronger than the latest bounds from XENONnT and LUX-ZEPLIN experiments, and comparable with those from stellar cooling considerations.

Cosmological observables such as the effective neutrino number, N eff , can also be used to constrain NMM. For Dirac neutrinos, NMM could flip the chirality and thermalize righthanded neutrinos, which would contribute to N eff significantly. Under this line of thought, previous studies derived cosmological bounds on NMM, µ ν < O(1) × 10 −11 µ B [34,35] and µ ν < 2.9 × 10 −10 µ B [36]. Alternatively, one may consider Majorana neutrinos which can only possess transition magnetic moments. In this case, no additional species are produced but large NMM could modify the SM neutrino decoupling at the MeV epoch. However, the cosmological constraint on this scenario is found to be weak [37].
Our work contains a careful calculation of the chirality-flipping rate, which was not treated consistently when an infrared (IR) divergence is involved in previous studies [34][35][36]. In Ref. [34], the chirality-flipping rate was obtained from a straightforward computation with a naive cut for the momentum transfer and the result exhibited a logarithmic dependence on the cut. In Ref. [35], the authors computed the rate in the real-time formalism of thermal quantum field theory (QFT), with a resummed photon propagator in the Hard-Thermal-Loop (HTL) approximation [57,58], but the imaginary part of the neutrino self-energy is time-ordered rather than retarded 2 . Later in Ref. [36], the chirality-flipping rate was computed with a retarded loop amplitude under the HTL approximation [57,58], and a momentum-cut approach was used to separate the hard-and soft-momentum transfers [60]. These different treatments are one of the reasons that lead to the different upper bounds of µ ν mentioned above.
In our more elaborated computation of the chirality-flipping rate, we adopt the realtime formalism of thermal QFT and take into account the resummed photon propagator with the damping rate not limited to the HTL approximation. We present the master integral for the collision rate that automatically contains contributions of both photonmediated scattering processes and plasmon decay γ * →ν L + ν R . Since the separation of hard-and soft-momentum contributions is known to be nontrivial [36,61], we will follow a numerical approach to compute the master integral, which allows us to obtain the collision rate more efficiently. Our calculation is focused on the Dirac neutrino case, but the collision rate computed in this work can be readily applied to Majorana neutrinos.
The paper is organized as follows. In Sec. 2, we start with a straightforward calculation Figure 1. The t-and s-channel scattering processes for ν R production in the early universe. The red blobs denote neutrino magnetic moments and ψ denotes a generic charged fermion. of the collision rate from the tree-level scattering amplitude, and show that the calculation relies on the infrared cut imposed on the momentum transfer. A more consistent treatment requires loop calculations in the thermal QFT, which will be elaborated in Sec. 3. Using the obtained collision rate, we compute the NMM correction to N eff and derive cosmological bounds on NMM in Sec. 4. Finally, we draw our conclusions in Sec. 5.
2 ν L → ν R flipping rate from tree-level scattering amplitudes The effective Lagrangian of a Dirac NMM is formulated as where ν denotes the Dirac spinor of a neutrino, σ αβ = i[γ α , γ β ]/2, and F αβ = ∂ α A β − ∂ β A α is the electromagnetic field tensor. In the chiral basis, ν = (ν L , ν R ) T , we can write it as which implies that the NMM operator flips the chirality of the neutrino. Neutrinos could also possess electric dipole moments, L ⊃ ν νσ αβ iγ 5 νF αβ /2, similar to Eq. (2.1) except for an additional iγ 5 . Our calculations for NMM can be applied to electric dipole moments by simply replacing µ 2 ν → 2 ν in the ν L → ν R rates because the squared amplitudes of all processes considered in this work are not affected by the additional iγ 5 , which can be seen manifestly in terms of Weyl spinors.
In the presence of such a chirality-flipping interaction, ν R can be produced in the thermal bath of the early universe via the tree-level scattering processes (see Fig. 1), ψ + ψ → ν L + ν R and ν L + ψ → ν R + ψ, where ψ denotes a generic charged fermion. In this paper, we assume that NMM are flavor universal and flavor diagonal for simplicity, which implies that each NMM is responsible for the thermal production of a single species of ν R .
Let us first calculate the chirality-flipping rate for the tree-level scattering processes in the zero-temperature limit and then show the necessity of considering finite-temperature corrections. For the s-channel process ψ + ψ → ν L + ν R , the squared amplitude reads where α ≡ e 2 /(4π) is the fine-structure constant, s, t, u are the Mandelstam variables, and all the fermion masses are neglected. The cross section is found to be constant in energy, The cross section is to be used in the Boltzmann equation, where n ν R denotes the number density of ν R , H is the Hubble parameter, and C (gain) ν R and C (loss) ν R are collision terms accounting for the gain and loss of ν R due to reactions. Since our discussions will be mainly on C for brevity. The s-channel contribution to C ν R , denoted by C ν R ,s , is computed as follows [62]: where T is the temperature of the SM thermal bath and K 1 is the modified Bessel function of order one. In Eq. (2.6) we have taken the same approximations adopted in Ref. [62] such as neglecting the Pauli-blocking effects and using the Boltzmann distributions.
For the purpose of studying when ν R can be in thermal equilibrium, we define the thermally averaged collision rate, where n eq ν R is the equilibrium value of n ν R . For the s-channel collision term, it reads With σvn defined, the condition of ν R in thermal equilibrium is formulated as The Hubble parameter is determined by GeV is the Planck mass and g is the number of effective degrees of freedom. Neglecting the temperature dependence of g , we see that H is proportional to T 2 while σvn is proportional to T 3 . Therefore, at a sufficiently high temperature, Eq. (2.9) is always satisfied. By solving H = C ν R /n eq ν R with respect to T , one can obtain the temperature of ν R decoupling, T dec . For T > T dec , ν R is in thermal equilibrium with the SM plasma. For T < T dec , ν R is decoupled and its temperature can be computed using entropy conservation.
Next, let us include the t-channel contribution to the collision rate. The squared amplitude of ν L + ψ → ν R + ψ is given by which is known to have an IR divergence in the soft-scattering limit, t → 0. In previous studies [34,63], an IR cut was imposed manually on the momentum transfer. This is qualitatively correct since the photon at finite temperatures has modified dispersion relations and acquires an effective thermal mass m γ ∼ eT in the relativistic QED plasma [57,58].
Keeping the photon thermal mass as an IR regulator, the cross section reads (2.11) Using m γ = eT / √ 6 to be derived in Sec. 3, we obtain the thermally averaged collision rate: where both ψ + ν L → ν R + ψ andψ + ν L →ψ + ν R have been taken into account. By comparing Eq. (2.12) to Eq. (2.8), one can see that σvn t / σvn s ≈ 55, which implies that the t-channel scattering dominates the ν R production, at least according to the above calculation with the simple IR regulator. A precise computation of ν R production rate that can consistently remove the IR divergence involves thermal QFT, as we will present in the next section.
3 ν L → ν R flip rate at finite temperatures

Computation method
The IR divergence in the t channel can be canceled by taking the finite-temperature effects into account. Previously, a momentum-cut approach was introduced to separate the hardand soft-momentum contributions [60]. The former is calculated in the tree-level scattering amplitude with a vacuum photon propagator, while the latter takes into account the resummed photon propagator at finite temperatures. Although it has been shown that the dependence on the momentum cut would be canceled after combining the hard-and soft-momentum contributions, this approach is based on some particular sum rules of the thermal propagators [36,61]. Sometimes finding the rules is a nontrivial task. In particular, the analytic extraction of the momentum-cut dependence in the soft regime is not so simple as that in the hard domain. If we are only concerned with the total rate, a full momentum integration that automatically combines the hard and soft regimes could be more efficient, as applied in Ref. [35]. Below we present the calculations in details. Readers who are not interested in the finite-temperature calculations are referred to Tab. 1 for the final results.

Collision rate in the Boltzmann equation
The evolution of phase-space distribution function of ν R is governed by the following Boltzmann equation, Figure 2. The self-energy diagram of ν R in the presence of neutrino magnetic moments (red blobs). The photon propagator is resummed to include the charged-fermion loop. In thermal QFT, the evaluation of this diagram can be used to obtain the total ν R production rate which automatically includes the contributions of plasmon decay and scattering processes. For instance, the dashed cut leads to a tree-level diagram corresponding to the on-shell scattering shown in Fig. 1. where Γ ν R ,gain/loss denotes the gain/loss rate of ν R , respectively. Their sum, Γ ν R ,tot ≡ Γ ν R ,gain + Γ ν R ,loss , can be physically interpreted as the rate of f ν R evolving towards equilibrium [64] 3 . At finite temperatures, it is related to the imaginary part of the retarded ν R self-energy as follows [64] (see also Appendix A for further details): where Σ R (p µ ) denotes the retarded self-energy of ν R and p µ = (E p , p) is the neutrino fourmomentum. In the real-time formalism of thermal QFT, the main task is to calculate the imaginary part of Σ R . The diagram of Σ R , as shown in Fig. 2, consists of a ν L propagator and a resummed γ propagator which includes the contribution of charged-fermion loops. The tree-level scattering amplitudes in Fig. 1 correspond to a half of the loop diagram after cutting it symmetrically along the blue dashed line in Fig. 2. According to the optical theorem, the squared amplitude of the tree-level diagram can be computed from the imaginary part of the loop diagram. Since in the resummed photon propagator an infinite number of one-particle-irreducible (1PI) loops are actually included, one can also add another 1PI loop to the photon propagator in Fig. 2 and then cut it symmetrically. This corresponds to the contribution of plasmon decay.
The collision term C ν R for ν R production previously introduced in Eq. (2.5) is related to Γ ν R ,gain via where f eq ν R is the distribution function in thermal equilibrium, and we have used the unitary condition Γ ν R ,gain = f eq ν R Γ ν R ,tot [64]. Figure 3. Two contributions, Σ −+ (p) and Σ +− (p), to the imaginary part of the retarded righthanded neutrino self-energy. The black blob denotes the charged-fermion loop in the resummed photon propagator, as shown in Fig. 2.

Neutrino self-energy from NMM interaction
Let us now compute the retarded self-energy of right-handed neutrino Σ R (p). The imaginary part of retarded self-energy can be evaluated by 4 where Σ −+ (p) and Σ +− (p) represent the two amplitudes in Fig. 3. The explicit forms are where k ≡ p + q and q denote the momentum of ν L and the photon. S ±∓ denote the free thermal propagators of the neutrino [57]: where f ν L is the thermal distribution function of ν L . The last parts of Eqs. (3.5) and (3.6), G ±∓,µν , denote the resummed photon propagators, which will be elucidated in Sec. 3.4.

Polarization tensors
For the retarded photon self-energy amplitude, Π R,µν , the most general tensor structure obeying the Ward identity (q µ Π R,µν = 0) is given by 5 where the longitudinal and transverse polarization tensors are given by 6 [65] (3.10) and u µ the 4-velocity of the plasma. In the rest frame, u µ = (1, 0) and the polarization tensors reduce to The retarded propagator for free photon is given by where the free 2 × 2 thermal propagator matrix elements, G ±±,µν (q), are given by and G −−,µν (q) = −G * ++,µν (q). Then, the Dyson-Schwinger equation in terms of Π L,T R reads where we have used the orthogonality, for A = L, T , and The modifications to the photon dispersion relation arising from the longitudinal and transverse parts are now encoded in the different thermal scalar functions Π L,T R (q). Given the free and resummed retarded propagators, we can compute any thermal components of resummedG AB,µν for A, B = ± via the diagonalization approach [66]. Explicitly, we havẽ where the diagonalization matrices U, V are given by [66] where f B/F are distribution functions for bosons and fermions respectively. Note that the physical result is independent of the unspecified scalar functions b q , c q . The diagonal matri-cesĜ µν ,Π µν are constructed by retarded/advanced propagators G R/A and loop amplitudes with G A,µν = G * R,µν and Π A,µν = Π * R,µν . We can then obtainG +−,µν ,G −+,µν in the retarded self-energy amplitude of ν R as follows: where the spectral densities ρ L,T (q) are given by (3.25) In the free limit, Π L,T R = 0,G +−,µν andG −+,µν reduce to the form given in Eqs. (3.15)-(3.16).

Self-energy at high temperatures
The remaining task towards determining the resummed photon propagator is to compute the transverse and longitudinal functions Π L,T R (q) from the retarded photon self-energy diagram, as shown in Fig. 4. For a given ψ in the photon self-energy loop, the amplitude is given by where S −+ and S +− are given in Eqs. (3.7) and (3.8) while the time-ordered propagator neglecting the mass of ψ is given by Note that in Eq.
For q 2 < 0, it leads to cos θ 1 : while for q 2 0, it gives We can also see from Eq. (3.28) that if we neglect the q 2 term, the conditions | cos θ 1,2 | 1 are only supported by spacelike photon propagation, i.e., −| q| < q 0 < | q|. This observation is equivalent to the results under the HTL approximation, where only the corrections from the q 2 < 0 regime are considered in the resummed photon propagator. Evaluating the longitudinal part ImΠ L R = q 2 ImΠ 00 R /| q| 2 in the q 2 < 0 and q 2 0 respectively with the full quantum statistics, we obtain where x 0 ≡ q 0 /T, x q ≡ | q|/T , and Li n is the polylogarithm function of order n. Similarly, the transverse part In the above equations, we have dropped the contributions of ImΠ L,T R in the x 0 > x q region. As will be found out in Sec. 3.5, this region is not supported by the integration of photon momentum q.
It should be pointed out that the above results for ImΠ L,T R are exact without the HTL approximation. If we go to the limit (| q| ± q 0 )/2 = 0 in Eqs. (3.29)-(3.30) and take the Fermi-Dirac statistics, ImΠ L,T R at q 2 < 0 would reduce to the known results in the HTL approximation [57], i.e., Since the imaginary parts can be analytically integrated even with the full quantum statistics due to the presence of two Dirac δ-functions, we will use ImΠ L,T R without the HTL approximation. Taking further into account the contributions in the q 2 0 region, the contributions of scattering with both hard-and soft-momentum transfers will be included simultaneously. On the other hand, an analytic integration with the full quantum statistics cannot be obtained for the real part of Π L,T R . Given that the real part ReΠ L,T R (q) ∼ αT 2 plays the role of IR regulator in the soft-momentum transfer q 2 T 2 but only serves as sub-leading corrections to the photon spectral density in the hard-momentum transfer q 2 αT 2 , we can apply the HTL-approximated results [57,67]: for the photon spectral density ρ L,T (q) in the full momentum space.
The computation thus far considers only the contribution of the electron. At temperatures above the QCD phase transition, one needs to include quarks into the charged-fermion loop. In general, for ν R decoupling above a few hundred GeV, all the SM charged fermions and the charged W boson can contribute to the photon dispersion relation in Fig. 2. Their contributions can be taken into account by replacing where Q ψ denotes the electric charge of ψ. Here the summation goes over all charged particles that are relativistic in the thermal bath. All quarks at a sufficiently high temperature contribute a combined factor of (1/3) 2 × 3 × 3 + (2/3) 2 × 3 × 3 = 5 and all charge leptons contribute a factor of 3. For the W boson, we assume its contribution is the same as the electron, though strictly speaking it would require a more dedicated calculation. When the temperature is not sufficiently high, we treat c ψ as a step function of T , including only charged particles with masses below T .

Full result of the collision rate
Combining results in the above calculations, we can write the imaginary retarded amplitude as Note that due to the presence of δ(k 2 ), which dictates k 2 = (p+q) 2 = 2(q 0 p 0 − p· q)+q 2 = 0, we obtain For q 2 > 0 (i.e. q 0 > | q| or q 0 < −| q|), the second inequality in Eq. (3.41) would not be satisfied if we take the positive branch q 0 > | q| > 0. In this case, we find that k 0 is always negative. Similarly, for q 2 < 0 we find that k 0 is always positive. Hence the sign function takes Assembling all the pieces in the previous calculation, the trace in Eq. (3.3) is given by where the Heaviside functions are defined as Finally, we obtain the master integral for the full collision rate as follows: where , and c ψ is the enhancement factor defined in Eq. (3.39). The three-dimensional integral can be integrated numerically, which yields 7 where we have used c ψ = 9 to include the contribution of all charged particles. If only ψ = e is included (i.e. c ψ = 1), we would have σvn ψ=e ≈ 1.53αµ 2 ν T 3 . If we use the HTL-approximate results given in Eq. (3.36), we obtain σvn ≈ 1.84αµ 2 ν T 3 and Γ tot ≡ Thermally averaged ν L → ν R rate σvn /αµ 2 ν T 3 Zero-T QFT + m γ cut + e ± plasma + Boltzmann statistics 2.26 Finite-T QFT + e ± plasma + HTL approximation 1.84 Finite-T QFT + e ± plasma + full quantum statistics 1.53 Finite-T QFT + ψψ plasma + full quantum statistics 6.47 Table 1. The thermally averaged chirality-flipping rates computed in different methods. "Zero-T QFT + e ± plasma + Boltzmann statistics" denotes the calculation in Sec. 2, where we use Feynman rules of zero-temperature QFT with a finite photon mass cut m γ = eT / √ 6 and the Boltzmann statistics in the relativistic electron-positron (e ± ) plasma. In the finite-T QFT approach from Sec. 3, finite-temperature effects and full quantum statistics (Fermi-Dirac/Bose-Einstein) are taken into account.
The latter should be compared with the rate Γ tot ≈ 1.81αµ 2 ν T 3 obtained in Ref. [35]. Note that the difference between σvn , which can be written as (1 − f ν R )Γ tot according to the definition, and Γ tot is that the former is slightly suppressed by the Pauli-blocking factor (1 − f ν R ). From Eq. (2.12), we can see that the result derived from the tree-level scattering amplitude with a photon thermal mass cut m γ = eT / √ 6 is overestimated. The comparison of our final results obtained in the finite-temperature approach (together with quantum statistics) to the result obtained in Sec. 2 is presented in Tab. 1.

Plasmon decay
At the end of this section, we would like to present an estimate of the plasmon decay (γ * →ν L + ν R ) effect, which would be important when the charged fermions become non-relativistic while their number densities remain high, such as in some stellar environments [29]. The contribution of plasmon decay is already included in the master formula, Eq. (3.46), in which it corresponds to q 2 = Re Π L,T R (q). For a simple estimate, the decay rate can be calculated by using the S-matrix formalism. With the effective Lagrangian Eq. (2.1), it is straightforward to obtain the decay width as follows: where q 2 = E 2 γ * − | q| 2 . The thermally averaged rate can be written as In the non-relativistic regime (T m ψ ) and the HTL approximation, we obtain from Eq. (3.26) the photon thermal mass which depends on the photon polarization: where n ψ denotes the number density of ψ. Here we would like to make a comparison to the results in Ref. [68] where the photon dispersion relation for T m e has also been calculated in the stellar plasma. The main difference in the stellar plasma is that, as an electrically neutral medium, its positively charged particles are protons (including protons in nuclei), whose contributions to the photon thermal masses are negligible due to the heavy proton/nucleus mass. In the thermal plasma of the early universe, we have equally high densities of both electrons and positrons. Therefore, the photon thermal masses in Eq. (3.50) are twice as large as the results for the stellar plasma.
Taking the approximate dispersion relation q 2 ≈ m 2 γ with m 2 γ ≈ 2e 2 n ψ /m ψ (i.e. only the transverse polarization mode is used) and the nonrelativistic limit of n ψ , we obtain the thermally averaged production rate We can see that the rate from plasmon decay is at O(α 2 ) and is exponentially suppressed when T m ψ . In the relativistic regime, on the other hand, the calculation is similar and we find Overall, the contribution of plasmon decay to ν R production in the early Universe is subdominant because the results in Eqs. (3.52) and (3.53) are proportional to α 2 .

NMM bounds from N eff constraints
With the ν R production rate obtained (see Tab. 1), we are ready to relate NMM to N eff and use cosmological measurements of N eff to constrain NMM.
As we have discussed in Sec. 2, due to σvn ∝ T 3 and H ∝ T 2 , we expect that ν R is in thermal equilibrium at sufficiently high temperatures and decouples from the thermal bath at low temperatures. The decoupling temperature, T dec , is determined by solving H = σvn , i.e., 1.66 In general, as g and σvn /T 3 both vary with the temperature, one has to solve Eq. The full numerical relation between T dec and µ ν is shown in Fig. 5, where the temperature dependence of g (T ) is taken from Ref. [69] and c ψ takes the step function as mentioned below Eq. (3.39). For comparison, the approximate relation in Eq. (4.2) is also shown in Fig. 5. The contribution of ν R to the effective neutrino number excess, ∆N eff , depends on T dec , and the relation can be derived from entropy conservation-see, e.g., [50,51,56]. Given T dec , ∆N eff is determined by where N R is the number of thermalized light ν R . Although N R = 3 is the most natural assumption because for three Dirac neutrinos with NMM all ν R 's are expected to thermalize at a sufficiently high T , it is still possible to have N R = 2 or 1. For instance, it could be that not all neutrinos are Dirac so that there are only one or two light ν R species present. Even if all neutrinos are Dirac, the actual temperature required by H = σvn could be too high so that the universe has never reached such a high temperature (e.g. if the temperature is above the reheating temperature after inflation). Therefore, we keep the possibilities of N R = 2 and 1 in our discussion below.

Results
Applying the numerical T dec (µ ν ) relation in Fig. 5 to Eq. (4.3), we plot ∆N eff as a function of µ ν for N R = 1, 2, 3 in Fig. 6, together with cosmological bounds on ∆N eff and other known bounds on µ ν . The bounds are explained as follows. Currently, the best measurements of N eff come from combinations of CMB, BAO, and BBN observations. We adopt results of the CMB+BAO combination from the Planck 2018  Figure 6. The dependence of ∆N eff on µ ν , assuming the number of ν R thermalized via NMM is N R = 1, 2, or 3. The red horizontal line represents the current best limit on ∆N eff from the combination of latest CMB (Planck 2018) and BBN data. For N R = 3, it can be recast to the corresponding bound on µ ν , as indicated by the red vertical line. Future experiments such as SO/SPT-3G and CMB-S4 will be able to probe the N R = 2 and N R = 1 scenarios. Other constraints are explained in the text. publication [70] and the CMB+BBN combination from Ref. [71]   Subtracting the standard value, N st. eff = 3.045 [43,72,73] 9 , we obtain the upper limits ∆N eff < 0.285 and ∆N eff < 0.163 at 95% (2σ) C.L., respectively.
Future CMB experiments such as the Simons Observatory (SO) and CMB Stage-IV (CMB-S4) can reach ∆N eff < 0.1 [40,77] and ∆N eff < 0.06 [38,78] at 2σ C.L., respectively. They are plotted in Fig. 6 as dashed lines. The South Pole Telescope (SPT-3G) [39] has approximately the same sensitivity reach as SO. So we refer to their common limit as SO/SPT-3G.
Low-energy elastic ν + e − scattering data in neutrino and dark matter detectors can be used to constrain NMM. In Tab. 2, we list such bounds from the very recent XENONnT [30] and LUX-ZEPLIN (LZ) [31] experiments as well as that from Borexino [16] which used to be the strongest laboratory constraint. Note that these bounds are all derived from solar neutrinos which change flavors when arriving at the Earth. At low energies, the probabilities of solar neutrinos appearing as ν e , ν µ , and ν τ are around 56%, 22%, and 22%,

Upper Limit
Ref.
respectively [79]. Here we assume flavor-universal and flavor-diagonal NMM, for which the reported bounds in these experiments can be directly compared to our results. If one focuses on a NMM of a specific flavor, then the above probabilities should be taken into account and would lead to weaker bounds. We refer to Ref. [80] for a more dedicated treatment of this issue.
In addition to laboratory bounds, we also include the astrophysical bounds derived from the red-giant branch, µ ν < 2.2 × 10 −12 µ B [32] and µ ν < 1.5 × 10 −12 µ B [33]. In Fig. 6, we only plot the latter to avoid cluttering. Due to potential astrophysical uncertainties, this bound is presented as a dash-dotted line.
As can be read from Fig. 6, if three ν R 's are thermalized via NMM in the early Universe the current BBN+CMB data would set an upper limit on µ ν with BBN+CMB : µ ν < 2.7 × 10 −12 µ B (for N R = 3) . (4.6) It should be noted that for N R = 3, the minimal value of ∆N eff is 0.14, as can be seen from the N R = 3 curve at lower µ ν values in Fig. 6. Future experiments like SO/SPT-3G and CMB-S4 are sensitive to ∆N eff below this level. If ∆N eff is found to be smaller than 0.14 in future measurements, the N R = 3 scenario (i.e. three ν R 's being thermalized via NMM) would be ruled out. There are multiple possibilities to go beyond this scenario: (i) the effective NMM vertex opens up at high temperatures, leading to the transition from freeze-out to freeze-in (see e.g. [51,54]); (ii) as T increases above the electroweak scale, g (T ) continues increasing due to new particles beyond the SM; (iii) the number of ν R that can be possibly thermalized via NMM is less than three. In either (i) or (ii), the curves in Fig. 6 will further decrease as µ ν decreases, but this part would be quite model dependent. We will investigate this in future work. As for (iii), one may consider N R = 2 or 1. For N R = 2, the SO/SPT-3G sensitivity is able to reach µ ν < 2.0 × 10 −12 µ B while for N R = 1 the future CMB-S4 sensitivity can set µ ν < 3.7 × 10 −12 µ B . We see from Fig. 6 that all these cosmological bounds are stronger than the laboratory ones from XENONnT, LZ, and Borexino, and comparable to astrophysical ones. Note that the strongest astrophysical bound cuts the N R = 3 curve at ∆N eff = 0.143 while the curve becomes flat at 0.141. So there is the possibility for upcoming precision measurements of N eff to make a discovery. If a small ∆N eff between 0.143 and 0.141 is probed, it can be attributed to NMM smaller than the strongest astrophysical bounds. What would be more interesting is that future experiments could measure an even smaller ∆N eff below 0.141. In this case, it could probe the freeze-in regime which is sensitive to how the effective vertex opens up at high energies. We leave this possibility for future exploration.

Discussions
We would like to make a few comments on the underlying assumptions of our results. Large NMM arising from new physics models typically involve new energy scales (denoted by Λ), as well as new heavy particles that directly couple to ν R . At temperatures around or above Λ, the effective NMM vertex might be invalid in the calculations of collision terms and the new particles could contribute to g significantly. So our results rely on the assumption that the new particle masses and Λ are well above the temperatures relevant to our calculations. More specifically, for the bounds presented in Fig. 6, we are only concerned with ν R decoupling at temperatures around or below the electroweak scale. Therefore, it is reasonable to assume that the new physics scales are well above the decoupling temperature so that our results are independent of the UV theories of NMM. Let us consider NMM generated at the one-loop level, for which the magnitude of µ ν is roughly given by [12] µ ν ∼ m X 16π 2 Λ 2 , (4.7) where 16π 2 is the one-loop suppression factor, m X denotes the chirality-flipping fermion mass, and Λ corresponds to the heaviest particle mass in the loop. Taking the SM prediction µ ν = 3eG F m ν /(8 √ 2π 2 ) as an example, Λ corresponds to the electroweak scale G −1/2 F and m X corresponds to m ν because the chirality has to flip only via the neutrino mass termchirality-flipping via a charged-lepton mass is not feasible because ν R does not participate in the SM gauge interactions. In left-right symmetric models, due to the presence of gauge interactions with ν R , chirality-flipping can be achieved by the charged fermion running in the loop [4]. In this case, we would have m X = m with = e, µ, or τ , leading to a great enhancement of the NMM. Taking µ ν ∼ 10 −12 µ B in Eq. (4.7), we have Λ ∼ 4.6 TeV · m X /GeV, which implies that for m X = m τ the new physics scale would be well above TeV. For new heavy fermions playing the role of chirality flipping, Λ could be even higher.

Conclusion
In this paper, we investigate cosmological constraints on Dirac NMM from current and future precision measurements of the effective number of neutrinos, N eff . As has been realized in previous studies, a straightforward calculation of the ν R production rate using tree-level scattering amplitudes is IR divergent. Therefore, in this work we compute the rate in a thermal QFT approach which, by taking into account the effects of modified photon dispersion relations in the ν R self-energy, is free from the IR divergence. Furthermore, this approach automatically includes the contributions of the s-and t-channel scattering processes and plasmon decay.
Using the refined ν R production rate, we obtain accurate relations between µ ν and N eff , in which the precision measurements of the latter can be recast to constrain the former. The main results we have obtained are presented in Fig. 6 and summarized in Tab. 2. For three ν R being thermalized via flavor-universal NMM, the combination of current CMB and BBN measurements of N eff puts a strong bound, µ ν < 2.7 × 10 −12 µ B at 2σ C.L. This is better than the latest laboratory bounds from XENONnT and LZ. Future measurements from SO/SPT-3G and CMB-S4 will even be able to exclude scenarios with three or two ν R being thermalized. Our results are applicable to a broad class of NMM models in which the new physic scale is much higher than the electroweak scale.

A Weldon's formula with chiral fermions.
In this appendix, we briefly review Weldon's formula [64] which relates the gain/loss rate of a particle in the thermal bath to its self-energy computed in finite-temperature field theories, as we have applied in Eq. (3.2). Although the formula has been elaborated in great detail in the original paper, when applying to chiral fermions, there could be a subtle issue regarding factors of two. Hence we would like to take this appendix to clear potential issues.
Let us first consider a non-chiral toy model: where ψ 1,2 are two Dirac fermions and φ is a real scalar. For simplicity, we ignore the masses of ψ 1 and φ, and concentrate on the production and decay of ψ 2 in the thermal bath. In this context, Weldon's formula (see Eq. (2.31) in Ref. [64]) reads where u 2 (u 2 ) denotes the initial (final) state of the ψ 2 particle, Σ denotes its self-energy shown in Fig. 7, E p and p are the particle energy and momentum, and Γ tot ≡ Γ gain + Γ loss with Γ gain and Γ loss defined similar to those in Eq.
where the momenta p, k, q have been specified in Fig. 7, S ψ 1 and D φ denote the propagators of ψ 1 and φ at a finite temperature, T . Following Ref. [64], here we use the imaginary-time formalism, in which the energy take discrete values i2πT b where b is an integer for a boson or a half-integer for a fermion. For this reason, the energy integral dk 0 /(2π) has been discretized to the summation of b. After computing the summation, one can get where "· · · " represents other terms proportional to δ( We refer to Eq. (2.22) in Ref. [64] for their explicit forms. Now let us turn to the right-hand side of Eq. (A.2). Recall that the collision term in the Boltzmann equation for f ψ 2 is given by where s 1 denotes the spin of ψ 1 and dΠ x ≡ d 3 p x / (2π) 3 2E x does not include internal degrees of freedom. Instead, we have s 1 to sum over the spin of ψ 1 explicitly. The squared amplitudes read After the spin summation, we have which is identical to the part between δ( where in the second step we have integrated out dΠ φ with a part of the delta function. By comparing Eq. (A.8) with Eq. (A.4), we see that Weldon's formula is explicitly verified, except for the "· · · " part in Eq. (A.4) accounting for contributions of other processes.
Note that so far we have not summed the spin of ψ 2 , nor have we included its internal degrees of freedom in the calculation. Hence both sides of Eq. (A.2) should be interpreted as quantities for only one of its degrees of freedom (if there are many). For the non-chiral model considered above, since the results for the two possible spin polarizations of a ψ 2 particle are equal, one can average over the spin polarizations: Now let us revisit the calculation for the following chiral toy model: i.e., we assume that only the left-handed part of ψ 2 participates in the yukawa interaction.
For this model, one can still explicitly verify Weldon's formula with a few modifications as follows. First, the two vertices in Fig. 7 are accompanied with the chiral projectors P L and P R . So we have S ψ 1 (k) → P R S ψ 1 (k)P L in Eq. (A.3) and u 2 / ku 2 → u 2 P R / kP L u 2 in Eq. (A.4). The squared amplitude is also modified: |M| 2 → y 2 |u 1 P L u 2 | 2 = y 2 |u 2 P R u 1 | 2 . (A.11) Here the projectors automatically guarantee that only left-handed ψ 2 and right-handed ψ 1 are involved in the amplitude. So one can still sum over the spin because the contribution of the wrong spin polarization automatically vanish due to P L/R . With the above details being noted, we can see that Eq. (A.2) for chiral fermions still holds. However, due to the absence of right-handed ψ 2 , Eq. (A.9) should be modified as where in principle the projectors P L/R should be included in Σ but here we prefer to write them out explicitly. For our application to neutrinos in Sec. 3, the projectors have been included implicitly in Eq. (3.2).