Neutrino oscillation constraints on U(1)′ models: from non-standard interactions to long-range forces

We quantify the effect of gauge bosons from a weakly coupled lepton flavor dependent U(1)′ 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 ∼ 10000 set of models characterized by the six relevant fermion U(1)′ charges and find that in all cases, constraints on the coupling and mass of the Z′ can be derived. We also find that about 5% of the U(1)′ 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 U1B−3Le\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{B-3{L}_e} $$\end{document}, U1B−3Lμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{B-3{L}_{\mu }} $$\end{document}, U1B−3Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{B-3{L}_{\tau }} $$\end{document}, U1B−32Lμ+Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{B-\frac{3}{2}\left({L}_{\mu }+{L}_{\tau}\right)} $$\end{document}, U1Le−Lμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{L_e-{L}_{\mu }} $$\end{document}, U1Le−Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{L_e-{L}_{\tau }} $$\end{document}, U1Le−12Lμ+Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{L_e-\frac{1}{2}\left({L}_{\mu }+{L}_{\tau}\right)} $$\end{document}. 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′ lighter than ∼ 10−8→ 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′ and MZ′ for which the oscillation analysis provides constraints beyond those imposed by laboratory experiments. Finally we illustrate the range of Z′ and couplings leading to a viable LMA-D solution for two sets of models.


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 [1,2], see Ref. [3] 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 [4,5].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 = √ 2G 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) [4,6,7] between neutrinos and matter (for recent reviews, see [8][9][10][11]), which can be parametrized as where G F is the Fermi constant, α, β are flavor indices, P ≡ 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, ε f,P αβ parametrizes the strength of the new interaction with respect to the Fermi constant, ε f,P αβ ∼ 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 ′ ≳ 10 −12 eV [12].
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 ∼ 1/M Z ′ around it [12][13][14][15][16].
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 [17][18][19].Analysis of early oscillation data was also used to impose constraints on flavor dependent long-range forces [12][13][14].
Straightforward constructions leading to Eq. (1.1) have an extended gauge sector with an additional U (1) ′ 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 ′ 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. [19] were adapted to place constraints on the coupling and mass of the new gauge boson both in the NSI limit [20] and in the long-range regime [16] for several U (1) ′ flavor symmetries.However, strictly speaking, the bounds derived in Ref. [19] cannot be directly used to constraint the U (1) ′ scenarios because in the latter case only flavor diagonal interactions (and only some of them depending on the U (1) ′ charge) are generated, while the bounds derived in Ref. [19] 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) ′ 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) ′ 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) ′ 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) [21] solution first observed in solar neutrinos, where for suitable NSI the data can be explained by a mixing angle θ 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 [9,10,[22][23][24][25].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.

Formalism
We are going to focus on U (1) ′ 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) ′ 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) ′ .In what respects the effect of the new U (1) ′ 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 [20].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 ′ 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) ′ Lagrangian relevant for propagation in ordinary matter has the most general form with arbitrary charges a u,d,e and b e,µ,τ . Model 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. [22] we have defined B y ≡ B 1 − yB 2 − (3 − y)B 3 where y is an arbitrary constant.The results presented in this article are independent of y.
These charges can be accommodated in generalized anomaly free UV-complete models including only the SM particles plus right-handed neutrinos [26].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 µ − L τ ), and (L µ − L e ) [20]: In particular models with charges as in Eq. (2.2) are not constrained by rare Z decays or flavor changing neutral current (FCNC) meson decays [27][28][29].For convenience, in table 1 we list the charges corresponding to some of the models discussed in the literature.
In general, the evolution of the neutrino and antineutrino flavor state during propagation is governed by the Hamiltonian: where H vac is the vacuum part which in the flavor basis (ν e , ν µ , ν τ ) reads 31 ) . (2.4) Here U vac denotes the three-lepton mixing matrix in vacuum [1,30,31].Following the convention of Ref. [32], we define U vac = R 23 (θ 23 )R 13 (θ 13 ) R12 (θ 12 , δ CP ), where R ij (θ ij ) is a rotation of angle θ ij in the ij plane and R12 (θ 12 , δ CP ) is a complex rotation by angle θ 12 and phase δ CP .
Concerning the matter part H mat of the Hamiltonian generated by the SM together with the U (1) ′ interactions in (2.1), its form depends on the new interaction length determined by the Z ′ mass as we discuss next.

2.1
The large (M Z ′ ≳ 10 −12 eV) M Z ′ limit: the NSI regime In the limit of large M Z ′ , the Z ′ field can be integrated out from the spectrum and Eq.(2.1) 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.1).The coefficients ε f,P αβ parametrizes the strength of the new interaction with respect to the Fermi constant, with where we have introduced the notation 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: These interactions lead to a flavor diagonal modification of the matter potential where the "+1" term in the ee entry accounts for the standard contribution, and describes the non-standard part.Here N f (⃗ x) is the number density of fermion f at the position ⃗ x along the neutrino trajectory.In Eq. (2.9) we have limited the sum to the charged fermions present in ordinary matter, f = e, u, d.
, and also that matter neutrality implies N p (⃗ x) = N e (⃗ x), Eq. (2.9) becomes: where and we have introduced the proton and neutron Z ′ couplings As discussed in Ref. [17], in the Earth the neutron/proton ratio Y n (⃗ x) which characterize the matter chemical composition can be taken to be constant to very good approximation.The PREM model [33] fixes Y n = 1.012 in the Mantle and Y n = 1.137 in the Core, with an average value Y ⊕ n = 1.051 all over the Earth.Setting therefore Y n (⃗ x) ≡ Y ⊕ n in Eqs.(2.9) and (2.10) we get E αα (⃗ x) ≡ ε ⊕ αα with: For what concerns the study of propagation of solar and KamLAND neutrinos one can work in the one mass dominance approximation, ∆m 2  31 → ∞ (which effectively means that 31 /E ν ).In this approximation the survival probability P ee can be written as [34,35] P ee = c and where Following Ref. [19] we can rewrite the Z ′ contribution as: where the angle η parametrizes the ratio of the charges of the matter particles as: and The neutrino oscillation phenomenology in this regime reduces to a special subclass of the general NSI interactions analyzed in Ref. [19]. 1 In particular, as a consequence of the CPT symmetry (see also Refs.[17,18,32] for a discussion in the context of NSI) the neutrino evolution is invariant if the relevant Hamiltonian is transformed as H → −H * .In vacuum this transformation can be realized by changing the oscillation parameters as where δ CP is the leptonic Dirac CP phase, and we are using here the parameterization conventions from Refs.[19,32].The symmetry is broken by the standard matter effect, which allows a determination of the octant of θ 12 and (in principle) of the sign of ∆m 2  31 .However, in the presence of the Z ′ -induced NSI, the symmetry can be restored if in addition to the transformation Eq. (2.23), the E αα (⃗ x) terms can be transformed as [18,32,36] (2.24) Eq. (2.23) shows that this degeneracy implies a change in the octant of θ 12 (as manifest in the LMA-D fit to solar neutrino data [21]) as well as a change in the neutrino mass ordering, i.e., the sign of ∆m 2  31 .For that reason it has been called "generalized mass ordering degeneracy" in Ref. [32].Because of the position dependence of the NSI hamiltonian described by E αα (⃗ 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 θ 12 we consider two distinct parts of the parameter space: one with θ 12 < 45 • , which we denote as LIGHT, and one with θ 12 > 45 • , 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ν oscillations in the LIGHT sector.Till recently this was indeed the case because for the last decade the value of ∆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 8 B measurements performed by SNO, SK and Borexino showed no evidence of the low energy spectrum turnup expected in the standard LMA-MSW [4,5] solution for the value of ∆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 ∆m 2  21 value indicated by KamLAND.With the data included in the analysis in Ref. [19,37] this resulted into a tension of ∆χ 2 ∼ 7.4 for 1 To be precise, the data analysis performed in Ref. [19] was restricted to NSI with quarks, ie ae = 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. [19], 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.
the standard 3ν oscillations.Such tension could be alleviated in presence of a non-standard matter potential, thus leading to a possible decrease in the minimum χ 2 .However, with the latest 2970-days SK4 results presented at the Neutrino2020 conference [38] in the form of total energy spectrum and updated day-night asymmetry, the tension between the best fit ∆m 2  21 of KamLAND and that of the solar results has decreased.Currently they are compatible within 1.1σ in the latest global analysis [39].

2.2
The finite M Z ′ case: the long-range interaction regime If M Z ′ is very light the four-fermion contact interaction approximation in Eq. (1.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 ∼ 1/M Z ′ around it.However, following Ref.[12] the generalized matter potential can still be written as Eq.(2.8) provided that Eq. (2.9) is modified as: where Taking into account that ordinary matter is neutral and only contains f = e, u, d, we can rewrite Eq. (2.25) in a way that generalizes Eq. (2.10): 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 ⃗ x ⊙ the center of the Sun and accounting for the spherical symmetry of the matter potential we can write: 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 ⃗ x ⊕ the center of the Earth and by X ⊖ = |⃗ x ⊙ − ⃗ x ⊕ | the Sun-Earth distance, we have: The solar-induced contribution becomes non-negligible when the range of the interactions, 1 M Z ′ , is comparable or larger than the Sun-Earth distance X ⊖ .The factors F ⊙ i (r, M Z ′ ) and F ⊕ i (r, M Z ′ ) represent the modification of the matter potential due to the finite range of the interaction mediated by the Z ′ 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 ⊙(⊕) .Hence: 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 (0, M Z ′ ).For the Earth matter potential the position dependence of the factor F ⊕ i (r, M Z ′ ) is very weak in the current experiments, so one can also scale the contact interaction potential with an approximate F ⊕ i (r, M Z ′ ) evaluated at a fix r which we take to be also r = 0.In Fig. 1 we plot these scale factors F ⊙(⊕) e (0, M Z ′ ).As seen in the figure for M Z ′ ≲ 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) ′ -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 ′ , 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 ′ 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 ′ coupled to L e − L µ , but the results are similar for models with other couplings) for each M Z ′ to those obtained in the large M Z ′ regime.As seen in the figure for M Z ′ ≳ 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 ′ ≲ 10 −13 eV all the matter potentials in the analysis show deviations from the contact interaction regime.
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.(2.23) and (2.24) because of the very different ⃗ x dependence of the SM matter potential and the one generated by the Z ′ .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 ⃗ 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.

Results of the global oscillation analysis
We have performed a global fit to neutrino oscillation data in the framework of 3ν massive neutrinos with new neutrino-matter interactions generated by U (1) ′ models and characterized by the Lagrangian in Eq. (2.1).For the detailed description of methodology and data included we refer to the comprehensive global fit in Ref. [19] 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 [39] which includes the previously cited solar neutrinos 2970-days SK4 results [38], the updated medium baseline reactor samples from RENO [40] and Double Chooz [41], and the latest long-baseline samples from T2K [42] and NOνA [43].Notice that in order to keep the fit manageable we proceed as in in Ref. [19] and restrict ourselves to the CP-conserving case and set δ CP ∈ {0, π}.Consequently the T2K and NOν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 χ 2 function: for each Z ′ model.In general each model belongs to a family characterized by a set of U (1) ′ charges; for each family χ 2 OSC+Z' depends on the two variables parametrizing the new interaction, g ′ and M Z ′ , plus the six oscillation parameters ⃗ ω ≡ (∆m 2 21 , ∆m 2 31 , θ 12 , θ 13 , θ 23 , δ CP ).Following the discussion in Sec. 2, we have performed the analysis in two physically distinctive domains of M Z ′ : • 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 ′ ≳ 10 −10 eV.
In this regime a possible conflict with the cosmological bound on ∆N eff may appear because of the contribution of either the Z ′ itself (if lighter than all active neutrinos) or by the extra contribution to the neutrino density produced by the Z ′ decay (if heavier than some ν mass eigenstate).These bounds can be evaded in two distinct ranges of the U (1) ′ interactions: -M Z ′ ≳ 5 MeV for which the contribution to the neutrino energy density due to the decay of the Z ′ into neutrinos is sufficiently suppressed by the Boltzmann factor ∼ exp(−M Z ′ /T ) [44]; -M Z ′ ≲ O(eV) but with very weak coupling g ′ < 10 −10 for which the Z ′ is produced through freeze-in.In this regime, even if Z ′ could decay to the lightest neutrinos this would happen after neutrino decoupling making the contribution to ∆N eff negligible or at most within the present allowed range [45].
• Long-Range Domain (DOM=LRI): If M Z ′ ≲ 10 −13 eV interactions in the Earth and Sun matter are long range, and for a given value of g ′ and M Z ′ the effects in the Earth are suppressed with respect to those in the Sun.As mentioned above, for g ′ < 10 −10 the contribution to ∆N eff is negligible.
Correspondingly we define In both domains we compare the results of the fit including the new U (1) ′ interaction with those obtained in the "standard" 3ν-mixing scenario, which we will denote as "OSC", and for which the present global fit yields χ 2 OSC,min = 718.5 . (3.3) We will classify the models according to the quality of the fit in the presence of the U (1) ′ interactions compared with that of OSC by defining where by | 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 where by | 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) ′ interactions characterized by the six couplings a u,d ∈ {−1, 0, 1} and a e , b e,µ,τ ∈ {−3, −2, −1, 0, 1, 2, 3}.In this way our survey covers a total of ∼ 10000 different sets of U (1) ′ charges which can produce effects in matter propagation in neutrino oscillation experiments.

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 ′ and M Z ′ .We find that in the NSI (LRI) domain 88% (90%) of the surveyed sets of charges lead to a decrease in the χ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 ′ models with gauging of SM global symmetries with SM plus right-handed neutrinos matter content of Eq. (2.2).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 ∆χ 2 LIGHT,DOM for U (1) ′ 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σ level.
Quantitatively we find that in the NSI regime none of the surveyed models yields an improvement beyond -1.9 units of χ 2 .This is the case, for example, of a Z ′ coupled to B − 3L e + 2L µ + 3L τ .Generically models in the LRI domain can provide a better fit with a reduction of up to 3.5 units of χ 2 .For example a model with charge L e + 2L µ − 3L τ and a Z ′ with M Z ′ ∼ 5 × 10 −15 eV provides a better fit than standard oscillations by 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σ level, this is min g ′ ,M Z ′ (∆χ 2 LIGHT,DOM ) was always larger than −4. 2 Consequently for all models sur-Model Table 3. Results for the models with charges in table 1 in the NSI regime.For the model in Ref. [22] we have defined The second column gives minimum ∆χ 2 defined w.r.t. the 3ν oscillation (see Eq. (3.4)).The last column gives the coefficient of the bound on the coupling over the mediator mass in units of 100 MeV.
veyed one can conclude that the analysis of neutrino oscillation experiments show no significant evidence of U (1) ′ 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 χ 2 (4 units for 95.45% CL), this is: 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 ′ and M Z ′ for which the U (1) ′ model characterized by such interaction strength and interaction length, gives a fit which is worse than standard oscillations by at least 4 units of χ 2 .As in the LIGHT sector the standard model is recovered for either g ′ → 0 or M Z ′ → ∞, the above condition yields also the 2σ excluded one-dimensional upper range of interaction coupling g ′ for each value of the interaction length (or, correspondingly, the 2σ excluded one-dimensional lower range of M Z ′ for each value of the interaction coupling), for all models characterized by a given set of charges.The corresponding excluded ranges for the Z ′ coupling and mass for the models with couplings listed in table 1 are shown in Figs. 2 and 3 for M Z ′ below 1 eV and in the O(MeV-GeV) range, respectively 3i.e., below and above the window strongly disfavored by the cosmological bound on ∆N eff .In particular in Fig. 2 we observe the slope change of the oscillation exclusion ranges for masses M Z ′ ∼ 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 ′ ≲ 10 −15 eV which 3 Tables with the numerical values of the bounds can be provided upon request to the authors.
Values of g ′ and M Z ′ ≤ O(eV) for which a U (1) ′ model coupled to the charges labeled in each panel gives a worse fit than standard oscillation by 2σ, Eq. (3.6) (hatched region).We also show the bounds on these models imposed by gravitational fifth force searches [50,51] and by equivalence principle tests [52].In the window corresponding to the model in Ref. [22] we have defined B y ≡ B 1 − yB 2 − (3 − y)B 3 .See text for details.
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 ′ mass and coupling derived at comparable confidence level under the minimal assumptions in Eq. (2.1).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 ν a → ν b Z ′ from Cosmic Microwave Background (CMB) data [46], bounds from Black-Hole superradiance [47], or from flavor composition of extra-galactic neutrinos [48].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 ′ in meson decays, neutrinoless ββ decay, or neutrino annihilation in Big-Bang Nucleosynthesis (BBN) and supernovae are relevant for Z ′ with lower masses than those shown in the figure (see, for example, Ref. [49]).
The bounds imposed by gravitational fifth force searches shown in Fig. 2 were obtained by rescaling the results shown in Ref. [50] (which, in turn, were recasted from Ref. [51]).Limits from equivalence principle tests are obtained rescaling the results from Ref. [52].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. (3.6).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 ′ lighter than ∼ 10 −11 eV.
Conversely the results shown in Fig. 3 for M Z ′ in the O(MeV-GeV) range correspond to U (1) ′ effects in oscillation experiments always in the NSI domain.In this case, the analysis of oscillation data results in a bound on g ′ 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) ′ 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 ′ and M Z ′ which are only constrained by the oscillation analysis.This is particularly the case for U (1) ′ coupled to B−3L τ .We have found no competitive bound from other experiments in the shown window with the only exception of the constraint reported in Ref. [53] using data from DONUT [54] and from NA62 (for bounds from other experiments relevant for larger couplings see for example Refs.[20,55]).to standard oscillations, irrespective of the U (1) ′ 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 ′ and M Z ′ .
• The excluded ranges for the Z ′ coupling and mass for the models with couplings listed in table 1 are shown in Figs. 2 and 3 for M Z ′ below 1 eV and in the O(MeV-GeV) range, respectivelyi.e., below and above the window strongly disfavored by the cosmological bound on ∆N 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 ′ lighter than ∼ 10 −11 eV.
• For mediators in the O(MeV-GeV) range we list in table 3 the derived constraints on g ′ versus M Z ′ .We find that for several of the considered models there are values of g ′ and M Z ′ 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 ∆χ 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.(2.2)).So, generically, Z ′ models for LMA-D with gauged SM global symmetries require additional states for anomaly cancellation.
where G is the gravitational constant, Ni 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. [52] 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 NBe = 0.99868 and NTi = 1.001077, as in Ref. [52].
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. [52] put a bound on which is used to set a constraint on α as a function of λ, α max , λ . (A.9) This can be used to set a bound on ∆V ′ (r) for a general model, noting that where u stands for the atomic mass unit in GeV (u ≃ 0.931 GeV).So the boundary in the g ′ vs M Z ′ plane will be

Bounds from white-dwarf cooling
We based the bounds shown in Figs. 3 and 4 on the study presented in Ref. [62].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) ′ scenarios here considered the minimum new contribution is due to Z ′ mediated decays into neutrinos where C e,V is the vector coupling to the electron current in the SM, Z s is the plasmon wavefunction renormalization and π s is the effective plasmon mass which enters in the dispersion relation ω 2 s − k 2 = π s (ω s , k).Here, s = T, L refers to the plasmon polarization (transverse or longitudinal).
Under the assumption that the mass of the Z ′ is much larger than the frequency of the plasmon we can write its rate into neutrinos of a given flavor β due to the new interactions as (a e b β ) 2 48π 2 α em Z s π 3 s ω s (A.13) • Constraints on invisible Z ′ decays in searches for π 0 → γ Z ′ at NA62 experiment [90] and in e N → e N Z ′ at NA64 experiment [91][92][93].There are bounds from NA64 for the Z ′ decaying into e + e − [94], 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) ′ models considered we have used the darkcast software [56], 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) ′ models coupled to charges including L µ can be constrained with data on production of µ + µ − in ν µ scattering off the Coulomb field of a nucleus (the so called neutrino trident production) at CHARM-II [95] and the CCFR experiment [96], see Ref. [97].They can also be bounded with data on µ + µ − production in e + e − collision at BaBar [98] and Belle [99].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
The impact of light Z ′ 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 ′ emission can be used to set very tight constraints on these models.Solar and stellar constraints are only relevant in the mass window eV ≲ M Z ′ ≲ 100 eV (see, e.g., Refs.[16,59,61,100]) 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 ′ case.As an example, in Ref. [101] a specific analysis carried out for supernovae emission for the L µ − L τ and B − L models showed large variations in the results with respect to the region constrained in the dark photon case [102,103].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 ′ onto e + e − in the early Universe, which would be applicable in the mass region above 1 MeV (see, e.g., Fig. 6 in [104]).However, again in this case it is uncertain how to recast these bounds to a general U (1) ′ model with arbitrary charges, and therefore we have not included them here.

Figure 1 .
Figure 1.Effective suppression factor of the matter potential due to the long range of the U (1) ′ interactions as a function of M Z ′ .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).

Figure 3 .
Figure 3. Values of g ′ and 5MeV ≤ M Z ′ ≤ 10 GeV for which a U (1) ′ model coupled to the charges labeled in each panel gives a worse fit than standard oscillation by 2σ, Eq. (3.6) (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.[22] we have defined B y ≡ B 1 − yB 2 − (3 − y)B 3 .See text for details.

Table 2 .
[22]lts for the models with charges in table 1 in the LRI regime.For the model in Ref.[22]we have defined B y ≡ B 1 − yB 2 − (3 − y)B 3 .The second column gives minimum ∆χ 2 defined w.r.t. the 3ν oscillation (see Eq. (3.4)).The last column gives the the upper bound for the coupling of asymptotically for ultra light mediators, M Z ′ ≲ 10 −15 eV