aainstitutetext: Instituto de Física Teórica UAM/CSIC, Calle de Nicolás
Cabrera 13–15, Universidad Autónoma de Madrid, Cantoblanco, E-28049 Madrid,
Spainbbinstitutetext: Departament de Física Quàntica i Astrofísica and
Institut de Ciències del Cosmos, Universitat de Barcelona, Diagonal 647,
E-08028 Barcelona, Spainccinstitutetext: Institució Catalana de Recerca i
Estudis Avançats (ICREA), Pg. Lluis Companys 23, E-08010 Barcelona,
Spainddinstitutetext: C.N. Yang Institute for Theoretical Physics, Stony Brook
University, Stony Brook, NY 11794-3840, USA

# Neutrino Oscillation Constraints on $U(1)^{\prime}$ Models: from Non-
Standard Interactions to Long-Range Forces

Pilar Coloma<EMAIL_ADDRESS>b,c,d M. C. Gonzalez-Garcia
<EMAIL_ADDRESS>a Michele Maltoni
<EMAIL_ADDRESS>
###### Abstract

We quantify the effect of gauge bosons from a weakly coupled lepton flavor
dependent $U(1)^{\prime}$ interaction on the matter background in the
evolution of solar, atmospheric, reactor and long-baseline accelerator
neutrinos in the global analysis of oscillation data. The analysis is
performed for interaction lengths ranging from the Sun-Earth distance to
effective contact neutrino interactions. We survey $\sim 10000$ set of models
characterized by the six relevant fermion $U(1)^{\prime}$ charges and find
that in all cases, constraints on the coupling and mass of the $Z^{\prime}$
can be derived. We also find that about 5% of the $U(1)^{\prime}$ model
charges lead to a viable LMA-D solution but this is only possible in the
contact interaction limit. We explicitly quantify the constraints for a
variety of models including $U(1)_{B-3L_{e}}$, $U(1)_{B-3L_{\mu}}$,
$U(1)_{B-3L_{\tau}}$, $U(1)_{B-\frac{3}{2}(L_{\mu}+L_{\tau})}$,
$U(1)_{L_{e}-L_{\mu}}$, $U(1)_{L_{e}-L_{\tau}}$,
$U(1)_{L_{e}-\frac{1}{2}(L_{\mu}+L_{\tau})}$. We compare the constraints
imposed by our oscillation analysis with the strongest bounds from fifth force
searches, violation of equivalence principle as well as bounds from scattering
experiments and white dwarf cooling. Our results show that generically, the
oscillation analysis improves over the existing bounds from gravity tests for
$Z^{\prime}$ lighter than $\sim 10^{-8}\to 10^{-11}$ eV depending on the
specific couplings. In the contact interaction limit, we find that for most
models listed above there are values of $g^{\prime}$ and $M_{Z^{\prime}}$ for
which the oscillation analysis provides constraints beyond those imposed by
laboratory experiments. Finally we illustrate the range of $Z^{\prime}$ and
couplings leading to a viable LMA-D solution for two sets of models.

###### Keywords:

neutrino physics, solar and atmospheric neutrinos

††preprint: IFT-UAM/CSIC-133, YITP-SB-2020-30

##  1 Introduction

Experiments measuring the flavor composition of neutrinos produced in the Sun,
in the Earth’s atmosphere, in nuclear reactors and in particle accelerators
have established that lepton flavor is not conserved in neutrino propagation,
but it oscillates with a wavelength which depends on distance and energy. This
demonstrates beyond doubt that neutrinos are massive and that the mass states
are non-trivial admixtures of flavor states Pontecorvo:1967fh ; Gribov:1968kq
, see Ref. GonzalezGarcia:2007ib  for an overview.

When traveling through matter, the flavor evolution of the neutrino ensemble
is affected by the difference in the effective potential induced by elastic
forward scattering of neutrino with matter, the so-called Mikheev-Smirnov-
Wolfenstein (MSW) mechanism Wolfenstein:1977ue ; Mikheev:1986gs . Within the
context of the Standard Model (SM) of particle interactions, this effect is
fully determined and leads to a matter potential which, for neutral matter, is
proportional to the number density of electrons at the neutrino position,
$V=\sqrt{2}G_{F}N_{e}(r)$, and which only affects electron neutrinos. New
flavor dependent interactions can modify the matter potential and consequently
alter the pattern of flavor transitions, thus leaving imprints in the
oscillation data involving neutrinos which have traveled through large regions
of matter, as is the case for solar and atmospheric neutrinos.

Forward elastic scattering takes place in the limit of zero momentum transfer,
so as long as the range of the interaction is shorter than the scale over
which the matter density extends, the effective matter potential can be
obtained in the contact interaction approximation between the neutrinos and
the matter particles. The paradigmatic example is provided by neutral current
non-standard interactions (NSI) Wolfenstein:1977ue ; Valle:1987gv ;
Guzzo:1991hi  between neutrinos and matter (for recent reviews, see
Ohlsson:2012kf ; Miranda:2015dra ; Farzan:2017xzy ; Dev:2019anc ), which can
be parametrized as

$\mathcal{L}_{\text{NSI}}=-2\sqrt{2}G_{F}\sum_{f,P,\alpha,\beta}\varepsilon_{\alpha\beta}^{f,P}(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\beta})(\bar{f}\gamma_{\mu}Pf)\,,$
(1)

where $G_{F}$ is the Fermi constant, $\alpha,\beta$ are flavor indices,
$P\equiv P_{L},P_{R}$ and $f$ is a SM charged fermion. These operators are
expected to arise generically from the exchange of some mediator state heavy
enough for the contact interaction approximation to hold. In this notation,
$\varepsilon_{\alpha\beta}^{f,P}$ parametrizes the strength of the new
interaction with respect to the Fermi constant,
$\varepsilon_{\alpha\beta}^{f,P}\sim\mathcal{O}(G_{X}/G_{F})$. Generically
they modify the matter potential in neutrino propagation, but – being local
interactions – the resulting potential is still proportional to the number
density of particles in the medium at the neutrino position. Since such
modifications arise from a coherent effect, oscillation bounds apply even to
NSI induced by ultra light mediators, as long as their interaction length is
shorter than the neutrino oscillation length. For the experiments considered
here, such condition is fulfilled as long as $M_{Z^{\prime}}\gtrsim 10^{-12}$
eV GonzalezGarcia:2006vp .

Conversely, if the mediator is too light then the contact interaction
approximation is no longer valid, and the flavor dependent forces between
neutrino and matter particles become long-range. In this case neutrino
propagation can still be described in terms of a matter potential, which
however is no longer simply determined by the number density of particles in
the medium at the neutrino position, but it depends instead on the average of
the matter density within a radius $\sim 1/M_{Z^{\prime}}$ around it
Grifols:2003gy ; Joshipura:2003jh ; GonzalezGarcia:2006vp ; Davoudiasl:2011sz
; Wise:2018rnb .

At present, the global analysis of data from oscillation experiments provides
some of the strongest constraints on the size of the NSI affecting neutrino
propagation GonzalezGarcia:2011my ; Gonzalez-Garcia:2013usa ; Esteban:2018ppq
. Analysis of early oscillation data was also used to impose constraints on
flavor dependent long-range forces Grifols:2003gy ; Joshipura:2003jh ;
GonzalezGarcia:2006vp .

Straightforward constructions leading to Eq. (1) have an extended gauge sector
with an additional $U(1)^{\prime}$ symmetry with charge involving some of the
lepton flavors and an heavy enough gauge boson. Conversely if the gauge boson
is light enough a long range force will be generated. Thus the analysis of
neutrino oscillation data can shed light on the valid range of $Z^{\prime}$
mass and coupling in both regimes. Following this approach, the marginalized
bounds on the NSI coefficients derived from the global analysis of
oscillations in presence of NSI performed in Ref. Esteban:2018ppq  were
adapted to place constraints on the coupling and mass of the new gauge boson
both in the NSI limit Heeck:2018nzc  and in the long-range regime Wise:2018rnb
for several $U(1)^{\prime}$ flavor symmetries. However, strictly speaking, the
bounds derived in Ref. Esteban:2018ppq  cannot be directly used to constraint
the $U(1)^{\prime}$ scenarios because in the latter case only flavor diagonal
interactions (and only some of them depending on the $U(1)^{\prime}$ charge)
are generated, while the bounds derived in Ref. Esteban:2018ppq  were obtained
in the most general parameter space with all relevant four-fermion
interactions (flavor conserving and flavor changing) being simultaneously non-
vanishing. In order to derive statistically consistent bounds on each
$U(1)^{\prime}$ scenario a dedicated analysis has to be performed in its
reduced parameter space.

With this motivation, in this work we perform such dedicated global analysis
of oscillation data in the framework of lepton flavor dependent
$U(1)^{\prime}$ interactions which affect the neutrino evolution in matter,
with interaction lengths ranging from the Sun-Earth distance to effective
contact neutrino interactions. In Sec. 2 we describe the models which will be
studied and derive the matter potential generated both in the contact
interaction limit (in Sec. 2.1) and in the case of finite interaction range
(in Sec. 2.2) as a function of the $U(1)^{\prime}$ charges. The results of the
global analysis are presented in Sec. 3. In particular the bounds imposed by
the analysis and how they compare with those from other experiments are
presented in Sec. 3.1. An additional consideration that we take into account
is that in the presence of NSI a degeneracy exists in oscillation data,
leading to the so called LMA-Dark (LMA-D) Miranda:2004nb  solution first
observed in solar neutrinos, where for suitable NSI the data can be explained
by a mixing angle $\theta_{12}$ in the second octant. For this new solution to
appear the new interactions must be such that the matter potential difference
for electron neutrinos reverses its sign with respect to that in the SM. It is
not trivial to generate such large effects without conflicting with bounds
from other experiments, though models with light mediators (i.e., below the
electroweak scale) have been proposed as viable candidates Farzan:2017xzy ;
Miranda:2015dra ; Farzan:2015doa ; Farzan:2015hkd ; Babu:2017olk ;
Denton:2018xmq . Section 3.2 contains our findings on viable models for LMA-D.
We summarize our conclusions in Sec. 4. We present some details of the
translation of the bounds from some experiments to the models studied in an
appendix.

##  2 Formalism

We are going to focus on $U(1)^{\prime}$ models which can be tested in
neutrino oscillation experiments via its contribution to the matter potential.
As a start this requires that the new gauge boson couples to the fermions of
the first generation.

An important issue when enlarging the Standard Model with a new
$U(1)^{\prime}$ gauge group is the possibility of mixing between the three
neutral gauge bosons of the model which can, in general, be induced in either
kinetic or mass terms. While kinetic mixing is fairly generic as it can be
generated at the loop level with the SM particle contents, matter mixing is
model dependent as it requires an extended scalar sector with a vacuum
expectation value charged both under the SM and the $U(1)^{\prime}$. In what
respects the effect of the new $U(1)^{\prime}$ in oscillation experiments an
important observation is that if the new interaction does not couple directly
to fermions of the first generation, no matter effects can be generated by
kinetic mixing Heeck:2018nzc . Thus neglecting mixing effects yields the most
model independent and conservative bounds from neutrino oscillation results.
So in what follows we are going to work under the assumption that the
$Z^{\prime}$ mixing with SM gauge bosons can be safely neglected.

In addition we notice that only vector interactions contribute to the matter
potential in neutrino propagation. Altogether the part of the $U(1)^{\prime}$
Lagrangian relevant for propagation in ordinary matter has the most general
form

$\mathcal{L}_{Z^{\prime}}^{\text{matter}}=-g^{\prime}\big{(}a_{u}\,\bar{u}\gamma^{\alpha}u+a_{d}\,\bar{d}\gamma^{\alpha}d+a_{e}\,\bar{e}\gamma^{\alpha}e\\\
+b_{e}\,\bar{\nu}_{e}\gamma^{\alpha}P_{L}\nu_{e}+b_{\mu}\,\bar{\nu}_{\mu}\gamma^{\alpha}P_{L}\nu_{\mu}+b_{\tau}\,\bar{\nu}_{\tau}\gamma^{\alpha}P_{L}\nu_{\tau}\big{)}Z^{\prime}_{\alpha}$
(2)

with arbitrary charges $a_{u,d,e}$ and $b_{e,\mu,\tau}$.

These charges can be accommodated in generalized anomaly free UV-complete
models including only the SM particles plus right-handed neutrinos
Allanach:2018vjg . If, in addition, one requires all couplings to be vector-
like and the quark couplings to be generation independent, the condition of
anomaly cancellation for models with only SM plus right handed neutrinos
imposes constrains over the six charges above and one ends with a subclass of
models characterized by three independent charges which can be chosen to be,
for example, $B-L$, $(L_{\mu}-L_{\tau})$, and $(L_{\mu}-L_{e})$ Heeck:2018nzc
:

$c_{\textsc{bl}}(B-L)+c_{\mu\tau}(L_{\mu}-L_{\tau})+c_{\mu
e}(L_{\mu}-L_{e})\,,$ (3)

so $a_{u}=a_{d}=c_{\textsc{bl}}/3$, $a_{e}=b_{e}=-(c_{\textsc{bl}}+c_{\mu
e})$, $b_{\mu}=-c_{\textsc{bl}}+c_{\mu e}+c_{\mu\tau}$, and
$b_{\tau}=-(c_{\textsc{bl}}+c_{\mu\tau})$. In particular models with charges
as in Eq. (3) are not constrained by rare $Z$ decays or flavor changing
neutral current (FCNC) meson decays Dror:2017ehi ; Dror:2017nsg ; Dror:2018wfl
. For convenience, in table 1 we list the charges corresponding to some of the
models discussed in the literature.

Model | $a_{u}$ | $a_{d}$ | $a_{e}$ | $b_{e}$ | $b_{\mu}$ | $b_{\tau}$
---|---|---|---|---|---|---
$B-3L_{e}$ | $\frac{1}{3}$ | $\frac{1}{3}$ | $-3$ | $-3$ | $\hphantom{-}0$ | $\hphantom{-}0$
$B-3L_{\mu}$ | $\frac{1}{3}$ | $\frac{1}{3}$ | $\hphantom{-}0$ | $\hphantom{-}0$ | $-3$ | $\hphantom{-}0$
$B-3L_{\tau}$ | $\frac{1}{3}$ | $\frac{1}{3}$ | $\hphantom{-}0$ | $\hphantom{-}0$ | $\hphantom{-}0$ | $-3$
$B-\frac{3}{2}(L_{\mu}+L_{\tau})$ | $\frac{1}{3}$ | $\frac{1}{3}$ | $\hphantom{-}0$ | $\hphantom{-}0$ | $-\frac{3}{2}$ | $-\frac{3}{2}$
$L_{e}-L_{\mu}$ | $0$ | $0$ | $\hphantom{-}1$ | $\hphantom{-}1$ | $-1$ | $\hphantom{-}0$
$L_{e}-L_{\tau}$ | $0$ | $0$ | $\hphantom{-}1$ | $\hphantom{-}1$ | $\hphantom{-}0$ | $-1$
$L_{e}-\frac{1}{2}(L_{\mu}+L_{\tau})$ | $0$ | $0$ | $\hphantom{-}1$ | $\hphantom{-}1$ | $-\frac{1}{2}$ | $-\frac{1}{2}$
$B_{y}+L_{\mu}+L_{\tau}$ Ref. Farzan:2015doa  | $\frac{1}{3}$ | $\frac{1}{3}$ | $\hphantom{-}0$ | $\hphantom{-}0$ | $\hphantom{-}1$ | $\hphantom{-}1$
$L_{e}+2L_{\mu}+2L_{\tau}$ | $0$ | $0$ | $\hphantom{-}1$ | $\hphantom{-}1$ | $\hphantom{-}2$ | $\hphantom{-}2$
Table 1: Relevant charges for the matter effects in neutrino oscillation
experiments corresponding to a selection of models studied in the literature.
For the model presented in Ref. Farzan:2015doa  we have defined $B_{y}\equiv
B_{1}-yB_{2}-(3-y)B_{3}$ where $y$ is an arbitrary constant. The results
presented in this article are independent of $y$.

In general, the evolution of the neutrino and antineutrino flavor state during
propagation is governed by the Hamiltonian:

$H^{\nu}=H_{\text{vac}}+H_{\text{mat}}\quad\text{and}\quad
H^{\bar{\nu}}=(H_{\text{vac}}-H_{\text{mat}})^{*}\,,$ (4)

where $H_{\text{vac}}$ is the vacuum part which in the flavor basis
$(\nu_{e},\nu_{\mu},\nu_{\tau})$ reads

$H_{\text{vac}}=U_{\text{vac}}D_{\text{vac}}U_{\text{vac}}^{\dagger}\quad\text{with}\quad
D_{\text{vac}}=\frac{1}{2E_{\nu}}\mathop{\mathrm{diag}}(0,\Delta
m^{2}_{21},\Delta m^{2}_{31})\,.$ (5)

Here $U_{\text{vac}}$ denotes the three-lepton mixing matrix in vacuum
Pontecorvo:1967fh ; Maki:1962mu ; Kobayashi:1973fv . Following the convention
of Ref. Coloma:2016gei , we define
$U_{\text{vac}}=R_{23}(\theta_{23})R_{13}(\theta_{13})\tilde{R}_{12}(\theta_{12},\delta_{\text{CP}})$,
where $R_{ij}(\theta_{ij})$ is a rotation of angle $\theta_{ij}$ in the $ij$
plane and $\tilde{R}_{12}(\theta_{12},\delta_{\text{CP}})$ is a complex
rotation by angle $\theta_{12}$ and phase $\delta_{\text{CP}}$.

Concerning the matter part $H_{\text{mat}}$ of the Hamiltonian generated by
the SM together with the $U(1)^{\prime}$ interactions in (2), its form depends
on the new interaction length determined by the $Z^{\prime}$ mass as we
discuss next.

###  2.1 The large ($M_{Z^{\prime}}\gtrsim 10^{-12}$ eV) $M_{Z^{\prime}}$
limit: the NSI regime

In the limit of large $M_{Z^{\prime}}$, the $Z^{\prime}$ field can be
integrated out from the spectrum and Eq. (2) generate effective dimension-six
four-fermion interactions leading to Neutral Current NSI between neutrinos and
matter which are usually parametrized in the form of Eq. (1). The coefficients
$\varepsilon_{\alpha\beta}^{f,P}$ parametrizes the strength of the new
interaction with respect to the Fermi constant, with

$\varepsilon_{\alpha\beta}^{f,L}=\varepsilon_{\alpha\beta}^{f,R}=\delta_{\alpha\beta}\,\frac{1}{2\sqrt{2}G_{F}}\,\frac{g^{\prime
2}}{M_{Z^{\prime}}^{2}}\,a_{f}\,b_{\alpha}\equiv\delta_{\alpha\beta}\frac{1}{2}\,a_{f}\,b_{\alpha}\,\varepsilon^{0}$
(6)

where we have introduced the notation

$\varepsilon^{0}\equiv\frac{1}{\sqrt{2}G_{F}}\frac{g^{\prime
2}}{M_{Z^{\prime}}^{2}}$ (7)

As it is well known, only vector NSI contribute to the matter potential in
neutrino oscillations. It is therefore convenient to define the parameters
relevant for neutrino oscillation experiments as:

$\varepsilon_{\alpha\beta}^{f}\equiv\varepsilon_{\alpha\beta}^{f,L}+\varepsilon_{\alpha\beta}^{f,R}=\delta_{\alpha\beta}\,a_{f}\,b_{\alpha}\,\varepsilon^{0}\,.$
(8)

These interactions lead to a flavor diagonal modification of the matter
potential

$H_{\text{mat}}=\sqrt{2}G_{F}N_{e}(\vec{x})\begin{pmatrix}1+\mathcal{E}_{ee}(\vec{x})&0&0\\\
0&\mathcal{E}_{\mu\mu}(\vec{x})&0\\\
0&0&\mathcal{E}_{\tau\tau}(\vec{x})\end{pmatrix}$ (9)

where the “$+1$” term in the $ee$ entry accounts for the standard
contribution, and

$\mathcal{E}_{\alpha\alpha}(\vec{x})=\sum_{f=e,u,d}\frac{N_{f}(\vec{x})}{N_{e}(\vec{x})}\varepsilon_{\alpha\alpha}^{f}$
(10)

describes the non-standard part. Here $N_{f}(\vec{x})$ is the number density
of fermion $f$ at the position $\vec{x}$ along the neutrino trajectory. In Eq.
(10) we have limited the sum to the charged fermions present in ordinary
matter, $f=e,u,d$. Taking into account that
$N_{u}(\vec{x})=2N_{p}(\vec{x})+N_{n}(\vec{x})$ and
$N_{d}(\vec{x})=N_{p}(\vec{x})+2N_{n}(\vec{x})$, and also that matter
neutrality implies $N_{p}(\vec{x})=N_{e}(\vec{x})$, Eq. (10) becomes:

$\mathcal{E}_{\alpha\alpha}(\vec{x})=\big{(}\varepsilon_{\alpha\alpha}^{e}+\varepsilon_{\alpha\alpha}^{p}\big{)}+Y_{n}(\vec{x})\varepsilon_{\alpha\alpha}^{n}\quad\text{with}\quad
Y_{n}(\vec{x})\equiv\frac{N_{n}(\vec{x})}{N_{e}(\vec{x})}$ (11)

where

$\varepsilon_{\alpha\alpha}^{p}\equiv
2\varepsilon_{\alpha\alpha}^{u}+\varepsilon_{\alpha\alpha}^{d}=a_{p}\,b_{\alpha}\,\varepsilon^{0}\,,\qquad\varepsilon_{\alpha\alpha}^{n}\equiv
2\varepsilon_{\alpha\alpha}^{d}+\varepsilon_{\alpha\alpha}^{u}=a_{n}\,b_{\alpha}\,\varepsilon^{0}\,.$
(12)

and we have introduced the proton and neutron $Z^{\prime}$ couplings

$a_{p}\equiv 2a_{u}+a_{d}\,,\qquad a_{n}\equiv 2a_{d}+a_{u}\,.$ (13)

As discussed in Ref. GonzalezGarcia:2011my , in the Earth the neutron/proton
ratio $Y_{n}(\vec{x})$ which characterize the matter chemical composition can
be taken to be constant to very good approximation. The PREM model
Dziewonski:1981xy  fixes $Y_{n}=1.012$ in the Mantle and $Y_{n}=1.137$ in the
Core, with an average value $Y_{n}^{\oplus}=1.051$ all over the Earth. Setting
therefore $Y_{n}(\vec{x})\equiv Y_{n}^{\oplus}$ in Eqs. (10) and (11) we get
$\mathcal{E}_{\alpha\alpha}(\vec{x})\equiv\varepsilon_{\alpha\alpha}^{\oplus}$
with:

$\begin{split}\varepsilon_{\alpha\alpha}^{\oplus}&=\varepsilon_{\alpha\alpha}^{e}+\big{(}2+Y_{n}^{\oplus}\big{)}\varepsilon_{\alpha\alpha}^{u}+\big{(}1+2Y_{n}^{\oplus}\big{)}\varepsilon_{\alpha\alpha}^{d}=\big{(}\varepsilon_{\alpha\alpha}^{e}+\varepsilon_{\alpha\alpha}^{p}\big{)}+Y_{n}^{\oplus}\varepsilon_{\alpha\alpha}^{n}\\\
&=\big{[}(a_{e}+a_{p})+Y_{n}^{\oplus}a_{n}\big{]}\,b_{\alpha}\,\varepsilon^{0}\,.\end{split}$
(14)

For what concerns the study of propagation of solar and KamLAND neutrinos one
can work in the one mass dominance approximation, $\Delta m^{2}_{31}\to\infty$
(which effectively means that
$G_{F}N_{e}(\vec{x})\mathcal{E}_{\alpha\alpha}(\vec{x})\ll\Delta
m^{2}_{31}/E_{\nu}$). In this approximation the survival probability $P_{ee}$
can be written as Kuo:1986sk ; Guzzo:2000kx

$P_{ee}=c_{13}^{4}P_{\text{eff}}+s_{13}^{4}$ (15)

The probability $P_{\text{eff}}$ can be calculated in an effective $2\times 2$
model described by the Hamiltonian
$H_{\text{eff}}=H_{\text{vac}}^{\text{eff}}+H_{\text{mat,SM}}^{\text{eff}}+H_{\text{mat,$Z^{\prime}$}}^{\text{eff}}$,
with:

$\displaystyle H_{\text{vac}}^{\text{eff}}$ $\displaystyle=\frac{\Delta
m^{2}_{21}}{4E_{\nu}}\begin{pmatrix}-\cos
2\theta_{12}\,\hphantom{e^{-i\delta_{\text{CP}}}}&~{}\sin
2\theta_{12}\,e^{i\delta_{\text{CP}}}\\\ \hphantom{-}\sin
2\theta_{12}\,e^{-i\delta_{\text{CP}}}&~{}\cos
2\theta_{12}\,\hphantom{e^{i\delta_{\text{CP}}}}\end{pmatrix},$ (16)
$\displaystyle H_{\text{mat,SM}}^{\text{eff}}$
$\displaystyle=\sqrt{2}G_{F}N_{e}(\vec{x})\begin{pmatrix}c_{13}^{2}&0\\\
0&0\end{pmatrix}$ (17)

and

$H_{\text{mat,Z'}}^{\text{eff}}=\sqrt{2}G_{F}N_{e}(\vec{x})\,\varepsilon^{0}\,\big{[}a_{e}+a_{p}+Y_{n}(\vec{x})a_{n}\big{]}\begin{pmatrix}-b_{D}&b_{N}\\\
\hphantom{+}b_{N}&b_{D}\end{pmatrix}$ (18)

where

$\displaystyle\begin{split}b_{D}&=-\frac{c_{13}^{2}}{2}\left(b_{e}-b_{\mu}\right)+\frac{s_{23}^{2}-s_{13}^{2}c_{23}^{2}}{2}\left(b_{\tau}-b_{\mu}\right)\,,\end{split}$
(19) $\displaystyle b_{N}$
$\displaystyle=s_{13}c_{23}s_{23}\left(b_{\tau}-b_{\mu}\right)\,.$ (20)

Following Ref. Esteban:2018ppq  we can rewrite the $Z^{\prime}$ contribution
as:

$\displaystyle H_{\text{mat,Z'}}^{\text{eff}}$
$\displaystyle=\sqrt{2}G_{F}N_{e}(\vec{x})\big{[}\cos\eta+Y_{n}(\vec{x})\sin\eta\big{]}\begin{pmatrix}-\varepsilon_{D}^{\eta\hphantom{*}}&\varepsilon_{N}^{\eta}\\\
\hphantom{+}\varepsilon_{N}^{\eta}&\varepsilon_{D}^{\eta}\end{pmatrix},$ (21)

where the angle $\eta$ parametrizes the ratio of the charges of the matter
particles as:

$\cos\eta=\frac{a_{e}+a_{p}}{\sqrt{(a_{e}+a_{p})^{2}+a_{n}^{2}}}\,,\qquad\sin\eta=\frac{a_{n}}{\sqrt{(a_{e}+a_{p})^{2}+a_{n}^{2}}}\,,$
(22)

and

$\varepsilon_{D,N}^{\eta}=\sqrt{(a_{e}+a_{p})^{2}+a_{n}^{2}}\,b_{D,N}\,\varepsilon^{0}\,.$
(23)

The neutrino oscillation phenomenology in this regime reduces to a special
subclass of the general NSI interactions analyzed in Ref. Esteban:2018ppq
.111To be precise, the data analysis performed in Ref. Esteban:2018ppq  was
restricted to NSI with quarks, ie $a_{e}=0$. The formalism for matter effects
can be trivially extended to NSI coupled to electrons as shown above. However,
NSI coupled to electrons would affect not only neutrino propagation in matter
as described, but also the neutrino-electron (ES) scattering cross-section in
experiments such as SK, SNO and Borexino. In order to keep the analysis
manageable, in Ref. Esteban:2018ppq , and in what follows, the NSI corrections
to the ES scattering cross section in SK, SNO, and Borexino are neglected. In
the absence of cancellations between propagation and interaction effects this
renders the results of the oscillation analysis conservative. In particular,
as a consequence of the CPT symmetry (see also Refs. GonzalezGarcia:2011my ;
Gonzalez-Garcia:2013usa ; Coloma:2016gei  for a discussion in the context of
NSI) the neutrino evolution is invariant if the relevant Hamiltonian is
transformed as $H\to-H^{*}$. In vacuum this transformation can be realized by
changing the oscillation parameters as

$\displaystyle\Delta m^{2}_{31}$ $\displaystyle\to-\Delta m^{2}_{31}+\Delta
m^{2}_{21}=-\Delta m^{2}_{32}\,,$ (24) $\displaystyle\theta_{12}$
$\displaystyle\to\pi-\theta_{12}\,,$ $\displaystyle\delta_{\text{CP}}$
$\displaystyle\to\pi-\delta_{\text{CP}}\,,$

where $\delta_{\text{CP}}$ is the leptonic Dirac CP phase, and we are using
here the parameterization conventions from Refs. Esteban:2018ppq ;
Coloma:2016gei . The symmetry is broken by the standard matter effect, which
allows a determination of the octant of $\theta_{12}$ and (in principle) of
the sign of $\Delta m^{2}_{31}$. However, in the presence of the
$Z^{\prime}$-induced NSI, the symmetry can be restored if in addition to the
transformation Eq. (24), the $\mathcal{E}_{\alpha\alpha}(\vec{x})$ terms can
be transformed as Gonzalez-Garcia:2013usa ; Bakhti:2014pva ; Coloma:2016gei

$\displaystyle\big{[}\mathcal{E}_{ee}(\vec{x})-\mathcal{E}_{\mu\mu}(\vec{x})\big{]}$
$\displaystyle\to-\big{[}\mathcal{E}_{ee}(\vec{x})-\mathcal{E}_{\mu\mu}(\vec{x})\big{]}-2\,,$
(25)
$\displaystyle\big{[}\mathcal{E}_{\tau\tau}(\vec{x})-\mathcal{E}_{\mu\mu}(\vec{x})\big{]}$
$\displaystyle\to-\big{[}\mathcal{E}_{\tau\tau}(\vec{x})-\mathcal{E}_{\mu\mu}(\vec{x})\big{]}\,.$

Eq. (24) shows that this degeneracy implies a change in the octant of
$\theta_{12}$ (as manifest in the LMA-D fit to solar neutrino data
Miranda:2004nb ) as well as a change in the neutrino mass ordering, i.e., the
sign of $\Delta m^{2}_{31}$. For that reason it has been called “generalized
mass ordering degeneracy” in Ref. Coloma:2016gei . Because of the position
dependence of the NSI hamiltonian described by
$\mathcal{E}_{\alpha\alpha}(\vec{x})$ this degeneracy is only approximate,
mostly due to the non-trivial neutron/proton ratio along the neutrino path
inside the Sun. In what follows when marginalizing over $\theta_{12}$ we
consider two distinct parts of the parameter space: one with
$\theta_{12}<45^{\circ}$, which we denote as LIGHT, and one with
$\theta_{12}>45^{\circ}$, which we denote by DARK.

Apart from the appearance of this degenerate solution, another feature to
consider in the global analysis of oscillation data in presence of NSI is the
possibility to further improve the quality of the fit with respect to that of
_standard_ $3\nu$ oscillations in the LIGHT sector. Till recently this was
indeed the case because for the last decade the value of $\Delta m^{2}_{21}$
preferred by KamLAND was somewhat higher than the one from solar experiments.
This tension appeared due to a combination of two effects: the fact that the
$\mbox{}^{8}$B measurements performed by SNO, SK and Borexino showed no
evidence of the low energy spectrum turn-up expected in the standard LMA-MSW
Wolfenstein:1977ue ; Mikheev:1986gs  solution for the value of $\Delta
m^{2}_{21}$ favored by KamLAND, and the observation of a non-vanishing day-
night asymmetry in SK, whose size was larger than the one predicted for the
$\Delta m^{2}_{21}$ value indicated by KamLAND. With the data included in the
analysis in Ref. Esteban:2018azc ; Esteban:2018ppq  this resulted into a
tension of $\Delta\chi^{2}\sim 7.4$ for the standard $3\nu$ oscillations. Such
tension could be alleviated in presence of a non-standard matter potential,
thus leading to a possible decrease in the minimum $\chi^{2}$. However, with
the latest 2970-days SK4 results presented at the Neutrino2020 conference
SK:nu2020  in the form of total energy spectrum and updated day-night
asymmetry, the tension between the best fit $\Delta m^{2}_{21}$ of KamLAND and
that of the solar results has decreased. Currently they are compatible within
$1.1\sigma$ in the latest global analysis Esteban:2020cvm .

###  2.2 The finite $M_{Z^{\prime}}$ case: the long-range interaction regime

If $M_{Z^{\prime}}$ is very light the four-fermion contact interaction
approximation in Eq. (1) does not hold and the potential encountered by the
neutrino in its trajectory depends on the integral of the source density
within a radius $\sim 1/M_{Z^{\prime}}$ around it. However, following Ref.
GonzalezGarcia:2006vp  the generalized matter potential can still be written
as Eq. (9) provided that Eq. (10) is modified as:

$\mathcal{E}_{\alpha\alpha}(\vec{x})=\sum_{f=e,u,d}\frac{\hat{N}_{f}(\vec{x},M_{Z^{\prime}})}{N_{e}(\vec{x})}\varepsilon_{\alpha\alpha}^{f}$
(26)

where

$\hat{N}_{f}(\vec{x},M_{Z^{\prime}})\equiv\frac{4\pi}{M_{Z^{\prime}}^{2}}\int
N_{f}(\vec{\rho})\,\frac{e^{-M_{Z^{\prime}}|\vec{\rho}-\vec{x}|}}{|\vec{\rho}-\vec{x}|}\,d^{3}\vec{\rho}\,.$
(27)

Taking into account that ordinary matter is neutral and only contains
$f=e,u,d$, we can rewrite Eq. (26) in a way that generalizes Eq. (11):

$\mathcal{E}_{\alpha\alpha}(\vec{x})=F_{e}(\vec{x},M_{Z^{\prime}})\big{(}\varepsilon_{\alpha\alpha}^{e}+\varepsilon_{\alpha\alpha}^{p}\big{)}+F_{n}(\vec{x},M_{Z^{\prime}})Y_{n}(\vec{x})\varepsilon_{\alpha\alpha}^{n}\\\
\text{with}\quad
F_{i}(\vec{x},M_{Z^{\prime}})\equiv\frac{\hat{N}_{i}(\vec{x},M_{Z^{\prime}})}{N_{i}(\vec{x})}\quad\text{and}\quad
i\in\\{e,n\\}\,.$ (28)

For what concerns neutrinos traveling inside the Sun, the propagation effects
induced by the new interactions are completely dominated by the solar matter
distribution. Denoting by $\vec{x}_{\odot}$ the center of the Sun and
accounting for the spherical symmetry of the matter potential we can write:

$F_{i}(\vec{x},M_{Z^{\prime}})\simeq
F_{i}^{\odot}(|\vec{x}-\vec{x}_{\odot}|,M_{Z^{\prime}})\\\ \text{with}\quad
F_{i}^{\odot}(r,M_{Z^{\prime}})=\frac{1}{N_{i}^{\odot}(r)}\cdot\frac{M_{Z^{\prime}}}{2\,r}\int_{0}^{R_{\odot}}\rho\,N_{i}^{\odot}(\rho)\big{[}e^{-M_{Z^{\prime}}|\rho-r|}-e^{-M_{Z^{\prime}}(\rho+r)}\big{]}\,d\rho\,.$
(29)

A similar formula can be derived for neutrinos traveling inside the Earth, but
in this case the effective potential has an extra term induced by the Sun
matter density. Concretely, denoting by $\vec{x}_{\oplus}$ the center of the
Earth and by $X_{\ominus}=|\vec{x}_{\odot}-\vec{x}_{\oplus}|$ the Sun-Earth
distance, we have:

$F_{i}(\vec{x},M_{Z^{\prime}})\simeq
F_{i}^{\oplus}(|\vec{x}-\vec{x}_{\oplus}|,M_{Z^{\prime}})\\\ \text{with}\quad
F_{i}^{\oplus}(r,M_{Z^{\prime}})=\frac{1}{N_{i}^{\oplus}(r)}M_{Z^{\prime}}\Bigg{\\{}\frac{1}{2\,r}\int_{0}^{R_{\oplus}}\rho\,N_{i}^{\oplus}(\rho)\big{[}e^{-M_{Z^{\prime}}|\rho-r|}-e^{-M_{Z^{\prime}}(\rho+r)}\big{]}\,d\rho\\\
+\frac{e^{-M_{Z^{\prime}}X_{\ominus}}}{X_{\ominus}}\int_{0}^{R_{\odot}}\rho\,N_{i}^{\odot}(\rho)\sinh(\rho\,M_{Z^{\prime}})\,d\rho\Bigg{\\}}\,.$
(30)

The solar-induced contribution becomes non-negligible when the range of the
interactions, $1\big{/}M_{Z^{\prime}}$, is comparable or larger than the Sun-
Earth distance $X_{\ominus}$.

The factors $F_{i}^{\odot}(r,M_{Z^{\prime}})$ and
$F_{i}^{\oplus}(r,M_{Z^{\prime}})$ represent the modification of the matter
potential due to the finite range of the interaction mediated by the
$Z^{\prime}$ with respect to that obtained in the contact interaction limit.
Such limit is recovered when the range of the new interactions become shorter
than the typical size of the matter distribution, i.e., $R_{\odot(\oplus)}$.
Hence:

$F_{i}^{\odot(\oplus)}(r,M_{Z^{\prime}})\to 1\quad\text{for}\quad
M_{Z^{\prime}}\gg 1/R_{\odot(\oplus)}\,.$ (31)

For solar neutrinos further simplification follows if one takes into account
that for adiabatic transitions the dominant matter effects is generated by the
potential at the neutrino production point which is is close to the Sun
center. So to a very good approximation one can scale the contact interaction
potential with a position independent factor
$F_{i}^{\odot}(0,M_{Z^{\prime}})$. For the Earth matter potential the position
dependence of the factor $F_{i}^{\oplus}(r,M_{Z^{\prime}})$ is very weak in
the current experiments, so one can also scale the contact interaction
potential with an approximate $F_{i}^{\oplus}(\bar{r},M_{Z^{\prime}})$
evaluated at a fix $\bar{r}$ which we take to be also $\bar{r}=0$.

In Fig. 1 we plot these scale factors
$F_{e}^{\odot(\oplus)}(0,M_{Z^{\prime}})$. As seen in the figure for
$M_{Z^{\prime}}\lesssim 10^{-13}$ eV the matter potential in the Earth is more
suppressed with respect to that in the Sun. In principle, this opens the
possibility of configurations for which the $U(1)^{\prime}$-induced matter
potential in the Sun is large enough without conflicting with bounds imposed
by atmospheric and long-baseline experiments. This also implies that in the
combined analysis of solar and KamLAND data, for a given value of
$M_{Z^{\prime}}$, the effective matter potential for solar neutrinos will be
suppressed by a different factor than that for KamLAND antineutrinos. To
illustrate the overall $M_{Z^{\prime}}$ dependence of the effect we show in
Fig. 1 the effective suppression factor in the combined solar+KamLAND analysis
calculated by scaling the results obtained for a specific model (concretely,
for a $Z^{\prime}$ coupled to $L_{e}-L_{\mu}$, but the results are similar for
models with other couplings) for each $M_{Z^{\prime}}$ to those obtained in
the large $M_{Z^{\prime}}$ regime. As seen in the figure for
$M_{Z^{\prime}}\gtrsim 10^{-10}$ eV both the effective potential in the
Solar+KamLAND analysis and the Earth matter potential relevant for atmospheric
and LBL neutrinos are well within the contact interaction regime. Conversely
for $M_{Z^{\prime}}\lesssim 10^{-13}$ eV all the matter potentials in the
analysis show deviations from the contact interaction regime.

Figure 1: Effective suppression factor of the matter potential due to the long
range of the $U(1)^{\prime}$ interactions as a function of $M_{Z^{\prime}}$.
The red (blue) curve corresponds to the potential in the Sun (Earth). The
purple line corresponds to the effective combined suppression factor in the
analysis of Solar+KamLAND data (see text for details).

Concerning the phenomenology of neutrino oscillations in the presence of the
modified matter potential in this long-range interaction regime, the main
difference with the NSI contact interaction case is the impossibility of
realizing the “generalized mass ordering degeneracy” in Eqs. (24) and (25)
because of the very different $\vec{x}$ dependence of the SM matter potential
and the one generated by the $Z^{\prime}$. In other words, one cannot “flip
the sign of the matter hamiltonian” by adding to the standard $N_{e}(x)$
something which has a completely different $\vec{x}$ profile. Hence, there is
no LMA-D solution for these models. On the other hand, it is still possible,
at least in principle, that a long-range potential leads to an improvement on
the fit to solar and KamLAND data with respect to the pure LMA solution.

##  3 Results of the global oscillation analysis

We have performed a global fit to neutrino oscillation data in the framework
of $3\nu$ massive neutrinos with new neutrino-matter interactions generated by
$U(1)^{\prime}$ models and characterized by the Lagrangian in Eq. (2). For the
detailed description of methodology and data included we refer to the
comprehensive global fit in Ref. Esteban:2018ppq  performed in the framework
of three-flavor oscillations plus NSI. In addition in the present analysis we
account for the latest LBL data samples included in NuFIT-5.0 Esteban:2020cvm
which includes the previously cited solar neutrinos 2970-days SK4 results
SK:nu2020 , the updated medium baseline reactor samples from RENO RENO:nu2020
and Double Chooz DoubleC:nu2020 , and the latest long-baseline samples from
T2K T2K:nu2020  and NO$\nu$A NOvA:nu2020 . Notice that in order to keep the
fit manageable we proceed as in in Ref. Esteban:2018ppq  and restrict
ourselves to the CP-conserving case and set
$\delta_{\text{CP}}\in\\{0,\pi\\}$. Consequently the T2K and NO$\nu$A
appearance data (which exhibit substantial dependence on the leptonic CP
phase) are not included in the fit. With these data we construct a global
$\chi^{2}$ function:

$\chi^{2}_{\text{OSC+Z'}}(g^{\prime},M_{Z^{\prime}}|\vec{\omega})$ (32)

for each $Z^{\prime}$ model. In general each model belongs to a family
characterized by a set of $U(1)^{\prime}$ charges; for each family
$\chi^{2}_{\text{OSC+Z'}}$ depends on the two variables parametrizing the new
interaction, $g^{\prime}$ and $M_{Z^{\prime}}$, plus the six oscillation
parameters $\vec{\omega}\equiv(\Delta m^{2}_{21},\Delta
m^{2}_{31},\theta_{12},\theta_{13},\theta_{23},\delta_{\text{CP}})$.

Following the discussion in Sec. 2, we have performed the analysis in two
physically distinctive domains of $M_{Z^{\prime}}$:

  * •

NSI Domain (DOM=NSI): In this case, the range of the induced non-standard
interactions in both the Sun and the Earth matter is short enough for the
four-fermion effective description to hold. As seen in Fig. 1 this happens for
$M_{Z^{\prime}}\gtrsim 10^{-10}$ eV. In this regime a possible conflict with
the cosmological bound on $\Delta N_{\text{eff}}$ may appear because of the
contribution of either the $Z^{\prime}$ itself (if lighter than all active
neutrinos) or by the extra contribution to the neutrino density produced by
the $Z^{\prime}$ decay (if heavier than some $\nu$ mass eigenstate). These
bounds can be evaded in two distinct ranges of the $U(1)^{\prime}$
interactions:

    * –

$M_{Z^{\prime}}\gtrsim 5$ MeV for which the contribution to the neutrino
energy density due to the decay of the $Z^{\prime}$ into neutrinos is
sufficiently suppressed by the Boltzmann factor $\sim\exp(-M_{Z^{\prime}}/T)$
Kamada:2015era ;

    * –

$M_{Z^{\prime}}\lesssim\mathcal{O}(\text{eV})$ but with very weak coupling
$g^{\prime}<10^{-10}$ for which the $Z^{\prime}$ is produced through freeze-
in. In this regime, even if $Z^{\prime}$ could decay to the lightest neutrinos
this would happen after neutrino decoupling making the contribution to $\Delta
N_{\text{eff}}$ negligible or at most within the present allowed range
Escudero:2019gzq .

  * •

Long-Range Domain (DOM=LRI): If $M_{Z^{\prime}}\lesssim 10^{-13}$ eV
interactions in the Earth and Sun matter are long range, and for a given value
of $g^{\prime}$ and $M_{Z^{\prime}}$ the effects in the Earth are suppressed
with respect to those in the Sun. As mentioned above, for
$g^{\prime}<10^{-10}$ the contribution to $\Delta N_{\text{eff}}$ is
negligible.

Correspondingly we define

$\chi^{2}_{\text{OSC+Z',DOM}}(g^{\prime},M_{Z^{\prime}}|\vec{\omega})\equiv\chi^{2}_{\text{OSC+Z'}}(g^{\prime},M_{Z^{\prime}}\in\text{DOM}|\vec{\omega})$
(33)

In both domains we compare the results of the fit including the new
$U(1)^{\prime}$ interaction with those obtained in the “standard”
$3\nu$-mixing scenario, which we will denote as “OSC”, and for which the
present global fit yields

$\chi^{2}_{\text{OSC,min}}=718.5\,.$ (34)

We will classify the models according to the quality of the fit in the
presence of the $U(1)^{\prime}$ interactions compared with that of OSC by
defining

$\Delta\chi^{2}_{\text{LIGHT,DOM}}(g^{\prime},M_{Z^{\prime}})\equiv\left.\chi^{2}_{\text{OSC+Z',DOM}}(g^{\prime},M_{Z^{\prime}}|\vec{\omega})\right|_{\text{marg,LIGHT}}-\chi^{2}_{\text{OSC,min}}\,,$
(35)

where by $\left.\right|_{\text{marg,LIGHT}}$ we imply that the minimization
over the oscillation parameters is done in the LIGHT sector of parameter
space. In addition the presence of a viable LMA-D solution can be quantified
in terms of

$\Delta\chi^{2}_{\text{DARK,DOM}}(g^{\prime},M_{Z^{\prime}})\equiv\left.\chi^{2}_{\text{OSC+Z',DOM}}(g^{\prime},M_{Z^{\prime}}|\vec{\omega})\right|_{\text{marg,
DARK}}-\chi^{2}_{\text{OSC,min}}\,,$ (36)

where by $\left.\right|_{\text{marg,DARK}}$ we imply that the minimization
over the oscillation parameters is done in the DARK sector.

We have surveyed the model space by performing the global oscillation analysis
for a grid of $U(1)^{\prime}$ interactions characterized by the six couplings
$a_{u,d}\in\\{-1,0,1\\}$ and
$a_{e},\,b_{e,\mu,\tau}\in\\{-3,-2,-1,0,1,2,3\\}$. In this way our survey
covers a total of $\sim 10000$ different sets of $U(1)^{\prime}$ charges which
can produce effects in matter propagation in neutrino oscillation experiments.

###  3.1 Bounds

We first search for models for which in the LIGHT sector the new interactions
lead to a significantly better fit of the oscillation data compared to pure
oscillations for some value of $g^{\prime}$ and $M_{Z^{\prime}}$. We find that
in the NSI (LRI) domain 88% (90%) of the surveyed sets of charges lead to a
decrease in the $\chi^{2}$ of the global analysis when compared to the
standard oscillation case. The percentages grow to 100% when restricting to
the subclass of anomaly free vector $Z^{\prime}$ models with gauging of SM
global symmetries with SM plus right-handed neutrinos matter content of Eq.
(3). However the improvement in the quality of the fit is never statistically
significant. As an illustration we show in the second column of tables 2 and 3
the minimum values of $\Delta\chi^{2}_{\text{LIGHT,DOM}}$ for $U(1)^{\prime}$
interactions characterized by the specific set of charges in table 1.
Comparing the two tables we notice that for some of the cases the fit can be
slightly better in the LRI domain than in the NSI domain but still below the
$2\sigma$ level.

Model | $(\Delta\chi_{\text{LIGHT,LRI}}^{2})_{\text{min}}$ |  $g^{\prime}\leq$ bound
---|---|---
$B-3L_{e}$ | $-1.4$ | $6.6\times 10^{-27}$
$B-3L_{\mu}$ | $-1.1$ | $7.0\times 10^{-27}$
$B-3L_{\tau}$ | $-1.8$ | $7.3\times 10^{-27}$
$B-\frac{3}{2}(L_{\mu}+L_{\tau})$ | $-1.2$ | $7.2\times 10^{-27}$
$L_{e}-L_{\mu}$ | $-1.3$ | $9.7\times 10^{-27}$
$L_{e}-L_{\tau}$ | $-1.7$ | $1.0\times 10^{-26}$
$L_{e}-\frac{1}{2}(L_{\mu}+L_{\tau})$ | $-1.4$ | $9.8\times 10^{-27}$
$B_{y}+L_{\mu}+L_{\tau}$ Ref. Farzan:2015doa  | $\hphantom{-}0\hphantom{.0}$ | $4.9\times 10^{-27}$
$L_{e}+2L_{\mu}+2L_{\tau}$ | $\hphantom{-}0\hphantom{.0}$ | $6.0\times 10^{-27}$
Table 2: Results for the models with charges in table 1 in the LRI regime. For
the model in Ref. Farzan:2015doa  we have defined $B_{y}\equiv
B_{1}-yB_{2}-(3-y)B_{3}$. The second column gives minimum $\Delta\chi^{2}$
defined w.r.t. the $3\nu$ oscillation (see Eq. (35)). The last column gives
the the upper bound for the coupling of asymptotically for ultra light
mediators, $M_{Z^{\prime}}\lesssim 10^{-15}$ eV

Quantitatively we find that in the NSI regime none of the surveyed models
yields an improvement beyond -1.9 units of $\chi^{2}$. This is the case, for
example, of a $Z^{\prime}$ coupled to $B-3L_{e}+2L_{\mu}+3L_{\tau}$.
Generically models in the LRI domain can provide a better fit with a reduction
of up to 3.5 units of $\chi^{2}$. For example a model with charge
$L_{e}+2L_{\mu}-3L_{\tau}$ and a $Z^{\prime}$ with $M_{Z^{\prime}}\sim 5\times
10^{-15}$ eV provides a better fit than standard oscillations by
$(\Delta\chi^{2}_{\text{LIGHT,LRI}})_{\text{min}}=-3.2$.

Model | $(\Delta\chi_{\text{LIGHT,NSI}}^{2})_{\text{min}}$ |  $g^{\prime}\leq$ bound $\bigg{(}\dfrac{M_{Z^{\prime}}}{\text{100~{}MeV}}\bigg{)}$
---|---|---
$B-3L_{e}$ | $-1.4$ | $2.0\times 10^{-4}$
$B-3L_{\mu}$ | $\hphantom{-}0\hphantom{.0}$ | $4.6\times 10^{-5}$
$B-3L_{\tau}$ | $-0.6$ | $4.7\times 10^{-5}$
$B-\frac{3}{2}(L_{\mu}+L_{\tau})$ | $-1.1$ | $2.2\times 10^{-4}$
$L_{e}-L_{\mu}$ | $-1.3$ | $1.2\times 10^{-4}$
$L_{e}-L_{\tau}$ | $-1.0$ | $1.2\times 10^{-4}$
$L_{e}-\frac{1}{2}(L_{\mu}+L_{\tau})$ | $-1.3$ | $3.0\times 10^{-4}$
$B_{y}+L_{\mu}+L_{\tau}$ Ref. Farzan:2015doa  | $\hphantom{-}0\hphantom{.0}$ | $1.5\times 10^{-4}$
$L_{e}+2L_{\mu}+2L_{\tau}$ | $-0.1$ | $1.8\times 10^{-4}$
Table 3: Results for the models with charges in table 1 in the NSI regime. For
the model in Ref. Farzan:2015doa  we have defined $B_{y}\equiv
B_{1}-yB_{2}-(3-y)B_{3}$. The second column gives minimum $\Delta\chi^{2}$
defined w.r.t. the $3\nu$ oscillation (see Eq. (35)). The last column gives
the coefficient of the bound on the coupling over the mediator mass in units
of 100 MeV.

But in summary, our analysis shows that none of the set of charges surveyed in
both NSI or LRI domains improved over standard oscillations at the $2\sigma$
level, this is
$\min_{g^{\prime},M_{Z^{\prime}}}(\Delta\chi^{2}_{\text{LIGHT,DOM}})$ was
always larger than $-4$.222We notice that before the new results from Super-
Kamiokande SK:nu2020  there were models which could improve the mismatch
between the best fit $\Delta m^{2}_{21}$ in Solar and KamLAND. For such
models, one could find values of $g^{\prime}$ and $M_{Z^{\prime}}$ for which
the fit was better than standard oscillations by more than 4 units of
$\chi^{2}$. Consequently for all models surveyed one can conclude that the
analysis of neutrino oscillation experiments show no significant evidence of
$U(1)^{\prime}$ interactions. Consequently one can exclude models at a certain
confidence level – which we have chosen to be 95.45% – by verifying that their
global fit is worse than in OSC by the corresponding units of $\chi^{2}$ (4
units for 95.45% CL), this is:

$\Delta\chi^{2}_{\text{LIGHT,DOM}}(g^{\prime},M_{Z^{\prime}})>4.$ (37)

Let us stress that with the above condition we are not “deriving two-
dimensional excluded regions in the parameter space”, but we are instead
determining the values of $g^{\prime}$ and $M_{Z^{\prime}}$ for which the
$U(1)^{\prime}$ model characterized by such interaction strength and
interaction length, gives a fit which is worse than standard oscillations by
at least 4 units of $\chi^{2}$. As in the LIGHT sector the standard model is
recovered for either $g^{\prime}\to 0$ or $M_{Z^{\prime}}\to\infty$, the above
condition yields also the $2\sigma$ excluded one-dimensional upper range of
interaction coupling $g^{\prime}$ for each value of the interaction length
(or, correspondingly, the $2\sigma$ excluded one-dimensional lower range of
$M_{Z^{\prime}}$ for each value of the interaction coupling), for all models
characterized by a given set of charges.

Figure 2: Values of $g^{\prime}$ and $M_{Z^{\prime}}\leq\mathcal{O}$(eV) for
which a $U(1)^{\prime}$ model coupled to the charges labeled in each panel
gives a worse fit than standard oscillation by $2\sigma$, Eq. (37) (hatched
region). We also show the bounds on these models imposed by gravitational
fifth force searches Salumbides:2013dua ; Adelberger:2009zz  and by
equivalence principle tests Schlamminger:2007ht . In the window corresponding
to the model in Ref. Farzan:2015doa  we have defined $B_{y}\equiv
B_{1}-yB_{2}-(3-y)B_{3}$. See text for details.

The corresponding excluded ranges for the $Z^{\prime}$ coupling and mass for
the models with couplings listed in table 1 are shown in Figs. 2 and 3 for
$M_{Z^{\prime}}$ below 1 eV and in the $\mathcal{O}(\text{MeV--GeV})$ range,
respectively333Tables with the numerical values of the bounds can be provided
upon request to the authors. – i.e., below and above the window strongly
disfavored by the cosmological bound on $\Delta N_{\text{eff}}$. In particular
in Fig. 2 we observe the slope change of the oscillation exclusion ranges for
masses $M_{Z^{\prime}}\sim 10^{-15}$–$10^{-13}$ GeV, for which the interaction
length is longer than the Earth and Sun radius and the matter potential in the
Earth and in the Sun becomes saturated (see Fig. 1). Quantitatively, for the
models with charges in table 1 we find that in the LRI domain the analysis of
oscillation data yields the upper bound for the coupling of asymptotically
ultra light mediators, $M_{Z^{\prime}}\lesssim 10^{-15}$ eV which we list in
table 2. For the sake of comparison, we also show in Fig. 2 the bounds on
these models imposed by gravitational fifth force searches, and by equivalence
principle tests. Those are the strongest model independent constraints in the
shown range of $Z^{\prime}$ mass and coupling derived at comparable confidence
level under the minimal assumptions in Eq. (2). For some of the models shown
in the figure, additional bounds on this range of masses can arise from
cosmological and astrophysical observations, including constraints on
invisible neutrino decay $\nu_{a}\to\nu_{b}Z^{\prime}$ from Cosmic Microwave
Background (CMB) data Escudero:2019gfk , bounds from Black-Hole superradiance
Baryakhtar:2017ngi , or from flavor composition of extra-galactic neutrinos
Bustamante:2018mzu . All of them, however, largely depend on the assumptions
made, and can be evaded in specific scenarios. In addition bounds from
production of the light $Z^{\prime}$ in meson decays, neutrinoless
$\beta\beta$ decay, or neutrino annihilation in Big-Bang Nucleosynthesis (BBN)
and supernovae are relevant for $Z^{\prime}$ with lower masses than those
shown in the figure (see, for example, Ref. Dror:2020fbh ).

The bounds imposed by gravitational fifth force searches shown in Fig. 2 were
obtained by rescaling the results shown in Ref. Salumbides:2013dua  (which, in
turn, were recasted from Ref. Adelberger:2009zz ). Limits from equivalence
principle tests are obtained rescaling the results from Ref.
Schlamminger:2007ht . Details of the rescaling of the published bounds applied
for the specific models can be found in the Appendix. It is important to
notice that these exclusion regions obtained by recasting the boundaries of
the published regions may not correspond to the statistical condition we
employed, Eq. (37). So the comparison has to be taken with a pinch of salt.
Still, from the figures we see that, generically, for all models shown the
oscillation analysis improves over the existing bounds for $Z^{\prime}$
lighter than $\sim 10^{-8}$ eV or $\sim 10^{-11}$ eV depending on whether the
$U(1)^{\prime}$ current involves coupling to electron lepton number.

Figure 3: Values of $g^{\prime}$ and $\text{5 MeV}\leq
M_{Z^{\prime}}\leq\text{10~{}GeV}$ for which a $U(1)^{\prime}$ model coupled
to the charges labeled in each panel gives a worse fit than standard
oscillation by $2\sigma$, Eq. (37) (hatched region). We also show the bounds
on these models imposed by a compilation of experiments as labeled in the
figure. In the window corresponding to the model in Ref. Farzan:2015doa  we
have defined $B_{y}\equiv B_{1}-yB_{2}-(3-y)B_{3}$. See text for details.

Conversely the results shown in Fig. 3 for $M_{Z^{\prime}}$ in the
$\mathcal{O}(\text{MeV--GeV})$ range correspond to $U(1)^{\prime}$ effects in
oscillation experiments always in the NSI domain. In this case, the analysis
of oscillation data results in a bound on $g^{\prime}$ which scales as the
inverse of the mediator mass with coefficients which we list in table 3. For
the sake of comparison, we also show in Fig. 3 a compilation of the most
relevant experimental bounds on these $U(1)^{\prime}$ models. These include
constraints from electron and proton fixed target experiments, neutrino
electron elastic scattering, coherent neutrino nucleus elastic scattering,
white dwarf cooling and collider constraints. Let us point out that we have
neglected kinetic mixing so far because it does not affect the oscillation
bounds as previously discussed. But it should be kept in mind that the
presence of kinetic mixing could either strengthen or weaken the bounds
reported from other experiments. Appendix A contains all relevant details on
the derivation of these bounds.

As seen in Fig. 3 for several models there are values of $g^{\prime}$ and
$M_{Z^{\prime}}$ which are only constrained by the oscillation analysis. This
is particularly the case for $U(1)^{\prime}$ coupled to $B-3L_{\tau}$. We have
found no competitive bound from other experiments in the shown window with the
only exception of the constraint reported in Ref. Kling:2020iar  using data
from DONUT Kodama:2007aa  and from NA62 (for bounds from other experiments
relevant for larger couplings see for example Refs. Farzan:2016wym ;
Heeck:2018nzc ).444Potentially, for the model in Ref. Farzan:2015doa , coupled
to $B_{1}-yB_{2}-(3-y)B_{3}+L_{\mu}+L_{\tau}$ there are additional strong
bounds associated to quark mixing effects due to the non-conservation of the
$U(1)^{\prime}$ current by the SM particles (see Refs. Dror:2017ehi ;
Dror:2017nsg ; Dror:2018wfl ). According to Ref. Farzan:2015doa  such bounds
could be mitigated for a particular choice of $y$ and with the inclusion of a
$U(1)^{\prime}$-charged scalar sector. But strictly speaking a careful
evaluation of these bounds has not been presented in the literature, and it is
beyond the scope of this paper.

###  3.2 Models for LMA-D

A subset of models can lead to an allowed region in the DARK sector with LMA-D
within 4 units of $\chi^{2}$ with respect to the standard OSC solution,

$\Delta\chi^{2}_{\text{DARK,DOM}}(g^{\prime},M_{Z^{\prime}})<4\,.$ (38)

As discussed in Sec. 2.2 this can only happen in the NSI domain. In our survey
we have found that 4.8% of the set of charges studied can have a best fit in
LMA-D, and 5.2% lead to LMA-D verifying Eq. (38). However none of these set of
charges correspond to the subclass of anomaly free vector $Z^{\prime}$ models
with gauging of SM global symmetries with SM plus right-handed neutrinos
matter content (Eq. (3)).

In particular for the first seven set of models in table 1 the LMA-D solutions
lies at more than $5\sigma$ from the standard oscillation fit. On the contrary
the models in the last two lines yield a viable LMA-D solution The first one
was proposed in Ref. Farzan:2015doa  precisely as a viable model for LMA-D.
Indeed in this case we find
$(\Delta\chi^{2}_{\text{DARK,NSI}})_{\text{min}}=1.2$, which is within 4 units
of $\chi^{2}$ from the pure oscillation result but the best fit for this model
charges still lies within the LIGHT sector. We also show the results for a
$U(1)^{\prime}$ model coupled to $L_{e}+2L_{\mu}+2L_{\tau}$ for which we find
that the best fit is LMA-D with
$(\Delta\chi^{2}_{\text{DARK,NSI}})_{\text{min}}=-1.3$.

In Fig. 4 we plot as black bands the range of coupling and masses for these
two set of charges verifying the condition (38), together with the compilation
of relevant bounds from other experiments.

Figure 4: Range of $g^{\prime}$ and $M_{Z^{\prime}}$ for which the global
analysis of oscillation data can be consistently described within the LMA-D
solution, Eq. (38), (black band) for two viable models. We also show the
bounds on these models imposed by a compilation of experiments as labeled in
the figure. In the window corresponding to the model in Ref. Farzan:2015doa
we have defined $B_{y}\equiv B_{1}-yB_{2}-(3-y)B_{3}$. See text for details.

As seen in the figure, for the model in Ref. Farzan:2015doa , there are
solutions for the $Z^{\prime}$ coupling and mass for which LMA-D is allowed
without conflict with bounds from other experiments both for very light
mediators as well as for an $\mathcal{O}(\text{MeV--GeV})$ $Z^{\prime}$. More
quantitatively for ultra light mediator part of the LMA-D allowed parameter
space for the model in Ref. Farzan:2015doa  is in conflict with the bounds
from fifth force tests which impose the stronger constraints for this model in
this regime. But the LMA-D is still a viable solution for

$10^{-13}~{}\text{eV}\leq M_{Z^{\prime}}\leq 2.5\times 10^{-9}~{}\text{eV}$
(39)

with couplings in a very narrow band and seen in the figure, for example

$g^{\prime}=(9.1\pm 0.2)\times 10^{-25}\quad\text{for}\quad
M_{Z^{\prime}}=10^{-13}~{}\text{eV}$ (40)

and

$g^{\prime}=(1.22\pm 0.04)\times 10^{-20}\quad\text{for}\quad
M_{Z^{\prime}}=2.5\times 10^{-9}~{}\text{eV}\,.$ (41)

In addition LMA-D is also a viable solution for this model with

$5~{}\text{MeV}\leq M_{Z^{\prime}}\leq
20~{}\text{MeV}\quad\text{with}\quad\frac{g^{\prime}}{M_{Z^{\prime}}}=\frac{(4.85\pm
0.15)\times 10^{-5}}{\text{MeV}}$ (42)

On the contrary, as seen in the figure, the model with charge
$L_{e}+2L_{\mu}+2L_{\tau}$ can only provide a viable LMA-D solution without
conflict with bounds from other experiments for and ultra light mediator with
$M_{Z^{\prime}}\lesssim 5\times 10^{-12}$ eV.

##  4 Conclusions

In this work we have performed dedicated global analysis of oscillation data
in the framework of lepton flavor dependent $U(1)^{\prime}$ interactions which
affect the neutrino evolution in matter. The analysis is performed for
interaction lengths ranging from larger than the Sun radius (covering what we
label as LRI domain) to effective contact neutrino interactions (NSI domain).
We survey $\sim 10000$ set of models characterized by the charges of the first
generation charged fermions and the three flavor neutrinos. We find that

  * •

In the LIGHT sector of the oscillation parameter space the introduction of new
interactions does not lead to a significantly better fit of the oscillation
data compared to standard oscillations, irrespective of the $U(1)^{\prime}$
coupling in either NSI or LRI domains. Thus for all cases the analysis of
oscillation data in the LIGHT sector results in excluded ranges of
$g^{\prime}$ and $M_{Z^{\prime}}$.

  * •

The excluded ranges for the $Z^{\prime}$ coupling and mass for the models with
couplings listed in table 1 are shown in Figs. 2 and 3 for $M_{Z^{\prime}}$
below 1 eV and in the $\mathcal{O}(\text{MeV--GeV})$ range, respectively –
i.e., below and above the window strongly disfavored by the cosmological bound
on $\Delta N_{\text{eff}}$.

  * •

In the regime of ultra-light mediators, for all models shown, the oscillation
analysis improves over the bounds from tests of fifth forces and of violation
of equivalence principle and for $Z^{\prime}$ lighter than $\sim 10^{-8}$ eV
or $\sim 10^{-11}$ eV depending on whether the $U(1)^{\prime}$ current
involves coupling to electron lepton number.

  * •

For mediators in the $\mathcal{O}(\text{MeV--GeV})$ range we list in table 3
the derived constraints on $g^{\prime}$ versus $M_{Z^{\prime}}$. We find that
for several of the considered models there are values of $g^{\prime}$ and
$M_{Z^{\prime}}$ for which the oscillation analysis provides constraints
extending beyond those from other experiments.

  * •

In what respects to LMA-D we find that it cannot be realized in the LRI
domain. In the NSI domain we have found that 4.8% of the set of charges
studied can have a best fit in LMA-D, and 5.2% lead to LMA-D as a valid
solution within 4 units of $\Delta\chi^{2}$ of standard oscillations. None of
these set of charges correspond to an anomaly free model based on gauging SM
global symmetries with SM plus right-handed neutrinos matter content (Eq.
(3)). So, generically, $Z^{\prime}$ models for LMA-D with gauged SM global
symmetries require additional states for anomaly cancellation.

### Acknowledgement

We want to thank Renata Zukanovich for discussions and her participation in
the early stages of this work. The authors also acknowledge use of the HPC
facilities at the IFT (Hydra cluster) This work was supported by the Spanish
grants FPA2016-76005-C2-1-P, FPA2016-78645-P, PID2019-105614GB-C21 and
PID2019-110058GB-C21, by EU Network FP10 ITN HIDDEN (H2020-MSCA-
ITN-2019-860881), by USA-NSF grant PHY-1915093, by AGAUR (Generalitat de
Catalunya) grant 2017-SGR-929. The authors acknowledge the support of the
Spanish Agencia Estatal de Investigación through the grant “IFT Centro de
Excelencia Severo Ochoa SEV-2016-0597”. PC acknowledges funding from the
Spanish MICINN through the “Ramón y Cajal” program under grant
RYC2018-024240-I.

##  Appendix A Bounds from non-oscillation experiments

Bounds on dark photons or light $Z^{\prime}$ bosons can be set using a variety
of laboratory experiments, collider searches and cosmological and
astrophysical probes. While in this work we focus on the most relevant ones
for the range of masses under consideration, we refer the interested reader to
Refs. Ilten:2018crw ; Fabbrichesi:2020wbt ; Bauer:2018onh ; Wise:2018rnb ;
Jaeckel:2010ni ; Essig:2013lka ; Harnik:2012ni  for a selection of works which
include a systematic compilation of bounds. In this appendix we summarize the
main details of the bounds shown by the colored regions in Figs. 2 and 3, as
well as the procedure used to rescale them for the different models shown in
each panel.

### Bounds from searches for gravitational fifth forces

Generically gravitational fifth force experiments look for deviation from the
standard Newton potential ($\propto 1/r$) between two objects.

We have taken the results from Ref. Salumbides:2013dua  (which, in turn, were
recasted from Fig. 10 in Ref. Adelberger:2009zz ) where they present the
constraints on the coupling $\alpha_{5}$ versus the interaction length
$\lambda$ (or, equivalently, versus $m_{5}=\frac{1}{\lambda}$) defined as the
constant entering the potential of the fifth force

$V_{5}(r)=\alpha_{5}\,\mathcal{N}_{1}\,\mathcal{N}_{2}\,\frac{e^{-r/\lambda}}{r}$
(43)

where $\mathcal{N}_{1,2}$ is the total charge of each object, which they take
as the total number of baryons.

In particular from Fig. 3 in Salumbides:2013dua  we read their boundary curve

$\left(\frac{\alpha^{\text{max}}_{5}}{\alpha_{em}}\,,\,m_{5}\right)\,,$ (44)

where $\alpha_{em}$ is the SM fine structure constant. In order to rescale it
to the different $U(1)^{\prime}$ models we notice that the corresponding
potential for the $U(1)^{\prime}$ interaction for the same objects is

$V^{\prime}(r)=C_{1}C_{2}\frac{{g^{\prime}}^{2}}{4\pi}\frac{e^{-r/\lambda^{\prime}}}{r}\,,$
(45)

where $C_{i}$ is the total $Z^{\prime}$ charge of object $i$ and
$\lambda^{\prime}$ is the $U(1)^{\prime}$ interaction length

$\lambda^{\prime}=\frac{1}{M_{Z^{\prime}}}\,,\quad
C_{i}\equiv\frac{\mathcal{N}_{i}}{A_{i}}c_{i}\equiv\frac{\mathcal{N}_{i}}{A_{i}}\left[Z_{i}(a_{e}+2a_{u}+a_{d})+(A_{i}-Z_{i})(a_{u}+2a_{d})\right]\,,$
(46)

where $Z_{i}$ and $A_{i}$ are the atomic number and mass number of the
material of which the object $i$ is made, so $c_{i}$ is the charge under
$Z^{\prime}$ for each atom of the material. Assuming that the number of
protons and neutrons in the material are not very different (that is,
$A_{i}/Z_{i}\sim 2$), we get

$V^{\prime}(r)=V_{5}(r)\times\frac{g^{\prime
2}}{4\pi\alpha_{5}}\frac{(a_{e}+3a_{u}+3a_{d})^{2}}{4}$ (47)

with $m_{5}=M_{Z^{\prime}}$. So the boundary in the $g^{\prime}$ vs
$M_{Z^{\prime}}$ plane will be

$\big{(}g^{\prime}_{\text{max}},M_{Z^{\prime}}\big{)}=\left(\frac{2}{a_{e}+3a_{u}+3a_{d}}\sqrt{4\pi\alpha_{em}}\times\sqrt{\frac{\alpha^{\text{max}}_{5}}{\alpha_{em}}}\,,\,m_{5}\right)\,.$
(48)

### Bounds from searches for violation of the equivalence principle

Similarly as in the case of fifth force searches, this bound comes from
precise measurements of the gravitational potential between two objects.
However, in this case one tests the differences of the potential for the same
total mass of two different test materials, using a pendulum which is
attracted by the same mass. The bounds are taken from Ref. Schlamminger:2007ht
, where they define the potential due to the new force as:

$V_{G}(r)=\alpha\,G\frac{m_{t}m_{s}}{r}\hat{\mathcal{N}}_{t}\hat{\mathcal{N}}_{s}e^{-r/\lambda}\,.$
(49)

where $G$ is the gravitational constant, $\hat{\mathcal{N}}_{i}$ refers to the
new interaction charge per mass unit, and the subindices $t$ and $s$ stand for
test (or pendulum) and source masses, respectively. In Ref.
Schlamminger:2007ht  they assume that the new interaction couples to the
number of baryons. Thus, they use beryllium and titanium as test materials,
chosen to maximize the difference in baryon number per unit mass. To be
specific, we use $\hat{\mathcal{N}}_{\text{Be}}=0.99868$ and
$\hat{\mathcal{N}}_{\text{Ti}}=1.001077$, as in Ref. Schlamminger:2007ht .

Noting that the total potential will be the sum of the standard gravitational
potential plus the contribution from the new interaction, the authors of Ref.
Schlamminger:2007ht  put a bound on

$\eta=2\frac{V_{\text{Be}}-V_{\text{Ti}}}{V_{\text{Be}}+V_{\text{Ti}}}\sim\frac{\Delta
V_{G}(r)}{\frac{Gm_{s}m_{t}}{r}}=\alpha(\hat{\mathcal{N}}_{\text{Ti}}-\hat{\mathcal{N}}_{\text{Be}})\hat{\mathcal{N}_{s}}e^{-r/\lambda}$
(50)

which is used to set a constraint on $\alpha$ as a function of $\lambda$,

$\big{(}\alpha^{\text{max}}\,,\,\lambda\big{)}\,.$ (51)

This can be used to set a bound on $\Delta V^{\prime}(r)$ for a general model,
noting that

$\begin{split}\Delta
V^{\prime}(r)&=\frac{{g^{\prime}}^{2}}{4\pi}\frac{e^{-rM_{Z^{\prime}}}}{r}\,C_{s}(C_{\text{Be}}-C_{\text{Ti}})\\\
&=\Delta V_{G}(r)\times\frac{g^{\prime 2}/4\pi}{\alpha G\cdot
u^{2}}\times\frac{c_{\text{Ti}}/A_{\text{Ti}}-c_{\text{Be}}/A_{\text{Be}}}{\hat{\mathcal{N}}_{\text{Ti}}-\hat{\mathcal{N}}_{\text{Be}}}\times\frac{c_{s}/A_{s}}{\hat{\mathcal{N}}_{s}}\\\
&\simeq\Delta V_{G}(r)\times\frac{g^{\prime 2}/4\pi}{\alpha G\cdot
u^{2}}\times\frac{(a_{e}+a_{u}-a_{d})(a_{e}+3a_{u}+3a_{d})}{144(\hat{\mathcal{N}}_{\text{Ti}}-\hat{\mathcal{N}}_{\text{Be}})}\end{split}$
(52)

where $u$ stands for the atomic mass unit in GeV ($u\simeq 0.931$ GeV), and in
the last line we have approximated for the source material $A/Z\simeq 2$ and
$\hat{\mathcal{N}}_{s}\simeq 1$. So the boundary in the $g^{\prime}$ vs
$M_{Z^{\prime}}$ plane will be

$\big{(}g^{\prime}_{\text{max}}\,,\,M_{Z^{\prime}}\big{)}=\left(\sqrt{\frac{144(\hat{\mathcal{N}}_{\text{Ti}}-\hat{\mathcal{N}}_{\text{Be}})}{(a_{e}+3a_{u}+3a_{d})(a_{e}+a_{u}-a_{d})}\,4\pi\,G\,u^{2}}\times\sqrt{\alpha^{\text{max}}}\,,\,\frac{1}{\lambda}\right)$
(53)

### Bounds from white-dwarf cooling

We based the bounds shown in Figs. 3 and 4 on the study presented in Ref.
Dreiner:2013tja . There, upper bounds on new interactions are imposed on the
basis that the energy losses from plasmon decays into particles that escape
the star is not larger than the energy losses due to neutrino emission in the
SM.

In the $U(1)^{\prime}$ scenarios here considered the minimum new contribution
is due to $Z^{\prime}$ mediated decays into neutrinos

$\Gamma_{\text{plasmon}\to\nu\bar{\nu},Z^{\prime}}^{s}\lesssim\Gamma_{\text{plasmon}\to\nu\bar{\nu},\text{SM}}^{s}=\frac{C_{e,V}^{2}G_{F}^{2}}{48\pi^{2}\alpha_{em}}\frac{Z_{s}\pi_{s}^{3}}{\omega_{s}}$
(54)

where $C_{e,V}$ is the vector coupling to the electron current in the SM,
$Z_{s}$ is the plasmon wavefunction renormalization and $\pi_{s}$ is the
effective plasmon mass which enters in the dispersion relation
$\omega^{2}_{s}-k^{2}=\pi_{s}(\omega_{s},k)$. Here, $s=T,L$ refers to the
plasmon polarization (transverse or longitudinal).

Under the assumption that the mass of the $Z^{\prime}$ is much larger than the
frequency of the plasmon we can write its rate into neutrinos of a given
flavor $\beta$ due to the new interactions as

$\Gamma_{\text{plasmon}\to\nu_{\beta}\bar{\nu}_{\beta},Z^{\prime}}^{s}=\frac{1}{3}\frac{{g^{\prime}}^{4}}{M_{Z^{\prime}}^{4}}\frac{(a_{e}\,b_{\beta})^{2}}{48\pi^{2}\alpha_{em}}\frac{Z_{s}\pi_{s}^{3}}{\omega_{s}}$
(55)

And the upper bound obtained in Ref. Dreiner:2013tja  translates into:

$\sqrt{\sum_{\beta}\frac{(a_{e}\,b_{\beta})^{2}}{3}}\cdot\frac{{g^{\prime}}^{2}}{M_{Z^{\prime}}^{2}}\leq
C_{e,V}\,G_{F}=1.12\times 10^{-5}~{}\text{GeV}^{-2}\,.$ (56)

### Bounds from coherent neutrino-nucleus scattering

We have performed our own reanalysis of coherent neutrino-nucleus scattering
data using the time and energy information from COHERENT experiment
Akimov:2017ade ; Akimov:2018vzs  on CsI based on our recent analysis in Ref.
Coloma:2019mbs  performed for NSI with a variety of nuclear form factors,
quenching factors and parametrization of the background. In particular the
results shown in Fig. 3 correspond to the analysis performed using the
quenching factor obtained with the fit to the calibration data of the Duke
(TUNL) group Akimov:2017ade  together with our data driven reevaluation of the
steady-state background (see Ref. Coloma:2019mbs  for details).

Model predictions are obtained exactly as in Coloma:2019mbs , but replacing
the cross section for coherent scattering in the presence of NSI with that
induced by the $Z^{\prime}$. In particular if we define

$\varepsilon^{\prime}(Q^{2},g^{\prime},M_{Z^{\prime}})\equiv\frac{1}{\sqrt{2}G_{F}}\frac{{g^{\prime}}^{2}}{M_{Z^{\prime}}^{2}+Q^{2}}$
(57)

one can use the same cross section expression for NSI simply replacing:

$\varepsilon_{\alpha\beta}^{q,V}\to\delta_{\alpha\beta}\,a_{q}\,b_{\alpha}\,\varepsilon^{\prime}(Q^{2},g_{X},M_{Z^{\prime}})$
(58)

where $Q^{2}=2MT$ is the momentum transfer. In this case, the differential
cross section for coherent scattering of a neutrino with flavor $\alpha$
reads:

$\frac{d\sigma_{\alpha}}{dT}=\frac{G_{F}^{2}}{2\pi}W_{\alpha}^{2}(Q^{2},g_{X},M_{Z^{\prime}})F^{2}(Q^{2})M\left(2-\frac{MT}{E_{\nu}^{2}}\right)$
(59)

where $M$ is the mass of the nucleus, $T$ is the nuclear recoil energy, and
$E_{\nu}$ is the incident neutrino energy. We have defined a modified weak
charge for the nucleus as:

$\begin{split}W_{\alpha}(Q^{2},g_{X},M_{Z^{\prime}})&=Z\big{[}g_{p}^{V}+(2a_{u}+a_{d})\,b_{\alpha}\,\varepsilon^{\prime}(Q^{2},g_{X},M_{Z^{\prime}})\big{]}+\\\
&+N\big{[}g_{n}^{V}+(a_{u}+2a_{d})\,b_{\alpha}\,\varepsilon^{\prime}(Q^{2},g_{X},M_{Z^{\prime}})\big{]}\,.\end{split}$
(60)

Thus, unlike in the case of NSI, the weak charge now depends on the momentum
transferred in the process.

The construction of $\chi^{2}_{\text{COH}}(g^{\prime},M_{Z^{\prime}})$ is
totally analogous to that of $\chi^{2}_{\text{COH}}(\vec{\varepsilon})$ in
Ref. Coloma:2019mbs . In consistency with the condition impose when deriving
the oscillation bounds the COHERENT regions shown in Fig. 3 corresponds to the
$Z^{\prime}$ coupling and mass for which the fit to COHERENT data is worse
than the one obtained in the SM by 4 units:

$\chi^{2}_{\text{COH}}(g^{\prime},M_{Z^{\prime}})-\chi^{2}_{\text{COH,SM}}>4.$
(61)

Although it is difficult to compare to other works in the literature due to
the very different implementation of systematics, quenching factor, form
factor and background treatment, we find that our results are roughly
consistent with those presented in Refs. Cadeddu:2020nbr ; Papoulias:2019txv ;
Denton:2018xmq .

### Bounds from measurements of neutrino scattering on electrons

For neutrino-electron scattering experiments TEXONO and GEMMA we performed our
own analysis following the procedure in Ref. Lindner:2018kjo  which explicitly
studied the bounds imposed by those experiments in some $Z^{\prime}$ models
and provide all the relevant cross section expressions.

For TEXONO we use the data from Fig. 16 in Ref. Deniz:2009mu . This
corresponds to 29882 (7369) kg-day of fiducial mass exposure during Reactor ON
(OFF), respectively. The adopted analysis window is 3–8 MeV, spread out
uniformly over $N_{\text{bin}}=10$ energy bins. An overall normalization
constant has been manually set to reproduce the SM prediction shown in Fig 16
in Deniz:2009mu . Cross section has been implemented following Lindner:2018kjo
. With this we get $\chi^{2}_{\text{SM,TEXONO}}=15.87$.

For GEMMA we take data from Fig. 8 in Ref. Beda:2010hk  and consider both
phase I and phase II data so we define
$\chi^{2}_{\text{GEMMA,TOT}}=\chi^{2}_{\text{GEMMA,I}}+\chi^{2}_{\text{GEMMA,II}}$
This corresponds to about 12000 ON-hours and 3000 OFF-hours of active time.
Data was taken for a detector mass of 1.5 kg. The energy window used is
$0.015~{}\text{MeV}\leq E_{\nu}\leq 8.0~{}\text{MeV}$. With this procedure we
get $\chi^{2}_{\text{SM,GEMMA}}=154.97$.

For both TEXONO and GEMMA we draw the contours with the equivalent condition
used for the oscillation analysis so for a given set of model charges the
contours are defined by
$\chi^{2}(g^{\prime},M_{Z^{\prime}})-\chi^{2}_{\text{SM}}=4$.

Additional bounds from neutrino scattering on electrons can be derived from
the analysis of LSND or CHARM II data (see, e.g., Ref. Bilmis:2015lja ).
However, they are less constraining than those obtained for TEXONO and GEMMA
in the shown range of masses. They can be stronger for masses above $\sim 50$
MeV but in such case they are weaker than bounds from other experiments, in
particular in $e^{+}e^{-}$ collisions in BaBaR Lees:2014xha  (see below). For
this reason they are not shown in the figures.

### Fixed target experiments and colliders

As mentioned in the text, the bounds from fixed target experiments and
colliders are taken directly from the literature. In particular, we consider
the following set of bounds:

  * •

Electron beam dump experiments: we take these from the compilation in Ref.
Andreas:2012mt , which were obtained using data from E137 Bjorken:1988as ,
E141 Riordan:1987aw , E774 Bross:1989mp , Orsay Davier:1989wz , and KEK-PF-000
Konaka:1986cb .

  * •

Proton beam dump experiments: for LSND we use the results obtained in Ref.
Essig:2010gu  which used data from Ref. Auerbach:2003fz  (see also Ref.
Batell:2009di ); for CHARM Bergsma:1985qz  we use the limit derived in Ref.
Gninenko:2012eq ; finally, for $\nu$-Cal Blumlein:1990ay  we use the lmit
derived in Ref. Tsai:2019mtm  (see also Ref. Blumlein:2011mv  for a similar
analysis).

  * •

We also consider constraints from $Z^{\prime}$ production in $e^{+}e^{-}$
collisions in BaBaR both in visible Lees:2014xha  and invisible Lees:2017lec
final states, as well as constraints from LHCb for $U(1)^{\prime}$ decaying
into $\mu^{+}\mu^{-}$, the most relevant ones being from prompt decay searches
Aaij:2019bvg ; Ilten:2016tkc .

  * •

Constraints on invisible $Z^{\prime}$ decays in searches for
$\pi^{0}\to\gamma\,Z^{\prime}$ at NA62 experiment CortinaGil:2019nuo  and in
$e\,N\to e\,N\,Z^{\prime}$ at NA64 experiment Banerjee:2016tad ;
Banerjee:2017hhz ; NA64:2019imj . There are bounds from NA64 for the
$Z^{\prime}$ decaying into $e^{+}e^{-}$ Banerjee:2019hmi , but for the models
in Fig. 3 they are weaker than those from other experiments.

All the bounds mentioned above were obtained for a dark photon coupled to the
SM fermions via kinetic mixing. In order to recast these to bounds on the
different $U(1)^{\prime}$ models considered we have used the darkcast software
Ilten:2018crw , which takes into account the difference in production
branching ratio and lifetime of the new boson, as well as its decay into a
given final state. In particular for invisible decays we include only the
decays into neutrinos.

In addition bounds on $U(1)^{\prime}$ models coupled to charges including
$L_{\mu}$ can be constrained with data on production of $\mu^{+}\mu^{-}$ in
$\nu_{\mu}$ scattering off the Coulomb field of a nucleus (the so called
neutrino trident production) at CHARM-II Geiregat:1990gz  and the CCFR
experiment Mishra:1991bv , see Ref. Altmannshofer:2014pba . They can also be
bounded with data on $\mu^{+}\mu^{-}$ production in $e^{+}e^{-}$ collision at
BaBar TheBABAR:2016rlg  and Belle Adachi:2019otg . We have verified that all
these bounds are always weaker than those imposed by either coherent neutrino-
nucleus scattering, or LHCb for the same models. For this reason they are not
shown in the figures.

### Astrophysical and cosmological bounds for
$M_{Z^{\prime}}\gtrsim\mathcal{O}(\text{MeV})$

The impact of light $Z^{\prime}$ on astrophysical and cosmological observables
can also be used to set strong constraints on these models.

For example, non-standard cooling mechanisms in the Sun, other stars or
supernovae from $Z^{\prime}$ emission can be used to set very tight
constraints on these models. Solar and stellar constraints are only relevant
in the mass window $\textrm{eV}\lesssim M_{Z^{\prime}}\lesssim
100~{}\textrm{eV}$ (see, e.g., Refs. Harnik:2012ni ; Wise:2018rnb ;
Jaeckel:2010ni ; Davidson:2000hf ) and therefore will not be considered here.
More relevant are the SN constraints, which also apply to masses above a few
MeV and therefore would be relevant in the high mass window considered in this
work. However, while these have been derived in the dark photon scenario, the
production mechanisms would be significantly affected in the $Z^{\prime}$
case. As an example, in Ref. Croon:2020lrf  a specific analysis carried out
for supernovae emission for the $L_{\mu}-L_{\tau}$ and $B-L$ models showed
large variations in the results with respect to the region constrained in the
dark photon case Raffelt:2006cw ; Chang:2016ntp . A dedicated analysis would
be required to adapt these bounds to models with arbitrary charges.

Finally additional constraints can be derived from cosmological observations,
from the energy injection of the $Z^{\prime}$ onto $e^{+}e^{-}$ in the early
Universe, which would be applicable in the mass region above 1 MeV (see, e.g.,
Fig. 6 in Coffey:2020oir ). However, again in this case it is uncertain how to
recast these bounds to a general $U(1)^{\prime}$ model with arbitrary charges,
and therefore we have not included them here.

## References

  * (1) B. Pontecorvo, _Neutrino Experiments and the Problem of Conservation of Leptonic Charge_ , _Sov. Phys. JETP_ 26 (1968) 984. [Zh. Eksp. Teor. Fiz.53,1717(1967)].
  * (2) V.N. Gribov and B. Pontecorvo, _Neutrino astronomy and lepton charge_ , _Phys. Lett._ B28 (1969) 493.
  * (3) M.C. Gonzalez-Garcia and M. Maltoni, _Phenomenology with Massive Neutrinos_ , _Phys. Rept._ 460 (2008) 1 [0704.1800].
  * (4) L. Wolfenstein, _Neutrino Oscillations in Matter_ , _Phys. Rev._ D17 (1978) 2369.
  * (5) S.P. Mikheev and A.Y. Smirnov, _Resonance enhancement of oscillations in matter and solar neutrino spectroscopy_ , _Sov. J. Nucl. Phys._ 42 (1985) 913.
  * (6) J.W.F. Valle, _Resonant Oscillations of Massless Neutrinos in Matter_ , _Phys. Lett._ B199 (1987) 432.
  * (7) M.M. Guzzo, A. Masiero and S.T. Petcov, _On the MSW effect with massless neutrinos and no mixing in the vacuum_ , _Phys. Lett._ B260 (1991) 154.
  * (8) T. Ohlsson, _Status of non-standard neutrino interactions_ , _Rept. Prog. Phys._ 76 (2013) 044201 [1209.2710].
  * (9) O.G. Miranda and H. Nunokawa, _Non standard neutrino interactions: current status and future prospects_ , _New J. Phys._ 17 (2015) 095002 [1505.06254].
  * (10) Y. Farzan and M. Tortola, _Neutrino oscillations and Non-Standard Interactions_ , _Front. in Phys._ 6 (2018) 10 [1710.09360].
  * (11) P.S. Bhupal Dev et al., _Neutrino Non-Standard Interactions: A Status Report_ , in _NTN Workshop on Neutrino Non-Standard Interactions St Louis, MO, USA, May 29-31, 2019_ , 2019, http://lss.fnal.gov/archive/2019/conf/fermilab-conf-19-299-t.pdf [1907.00991].
  * (12) M.C. Gonzalez-Garcia, P.C. de Holanda, E. Masso and R. Zukanovich Funchal, _Probing long-range leptonic forces with solar and reactor neutrinos_ , _JCAP_ 0701 (2007) 005 [hep-ph/0609094].
  * (13) J. Grifols and E. Masso, _Neutrino oscillations in the sun probe long range leptonic forces_ , _Phys. Lett. B_ 579 (2004) 123 [hep-ph/0311141].
  * (14) A.S. Joshipura and S. Mohanty, _Constraints on flavor dependent long range forces from atmospheric neutrino observations at super-Kamiokande_ , _Phys. Lett. B_ 584 (2004) 103 [hep-ph/0310210].
  * (15) H. Davoudiasl, H.-S. Lee and W.J. Marciano, _Long-Range Lepton Flavor Interactions and Neutrino Oscillations_ , _Phys. Rev. D_ 84 (2011) 013009 [1102.5352].
  * (16) M.B. Wise and Y. Zhang, _Lepton Flavorful Fifth Force and Depth-dependent Neutrino Matter Interactions_ , _JHEP_ 06 (2018) 053 [1803.00591].
  * (17) M.C. Gonzalez-Garcia, M. Maltoni and J. Salvado, _Testing matter effects in propagation of atmospheric and long-baseline neutrinos_ , _JHEP_ 05 (2011) 075 [1103.4365].
  * (18) M.C. Gonzalez-Garcia and M. Maltoni, _Determination of matter potential from global analysis of neutrino oscillation data_ , _JHEP_ 09 (2013) 152 [1307.3092].
  * (19) I. Esteban, M.C. Gonzalez-Garcia, M. Maltoni, I. Martinez-Soler and J. Salvado, _Updated Constraints on Non-Standard Interactions from Global Analysis of Oscillation Data_ , _JHEP_ 08 (2018) 180 [1805.04530].
  * (20) J. Heeck, M. Lindner, W. Rodejohann and S. Vogl, _Non-Standard Neutrino Interactions and Neutral Gauge Bosons_ , _SciPost Phys._ 6 (2019) 038 [1812.04067].
  * (21) O.G. Miranda, M.A. Tortola and J.W.F. Valle, _Are solar neutrino oscillations robust?_ , _JHEP_ 10 (2006) 008 [hep-ph/0406280].
  * (22) Y. Farzan, _A model for large non-standard interactions of neutrinos leading to the LMA-Dark solution_ , _Phys. Lett._ B748 (2015) 311 [1505.06906].
  * (23) Y. Farzan and I.M. Shoemaker, _Lepton Flavor Violating Non-Standard Interactions via Light Mediators_ , _JHEP_ 07 (2016) 033 [1512.09147].
  * (24) K.S. Babu, A. Friedland, P.A.N. Machado and I. Mocioiu, _Flavor Gauge Models Below the Fermi Scale_ , 1705.01822.
  * (25) P.B. Denton, Y. Farzan and I.M. Shoemaker, _Testing large non-standard neutrino interactions with arbitrary mediator mass after COHERENT data_ , _JHEP_ 07 (2018) 037 [1804.03660].
  * (26) B.C. Allanach, J. Davighi and S. Melville, _An Anomaly-free Atlas: charting the space of flavour-dependent gauged $U(1)$ extensions of the Standard Model_, _JHEP_ 02 (2019) 082 [1812.04602]. [Erratum: JHEP08,064(2019)].
  * (27) J.A. Dror, R. Lasenby and M. Pospelov, _New constraints on light vectors coupled to anomalous currents_ , _Phys. Rev. Lett._ 119 (2017) 141803 [1705.06726].
  * (28) J.A. Dror, R. Lasenby and M. Pospelov, _Dark forces coupled to nonconserved currents_ , _Phys. Rev. D_ 96 (2017) 075036 [1707.01503].
  * (29) J.A. Dror, R. Lasenby and M. Pospelov, _Light vectors coupled to bosonic currents_ , _Phys. Rev. D_ 99 (2019) 055016 [1811.00595].
  * (30) Z. Maki, M. Nakagawa and S. Sakata, _Remarks on the unified model of elementary particles_ , _Prog. Theor. Phys._ 28 (1962) 870.
  * (31) M. Kobayashi and T. Maskawa, _CP Violation in the Renormalizable Theory of Weak Interaction_ , _Prog. Theor. Phys._ 49 (1973) 652.
  * (32) P. Coloma and T. Schwetz, _Generalized mass ordering degeneracy in neutrino oscillation experiments_ , _Phys. Rev._ D94 (2016) 055005 [1604.05772].
  * (33) A. Dziewonski and D. Anderson, _Preliminary reference earth model_ , _Phys.Earth Planet.Interiors_ 25 (1981) 297.
  * (34) T.-K. Kuo and J.T. Pantaleone, _THE SOLAR NEUTRINO PROBLEM AND THREE NEUTRINO OSCILLATIONS_ , _Phys.Rev.Lett._ 57 (1986) 1805.
  * (35) M. Guzzo, H. Nunokawa, P. de Holanda and O. Peres, _On the massless ’just-so’ solution to the solar neutrino problem_ , _Phys.Rev._ D64 (2001) 097301 [hep-ph/0012089].
  * (36) P. Bakhti and Y. Farzan, _Shedding Light on Lma-Dark Solar Neutrino Solution by Medium Baseline Reactor Experiments: Juno and Reno-50_ , _JHEP_ 07 (2014) 064 [1403.0744].
  * (37) I. Esteban, M. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, _Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of $\theta_{23}$, $\delta_{CP}$, and the mass ordering_, _JHEP_ 01 (2019) 106 [1811.05487].
  * (38) Y. Nakajima, “SuperKamiokande.” Talk given at the XXIX International Conference on Neutrino Physics and Astrophysics, Chicago, USA, June 22–July 2, 2020 (online conference) doi.org/10.5281/zenodo.3959640.
  * (39) I. Esteban, M. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, _The fate of hints: updated global analysis of three-flavor neutrino oscillations_ , 2007.14792.
  * (40) J. Yoo, “RENO.” Talk given at the XXIX International Conference on Neutrino Physics and Astrophysics, Chicago, USA, June 22–July 2, 2020 (online conference) doi.org/10.5281/zenodo.3959698.
  * (41) T. Bezerra, “New Results from the Double Chooz Experiment.” Talk given at the XXIX International Conference on Neutrino Physics and Astrophysics, Chicago, USA, June 22–July 2, 2020 (online conference) doi.org/10.5281/zenodo.3959542.
  * (42) P. Dunne, “Latest Neutrino Oscillation Results from T2K.” Talk given at the XXIX International Conference on Neutrino Physics and Astrophysics, Chicago, USA, June 22–July 2, 2020 (online conference) doi.org/10.5281/zenodo.3959558.
  * (43) A. Himmel, “New Oscillation Results from the NOvA Experiment.” Talk given at the XXIX International Conference on Neutrino Physics and Astrophysics, Chicago, USA, June 22–July 2, 2020 (online conference) doi.org/10.5281/zenodo.3959581.
  * (44) A. Kamada and H.-B. Yu, _Coherent Propagation of PeV Neutrinos and the Dip in theNeutrino Spectrum at IceCube_ , _Phys. Rev._ D92 (2015) 113004 [1504.00711].
  * (45) M. Escudero, D. Hooper, G. Krnjaic and M. Pierre, _Cosmology with a Very Light L μ-Lτ Gauge Boson_, _JHEP_ 03 (2019) 071 [1901.02010].
  * (46) M. Escudero and M. Fairbairn, _Cosmological Constraints on Invisible Neutrino Decays Revisited_ , _Phys. Rev. D_ 100 (2019) 103531 [1907.05425].
  * (47) M. Baryakhtar, R. Lasenby and M. Teo, _Black Hole Superradiance Signatures of Ultralight Vectors_ , _Phys. Rev. D_ 96 (2017) 035019 [1704.05081].
  * (48) M. Bustamante and S.K. Agarwalla, _Universe’s Worth of Electrons to Probe Long-Range Interactions of High-Energy Astrophysical Neutrinos_ , _Phys. Rev. Lett._ 122 (2019) 061103 [1808.02042].
  * (49) J.A. Dror, _Discovering leptonic forces using nonconserved currents_ , _Phys. Rev. D_ 101 (2020) 095013 [2004.04750].
  * (50) E.J. Salumbides, W. Ubachs and V.I. Korobov, _Bounds on fifth forces at the sub-Angstrom length scale_ , _J. Molec. Spectrosc._ 300 (2014) 65 [1308.1711].
  * (51) E. Adelberger, J. Gundlach, B. Heckel, S. Hoedl and S. Schlamminger, _Torsion balance experiments: A low-energy frontier of particle physics_ , _Prog. Part. Nucl. Phys._ 62 (2009) 102.
  * (52) S. Schlamminger, K.Y. Choi, T.A. Wagner, J.H. Gundlach and E.G. Adelberger, _Test of the equivalence principle using a rotating torsion balance_ , _Phys. Rev. Lett._ 100 (2008) 041101 [0712.0607].
  * (53) F. Kling, _Probing light gauge bosons in tau neutrino experiments_ , _Phys. Rev. D_ 102 (2020) 015007 [2005.03594].
  * (54) DONuT collaboration, _Final tau-neutrino results from the DONuT experiment_ , _Phys. Rev. D_ 78 (2008) 052002 [0711.0728].
  * (55) Y. Farzan and J. Heeck, _Neutrinophilic nonstandard interactions_ , _Phys. Rev._ D94 (2016) 053010 [1607.07616].
  * (56) P. Ilten, Y. Soreq, M. Williams and W. Xue, _Serendipity in dark photon searches_ , _JHEP_ 06 (2018) 004 [1801.04847].
  * (57) M. Fabbrichesi, E. Gabrielli and G. Lanfranchi, _The Dark Photon_ , 2005.01515.
  * (58) M. Bauer, P. Foldenauer and J. Jaeckel, _Hunting All the Hidden Photons_ , _JHEP_ 18 (2020) 094 [1803.05466].
  * (59) J. Jaeckel and A. Ringwald, _The Low-Energy Frontier of Particle Physics_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 405 [1002.0329].
  * (60) R. Essig et al., _Working Group Report: New Light Weakly Coupled Particles_ , in _Proceedings, 2013 Community Summer Study on the Future of U.S. Particle Physics: Snowmass on the Mississippi (CSS2013): Minneapolis, MN, USA, July 29-August 6, 2013_ , 2013, http://www.slac.stanford.edu/econf/C1307292/docs/IntensityFrontier/NewLight-17.pdf [1311.0029].
  * (61) R. Harnik, J. Kopp and P.A.N. Machado, _Exploring nu Signals in Dark Matter Detectors_ , _JCAP_ 1207 (2012) 026 [1202.6073].
  * (62) H.K. Dreiner, J.-F. Fortin, J. Isern and L. Ubaldi, _White Dwarfs constrain Dark Forces_ , _Phys. Rev._ D88 (2013) 043517 [1303.7232].
  * (63) COHERENT collaboration, _Observation of Coherent Elastic Neutrino-Nucleus Scattering_ , _Science_ 357 (2017) 1123 [1708.01294].
  * (64) COHERENT collaboration, _COHERENT Collaboration data release from the first observation of coherent elastic neutrino-nucleus scattering_ , 1804.09459.
  * (65) P. Coloma, I. Esteban, M. Gonzalez-Garcia and M. Maltoni, _Improved global fit to Non-Standard neutrino Interactions using COHERENT energy and timing data_ , _JHEP_ 02 (2020) 023 [1911.09109].
  * (66) M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. Li, E. Picciau et al., _Constraints on light vector mediators through coherent elastic neutrino nucleus scattering data from COHERENT_ , 2008.05022.
  * (67) D.K. Papoulias, _COHERENT constraints after the Chicago-3 quenching factor measurement_ , 1907.11644.
  * (68) M. Lindner, F.S. Queiroz, W. Rodejohann and X.-J. Xu, _Neutrino-electron scattering: general constraints on $Z^{\prime}$ and dark photon models_, _JHEP_ 05 (2018) 098 [1803.00060].
  * (69) TEXONO collaboration, _Measurement of Nu(e)-bar -Electron Scattering Cross-Section with a CsI(Tl) Scintillating Crystal Array at the Kuo-Sheng Nuclear Power Reactor_ , _Phys. Rev. D_ 81 (2010) 072001 [0911.1597].
  * (70) A. Beda, V. Brudanin, V. Egorov, D. Medvedev, V. Pogosov, M. Shirchenko et al., _Upper limit on the neutrino magnetic moment from three years of data from the GEMMA spectrometer_ , 1005.2736.
  * (71) S. Bilmis, I. Turan, T. Aliev, M. Deniz, L. Singh and H. Wong, _Constraints on Dark Photon from Neutrino-Electron Scattering Experiments_ , _Phys. Rev. D_ 92 (2015) 033009 [1502.07763].
  * (72) BaBar collaboration, _Search for a Dark Photon in $e^{+}e^{-}$ Collisions at BaBar_, _Phys. Rev. Lett._ 113 (2014) 201801 [1406.2980].
  * (73) S. Andreas, C. Niebuhr and A. Ringwald, _New Limits on Hidden Photons from Past Electron Beam Dumps_ , _Phys. Rev._ D86 (2012) 095019 [1209.6083].
  * (74) J.D. Bjorken, S. Ecklund, W.R. Nelson, A. Abashian, C. Church, B. Lu et al., _Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump_ , _Phys. Rev._ D38 (1988) 3375.
  * (75) E.M. Riordan et al., _A Search for Short Lived Axions in an Electron Beam Dump Experiment_ , _Phys. Rev. Lett._ 59 (1987) 755.
  * (76) A. Bross, M. Crisler, S.H. Pordes, J. Volk, S. Errede and J. Wrbanek, _A Search for Shortlived Particles Produced in an Electron Beam Dump_ , _Phys. Rev. Lett._ 67 (1991) 2942.
  * (77) M. Davier and H. Nguyen Ngoc, _An Unambiguous Search for a Light Higgs Boson_ , _Phys. Lett._ B229 (1989) 150.
  * (78) A. Konaka et al., _Search for Neutral Particles in Electron Beam Dump Experiment_ , _Phys. Rev. Lett._ 57 (1986) 659.
  * (79) R. Essig, R. Harnik, J. Kaplan and N. Toro, _Discovering New Light States at Neutrino Experiments_ , _Phys. Rev._ D82 (2010) 113008 [1008.0636].
  * (80) LSND collaboration, _Search for pi0 — $>$ nu(mu) anti-nu(mu) decay in LSND_, _Phys. Rev. Lett._ 92 (2004) 091801 [hep-ex/0310060].
  * (81) B. Batell, M. Pospelov and A. Ritz, _Exploring Portals to a Hidden Sector Through Fixed Targets_ , _Phys. Rev._ D80 (2009) 095024 [0906.5614].
  * (82) CHARM collaboration, _Search for Axion Like Particle Production in 400-GeV Proton - Copper Interactions_ , _Phys. Lett._ 157B (1985) 458.
  * (83) S.N. Gninenko, _Constraints on sub-GeV hidden sector gauge bosons from a search for heavy neutrino decays_ , _Phys. Lett._ B713 (2012) 244 [1204.3583].
  * (84) J. Blumlein et al., _Limits on neutral light scalar and pseudoscalar particles in a proton beam dump experiment_ , _Z. Phys. C_ 51 (1991) 341.
  * (85) Y.-D. Tsai, P. deNiverville and M.X. Liu, _The High-Energy Frontier of the Intensity Frontier: Closing the Dark Photon, Inelastic Dark Matter, and Muon g-2 Windows_ , 1908.07525.
  * (86) J. Blumlein and J. Brunner, _New Exclusion Limits for Dark Gauge Forces from Beam-Dump Data_ , _Phys. Lett._ B701 (2011) 155 [1104.2747].
  * (87) BaBar collaboration, _Search for Invisible Decays of a Dark Photon Produced in ${e}^{+}{e}^{-}$ Collisions at BaBar_, _Phys. Rev. Lett._ 119 (2017) 131804 [1702.03327].
  * (88) LHCb collaboration, _Search for $A^{\prime}\to\mu^{+}\mu^{-}$ Decays_, _Phys. Rev. Lett._ 124 (2020) 041801 [1910.06926].
  * (89) P. Ilten, Y. Soreq, J. Thaler, M. Williams and W. Xue, _Proposed Inclusive Dark Photon Search at LHCb_ , _Phys. Rev. Lett._ 116 (2016) 251803 [1603.08926].
  * (90) NA62 collaboration, _Search for production of an invisible dark photon in $\pi^{0}$ decays_, _JHEP_ 05 (2019) 182 [1903.08767].
  * (91) NA64 collaboration, _Search for invisible decays of sub-GeV dark photons in missing-energy events at the CERN SPS_ , _Phys. Rev. Lett._ 118 (2017) 011802 [1610.02988].
  * (92) NA64 collaboration, _Search for vector mediator of Dark Matter production in invisible decay mode_ , _Phys. Rev. D_ 97 (2018) 072002 [1710.00971].
  * (93) D. Banerjee et al., _Dark matter search in missing energy events with NA64_ , _Phys. Rev. Lett._ 123 (2019) 121801 [1906.00176].
  * (94) NA64 collaboration, _Improved limits on a hypothetical X(16.7) boson and a dark photon decaying into $e^{+}e^{-}$ pairs_, _Phys. Rev. D_ 101 (2020) 071101 [1912.11389].
  * (95) CHARM-II collaboration, _First observation of neutrino trident production_ , _Phys. Lett. B_ 245 (1990) 271.
  * (96) CCFR collaboration, _Neutrino tridents and W Z interference_ , _Phys. Rev. Lett._ 66 (1991) 3117.
  * (97) 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 [1406.2332].
  * (98) BaBar collaboration, _Search for a muonic dark force at BABAR_ , _Phys. Rev. D_ 94 (2016) 011102 [1606.03501].
  * (99) Belle II collaboration, _Search for an Invisibly Decaying $Z^{\prime}$ Boson at Belle II in $e^{+}e^{-}\to\mu^{+}\mu^{-}(e^{\pm}\mu^{\mp})$ Plus Missing Energy Final States_, _Phys. Rev. Lett._ 124 (2020) 141801 [1912.11276].
  * (100) S. Davidson, S. Hannestad and G. Raffelt, _Updated bounds on millicharged particles_ , _JHEP_ 05 (2000) 003 [hep-ph/0001179].
  * (101) D. Croon, G. Elor, R.K. Leane and S.D. McDermott, _Supernova Muons: New Constraints on Z’ Bosons, Axions, and ALPs_ , 2006.13942.
  * (102) G.G. Raffelt, _Astrophysical axion bounds_ , _Lect. Notes Phys._ 741 (2008) 51 [hep-ph/0611350].
  * (103) J.H. Chang, R. Essig and S.D. McDermott, _Revisiting Supernova 1987A Constraints on Dark Photons_ , _JHEP_ 01 (2017) 107 [1611.03864].
  * (104) J. Coffey, L. Forestell, D.E. Morrissey and G. White, _Cosmological Bounds on sub-GeV Dark Vector Bosons from Electromagnetic Energy Injection_ , 2003.02273.


mixing weight parameter estimate $0.77$, and it covers approximately $63\%$ of
the series (approximated as the mean of the estimated mixing weights), whereas
the second regime has the implied mixing weight parameter estimate $0.23$ and
it covers approximately $37\%$ of the series. Therefore, they are somewhat
disproportionate, but seem reasonable enough not to distort the generalized
impulse response functions too much.

|  |  | $IPI$ | $HCPI$ | $OIL$ | $RATE$
---|---|---|---|---|---|---
| $\hat{\alpha}_{m}$ | $\hat{\nu}_{m}$ | $\hat{\mu}_{m,1}$ | $\hat{\sigma}_{m,1}^{2}$ | $\hat{\mu}_{m,2}$ | $\hat{\sigma}_{m,2}^{2}$ | $\hat{\mu}_{m,3}$ | $\hat{\sigma}_{m,3}^{2}$ | $\hat{\mu}_{m,4}$ | $\hat{\sigma}_{m,4}^{2}$
Regime 1 | $0.77$ | $\phantom{1}3.40$ | $-2.89$ | $111.22$ | $0.12$ | $0.08$ | $0.02$ | $1.79$ | $0.24$ | $78.96$
Regime 2 | $0.23$ | $12.89$ | $\phantom{-}2.79$ | $\phantom{1}10.25$ | $0.19$ | $0.02$ | $0.33$ | $0.77$ | $2.71$ | $\phantom{7}0.47$
Table 2: Mixing weight parameter estimates ($\hat{\alpha}_{m}$), degrees of
freedom parameter estimates ($\hat{\nu}_{m}$), and marginal stationary means
($\hat{\mu}_{m,i}$) and variances ($\hat{\sigma}_{m,i}^{2}$) of the component
series implied by the fitted StMVAR($2,2$) model for each of the regimes.

The degrees of freedom parameter estimates in Table 2 show that the first
regime has fatter tailed distribution than the second one. The first regime
also has negative and volatile long-run output gap, while the second regime
has positive and less volatile long-run output gap. Long-run inflation, is low
in the first regime (roughly $1.4\%$ yearly), whereas it is moderate in the
second regime (roughly $2.3\%$ yearly). Also oil price inflation is relatively
low in first regime, whereas it is high in the second one.181818The log-
difference of oil price was multiplied by $10$ and not $100$ for numerical
reasons, so the unconditional means should be multiplied $10$ to obtain
estimates for the (approximate) monthly long-run oil price inflation in
percentage units. The interest rate variable has low mean in first regime, but
the variance is high, which reflects wandering movements of the shadow rate
(Wu and Xia, 2016) after the early 2010’s recession. In the second regime, the
interest rate variable has moderate mean and low variance. According to the
unconditional variances of the observable variables in Table 2, the first
regime appears overall more volatile than the second.

According to the estimated mixing weights presented in Figure 1, the first
regime (blue) mainly prevails after the collapse of Lehman Brothers in the
Financial crisis in September 2008. The first regime also obtains large mixing
weights during and before the early 2000’s recession, however. The second
regime (red) dominates when the first one does not; that is, mainly before the
Financial crisis, but excluding the aforementioned periods when the first
regime obtains large mixing weights. After the Financial crisis, the second
regime’s mixing weights stay close to zero, excluding a short period before
the early 2010’s recession. Because the first regime is characterized by
negative output gap and it mainly prevails after the Financial crisis, we
refer to it as the low-growth post-Financial crisis regime. Accordingly,
because the second regime is characterized by positive output gap and it
mainly prevails before the Financial crisis, we refer to it as the high-growth
pre-Financial crisis regime.

###  C.3 Estimates of the structural parameters

The constraints imposed on the B-matrix of the StMVAR($2,2$) model in Section
7.1 produced the following estimates:

$\hat{W}=\begin{bmatrix}\phantom{-}0.88\ (0.692)&-1.16\ (0.812)&2.00\
(1.164)&-0.40\ (0.351)\\\ -0.13\ (0.091)&\phantom{-}0.14\ (0.099)&0.20\
(0.122)&\phantom{-}0\phantom{.00\ (0.000)}\\\ -1.20\ (0.686)&-0.42\
(0.561)&0.65\ (0.435)&\phantom{-}0\phantom{.00\ (0.000)}\\\ \phantom{-}0.01\
(0.021)&-0.03\ (0.029)&0.06\ (0.045)&\phantom{-}0.56\ (0.313)\\\
\end{bmatrix},\ \ \hat{\lambda}_{2}=\begin{bmatrix}0.37\ (0.423)\\\ 0.27\
(0.303)\\\ 0.14\ (0.161)\\\ 0.02\ (0.027)\\\ \end{bmatrix},$ (C.1)

where the identified monetary policy shock is ordered last. The approximate
standard errors are presented in parentheses next to the estimates.