#  Hubble Tension and Cosmological Imprints of $U(1)_{X}$ Gauge Symmetry:
$U(1)_{B_{3}-3L_{i}}$ as a case study

Dilip Kumar Ghosh<EMAIL_ADDRESS>School of Physical Sciences, Indian
Association for Cultivation of Science,
2A $\&$ 2B Raja S C Mullick Road, Kolkata 700032, India  Purusottam Ghosh
<EMAIL_ADDRESS>School of Physical Sciences, Indian Association for
Cultivation of Science,
2A $\&$ 2B Raja S C Mullick Road, Kolkata 700032, India  Sk Jeesun
<EMAIL_ADDRESS>School of Physical Sciences, Indian Association for
Cultivation of Science,
2A $\&$ 2B Raja S C Mullick Road, Kolkata 700032, India  Rahul Srivastava
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Science
Education and Research - Bhopal
Bhopal Bypass Road, Bhauri, Bhopal, India

###### Abstract

The current upper limit on $N_{\rm eff}$ at the time of CMB by Planck 2018 can
place stringent constraints in the parameter space of BSM paradigms where
their additional interactions may affect neutrino decoupling. Motivated by
this fact in this paper we explore the consequences of light gauge boson
($Z^{\prime}$) emerging from local $U(1)_{X}$ symmetry in $N_{\rm eff}$ at the
time of CMB. First, we analyze the generic $U(1)_{X}$ models with arbitrary
charge assignments for the SM fermions and show that, in the context of
$N_{\rm eff}$ the generic $U(1)_{X}$ gauged models can be broadly classified
into two categories, depending on the charge assignments of first generation
leptons. We then perform a detailed analysis with two specific $U(1)_{X}$
models: $U(1)_{B_{3}-3L_{e}}$ and $U(1)_{B_{3}-3L_{\mu}}$ and explore the
contribution in $N_{\rm eff}$ due to the presence of $Z^{\prime}$ realized in
those models. For comparison, we also showcase the constraints from low energy
experiments like: Borexino, Xenon 1T, neutrino trident, etc. We show that in a
specific parameter space, particularly in the low mass region of $Z^{\prime}$,
the bound from $N_{\rm eff}$ (Planck 2018) is more stringent than the
experimental constraints. Additionally, a part of the regions of the same
parameter space may also relax the $H_{0}$ tension.

##  1 INTRODUCTION

The cosmological parameter $N_{\rm eff}$, associated with the number of
relativistic degrees of freedom, is crucial in describing the dynamics of the
thermal history of the early universe. At very high temperature of the
universe, the photon and electron bath were coupled , whereas at low
temperature the interaction rate drops below the Hubble expansion rate and the
two baths decouple [1]. $N_{\rm eff}$ is parameterised in terms of the ratio
of energy densities of photon and neutrino bath. Within the Standard Model
(SM) particle contents, the two aforementioned baths were coupled through weak
interactions at high temperature and as temperature drops ($T\sim 2$ MeV) they
decouple. Assuming such scenario the predicted value of $N_{\rm eff}^{\rm SM}$
turns out to be $3.046$ [2, 3]. This value of $N_{\rm eff}$ deviates from the
number of neutrinos ($3$) in the SM particle content is due to various non-
trivial effects like non-instantaneous neutrino decoupling, finite temperature
QED corrections and flavour oscillations of neutrinos [2, 3]. However, in
observational cosmology also, measurements of the cosmic microwave background
(CMB), baryon acoustic oscillations (BAO), and other cosmological probes
provide constraints on $N_{\rm eff}$. The current Planck 2018 data has precise
measurement of $N_{\rm eff}$ at the time of CMB with $95\%$ confidence level,
$N_{\rm eff}=2.99^{+0.34}_{-0.33}$ [4]. Thus the current upper limit from
Planck 2018 data shows that there can be additional contribution (apart from
SM predicted value) to $N_{\rm eff}$, indicating the scope for new physics.

It is evident that $N_{\rm eff}$ will change in the presence of any beyond
standard model (BSM) particles with sufficient interactions with either of the
photon or neutrino bath at temperatures relevant for neutrino decoupling [5,
6, 7, 8] or in presence of any extra radiation [9, 10, 11, 12]. Thus the upper
limit on $N_{\rm eff}$ from CMB can be used to constrain such BSM paradigms
dealing with any extra energy injection. Several studies have been performed
to explore the imprints of BSM models in $N_{\rm eff}$ like models with early
dark energy [10, 13], relativistic decaying dark matter [14] and non-standard
neutrino interactions (NSI) [5, 6, 15, 16, 17, 18]. From the perspective of
neutrino physics, the last one is very interesting since the non-standard
interactions of light neutrinos may alter the late-time dynamics between
photon and neutrino baths, contributing to $N_{\rm eff}$ and the same NSI
interactions can be probed from ground based neutrino experiments as well [19,
20].

On the other hand, the anomaly-free $U(1)_{X}$ gauge extended BSM models are
well-motivated from several aspects like non zero neutrino masses [21], flavor
anomalies [22] etc. The anomaly condition allows the introduction of right-
handed neutrinos in the theory and thus it can also explain non-zero neutrino
masses via the Type-I seesaw mechanism [21]. These scenarios naturally involve
a gauge boson (${Z^{\prime}}$) that originates from the $U(1)_{X}$ abelian
gauge symmetry and has neutral current interactions with neutrinos and
electrons which may have some nontrivial role in neutrino decoupling and hence
in deciding $N_{\rm eff}$ [6, 23, 7, 15]. The detectability of gauge boson
throughout the mass scale ($M_{Z^{\prime}}$) also motivates such scenarios.
For TeV-scale $Z^{\prime}$, constraints arise from collider experiments [24,
25, 26, 27], whereas in the sub-GeV mass region, low energy scattering
experiments (neutrino electron scattering [28], neutrino-nucleus scattering
[29] etc.) are relevant to constrain the parameter space. However, both types
of direct searches become less sensitive in the mass
$M_{Z^{\prime}}\lesssim\mathcal{O}(10~{}{\rm MeV})$ and lower coupling
($g_{X}\lesssim 10^{-5}$) region. In that case, the CMB observation on $N_{\rm
eff}$ plays a crucial role and can impose severe constraints on the
$M_{Z^{\prime}}-g_{X}$ plane, which stands as the primary focus of our study.

In cosmology, there exists some discrepancy between the values of expansion
rate $H_{0}$ obtained from CMB and local measurement [30, 31, 32, 33, 4].
Using direct observations of celestial body distances and velocities, the
SHOES collaboration calculated $H_{0}=73.04\pm 1.04~{}{\rm
Km~{}s^{-1}~{}Mpc^{-1}}$ [33]. However, the CMB measurements like Planck 2015
TT data predict $H_{0}=68.0^{+2.6}_{-3.0}~{}{\rm Km~{}s}^{-1}\rm{Mpc}^{-1}$ at
$1\sigma$ [34] and from recent Planck 2018 collaboration it turns out to be
$H_{0}=67.36\pm 0.54{~{}\rm Km~{}s}^{-1}~{}\rm{Mpc}^{-1}$ [4], where both the
collaboration analyzes the CMB data under the assumption of $\Lambda$CDM
cosmology. So there is a disagreement between the local and CMB measurements
roughly at the level of $4\sigma-6\sigma$ [30]. Although such discrepancy may
arise from systematic error in measurements [35, 36], it also provides a hint
for BSM scenarios affecting the dynamics of the early universe. One possible
way out to relax this so called “$H_{0}$ tension” involves increasing $N_{\rm
eff}$ in the approximate range $\sim 3.2-3.5$ [37, 6, 38]111However, it is
essential to acknowledge that increasing $N_{\rm eff}$ also provokes the
tension related to the other parameter of $\Lambda$CDM model, $\sigma_{8}$ [4,
39].. As discussed earlier $U(1)_{X}$ local gauge extension leads to new gauge
boson and new interaction with SM leptons which might affect neutrino
decoupling and hence contribute to $N_{\rm eff}$. Thus $U(1)_{X}$ gauge
extension can be one possible resolution to the $H_{0}$ tension also. Two
kinds of $U(1)_{X}$ scenarios have been considered so far in the literature to
address the Hubble tension problem and/or the excess in $N_{\rm eff}$ at CMB:
$U(1)_{\mu-\tau}$ [6] and $U(1)_{B-L}$ [7]. In both these two cases, the
corresponding gauge boson mass lies below 100 MeV, which attracts severe
constraints from scattering experiments [19, 20], energy loss in supernovae,
etc [40]. Note that in the second scenario, $U(1)_{\mu-\tau}$, the light gauge
boson couples to the $\mu$ and $\tau$ leptons but does not to the electron and
quarks at the tree level. As a result, it is subject to less existing
constraints than $U(1)_{B-L}$.

In this work, we begin with a generic $U(1)_{X}$ extension and explore the
aforementioned cosmological phenomena more comprehensively. For simplicity, we
assume the light gauge boson was in a thermal bath, and the right-handed
neutrinos are heavy enough ($\gtrsim 100$ MeV) that they hardly play any role
in neutrino decoupling. We study the dynamics of neutrino decoupling in the
presence of the light gauge boson ($Z^{\prime}$) in a generic $U(1)_{X}$
scenario with arbitrary $U(1)_{X}$ charge assigned to SM fermions. We evaluate
$N_{\rm eff}$ by solving a set of coupled Boltzmann equations that describe
the evolution of light particles i.e. electron, SM neutrino, and light gauge
boson ( $M_{Z^{\prime}}\lesssim\mathcal{O}(10~{}{\rm MeV})$ ). Thus we
identify the region of parameter space that is consistent with Planck 2018
observations. We also illustrate the parameter space in which Hubble tension
can be relaxed where the value of $N_{\rm eff}$ falls within the range of
$3.2$ to $3.5$ [37]. Our analysis show that the depending on the $U(1)_{X}$
charge assignment of the first generation lepton, the generic $U(1)_{X}$
gauged models can be classified in context of $N_{\rm eff}$. We incorporate
these cosmological observations into specific gauged $U(1)_{X}$ models,
$U(1)_{B_{3}-3L_{i}}$ ($i=e,\mu,\tau$). In comparison to others, these
extensions encounter fewer constraints, as they are exclusively linked only to
the third generation of quarks and one generation of leptons. Our results show
that the upper limit on $N_{\rm eff}$ from Planck 2018 [4] can put stringent
bounds on the parameter space of light $Z^{\prime}$ in the mass region
$\lesssim\mathcal{O}(10~{}{\rm MeV})$ where the other experimental bounds are
comparatively relaxed.

The paper is organized as follows. In sec.2, we begin with a brief model-
independent discussion on the generic $U(1)_{X}$ model. Sec.3 is dedicated to
a comprehensive analysis and discussion of the dynamics governing neutrino
decoupling in terms of the cosmological parameter $N_{\rm eff}$ in the
presence of the light $Z^{\prime}$ originated from the generic $U(1)_{X}$.
Subsequently, in sec.4, we present the numerical results in terms of $N_{\rm
eff}$ for this generic $U(1)_{X}$ model. In sec.5, we discuss the analysis
presented earlier, which applies to the specific $U(1)_{X}$ model,
$U(1)_{B_{3}-3L_{i}}(i=e,\mu,\tau)$. We present a brief discussion on $H_{0}$
tension in sec.6. Finally, in sec.7, we summarise and conclude with the
outcomes of our analysis. In Appendices A to C, we provide various technical
details for the calculation of $N_{\rm eff}$.

##  2 Model Independent Discussion: Effective $Z^{\prime}$ Models

In this section, we look at effective light $Z^{\prime}$ models. We take a
model-independent but minimalist approach and only assume that the
$Z^{\prime}$ originates from the breaking of a new local $U(1)_{X}$ abelian
gauge symmetry under which the SM particles are charged as shown in Table 1.
The anomaly cancellation conditions for a typical $U(1)_{X}$ local gauge
symmetry requires the presence of additional chiral fermions beyond the SM
fermions222Note that for one generation of SM fermions, the only $U(1)$
symmetry which is anomaly free is the SM hypercharge $U(1)_{Y}$ symmetry. For
the full three generations of SM fermions, it is possible to have additional
$U(1)_{X}$ symmetries which can be made anomaly free without adding any new
chiral fermions such as $U(1)_{q_{i}-q_{j}}$ ; $i,j=1,2,3$ with $q_{i},q_{j}$
denoting charges of a given generation of SM quarks but such symmetries have
other phenomenological problems and are not of interest to the current study..
To keep the minimal scenario we assume that the only BSM fermions needed for
anomaly cancellation are the (three) right handed neutrinos $\nu_{R}$ which
are taken to be charged under the $U(1)_{X}$ symmetry. With the addition of
$\nu_{R}$, several different type of gauged $U(1)_{X}$ models including the
popular $U(1)_{X}$ models such as the $U(1)_{B-L}$ or $U(1)_{\mu-\tau}$ can be
constructed. We will discuss some of the models in the next Section 5.

For this purpose within this section, we will not specify the nature of
$U(1)_{X}$ symmetry nor the charges of SM particles and $\nu_{R}$ under it and
will treat them as free parameters and proceed with a generic discussion. We
will also not go into details of anomaly cancellation constraints. All these
things will be clarified in the following sections during our discussions on
some well motivated $U(1)_{X}$ models.

Fields |  $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ |  $U(1)_{{X}}$
---|---|---
$Q_{i}$ | $(3,2,\frac{1}{3})$ | $\mathbb{X}_{Q_{i}}$
$u_{i}$ | $(3,1,\frac{4}{3})$ | $\mathbb{X}_{u_{i}}$
$d_{i}$ | $(3,1,-\frac{2}{3})$ | $\mathbb{X}_{d_{i}}$
$L_{i}$ | $(1,2,-1)$ | $\mathbb{X}_{L_{i}}$
$\ell_{i}$ | $(1,1,-2)$ | $\mathbb{X}_{\ell_{i}}$
$\nu_{R_{i}}$ | $(1,1,0)$ | $\mathbb{X}_{\nu_{i}}$
$\Phi$ | $(1,2,1)$ | $\mathbb{X}_{\Phi}$
$\sigma$ | $(1,1,0)$ | $\mathbb{X}_{\sigma}$
Table 1: Particle content and $U(1)_{X}$ gauge charge assignments of Standard
Model and new particles. In this section, we will take $\mathbb{X}_{\Phi}=0$
as we are interested in a light $Z^{\prime}$ gauge boson.

In Table 1, in addition to $\nu_{R}$ (typically required for anomaly
cancellation) we also add an SM singlet scalar $\sigma$ carrying $U(1)_{X}$
charge whose vacuum expectation value (VEV) will break the $U(1)_{X}$
symmetry. Note that apart from the minimal particle content of Table 1, most
of the $U(1)_{X}$ models available in literature may also contain additional
BSM particles333Such new particles are typically much heavier than MeV scale
and will not change our analysis [22].. These additional particles are model
dependent and we refrain from adding them to proceed with a minimal setup.

Now some general model independent simplifications and conclusions can be
immediately drawn for the charges of the particles under $U(1)_{X}$ symmetry
listed in the aforementioned Table 1.

  1. 1.

Light $Z^{\prime}$: In the upcoming sections, we delve into an effective
resolution of the well-known Hubble parameter tension [37] by increasing
$N_{\rm eff}$ considering a light $Z^{\prime}$ with mass
$M_{Z^{\prime}}\sim\mathcal{O}(\rm MeV)$ mass range[6]. To avoid any fine-
tuning in keeping the $M_{Z^{\prime}}$ very light compared to the SM $Z$ gauge
boson mass i.e. $M^{2}_{Z^{\prime}}\ll M^{2}_{Z}$, we consider the
corresponding $U(1)_{X}$ charge of the SM Higgs doublet $\mathbb{X}_{\Phi}=0$,
so that $M^{2}_{Z^{\prime}}$ does not receive any contribution from the SM
vev.

  2. 2.

Mass generation for quarks and charged leptons: For the choice of
$\mathbb{X}_{\Phi}=0$, to generate the SM quark and charged lepton masses, we
must take $\mathbb{X}_{Q_{i}}=\mathbb{X}_{u_{i}}=\mathbb{X}_{d_{i}}$ and
$\mathbb{X}_{L_{i}}=\mathbb{X}_{\ell_{i}}$ such that the standard Yukawa term
$y^{u}_{ij}\bar{Q}_{i}\tilde{\Phi}u_{j}$ involving only SM fields can be
written in the canonical form as shown in eq.(1). It is important to note that
although within each generation, the $U(1)_{X}$ charges of quark doublet
should be the same as that of the up and down quark singlets, the charges may
differ when comparing across different generations i.e.
$\mathbb{X}_{Q_{1}}\neq\mathbb{X}_{Q_{2}}\neq\mathbb{X}_{Q_{3}}$. The same
applies to charged leptons. In fact, in later sections, we will indeed
consider flavour dependent $U(1)_{X}$ symmetries. To simplify our notation
throughout the remaining part of this paper we will denote
$\mathbb{X}_{L_{i}}=\mathbb{X}_{\ell_{i}}=X_{i}$.

  3. 3.

Quark mixing: As a follow up point note that if the charges of all three
generations of quarks are unequal i.e. if we take
$\mathbb{X}_{Q_{1}}\neq\mathbb{X}_{Q_{2}}\neq\mathbb{X}_{Q_{3}}$ then the
Yukawa matrices $Y^{u}$ and $Y^{d}$ and hence the resulting mass matrices will
only have diagonal entries and we will not be able to generate CKM mixing.
Thus, the charges of some (but not all) generations should match with each
other in order to allow the generation of quark mixing. The same is true for
charged lepton mass matrices. Again we will elaborate it further with specific
examples in the coming sections.

  4. 4.

BSM fermions: The $U(1)_{X}$ charge of right handed neutrinos is typically
fixed by anomaly cancellation conditions as we will discuss in a later section
with specific examples. As mentioned earlier, for the sake of simplicity, we
will not consider $U(1)_{X}$ symmetries involving chiral fermions beyond the
fermion content of Table 1.

Based on the aforementioned assumptions one can write the Yukawa and scalar
potential for the general $U(1)_{X}$ model as:

$\displaystyle\mathcal{L}_{Yuk}\,=\,y^{u}_{ij}\bar{Q}_{i}\tilde{\Phi}u_{j}\,+\,y^{d}_{ij}\bar{Q}_{i}\Phi
d_{j}\,+\,y^{e}_{ij}\bar{L}_{i}\Phi\ell_{j}\,+\,\mathcal{L}_{\nu}+h.c.$ (1)

and

$\displaystyle
V(\Phi,\sigma)\,=\,\mu^{2}_{\Phi}\Phi^{\dagger}\Phi\,+\,\mu^{2}_{\sigma}\sigma^{\dagger}\sigma\,+\,\lambda_{\Phi}\Phi^{\dagger}\Phi\Phi^{\dagger}\Phi\,+\,\lambda_{\sigma}\sigma^{\dagger}\sigma\sigma^{\dagger}\sigma\,+\,\lambda_{\Phi\sigma}\Phi^{\dagger}\Phi\sigma^{\dagger}\sigma~{}~{}.$
(2)

It is worth highlighting a couple of salient features of the Yukawa terms and
the scalar potential of this scenario.

  1. 1.

In eq.(1), $\mathcal{L}_{\nu}$ refers to the terms needed for light active
neutrino mass generation. These terms depend on the details of the $U(1)_{X}$
symmetry, the charges of leptons under it as well as the nature of neutrinos
(Dirac or Majorana) and the mechanism involved for mass generation.
Furthermore, this typically requires the presence of additional scalars or
fermions or both, beyond the particle content listed in Tab. 1. Note that even
in the massless limit and absence of any other interactions, the $\nu_{R}$ is
still interacting with the rest of the particles through its $U(1)_{X}$ gauge
interactions and depending on the charges and strength of the $U(1)_{X}$ gauge
coupling ($g_{X}$), they can be in thermal equilibrium with the rest of the
plasma at a given epoch in the evolution phase of the early universe.

  2. 2.

Since we want a very light $Z^{\prime}$, therefore the vev of $\sigma$ field
$\langle\sigma\rangle=v_{\sigma}$ which is responsible for $Z^{\prime}$ mass
generation should be small. Furthermore, after SSB, if we want the real
physical scalar $\sigma_{R}$ to be light as well, we should have
$\lambda_{\Phi\sigma}\ll 1$. The condition $\lambda_{\Phi\sigma}\ll 1$ will
also imply that the 125 GeV scalar ($h$) is primarily composed of the real
part of $\Phi^{0}$ (neutral component of $\Phi$) and hence the LHC constraints
on it can be trivially satisfied.

With these general assumptions, one can examine the potential range of a
$Z^{\prime}$ parameter that may account for the observed excess of $N_{\rm
eff}$ in the Planck 2018 data, leaving the charge of the leptons and the gauge
coupling as free parameters.

##  3 $\nu_{L}$ decoupling in presence of light $Z^{\prime}$ and ${N_{\rm
eff}}$

It is a well known fact that one of the most precisely measured quantities
from cosmology is the effective neutrino degrees of freedom ($N_{\rm eff}$) ,
which may get altered by the presence of light BSM particles. This change in
$N_{\rm eff}$ from the SM value can in principle provide one of the solutions
to relax the Hubble tension [37]. In this work, we focus on a scenario with
light $Z^{\prime}$ interacting with SM neutrinos ($\nu_{i}~{},(i\equiv
e,\mu,\tau)$) that may lead to change in $N_{\rm{eff}}$. After having a brief
description of generic features of light $Z^{\prime}$ models emerging from
$U(1)_{X}$ symmetries in the previous section, we now move to the most crucial
part of our paper which is the cosmological implications of those scenarios.
Before delving into the analysis of $N_{\rm eff}$ in the presence of such
light $Z^{\prime}$, we would like to mention some key aspects of $N_{\rm eff}$
within the SM framework.

In a standard cosmological scenario at temperature $T\sim 20$ MeV 444 By the
temperature of the universe ($T$), we mean it as the temperature
($T_{\gamma}$) of the photon bath. , only $e^{\pm},\nu_{i}$ are particles
coupled to thermal (photon) bath as the energy densities of other heavier SM
particles are already Boltzmann suppressed. As the universe cools, once the
weak interactions involving $e^{\pm}$ and $\nu_{i}$ drop below the Hubble
expansion rate $H(T)$, neutrinos decouple from the photon bath. Considering
only the SM weak interactions, the neutrino decoupling555by $\nu_{L}$ we refer
the whole SM neutrinos ($\nu_{i}$) temperature turns out to be $2$ MeV [41,
1]. After that, there exist two separate baths of photon ($+e$) and $\nu_{i}$,
each with different temperatures; $T_{\gamma}$ and $T_{\nu}$ respectively.
Approximately, at a temperature below $T_{\gamma}\lesssim 0.5$ MeV, $e^{\pm}$
annihilate, and the entropy is transferred entirely to photon bath leading to
a increment in $T_{\gamma}$ compared to neutrino bath. This difference in
temperature is parameterised in terms of $N_{\rm eff}$ which is given by[1],

$\displaystyle N_{\rm
eff}=\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\left(\frac{\rho_{\nu_{L}}(T_{\nu})}{\rho_{\gamma}(T_{\gamma})}\right)=3\times\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\left(\frac{T_{\nu}^{4}}{T_{\gamma}^{4}}\right)$
(3)

At the time of CMB formation, the predicted value of $N_{\rm eff}$ within the
SM particle content is $N_{\rm eff}^{\rm CMB}=3.046$ [2], whereas the recent
Planck 2018 data [4] estimates it to be $N_{\rm eff}^{\rm
CMB}=2.99^{+0.33}_{-0.34}$ at $95\%$ confidence level (C.L.).

(a) (b)

Figure 1: Cartoon diagram of (a) particle spectrum and (b) interactions
between three sectors. The vertical axis in (a) denotes the mass scale (or,
temperature) and the dotted line indicates the temperature where BBN started.
The blue lines in (a) signify the interactions between particles that affect
$N_{\rm eff}$. Here $\nu$ in the figure signify SM neutrinos ($\nu_{i}$) only.

In this scenario, the $Z^{\prime}$ arising from the aforementioned $U(1)_{X}$
models introduces new interactions with both $\nu_{i}$ and $e^{\pm}$, can
potentially impact the neutrino decoupling consequently altering $N_{\rm
eff}^{\rm CMB}$. As pointed out previously, the SM particles relevant for
$\nu_{L}$ decoupling are only $e^{\pm}$ and $\nu_{i}$, whereas the energy
densities of heavy SM particles are negligible due to Boltzmann suppression.
The light quarks also do not take part due to the QCD confinement at a much
higher temperature around $\sim 150$ MeV. Following the same argument in BSM
$U(1)_{X}$ scenario, $Z^{\prime}$ must be light enough
($M_{Z^{\prime}}\lesssim\mathcal{O}(10~{}{\rm MeV})$) to affect $\nu_{L}$
decoupling which will be shown in the later part of this section. There are
two other BSM particles present in our model: the BSM scalar $\sigma$ and the
RHN $\nu_{R}$. As we elaborated in sec.2, the BSM scalar $\sigma$ can be taken
as heavy enough so that they are irrelevant for phenomenology at the MeV scale
temperature. Hence, we integrate out $\sigma$ to perceive the sole effect of
light $Z^{\prime}$ in late time cosmology. And in the Majorana type mass
models, $\nu_{R}$ are also too heavy to affect $\nu_{L}$ decoupling [42, 43]
and we can neglect them also. However, in Dirac-type mass models, $\nu_{R}$
are relativistic at MeV temperature and can significantly alter $N_{\rm eff}$
[9]. In this work, we only consider heavy Majorana RHN and ignore their
contribution to temperature evaluation. In Fig.1(a) we present the relevant
mass scales by a schematic diagram.

So, in our proposed scenario, we have to trace the interactions between only
three baths i.e. $e,\nu_{L}$ and $Z^{\prime}$ to evaluate $T_{\nu}$ or $N_{\rm
eff}$ (eq.(3)). We describe the scenario using a cartoon diagram in Fig.1(b).
It is essential to take care of various energy transfers among these $3$
particles as they will play a key role in computing temperature evolution
equations (see Appendix A). The energy transfer rates are dictated by the
various collision processes and the distribution functions of respective
particles [5, 44, 45]. Here we enlist the relevant processes to consider for
the successful evaluation of $N_{\rm eff}$.

  1. 1.

SM contributions: The SM weak interactions are active at temperature $T\sim$
MeV. At this point, active neutrino annihilations
($\nu_{i}\bar{\nu}_{i}\leftrightarrow e^{+}e^{-}$) as well as elastic
scatterings ($\nu_{i}e^{\pm}\leftrightarrow\nu_{i}e^{\pm}$) mediated by SM $Z$
or $W^{\pm}$ take place to maintain the required thermal equilibrium of the
early universe.

  2. 2.

BSM contributions to $\gamma$ bath: For a light $Z^{\prime}$ with mass
$M_{Z^{\prime}}\sim\mathcal{O}$(MeV) sufficient energy density can be pumped
into the thermal bath via the decay and inverse decay between $Z^{\prime}$ and
electrons ($Z^{\prime}\leftrightarrow e^{+}e^{-}$) at temperature around
${\cal O}({\rm MeV})$. Additional contributions to the thermal bath may in
principle come from scattering processes like
$Z^{\prime}Z^{\prime}\leftrightarrow e^{+}e^{-}$,
$Z^{\prime}\gamma\leftrightarrow e^{+}e^{-}$. However, it turns out that for a
very light $Z^{\prime}$, with mass $M_{Z^{\prime}}\sim{\cal O}$ MeV, its
coupling $g_{X}$ with SM fermions is highly constrained from various
experimental data, $g_{X}\lesssim(10^{-3}-10^{-5})$. For such a small
coupling, the aforementioned decay process of $Z^{\prime}$ significantly
dominates over the scattering processes. Thus we ignore the scattering
contributions in our numerical calculations.

  3. 3.

BSM contributions to $\nu$ bath: Similarly, $Z^{\prime}$ can transfer energy
to $\nu_{i}$ bath through decay and inverse decay
($Z^{\prime}\leftrightarrow\nu_{i}\bar{\nu}_{i})$ processes. Moreover,
($\nu_{i}\bar{\nu}_{i}\leftrightarrow e^{+}e^{-}$) scattering process can play
an important role through one loop coupling of $Z^{\prime}$ with electron.
This one-loop coupling of $Z^{\prime}$ with electrons is responsible for
connecting two separate thermal baths containing $\nu_{i}$ and electrons. This
feature can be seen in certain scenarios of the $U(1)_{X}$ models and we will
elaborate on this issue in great detail in a later section.

  4. 4.

Within $\nu$ bath: In this case, if we assume that different $\nu_{i}$
flavours have different temperatures then
$\nu_{i}\bar{\nu}_{i}\leftrightarrow\nu_{j}\bar{\nu}_{j},~{}(i\neq j)$
mediated by both $Z$ and $Z^{\prime}$ will have significant impact on the
overall $\nu$ thermal bath. We will address this point in detail in the last
part of this section.

To construct the temperature equations and compute the energy transfer rates
we adopt the formalism already developed in ref.[5, 46, 6, 47]. It is worth
highlighting the approximations made in the formalism prescribed in ref.[5]
before the description of temperature equations. Firstly, Maxwell Boltzmann
distributions were considered to characterize the phase space distribution of
all particles in equilibrium, aiding in simplifying the collision term
integral. The use of the Fermi Dirac distribution in the collision term does
not alter the energy transfer rates substantially [48, 16]. Additionally, the
electron’s mass was ignored to simplify the collision terms as non-zero
electron mass would have resulted in a minimal modification of the energy
transfer rate, typically less than a few percent [49]. It is to be noted that
the $\nu_{L}$ masses can be easily neglected as the relevant temperatures for
neutrino decoupling is sufficiently higher than $\nu_{L}$-mass.

After successfully demonstrating all relevant processes and stating the
assumptions, we are now set to construct the temperature evolution equations.
The temperature equations are derived from the Liouville equation for phase
space distribution of particles in a thermal bath (eq.(15)) and the collision
terms take care of the energy transfers among involved particles through the
processes discussed before [50]. Following the detailed calculations of the
temperature evolution for the SM and BSM scenarios as displayed in appendices
A and B, here we quote the final results of aforementioned temperature
evolution equations [5, 6]:

$\displaystyle\frac{dT_{\nu_{L}}}{dt}$ $\displaystyle=$
$\displaystyle-\left(4H\rho_{\nu_{L}}-\frac{\delta\rho_{\nu_{L}\to
e^{\pm}}}{\delta t}+\frac{\delta\rho_{Z^{\prime}\to\nu_{L}}}{\delta
t}\right)\left(\frac{\partial\rho_{\nu_{L}}}{\partial
T_{\nu_{L}}}\right)^{-1}$ (4) $\displaystyle\frac{dT_{Z^{\prime}}}{dt}$
$\displaystyle=$
$\displaystyle-\left(3H\left(\rho_{Z^{\prime}}+P_{Z^{\prime}}\right)-\frac{\delta\rho_{Z^{\prime}\to\nu_{L}}}{\delta
t}-\frac{\delta\rho_{Z^{\prime}\to e^{\pm}}}{\delta
t}\right)\left(\frac{\partial\rho_{Z^{\prime}}}{\partial
T_{Z^{\prime}}}\right)^{-1}$ (5) $\displaystyle\frac{dT_{\gamma}}{dt}$
$\displaystyle=$
$\displaystyle-\left(4H\rho_{\gamma}+3H\left(\rho_{e}+p_{e}\right)+\frac{\delta\rho_{\nu_{L}\to
e^{\pm}}}{\delta t}+\frac{\delta\rho_{Z^{\prime}\to e^{\pm}}}{\delta
t}\right)\left(\frac{\partial\rho_{\gamma}}{\partial
T_{\gamma}}+\frac{\partial\rho_{e}}{\partial T_{\gamma}}\right)^{-1},$ (6)

where, $\rho_{r},P_{r}$ and $T_{r}$ signify the energy density, pressure
density, and temperature of species $r$. Terms like $\frac{\delta\rho_{a\to
b}}{\delta t}$ indicate the energy transfer rate from bath $a$ to $b$ and is
determined by integrating the collision terms (see eq.(16)). The energy
transfer rates are discussed in great detail in Appendix B. Here we assume all
three $\nu_{i}$ generations share the same temperature $T_{\nu_{L}}$ 666as for
our analysis RHNs are not relevant we will often denote $T_{\nu_{L}}$ as
$T_{\nu}$ and $\rho_{\nu_{L}}$ is the summation over the energy densities of
three generations of $\nu_{i}$.

Note that, the equations in eq.(4-6) are dependent on the thermal history of
heavy $Z^{\prime}$ when it remains in equilibrium with thermal bath in the
early universe ($T_{Z^{\prime}}\gtrsim M_{Z^{\prime}}$) through its
interaction with fermions ($f$). $Z^{\prime}$ preserves its thermal
equilibrium via decay, inverse decay ($Z^{\prime}\leftrightarrow f\bar{f}$)
and also through scattering ($f\bar{f}\leftrightarrow
Z^{\prime}Z^{\prime}(\gamma)$) process. Thanks to processes like
$f\bar{f}\leftrightarrow Z^{\prime}Z^{\prime}(\gamma)$ the chemical potentials
($\mu_{i}(T)$) are suppressed and $Z^{\prime}$ remain in chemical equilibrium
with the SM bath. The condition for the thermal equilibrium of $Z^{\prime}$
incorporates a lower bound on $g_{X}(\gtrsim 10^{-9})$ and the lower bound may
shift slightly depending on the specific choice of $U(1)_{X}$ charge and the
value of $M_{Z^{\prime}}$. Alternatively, it is highly plausible that
$Z^{\prime}$ was initially not in thermal equilibrium in the early universe,
but produced from other SM particles via freeze in process. In these
scenarios, it is necessary to solve the coupled equations for $\mu_{i}(T)$ to
compute $N_{\rm eff}$ [7]. This alternate scenario requires a detailed
complementary study which will be reported elsewhere.

As mentioned above we affix to the simplest BSM scenario assuming $Z^{\prime}$
in thermal equilibrium and set the initial condition for the set of equations
eq.(4-6) as $T_{\nu_{L}}=T_{Z^{\prime}}=T_{\gamma}$ at $T_{\gamma}\gtrsim
M_{Z^{\prime}}$ for the reasons already discussed before. In such case one can
further simplify the scenario assuming $T_{Z^{\prime}}=T_{\nu_{L}}$ as
$Z^{\prime}$ remains coupled to $\nu$ bath for longer time than with $\gamma$
bath [6]. This makes the term ${\delta\rho_{Z^{\prime}\to\nu_{L}}}/{\delta
t}=0$ and reduces the three equations in eq.(4)-(6) into two. However, the
equations described in the above format are useful for generic scenarios.

So far, we have assumed that all three $\nu_{L}$ share common temperature
$T_{\nu}$ as mentioned earlier. The approximation is valid as neutrino
oscillations are active around MeV temperature leading to all three
generations of neutrinos equilibriating with each other[51, 41, 52, 53].
However one can also evaluate the temperatures assuming different temperatures
$T_{\nu_{i}},(i=e,\mu,\tau)$ for all three generations and the relevant
equation reads as (see eq.(52)),

$\displaystyle\frac{dT_{\nu_{i}}}{dt}$ $\displaystyle=$
$\displaystyle-\left(4H\rho_{\nu_{i}}-\left(\frac{\delta\rho_{\nu_{i}\to
e}}{\delta t}\right)_{\rm tot}-\sum_{j\neq
i}\left(\frac{\delta\rho_{\nu_{i}\to\nu_{j}}}{\delta t}\right)_{\rm
tot}+\frac{\delta\rho_{Z^{\prime}\to\nu_{i}}}{\delta
t}\right)\left(\frac{\partial\rho_{\nu_{i}}}{\partial
T_{\nu_{i}}}\right)^{-1}~{}.$ (7)

Note that earlier eq.(4) differs from the one described in eq.(7) as the later
one also contains the energy transfer rate between different generations of
$\nu_{i}$. However the value of $N_{\rm eff}$ does not change significantly
($\lesssim 10\%$) if one solves the temperature equation without considering
different $T_{\nu_{i}}$ (see Appendix C) [6]. For the same reason we stick to
the simpler scenario assuming a common temperature of $\nu_{L}$ bath and solve
eq.(4)-(6) for numerical estimation of $N_{\rm eff}$ throughout this paper.

(a) (b) (c) (d)

Figure 2: Evolution of $T_{\gamma}/T_{\nu}$ with photon bath temperature
$T_{\gamma}$ assuming all 3 $\nu_{L}$ share same temperature for different
$U(1)_{X}$ charge combinations. We chose a benchmark parameter value
$M_{Z^{\prime}}=10$ MeV for both the plots. The couplings are taken
$g_{X}=10^{-7}$ for (a) and $g_{X}=10^{-8}$ for (b). The lines corresponding
to different charge combinations are indicated by the plot legends in the
figure. The legends named $N_{\rm eff}$ in the plots refer to $N_{\rm
eff}^{\rm CMB}$.

After the detailed discussion of the basic framework, we are now set to
perform an exhaustive numerical analysis. In Fig.2 we show case the evolution
of temperature ratio ($T_{\gamma}/T_{\nu}$) as well as $\Delta N_{\rm
eff}~{}(\equiv N_{\rm eff}^{\rm CMB}-3.046)$ with $T_{\gamma}$ for different
$U(1)_{X}$ charge combinations. Here we denote $T_{\nu}\equiv T_{\nu_{i}}$
with the assumption all three neutrino share a common temperature. Following
the discussion in the previous paragraphs, we stress that the relevant charges
for $\nu_{L}$ decoupling are $X_{L_{i}}=X_{l_{i}}\equiv X_{i}$ or more
precisely their modulus values as their squared value will enter in the
collision terms (see Appendix B). In the aforementioned plot, we present the
simplest scenario where all 3 $\nu_{L}$ share the same temperature. We
consider benchmark parameter (BP) values $M_{Z^{\prime}}=10$ MeV with
$g_{X}=10^{-7}$ (Fig.2(a) & 2(c)) and $g_{X}=10^{-8}$ (Fig.2(b) &2(d)). We
will justify the importance of such a light $Z^{\prime}$ in this process at
the end of this section. For such a light $M_{Z^{\prime}}$ and
$g_{X}=10^{-7}$, we show the variation of $T_{\gamma}/T_{\nu}$ with
$T_{\gamma}$ and $\Delta N_{\rm eff}$ with $T_{\gamma}$ in Fig.2(a) and
Fig.2(c) respectively. From the above figure, it is very clear that around
$T_{\gamma}\gtrsim 10$MeV, $T_{\gamma}/T_{\nu}=1$ as both $\nu_{L}$ and
$Z^{\prime}$ was coupled to photon bath at that time. However, the ratio
starts to increase after $T_{\gamma}\sim 0.5$ MeV and then saturates at low
temperature ($T_{\gamma}\sim 10^{-2}$ MeV) as almost all the processes
mentioned earlier gradually become inefficient at low $T_{\gamma}$. In
Fig.2(c) we portray the corresponding variation in $\Delta N_{\rm eff}$ for
all the $U(1)_{X}$ charge combinations in Fig.2(a). As around $T_{\gamma}\sim
10^{-2}$ MeV the values of $\Delta N_{\rm eff}$ become saturated we can
surmise that it will remain unchanged till recombination epoch
($T_{\gamma}\sim 0.1$ eV) and say $\Delta N_{\rm eff}(T_{\gamma}\sim
10^{-2}{\rm MeV})\equiv\Delta N_{\rm eff}^{\rm CMB}$. In a similar way, for
$g_{X}=10^{-8}$ also we display the evolution of $T_{\gamma}/T_{\nu}$ and
$\Delta N_{\rm eff}$ in Fig.2(b) and Fig.2(d) respectively. In Fig.2(c) and
Fig.2(d) we also showcase the $2\sigma$ exclusion limit $\Delta N_{\rm
eff}^{\rm CMB}=0.28$ from Planck 2018 [4] and shown in black dotted line. Note
that this limit is only valid at CMB and it excludes any values of $\Delta
N_{\rm eff}$ at that time above the black dotted line. It is easy to infer
that in the presence of the light $Z^{\prime}$ the values differ from SM
prediction. Before spelling out the physical implications of the $U(1)_{X}$
scenario we tabulate our findings from Fig.2 for the ease of understanding in
Table2.

$U(1)_{\mathbb{X}}$ coupling |  $U(1)_{\mathbb{X}}$ charge |  Notation | $N_{\rm eff}^{\rm CMB}$
---|---|---|---
$|X_{1}|$ | $|X_{2}|$ | $|X_{3}|$ |  |
| 1 | 1 | 1 |  blue solid | 3.38
| 0 | 1 | 1 |  green dashed | 3.34
$g_{X}=10^{-7}$ | 1 | 0 | 0 |  brown dotted | 3.38
| 0 | 1 | 0 |  magenta dashed dot | 3.34
Fig.2(a,c) | 0 | 3 | 0 |  thick yellow | 3.34
| 3 | 0 | 0 |  red dashed | 3.47
| 1 | 1 | 1 |  blue solid | 3.21
| 0 | 1 | 1 |  green dashed | 3.34
$g_{X}=10^{-8}$ | 1 | 0 | 0 |  brown dotted | 3.21
| 0 | 1 | 0 |  magenta dashed dot | 3.34
Fig.2(b,d) | 0 | 3 | 0 |  thick yellow | 3.34
| 3 | 0 | 0 |  red dashed | 3.29
Table 2: Values of $N^{\rm CMB}_{\rm eff}$ for different $U(1)_{X}$ charges
inferred from Fig.2.

Both from Fig.2 and Table2 it is evident that in our proposed $U(1)_{X}$
scenario the value of temperature ratio or $N^{\rm CMB}_{\rm eff}$ differs
from the values predicted by SM only. One can interpret this feature from the
new processes involved in $\nu_{L}$ decoupling apart from the SM weak
interactions (see Fig.1). Thus the light $Z^{\prime}$ acts as the bridge
between photon and $\nu_{L}$ bath and tries to balance their energy densities
(through decays and scatterings) and hence reduces the temperature ratio from
the value predicted by the SM only. From both Fig.2(a) and Fig.2(b) we notice
that for a fixed value of $g_{X}$ and $X_{L_{1}}\neq 0$ the ratio
$T_{\gamma}/T_{\nu}$ (at very late time) and $\Delta N^{\rm CMB}_{\rm eff}$
grow with an increase in $|X_{1}|$. For $X_{1}\neq 0$ the following BSM
processes affect $\nu_{L}$ decoupling: (i) $Z^{\prime}$ decaying to both
$e^{+}e^{-}$ and $\nu_{i}\bar{\nu_{i}}$ and (ii) scattering process
$\nu_{i}\bar{\nu_{i}}\to e^{+}e^{-}$ mediated by $Z^{\prime}$. Thus with an
increase in $|X_{1}|$, the effective coupling ($X_{1}g_{X}$) governing these
BSM processes increases and hence boosts the BSM contribution (see eq.(43) and
eq.(46)) . As a result, picking higher values of $X_{1}$ leads to a higher
interaction rate between $\nu$ bath and photon bath leading to an enhancement
in $\Delta N^{\rm CMB}_{\rm eff}$ or, more precisely a diminution in
$T_{\gamma}/T_{\nu}$ .

On the contrary when $X_{1}=0$ the only BSM process relevant for $\nu_{L}$
decoupling is $Z^{\prime}\leftrightarrow\nu_{\mu,\tau}\bar{\nu_{\mu,\tau}}$ as
there is no tree level coupling of $Z^{\prime}$ with electrons. At
$T_{\nu}<M_{Z^{\prime}}$ eventually all $Z^{\prime}$ decay to $\nu_{L}$
transferring all their energy density to $\nu$ bath only. As all the
equilibrium number density of $Z^{\prime}$ finally gets diluted to $\nu_{L}$
bath (with $100\%$ branching ratio) it does not depend on the coupling
strength ($X_{2/3}g_{X}$). There is no change in $N_{\rm eff}$ with the change
in charge assignments ($X_{2},X_{3}$) for $X_{1}=0$. For the same reason
described above, we infer that for $X_{1}\neq 0$, $N_{\rm eff}$ increases with
an increase in $g_{X}$ whereas with $X_{1}=0$ it does not change at all with
change in $g_{X}$ (comparing Fig.2(c) and Fig.2(d) ). Due to the fact that for
$X_{1}=0$, $Z^{\prime}$ has only decay mode to $\nu_{\mu,\tau}$, $T_{\nu}$
starts to increase before $e^{\pm}$ decouples ($T_{\gamma}\sim 0.5$ MeV). This
causes a slight dip in the $T_{\gamma}/T_{\nu}$ evolution line at higher
temperature for $X_{1}=0$ in Fig.2(a) and Fig.2(b).

Before concluding the discussion about Fig.2, we point out the fact that we
compute $N_{\rm eff}$ with a basic assumption that all $\nu_{L}$ share same
temperature and hence all $\nu_{L}$ equilibrate with each other even if
$Z^{\prime}$ decays (transfers energy) to any one of them. For this reason we
observe similar behaviour of $T_{\gamma}/T_{\nu}$ for all $U(1)_{X}$ charge
combinations with same value of $X_{1}$ and the same phenomenology elaborated
in the earlier paragraph accounts for that. One should note that the
$U(1)_{X}$ charge combination used in Fig.2 is not an exhaustive one, yet one
can easily deduce the outcome of other combinations from the reasoning we made
above. For example all charge combinations with $|X_{1}|=1$ and $|X_{1}|=0$
will lead to $N_{\rm eff}^{\rm CMB}=3.38$ and $3.34$ respectively for
$g_{X}=10^{-7}$,$M_{Z^{\prime}}=10$ MeV. This is the most pivotal point of our
analysis and we will revisit this in the next section. At this juncture, it is
worth pointing out that the dependence of $N^{\rm CMB}_{\rm eff}$ on $g_{X}$
for $\mid X_{i}\mid=0$, indicating the absence of tree-level
$Z^{\prime}e^{+}e^{-}$ coupling. However, some effective coupling between
$Z^{\prime}$ and $e^{\pm}$ can be generated depending on the specific
$U(1)_{X}$ model [6]. In the subsequent section, we will extensively discuss
such a scenario where the induced $Z^{\prime}e^{+}e^{-}$ coupling plays a
crucial role in deciding $N^{\rm CMB}_{\rm eff}$. For the ease of our notation
from now on whenever we say $N_{\rm eff}$, we refer to $N_{\rm eff}^{\rm CMB}$
only.

Figure 3: Variation of $\Delta N^{\rm CMB}_{\rm eff}$ with $M_{Z^{\prime}}$
assuming all three $\nu_{L}$ share same temperature for different $U(1)_{X}$
charge combinations. We choose a fixed value of $g_{X}=10^{-7}$.

Having discussed the dependence of $N_{\rm eff}$ on the two most important
parameters of the BSM model: the $Z^{\prime}$ universal gauge coupling $g_{X}$
and $U(1)_{X}$ charge combinations, now we turn our attention to investigate
its dependence on the light $M_{Z^{\prime}}$. In Fig.3 we show variation of
$\Delta N_{\rm eff}^{\rm CMB}$ 777To amplify the change in $N_{\rm eff}$ with
$M_{Z^{\prime}}$ and portray more lucidly, for this particular plot we switch
to $\Delta N_{\rm eff}^{\rm CMB}\equiv N_{\rm eff}-3.046$. The variation in
$N_{\rm eff}$ will immediately follow from it.  with $M_{Z^{\prime}}$ for a
fixed coupling $g_{X}=10^{-7}$ with different $U(1)_{X}$ charge combinations.
Rather than showing all the $U(1)_{X}$ charge combinations used in Fig.2, we
just portray only three distinct combinations of them in Fig.3. We indicate
the combinations ($|X_{1,2,3}|=1$), ($|X_{1}|=0,|X_{2,3}|=1$), and
($|X_{1}|=3,|X_{2,3}|=0$) by blue, green and red lines respectively. From the
aforementioned figure, we observe that $\Delta N_{\rm eff}^{\rm CMB}$
decreases as $M_{Z^{\prime}}$ increases and eventually it becomes almost zero,
reproducing the SM value of $N_{\rm eff}$ when $M_{Z^{\prime}}\gtrsim 30$ MeV.
This feature can be interpreted from our previous discussion in the context of
Fig.1. The energy density of heavier $Z^{\prime}$ gets Boltzmann suppressed at
$\nu_{L}$ decoupling temperature ($T_{\gamma}\sim 2$ MeV) making the BSM
contribution less significant in deciding $\nu_{L}$ decoupling. On the other
hand $Z^{\prime}$ mediated scattering processes between $\nu$ and $e$, also
propagator suppressed for higher $M_{Z^{\prime}}$. At very high
$M_{Z^{\prime}}$ ($\gtrsim 30$ MeV), the BSM contribution hardly plays any
role in deciding $\nu_{L}$ decoupling resulting $\Delta N_{\rm eff}^{\rm
CMB}=0$. Following the same argument we infer that a lower value of
$M_{Z^{\prime}}$ will enhance the BSM contribution and hence $\Delta N_{\rm
eff}^{\rm CMB}$. For a fixed value of $M_{Z^{\prime}}$, the difference in
$\Delta N_{\rm eff}^{\rm CMB}$ for different $U(1)_{X}$ charge combinations
can be easily apprehended from the discussion in the context of Fig.2. So far
we explored the dependence of $\Delta N_{\rm eff}^{\rm CMB}$ on various model
parameters, and the dependence of $N_{\rm eff}$ also can be easily understood
from that. Keeping in mind the key findings from $U(1)_{X}$ $Z^{\prime}$
models in the context of $N_{\rm eff}$, in the following section we will
explore their contribution to alleviating the Hubble tension.

##  4 Numerical results

In the previous section, we pinned down the key aspects of light $Z^{\prime}$
from generic $U(1)_{X}$ extension and showed that the $U(1)_{X}$ charge
assignments play a key role in deciding $N_{\rm eff}$ as well as its
dependence on coupling $g_{X}$. In this section, we will take a closer look at
the model parameters and explore their cosmological implication through
exhaustive numerical scans. Though one can have numerous $U(1)_{X}$ charge
assignments as suggested in the Table 1., in the context of $N_{\rm eff}$ the
arbitrary charge assignments can be categorised into only two classes when we
assume all $\nu_{L}$ share the same temperature (through oscillation).
Needless to say, it is the coupling of $Z^{\prime}$ with electron (apart from
$\nu_{i}$) that affects the $N_{\rm eff}$ and not the couplings with
$\tau,\mu$ or quarks, for the reasons already discussed earlier. Thus the
light $Z^{\prime}$ models can be broadly classified into two pictures:
${X}_{1}\neq 0$ and ${X}_{1}=0$, more precisely, whether $Z^{\prime}$ has
coupling with electron or not (see Fig.4).

Figure 4: Classification of $Z^{\prime}$ models depending on $X_{1}=0$ and
$X_{1}\neq 0$

In our discussions on the light $Z^{\prime}$ phenomenology so far we have
mainly focused on its tree level couplings with fermions. Yet, it is important
to highlight that even in the absence of tree level $Z^{\prime}e^{+}e^{-}$
interaction ($X_{1}=0$), $Z^{\prime}$ can develop induced coupling888The
source of the induced coupling is model dependent. It can originate from
kinetic mixing, can be loop induced or can originate from flavor violating
couplings. In this model independent section we will remain agnostic about its
origin. with $e^{\pm}$ [54]. As a result, for $\mid X_{1}\mid=0$, while
computing $N_{\rm eff}$, we need to consider the following effective
$Z^{\prime}e^{+}e^{-}$ interaction Lagrangian:

$\mathcal{L}_{int}=(\epsilon e)\bar{e}\gamma^{\mu}eZ^{\prime}_{\mu},$ (8)

where $\epsilon$ is the induced effective coupling. If, for instance, the
induced coupling is generated from the $\gamma-Z^{\prime}$ kinetic mixing at
one loop level, its expression is given as (see Appendix A of ref. [54]).

$\epsilon=\frac{g_{X}}{2\pi^{2}}\sum_{\ell}X_{\ell}\int_{0}^{1}dx~{}x(1-x)\log\left(\frac{\Lambda^{2}}{\Delta_{\ell}}\right),$
(9)

where, $\Delta_{\ell}=m_{\ell}^{2}-x(1-x)q^{2}$, $\Lambda$ denotes an
arbitrary mass scale and the summation includes all $U(1)_{{X}}$ charged
fermions with mass $m_{\ell}$. This induced coupling will have a significant
impact in scenarios where $X_{1}=0$, as we’ll shortly discover. While we aim
to maintain a model-independent discussion, the effective coupling described
in eq.(8) relies on the characteristics of particular $U(1)_{{X}}$ models.
Therefore, to compare our findings with existing literature, we opt for the
benchmark value of the effective coupling, setting
$\epsilon=-\frac{g_{X}}{70}$, which is very commonly used for
$L_{\mu}-L_{\tau}$ models (with $X_{1}=0$)[55, 6]. Henceforth, in all our
numerical results, we will use this particular value of $\epsilon$.

(a) (b)

Figure 5: Variation of $\Delta N_{\rm eff}^{\rm CMB}$ with $g_{X}$ in (a) with
$\epsilon=0$ and in (b) with $\epsilon=-\frac{g_{X}}{70}$ assuming all 3
$\nu_{L}$ share same temperature for different $U(1)_{X}$ charge combinations.
We choose a fixed value of $M_{Z^{\prime}}=10$ MeV for both (a) and (b).

In Fig.5, we aim to investigate the role of $g_{X}$ on $\Delta N_{\rm
eff}^{\rm CMB}$ while maintaining a constant $M_{Z^{\prime}}=10$ MeV and as
stated earlier the variation of $N_{\rm eff}$ also follows from it. This
scrutiny involves two distinct scenarios for the $Z^{\prime}e^{+}e^{-}$
coupling: (a) solely with tree-level couplings ($\epsilon=0$) depicted in
Fig.5(a), and (b) considering induced coupling as well ($\epsilon\neq 0$)
shown in Fig.5(b). We display the variation of $\Delta N_{\rm eff}^{\rm CMB}$
with $g_{X}$ for three $U(1)_{\mathbb{X}}$ charge combinations
($|X_{1,2,3}|=1$), ($|X_{1}|=0,|X_{2,3}|=1$) and ($|X_{1}|=3,|X_{2,3}|=0$) by
blue, green and red lines respectively. In both Fig.5(a) and Fig.5(b), the
grey band corresponds to the value of $\Delta N^{\rm CMB}_{\rm eff}$ that is
excluded by the $2\sigma$ upper limit obtained from the Planck 2018
measurement [4]. Note that the dependence of $\Delta N_{\rm eff}^{\rm CMB}$
(and hence $N_{\rm eff}$) on $g_{X}$ for the cases with $X_{1}\neq 0$ remains
the same even after including mixing. This feature is easy to realize as the
$Z^{\prime}e^{+}e^{-}$ induced coupling ($\epsilon$) is suppressed by an order
of magnitude compared to the corresponding $Ze^{+}e^{-}$ tree level
interaction. Hence, in the presence of the tree-level $Z^{\prime}e^{+}e^{-}$
interaction $(X_{1}\neq 0)$, one can easily ignore the contribution of
$Z^{\prime}$ arising due to the induced coupling $\epsilon$. For the reasons
already discussed in sec.3, the BSM contribution increases as the values of
$g_{X}$ rise, leading to an increase in $\Delta N_{\rm eff}^{\rm CMB}$ (also
$N_{\rm eff}$) also, as shown in the aforementioned figure.

However, comparing Fig.5(a) and Fig.5(b) one notices that the dependence
$\Delta N_{\rm eff}^{\rm CMB}$ on $g_{X}$ for the case with $X_{1}=0$ changes
drastically when the induced coupling $(\epsilon\neq 0)$ is present. In
Fig.5(a), we discern that $\Delta N_{\rm eff}^{\rm CMB}$ remains unchanged
with variations in $g_{X}$ for $X_{1}=0$ in the absence of induced coupling
($\epsilon=0$). In this scenario, regardless of the associated coupling, the
$Z^{\prime}$ decay is limited to $\nu_{L}$ exclusively, transferring its
entire energy density, as reasoned in sec.3. This observation upholds our
earlier analysis discussed in the preceding section. On the other hand, when
$|X_{1}|=0$, in the presence of the induced coupling $(\epsilon\neq 0)$,
$Z^{\prime}$ can couple to $e^{\pm}$. Thus for a nonzero $\epsilon$ the decay
$Z^{\prime}\to e^{+}e^{-}$ and $\nu_{i}\bar{\nu_{i}}\to e^{+}e^{-}$ scattering
processes mediated by $Z^{\prime}$ continue during $\nu_{L}$ decoupling
temperature ($T_{\gamma}\sim 1$ MeV). Among these two processes, the
scattering process tries to balance $e$ and $\nu_{L}$ bath by increasing
$T_{\nu}$ or increasing $N_{\rm eff}$. Therefore, as $g_{X}$ increases, the
contribution of BSM scenarios also increases, resulting in an overall rise in
$\Delta N_{\rm eff}^{\rm CMB}$. This feature is reflected in Fig.5(b) (red
line) for $g_{X}\gtrsim 4\times 10^{-8}$. On the other extreme, for
$g_{X}\lesssim 4\times 10^{-8}$, the scattering process
($\propto\epsilon^{2}g_{X}^{2}$) is unable to compete with the tree level
decay $Z^{\prime}\to\nu_{L}\bar{\nu_{L}}$ ($\propto g_{X}^{2}$) process.
Hence, for lower values of $g_{X}(\lesssim 4\times 10^{-8})$, $Z^{\prime}$
promptly decays to $\nu$ bath transferring all its energy to $\nu$ sector and
the scattering processes become inefficient to dilute this extra energy
density to $e$ bath. So, when $g_{X}(\lesssim 4\times 10^{-8})$ we see a
distinct rise in $\Delta N_{\rm eff}^{\rm CMB}$ for $X_{1}=0$ in Fig.5(b).
When $g_{X}(\lesssim 7\times 10^{-8})$ is very small, the BSM contribution
that affects the evolution of the energy density of $\nu_{L}$ is primarily
dominated by $Z^{\prime}$ decay processes. Irrespective of the coupling, at
this level, the $Z^{\prime}$ particle transfers all its energy density to
$\nu_{L}$. It is essential to note that at this point, the value of $\Delta
N_{\rm eff}^{\rm CMB}$ for $X_{1}=0$ turns out to be identical for both
scenarios with and without induced coupling when comparing Fig.5(a) and
Fig.5(b). This outcome serves as validation for our previous argument. At this
point, it is worth pointing out another aspect of this tree-level vs induced
$Z^{\prime}$ coupling in generating non-zero contributions of BSM physics to
$\Delta N_{\rm eff}^{\rm CMB}$. When $g_{X}$ exceeds $10^{-8}$ and $X_{1}=0$,
the other two $U(1)_{X}$ charge combinations with non-zero $\mid X_{1}\mid$
result in the tree-level $Z^{\prime}e^{+}e^{-}$ interaction yielding greater
contributions to $\Delta N_{\rm eff}^{\rm CMB}$ compared to the aforementioned
induced $Z^{\prime}$ coupling. This distinction is particularly noticeable in
Fig. 5(b).

Now, we put forward the main thrust of this paper and split the analysis of
light $Z^{\prime}$ realized in different $U(1)_{{X}}$ models in the context of
$N_{\rm eff}$ into two categories as shown in Fig.4. The consequences of all
other $U(1)_{{X}}$ charge combinations can be easily anticipated from the
broad classifications. Following this line of thought we will now perform
numerical scans to explore the imprints in $N_{\rm eff}$ in the presence of
light $Z^{\prime}$ emerging from two broad classes of generic $U(1)_{X}$
models in the following two subsections.

###  4.1 $Z^{\prime}$ having tree level coupling with $e^{\pm}$

Figure 6: The parameter space for $N_{\rm eff}$ is depicted in the
$M_{Z^{\prime}}$ vs. $g_{X}$ plane for a generic $U(1)_{{X}}$ gauge extension
with the charge assignment $|X_{1,2,3}|=1$. The 2$\sigma$ upper bound from
Planck 2018 [4] is depicted by the grey dashed line, excluding the parameter
space to the left of that line shown by the grey region. The blue and red
dashed lines indicate the values of $N_{\rm eff}^{\rm CMB}=3.2$ and $3.5$,
respectively, where the $H_{0}$ tension can be relaxed as pointed out in
ref.[30].

Here we consider the charge assignment $|X_{1,2,3}|=1$ and show the contours
of constant $N_{\rm eff}$999As stated earlier by $N_{\rm eff}$ we refer
$N_{\rm eff}$ at the time of CMB. in $M_{Z^{\prime}}$ vs. $g_{X}$ plane in
Fig.6. In the same plot, we also showcase the $2\sigma$ bound from Planck 2018
data [4] shown by the grey dashed line. The parameter space to the left of the
grey dashed line, as shown by the grey region, is excluded by Planck 2028 data
at $2\sigma$ [4]. From the figure we observe that for a fixed $g_{X}$, $N_{\rm
eff}$ decreases with an increase in $M_{Z^{\prime}}$ as we explained in the
context of Fig.3. For very high $M_{Z^{\prime}}$ the contribution to $N_{\rm
eff}$ becomes negligible. We consider the lowest value of $g_{X}=10^{-9}$ as
below that, $Z^{\prime}$ fails to thermalize in the early universe. Following
the discussion made earlier in this section, it is evident that increasing the
value of $|X_{1}|$ will lead to a gradual shift of the contour for Planck 2018
upper limit (grey dashed line) towards right i.e. towards higher values of
$M_{Z^{\prime}}$ .

###  4.2 $Z^{\prime}$ having induced coupling with $e^{\pm}$

(a) (b)

Figure 7: Parameter space space for $N_{\rm eff}$ in $M_{Z^{\prime}}$ vs.
$g_{X}$ plane for a generic $U(1)_{{X}}$ gauge extension with the charge
assignment $|X_{1}|=0,|X_{2,3}|=1$. We show the results without kinetic mixing
($\epsilon=0$) in (a) and with kinetic mixing ($\epsilon\neq 0$) in (b). The
grey dashed line in each figure represents the 2$\sigma$ upper bound from
Planck 2018 [4], excluding the parameter space to its left, as depicted by the
grey region. The blue and red dashed lines in each figure denote $N_{\rm
eff}^{\rm CMB}$ values at $3.2$ and $3.5$, respectively, indicating the
intermediate regions where the $H_{0}$ tension can be relaxed [30].

Here we consider the following charge assignment $|X_{1}|=0,|X_{2,3}|=1$ and
show the contours of constant $N_{\rm eff}$ in $M_{Z^{\prime}}$ vs. $g_{X}$
plane in Fig.7. Similar to Fig.6 we portray the contour lines for $N_{\rm
eff}=3.2$, $N_{\rm eff}=3.5$ and $2\sigma$ upper limit from Planck 2018 [4]
depicted by blue, red and grey dashed lines. The grey region to the left of
the grey dashed line in each figure is excluded by the 2$\sigma$ upper bound
from Planck 2018 data [4].

We show the numerical results for the case $|X_{1}|=0$ without induced
coupling ($\epsilon=0$) in Fig.7(a). The dependence of $N_{\rm eff}$ on the
model parameter $g_{X}$ is pretty straight forward as we noticed in earlier
Fig.5(a). For $|X_{1}|=0$, in the absence of induced coupling, $N_{\rm eff}$
stays unchanged despite the variation in $g_{X}$, resulting in distinct
vertical lines of $N_{\rm eff}$ contours in $M_{Z^{\prime}}-g_{X}$ plane in
Fig.7(a). Also, from the same figure we note a decrease in $N_{\rm eff}$ with
an increase in $M_{Z^{\prime}}$ for a fixed $g_{X}$, as elaborated earlier.
Note that, in the absence of induced coupling ($\epsilon=0$)the contour for
Planck 2018 upper limit (grey dashed line) will not change with increasing
$X_{2,3}$ for such kind ($X_{1}=0$) of $U(1)_{X}$ models as argued before in
context of Fig.5(a).

However, the results are quite different after including the induced coupling
($\epsilon\neq 0$) of $Z^{\prime}$ with electrons as shown in Fig.7(b). Due to
the induced coupling of $Z^{\prime}$ with electrons, the $N_{\rm eff}$
contours replicate the feature of the $|X_{1}|=1$ case (as shown in Fig.6) for
higher values of $g_{X}~{}(\gtrsim 10^{-8})$ with a knee like pattern around
$g_{X}\sim 10^{-8}$. Such non-trivial dependence is the consequence of the
interplay between tree-level decay ( $Z^{\prime}\to\nu_{L}\bar{\nu_{L}}$) and
the induced $Z^{\prime}$ mediated scattering as explained earlier in detail.
We notice the bend in the contour lines in Fig.7(b) since the collision term
accounting for the BSM contribution in computing $N_{\rm eff}$ is decay
dominated in lower coupling region ($g_{X}\lesssim 4\times 10^{-8}$). It is
worth mentioning that for $\epsilon\neq 0$ the contour for Planck 2018 upper
limit (grey dashed line) will shift towards right with increasing $X_{2,3}$
for $g_{X}\gtrsim 4\times 10^{-8}$ and will remain unchanged for
$g_{X}\lesssim 4\times 10^{-8}$.

Thus from this section, we propound that the cosmological imprints of light
$Z^{\prime}$ models due to different charge assignments leading to different
$U(1)_{X}$ models can be put under the same roof following our prescription.
However, in the model independent analysis, we ignored the experimental
constraints which are inevitably relevant for our parameter space. The detail
analysis of those constraints on light $Z^{\prime}$ depends on specific
$U(1)_{{X}}$ choice and will be reported elsewhere [56]. For completeness, we
will show the numerical results for some specific $U(1)_{{X}}$ models along
with experimental constraints in the next section.

##  5 Specific $U(1)_{X}$ Symmetries:$U(1)_{B_{3}-3L_{j}}$

Let’s now look at the type of $U(1)_{X}$ symmetry where the symmetry is
flavour dependent in both quarks as well as the lepton sector namely the
$U(1)_{B_{3}-3L_{j}}$ gauge symmetries [22]. In this case, we will consider
the following three gauged $U(1)$ symmetries: $U(1)_{B_{3}-3L_{e}}$,
$U(1)_{B_{3}-3L_{\mu}}$ and $U(1)_{B_{3}-3L_{\tau}}$.

Fields | $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ | $U(1)_{B_{3}-3L_{1}}$ | $U(1)_{B_{3}-3L_{2}}$ | $U(1)_{B_{3}-3L_{3}}$
---|---|---|---|---
$Q_{i}$ | $(3,2,\frac{1}{3})$ | $(0,0,1)$ | $(0,0,1)$ | $(0,0,1)$
$u_{i}$ | $(3,1,\frac{4}{3})$ | $(0,0,1)$ | $(0,0,1)$ | $(0,0,1)$
$d_{i}$ | $(3,1,-\frac{2}{3})$ | $(0,0,1)$ | $(0,0,1)$ | $(0,0,1)$
$L_{i}$ | $(1,2,-1)$ | $(-3,0,0)$ | $(0,-3,0)$ | $(0,0,-3)$
$l_{i}$ | $(1,1,-2)$ | $(-3,0,0)$ | $(0,-3,0)$ | $(0,0,-3)$
$\nu_{R_{i}}$ | $(1,1,0)$ | $(x_{1},x_{2},x_{3})$ | $(x_{1},x_{2},x_{3})$ | $(x_{1},x_{2},x_{3})$
$\Phi$ | $(1,2,1)$ | $0$ | $0$ | $0$
$\sigma$ | $(1,1,0)$ | $y$ | $y$ | $y$
Table 3: Particle content and $U(1)_{B_{3}-3L_{j}}$ gauge charge assignments
of Standard Model and new particles. The charges of $\nu_{R_{i}}$ can be fixed
by anomaly cancellation condition while the charge of $\sigma$ depends on the
details of the nature and model for neutrino mass generation.

For all three cases, the charges of $\nu_{R_{i}}$ can be fixed by the anomaly
cancellation conditions. A convenient anomaly free charge assignment for
$\nu_{R_{i}}$ is $\nu_{R_{i}}=-3,\nu_{R_{j}}=0$; $j\neq i$ for the
$U(1)_{B_{i}-3L_{j}}$ symmetry e.g. for say $U(1)_{B_{3}-3L_{2}}$ gauge
symmetry the charges can be $\nu_{R}\sim(0-3,0)$ [22]. The charges of the
$\sigma$ field depend on the details of the model and we will not go into
details of the model building unless needed 101010see for example ref.[22].
And as we mentioned earlier in sec.3, we assume this $\nu_{R}$ and $\sigma$ to
be heavy enough that they are irrelevant for the analysis of $N_{\rm eff}$.

Before exploring the phenomenology of the specific $U(1)_{X}$, we would like
to outline an overview of the exclusion bounds on the mass of the light-gauge
boson ($M_{Z^{\prime}}$) and corresponding gauge coupling ($g_{X}$) within the
mass range $M_{Z^{\prime}}\sim\mathcal{O}(10~{}{\rm MeV})$, which is our point
of interest. Here, we will briefly review various types of low-energy
experimental observations that can be used to constrain the scenarios
involving any $U(1)_{X}$. In the following subsections, we will present the
exclusion bound identified by each low-energy experiment for a specific
$U(1)_{B_{3}-3L_{i}}$ scenario ($i=e,\mu,\tau$), along with our cosmological
findings.

Elastic electron-neutrino scattering (E$\nu$ES): The elastic scatterings of
neutrinos with electrons
($\nu_{\alpha}~{}e\to\nu_{\alpha}~{}e,~{}\alpha=e,\mu,\tau$) in laboratory
experiments serve as one of the probes for non-standard interaction of
neutrino and electron with the light-gauged boson ($Z^{\prime}$). In SM, the
$e-\nu_{e}$ scattering involves both charged-current (CC) and neutral-current
(NC) weak interactions, whereas $e-\nu_{\mu,\tau}$ scattering is solely
governed CC interaction [28]. The elastic $e-\nu_{e,\mu,\tau}$ scattering can
be altered in the presence of an additional NC interaction mediated by the
$U(1)_{X}$ gauge boson $Z^{\prime}$, which can be probed in the low-energy
scattering experiments. Note that the interference terms in the matrix
amplitude between the SM (W and Z-mediated ) and BSM ($Z^{\prime}$-mediated),
play a critical role in altering the scattering rate. The interference term in
differential cross section is given by [54],

$\displaystyle\left[\frac{d\sigma_{i-e}}{dE_{R}}\right]_{\rm Interference}$
$\displaystyle=$
$\displaystyle\frac{\sqrt{2}G_{F}m_{e}}{\pi}\left(\frac{(g_{X}X_{1})(g_{X}X_{i})}{(2E_{R}m_{e}+M_{Z^{\prime}}^{2})}\right)$
(10)
$\displaystyle\times\left[(g_{L}^{i}+g_{R}^{i})\left(1-\frac{m_{e}E_{R}}{2E_{\nu}^{2}}\right)-g_{R}^{i}\frac{E_{R}}{E_{\nu}}\left(2-\frac{E_{R}}{E_{\nu}}\right)\right],$

where $E_{R}$ and $E_{\nu}$ denote electron recoil energy and incoming
neutrino energy respectively. $g_{L,R}$ is defined in appendix A. For the
cases where $X_{1}=0$ the term $(g_{X}X_{1})$ will be replaced by $\epsilon e$
in eq.(10). Thus observing the electron recoil rate imposes constraints on the
$M_{Z^{\prime}}-g_{X}$ plane for a model-specific scenario. The coupling
strength of the light gauge boson with leptons varies over gauge extensions.
As a result, the constraints will vary from model to model.

Borexino [19] and dark matter experiments like XENON 1T [57] dedicated to
measuring the electron recoil rate, can be relevant in the context of
$e-\nu_{\alpha}$ elastic scattering [20]. The Borexino experiment is designed
to study the low-energy solar neutrinos ($\nu_{e}$)111111E$\nu$ES event rates
for atmospheric neutrinos are negligible compared to solar ones [20]. produced
via decay of 7Be by observing the electron recoil rate through the neutrino-
electron elastic scattering process [20]. These solar neutrinos undergo
flavour change as they travel from the sun to the detector. This flavor
changing phenomena $\nu_{e}\to\nu_{\alpha}(\alpha=e,\mu,\tau)$ can be
accounted using the transition probability $P_{e\alpha}$ [58] . Therefore the
number of events for solar neutrinos ($\nu_{e}$) interacting with the
electrons in the Borexino will be
$N_{\nu_{e}}\propto\sum_{\alpha=e,\mu,\tau}P_{e\alpha}~{}\sigma_{e-\nu_{\alpha}}$
[58, 28]. Note that in Solar neutrino flux $P_{ee}(\sim
50\%)>P_{e\mu}=P_{e\tau}$ and the dominant contribution in interference term
is due to $\nu_{e}$ as it has both CC and NC interaction with electron[28].
Hence for $X_{1}\neq 0$ one can approximate the differential recoil rate [28],

$\left[\frac{dR}{dE_{R}}\right]_{\rm BSM}\propto
P_{ee}\left[\frac{d\sigma_{ee}}{dE_{R}}\right]_{\rm Interference}$ (11)

Note that when $X_{1}=0$ i.e. in the absence of tree-level coupling of
$Z^{\prime}$ with electron one has to consider the elastic scattering (via
loop induced coupling with $Z^{\prime}$) of an electron with $\nu_{\mu/\tau}$
with specific transition probabilities.

Coherent elastic neutrino-nucleus scattering (CE$\nu$NS): The COHERENT
experiment investigates coherent elastic scattering between neutrinos and
nucleus in CsI material ($\nu~{}N\to\nu~{}N$) [59, 29, 55]. This mode of
interaction opens up new opportunities for studying neutrino interactions,
including the introduction of a new light gauge boson in this case. In SM, the
elastic neutrino-nucleus scattering ($\nu-N$) takes place via the NC
interactions between neutrinos ($\nu_{e,\mu,\tau}$) and nucleons or more
precisely with the first generation of quarks ($q=\\{u,d\\}$). The additional
NC interaction introduced by the light gauge boson $Z^{\prime}$ of $U(1)_{X}$
can contribute to the CE$\nu$NS process. The modification to SM CE$\nu$NS, due
to the $Z^{\prime}$ can be utilized to constrain the parameter space of
$M_{Z^{\prime}}$ vs $g_{X}$ plane. Similar to E$\nu$ES, the exclusion limit is
also dependent on the specific gauged scenarios. This is because the coupling
strength of quarks and neutrinos with the light-gauged boson is influenced by
the specific gauge choice. However, in contrast to E$\nu$ES experiment,
CE$\nu$NS is lepton flavor independent and one has to consider interaction
with all 3 $\nu$. For a detailed discussion see ref.[20].

Supernova 1987A (SN1987A): Non-standard neutrino interactions with a light
gauge boson can also affect the cooling of core-collapse supernovae (SN) which
is a powerful source of neutrinos ($\nu_{e}$) [60, 40]. The non-standard
neutrino interaction gives rise to the initial production of light-gauged
bosons in the core of a supernova (SN),
$\nu_{\alpha}~{}\overline{\nu_{\alpha}}\to Z^{\prime}$, resulting in energy
loss in the SN core [6]. After production the late decay of these $Z^{\prime}$
into neutrinos ($Z^{\prime}\to\nu_{i}~{}\overline{\nu_{i}}$) has the potential
to modify the neutrino flux emitted by the SN [40]. Additionally, the neutrino
self-interactions mediated by the light $Z^{\prime}$ may also affect energy
loss in the SN core [61]. Thus the observation of the SN 1987A neutrino burst
imposes constraints on these interactions [62]. However, there can be
additional processes that might place stronger constraints from SN 1987A on
light $Z^{\prime}$ models [63] and explicit derivation of that is beyond the
scope of this work.

Neutrino trident: The neutrino trident production mode, in which neutrinos
interact with a target nucleus to produce a pair of charged leptons without
changing the neutrino flavor ($\nu~{}N\to\nu~{}N~{}\mu^{+}~{}\mu^{-}$), is a
powerful tool for probing new physics [64]. The gauge extension introduces new
interactions between leptons and light gauge bosons, increasing the rate of
neutrino trident production, as predicted by the SM. Therefore any gauge
extended scenarios associated with muon ($X_{2}\neq 0$) face constraints from
existing neutrino trident production experimental results [22].

After having a generic discussion on the constraints relevant to our analysis,
in the following two subsections we present our numerical results for specific
$U(1)_{X}$ gauge choice. We portray the bound from $N_{\rm eff}$ on these
models along with the constraints discussed above. In the same plane, we also
indicate the parameter space that can relax the $H_{0}$ tension.

###  5.1 Gauged $U(1)_{B_{3}-3L_{e}}$ symmetry

In this subsection we show the values of $N_{\rm eff}$ for light $Z^{\prime}$
in $U(1)_{B_{3}-3L_{e}}$ gauge extension. We show our numerical results in
$M_{Z^{\prime}}$ vs. $g_{X}$ plane in Fig.8. In the same plane, we also
portray other relevant astrophysical and experimental constraints. The grey
dashed line in the same plot signify the $2\sigma$ upper limit on $N_{\rm
eff}$ from Planck 2018 data [4]. It excludes the parameter space to the left
of that line as shown by the grey shaded region. One can easily infer from the
figure that for a fixed $g_{X}$ with increasing $M_{Z^{\prime}}$ the value of
$N_{\rm eff}$ decreases which is explained in detail in the previous section.
Note that for $U(1)_{B_{3}-3L_{e}}$, $Z^{\prime}$ has tree level coupling with
$e^{\pm}$ ($X_{1}=-3$) and hence it belongs to the group of $U(1)_{X}$ models
with $|X_{1}|\neq 0$ following our discussion in previous section 4.1. For
this, we notice that the behavior of $N_{\rm eff}$ with $g_{X}$ is similar to
what we observed in the context of Fig.6. However, for $U(1)_{B_{3}-3L_{e}}$
model $|X_{1}|=3$ and hence for fixed $g_{X}$ and $M_{Z^{\prime}}$ the value
of $N_{\rm eff}$ will be higher compared to the case considered in Fig.6 for
the reason already elaborated before. As a result we notice the $N_{\rm eff}$
lines in Fig. 8 has moved rightwards compared to Fig.6.

Figure 8: Constraints from cosmology and low energy experiments in
$M_{Z^{\prime}}$ vs. $g_{X}$ plane for a generic $U(1)_{B_{3}-3L_{e}}$ gauge
extension. The $2\sigma$ upper bound from Planck 2018 [4] is shown by the grey
dashed line and it excludes the parameter space to the left side of that line
as shown by the grey shaded region. The region between the $N_{\rm eff}=3.2$
(blue dashed) and $N_{\rm eff}=3.5$ (red dashed) lines helps relax the $H_{0}$
tension, as mentioned in the ref. [30].

In the same plane in Fig.8 we showcase the relevant phenomenological
constraints that we mentioned earlier. Since $Z^{\prime}$ has tree level
coupling with electron in $U(1)_{B_{3}-3L_{e}}$ model, it attracts strong
constraint from E$\nu$ES [20]. The exclusion region from E$\nu$ES with XENON1T
experiment is shown by the region shaded by dark cyan diagonal lines. In a
similar fashion, Borexino also constrains the parameter space shown by the
magenta shaded region [20]. In $U(1)_{B_{3}-3L_{e}}$ model the light
$Z^{\prime}$ has no tree-level coupling with the first two generations of
quarks $(u,d),(c,s)$ which is important for CE$\nu$NS. However, $Z^{\prime}$
can have induced coupling with first-generation quarks via CKM mixing [22] and
can contribute to CE$\nu$NS. In SM the CKM matrix $V_{\rm CKM}$ is constructed
from charge current interaction betweeen quarks [65],

$\displaystyle J^{\mu}_{W}$ $\displaystyle=$
$\displaystyle(\bar{q_{u}})_{L}~{}\mathcal{U}_{L}~{}\gamma^{\mu}~{}\mathbb{I}_{3\times
3}~{}\mathcal{D}_{L}^{\dagger}~{}(q_{d})_{L}$ (12)

where, $q_{u}\equiv(u~{}c~{}t)^{\dagger}$ and
$q_{d}\equiv(d~{}s~{}b)^{\dagger}$. $\mathcal{U}_{L}$ and $\mathcal{D}_{L}$
are rotation matrices that diagonalize mass matrices for $q_{u}$ and $q_{d}$
respectively. The CKM matrix can be constructed as $V_{\rm
CKM}=\mathcal{U}_{L}\mathcal{D}_{L}^{\dagger}$ [65]. For our estimation we
take $\mathcal{U}_{L}=\mathbb{I}$ and hence $\mathcal{D}_{L}=V_{\rm
CKM}^{\dagger}$[22]. In our scenario as $Z^{\prime}$ couples to only third
generation of quarks (with $\mathcal{L}_{int}\supset
g_{X}J^{\mu}_{Z^{\prime}}Z^{\prime}_{\mu}$) it may induce flavor changing
neutral current (FCNC) and the Noether current ($J^{\mu}_{Z^{\prime}}$) is
given by,

$\displaystyle J^{\mu}_{Z^{\prime}}$ $\displaystyle=$
$\displaystyle(\bar{q_{d}})_{L}~{}\mathcal{D}_{L}~{}\gamma^{\mu}~{}\begin{pmatrix}0&0&0\\\
0&0&0\\\ 0&0&1\end{pmatrix}~{}\mathcal{D}_{L}^{\dagger}~{}(q_{d})_{L}.$ (13)

After inserting the values of the elements of $\mathcal{D}_{L}$ matrix [66],
in eq.(13) the $Z^{\prime}$ effective coupling with $d$ quark takes the
following form:

$\mathcal{L}\supset(0.006~{}g_{X})~{}Z^{\prime}_{\mu}\bar{d}_{L}\gamma^{\mu}d_{L}$
(14)

Using this we translate the bound from CE$\nu$NS experiment in our scenario
and obtain that it excludes the parameter space for $g_{X}\gtrsim 10^{-2}$.
This limit is weaker than E$\nu$ES and hence not shown in Fig.8. As
$Z^{\prime}$ has coupling with $b$ quark it generates $B_{s}-\overline{B_{s}}$
mixing and hence it attracts constraint indicated by orange shaded region
[22]. Finally, the bound from supernova cooling is shown by the region shaded
with grey vertical lines [6].

###  5.2 Gauged $U(1)_{B_{3}-3L_{\mu}}$ symmetry

(a) (b)

Figure 9: Constraints from cosmology and low energy experiments in
$M_{Z^{\prime}}$ vs. $g_{X}$ plane for a generic $U(1)_{B_{3}-3L_{\mu}}$ gauge
extension. We present the results without loop induced mixing ($\epsilon=0$)
in (a) and with loop induced mixing ($\epsilon\neq 0$) in (b). The grey dashed
line represents the $2\sigma$ upper bound from Planck 2018 [4], excluding the
parameter space on the left side of the line as shown by the grey shaded
region. The region between the blue dashed line ($N_{\rm eff}=3.2$) and the
red dashed line ($N_{\rm eff}=3.5$) indicates the parameter space where the
$H_{0}$ tension can be relaxed as mentioned in ref.[30].

In this subsection, we present our analysis with $N_{\rm eff}$ for a light
$Z^{\prime}$ in the gauged extension of $U(1)_{B_{3}-3L\mu}$. Analogous to the
previous subsection we show our numerical results in $M_{Z^{\prime}}$ vs.
$g_{X}$ plane in Fig.9 along with all relevant astrophysical and experimental
constraints. The grey dashed line in the same plot signify the $2\sigma$ upper
limit on $N_{\rm eff}$ from Planck 2018 data [4] and likewise before. The grey
shaded region depicts the exclusion of the parameter space to the left of the
grey dashed line.

However, unlike previous scenario, for $U(1)_{B_{3}-3L_{\mu}}$, $Z^{\prime}$
has no tree level coupling with $e^{\pm}$ ($X_{1}=0$) and hence it belongs to
other class of $U(1)_{X}$ models with $|X_{1}|=0$ as discussed in previous
sub-section 4.2. We present our numerical results with the two values of
effective coupling, $\epsilon=0$ and $\epsilon=-g_{X}/70$ in Fig.9(a) and
Fig.9(b) respectively. For $\epsilon=0$ we observe no change in $N_{\rm eff}$
with $g_{X}$ in Fig.9(a) as argued in the previous sections in the context of
Fig.7(a). On the other hand, for $\epsilon=-g_{X}/70$ we observe the non-
trivial dependence of $N_{\rm eff}$ with $g_{X}$ in Fig.9(a) due to the
interplay between scattering and decays as elaborated in detail in the context
of Fig.7(b).

In the same plane for both the plots Fig.9(a) and Fig.9(b) we portray the
constraints arising from E$\nu$ES experiment with XENON1T[20] shown by the
regions shaded by dark cyan diagonal lines. Similarly, we showcase the
exclusion region from Borexino depicted by magenta shaded region. Note that
these bounds are weaker for $U(1)_{B_{3}-3L_{\mu}}$ compared to Fig.8 as in
this case $Z^{\prime}$ has no tree level coupling with $e^{\pm}$ and the
recoil rate gets smaller here (see eq.(10)). Because, in this scenario
$Z^{\prime}$ couples directly with $\mu$, the Neutrino trident experiment
significantly constrains the parameter space shown by the dark blue dashed dot
line [22]. The parameter space above that line is excluded. In this model also
$Z^{\prime}$ can generate $B_{s}-\overline{B_{s}}$ mixing and attracts
constraint depicted by orange shaded region [22]. Similar to the previous
subsection the bound from CE$\nu$NS experiment is weaker than E$\nu$ES and
hence not shown in Fig.9. The bound from supernova cooling is shown by the
region shaded with grey vertical lines [6]. We do not show the results for
$U(1)_{B_{3}-3L\tau}$ explicitly, though the imprint in $N_{\rm eff}$ for this
model will be similar like the one in $U(1)_{B_{3}-3L\tau}$ as both these
models have $|X_{1}|=0$.

##  6 The Neutrino Portal Solution to the Hubble Tension

The change in $N_{\rm eff}$ due to the presence of light $Z^{\prime}$ gauge
bosons can also lead to an interesting consequence namely a simple and elegant
resolution to the cosmological Hubble constant problem. Recall that the Hubble
constant problem is the discrepancy between the obtained value of the Hubble
constant using early and late time probes [30, 33]. As mentioned in the
introduction the discrepancy in $H_{0}$ value can be ameliorated with
increasing $N_{\rm eff}$ [30]. Several different resolutions have been
proposed to address this discrepancy. One of the simplest possible resolution
is the presence of light degrees of freedom changing the $N_{\rm eff}$ [6, 9,
17, 18]. There are already a few model specific works successfully addressing
this problem for particular $U(1)_{X}$ extensions such as $U(1)_{\mu-\tau}$
[6] and $U(1)_{B-L}$ [46].

The general formalism developed in this work also can be used to show that the
$U(1)_{X}$ symmetries with light $Z^{\prime}$ gauge bosons, in general, can
address the Hubble constant problem reducing the tension. According to the
Planck 2015 TT data [37, 30], the $H_{0}$ tension issue can be relaxed upto
$1.8\sigma$ with the $N_{\rm eff}$ value between $3.2$ to $3.5$. However, the
Planck 2018 polarization measurements provide a more stringent bound on
$N_{\rm eff}$ [4] and it is difficult to reach $H_{0}>70~{}{\rm
Km~{}s}^{-1}\rm{Mpc}^{-1}$ [67] within the $\Lambda$CDM. Hence there still
exists a disagreement of $3.6\sigma$ in the $H_{0}$ value predicted from CMB
and local measurement. In our discussion, we only highlight the parameter
space where $N_{\rm eff}$ lies between $3.2$ to $3.5$ in
$M_{Z^{\prime}}-g_{X}$ plane, that may relax $H_{0}$ tension followed by the
analysis in ref. [37, 30]. In sec.4 we portray the parameter space for generic
$U(1)_{X}$ models in Fig.6 (for $|X_{1}|\neq 0$)and Fig.7 ($|X_{1}|=0$) where
the blue and red dashed lines indicate the values of $N_{\rm eff}=3.2$ and
$3.5$ respectively. Similarly, in sec.5 we also showcase the same lines with
constant $N_{\rm eff}=3.2$ (blue dashed line) and $3.5$ (red dashed line)
contours for the specific models: $U_{B_{3}-3L_{e}}$ (in Fig.8) and
$U_{B_{3}-3L_{\mu}}$ (in Fig. 9). The region between these two contour lines
may ameliorate $H_{0}$ tension. However claiming any possible resolution of
$H_{0}$ problem with quantitative effectiveness, itself requires a dedicated
analysis which is beyond the scope of this work [30].

##  7 Summary and conclusion

The Planck experiment has a very accurate measurement of the CMB which has
established stringent limits on $\Delta N_{\rm eff}$, representing the number
of effective relativistic degrees of freedom during the early epoch of the
universe [4]. This sensitivity makes $\Delta N_{\rm eff}$ a useful probe for
investigating various BSM scenarios that affect the neutrino decoupling at the
time of the CMB. In this work, we have studied the impact of light
$Z^{\prime}$ particle, arising from generic $U(1)_{X}$ model, on $N_{\rm
eff}$. In the presence of this light $Z^{\prime}$, $N_{\rm eff}$ receives two-
fold contributions: the decay of
$Z^{\prime}\to\nu_{L}{\bar{\nu}_{L}},~{}e^{+}e^{-}$ and scattering of light SM
leptons ($\nu_{L},e$) via this $Z^{\prime}$. At first, we considered a generic
$U(1)_{X}$ model with arbitrary charge assignments. Apart from the light
$Z^{\prime}$ gauge boson, the model also contains a BSM scalar $\sigma$ and
RHN $\nu_{R}$ which are necessary for the model construction. To understand
the sole effect of the light $Z^{\prime}$ in early universe temperature
evolution, we assume the other BSM particles are sufficiently heavy ($\gtrsim
100$ MeV) that they decouple at the time of neutrino decoupling. In this work,
we only consider the scenario where $Z^{\prime}$ was initially ($T\gtrsim
M_{Z^{\prime}}$) in a thermal bath and leave the other scenario with non-
thermal $Z^{\prime}$ for future work [56]. We adopt the formalism developed in
ref.[6] and solve the temperature equations to evaluate $N_{\rm eff}$ in
sec.3. We enumerate our findings below.

  * •

Firstly, we noticed that the $U(1)_{X}$ charges for $e,\nu_{e,\mu,\tau}$ play
a crucial role in $N_{\rm eff}$ for our proposed generic $U(1)_{X}$ model in
sec.2 as these are the only SM particles relevant for $\nu$ decoupling. In a
similar vein, we found that the mass of the BSM gauge boson $Z^{\prime}$
should be around $\lesssim 30$ MeV to affect neutrino decoupling.

  * •

After a careful analysis in sec.3 we observed that the value of $N_{\rm eff}$
depends differently on the coupling $g_{X}$ for two distinct scenarios those
are $|X_{1}|=0$ and $|X_{1}|\neq 0$. Based on this fact, in sec.4 we
categorise $U(1)_{X}$ models into two classes: (a) electrons having tree level
coupling with $Z^{\prime}$ ($|X_{1}|=0$) and (b) without tree level coupling
of an electron with $Z^{\prime}$ ($|X_{1}|\neq 0$). However, in the second
scenario electron may have some induced coupling with $Z^{\prime}$ even if
$|X_{1}|=0$ and we explored that scenario too. For both the two classes we
present the contours from the upper limit on $N_{\rm eff}$ from Planck 2018
[4] in $M_{Z^{\prime}}$ vs. $g_{X}$ plane as shown in Fig.6 and Fig.7. In the
same plane, we also highlight the parameter space favoured to relax $H_{0}$
tension shown in earlier analyses [37, 6].

  * •

As an example model, in sec.5 we considered a specific $U(1)_{X}$ model and
discuss the cosmological implications in context of $N_{\rm eff}$. We present
the numerical results for $U(1)_{B_{3}-3L_{e}}~{}(|X_{1}|=3)$ and
$U(1)_{B_{3}-3L_{\mu}}~{}(|X_{1}|=0)$ models in sub-sec.5.1 and sub-sec.5.2
respectively. Depending on the coupling of an electron with $Z^{\prime}$ these
two models also lead to distinguishable $N_{\rm eff}$ contours in
$M_{Z^{\prime}}$ vs. $g_{X}$ plane as portrayed in Fig.8 and Fig.9.

  * •

For both $U(1)_{B_{3}-3L_{e}}$ and $U(1)_{B_{3}-3L_{\mu}}$ models we also
showcase the relevant astrophysical and laboratory constraints in Fig.8 and
Fig.9 respectively. We displayed the constraints from E$\nu$ES with XENON1T
[20] and Borexino [58], neutrino trident [22], $B_{s}-\overline{B_{s}}$ mixing
[22] as well as from SN1987A [6]. We checked the constraints from CE$\nu$NS
[20] is weaker than the other constraints as our chosen $U(1)_{X}$ models do
not have $Z^{\prime}$ coupling with first generation quarks. For
$U(1)_{B_{3}-3L_{\mu}}$ model electron has no BSM gauge coupling and hence,
the bounds from E$\nu$ES for this model are weaker than the other one. In the
same plane, we also indicate the parameter space that can relax the $H_{0}$
tension [37, 6].

For both $U(1)_{B_{3}-3L_{e}}$ and $U(1)_{B_{3}-3L_{\mu}}$ models we have
shown that the bounds on $N_{\rm eff}$ from Planck 2018 data [4] can provide
more stringent bound on the parameter space than the laboratory searches. The
future generation experiments like CMB-S4 ( $\Delta N_{\rm eff}=0.06$ at
$2\sigma$) [68] can even probe more parameter space for such models. Our
analysis also shows that there is a certain parameter space still left to
relax $H_{0}$ tension [37, 6] allowed from all kinds of constraints. On the
other hand, as $U(1)_{X}$ models are well motivated BSM scenarios and widely
studied in several aspects, this analysis may enhance insights to explore
their connection with the $H_{0}$ problem too. Our generalised prescription
for $N_{\rm eff}$ analysis with generic $U(1)_{X}$ models is extremely helpful
to put stringent constraints on the parameter space from cosmology
complementary to the bounds obtained from gorund based experiments. In this
analysis we considered the thermal $Z^{\prime}$ scenario only and hence the
bounds from $N_{\rm eff}$ for extremely small couplings are not shown.
However, in future, we will explore the non-thermal $Z^{\prime}$ scenario also
to constrain the parameter space from $N_{\rm eff}$ even for the smaller
$Z^{\prime}$ mass and couplings [56].

###### Acknowledgements.

SJ thanks Sougata Ganguly for the helpful discussion and Miguel Escudero for
the email conversations. The work of SJ is funded by CSIR, Government of
India, under the NET JRF fellowship scheme with Award file No.
09/080(1172)/2020-EMR-I. PG would like to acknowledge the Indian Association
for the Cultivation of Science, Kolkata for the financial support.

##  Appendix A Neutrino decoupling in SM scenario

In this section, we recapitulate the formalism for $\nu_{L}$ decoupling in the
SM scenario which is already developed in ref.[5]. However, for better
understanding and smooth transition to the extra $U(1)_{X}$ scenario we just
note down the key points for the SM case here. As mentioned in sec.3 the only
particles relevant at the temperature scale of $\nu_{L}$ decoupling are
$\nu_{i}~{}(i=e,\mu,\tau)$,$e^{\pm}$ and $\gamma$. So in the SM scenario, the
temperature evaluation of only these $3$ particles is needed to evaluate
$T_{\nu}/T_{\gamma}$. Initially all $\nu_{i}~{}(i=e,\mu,\tau)$ were coupled to
$\gamma$ and $e^{\pm}$ bath at high temperature via weak interaction
processes. To track the exact temperature of photon bath ($T_{\gamma}$) and
neutrino bath ($T_{\nu}$) we start from the Liouville equation [50]:

$\displaystyle\frac{\partial f(p,t)}{\partial t}-Hf(p,t)\frac{\partial
f(p,t)}{\partial p}=\mathcal{C}[f],$ (15)

where,$f(p,t)$ is the distribution function in a homogeneous and isotropic
universe and $p,t$ signifies the momentum and time respectively.
$\mathcal{C}[f]$ is the collision term and it describes the total change in
$f(p,t)$ with time.

For SM $\nu_{L}$ decoupling the involved processes are $2\to 2$ scattering
between $\nu,e$ as shown in Table 4. So, we write the collision terms for a
particle $\chi$ undergoing $2\to 2$ process say, $\chi+i\to j+k$ for which
$\mathcal{C}[f]$ can be written as (assuming CP invariance),

$\displaystyle\mathcal{C}[f_{\chi}]$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{2E_{\chi}}\int
d\Pi_{i}d\Pi_{j}d\Pi_{k}(2\pi)^{4}\delta^{4}(p_{\chi}+p_{i}-p_{j}-p_{k})\times$
(16) $\displaystyle|\mathcal{M}|^{2}_{\chi+i\to
j+k}\left[f_{\chi}f_{i}\left[1\pm f_{j}\right]\left[1\pm
f_{k}\right]-f_{j}f_{k}\left[1\pm f_{\chi}\right]\left[1\pm
f_{i}\right]\right]\,$

where, $f_{\ell}$,$E_{\ell}$, and $p_{\ell}$ denote the distribution function,
energy, and 4-momentum respectively for particle $\ell$.
$|\mathcal{M}|^{2}_{\chi+i\to j+k}$ is the amplitude square for the above
mentioned process.

Multiplying eq.(15) with $g_{\chi}\,E\,d^{3}p_{\chi}/(2\pi)^{3}$ and
integrating over momenta will give the following energy density equation:

$\displaystyle\frac{d\rho_{\chi}}{dt}+3H(\rho_{\chi}+P_{\chi})=\int
g_{\chi}\,E\,\frac{d^{3}p_{\chi}}{(2\pi)^{3}}\,\mathcal{C}[f_{\chi}]=\frac{\delta\rho_{\chi\to
j}}{\delta t}$ (17)

$\frac{\delta\rho_{\chi\to j}}{\delta t}$ is the energy density transfer rate
from $\chi\to j({\rm or}~{}k)$. $P_{\chi}$ is the pressure density. So, we can
rewrite energy density equation eq.(17) for $\chi$,

$\displaystyle\frac{dT_{\chi}}{dt}$ $\displaystyle=$
$\displaystyle\left(-3H(\rho_{\chi}+P_{\chi})+\frac{\delta\rho_{\chi}}{\delta
t}\right)\left({\frac{\partial\rho_{\chi}}{\partial T_{\chi}}}\right)^{-1}$
(18)

Now for SM $\nu_{L}$ decoupling we consider the energy transfers between
$\nu_{i}~{}(i=e,\mu,\tau)$ to $e^{\pm},\gamma$ photon bath. For photon bath
containing $\gamma,e^{\pm}$ we consider a common temperature $T_{\gamma}$ and
add their energy density equations and rearrange accordingly,

$\displaystyle\frac{d\rho_{\gamma}}{dt}+\frac{d\rho_{e}}{dt}$ $\displaystyle=$
$\displaystyle-4H\rho_{\gamma}-3H(\rho_{e}+P_{e})+\sum_{i=e,\mu,\tau}\frac{\delta\rho_{e\to\nu_{i}}}{\delta
t}$ (19) $\displaystyle\frac{dT_{\gamma}}{dt}$ $\displaystyle=$
$\displaystyle\left(-4H\rho_{\gamma}-3H(\rho_{e}+P_{e})-\sum_{i=e,\mu,\tau}\frac{\delta\rho_{\nu_{i}\to
e}}{\delta t}\right)\left({\frac{\partial\rho_{\gamma}}{\partial
T_{\gamma}}+\frac{\partial\rho_{e}}{\partial T_{\gamma}}}\right)^{-1},$ (20)

where, in the last step we have used $\frac{\delta\rho_{\nu_{\chi}\to
j}}{\delta t}=-\frac{\delta\rho_{\nu_{j}\to\chi}}{\delta t}$. Similarly,
following eq.(17) we can write the respective temperature ($T_{\nu_{i}}$)
equations for each type of $\nu$,

$\displaystyle\frac{dT_{\nu_{i}}}{dt}$ $\displaystyle=$
$\displaystyle\left(-4H\rho_{\nu_{i}}+\frac{\delta\rho_{\nu_{i}}}{\delta
t}\right)\left(\frac{\partial\rho_{\nu_{i}}}{\partial
T_{\nu_{i}}}~{}\right)^{-1}$ (21)

Once we have the temperature equations the remaining step is to evaluate the
energy transfer rates i.e. the collision terms. The relevant processes in the
SM scenario are given in Table 4.

Process | $2^{-5}G_{F}^{2}\overline{|M|^{2}}$
---|---
$\nu_{i}+\bar{\nu_{i}}\longleftrightarrow e^{+}+e^{-}$ |  $4[g_{L_{i}}^{2}(p_{1}.p_{4})(p_{2}.p_{3})+g_{R_{i}}^{2}(p_{1}.p_{3})(p_{2}.p_{4})]$
$\nu_{i}+e^{-}\longleftrightarrow\nu_{i}+e^{-}$ | $4[g_{L_{i}}^{2}(p_{1}.p_{2})(p_{3}.p_{4})+g_{R_{i}}^{2}(p_{1}.p_{4})(p_{2}.p_{3})]$
$\nu_{i}+e^{+}\longleftrightarrow\nu_{i}+e^{+}$ | $4[g_{L_{i}}^{2}(p_{1}.p_{4})(p_{2}.p_{3})+g_{R_{i}}^{2}(p_{1}.p_{2})(p_{3}.p_{4})]$
$\nu_{i}+\nu_{j}\longleftrightarrow\nu_{i}+\nu_{j}$ | $(p_{1}.p_{2})(p_{3}.p_{4})$
$\nu_{i}+\bar{\nu_{j}}\longleftrightarrow\nu_{i}+\bar{\nu_{j}}$ | $(p_{1}.p_{2})(p_{3}.p_{4})$
Table 4: Squared amplitudes for interactions relevant for $\nu$ decoupling in
SM. The $i$ in subscripts stand for different generations of $\nu_{L}$ [52].

All the processes are either elastic scatterings between $\nu_{i}$ and $e$ or
annihilations and all are mediated by $W^{\pm},Z$. As we are focusing on the
MeV scale we can integrate out the mediators and write in terms of Fermi
constant $G_{F}$. In Table 4, $p_{1,2}$ and $p_{3,4}$ denote the momentum of
incoming and outgoing particles. Here, $g_{L_{e}}=\frac{1}{2}+s_{W}^{2}$,
$g_{R_{e}}=s_{W}^{2}$, $g_{L_{\mu,\,\tau}}=-\frac{1}{2}+s_{W}^{2}$,
$g_{R_{\mu,\,\tau}}=s_{W}^{2}$ and $s_{W}=\sin\theta_{W}$ where $\theta_{W}$
is Weinberg angle. The reason behind the difference of $1$ between $g_{L_{e}}$
and $g_{R_{\mu,\,\tau}}$ is the additional charged current process available
for $\nu_{e}$.

Collision terms for $\nu,e^{\pm}$ scattering

Process | $\pi^{5}G_{F}^{-2}\frac{\delta\rho}{\delta t}$
---|---
$\nu_{i}+\bar{\nu}_{i}\longleftrightarrow e^{+}+e^{-}$ |  $64(g_{L_{i}}^{2}+g_{R_{i}}^{2})\left[T_{\gamma}^{9}-T_{\nu_{i}}^{9}\right]$
$\nu_{i}+e^{\pm}\longleftrightarrow\nu_{i}+e^{\pm}$ | $112(g_{L_{i}}^{2}+g_{R_{i}}^{2})~{}T_{\gamma}^{4}T_{\nu_{i}}^{4}\left[T_{\gamma}-T_{\nu_{i}}\right]$
$\nu_{i}+\bar{\nu}_{i}\longleftrightarrow\nu_{j}+\bar{\nu}_{j}$ | $32\left[T_{\nu_{j}}^{9}-T_{\nu_{i}}^{9}\right]$
$\nu_{i}+\nu_{j}\longleftrightarrow\nu_{i}+\nu_{j}$ | $48T_{\nu_{j}}^{4}T_{\nu_{i}}^{4}\left[T_{\nu_{j}}-T_{\nu_{i}}\right]$
$\nu_{i}+\bar{\nu}_{j}\longleftrightarrow\nu_{i}+\bar{\nu}_{j}$ | $8T_{\nu_{j}}^{4}T_{\nu_{i}}^{4}\left[T_{\nu_{j}}-T_{\nu_{i}}\right]$
Table 5: Collision terms for $\nu,e^{\pm}$ scattering assuming MB distribution
function [47, 5]. The $i$ in subscripts stand for different generations of
$\nu_{L}$.

The first 2 processes are responsible for energy transfer from $\nu_{e}$ to
$e$ sector (photon bath) and the rest are responsible for energy transfer from
$\nu_{e}$ to $\nu_{\mu/\tau}$ sector. Finally, we enumerate below all the
energy transfer rates for the SM scenario following the formalism in ref.[47,
5]:

$\displaystyle\left(\frac{\delta\rho_{\nu_{e}\to e}}{\delta t}\right)_{\rm
SM}$ $\displaystyle=$
$\displaystyle\frac{G_{F}^{2}}{\pi^{5}}\left\\{(1+4s_{W}^{2}+8s_{W}^{4})F_{\rm
MB}(T_{\gamma},T_{\nu_{e}})\right\\}$ (22)
$\displaystyle\left(\frac{\delta\rho_{\nu_{\mu/\tau}\to e}}{\delta
t}\right)_{\rm SM}$ $\displaystyle=$
$\displaystyle\frac{G_{F}^{2}}{\pi^{5}}\left\\{(1-4s_{W}^{2}+8s_{W}^{4})F_{\rm
MB}(T_{\gamma},T_{\nu_{\mu}})\right\\}$ (23)
$\displaystyle\left(\frac{\delta\rho_{\nu_{e}\to\mu/\tau}}{\delta
t}\right)_{\rm SM}$ $\displaystyle=$
$\displaystyle\frac{G_{F}^{2}}{\pi^{5}}F_{\rm MB}(T_{\nu_{\mu}},T_{\nu_{e}}),$
(24)

where,

$\displaystyle F_{\rm
MB}(T_{1},T_{2})=32\,(T_{1}^{9}-T_{2}^{9})+56\,T_{1}^{4}\,T_{2}^{4}\,(T_{1}-T_{2})\,.$
(25)

Plugging these collision terms in eq.(20) and eq.(21) we can track the
temperature evolution.

##  Appendix B Neutrino decoupling in presence of a $U(1)_{X}$ gauge boson

The relevant lagrangian for the calculation of $N_{\rm eff}$ in $U(1)_{X}$
model has the following generic form,

$\mathcal{L}_{\rm int}\supset Z^{\prime}_{\alpha}J^{\alpha}_{\mathbb{X}},$
(26)

where, $J^{\alpha}_{\mathbb{X}}$ is the Noether current associated with
$U(1)_{\mathbb{X}}$ and given by,

$\displaystyle J^{\alpha}_{\mathbb{X}}$ $\displaystyle=$ $\displaystyle
g_{X}\left(X_{3}\bar{\tau}\gamma^{\alpha}\tau+X_{3}\bar{\nu}_{\tau}\gamma^{\alpha}P_{L}\nu_{\tau}+X_{2}\bar{\mu}\gamma^{\alpha}\mu+X_{2}\bar{\nu}_{\mu}\gamma^{\alpha}P_{L}\nu_{\mu}\right)$
(27)
$\displaystyle+g_{X}\left(X_{1}\bar{e}\gamma^{\alpha}e+X_{1}\bar{\nu}_{e}\gamma^{\alpha}P_{L}\nu_{e}\right)$

We have not written the corresponding interactions with quark sectors since
around $\sim 2$ MeV quarks are already confined [69]. The $\mu,\tau$ particles
also do not take part as they have suppressed energy density around MeV
temperature. In the presence of these $Z^{\prime}$ the following things will
affect neutrino decoupling (see Fig.1).

  1. 1.

Depending on charge assignments ($X_{1}\neq 0$) $Z^{\prime}$ can decay to
$e^{+}e^{-}$ and the inverse decay can also happen upto
$T_{\gamma}=M_{Z^{\prime}}$. So there will be energy transfer from
$Z^{\prime}$ to $e^{\pm}$ sector($\gamma$ bath) and vice versa. Similarly,
$Z^{\prime}$ can transfer energy to $\nu$ bath with the decays and inverse
decays to $\nu_{i}\bar{\nu_{i}}$.

  2. 2.

Through $Z^{\prime}$ portal $\nu_{i}+\overline{\nu_{i}}\longleftrightarrow
e^{+}e^{-}$,
$\nu_{i}+e^{\pm}\longleftrightarrow{\nu_{i}}+e^{\pm}~{},i=(e,\mu,\tau)$
processes can happen (depending on $X_{1,2,3}$) leading to an additional
energy transfer between $\nu_{i}$ sector and $e^{\pm}$.

We will discuss the above points in the following subsections. As discussed in
Sec.3, in this work we will only focus on the scenarios where $Z^{\prime}$ was
in thermal equilibrium with $e,\nu_{i}$ at $T>T^{\nu}_{\rm dec}$.

###  2.1 $Z^{\prime}$ decay to $e^{+}e^{-}$

We consider the process

$Z^{\prime}(k)\leftrightarrow e^{-}(p_{1})+e^{+}(p_{2}),$ (28)

where, $k,p_{1},p_{2}$ denote the four momenta of respective particles.
Following eq.(16) and eq.(17) we calculate the energy transfer rate from
$Z^{\prime}$ to $e$ bath,

$\displaystyle\int
g_{Z^{\prime}}E_{Z^{\prime}}\frac{d^{3}k}{(2\pi)^{3}2E_{Z^{\prime}}}\mathcal{C}[f]$
$\displaystyle=$ $\displaystyle-\int
g_{Z^{\prime}}\frac{d^{3}k}{(2\pi)^{3}2E_{Z^{\prime}}}g_{1}\frac{d^{3}p_{1}}{(2\pi)^{3}2E_{1}}g_{2}\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}\,E_{Z^{\prime}}\frac{1}{g_{Z^{\prime}}g_{1}g_{2}}|{\cal
M}_{e^{+}e^{-}\to Z^{\prime}}|^{2}$ (29)
$\displaystyle(f^{\rm(eq)}_{z^{\prime}}(k)-f^{\rm(eq)}_{e}(p_{1})f^{\rm(eq)}_{e}(p_{2}))(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-k)$

Here the degrees of freedom and energy of the corresponding particle are
denoted as $g_{\ell}\&E_{\ell},~{}(\ell\equiv Z^{\prime},1,2)$ respectively.
$f^{eq}_{i},~{}(i\equiv Z^{\prime},e)$ signify the equilibrium distribution
function of $i$ particle

The decay width of $Z^{\prime}$ to $e^{+}e^{-}$ is given by,

$\displaystyle\Gamma_{Z^{\prime}\to e^{+}e^{-}}$ $\displaystyle=$
$\displaystyle\frac{1}{2M_{Z^{\prime}}}\int\frac{d^{3}p_{1}}{(2\pi)^{3}2E_{1}}\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}\frac{1}{g_{Z^{\prime}}}|{\cal
M}_{Z^{\prime}\to e^{+}e^{-}}|^{2}(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-p_{k})$
(30) $\displaystyle=$ $\displaystyle
X_{1}^{2}\frac{g_{X}^{2}M_{Z^{\prime}}}{12\pi}(1+\frac{2m_{e}^{2}}{M^{2}_{Z^{\prime}}})\sqrt{1-\frac{4m_{e}^{2}}{M^{2}_{Z^{\prime}}}}$
(31)

where $m_{e}$ is the electron mass.

Using eq.(31) we can simplify above eq.(29) as

$\displaystyle\int
g_{z^{\prime}}E_{k}\frac{d^{3}k}{(2\pi)^{3}2E_{k}}\mathcal{C}[f]$
$\displaystyle=$ $\displaystyle\frac{3M_{Z^{\prime}}^{3}\Gamma_{Z^{\prime}\to
e^{+}e^{-}}}{2\pi^{2}}\left[T_{\gamma}~{}K_{2}\left(\frac{M_{Z^{\prime}}}{T_{\gamma}}\right)-T_{Z^{\prime}}K_{2}\left(\frac{M_{Z^{\prime}}}{T_{Z^{\prime}}}\right)\right]$
(32)

Therefore, from eq.(16) we write

$\frac{\delta\rho_{Z^{\prime}\to e}}{\delta
t}=\frac{3m_{Z^{\prime}}^{3}}{2\pi^{2}}\left[T_{\gamma}K_{2}\left(\frac{M_{Z^{\prime}}}{T_{\gamma}}\right)-T_{Z^{\prime}}K_{2}\left(\frac{M_{Z^{\prime}}}{T_{Z^{\prime}}}\right)\right]\Gamma_{Z^{\prime}\to
e^{+}e^{-}},~{}$ (33)

$T_{Z^{\prime}}$ denotes temperature of $Z^{\prime}$ bath.

###  2.2 $Z^{\prime}$ decay to $\nu_{i}\bar{\nu_{i}}$

Here we consider the processes

$Z^{\prime}(k)\leftrightarrow\nu_{i}(p_{1})+\bar{\nu_{i}}(p_{2}),~{}{\rm
where~{}}i=e,\mu,\tau.$ (34)

The corresponding decay widths of $Z^{\prime}$ to $i$ type $\nu$ are given as,

$\displaystyle\Gamma_{Z^{\prime}\to\nu_{i}\bar{\nu_{i}}}$ $\displaystyle=$
$\displaystyle X_{i}^{2}\frac{g_{X}^{2}M_{Z^{\prime}}}{24\pi}.$ (35)

Similar to the previous section we can write the energy density transfer rates
from $Z^{\prime}$ to each generation of $\nu$ as,

$\frac{\delta\rho_{Z^{\prime}\to\nu_{i}}}{\delta
t}=\frac{3m_{Z^{\prime}}^{3}}{2\pi^{2}}\left[T_{\nu_{i}}K_{2}\left(\frac{M_{Z^{\prime}}}{T_{\nu_{i}}}\right)-T_{Z^{\prime}}K_{2}\left(\frac{M_{Z^{\prime}}}{T_{Z^{\prime}}}\right)\right]\Gamma_{Z^{\prime}\to\nu_{i}\bar{\nu_{i}}},~{}$
(36)

###  2.3  Energy transfer from $\nu_{i}$ to $e^{+}e^{-}$ mediated by BSM
$Z^{\prime}$

In the presence of light $Z^{\prime}$ there will be additional energy transfer
between $e^{\pm}$ and $\nu_{i}$ apart from the SM weak processes (mediated by
SM $W\&Z$). The two relevant processes are:

$\displaystyle\nu(p_{1})+\bar{\nu}(p_{2})$ $\displaystyle\leftrightarrow$
$\displaystyle e^{-}(p_{3})+e^{+}(p_{4})$ (37)
$\displaystyle\nu(p_{1})+e^{\pm}(p_{2})$ $\displaystyle\leftrightarrow$
$\displaystyle\nu(p_{3})+e^{\pm}(p_{4})$ (38)

We denote the amplitudes for SM and BSM processes as $\mathcal{M}_{\rm SM}$
and $\mathcal{M}_{\rm BSM}$ respectively. Then the total amplitude
$\mathcal{M}_{\rm tot}$ can be written as

$\displaystyle|\mathcal{M}_{\rm tot}|^{2}=|\mathcal{M}_{\rm
SM}|^{2}+2Re(\mathcal{M}_{\rm SM}\mathcal{M}^{\dagger}_{\rm
BSM})+|\mathcal{M}_{\rm BSM}|^{2}$ (39)

Plugging this in eq.(16) we will get the total energy transfer rate
($\nu\leftrightarrow e$) Three sources contribute to the total energy transfer
rate; pure SM, SM-BSM interference, and pure BSM respectively. For simplicity,
we calculate the energy transfer rates for each of them separately. The
amplitudes for pure SM contributions can be found in sec.A. And the energy
transfer rates for pure SM contributions are given in eq.(22\- 24) [5].

The amplitude for BSM processes (with proper momentum assignment) will be
($s\ll M^{2}_{Z^{\prime}}$)

$\displaystyle\mathcal{M}_{BSM}=-\frac{(X_{i}g_{X})(X_{1}g_{X})}{M_{Z^{\prime}}^{2}}[\bar{\nu}_{i}\gamma^{\mu}P_{L}\nu_{i}]~{}[\bar{e}\gamma_{\mu}e]$
(40)

The interference amplitudes in eq.(39) are tabulated in Tab.-6.

Process | $2^{-3}/\sqrt{2}~{}G_{F}~{}(X_{i}g_{X})^{-1}~{}(X_{1}g_{X})^{-1}M_{Z^{\prime}}^{2}~{}\overline{|M|^{2}}$
---|---
$\nu_{i}+e^{-}\longleftrightarrow\nu_{i}+e^{-}$ | $~{}1/2~{}[g_{L}(p_{1}.p_{2})(p_{3}.p_{4})+g_{R}(p_{1}.p_{4})(p_{2}.p_{3})]$
$\nu_{i}+e^{+}\longleftrightarrow\nu_{i}+e^{+}$ | $~{}1/2~{}[g_{L}(p_{1}.p_{4})(p_{2}.p_{3})+g_{R}(p_{1}.p_{2})(p_{3}.p_{4})]$
$\nu_{i}+\bar{\nu_{i}}\longleftrightarrow e^{+}+e^{-}$ |  $[g_{L}(p_{1}.p_{4})(p_{2}.p_{3})+g_{R}(p_{1}.p_{3})(p_{2}.p_{4})]$
Table 6: Relevant SM-BSM interference terms for $\nu$ decoupling in SM. The
interactions are written with $i~{}(i=e,\mu,\tau)$ type $\nu$.

Following the formalism developed in ref.[47] we write the collision terms for
the interference terms:

$\displaystyle\int\frac{d^{3}p}{2\pi^{3}}p~{}\mathcal{C}[f]_{\nu_{i}+e\leftrightarrow\nu_{i}+e}$
$\displaystyle=$
$\displaystyle(X_{i}g_{X})(X_{1}g_{X})~{}\frac{7\sqrt{2}G_{F}}{\pi^{5}M_{Z^{\prime}}^{2}}(g_{L}+g_{R})~{}T_{\gamma}^{4}T_{\nu_{i}}^{4}\left[T_{\gamma}-T_{\nu_{i}}\right]$
(41)
$\displaystyle\int\frac{d^{3}p}{2\pi^{3}}p~{}\mathcal{C}[f]_{\nu_{i}\bar{\nu_{i}}\leftrightarrow
ee}$ $\displaystyle=$
$\displaystyle(X_{i}g_{X})(X_{1}g_{X})~{}\frac{4\sqrt{2}G_{F}}{\pi^{5}M_{Z^{\prime}}^{2}}(g_{L}+g_{R})~{}\left[T^{9}_{\gamma}-T^{9}_{\nu_{i}}\right]$
(42)

The first one in the above two collision terms accounts for elastic scattering
of $\nu_{i}$ with both $e^{\pm}$. Here we have multiplied with an additional
factor $2$ to count the effect of $\bar{\nu_{i}}$. The second collision term
in eq.(42) is due to annihilations. So, the total energy transfer rate from
$\nu$ to $e$ bath due to the interference with $U(1)_{X}$ gauge boson,

$\displaystyle\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta t}\right)_{\rm
int}$ $\displaystyle=$
$\displaystyle(X_{i}g_{X})(X_{1}g_{X})~{}\frac{\sqrt{2}G_{F}}{8\pi^{5}M_{Z^{\prime}}^{2}}\left\\{(g_{L}+g_{R})F_{\rm
MB}(T_{\gamma},T_{\nu_{i}})\right\\},$ (43)

In a similar fashion, we now move to calculate the collision term due to pure
BSM amplitudes (third term in eq.(39)). The pure BSM amplitudes in eq.(39) are
tabulated in Tab.-7.

Process | $2^{-4}~{}(X_{i}g_{X})^{-2}~{}(X_{1}g_{X})^{-2}~{}M_{Z^{\prime}}^{2}~{}\overline{|M|^{2}}$
---|---
$\nu_{i}+e^{-}\longleftrightarrow\nu_{i}+e^{-}$ | $~{}1/2~{}[(p_{1}.p_{2})(p_{3}.p_{4})+(p_{1}.p_{4})(p_{2}.p_{3})]$
$\nu_{i}+e^{+}\longleftrightarrow\nu_{i}+e^{+}$ | $~{}1/2~{}[(p_{1}.p_{4})(p_{2}.p_{3})+(p_{1}.p_{2})(p_{3}.p_{4})]$
$\nu_{i}+\bar{\nu_{i}}\longleftrightarrow e^{+}+e^{-}$ |  $[(p_{1}.p_{4})(p_{2}.p_{3})+(p_{1}.p_{3})(p_{2}.p_{4})]$
Table 7: Relevant BSM terms for $\nu$ decoupling in the SM. The interactions
are written with $i~{}(i=e,\mu,\tau)$ type $\nu$.

The collision terms will be:

$\displaystyle\int\frac{d^{3}p}{2\pi^{3}}p~{}\mathcal{C}[f]_{\nu_{i}+e\leftrightarrow\nu_{i}+e}$
$\displaystyle=$
$\displaystyle(X_{i}g_{X})^{2}(X_{1}g_{X})^{2}~{}\frac{14}{\pi^{5}M_{Z^{\prime}}^{4}}~{}T_{\gamma}^{4}T_{\nu_{i}}^{4}\left[T_{\gamma}-T_{\nu_{i}}\right]$
(44)
$\displaystyle\int\frac{d^{3}p}{2\pi^{3}}p~{}\mathcal{C}[f]_{\nu_{i}\bar{\nu_{i}}\leftrightarrow
ee}$ $\displaystyle=$
$\displaystyle(X_{i}g_{X})^{2}(X_{1}g_{X})^{2}~{}\frac{8}{\pi^{5}M_{Z^{\prime}}^{4}}~{}\left[T^{9}_{\gamma}-T^{9}_{\nu_{i}}\right]$
(45)

So, the total energy transfer rate from $\nu$ to $e$ bath due to the pure BSM
matrix element of $U(1)_{X}$ gauge boson,

$\displaystyle\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta t}\right)_{\rm
BSM}$ $\displaystyle=$
$\displaystyle(X_{i}g_{X})^{2}(X_{1}g_{X})^{2}~{}\frac{1}{4\pi^{5}M_{Z^{\prime}}^{4}}F_{\rm
MB}(T_{\gamma},T_{\nu_{i}}),$ (46)

So the total energy transfer rate from $\nu_{i}$ to $e$ bath will be

$\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta t}\right)_{\rm
tot}=\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta t}\right)_{\rm
SM}+\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta t}\right)_{\rm
int}+\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta t}\right)_{\rm BSM}$ (47)

$X_{1}=0$: Here we would like to highlight the scenarios where $X_{1}=0$. In
such case apparently, it seems that the BSM contributions in
eq.(29),eq.(45),eq.(46) vanish. However, in such cases even though
$Z^{\prime}$ doesn’t have any tree level coupling with $e^{\pm}$, it can
couple to $e^{\pm}$ through induced couplings which may generate from kinetic
mixing, fermion loops or CKM mixing. If the loop factor induced coupling is
denoted as $\epsilon$, then in such cases we have to just replace the term
$(X_{1}g_{X})$ in eq.(29),eq.(45),eq.(46) with $\epsilon~{}e$.

$(X_{1}g_{X})\longrightarrow(\epsilon e)~{}~{}~{}({\rm for}~{}~{}X_{1}=0)$
(48)

So, even if $X_{1}=0$, there can be energy transfer from $Z^{\prime}$ to $e$
bath as well as $\nu_{i}$ to $e$ bath by $Z^{\prime}$ mediation.

Note that eq.(47) differs from ref.[6] where the authors ignored the SM-BSM
interference term in the calculation of energy transfer rates. However, our
BSM contribution term $\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta
t}\right)_{\rm BSM}$ in eq.(46) agrees with ref.[6] for a specific choice of
$U(1)_{X}~{}(~{}U(1)_{L_{\mu}-L_{\tau}})$.

###  2.4 Temperature evolution

Now with all the collision terms, we are set to formulate the energy
evaluation equations. Following the prescription in appendix-A we can write
the following energy density evaluation equations:

$\displaystyle\frac{d\rho_{\nu_{i}}}{dt}$ $\displaystyle=$
$\displaystyle-\left(4H\rho_{\nu_{i}}-\left(\frac{\delta\rho_{\nu_{i}\to
e}}{\delta t}\right)_{\rm tot}-\sum_{j\neq
i}\left(\frac{\delta\rho_{\nu_{i}\to\nu_{j}}}{\delta t}\right)_{\rm
tot}+\frac{\delta\rho_{Z^{\prime}\to\nu_{i}}}{\delta t}\right)$ (49)
$\displaystyle\frac{d\rho_{Z^{\prime}}}{dt}$ $\displaystyle=$
$\displaystyle-\left(3H(\rho_{Z^{\prime}}+P_{Z^{\prime}})-\sum_{i}\frac{\delta\rho_{Z^{\prime}\to\nu_{i}}}{\delta
t}-\frac{\delta\rho_{Z^{\prime}\to e}}{\delta t}\right)$ (50)
$\displaystyle\frac{d\rho_{\gamma}}{dt}$ $\displaystyle=$
$\displaystyle-\left(4H\rho_{\gamma}+3H\left(\rho_{e}+p_{e}\right)+\sum_{i}\left(\frac{\delta\rho_{\nu_{i}\to
e}}{\delta t}\right)_{\rm tot}+\frac{\delta\rho_{Z^{\prime}\to e}}{\delta
t}\right),$ (51)

where the $\sum_{i}$ indicates the summation of energy density transfer rates
over all $3$ generatation of $\nu~{}(e,\mu,\tau)$. And
$\left(\frac{\delta\rho_{\nu_{i}\to\nu_{j}}}{\delta t}\right)_{\rm tot}$ bears
a similar kind of form as in the SM case. However, in the presence of the
light $Z^{\prime}$, it will get some additional contribution from $Z^{\prime}$
mediation. So the $G_{F}^{2}$ in eq.(24) will be simply replaced by
$G_{F}\to(G_{F}+(X_{i}g_{X})(X_{j}g_{X})/M_{Z^{\prime}}^{2})$. Similar to the
previous section we can write the temperature evaluation equations using the
partial derivatives and we get,

$\displaystyle\frac{dT_{\nu_{i}}}{dt}$ $\displaystyle=$
$\displaystyle-\left(4H\rho_{\nu_{i}}-\left(\frac{\delta\rho_{\nu_{i}\to
e}}{\delta t}\right)_{\rm tot}-\sum_{j\neq
i}\left(\frac{\delta\rho_{\nu_{i}\to\nu_{j}}}{\delta t}\right)_{\rm
tot}+\frac{\delta\rho_{Z^{\prime}\to\nu_{i}}}{\delta
t}\right)\left(\frac{\partial\rho_{\nu_{i}}}{\partial
T_{\nu_{i}}}\right)^{-1}$ (52) $\displaystyle\frac{dT_{Z^{\prime}}}{dt}$
$\displaystyle=$
$\displaystyle-\left(3H(\rho_{Z^{\prime}}+P_{Z^{\prime}})-\sum_{i}\frac{\delta\rho_{Z^{\prime}\to\nu_{i}}}{\delta
t}-\frac{\delta\rho_{Z^{\prime}\to e}}{\delta
t}\right)\left(\frac{\partial\rho_{Z^{\prime}}}{\partial
T_{Z^{\prime}}}\right)^{-1}$ (53) $\displaystyle\frac{dT_{\gamma}}{dt}$
$\displaystyle=$
$\displaystyle-\left(4H\rho_{\gamma}+3H\left(\rho_{e}+p_{e}\right)+\sum_{i}\left(\frac{\delta\rho_{\nu_{i}\to
e}}{\delta t}\right)_{\rm tot}+\frac{\delta\rho_{Z^{\prime}\to e}}{\delta
t}\right)\left(\frac{\partial\rho_{\gamma}}{\partial
T_{\gamma}}+\frac{\partial\rho_{e}}{\partial T_{\gamma}}\right)^{-1}$ (54)

Hence solving the aforementioned five Boltzmann equations on can track the
evolution of $\nu_{L}$ decoupling. By exploiting the neutrino oscillations
which are active around $\sim$ MeV temperature [51, 41, 52, 53], one can
further simplify the scenario by assuming all 3 $\nu_{i}$ equilibrate with
each other and acquire a common temperature i.e.
$T_{\nu_{e}}=T_{\nu_{\mu}}=T_{\nu_{\tau}}\equiv T_{\nu}$. The benefit of such
an assumption is that we no longer have to keep track of the energy transfer
rates within different $\nu$ sectors i.e.
$\left(\frac{\delta\rho_{\nu_{i}\to\nu_{j}}}{\delta t}\right)_{\rm tot}=0$ in
eq.(52). In such case the three $T_{\nu_{i}}$ equations (eq.(52)) simply
reduces to a single temperature ($T_{\nu}$) evaluation equation:

$\displaystyle\frac{dT_{\nu}}{dt}$ $\displaystyle=$
$\displaystyle-\left(4H\rho_{\nu}-\left(\frac{\delta\rho_{\nu\to e}}{\delta
t}\right)_{\rm tot}+\frac{\delta\rho_{Z^{\prime}\to\nu}}{\delta
t}\right)\left(\frac{\partial\rho_{\nu}}{\partial T_{\nu}}\right)^{-1},$ (55)

where, $\rho_{\nu}=\sum_{i=e,\mu,\tau}\rho_{\nu_{i}}$,
$\left(\frac{\delta\rho_{\nu\to e}}{\delta t}\right)_{\rm
tot}=\sum_{i=e,\mu,\tau}\left(\frac{\delta\rho_{\nu_{i}\to e}}{\delta
t}\right)_{\rm tot}$ and $\frac{\delta\rho_{Z^{\prime}\to\nu}}{\delta
t}=\sum_{i=e,\mu,\tau}\frac{\delta\rho_{Z^{\prime}\to\nu_{i}}}{\delta t}$

##  Appendix C Evaluation of $N_{\rm eff}$ with different $T_{\nu_{i}}$

Throughout the paper we evaluated $N_{\rm eff}$ assuming all three
$\nu_{i},~{}(i\equiv e,\mu,\tau)$ share the same temperature. However, the
values of $N_{\rm eff}$ change if we allow all $\nu_{i}$ to develop different
temperatures as prescribed in context of eq.(52-54). We show our results for
$g_{X}=10^{-8},~{}M_{Z^{\prime}}=10$ MeV for three different $U(1)_{X}$ charge
combinations in Fig.10 assuming only tree level couplings (i.e. $\epsilon=0$).

(a) (b) (c)

Figure 10: Evolution of $T_{\gamma}/T_{\nu}$ with photon bath temperature
$T_{\gamma}$ assuming all 3 $\nu_{i}$ with different temperatures for
different $U(1)_{X}$ charge combinations. We chose a benchmark parameter value
$M_{Z^{\prime}}=10$ MeV and $g_{X}=10^{-7}$ for all 3 plots. We consider 3
different charge combination (a)$X_{1}=-1,X_{2,3}=0$ (b)$X_{2}=-1,X_{1,3}=0$
and (a)$X_{3}=-1,X_{1,2}=0$. Here even if $X_{1}=0$ $N_{eff}$ changes slightly
depending on $X_{2,3}$ (in contrast to the case considered in the papaer.)

We consider $|X_{1}|=1,|X_{2,3}|=0$, $|X_{1}|=0,|X_{2,3}|=1$ and
$|X_{3}|=-1,|X_{1,2}|=0$ in Fig.10(a),Fig.10(b) and Fig.10(c) respectively.
The red, blue and green lines indicate the temperature ratio
$T_{\gamma}/T_{\nu_{i}}$ for $\nu_{e},\nu_{\mu}$ and $\nu_{\tau}$
respectively. Note that for $|X_{1}|=1,|X_{2,3}|=0$ only $\nu_{e}$ has BSM
interactions which lead to enhancing $T_{\nu_{e}}$, whereas
$T_{\nu_{\mu/\tau}}$ reproduces the same value as predicted in SM scenario
(where, $(T_{\gamma}/T_{\nu_{i}})_{\rm SM}\sim 1.4$ [5]) as portrayed in
Fig.10(a). It is interesting to note that this value of $N_{\rm eff}=3.09$ is
less ($\sim 8\%$) than the one predicted for the same charge configuration in
sec.3. This is due to the fact we overestimated the neutrino energy increment
in presence of $Z^{\prime}$ while assuming all 3 $\nu_{i}$ share the same
temperature.

For Fig.10(b) and Fig.10(c) also we notice the calculated values of $N_{\rm
eff}$ is slightly less than what was predicted in sec.3. For
$|X_{1}|=0,|X_{2,3}|=1$, $Z^{\prime}$ promptly decays to $\nu_{\mu}$ and
$\nu_{\tau}$ and not in electron. This causes a sudden dip (before
$T_{\gamma}\approx 0.5$ MeV ) in the respective temperature ratio curves shown
by the blue and green line in Fig.10(b). For the same reason we notice similar
dip in the temperature ratio curve (green line) for $\nu_{\tau}$ in Fig.10(c).
The value of $N_{\rm eff}$ is slightly higher for $|X_{1}|=0,|X_{2,3}|=1$ in
Fig.10(b) compared to the one obtained for Fig.10(c) as the former one
contains two type of $\nu_{i}$ having BSM interactions. Also it is worth
pointing that the value of $N_{\rm eff}$ obtained for $|X_{1}|=1,|X_{2,3}|=0$
is less compared to the other two as here $Z^{\prime}$ decays to both
$\nu_{e}$ and $e$ (i.e. photon bath) whereas in the other two cases
$Z^{\prime}$ decays only to $\nu_{i}$ bath.

## References

  * [1] S. Dodelson, Modern Cosmology.  Academic Press, Amsterdam, 2003.
  * [2] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico, “Relic neutrino decoupling including flavor oscillations,” Nucl. Phys. B 729 (2005) 221–234, arXiv:hep-ph/0506164.
  * [3] E. Grohs, G. M. Fuller, C. T. Kishimoto, M. W. Paris, and A. Vlasenko, “Neutrino energy transport in weak decoupling and big bang nucleosynthesis,” Phys. Rev. D 93 no. 8, (2016) 083522, arXiv:1512.02205 [astro-ph.CO].
  * [4] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [5] M. Escudero, “Neutrino decoupling beyond the Standard Model: CMB constraints on the Dark Matter mass with a fast and precise $N_{\rm eff}$ evaluation,” JCAP 02 (2019) 007, arXiv:1812.05605 [hep-ph].
  * [6] M. Escudero, D. Hooper, G. Krnjaic, and M. Pierre, “Cosmology with A Very Light Lμ $-$ Lτ Gauge Boson,” JHEP 03 (2019) 071, arXiv:1901.02010 [hep-ph].
  * [7] H. Esseili and G. D. Kribs, “Cosmological Implications of Gauged $U(1)_{B-L}$ on $\Delta N_{\rm eff}$ in the CMB and BBN,” arXiv:2308.07955 [hep-ph].
  * [8] S. Ganguly, S. Roy, and A. K. Saha, “Imprints of MeV scale hidden dark sector at Planck data,” Phys. Lett. B 834 (2022) 137463, arXiv:2201.00854 [hep-ph].
  * [9] K. N. Abazajian and J. Heeck, “Observing Dirac neutrinos in the cosmic microwave background,” Phys. Rev. D 100 (2019) 075027, arXiv:1908.03286 [hep-ph].
  * [10] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, “Early Dark Energy Can Resolve The Hubble Tension,” Phys. Rev. Lett. 122 no. 22, (2019) 221301, arXiv:1811.04083 [astro-ph.CO].
  * [11] D. K. Ghosh, S. Jeesun, and D. Nanda, “Long-lived inert Higgs boson in a fast expanding universe and its imprint on the cosmic microwave background,” Phys. Rev. D 106 no. 11, (2022) 115001, arXiv:2206.04940 [hep-ph].
  * [12] D. K. Ghosh, P. Ghosh, and S. Jeesun, “CMB signature of non-thermal Dark Matter produced from self-interacting dark sector,” JCAP 07 (2023) 012, arXiv:2301.13754 [hep-ph].
  * [13] V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, “Cosmological implications of ultralight axionlike fields,” Phys. Rev. D 98 no. 8, (2018) 083525, arXiv:1806.10608 [astro-ph.CO].
  * [14] T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg, and P. Walia, “Converting nonrelativistic dark matter to radiation,” Phys. Rev. D 98 no. 2, (2018) 023543, arXiv:1803.03644 [astro-ph.CO].
  * [15] S.-P. Li and X.-J. Xu, “Neff constraints on light mediators coupled to neutrinos: the dilution-resistant effect,” JHEP 10 (2023) 012, arXiv:2307.13967 [hep-ph].
  * [16] A. Biswas, D. K. Ghosh, and D. Nanda, “Concealing Dirac neutrinos from cosmic microwave background,” JCAP 10 (2022) 006, arXiv:2206.13710 [hep-ph].
  * [17] M. Berbig, S. Jana, and A. Trautner, “The Hubble tension and a renormalizable model of gauged neutrino self-interactions,” Phys. Rev. D 102 no. 11, (2020) 115008, arXiv:2004.13039 [hep-ph].
  * [18] M. Fabbrichesi, E. Gabrielli, and G. Lanfranchi, “The Dark Photon,” arXiv:2005.01515 [hep-ph].
  * [19] Borexino Collaboration, T. Shutt, “Borexino: A status report,” Nucl. Phys. B Proc. Suppl. 110 (2002) 323–325.
  * [20] A. Majumdar, D. K. Papoulias, and R. Srivastava, “Dark matter detectors as a novel probe for light new physics,” Phys. Rev. D 106 no. 1, (2022) 013001, arXiv:2112.03309 [hep-ph].
  * [21] E. Ma and R. Srivastava, “Dirac or inverse seesaw neutrino masses from gauged $B–L$ symmetry,” Mod. Phys. Lett. A 30 no. 26, (2015) 1530020, arXiv:1504.00111 [hep-ph].
  * [22] C. Bonilla, T. Modak, R. Srivastava, and J. W. F. Valle, “$U(1)_{B_{3}-3L_{\mu}}$ gauge symmetry as a simple description of $b\to s$ anomalies,” Phys. Rev. D 98 no. 9, (2018) 095002, arXiv:1705.00915 [hep-ph].
  * [23] B. Dutta, S. Ghosh, and J. Kumar, “Contributions to $\Delta N_{eff}$ from the dark photon of $U(1)_{T3R}$,” Phys. Rev. D 102 no. 1, (2020) 015013, arXiv:2002.01137 [hep-ph].
  * [24] CMS Collaboration, V. Khachatryan et al., “Search for narrow resonances in dilepton mass spectra in proton-proton collisions at $\sqrt{s}$ = 13 TeV and combination with 8 TeV data,” Phys. Lett. B 768 (2017) 57–80, arXiv:1609.05391 [hep-ex].
  * [25] ATLAS Collaboration, G. Aad et al., “Search for high-mass dilepton resonances using 139 fb-1 of $pp$ collision data collected at $\sqrt{s}=$13 TeV with the ATLAS detector,” Phys. Lett. B 796 (2019) 68–87, arXiv:1903.06248 [hep-ex].
  * [26] A. Das, S. Oda, N. Okada, and D.-s. Takahashi, “Classically conformal U(1)’ extended standard model, electroweak vacuum stability, and LHC Run-2 bounds,” Phys. Rev. D 93 no. 11, (2016) 115038, arXiv:1605.01157 [hep-ph].
  * [27] E. Accomando, L. Delle Rose, S. Moretti, E. Olaiya, and C. H. Shepherd-Themistocleous, “Extra Higgs boson and Z′ as portals to signatures of heavy neutrinos at the LHC,” JHEP 02 (2018) 109, arXiv:1708.03650 [hep-ph].
  * [28] P. Coloma, P. Coloma, M. C. Gonzalez-Garcia, M. C. Gonzalez-Garcia, M. Maltoni, M. Maltoni, J. a. P. Pinheiro, J. a. P. Pinheiro, S. Urrea, and S. Urrea, “Constraining new physics with Borexino Phase-II spectral data,” JHEP 07 (2022) 138, arXiv:2204.03011 [hep-ph]. [Erratum: JHEP 11, 138 (2022)].
  * [29] COHERENT Collaboration, D. Akimov et al., “Measurement of the Coherent Elastic Neutrino-Nucleus Scattering Cross Section on CsI by COHERENT,” Phys. Rev. Lett. 129 no. 8, (2022) 081801, arXiv:2110.07730 [hep-ex].
  * [30] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, “In the realm of the Hubble tension—a review of solutions,” Class. Quant. Grav. 38 no. 15, (2021) 153001, arXiv:2103.01183 [astro-ph.CO].
  * [31] A. G. Riess et al., “Milky Way Cepheid Standards for Measuring Cosmic Distances and Application to Gaia DR2: Implications for the Hubble Constant,” Astrophys. J. 861 no. 2, (2018) 126, arXiv:1804.10655 [astro-ph.CO].
  * [32] A. G. Riess et al., “A 2.4 $\%$ Determination of the Local Value of the Hubble Constant,” Astrophys. J. 826 no. 1, (2016) 56, arXiv:1604.01424 [astro-ph.CO].
  * [33] A. G. Riess et al., “A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km/s/Mpc Uncertainty from the Hubble Space Telescope and the SH0ES Team,” Astrophys. J. Lett. 934 no. 1, (2022) L7, arXiv:2112.04510 [astro-ph.CO].
  * [34] Planck Collaboration, P. A. R. Ade et al., “Planck 2015 results. XIII. Cosmological parameters,” Astron. Astrophys. 594 (2016) A13, arXiv:1502.01589 [astro-ph.CO].
  * [35] G. Efstathiou, “H0 Revisited,” Mon. Not. Roy. Astron. Soc. 440 no. 2, (2014) 1138–1152, arXiv:1311.3461 [astro-ph.CO].
  * [36] W. L. Freedman, “Cosmology at a Crossroads,” Nature Astron. 1 (2017) 0121, arXiv:1706.02739 [astro-ph.CO].
  * [37] J. L. Bernal, L. Verde, and A. G. Riess, “The trouble with $H_{0}$,” JCAP 10 (2016) 019, arXiv:1607.05617 [astro-ph.CO].
  * [38] T. Araki, K. Asai, K. Honda, R. Kasuya, J. Sato, T. Shimomura, and M. J. S. Yang, “Resolving the Hubble tension in a U(1)${}_{L_{\mu}-L_{\tau}}$ model with the Majoron,” PTEP 2021 no. 10, (2021) 103B05, arXiv:2103.07167 [hep-ph].
  * [39] H. Hildebrandt et al., “KiDS+VIKING-450: Cosmic shear tomography with optical and infrared data,” Astron. Astrophys. 633 (2020) A69, arXiv:1812.06076 [astro-ph.CO].
  * [40] K. Akita, S. H. Im, and M. Masud, “Probing non-standard neutrino interactions with a light boson from next galactic and diffuse supernova neutrinos,” JHEP 12 (2022) 050, arXiv:2206.06852 [hep-ph].
  * [41] S. Hannestad, “Oscillation effects on neutrino decoupling in the early universe,” Phys. Rev. D 65 (2002) 083006, arXiv:astro-ph/0111423.
  * [42] R. N. Mohapatra and G. Senjanovic, “Neutrino Mass and Spontaneous Parity Nonconservation,” Phys. Rev. Lett. 44 (1980) 912.
  * [43] J. Schechter and J. W. F. Valle, “Neutrino Masses in SU(2) x U(1) Theories,” Phys. Rev. D 22 (1980) 2227.
  * [44] A. Fradette, M. Pospelov, J. Pradler, and A. Ritz, “Cosmological beam dump: constraints on dark scalars mixed with the Higgs boson,” Phys. Rev. D 99 no. 7, (2019) 075004, arXiv:1812.07585 [hep-ph].
  * [45] C. D. Kreisch, F.-Y. Cyr-Racine, and O. Doré, “Neutrino puzzle: Anomalies, interactions, and cosmological tensions,” Phys. Rev. D 101 no. 12, (2020) 123505, arXiv:1902.00534 [astro-ph.CO].
  * [46] M. Escudero Abenza, “Precision early universe thermodynamics made simple: $N_{\rm eff}$ and neutrino decoupling in the Standard Model and beyond,” JCAP 05 (2020) 048, arXiv:2001.04466 [hep-ph].
  * [47] X. Luo, W. Rodejohann, and X.-J. Xu, “Dirac neutrinos and $N_{{\rm eff}}$,” JCAP 06 (2020) 058, arXiv:2005.01629 [hep-ph].
  * [48] A. D. Dolgov and K. Kainulainen, “Fermi-Dirac corrections to the relic abundances,” Nucl. Phys. B 402 (1993) 349–359, arXiv:hep-ph/9211231.
  * [49] M. Kawasaki, K. Kohri, and N. Sugiyama, “MeV scale reheating temperature and thermalization of neutrino background,” Phys. Rev. D 62 (2000) 023506, arXiv:astro-ph/0002127.
  * [50] E. W. Kolb and M. S. Turner, The Early Universe, vol. 69.  1990\.
  * [51] E. K. Akhmedov, “Neutrino physics,” in ICTP Summer School in Particle Physics, pp. 103–164.  6, 1999.  arXiv:hep-ph/0001264.
  * [52] A. D. Dolgov, “Neutrinos in cosmology,” Phys. Rept. 370 (2002) 333–535, arXiv:hep-ph/0202122.
  * [53] A. D. Dolgov, S. H. Hansen, S. Pastor, S. T. Petcov, G. G. Raffelt, and D. V. Semikoz, “Cosmological bounds on neutrino degeneracy improved by flavor oscillations,” Nucl. Phys. B 632 (2002) 363–382, arXiv:hep-ph/0201287.
  * [54] D. W. P. d. Amaral, D. G. Cerdeno, P. Foldenauer, and E. Reid, “Solar neutrino probes of the muon anomalous magnetic moment in the gauged $\mathrm{U}{(1)}_{L_{\mu}-{L}_{\tau}}$,” JHEP 12 (2020) 155, arXiv:2006.11225 [hep-ph].
  * [55] H. Banerjee, B. Dutta, and S. Roy, “Probing L$\mu$-L$\tau$ models with CE$\nu$NS: A new look at the combined COHERENT CsI and Ar data,” Phys. Rev. D 104 no. 1, (2021) 015015, arXiv:2103.10196 [hep-ph].
  * [56] D. K. Ghosh, P. Ghosh, S. Jeesun, and R. Srivastava, “Probing $U(1)_{\mathbb{X}}$ models from cosmology and ground based experiments,” arXiv:2401.xxxx [hep-ph].
  * [57] XENON Collaboration, E. Aprile et al., “The XENON1T Dark Matter Experiment,” Eur. Phys. J. C 77 no. 12, (2017) 881, arXiv:1708.07051 [astro-ph.IM].
  * [58] K. Chakraborty, A. Das, S. Goswami, and S. Roy, “Constraining general U(1) interactions from neutrino-electron scattering measurements at DUNE near detector,” JHEP 04 (2022) 008, arXiv:2111.08767 [hep-ph].
  * [59] M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, and Y. Y. Zhang, “Constraints on light vector mediators through coherent elastic neutrino nucleus scattering data from COHERENT,” JHEP 01 (2021) 116, arXiv:2008.05022 [hep-ph].
  * [60] E. W. Kolb and M. S. Turner, “Supernova SN 1987a and the Secret Interactions of Neutrinos,” Phys. Rev. D 36 (1987) 2895.
  * [61] P.-W. Chang, I. Esteban, J. F. Beacom, T. A. Thompson, and C. M. Hirata, “Toward Powerful Probes of Neutrino Self-Interactions in Supernovae,” Phys. Rev. Lett. 131 no. 7, (2023) 071002, arXiv:2206.12426 [hep-ph].
  * [62] H. T. Janka, “Neutrino Emission from Supernovae,” arXiv:1702.08713 [astro-ph.HE].
  * [63] D. Croon, G. Elor, R. K. Leane, and S. D. McDermott, “Supernova Muons: New Constraints on $Z$’ Bosons, Axions and ALPs,” JHEP 01 (2021) 107, arXiv:2006.13942 [hep-ph].
  * [64] W. Altmannshofer, S. Gori, M. Pospelov, and I. Yavin, “Neutrino Trident Production: A Powerful Probe of New Physics with Neutrino Beams,” Phys. Rev. Lett. 113 (2014) 091801, arXiv:1406.2332 [hep-ph].
  * [65] M. D. Schwartz, Quantum Field Theory and the Standard Model.  Cambridge University Press, 3, 2014.
  * [66] Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [67] A. S. de Jesus, N. Pinto-Neto, F. S. Queiroz, J. Silk, and D. R. da Silva, “The hubble rate trouble: an effective field theory of dark matter,” Eur. Phys. J. C 83 no. 3, (2023) 203, arXiv:2212.13272 [hep-ph].
  * [68] CMB-S4 Collaboration, K. Abazajian et al., “Snowmass 2021 CMB-S4 White Paper,” arXiv:2203.08024 [astro-ph.CO].
  * [69] A. R. Zhitnitsky, “Confinement Deconfinement Phase Transition in Hot and Dense QCD at Large N(c),” Nucl. Phys. A 813 (2008) 279–292, arXiv:0808.1447 [hep-ph].