Dark sector-photon interactions in proton-beam experiments

We consider electromagnetically neutral dark states that couple to the photon through higher dimensional effective operators, such as electric and magnetic dipole moment, anapole moment and charge radius operators. We investigate the possibility of probing the existence of such dark states, taking a Dirac fermion χ as an example, at several representative proton-beam experiments. As no positive signal has been reported, we obtain upper limits (or projected sensitivities) on the corresponding electromagnetic form factors for dark states lighter than several GeV. We demonstrate that while the current limits from proton-beam experiments are at most comparable with those from high-energy electron colliders, future experiments, such as DUNE and SHiP, will be able to improve the sensitivities to electric and magnetic dipole moment interactions, owing to their high intensity.


I. INTRODUCTION
The operation and development of high-intensity proton facilities are the backbone of the world-wide shortand long-baseline neutrino program. The collisions of high-energy proton beams on fixed targets deliver the neutrino fluxes that are registered in near [O(10-100 m)] and far [O(100-1000 km)] detectors through charged and neutral current interactions. In addition to mapping out the still elusive neutrino sector of the Standard Model (SM), the near detectors increasingly serve a second purpose: they become instruments to test new physics beyond SM. Dark sector particles with masses in the GeVrange and below can be produced and lead to observable signals in many previous, existing and upcoming neutrino experiments, such as LSND [1], MiniBooNE [2], COHER-ENT [3], DUNE [4], among others [5]. This dual purpose is further supported by dedicated experiments that aim to probe dark sector states, such as the proposed SHiP experiment or various beam-dump searches in the past; for an overview see [6] and references therein.
The interest in light dark states not only concerns DM detection, but also generally aims to test the presence of new sub-GeV particles in nature. Taking this broader point of view, in [41] we studied in-depth the possibility that χ particles-not necessarily the main component of DM-carry electromagnetic (EM) form factor interactions. In particular, we studied χ pair-production in electron beams on fixed targets at NA64 [42], LDMX [43], BDX [44], and mQ [9], in e + e − colliders BaBar [45] and Belle-II [46] and in proton-proton collisions at LHC [47]. In addition, flavor constraints from rare meson decays such as K ± → π ± + inv. [48][49][50][51][52] and precision observables such as (g − 2) of the muon [53][54][55][56] or constraints on the running of the fine-structure constant [57] were worked out. These studies of electromagnetic moment dark states in the MeV-GeV mass bracket were then further complemented by a detailed astrophysical study of stellar cooling constraints, once the mass drops below the MeV case [58]; see also related [59].
In this work we bring the above phenomenological studies to a closure by working out the current limits on and detection prospects of electromagnetically interacting χ particles utilizing high-intensity proton beams. Concretely, we set limits from LSND [1], MiniBooNE-DM [60], CHARM-II [61,62], and E613 [63] and forecast the sensitivity of SHiP [64] and DUNE [4,65]. In this process we take into account the most important χ pairproducing reactions from the Drell-Yan (DY) process, as well as from scalar-and vector-meson decays. With a detailed prediction of the energy-and angular-differential flux at hand, we subsequently derive the observable signals from χ-electron scattering and χ-nucleon deep inelastic scattering.
The paper is organized as follows: in Sec. II the electromagnetic form-factor interactions of χ are introduced. Sec. III contains a detailed discussion of dark state production in proton-beam experiments, and Sec. IV calcu-lates the generic signal-generation in the detectors. We list parameters of the considered experiments in Sec. V and present the derived results in the parameter space of dark state coupling vs. mass in Sec. VI. We conclude in Sec. VII and provide further details on our calculations in several appendices.

II. ELECTROMAGNETIC FORM FACTORS
We consider a neutral Dirac fermion χ as the dark state which interacts with the photon, A µ , or its field strength tensor, F µν . At mass dimension-5, the Lagrangian reads where µ χ and d χ are the MDM and EDM couplings, expressed in units of the Bohr magneton µ B ≡ e/(2m e ) below; m e is the electron mass and σ µν ≡ i[γ µ , γ ν ]/2. At mass dimension-6, the Lagrangian reads in which a χ and b χ are the AM and CR couplings. We take the values of all the couplings as real numbers. The possible UV completion for the effective interaction in Eqs. (1) and (2) could come from the compositeness of the dark states [66][67][68] or radiatively, from extra U(1) charged particles at high energy scales. For simplicity, we assume the effective operator approach holds for the beam energy of those experiments considered throughout the paper. The above interactions are then assembled in the matrix element of the dark current, where p i,f and q = p i − p f are both four-momenta. The vertex factor is The independence of momentum-transfer in the couplings follows from the assumption that the UV completion scale is much higher than the center-of-mass energies considered in this work.

III. PRODUCTION OF DARK STATES
At proton-beam experiments, dark states coupled to the photon can be produced via prompt processes (e.g. DY process or proton bremsstrahlung) and secondary processes (e.g. in meson decays or secondary collisions). In this section, we discuss these production processes and provide the calculations of dominant channels. Numerical results, taking the SHiP experiment as an example, are shown in Fig. 1. The relative importance of the individual contributions does not change significantly from experiment to experiment.

A. Drell-Yan Production
Dark states with effective couplings to the photon can be pair-produced directly through quark-antiquark annihilation. To correctly estimate the χ production from proton-proton collision, we utilize the event generator MadGraph 5 [69], to obtain the energy spectrum and angular distribution of dark states per collision, denoted as d 2N DY χ /(dE χ d cos θ χ ), as a function of χ energy E χ and the angle between their momentum and the beam axis, θ χ .
We then take the thick target limit, and calculate the total yield of dark states from the DY process via where the proton on target (POT) number is known for each experiment and A is the atomic mass number of the target; α 1 , α 2 are scaling-indices induced by scattering off a nucleus instead of a proton for the DY cross section, and the total scattering cross section, respectively. DY processes can be treated as incoherent and thus α 1 1. The value of α 2 , for inclusive proton-nucleus scattering, typically of the order O(0.8), depends on the exact target material, and only mildly affects the final results [70].

B. Meson Decay
Another important process is the secondary production of a χ-pair in the decays of scalar/vector mesons through an off-shell photon. Here we consider the scalar mesons π 0 , η and η , as well as vector mesons ρ, ω, φ and J/Ψ.
Typically, if the decays of meson into dark states are kinematically allowed, they tend to dominate the production rate. For example, [71] shows that the production of milli-charged particles from meson decay is several orders of magnitude larger than that from DY. Among them, the π 0 decay contribution is the most important. However, this picture changes when one considers higherdimensional operators. This is because the decay rate of light mesons into χ-pairs will receive additional suppression from their masses, as shown below.
For scalar mesons, the dominant decay channel producing dark states is a three-body decay with final states γχχ. By factorizing out the dark current part, we infer that where the subscript "sm" denotes "scalar meson". The branching ratios, Br(sm → γγ), are taken from the PDG [52]. It is worthwhile pointing out that in this step we neglect the mild q 2 -dependence induced by the meson  transition form factors F smγγ * (q 2 , 0). Such approximation is particularly justified for the lighter mesons: the photon virtuality is limited by kinematics, q 2 ≤ m 2 sm and corrections enter at the level of q 2 /m 2 ρ where m ρ is the ρ-meson mass; see e.g. [72][73][74] and Fig. 7 in App. A. To calculate Γ sm→γχχ and thus the ratio of the two channels, we follow our previous methodology in [41,58] and provide the corresponding expressions in App. A.
A vector meson, in turn, can decay into a χ pair directly. Thus we compute the branching ratio Br(vm → χχ), where the subscript "vm" denotes "vector meson". For two-body decays, one can separate the decay amplitude and phase space factors to obtain where the last two factors count the differences induced by the interaction type and the phase space, respectively. The expression of f (m 2 vm ) for each interaction type is given in App. A, and has previously been derived in [41,58]. In contrast to the (milli-)charged case (f e ), the function f χ heavily relies on the meson mass for higher-dimensional operators, and the χ production rate becomes enhanced for heavier meson decay; for more details see the appendix.
To calculate the χ production rate, the energy and angular distributions of the produced scalar/vector mesons are required. However, the latter are still poorly understood. One reasonable method is to estimate the normalized neutral meson distribution using charged meson distributions. Taking the neutral pion π 0 as an example, we follow [75] and stipulate that N π 0 ∼ (N π − + N π + )/2 and write the π 0 distribution as where E π 0,−,+ is the energy of the pion and θ π 0,−,+ is its respective emission angle relative to the beam axis. For charged meson distributions, we follow the literature and use the Burman-Smith parameterization [76] for sub-GeV kinetic energy proton beams such as LSND, use the Sanford-Wang distribution [77] for moderate beam energies (several GeV) such as MiniBooNE, and use the so-called BMPT distribution [78] for larger beam energies (from tens of GeV to hundreds of GeV), such as for DUNE and SHiP.
Besides the normalized meson distribution discussed above, we also require the total number of produced mesons in each experiment. For this, we use PYTHIA 8.2 [79,80] to simulate pp collisions, and list the average number of mesons produced per POT for each experiment in Tab. I. We assume that these meson production rates per POT remain the same for pN collisions; for the latter, current detailed simulations yield differing results, see, e.g., [81] for a recent discussion. 1 The meson multiplicities of our Tab. I lie within the range of their adopted values in previous works, e.g. [81,[83][84][85][86][87][88], and we estimate the uncertainties only affect the final bounds by a factor of 1.2 at most. Finally, we have also extracted the information on their momentum and angular distributions from PYTHIA 8.2, which is consistent with the fitted distributions mentioned above [88].
For the final distribution function of χ particles from meson decay in the lab frame we find where E * χ , θ * and φ * are defined in the meson rest frame and denote respectively the χ energy, as well as the polar and azimuthal angles of the χ momentum w.r.t. the direction of the boosted meson. In contrast, E χ and θ χ are defined in lab frame, and represent the energy of χ and the angle of the χ momentum w.r.t. the beam axis. At last, E m and θ m , the energy of the meson and the angle of the meson momentum w.r.t. the beam axis in lab frame, are functions of θ * , φ * , E * χ , E χ and θ χ . The dark state spectrum from each meson decay, dN m where the factor 2 accounts for the pair production of dark states and Br χ is the aforementioned Br(sm → γχχ) or Br(vm → χχ), depending on the spin of the meson. Their exact expressions are given in App. A. Note that to obtain Eq. (9) we have used the fact that the meson decay at rest is isotropic. In practice, we perform Monte Carlo simulations to numerically obtain the distribution function of χ from meson decay, instead of integrating Eq. (9) directly, as the latter is prohibitively time-consuming.

C. Other production mechanisms
Here we discuss additional channels of χ-pair production. Prominently, proton-nucleus bremsstrahlung contributes to the production of χ particles. The process can e.g. be estimated using the Fermi-Weizsäcker-Williams method [89][90][91], as has been done in [75,[92][93][94]. However, for the higher-dimensional interactions studied here, the production of χ-pairs through bremsstrahlung is generally dominated by the contribution of the vector meson resonance at s χχ m 2 ρ,ω [95]. Since we have already taken into account the resonant contribution through the vector meson decay processes above, we will not consider the proton bremsstrahlung any further; thereby we also avoid any double-counting.
Another source of χ-pair production is the capture of pions onto nuclei or protons via pπ − → nγ * → nχχ. This process will mostly result in low-energy χ-particles and is not considered further here. At last, secondary collisions, e.g. between secondary electron/photons and the target, should not appreciably contribute to the χ yield in our framework. We always neglect the latter contributions in this work.

IV. DETECTION OF DARK STATES
The dark states, produced in proton-nucleus collisions, travel relativistically through the shield into the downstream detector, leading to observable signals. In this work, we focus on their elastic scattering with electrons in the detector (LSND, MiniBooNE-DM, CHARM II, DUNE, SHiP) and hadronic shower signals caused by nuclear deep inelastic scattering (DIS) in E613.
For simplicity, we will approximate the detector-shapes as cylinders with a constant transverse cross-sectional area and a certain depth. Thus, the geometric acceptance of the dark states is determined by the target-detector distance and an effective size. For the nearly spherical detector in MiniBooNE-DM, we take the geometry into account in deriving the signal rate.

A. Scattering on electrons
When entering the detector, χ particles may scatter with electrons and cause detectable recoil signals. Following [41], the master formula to calculate the number of signal events reads where n e is the electron number density of the target, L det is the depth of the detector, E R is the electron recoil energy with respective experimental threshold and cutoff energies E min R and E max R , E χ is the initial χ energy in lab frame, and eff is the detection efficiency. The minimal energy of dark states E min χ can be expressed in terms of where m i is target mass, i.e. the electron mass in this case. The differential scattering cross section, dσ χe /dE R , is found in App. E of [41]. The spectrum of dark states that have entered the detector, dN χ /dE χ , is obtained by summing up all production processes in the previous section, and applying the detector geometric cut, The maximum opening angle θ max χ is obtained from the target-detector distance and the effective size of the detector. This is illustrated in Fig. 2 for the SHiP experiment (400 GeV proton) and which is obviously independent of the values of form factor couplings. As shown by the figures, only about 0.1-10 −5 of the total number of χ particles produced reach the detectors, and this strongly suppresses the number of final events at low energy experiments, such as at MiniBooNE-DM. 2 Moreover, such reduction becomes more severe for dark particles generated from heavy meson decay, and is largely insensitive to m χ for χ particles from DY processes. Besides, for higher-dimensional operators, a preference for more energetic χ particles can also be observed by comparing the left and right panels in Fig. 2 (also in Fig. 3). This is due to their different energy-dependence in the production rate, and will be further discussed in Sec. VI A.
Several experiments also make cuts on the electron recoil angle, θ R , in order to reduce backgrounds. From kinematics, the recoil energy E R can be expressed in terms of E χ and θ R as Hence, we take the cuts on θ R as a further requirement on the boundaries of E R , where θ min For the spherical detector in the MiniBooNE-DM experiment, we use an incoming angle-dependent depth for the detector, which reads where the radius of the detector R = 6.1 m, and the distance between the collision point and the detector center D = 490 m. Finally, we note in passing that while we focus on electron recoil signals, which in general provide better bounds, the expressions above are easily generalized to describe nucleon recoils.

B. Hadronic showers
The dark states may also cause hadronic showers, which is relevant for E613. Following [11] we consider the deep inelastic scattering of χ with nucleons as the energy deposition process, while neglecting any coherence effects since the typical momentum transfer is larger than the QCD confinement scale. It is worth pointing out that we do not consider the possibility of multiple scatterings in the detector, since the coupling between the χ particle and the photon is assumed to be weak; see Sec. IV C.
To derive the expected number of signal events, we first compute the differential cross section of χ-N deep inelastic scattering. The 4-momentum of χ before (after) scattering is denoted as p χ (p χ ). The momentum transfer carried by the intermediate photon is defined as q = p χ − p χ , which is spacelike. Following the DIS formalism for leptons, we introduce the Bjorken variable x ≡ Q 2 /(2m N ν), with m N m p being the nucleon mass, Q 2 ≡ −q 2 > 0 and ν being the energy transfer in the rest frame of the nucleons. The differential cross section is then written as where the dark current L µν can be written in terms of the vertex factors of Sec. II, with the factor 1/2 coming from average over initial state χ-spins. The hadronic tensor W µν may be expressed as in which with p N being the 4-momentum of the nucleon before the scattering. We adopt the results for the two structure functions as where e q is the charge of quarks in unit of electron charge. We sum over flavors of light quarks/antiquarks, q = u,ū, d,d, s,s, and use the values of parton distribution function xf q (x, Q 2 ) averaged over nucleons for each corresponding nucleus from [99]. 3 The expected number of signal events is given by where n N is the number density of nucleons in the detector. The integration boundaries for ν and Q 2 are derived from kinematics as where E cut is the experiment-specific threshold energy. The squared momentum transfer Q 2 lies in the range Finally, there is the general requirement x < 1.
C. Mean-free-path of dark states As already mentioned above, our calculations are based on the assumption that χ particles travel freely, both in the shield and in the detector. This may be validated by estimating the mean free path of χ, using transport cross section σ T χp of χ-proton scatterings, where n p is the proton number density. The transport cross section is used as it removes the influence of soft scatterings that would not attenuate the flux of dark particles.
To obtain an estimate, we use the elastic scattering processes for which the formulae can be found in App. E of [41]. Here we take the typical distance between the collision point and the detector to be 100 m and the dump/shield mass density to be 10 g/cm 3 . By requiring λ χ ≥ 100 m, one can see that these proton-beam experiments are sensitive to the EM form factor parameters for sub-GeV χ particles with E χ = 5 GeV. As parameters larger than these values above are already excluded by other probes, we may always assume that χ particles scatter at best once inside the entire experimental setup. 3 Such parameterization is numerically equivalent to the one of [11] in the limit of ν 2 Q 2 , which is the case in E613.

V. EXPERIMENTS
In this section, we briefly review the relevant details of each experiment under consideration. In order to derive the ensuing 90% C.L. limits, we require that the number of events generated by the dark states, where N obs is the number of actual observed events and N bkg is the expected number of background events. When making forecasts for future experiments, we assume N obs = N bkg . The standard criterion N sig ≤ 2.3 is adopted if no events were observed. For each experiment, the summary of relevant parameters can be found in Table. II.

A. LSND
At the Liquid Scintillator Neutrino Detector (LSND) experiment, a proton beam of 800 MeV kinetic energy was conducted onto water or a high-Z target such as copper [1]. The detector was located at a distance of 35 m from the beam dump, with an off-axis angle of 31 • , and an active volume comprised of an 8.3 m long cylinder with a diameter of 5.7 m, filled with 167 tonnes of mineral oil CH 2 [103].
Due to the low beam energy, we consider π 0 decay as the only χ production channel in LSND as other heavier mesons decay and DY channels are suppressed. As it is difficult to generate the total production rate of π 0 in PYTHIA 8.2 at such low energy, we instead estimate it via the ratio (σ pp→X+π 0 +2σ pp→X+2π 0 )/σ pp , which measurements put at a value of approximately 0.1 [104,105]. Under the assumption that this ratio remains unchanged for proton-nuclear scattering, we adopt the value 0.1π 0 /POT as our fiducial value in the calculation. This is close to the production rate of positively charged mesons in LSND, about 0.08π + /POT [106], as well as the value used in COHERENT experiment, 0.09π 0 /POT [107].
In the MDM case with m χ m π , the χ flux entering the detector is then approximately yielding the constraint µ χ ≤ 10 −5 µ B . This can be rescaled to compare with the LSND results [108], which estimates that the ν e flux entering the detector, Φ νe , is about 1.2 × 10 14 cm −2 , leading to a bound on ν e 's MDM at µ νe ≤ 10 −9 µ B [108]. One can see the equality is approximately satisfied, suggesting that our treatment of the detector works well. is the maximal angle between χ's momentum and the beam axis in order for χ to pass through the detector, ER is the recoil energy of the target, θR is the recoil angle of the target with respect to the χ momentum and eff is the detection efficiency of considered signal.

B. MiniBooNE-DM
The Booster Neutrino Experiment, MiniBooNE, operates at the Fermi National Accelerator [109]. The Booster delivers a proton beam with kinetic energy E beam = 8 GeV ( √ s ∼ 4.3 GeV) on a beryllium (A Be = 9) target. The center of the spherical on-axis detector is placed 490 m downstream from the beam dump with a diameter of 12.2 m filled with 818 tonnes of mineral oil C n H 2n+2 (n ∼ 20). In practice, we are more interested in the off-target mode of MiniBooNE, where the proton beam hits directly the steel beam dump, with an ensuing smaller high-energy neutrino background. This is referred to as MiniBooNE-DM, which has data with 1.86 × 10 20 POT [60]. By only focusing on electrons with extremely small recoil angles, the background was effectively reduced to zero in this off-target mode [60]. That is, we derive the 90% C.L. limits on the couplings of dark states to the photon by requiring N sig ≤ 2.3. It is well known that in the on-target mode with 1.3 × 10 21 POT, MiniBooNE reported a significant excess of electron-like events [2]. In addition, the background event of a single electron recoil is estimated to be about one hundred, after the same cuts as above [110]. Substituting these values into Eq. (22) in turn suggests that the on-target mode should lead to slightly weaker limits than those from MiniBooNE-DM, despite its larger POT number.
C. CHARM II CERN High energy AcceleRator Mixed field facility II (CHARM II) was a fixed-target experiment designed for a precision measurement of the weak angle. It utilized a 450 GeV proton beam on a Be target, and collected data with 2.5 × 10 19 POT during 1987-1991 [62]. The main detector is a 692 t glass calorimeter (SiO 2 , on average A 20.7 per nucleus), and has an active area of 3.7 × 3.7 m 2 , about 870 m away from the target along the beam axis [61]. In this study, we focus on the single electron recoil signals, as the detector has an almost 100% efficiency to record electromagnetic showers for recoil energy E R ∈ [3,24] GeV.
To estimate the number of background events, we take N obs = 5429 reported in [62], largely induced by electron scattering with energetic ν µ +ν µ particles. This estimation is conservative, as CHARM II was able to determine the value of the Weinberg angle with the uncertainty below several percents.

D. DUNE
The Deep Underground Neutrino Experiment (DUNE) is proposed to be performed at the Long-Baseline Neutrino Facility (LBNF), and can be used to probe light dark particles [4,65]. At DUNE, a graphite (A C = 12) target is hit by a proton beam with an initial energy E beam = 120 GeV. The near detector (75 t fiducial mass) will be placed 574 m downstream from the target. It is on-axis and a parallelepiped with a size 4 × 3 × 5 m 3 and we use 5 m as its effective depth [71]. The detector is filled with liquid Argon (LAr).
We take a 10-year run of the DUNE experiment, with a total POT of 1.1 × 10 22 . The observable signals we consider for DUNE are single electron events caused by χ-e scatterings. The detection efficiency is assumed to be eff = 0.5 for the LAr time projection chamber. Following [87,101], we require the cut on the electron recoil angle to satisfy E R θ 2 R ≤ 1 MeV, which significantly reduces the number of background events from charged-current ν e -n scattering; see Tab. II for details of the parameters.

E. SHiP
A fixed-target facility to Search for Hidden Particles (SHiP) is proposed at the CERN super proton synchrotron (SPS) accelerator [64]. At the SPS facility, a proton beam with E beam = 400 GeV ( √ s ∼ 27.4 GeV) is deployed to collide with the titanium-zirconium doped molybdenum target (A Mo = 95.95). An emulsion cloud chamber detector will be located 56.5 m downstream from the target, and it will be filled with layers of nuclear emulsion films. Following the latest SHiP report [111], the size of the detector (∼ 8 tonnes) is set to be 80 × 80 × 100 cm 3 . We assume a 100% detection efficiency for simplicity. 4 The detection process we consider for SHiP is also χ-e scatterings. With 2×10 20 POT after 5-years of operation the number of background events is estimated to be 846, which is dominated by ν e quasi-elastic scattering with a soft final state proton [111].

F. E613
E613 was a beam dump experiment at Fermilab, set up to study neutrino production, with a 400 GeV proton beam hitting a tungsten target [63]. The detector, 55.8 m away from the target, consisted of 200 tonnes lead plus liquid scintillator. Its size was 1.5 × 3 × 3 m 3 , with a mass density of about 10 g/cm 3 . In order to compare with the previous results [11,113], we only consider a circular region of the detector with a radius of 0.75 m along the beam axis. Moreover, for nucleon-recoil events in E613, the energy deposit needs to be larger than 20 GeV, in order to be recorded. We require the number of such events to be below 180 during its 1.8 × 10 17 POT run to obtain the constraints.
We assume a thick target so that each incident proton scatters once. This is different from the treatment by [11,113], which estimated the number of scatter events per POT following the scaling L T × n T σ pT (30) with L T (n T ) being the total length (the nucleon number density) of the target, and σ pT the scattering cross section between proton and target. For E613, where L T is much larger than the mean-free-path of a 400 GeV proton (a few cm in tungsten), Eq. (30) significantly overestimates the total number of produced χ particles. As a result, our limits are weaker than those derived in [113]. We revise the previous results in the next section.

G. Other experiments
There also exist many other proton-beam experiments which adopt similar setups to those we have studied above, such as COHERENT with a 1 GeV proton beam [3], JSNS 2 with a 3 GeV proton beam [114] and NOνA with a 120 GeV proton beam [115]. Nevertheless, these experiments are in general not expected to provide noticeably stronger (projected) bounds than those obtained above (see e.g. [86,107,[116][117][118]), and are thus not further studied in this work. 4 A unity efficiency was also used in [102,112].
A different new bound on light dark states comes from the NA62 experiment, which has recently improved the constraint on π 0 → γ + inv. by three orders of magnitude [119]. This puts upper bounds on the MDM/EDM interactions of our interest as for m χ m π /2. They are weaker than the bounds obtained above, and become even weaker for higherdimensional operators, i.e. the AM/CR interactions.
High-energy colliders become more important for χ particles heavier than pions. For instance, at LHC, the upgrade of the MoEDAL experiment will be equipped with three deep liquid scintillator layers [120]. In addition, there will be the milliQan detector [121,122] which will be composed of three stacks of plastic scintillators. Both experiments are designed to be sensitive to milli-charged dark particles, of which the scattering cross section with electron/nucleus is dramatically enhanced at low momentum-transfer. As suggested in [123,124], such experiments will constrain the EDM form factor of dark states, where there also exists an enhancementalthough milder-in low momentum-transfer χ-e (χ-N ) region of elastic scattering. Moreover, proposed future colliders, such as HL-LHC and ILC, will be able to further improve the experimental sensitivity on all the EM form factors studied here; see e.g. [125][126][127].

VI. RESULTS
In this section, we first compare the production efficiency of various production channels, and then summarize our bounds on the EM form factors of dark states.

A. Comparison of Production Channels
In contrast to dark state-photon interactions through milli-charge, higher-dimensional operators are considered in this work. Therefore, dimensional analysis demands an extra energy scale E to compensate for the presence of the dimensionful coupling (E for dimension-5 operators and E 2 for dimension-6 operators) in cross sections and branching ratios, in comparison to the dimension-4 case. This typically suppresses the yield of dark states.
For DY, the relevant energy scale is of the order of the pp collision energy, √ s. We can then infer that for dimension-5 (dimension-6) operators the resulting cross section will contain a dimensionless factor µ 2 χ s and d 2 χ s (a 2 χ s 2 and b 2 χ s 2 ). Thus, the cross sections involving dimension-5 and 6 operators are further suppressed relative to dimension-4 interactions for d −1 χ , µ −1 χ √ s and a −1 χ , b −1 χ s, which incidentally are also required for the treatment of Eqs. (1) and (2) as effective operators. As a result, the DY process gains in relevance relative to the meson decay in the production of χ particles, especially  The SN1987A bound is taken from [58].
for dimension-6 operators as for the latter, the relevant energy scale is roughly the meson mass.
In addition, because of the mass-scaling, the relative importance of decaying meson contributions is also modified. The branching ratios into χ-pairs from light mesons become suppressed. Therefore, we can see that although heavier mesons are produced at lower rates, as shown in Tab. I, the final yields of dark states from their decay are comparable to (dominate over) those from light mesons for dimension-5 (dimension-6) operators.
The χ production rate of each channel, after applying the geometric cut, is demonstrated in Fig. 1. One can see that due to the reasons above, the overall pattern in our χ production rate becomes very different from those of milli-charged particles (see e.g. [71]) and dark photons (see e.g. [75]), where light meson decay is the most important production channel unless it is kinematically suppressed. 5

B. Constraints
The 90% C.L. constraints on the EM form factors derived above are shown by the colored regions in Fig. 4 (MDM and EDM) and Fig. 5 (AM and CR), together with our previous constraints (gray regions) [41,58]. As explained above, the strengths of higher-dimensional interactions are energy-sensitive, and constraints derived from current proton-beam experiments, with √ s below several to tens of GeV, turn out not to be competitive with the constraint from LEP [41]. For dimension-5 operators, future experiments such as DUNE (10-year) and SHiP will improve the sensitivity by a factor of 2-3, and become stronger than LEP for m χ < 1 GeV due to their high intensity. It is worth pointing out that the astrophysical bound from SN1987A constrains the MeV-region below 10 −8 µ B [58], well below the current and projected experimental sensitivity. For dimension-6 operators, the production and detection rates of light dark states are even more sensitive to the center-of-mass energy, suggesting it is unlikely for low-energy experiments to play any role in the foreseeable future. In E613 the initial energy of χ needs to be above 20 GeV to trigger an observable signal, but such large E χ also enhances the χ-proton scattering, making it difficult for χ particles to travel through the shield unless a χ , b χ 10 −2 GeV −2 . 6 Thus, future high energy colliders have better potential to probe dimension-6 dark state interactions.
At last, due to the consideration given in Sec. V F, we also revise the E613 bound on milli-charged particles from [11], although it has been surpassed by bounds derived from later experiments [9,13,41,128,129]. Our derivation also improves w.r.t. a much earlier work [8], by adding the production through decays of scalar mesons and by imposing the BMPT distribution for mesons. As shown in Fig. 6, if only DY processes are taken into account, our bound is weaker than that from [11] by about a factor of 7. By adding contributions from vector meson decay, the bound becomes stronger, approximately in agreement with [8] (dashed lines). Our final exclusion limit, taking into account all these contributions, is shown as the pink shaded region in the figure.

VII. CONCLUSIONS
In this work we study the production and detection of neutral fermionic dark states χ that carry EM form factors in proton-beam experiments. We consider the production of χχ-pairs in the collision of high-intensity protons on nuclear targets through prompt Drell-Yan scattering and in secondary meson decays. The detectable signals considered are single electron recoil events at LSND and MiniBooNE-DM, CHARM II, as well as at the proposed DUNE and SHiP experiments, and hadronic showers caused by nuclear deep inelastic scattering at E613.
Owing to the higher dimensionality of the considered operators (dimension 5 and 6), the relative importance of production channels is biased towards processes with larger intrinsic energy. As a consequence, Drell-Yan production and production in heavy meson decays gain prominence when compared to the milli-charged and dark photon cases, for which pion decays dominate the dark state yield. . This corrects a previously derived limit from DY production (blue dotted) [11] and improves previous work utilizing DY + vector meson decay only (blue dashed) [8]. Other bounds shown are from the MiniBooNE [13] (solid gray), SLAC beam dump (dash-dotted gray) and mQ [9] (dashed gray) experiments. See [13,14] for the sensitivity reach of Babar, BESIII, and future experiments.
We compute in detail the energy and angular distribution of the produced dark state flux and set the strongest constraints on the existence of χ-particles with MDM and EDM interactions in the MeV-GeV mass bracket, excluding dimensionful coefficients µ χ , d χ 8 × 10 −6 µ B , corresponding to an effective scale Λ 5 < 0.4 TeV. For the dimension-6 AM and CR interactions, we find a χ , b χ 3 × 10 −3 GeV −2 are excluded, pointing towards a comparably lower effective scale of Λ 6 < 20 GeV. In the latter case, the constraint is superseded by LEP. Finally, as a by-product of our study, we also revise previously obtained proton-beam dump bounds on milli-charged particles.
With a strong connection to the neutrino program, proton beam experiments constitute an active and diverse field, with a number of new experiments proposed such as SHiP and DUNE. However, because the interactions considered here are higher-dimensional, we find that the prospects of significantly improving the direct sensitivity on EM form factor couplings rather hinges on the future of high-energy collider experiments and their ability to produce collisions with an ever increased centerof-mass energy.
Acknowledgments The authors are supported by the New Frontiers program of the Austrian Academy of Sciences. JLK is supported by the Austrian Science Fund FWF under the Doctoral Program W1252-N27 Particles and Interactions. We acknowledge the use of computer packages for algebraic calculations [130,131].
Appendix A: Decay rates of scalar mesons The decay rate of scalar mesons into a photon plus a χ-pair, Γ χ ≡ Γ sm→γχχ , is given by where Γ sm→γγ * is the decay rate with an off-shell photon, with F sm being the decay constant of the meson. Since we are allowed to neglect the momentum-dependence of the EM transition form factors of the scalar mesons (shown below), F sm will drop out in the branching ratio, Eq. (6). In Eq. (A1), the function f χ (s χχ ) is defined by the phase space integral of the squared amplitude of γ * → χχ The explicit expressions of f χ (s χχ ) were already obtained in [41], and are also listed below: The expression for f e (m 2 vm ), used for vector meson decay in Sec. III B, is then given by To infer the energy spectrum of χ from scalar meson decay, we also need to know the differential decay rate dΓ χ /dE * χ in the rest frame of the meson ( p 1 = 0). To this end, we first compute the amplitude of the process sm(p 1 ) → γ(p 2 ) + χ(p 3 ) +χ(p 4 ). and define the two Lorentz-invariant variables s 23 = (p 2 + p 3 ) 2 and s 42 = (p 4 + p 2 ) 2 , so that s χχ becomes s χχ = m 2 sm + 2m Here, s 42 is The allowed kinematic range of E * χ is m χ ≤ E * χ ≤ m sm /2. At last we arrive at the differential branching ratio via from which one can directly see that the meson decay constant, F sm , cancels out in the ratio of Γ χ and Γ sm→γγ . At the end of this section, we comment on the assumption of using constant transition form factors for scalar meson decay. Vector meson dominance suggests that the assumption holds well for m 2 sm m 2 ρ , which is the case for π 0 decays. For the heavier scalar mesons considered in this work, η/η , we have numerically evaluated the differential decay rate using the EM transition form factor. For the η meson [132], the results are given in Fig. 7, which shows that the shape of dBr χ /dE * χ is only affected mildly by the (kinematically limited) virtuality of the intermediate photon. In the total decay rate for m χ = 1 (100) MeV, Br χ would increase by a factor of 1.3 (1.7) in the case of MDM/EDM, and by a factor of 1.8 (1.9) in the case of AM/CR. Hence, neglecting the momentum-dependence of the transition form factors leads to slightly weaker bounds, and is hence conservative.
Appendix B: χ distribution from meson decay In this section, the derivation of Eq. (9) is provided. In general, the number of χ particles produced from a certain meson distribution is given by