Light, Long-Lived $B-L$ Gauge and Higgs Bosons at the DUNE Near Detector

The low-energy $U(1)_{B-L}$ gauge symmetry is well-motivated as part of beyond Standard Model physics related to neutrino mass generation. We show that a light $B-L$ gauge boson $Z{'}$ and the associated $U(1)_{B-L}$-breaking scalar $\varphi$ can both be effectively searched for at high-intensity facilities such as the near detector complex of the Deep Underground Neutrino Experiment (DUNE). Without the scalar $\varphi$, the $Z{'}$ can be probed at DUNE up to mass of 1 GeV, with the corresponding gauge coupling $g_{BL}$ as low as $10^{-9}$. In the presence of the scalar $\varphi$ with gauge coupling to $Z{'}$, the DUNE capability of discovering the gauge boson $Z{'}$ can be significantly improved, even by one order of magnitude in $g_{BL}$, due to additional production from the decay $\varphi \to Z{'}Z{'}$. The DUNE sensitivity is largely complementary to other long-lived $Z{'}$ searches at beam-dump facilities such as FASER and SHiP, as well as astrophysical and cosmological probes. On the other hand, the prospects of detecting $\varphi$ itself at DUNE are to some extent weakened in presence of $Z{'}$, compared to the case without the gauge interaction.


Introduction
The U (1) B−L gauge symmetry is a well-motivated extension [1,2] of the Standard Model (SM), which provides a natural framework to account for the tiny neutrino masses via the type-I seesaw mechanism [3][4][5][6][7]. It may also accommodate a light dark matter (DM) particle which interacts with the SM particles through the scalar or gauge portal of U (1) B−L [8]. In this paper, we focus on the class of B − L extensions of SM, where B − L does not contribute to electric charge so that its gauge coupling g BL can be arbitrarily small making the associated gauge boson Z very light with mass in the sub-GeV range [9][10][11]. This is especially true if the U (1) B−L -breaking scale is low, e.g., below the GeV-scale. In this case, we can expect the corresponding U (1) B−L -breaking scalar (denoted here by ϕ) to be also light. Both Z and ϕ are essential ingredients of the U (1) B−L extension, and play important role in connecting DM to the SM, in DM phenomenology [8,[12][13][14][15][16][17], and in explaining the observed light neutrino masses [18][19][20][21][22][23][24]. Since the seesaw mechanism does not depend directly on the gauge coupling value, such low gauge coupling models can still accommodate the seesaw mechanism. For phenomenological exploration of this class of B − L extensions at a high mass scale and larger gauge coupling range, see, e.g., Refs. [25][26][27][28][29][30].
If kinematically allowed, both ϕ and Z can be produced and detected in the high-intensity experiments such as the Deep Underground Neutrino Experiment (DUNE) [31]. One or both of these particles may be relatively long-lived and be able to travel from the DUNE target to its near detector complex where it decays in a striking fashion. It may then be possible to search for signatures of this class of beyond Standard Model (BSM) physics in such experiments. Our goal in this work is to determine how these gauge bosons, scalar bosons, and their interplay can be searched for in the DUNE experiment at Fermilab. In particular, we will focus on the following three scenarios: (i) pure Z case with the scalar ϕ decoupled, (ii) the effect of ϕ → Z Z on the prospects of ϕ at DUNE, and (iii) improvement of the DUNE sensitivity of Z boson due to the scalar source ϕ → Z Z .
The Z boson in the U (1) B−L model couples to all the SM quarks and leptons and thus can be produced from meson decays such as π → γ + Z and η → γ + Z as well as other hadronic decays where the π 0 and η etc are final state particles in the pp collision [32]. The proton-proton bremsstrahlung process is also very important, in particular when the Z boson is heavier than the η meson [33,34] (cf. Fig. 3). After being produced, the Z boson may decay via Z → e + e − , µ + µ − , νν and hadrons such as πππ and KK [35][36][37], depending on its mass, and in the small gauge coupling limit can lead to displaced vertices with decays in the DUNE near detector hundreds of meters away. The observation of such displaced vertices may provide the key signatures of the light Z boson. An analogous analysis of a leptonic gauge boson (e.g. one of U (1) Lµ−Lτ ) has been performed in Ref. [31]. Such gauge bosons can lead to anomalous neutrino-electron scattering and neutrino trident events, which have been explored in Refs. [38,39]. We will perform a thorough study of the B − L Z gauge boson and determine the DUNE near detector sensitivity. It turns out to be qualitatively different from the leptonic Z boson, as presented in Fig. 4.
As for the light scalar ϕ, if it mixes with the SM Higgs boson, it will obtain loop-level flavor-changing neutral current (FCNC) couplings to the SM quarks. As a result, it can be produced from the loop-level FCNC meson decays such as K → π + ϕ at DUNE [40], and then decay into the light SM particles ϕ → e + e − , µ + µ − , ππ, γγ, with the decay also induced by the mixing of ϕ with the SM Higgs. The sensitivity of a light scalar ϕ without the Z boson at DUNE has been studied in Ref. [31].
In the U (1) B−L model, there is the gauge coupling ϕZ Z , This will induce the extra decay channel ϕ → Z Z , which could have multiple effects on the prospects of ϕ and Z at DUNE: • This new decay channel of ϕ will produce more Z bosons at DUNE, compared to the pure Z case, even by a factor of 10 5 (see Fig. 6). As a result, the Z can be probed at DUNE with a smaller gauge coupling g BL , with an improvement factor up to 45 (see Fig. 7). This is only possible for a Z mass m Z 200 MeV. This is because we require m ϕ > 2m Z , and ϕ heavier than ∼400 MeV cannot be produced from the meson decay K → π + ϕ at DUNE.
• The gauge portal scalar decay will not only change the branching fractions of ϕ (cf. Fig. 2), but also shorten its lifetime. This will have some effects on the search for ϕ at DUNE, depending on the BR of ϕ decaying into SM particles, i.e. BR(ϕ → SM). As presented in Fig. 5, the effect of ϕ → Z Z is most significant when ϕ is relatively heavy and the mixing angle sin ϑ is relatively large.
• The simultaneous existence of ϕ and Z will also induce the associated production of ϕ and Z at DUNE, for instance from the meson decay K → π + Z + ϕ. However, no matter whether the scalar ϕ is emitted from the Z boson line or from couplings to mesons, such decays will always be highly suppressed by either g 4 BL or sin 2 ϑg 2 BL , and can thus be neglected at DUNE.
The rest of this paper is organized as follows. First, we provide a brief review of the model in Section 2. Then, we begin the discussion of search strategies in Section 3 with a brief discussion of the experimental setup and how proton-beam experiments (including the DUNE beamline) are well-suited for searches of these types of models. We divide the experimental search for this model by taking the following approach: first, we discuss in Section 3.1 how the Z gauge boson of U (1) B−L can be searched for in DUNE making agnostic assumptions about the existence of the associated scalar ϕ. Then, in Section 3.2, we will demonstrate that searches for ϕ in accelerator beam environments, including DUNE, are weakened in the presence of Z , where the produced ϕ can decay into Z instead of SM particles. We will also revisit the Z search strategy and discuss how additional fluxes sourced by ϕ (including a brief discussion of associated production of both ϕ and Z ) decays can be detected in Section 3.3. This improves the capability of discovering the gauge boson Z . Finally, we conclude in Section 4. Some details of the ϕ decay calculations are relegated to Appendix A. Subdominant contributions to Z flux at DUNE are collected in Appendix B, and more Z limits can be found in Appendix C.

The Model
We consider the minimal U (1) B−L model based on the gauge group SU (2) L × U (1) Y × U (1) B−L . For the purpose of anomaly cancellation, three right-handed neutrinos (RHNs) N i (with i = 1, 2, 3) are naturally introduced. To break the U (1) B−L gauge symmetry, a complex singlet scalar field φ is introduced, which carries two units of B − L charge. When the φ-field develops a non-vanishing vacuum expectation value After the symmetry breaking, expanding the field around its VEV, we obtain The CP-odd component χ is "eaten" by the Z boson, while the CP-even component ϕ is left as a physical scalar field. In presence of the (H † H)(φ † φ) term in the scalar potential (H is the SM Higgs doublet), the scalar ϕ mixes with the SM Higgs h, with a mixing angle sin ϑ. Given the Yukawa coupling with C standing here for charge conjugate, the RHNs obtain masses m N = y N v BL / √ 2, which can be used to generate the tiny active neutrino masses via type-I seesaw mechanism [3][4][5][6][7]. For simplicity, we will assume the RHN masses m N > m Z /2 and m N > m ϕ /2 such that the decays Z → N N and ϕ → N N are kinematically forbidden, and we can neglect the contributions of RHNs to the DUNE prospects of Z and ϕ in this paper. We also assume for simplicity that there is no Z − Z mixing throughout this paper, which would otherwise potentially affect the decay branching ratios (BRs) of the Z boson. A small, loop-induced kinetic mixing between Z and the SM Z will arise, on the order of g BL g/(16π 2 ). Since we are interested in g BL 1, this mixing is always very small, and we will disregard it for the remainder of this work. We are also going to focus on the mass regime MeV m Z GeV -whether such small gauge couplings and light gauge boson masses are natural or fine-tuned requires a detailed analysis of the scalar potential, which is beyond the scope of this manuscript, and we just take a more phenomenological approach.
Then there are only four free phenomenological parameters in the minimal which we study to investigate the DUNE prospects. To produce the Z boson and/or the ϕ scalar at the DUNE experiment, their masses m Z and m ϕ have to be at or below the GeV-scale. As we will show in Figs. 4 and 5 respectively, the gauge coupling g BL and the mixing sin ϑ have to be sufficiently smaller than one to satisfy the current experimental constraints.  Neglecting the RHN channel, the gauge boson Z decays predominantly into pairs of SM fermions (and mesons), and we can express the decay widths of Z into SM fermions as [36] Γ with N C f being the color factor (3 for quarks and 1 for leptons), and the symmetry factor S f = 1 for quarks and charged leptons and 1/2 for neutrinos. Decay widths into hadrons from the quark hadronization are more complicated, but can be expressed using the R-ratio (see, e.g., Refs. [35][36][37] for further discussion). In practice, we use the darkcast code from Ref. [36] to determine the Z lifetime and BR into visible final states accessible in the DUNE near detector. These final states include e + e − , µ + µ − , and hadronic ones, dominated by π + π − π 0 and K + K − states (but not the two-pion final state, which is forbidden by G-parity conservation in the absence of isospin-breaking terms); Fig. 1 presents these BRs as a function of m Z .
Through mixing with the SM Higgs, the scalar ϕ decays into the SM leptons, pion pairs and two photons via the SM W boson and charged fermion loops. The corresponding partial widths are respectively where G F is the Fermi constant, α ≡ e 2 /4π is the fine-structure constant, and Q f is the charge of the fermion f in units of the proton charge e. In Eq. (2.6), G(m 2 ϕ ) is a dimensionless transition amplitude defined in Appendix A. Similarly in Eq. (2.7), A 1/2 (τ f ) and A 1 (τ W ) respectively (with τ X = m 2 ϕ /4m 2 X ) are standard fermion and W loop functions for the Higgs decay, also given in Appendix A for completeness.
As a result of the gauge coupling of ϕ to Z bosons, if kinematically allowed, we have also the decay channel As shown in Eqs. (2.5) -(2.7), all the decay channels of ϕ into the SM particles are universally proportional to the mixing angle sin ϑ, while the Z boson channel in Eq. (2.8) is dictated by the gauge coupling g BL . In the limit of g BL → 0 or m ϕ < 2m Z , the scalar ϕ only decays into the SM particles, which is equivalent to the case of singlet extension of the SM. In presence of the U (1) B−L gauge couplings and the Z boson, the decay products of ϕ are changed dramatically, depending on the values of m Z , sin ϑ and g BL . The BRs of ϕ for the following two Benchmark Points (BP) of m Z and g BL are shown in Fig. 2:

Searches for Scalars and Gauge Bosons at DUNE
The Fermilab DUNE facility [41,42] uses a 120 GeV proton beam 1 striking a graphite target to generate the high-intensity neutrino beam. Of interest for this work, the proton collisions produce various hadronic final states such as pions, Kaons, η mesons, etc. Specifically, we are focused on neutral and charged pions, Kaons and η mesons. In their decays, these can produce the new BSM particles Z and ϕ. The proton beam is expected to deliver at least 1.47 × 10 21 protons on target (POT) per year, with potential upgrades throughout the lifecycle of the experiment. We will assume a nominal exposure of 1.47 × 10 22 POT for our analyses, conservatively corresponding to at most ten years of data collection. The DUNE near detector complex [43], including several detector components, is located 574 meters away from the target and therefore when the decays of the new particles Z and ϕ occur far away from the production vertex, the final states can give a signal in the DUNE near detector. This setup is ideal in probing the parameter space in which g BL and the new particle masses are small, and therefore they are long-lived.
As detailed in Section 3.1.1, in our U (1) B−L model, π 0 , η → γ + Z provides the source for the new BSM particle Z . This decay will take place promptly, effectively in the DUNE target, and will be our source for Z . There can also be direct bremsstrahlung production of Z from pp → pp + Z since protons have nonzero B − L quantum number. The ϕ production arises mostly from the decays K + → π + + ϕ and K L → π 0 + ϕ. Whether we are looking for ϕ or Z signatures in the near detector, we will be interested in final states including opposite-sign pairs of charged particles. This includes e + e − , µ + µ − , and pionic final states. The DUNE capability, specifically using its gaseous argon near detector (immediately downstream of the liquid argon near detector) to identify final-states of particles decaying in flight has been detailed in Refs. [31,[44][45][46][47][48][49]. We will use the results of these works and consider that a background-free search for these final states is possible.
As discussed in Section 1, our search strategy naturally divides into three different categories based on the (non)existence of either the scalar ϕ and gauge boson Z . Phenomenologically speaking, we are interested in the four free parameters in Eq. (2.3). We divide the search strategies by considering the following cases: Figure 2. BRs of the scalar ϕ in the U (1) B−L model into ee, µµ, γγ, ππ (including both π 0 π 0 and π + π − ) and Z Z , as a function of the scalar mass. In both panels, we take g BL = 5 × 10 −9 and sin 2 ϑ = 10 −8 . The top (bottom) panel assumes m Z = 30 MeV (140 MeV).
• Pure Z case: sin ϑ = 0 -if a scalar boson exists at all, it is decoupled and does not have any impact on the phenomenology of Z . This is explored in Section 3.1.
• Pure ϕ case: g BL = 0 -if a gauge boson exists at all, it is decoupled and does not impact the phenomenology of ϕ. This has been studied in Ref. [31], and will not be detailed in this paper any more.
• Combined Z and ϕ case (effects of ϕ → Z Z ): If both particles are relevant at DUNE and m ϕ > 2m Z then the decay channel ϕ → Z Z can impact the prospects of the DUNE search for ϕwe investigate this effect for DUNE and other ϕ searches in Section 3.2.2. Furthermore, this can also provide additional flux of Z at the DUNE near detector. This additional Z flux allows for increased sensitivity as a function of m Z and g BL , subject to current constraints on m ϕ and sin ϑ. This is discussed in Section 3.3 and is the main emphasis of this paper. As mentioned in the introduction, the associated production of both Z and ϕ in some decay channels will be highly suppressed either by g 4 BL or g 2 BL sin 2 ϑ, and can be neglected.

Pure Gauge Boson Search
Here we discuss the case in which the only new physics particle is the gauge boson Z with a mass m Z and gauge coupling g BL . In the context of a search in a proton beam-dump environment, this scenario is very similar to the study performed for a U (1) dark photon in Ref. [31], with a mapping of the kinetic mixing parameter ε onto this gauge coupling. This remapping is discussed for proton beam-dump experiments, as well as other experimental contexts, in Ref. [37]. In the following subsections, we discuss Z production in Section 3.1.1 and DUNE sensitivity in Section 3.1.2; expressions for the decay of Z can be found in Section 2.

Production of Gauge Bosons
The new gauge boson Z can be produced in the same ways as the U (1) dark photon -decays of neutral pseudoscalar mesons m → γZ , and proton bremsstrahlung pp → ppZ . We give expressions for the flux of these Z at the DUNE near detector here. The neutral pseudoscalar meson production flux is determined using The quantities in Eq. (3.1) are as follows: c m is the average number of mesons m produced in a given POT collision; N POT is the POT number considered in the experimental exposure; A Det. is the detector area as viewed by an incoming particle; (m Z ) is the geometrical acceptance factor of the Z particles emerging from the decay m → γZ , determined using Monte Carlo simulations; e is the electric charge, m m the mass of m, and BR (m → γγ) is the SM BR of the meson m into two photons. We find that the dominant contribution to Φ mZ comes from π 0 and η mesons. For proton-proton bremsstrahlung, we use the calculations of Refs. [33,34] with specifics for DUNE discussed in Ref. [31]. The total flux from this bremsstrahlung process is Here, σ pN is the total proton-target cross section evaluated at s = 2E p m p with E p = 120 GeV being the DUNE beam energy. F 1,N is a form factor allowing for mixing between Z and the SM vector mesons, specifically when m Z approaches m ρ . The two integrals are performed over the variables p 2 T (the transverse momentum squared of the outgoing Z ) and z, the fraction of the incoming proton's initial momentum that is transferred to the longitudinal momentum of the outgoing Z , where we label "Det." to indicate that the range of integration of these two is such that the outgoing Z is pointing from its production point to the front face of the near detector. The photon splitting function w ba is given in Refs. [31,33,34] and, for the case of a U (1) B−L gauge boson, is proportional to g 2 BL /4π. Other production mechanisms for Z are possible, however we find that all of them are subdominant to the neutral pseudoscalar meson decays and proton-proton bremsstrahlung. These include K L → γZ , three-body decays of charged mesons π + , K + → + νZ (with = e, µ), loop-level flavor-changing decays such as K → πZ (including both neutral and charged Kaon decays), and heavy baryon decay ∆ → N Z . For completeness, we give the BR expressions in these channels in Appendix B.  As discussed earlier in this Section, we will assume ten years of data collection. 2 Assuming only Z exists (either ϕ does not exist, sin ϑ → 0, or m ϕ < 2m Z ), we can determine the flux of Z at the DUNE Near Detector complex using the above. This is shown, divided by different contributions, in Fig. 3.
Given the flux of Z , we may determine the number of decay events in the DUNE Near Detector volume using The outermost integral is over the energy of the Z bosons and can account for detector thresholds, efficiencies, etc. (with E min and E max respectively the minimal and maximal energy of Z bosons) -in our simulation we allow for the full range of E Z and assume 100% signal identification efficiency. The second integral is over the surface area of the detector, 2.5 m in radius for the gaseous argon time projection chamber. The innermost integral is over the depth of the detector L Det. = 5 m with D Det. being the distance to the detector, 579 m 3 . The probability of a decay occurring at position x is In this decay expression, γ Z = E Z /m Z is the Lorentz boost factor and τ Z is the proper lifetime of Z , which can be obtained using the expressions for its width given in Section 2.
The form of Eq. (3.3) assumes that the entire Z flux dΦ Z /dE Z originates at a common location, assumed to be at the proton/target interaction location. If the Z flux is produced continuously along the beam pipe, Eq. (3.3) will require an additional integral over the Z production distribution. This will be relevant when we consider Z coming from ϕ → Z Z decay in Section 3.3. With Eq. (3.3), we now can obtain the sensitivity of DUNE to the pure-Z scenario.

Sensitivity to Pure Gauge Boson Scenario
Here, we detail the DUNE sensitivity to the pure-Z scenario, in which the scalar ϕ either does not exist or its decay into Z pairs is rare. The sensitivity is depicted in Fig. 4, where the solid black line demonstrates the DUNE discovery capability: in the region within the black line, at least ten signal events of Z → + − are expected in ten years of data collection. As argued in Ref. [31], this constitutes a statistically significant signal, especially since the invariant mass of m Z can be reconstructed by the fully visible decay. The DUNE sensitivity in the visible channel dies off for m Z < 2m e , because the Z in that case dominantly decays to invisible final states of light neutrino pairs. 4 In the pure-Z scenario, we see that DUNE can extend on existing beam-dump searches (shown as a purple region, collected from Ref. [37]). This includes the E141 [63][64][65][66], Orsay [67], NuCal [33,68,69], E137 [64], and LSND [70] experiments. In contrast, DUNE can reach heavier m Z (due to large η meson production and bremsstrahlung contributions) and lower g BL (due to the large target-detector distance). We also show the projected sensitivity of the FASER experiment near the LHC in its first run (dashed purple) and proposed FASER2 run (dot-dashed purple) [62]. FASER is sensitive to similar m Z with larger g BL , since the Z produced near the LHC collision point have larger boosts (and therefore longer lab-frame lifetime) than in the DUNE environment. In this sense, DUNE as a detector searching for decaying new-physics particles is highly complementary to FASER. In the time between now and DUNE data collection, other Fermilab-based neutrino experiments could reach into this parameter space as well, including the Short-Baseline Neutrino facility detectors. See, e.g., Ref. [71] for a study on searches for dark scalars in this setup.
For relatively larger gauge coupling and Z masses, there are limits on m Z and g BL from the proton beam-dump experiment CHARM [53], the searches of Z in the channel e + e − → γZ → γ + − (with = e, µ) at BaBar [54], and the neutrino scattering experiment TEXONO [51,52], which are shown respectively as the green, cyan and red shaded regions in Fig. 4. Belle II can improve on the BaBar limits [72], extending to ever lower g BL as it collects data [73]. More limits can be found in Refs. [37,74]. For completeness, we also present our calculated limits on K + → π + Z → πνν from E949 [55] and NA62 [56], which are however weaker than the limits from TEXONO, CHARM and NA62, as indicated by the purple and orange lines in Fig. 4. The future NA62 data can improve the measurement of BR(K + → π + νν) by roughly one order of magnitude [61], thus extend on the existing limits from CHARM and BaBar, as denoted by the dashed orange line in Fig. 4. More details of the E949 and NA62 limits and the prospects at future NA62 can be found in Appendix C.
We note in passing that although a light Z boson of U (1) B−L generates a positive contribution to the muon anomalous magnetic moment [75], the 2σ parameter space preferred by the (g − 2) µ anomaly [59,60] is deep inside the exclusion region in Fig. 4, as shown by the gray band. Thus, if the (g − 2) µ anomaly turns out to be a true signal of BSM physics 5 , the minimal U (1) B−L model with a flavor-universal Z coupling needs to be extended to explain this result; see, e.g., Refs. [78][79][80][81][82][83][84][85][86].
Finally, at low masses and small couplings, constraints from Big Bang Nucleosynthesis (BBN) and the observation of Supernova 1987A (SN1987A) apply, which are shown respectively as the blue and orange  Figure 4. Prospects of the Pure Z search without the scalar ϕ using the DUNE Near Detector (black) where the black lines encompass parameter space for which at least ten signal events (zero background assumed) are expected. Current limits from beam-dump experiments (shaded purple) from Ref. [37], TEXONO (shaded red) [51,52], CHARM (shaded green) [53], BaBar (shaded cyan) [54], E949 (purple line) [55] and NA62 (orange line) [56], and BBN (shaded blue) [57] are also shown. The dark (solid) and faint (dashed) orange regions are respectively the conservative and aggressive SN1987A limits [58]. The gray band denotes the preferred region for muon g − 2 anomaly at the 2σ C.L. [59,60]. Also shown are the prospects from future NA62 data (dashed orange line) [61], FASER (dashed purple line) and FASER2 (dot-dashed purple line) [62]. See text for more details.
shaded regions in Fig. 4. In the early universe, the Z boson may decay into neutrinos and electrons (if kinematically allowed) and thus contribute to the evolution of neutrinos and photons at the MeV scale. This could spoil the success of BBN and thus is constrained by the precise measurement of the light degrees of freedom N eff by Planck [87]. However, the N eff limits on the B − L Z boson depend largely on the Z mass and the coupling g BL , which dictates whether the Z boson could reach equilibrium with the SM particles at the MeV scale. Moreover, neutrino oscillations are also important for the N eff limits [88,89]. The dedicated calculations of the cosmological limits are beyond the main scope of this paper. For simplicity, we take the BBN limit from Ref. [57]. Other relevant studies can be found, e.g., in Refs. [24,[90][91][92]. The supernova constraints (dark orange region/solid lines corresponding to conservative constraints, and faint orange region/dashed lines corresponding to less conservative) come from Ref.
[58] -we direct the reader to this reference, as well as Ref. [57,93,94] for more discussion on supernova luminosity constraints on U (1) B−L vector bosons and associated models. We conclude this discussion by noting that both the BBN and SN1987A limits are model dependent -for example, the chameleon effect (due to the environmental matter density) can weaken astrophysical limits [95], and late reheating can open up parameter space with respect to BBN constraints [96]. This model-dependence will be more relevant when we revisit this parameter space in Section 3.3 where we perform a combined search for Z and the scalar boson ϕ -the interplay of these two new particles and their associated mass scales complicates the calculations of the references mentioned above. We leave a detailed analysis of astrophysical constraints on the combined Z , ϕ scenario to future work.

Effects of Gauge Coupling on Searches for Scalar ϕ
Now we shift focus to the prospects of detecting the scalar boson ϕ responsible for breaking the U (1) B−L symmetry. Its couplings to SM particles are controlled by the mixing angle which we refer to as sin ϑ. Section 2 detailed the different decay channels and widths of ϕ. In Section 3.2.1 we discuss the production of ϕ in the DUNE environment. Section 3.2.2 then explores how the existence of both ϕ and Z (and the decay channel ϕ → Z Z ) can hinder current and future searches for the scalar boson ϕ. The improvement of Z prospects at DUNE as a result of the ϕ → Z Z source will be investigated in Section 3.3.

Production of Scalar Boson
From mixing with the SM Higgs, the light scalar ϕ can be produced in/near the DUNE target in the following channels: • Loop-level flavor-changing meson decays K ± → π ± ϕ and K L → π 0 ϕ. This is identical to the production in Fig. 6.1 of Ref. [31], with the widths [97] Γ(K ± → π ± ϕ) m K ± |y sd | 2 sin 2 ϑ 64π where the loop-level FCNC coupling [40] which is induced by the W −top loop, and v EW = ( √ 2G F ) −1/2 / √ 2 174 GeV is the electroweak VEV. Note that the partial decay widths in Eqs. (3.5) and (3.6) are almost identical, except for the crucial difference that the decay K L → π 0 ϕ depends only on the real part of the coupling y sd .
• The bremsstrahlung process pp → ppϕ is analogous to the pp → ppZ process discussed in Eq. (3.2).
Ref. [98] found that this production process is dominant for such a search at LSND (proton beam energy of 800 MeV). In contrast to the Z bremsstrahlung process, the scalar coupling is suppressed by the effective proton Yukawa coupling. Ref. [99] found that this coupling is small, O(10 −3 ), and Ref. [71] used this result to find that the proton bremsstrahlung process is highly suppressed relative to the Kaon decay processes discussed above for scalar production in the NuMI beam environment (proton beam energy of 120 GeV). We have verified that it is also highly suppressed in the DUNE environment.
In Eq. (3.3) we provided the expected number of decay signal events using some flux of (potentially) decaying particles. The same formalism applies here with the ϕ flux instead of Z , however, as alluded to previously, an extra integral must be applied. Because the Kaons K ± and K L are long-lived, the production of ϕ is not all in/near the DUNE target -the distribution of ϕ production points spans hundreds of meters 6 .

Effects on ϕ prospects at DUNE
The decay ϕ → Z Z not only opens up a new decay channel, but also contributes to the total width of ϕ. Therefore, this extra channel will affect any existing search for ϕ by (a) modifying its branching fraction into visible final-state particles and (b) shortening its lifetime. In general, the constraints on ϕ from searches for decays like K + → π + νν (which could be mimicked by K + → π + ϕ and ϕ decaying invisibly/outside a detector) and beam-dump searches like that of DUNE can be understood by two pieces. At the bottom of these exclusion regions, the reach of an experiment is limited by exposure (and/or backgrounds) -for lower values of sin 2 ϑ, not enough ϕ can be produced to create a signal. At the top of these exclusion regions, the lifetime of ϕ is so short that it decays away before leaving an imprint in a detector. This is particularly relevant for beam-dump searches like CHARM, where the target (production point of ϕ) and detector are well separated.
We reproduce existing limits on m ϕ and sin 2 ϑ coming from a variety of experiments -CHARM [53,100], KOTO [101], NA62 [56], and E949 [55] -in light of this modification. We also reproduce the DUNE sensitivity from Ref. [31], including this new decay channel. We do so by allowing the branching fraction of ϕ into SM particles, denoted as BR(ϕ → SM), to vary as a free parameter. In actuality, for a given (m ϕ , sin 2 ϑ), this BR is dependent on those parameters as well as m Z and g BL . The limits from these experiments, including the modifications, are shown in Fig. 5 for BR(ϕ → SM) = 100% (top), BR(ϕ → SM) = 10% (center), and BR(ϕ → SM) = 1% (bottom). Changes to the no-Z case (effectively, the left panel), are most apparent at higher m ϕ and larger sin 2 ϑ, e.g. the "lobe" in the CHARM exclusion closes up.
In Fig. 5, we also include several other results -in gray, we include constraints from rare B meson decays at LHCb [102,103]. We do not perform our own simulation of how these limits change for BR(ϕ → SM) = 1, however we expect that their behavior in the mass/angle range of interest will not vary much. In our simulations in Section 3.3, we do not allow sin 2 ϑ > 10 −6 for any m ϕ , in addition to requiring it to be consistent with all the experimental limits shown in Fig. 5. Also, in the top panel of Fig. 5, we include a limit from a reinterpretation of the LSND experiment's results [98] and future projections of the Fermilab SBN detectors, a combination of using the SBND and ICARUS experiments with the BNB and NuMI beams as proposed in Ref. [71]. Due to the complicated nature of simulating the LSND and SBN setups, we do not provide these limits/projections in the center/bottom panels, where they should in principle exist as well.
As with the light gauge boson providing a positive contribution to the muon anomalous magnetic moment, so does a light neutral scalar [75]. However, it requires sin 2 ϑ 1 mixing with Higgs for our light scalar ϕ to explain the (g − 2) µ anomaly [60], which is certainly excluded in Fig. 5 (well above the range we present). In general, a light scalar explanation of the (g − 2) µ anomaly must involve lepton-flavor-violating couplings to be consistent with the existing constraints; see, e.g. Refs. [106][107][108][109][110].
BBN can also provide constraints on the existence of ϕ if it is in thermal equilibrium and has yet to decay before the formation of light elements. If the lifetime τ ϕ 1 sec, the scalar ϕ will contribute a factor For each panel, we make a different assumption about the BR of ϕ into SM particles. The left panel assumes 100%, the center 10%, and the right 1%. We include constraints from CHARM (blue) [53,100], KOTO (green) [101], NA62 (red) [56], E949 (purple) [55], LHCb (gray) [102,103], projections from DUNE (black) in all panels, as well as constraints from LSND (orange) [98] and SBN projections from SBND and ICARUS (dashed blue) [71] in the left panel. The supernova constraints assuming the luminosities of 3 × 10 53 erg/sec and 5 × 10 53 erg/sec are shown respectively by the dot-dashed purple contours and purple shaded regions [104,105]. of 4/7 to the light degrees of freedom N eff [111,112]. Although the scalar contribution to N eff is smaller by a factor of three compared to Z , it is still larger than the allowed range of ∆N eff by Planck at the time of CMB formation, which is roughly 0.2 − 0.5 depending on the data sets adopted [87]. However, to implement the ∆N eff limits, ϕ has to be in thermal equilibrium with the SM particles at the BBN temperature of MeV scale [113]. The decaying of ϕ into e + e − and γγ are both very small, therefore there is almost no parameter space of m ϕ and sin θ to satisfy both the conditions above. The supernova limits on the scalar ϕ have been investigated in Refs. [104,[114][115][116][117][118][119]. In the supernova core, the scalar ϕ can be produced via the bremsstrahlung process N N → N N ϕ (with N being nucleons). Taking into account the partial cancellation between different Feynman diagrams [104], the observed neutrino luminosity L ν can be used to set limits on m ϕ and sin 2 ϑ. Setting L ν = 3 × 10 53 erg/sec and 5 × 10 53 erg/sec respectively, the regions inside the dot-dashed light-purple contours and shaded regions in Fig. 5 can be excluded. The production amplitudes for the diagrams with the scalar ϕ coupling to nucleons in Ref. [104] have been corrected by a factor of 1/2 [105]. One should note that in the supernova limits here we do not include the possible decay channel ϕ → Z Z . If this channel is kinematically allowed, the supernova limits on ϕ (and also Z ) could be potentially changed, but this requires a dedicated calculation beyond the scope of this paper, and is left for future work.

Improving Z prospects at DUNE
As discussed in Eq. (2.8), if both ϕ and Z exist and m ϕ > 2m Z , then the decay ϕ → Z Z may occur. In Section 3.2.2 we discussed how this decay channel modifies limits on the scalar boson in the (m ϕ , sin 2 ϑ) plane. In this subsection, we will discuss how the additional flux of Z from this channel can improve search capabilities in the (m Z , g BL ) plane.
Before analyzing this flux contribution, we briefly note that additional associated production of both ϕ and Z may occur from some meson decays. In the presence of the gauge coupling of ϕ to Z bosons, we can in principle have the associated production of Z and ϕ from meson decays, for instance π 0 → γ + Z + ϕ , K → π + Z + ϕ , ρ → Z + ϕ . (3.8) If the scalar ϕ is emitted from the Z boson line, the partial widths will be proportional to g 4 BL [32]. If the scalar ϕ is emitted from the quark fermion lines, the widths will be proportional to sin 2 ϑg 2 BL . Therefore for small sin ϑ and g BL , these associated production channels are highly suppressed compared to the channels above. Therefore, we do not consider these channels in this paper.
Since the main source of ϕ arises from decays of long-lived Kaons (compared to the prompt direct production of Z -neutral, short-lived meson decays and proton-proton bremsstrahlung), the Z coming from ϕ decays will originate all along the distance between the beam target and the end of the decay pipe. We take this into account by both allowing the Kaons to be long lived and the ϕ to have a lifetime according to the relevant parameters -m ϕ , sin 2 ϑ, g BL , and m Z . With these ingredients, we can determine the resulting Z flux, as well as the number of Z → e + e − and Z → µ + µ − events in our detector. Fig. 6 demonstrates how we determine this. For a given combination of (m Z , g BL ) -for example the BPs in Eqs. (2.9) and (2.10) -we determine the flux of Z at the DUNE Near Detector complex from "direct" production, i.e. the available neutral meson decays and proton-proton bremsstrahlung. That, combined with the Z lifetime (determined also by m Z and g BL ) and BR into visible final states, determines the number of signal events we observe. For these two combinations, we would expect ∼0.4 or ∼0.2 events in ten years of data collection at DUNE -not enough to warrant a significant discovery. However, if there exists ϕ with mass greater than 2m Z , then an additional Z flux can be generated. The total flux, and the number of The different colored contours represent increasing signal events (from darker to lighter colors) by an order of magnitude each, purple corresponding to N sig. < 1. In both panels, when m ϕ < 2m Z , the contribution is purely from direct production (meson decays and proton-proton bremsstrahlung), where ϕ → Z Z production increases the signal rate when m ϕ > 2m Z . In white solid lines we show a collection of current constraints from Section 3.2.2 on m ϕ and sin 2 ϑ, and in dashed white we show the expected DUNE sensitivity to search for the ϕ decays. In the vertical gray shaded regions, the production of ϕ from Kaon decay is kinematically forbidden.
signal events, would depend on the mass of ϕ and the mixing sin 2 ϑ, however, we cannot take sin 2 ϑ to be arbitrarily large, given existing constraints discussed in Section 3.2.2 and Fig. 5. In the left panel of Fig. 6, the colored contours display the number of expected signal events for the BP in Eq. (2.9) as we vary m ϕ and sin 2 ϑ, each color change representing an order-of-magnitude increase in N sig. from bottom to top (the purple regions correspond to N sig. < 1 and each subsequent color change goes to 1 ≤ N sig. ≤ 10 and so on). We scan over this parameter space, determining the maximum N sig. subject to the existing constraints in this space (solid white line) as well as the future DUNE sensitivity (dashed white line). For this combination of m Z and g BL , we obtain N sig. as large as ∼10 4 in ten years of data collection (this occurs if m ϕ ≈ 150 MeV and for sin 2 ϑ ≈ 3 × 10 −7 , where the colored contours in Fig. 6 left panel are lightest), which would be detectable in DUNE. In contrast, we show the result for the BP in Eq. (2.10) with a heavier Z in the right panel of Fig. 6, where the color changing denotes also the one-order-of-magnitude higher from bottom to top. As a result of the heavier Z mass, there is only a narrow window for m ϕ with the enhanced Z production rate at the DUNE, i.e. 2m Z < m ϕ < m K − m π , or equivalently 280 MeV < m ϕ 363 MeV. Beyond the mass threshold of m K − m π , the scalar ϕ can not be produced from the flavor-changing Kaon decays, which is shown as the vertical gray shaded regions in both panels of Fig. 6.
As a result of the extra production channel of ϕ → Z Z for Z , the improved DUNE sensitivities to the parameter space of m Z and g BL are shown in Fig. 7 as the dashed and dot-dashed black lines. As in Fig. 4, we use ten signal events as the required amount to be statistically significant. The dashed line (reaching lower g BL ) denotes the extended reach subject to the current constraints on m ϕ and sin 2 ϑ from CHARM, KOTO, NA62, and E949 discussed in Section 3.2.2. The dot-dashed line (extending below the solid black one but not as low as the dashed one) indicates the improvement subject to DUNE's direct ϕ decay search. In effect, if a signal of Z decay is observed consistent with the parameter space between the dot-dashed and dashed lines, then a similar ϕ → + − or ϕ → ππ signature should also be observable with the same data. Unlike in Fig. 4, here we do not show the prospects at future FASER and FASER2, because the presence of ϕ could in principle modify these projections and it requires a dedicated FASER simulation 7 , which is beyond the scope of our current work. The other labels are the same as in Fig. 4. In principle all other limits on Z such as those from the beam-dump experiments could also be affected by the scalar ϕ. However, for simplicity we do not include these corrections and assume naïvely that they are small.
The two BPs in Eqs. (2.9) and (2.10) are indicated by the stars in Fig. 7. The BP with m Z = 140 MeV is out of all current existing limits and even the Z prospects at DUNE in the pure gauge boson case in Section 3.1. However, it can be probed at DUNE in presence of the channel ϕ → Z Z . As shown in Fig. 7, the BP with m Z = 30 MeV is out of the reach of DUNE prospects in the pure Z case, but can be probed in presence of the scalar ϕ. It should be noted that this BP is nearly precluded by the supernova limits on Z [58]. However, as mentioned in Section 3.2.2, in presence of ϕ the supernova limits on Z could be very different. Furthermore, the supernova limits on Z could also change dramatically and thus be possibly avoided if Z couples to a dark matter particle [121]. In addition, the interplay of ϕ and Z would also be important for the evolution of Z in the early universe and the resultant BBN limits in Fig. 7.
It is clear in Fig. 7 that with the extra contribution from the scalar ϕ, DUNE can extend its search for Z to smaller g BL by a factor of 7 when subject to DUNE sensitivity to direct ϕ decays and by a factor of 45 when considering current ϕ constraints. It is striking that the scalar ϕ can improve the Z prospects at DUNE far beyond the current beam-dump limits for m Z 100 MeV, which would be otherwise very challenging for the pure Z boson case.

Conclusion
We have explored the possibility of experimentally searching for a class of BSM scenarios where the SM is extended by a minimal U (1) B−L gauge symmetry to explain neutrino masses and possibly provide a portal to DM. We focused on the light Z boson and light associated U (1) B−L -breaking scalar ϕ in the low gauge coupling regime. We discussed how the presence of Z and ϕ affects the search range at high-intensity facilities such as the near detector complex at DUNE. Focusing primarily on the sub-GeV mass region for Z and ϕ, the main results of this paper can be summarized as follows: • If the scalar ϕ decouples, the DUNE prospects of pure Z boson is presented in Fig. 4. Compared to the FASER sensitivities, the DUNE near detector complex can probe a smaller coupling g BL as a result of the lower center-of-mass proton-beam/target energy at DUNE (compared to the LHC).
• As far as the B − L breaking scalar ϕ is concerned, the presence of Z somewhat weakens the search range of its parameters at DUNE, compared to the case without the gauge interaction, as shown in Fig. 5.
• With the extra source of Z from scalar decay ϕ → Z Z , the number of Z signal events can be greatly increased for m Z 200 MeV (cf. Fig. 6), and the DUNE prospects of g BL can be improved by up to a factor of 45, as demonstrated in Fig. 7.

A Details of ϕ Decay
The transition amplitude G(s) appearing in Eq. (2.6) (with s = m 2 ϕ ) consists of three form factors. In the limit of isospin conservation, G(s) can be written in the form of [122] G(s) = 2 9 θ π (s) + 7 9 [Γ π (s) + ∆ π (s)] , where in chiral perturbation theory the three form factors are respectively .09, and the function where f π = 130 MeV is the pion decay constant, and κ = 1 − 4m 2 π /s. In Eq. (2.7) the loop functions are defined as [123] with τ X = m 2 ϕ /4m 2 X , and the function

B Additional Subdominant Contributions to the Gauge Boson Flux
For completeness, we list here possible additional production channels of Z in the U (1) B−L model at DUNE, although they are small with respect to the dominant channels in Section 3.1.1.
• Three-body charged pion/Kaon decays π + → + νZ and K + → + νZ with = e, µ. Because the incoming mesons have zero baryon or lepton number, this decay width should be identical to the leptonic Z case derived in Appendix B of Ref. [31] (see also Ref. [125]). This process can be calculated numerically -we find that the production of B − L Z gauge bosons from charged meson decays is suppressed relative to neutral meson decay/bremsstrahlung production, which we considered in the main text.
• Loop-level flavor-changing two-body decays K + → π + Z and K S → π 0 Z . Based on the low-energy chiral theory [126,127], the width is [128,129] Γ(K → πZ ) = g 2 BL m 2 Z W 2 (m 2 Z ) 2 12 π 5 m K λ(m K , m π , m Z ) where W 2 (q 2 ) incorporates both the polarization and ππ contributions [127], and Note that in the limit of m Z → 0, the width vanishes, as required by angular momentum conservation. The decay width in Eq. (B.2) applies to both K + → π + Z and K S → π 0 Z , which, however, are highly suppressed by the loop factor: BR(K S → π 0 Z ) 6 × 10 −6 × m Z 100 MeV 2 λ(m K 0 , m π 0 , m Z ) . (B.5) On the contrary, the process K L → πγ * is CP-violating, and as a result the decay K L → π 0 Z is highly suppressed [126,127]. There are also the flavor-violating decays B → πZ and B → KZ . However, the production rate of B mesons at the DUNE is highly suppressed [130], and we will not consider these B meson decay channels any more.

C Additional Z Boson Constraints
The current limits on the light Z boson of the U (1) B−L model are summarized in Fig. 4. Here we discuss a couple more limits from rare meson decays, which are not included in Ref. [37]. The most relevant ones are those from the loop-level flavor-changing K + decay (see Eq. (B.2) for the width) K + → π + Z , Z → νν . (C.1) The current most stringent limits are from NA62 [56] and E949 [55]. The NA62 experiment has obtained a 95% C.L. upper limit BR(K + → π + νν) < 2.44 × 10 −10 , which can be translated into the limit on m Z and g BL plane, as shown by the orange line in Fig. 4. The E949 limit on BR(K + → π + + X) (with X a long-lived particle) can reach up to 5.4 × 10 −11 , which is shown as the purple line in Fig. 4. The future NA62 data can improve the BR(K + → π + νν) constraint down to 2.35 × 10 −11 [61], and the resulting limit on g BL can be enhanced by a factor of 3, as indicated by the dashed orange line in Fig. 4. The gaps in the NA62 and E949 limits are due to the significant background due to the decays K + → π + π 0 with π 0 → νν. Note that because the width of K + → π + Z is proportional to Z mass (cf. Eq. (B.2)), both the NA62 and E949 limits get weaker when Z is lighter. There are in principle also the limits for a long-lived Z , with Z decaying outside the detectors. However, to have a long-lived Z the coupling g BL needs to be very small, which will in turn highly suppress the production of Z . It turns out that no parameter space of m Z and g BL is excluded due to the long-lived Z at NA62 and E949. The KOTO experiment has put an upper bound of BR(K L → πνν) < 3.0 × 10 −9 [101,131]. However, this limit can not be used to set constraints on the decay K S → π 0 Z .
There are more limits from the meson and baryon decays which can be used to constrain m Z and g BL , which are however much weaker, thus are not shown in Fig. 4: • The two-body meson decay K L → γZ with the subsequent decays Z → e + e − , µ + µ − contributes respectively to the decays K L → e + e − γ [132][133][134] and K L → µ + µ − γ [135]. However, as stated in Appendix B, the meson decay K L → γZ is highly suppressed by the corresponding SM BR(K L → γγ).
• The baryon decay ∆ → pZ with Z → e + e − contributes to ∆ → pe + e − , which is highly suppressed by the small SM BR(∆ → N γ), as mentioned in Appendix B.