LSND Constraints on the Higgs Portal

High-luminosity fixed target experiments provide impressive sensitivity to new light weakly coupled degrees of freedom. We revisit the minimal case of a scalar singlet $S$ coupled to the Standard Model through the Higgs portal, that decays visibly to leptons for scalar masses below the di-pion threshold. The dataset from the LSND experiment is found to impose the leading constraints within two mass windows between $m_S \sim 100$ and 350 MeV. In the process, we analyze a number of scalar production channels in the target, finding that proton bremsstrahlung provides the dominant channel at LSND beam energies.

From an effective field theory perspective, classifying the interactions of new neutral states with the Standard Model (SM) according to their dimensionality, there are only three relevant or marginal 'portal' operators that are not suppressed by a new energy scale. The Higgs, vector and neutrino portals therefore comprise the leading couplings of the SM to a hidden or dark sector. Motivated in part by the phenomenology of light dark matter (DM), much theoretical and experimental effort has recently been focussed on these portals [39,40].
In this paper, we will consider the minimal Higgs portal [42], the unique relevant operator that can couple the SM model to a dark sector, where S is a new scalar singlet, H is the SM Higgs doublet, and A is a dimensional portal coupling. Along with being one of the few renormalizable portal couplings to a dark sector, and a potential force mediator for thermal relic models of light dark matter, this interaction is of intrinsic interest as an extension of the SM Higgs sector. The strongest existing constraints on the Higgs portal, in the low mass range where Br(S → l + l − ) ∼ 1, arise from searches for leptonic decays at the CHARM fixed target experiment at CERN [43][44][45], and analysis of K + → π + S signatures at the Brookhaven E949 experiment, with S escaping the detector before decaying and thus being counted as missing energy in the search for K + → π + νν [44][45][46]. The latter constraint is the most stringent, except in an S mass range relatively close to m π where significant backgrounds limit the reach of E949. For higher S masses, a range of accelerator and meson decay constraints apply to the Higgs portal [6,44,45,[47][48][49][50][51][52].
In this paper, we revisit the limits on the Higgs portal in the low mass m S < 350 MeV region, by studying the sensitivity of the LSND experiment, that is known to provide important constraints on the dark photon [18,19]. In particular, by analyzing a range of production channels in the interaction of the 800 MeV proton beam with the target, particularly proton bremsstrahlung [51], we find that existing LSND analyses with final state electrons and muons already impose the leading constraint on the Higgs portal in two mass windows between 100 and 350 MeV. Our final results are presented in Fig. 1, where we show LSND exclusions compared to existing limits e.g. from CHARM and E949, and recent projections for sensitivity at the Fermilab SBN facility [53].
The rest of this paper is organized as follows. In the next Section we define the Higgs portal, and summarize some of the relevant couplings and decays rates. In Section 3, the production of light scalars at LSND is discussed in some detail, and in Section 4 we present the sensitivity reach due to light scalars decaying to electrons and muons in the detector. Section 5 contains our concluding remarks.

HIGGS PORTAL
We extend the SM by adding a scalar singlet S, for which the leading relevant or marginal couplings to the Higgs doublet H comprise the Higgs portal [42], After electroweak symmetry breaking, and rediagonalizing by shifting the physical Higgs field as h → h + θS to remove hS mixing, the induced linear couplings take the form (Av/m 2 h ) 2 versus dark scalar mass m S . Exclusions from other sources (in gray) including LHCb [47], E949 K → π + invisible [44][45][46], CHARM S → e + e − , µ + µ − [43][44][45] and SN 1987A [54], as well as the projections for the on-axis SBND (orange) and off-axis ICARUS (purple) experiments at Fermilab [53], are shown for comparison (see the text for further details).
where the mixing angle θ Av/m 2 h 1 for the parameters of interest in this paper.
Integrating out the electroweak-scale degrees of freedom induces further couplings of S to light hadronic states. For the sub-GeV mass range of interest here, the relevant interactions take the form, Well-known 1-loop triangle diagrams generate the effective diphoton coupling [55], where F γ (m S GeV) ∼ O(1) is a loop function [51,55,56]. The coupling to nucleons can in turn be obtained through the use of low energy theorems [51], In principle this coupling should be extended to a formfactor, but for the kinematic regime of interest in this paper, there is no significant impact from hadronic scalar resonances, and the assumption that g SN N is a constant will be sufficient.
In analyzing the fixed target detection signatures of S decays, we will also require the leptonic decay width of S, which is given by [6], We have Br(S → e + e − ) 1 for 2m e < m S < 2m µ , which is the dominant decay channel over much of the mass range of interest here, while Br(S → µ + µ − ) 1 for 2m µ < m S < 2m π . Just above the pion threshold, Br(S → µ + µ − ) 0.15 − 0.2 [49], which will also be relevant below.

LIGHT SCALAR PRODUCTION AT LSND
The LSND experiment comprises an 800 MeV proton beam impacting a thick target, that was either water or a high Z metal at various stages of the experimental program. Over its lifetime LSND accumulated one of the largest proton on target (POT) datasets of any fixed target experiment, with over 10 23 POT in total [57,58]. The relevance of LSND for Higgs portal phenomenology was briefly addressed in [48]. In this section, we will revisit the production rate of scalars for m S < m π from a variety of channels.
Before we examine specific production modes, it is useful to compare this case to the scenario with a dark photon A µ kinetically mixed with the photon via the interaction 2 F µν F µν . This induces a low energy coupling of A to the electromagnetic current with strength e , with the kinetic mixing parameter. The leading production mode at LSND for low mass dark photons is pseudoscalar meson decay, e.g. Br(π 0 → A γ) ∼ 2 . Thus, for sufficiently light dark photons, we can estimate the number of dark photons produced as N (π) A ∼ 2 N π . The large pion (and eta) production rate, combined with the large radiative branching for pseudoscalar mesons makes this channel by far the most efficient. For scalars coupled through the Higgs portal, the situation is somewhat different, as the only mesons with substantial scalar branching rates are kaons and B mesons, which are not kinematically accessible at LSND. We find instead that the dominant production mode in this case is proton bremsstrahlung, p + N → X + S, with X the inclusive hadronic final state and N S ∼ 0.5θ 2 g 2 SN N N π ∼ 10 −6 θ 2 N π , which is substantially lower than the dark photon production rate due to the reduced scalar coupling to hadrons. In the rest of this section, we will discuss this production mode in more detail, along with further production channels via ∆ decay and the Primakov process for comparison.

A. Proton bremsstrahlung
Scalars can be produced through the SN N vertex via the proton-proton bremsstrahlung process p+p → S +X, where we focus on pp scattering due to its resonantly enhanced rate, proceeding via the ∆ ++ intermediate state.
At LSND beam energies, the beam protons are only moderately relativistic, and thus we will utilize two different procedures for the calculation adapted respectively to either sub-relativistic or highly relativistic beams. Comparing the scalar production rate using both techniques at LSND will allow for an assessment of the precision of the rate calculation.
Splitting function:-We will first follow the approach of Altarelli-Parisi and formulate the bremsstrahlung calculation in quantum mechanical perturbation theory, as recently discussed in this context in [51]. Since the beam protons are not ultra-relativistic at LSND, this is an extension of the conventional Weizsacker-Williams (WW) approximation in which the beam protons are often considered in the infinite-momentum frame. Nonetheless, we find that the kinematic range at LSND will still allow us to approximate the required rate in terms of the proton-proton cross-section and a calculable sub-process [59][60][61]. In this formalism, all states are on-shell and while 3-momentum is conserved, energy is not automatically conserved at each vertex, but only after summing all contributions.
The relevant diagram for the process is shown in Fig. 2(a).
We denote the momentum of the incoming proton and emitted S in the target rest frame by p µ p =(E p , 0, p p ) and where p T is the S transverse momentum with respect to the beam, and z is the fraction of longitudinal momentum carried by S.
The second order contribution to the matrix element, generically of the form V f j V ji /(E f − E i ) for a perturbation V , has two possible time orderings in this case for the process p+p t →S+X exchanging the intermediate state p . The two amplitudes can be written as [51,60]  where the intermediate proton's 3-momentum is fixed by p p = p p − p S , while the energy is not automatically conserved at the pp S vertex. Denoting the energy denominators as ∆E emit =E p +E S −E p and ∆E absorb =E p −E S +E p , then under the condition we can neglect the matrix element M absorb . This can be interpreted as the dominant contribution coming from initial state radiation. We have verified that this condition is satisfied to a few percent for LSND kinematics. Imposing a second condition, it is possible to write the differential cross section of the process p+p t →S+X in the approximate form [61], where σ pp is the total proton-proton scattering cross section, which varies between ∼ 30−45mb over the relevant energy range (see Fig. 4) [62], with s =2m p (E p −E S +m p ) the center of mass energy. Denoting the momentum transfer as q µ = (E p −E S , p p − p S ), the differential splitting probability of the proton to emit a scalar P split S can be represented in the form, The integration range for p T and z is determined by the kinematic conditions (8,9), where we require the kinematic variable on the left of each inequality to be at most 10% of the right hand side. The conditions are satisfied for z ∈ [0, 0.5] and p T < 300 MeV at LSND. The resulting distribution of scalars is shown in Fig. 3, which we see reaches above E S ∼ 300 MeV.
One pion exchange:-We will now consider a complementary approach, modelling proton-proton scattering via one pion exchange, which is expected to provide the dominant hadronic (as opposed to electromagnetic) contribution to bremsstrahlung at sub-relativistic beam energies. Using L = g πN NN γ 5 τ · πN , with g 2 πN N /(4π) ≈ 13.5, we first verify that the tree-level one pion exchange contribution to pp elastic scattering does provide a relatively good fit, after accounting for the electromagnetic component, as shown in Fig. 4. We utilize a dipole form for the pion-nucleon form factor ∼ 1/(1 + Q 2 /m 2 A ) 2 where m A ∼ 1 GeV is the axial mass [64,65], and similarly the proton electromagnetic form factor F 1 (Q 2 ) ∼ 1/(1 + Q 2 /(0.71GeV) 2 ) 2 . The contribution from one pion exchange is significant in a narrow energy range, and it is known that additional processes, such as two pion exchange, become important for beam momenta above 600-700 MeV [66,67]. Retaining just the one pion exchange contribution will nonetheless be sufficient in our case, as we are interested in the ratio of two-to three-body final states, in which the overall normalization of the pp cross section drops out as for the splitting function calculation above. Note that above a beam momentum of about a GeV, the inelastic channel via the ∆-resonance contributes at a comparable level to elastic scattering, but is not accounted for in this approximation.
We now compute the rate for initial state radiation of S, pp → ppS via one pion exchange, according to Fig. 2(b). For the analysis below, we use the full tree-level calculation of the 2-body and 3-body final states. How- ever, we can gain some intuition in the limit where Mandelstam s m 2 S , where the cross section takes the form s + · · · exhibiting the Sudakov double logarithim. For the finitem S kinematics of interest here, there are no sizeable IR/collinear effects, and so we will not need to include the corresponding loop contribution that is relevant in the m S → 0 limit.
To compare with the splitting function calculation above, we define the differential splitting probability of the proton to emit a scalar via one pion exchange in the form, The plot in Fig. (5) compares the two different methods of calculating the splitting probability as a function of the scalar energy E S . Similar results hold for other choices of m S . We observe that the ratio is O(1), with one pion exchange providing a rate that is slightly larger than the relativistic splitting function for LSND beam energies. This comparison nonetheless provides confidence in the rate calculation at the O(1) level.
Scalar production rate:-Utilizing only the splitting function calculation (10) as a conservative approximation for the total rate, the total number of scalars N S produced through the bremsstrahlung channel can be estimated numerically, where we normalize the rate to the number of π + produced, N π , which is given at LSND energies by the Burman-Smith distribution [68]. For m S = 1 MeV, we obtain N S ∼ 0.5θ 2 g 2 SN N N π . This calculational approach should capture part of the primary production channel, but as is apparent from the discussion above, it is only anticipated to be accurate up to O(1) factors.

B. Other production channels
In this subsection, we comment on a number of additional sub-leading scalar production channels.
∆ decay:-At LSND beam energies, roughly half the total proton-proton scattering cross section involves an inelastic process with resonant production, e.g. of ∆ ++ , which subsequently decays to p + π + . Indeed, the resonant excitation of ∆ (and Σ) hadronic resonances is the primary channel for pion production at LSND. This is partially incorporated in the analysis of bremsstrahlung above, in that it contributes to the total cross section, but there are additional channels involving S radiation from final states which are more probematic to calculate. A tractable contribution of this type involves 3-body ∆ decay, ∆ → π + p + S [48], as shown in Fig. 2(c). Computing the 3-body decay rate, using a phenomenological pion-Delta-nucleon vertex at the low-energy given by L int = g π∆N∆ µ N ∂ µ π, and assuming that the 2-body decay of ∆'s saturates pion production inside the target, we can estimate the number of scalars from the 3-body decay via the following ratio, N S ∼ N π × Γ ∆→pπS Γ∆→pπ . Evaluating the 3-body phase space integral numerically for m S = 1 MeV, we find N S ∼ 0.04θ 2 g 2 SN N N π , which is consistent with the estimate in [48] and about an order of magnitude below the bremsstrahlung rate. Note that in the collinear limit, scalars are produced isotropically in the ∆ rest frame. We have transformed the energy-angle distribution to the lab frame, using a Monte Carlo simulation, in which the energy-angle distribution of ∆ baryons in the lab frame was reconstructed from the Burman-Smith parameterization of the pion distribution. As expected this distribution is almost isotropic, reflecting the fact that the ∆'s are produced almost at rest, in comparison to the more forwarded-peaked distribution from bremsstrahlung. This further suppresses the event rate in the detector.
Primakov conversion:-There are several additional decay channels that will contribute to S production, as discussed in [48], but none are estimated to be larger than the ∆-decay channel discussed above. We have also considered a different topology that utilizes the effective diphoton coupling (3), via which scalars can be produced via the Primakov conversion of photons γ+N →S+N in the presence of nuclei with atomic (number) mass Z(A), as shown in Fig. 2(d). The dominant source of photons is provided by π 0 decays in the target [69]. The neutral pion decay length is roughly 0.1 µm at the LSND beam energy, and thus π 0 → γγ decays can effectively be treated as a distribution of real photons in the target. Applying the analysis of [70] to scalar production, the total cross section can be written as σ S = dk γ dΩ γ f γ (k γ , Ω γ )×σ γN →N S (k γ ), where f γ (k γ , θ γ )dk γ dΩ γ is the photon energy and angular distribution with angles (θ γ , φ γ ) defined with respect to the beam direction. The two-body cross section γ + N → S + N incorporates a Helm form factor [71], which is exponentially suppressed for momentum transfer above 200 MeV once coherence is lost [72]. Taking the Burman-Smith model of the pion distributions as an input, the dependence of the photon distribution on the energy and angle with respect to the beam axis was determined using a Monte Carlo simulation, and numerically evaluating the integrals for m S =1 MeV, we estimate the number of scalars produced as N S ∼10 −4 θ 2 g 2 SN N N π where the factor of g 2 SN N ∼ 10 −6 has been inserted purely for comparison. This S-production process is forward-peaked, but is subleading at LSND.

SENSITIVITY AT LSND
In this section, we focus on the dominant proton bremsstrahlung production mode and combine the rate and distribution of the last section with the experimental geometry and detection probability, in order to determine the LSND constraints on the Higgs portal. The LSND detector was a shielded 5.7m diameter cylinder of length 8.3m filled with 167 tons of mineral oil, that was on average at an angle of 14 degrees to the beamline, and at a distance of 30m from the target. Charged particles, such as electrons and muons, were detected via a combination of Cerenkov and scintillation light.
Once produced, the probability that an S particle decays inside the detector is where L i (and L f ) denote the distances from production at which the scalar will enter (and exit) the detector, while τ is the lifetime and β the velocity. This probability therefore depends on the scalar's energy as well as its direction with respect to beam axis. Due to the low beam energy at LSND, we do not consider scattering or Compton absorption signatures inside the detector, since the decay reach dominates the scattering reach by several orders of magnitude [48,73].
To normalize the overall event rate at LSND, we have used N π 0 , the total number of neutral pions produced. In practice, the π 0 distribution is taken to be an average of the measured π + and π − production rates in protonnucleon collisions, which differ by O(1) factors. For the LSND beam energy, we use the parameterization of the production cross-section given by Burman and Smith [68], and denote the total cross section as σ BS π . With this normalization, the number of scalars produced via proton bremsstrahlung, that subsequently deposit their energy in the LSND detector, can be schematically represented as follows, where ϑ(E S ,θ S ) summarizes the experimental cut conditions and ε eff is the corresponding detection efficiency.
To determine the sensitivity to scalar decays to electrons, we use the analysis [57,58], in which ν e were detected via the inclusive charged-current reaction ν e + 12 C→e − +X. Following [9], we make the assumption, based on the primary use of the scintillation to Cerenkov light ratio, that the e + e − pairs would be registered as indistinguishable from single electrons. Therefore, we assume that the scalar's energy would have been measured as the energy of a single-electron in the energy range 60 MeV to 200 MeV with the e + e − pair detection efficiency as for a single electron, i.e. ε eff ∼ 0.1. A similiar analysis [74] uses an energy cut between 160 MeV and 600 MeV on muons produced through the reactions ν µ (ν µ )+ 12 C→µ −(+) +p(n)+X in order to identify muon neutrino-like beam excess events inside the detector. We can use this analysis to find the sensitivity to S decays to muon pairs, although the efficiency is harder to estimate in this case given that the µ + µ − pair will have a lower boost than the corresponding electron decay. We will assume these events are also reconstructed as singlemuon events with efficiency ∼ 0.1 similar to the electron case, but show the results with hatching to indicate that the detection assumptions are distinct. In this case, we also account for the reduced branching fraction to muons when m S > 2m π . In both analyses, the number of beamexcess events does not exceed 20, which we take as the limit for both electron and muon decay channels. Using the energy-angle distribution of scalars produced dominantly through proton bremsstrahlung, as outlined in Sec. 3, and considering the geometric acceptance of the LSND detector as well as kinematic cuts and detection efficiencies for the final state particles, we numerically determined the event yields at LSND. The resulting event number contours are shown in Fig. 6, while our final 20 event limit contour is shown in Fig. 1, which also summarizes the results in comparison to a number of existing constraints as detailed in the Figure caption. We see that the LSND sensitivity to electron decays provides the leading constraint in a small window in scalar mass from 130 to 180 MeV, while the sensitivity to muon decays provides the leading constraint from 2m µ up to 350 MeV.

CONCLUDING REMARKS
In this paper, we have revisited the existing limits on one of the three UV-complete portals from the SM to a dark sector, namely the Higgs portal coupling to a singlet scalar. This portal is of particular interest as one of the generic mediation channels for the interaction with dark matter. We have shown that existing data from LSND, when combined with the dominant low energy production mode through proton bremsstrahlung, already excludes additional regions of parameter space for m S between 100 and 350 MeV. Future analyses are possible, which can extend this reach further. For example, NA62 at CERN provides greater sensitivity to K + → π + νν, and so the exclusion from E949 can be extended [52]. Simi- three blue-shaded contour regions corresponding to 1 event (light), 10 events (medium) and 1000 events (dark). Solid shading indicates event rates from electron decays, while hatched shading indicates event rates from muon decays. Existing exclusions from other sources (in gray) include LHCb [47], E949 K → π + invisible [44][45][46], and CHARM S → e + e − , µ + µ − [43][44][45] analyses.
larly, KOTO provides sensitivity through the neutral decay channel K L → π 0 νν (see e.g. the recent discussion in [75]). The short baseline neutrino (SBN) program at Fermilab will also provide new sensitivity to the Higgs portal, as recently analyzed in [53], and we exhibit the projected sensitivity for SBND and ICARUS from that reference in Fig. 1.