Probing Leptophilic Dark Sectors at Electron Beam-Dump Facilities

Medium-energy electron beam-dump experiments provide an intense sources of secondary muons. These particles can be used to search for muon-coupling light dark scalars that may explain the $(g-2)_\mu$ anomaly. We applied this idea to SLAC E137 experiment deriving new exclusion limits and evaluated the expected sensitivity for the planned Jefferson Lab BDX experiment (in case of null result report). The calculation is based on numerical simulations that include a realistic description of secondary muons generation in the dump, dark scalar production, propagation, and decay, and, finally, the decay products (electrons, positrons, or photons) interaction with the detector. For both experiments, exclusion limits were extended to cover a broader area in the scalar-to-muon coupling vs. scalar mass parameter space. This study demonstrates that electron beam-dump experiments have an enhanced sensitivity to new physics in processes that are usually studied using proton beams.


INTRODUCTION
Fifteen years since it was firstly proposed [1], the discrepancy between the measured value of the anomalous magnetic moment of muon and its Standard Model (SM) prediction remains unexplained. The so-called "(g − 2) µ anomaly" has provided a strong motivation for light hidden particles searches, opening a window to new physics beyond the Standard Model (SM). Popular candidates such as dark photons and dark Higgs have been tightly scrutinized as a possible explanation for such anomaly. However, this hypothesis has been excluded by existing measurements (see e.g. [2,3]). Nevertheless, other solutions, such as a new light particle that dominantly couples to muons, remain viable and deserve attention. In this paper, we examine models with a new leptophilic dark scalar, tuned to explain the (g − 2) µ anomaly. We employ a simplified model framework, with the effective Lagrangian for the dark scalar field S written as: where S is a real scalar field and the coupling between S and the SM leptons, g , is restricted to be (Lepton-specific) g e : g µ : g τ = m e : m µ : m τ (Muon-specific) g µ = 0, g e = g τ = 0.
Such effective Lagrangian may originate from an effective gauge-invariant dimension 5 operator [4][5][6][7]: where H, L, E are the SM Higgs doublet, lepton doublet, and lepton singlet respectively. Λ is the new physics scale and c i are the Wilson coefficients for flavor i. Although the relative size of c i and the effective coupling g i are a priori undetermined, one may naturally expect that their values are regulated by the principle of Minimal Flavor Violation (MFV) [8] and therefore proportional to the Yukawa couplings y i . This leads to the coupling relations dubbed as "lepton-specific". On the other end, one may go beyond the MFV principle and take the alignment hypothesis [7], i.e., the scalar only couples to a single fermion of SM and both the scalar and the Yukawa interactions are simultaneously diagonalizable in a single basis. As emphasized in [7], although the UV origin of the alignment hypothesis remains mysterious, it could give rise to couplings of S predominantly to one flavor in a technically natural way that suppresses lepton flavor violations. Here, we choose S to be dominantly coupled to muons and we refer to this scenario as "muon-specific". Given that scalar-to-quark coupling is absent and its electron coupling is suppressed in lepton-specific model (or absent in muon-specific model), the effective way to probe the leptonphilic dark scalar is via accelerator experiments with muon beams [9]. Such experiments have been proposed at Fermilab (FNAL-µ [6], M 3 [10]), and CERN (NA64-µ [6,11]). Alternatively, proposed-proton beam experiments such as SeaQuest [12], NA62 [13], or SHiP [14] can also probe the dark scalars via a significant amount of secondary muons or photons. In this paper, we propose a new and alternative way to prove such a model, making use of secondary muons generated in electron beam-dump facilities. Such secondary muons have much larger energy and spatial spread in comparison to muon beams originated by protons. Nevertheless, they are abundantly produced and this search can be performed with a minimal upgrade of existing or proposed electron beam-dump experiments. The novel approach proposed in this paper exploits the fact that, when a high-energy electron impinges on a fixed target, a large variety of secondary particles is produced. This significantly extends the number of physics cases accessible to electron beam dump experiments.
The paper is organized as follows. Section II describes in general the muon scalar production on a fixed tar-2 get. Section III describes the main properties of electron beam-dump experiments searching for S, including: secondary muons production; scalar particle emission, propagation and decay; and, finally, decay products detection. Section IV describes the Montecarlo-based numerical procedure we developed to evaluate the sensitivity of electron beam-dump experiments. Section V presents results for E137 [15] and BDX [16] experiments and, finally, Sec. VI discusses limits obtained by complementary probes relevant for the parameter space explored in this paper. Appendixes A and B report technical details of the Montecarlo calculation.

II. S PRODUCTION BY MUONS AND SUBSEQUENT DECAY
The main process responsible for scalar emission by an impinging muon on a fixed target is the so-called "radiative" production, i.e., the reaction µ + N → µ + N + S, with N being a target nucleus. The corresponding Feynman diagrams are reported in Fig. 1. In case of large muon energy, compared to the muon and to the scalar mass, the corresponding cross-section can be computed at first order using the improved Weizsacker-Williams approximation [6]. The differential cross-section, integrated over all the emission angles, reads [17] dσ dx where x is the scalar energy to muon energy ratio E S /E µ , 1 is the muon boost factor, β S = 1 − m 2 S /E 2 S is the scalar boost factor, and α is the finestructure constant. The effective photon flux χ is where −t is the square of the four-momentum transferred from the initial to the final state nucleus, t min and t max are the kinematic limits on t, and G 2 (t) is the combined atomic and nuclear form factor (see App. A). The energy distribution of the emitted scalar follows directly from the previous expressions. In case of a light scalar (m S m µ ), the x distribution is concentrated in the low-x region, with a maximum value at x max 1.4 m S /m µ . For a heavy scalar, instead, the outgoing S takes more significant portion of the muon energy and the corresponding x distribution is peaked close to x peak 1 − m S /E µ . The kinematic of the produced S is strongly peaked in the forward direction, with a typical emission angle wherex is the typical x value discussed before.
Under our simplified dark scalar model (Eq. 1), a light scalar S only decays to pairs of photons or leptons once the decay channel is kinetically accessible. The total decay width is given by where F 1/2 (τ ) is the fermionic loop-function for on-shell scalar and photons. It is given by In this paper, we primarily focus on 2m e < m S < 2m µ case, resulting in a dominant decay into an e + e − pair for the lepton-specific model or to a γγ pair for the muon-specific case. This is motivated by the fact that the beam dump experiments here considered have maximum sensitivity to these channels. The corresponding partial widths are where in the second equality of Eq. (9) we used the identity g e /m e = g µ /m µ valid for the lepton-specific model. We note that, for same values of model parameters m S and g µ , Γ e + e − Γ γγ , resulting in longer-lived S particles for the muon-specific model than for the lepton-specific one.
Finally, we observe that, for the case m S > 2m µ , S decays always to a µ + µ − pair, with larger decay width This results in a reduced sensitivity for a typical beamdump experiment, since the large decay width required to cope with the experimental acceptance (see next Sec.) results to a g µ value, and hence to a production crosssection, too small to have an appreciable production yield.  L sh is the total length of target and shielding, while L dec is the length of the downstream decay region, preceding the detector. The three insets show schematically the production of a µ + µ − pair by a photon in the dump, the production of S by a secondary muon, and the S decay to an e + e − pair. For a "BDX-like" experiment, instead, the shielding extends to the detector position, and the decay region overlaps with the detector volume.

III. SEARCHING FOR S WITH SECONDARY MUONS AT ELECTRON BEAM DUMPS
The typical setup of an electron beam dump experiment searching for muonphilic scalar particles is shown in Fig. 2. Scalar particles are produced by the primary electron beam impinging on the fixed target through a tertiary process involving secondary muons. These propagate in the target and in the surrounding materials, undergoing energy loss and multiple scattering, and radiatively emit the scalar particles. The resulting S kinematic distribution, including the production vertex, is thus given by the convolution of the muon distribution, altered by energy loss and multiple scattering, and the differential production cross-section.
After being produced, S particles propagate straight, until decaying to an e + e − or γγ pair. These particles are measured by a detector placed behind the beam-dump. A sizable amount of shielding material is placed between the dump and the detector to range out all other particles produced by the primary beam. Two detection setups are hence possible. In the first case (E137-like), the detector is sensitive to S decay products (e + , e − , or γ) produced in a free decay region downstream of the shielding. In the second scenario (BDX-like), the shielding extends up to the detector location, and the decay region overlaps with the detector volume. Production of high-energy muons by an impinging e − on a heavy thick target predominantly happens via two different mechanisms, the radiative emission of a µ + µ − pair and the decay-in-flight of photo-produced pions and kaons. The contribution of muons produced by the decay-at-rest of π's and K's can be neglected due to the very low emission energy, typically lower than the scalar production and detection thresholds.
Muons with energy of the order of the beam energy are mainly produced by pair emission, predominately happening via a two-step process [18]. First, an e − radiate a bremsstrahlung photon in a nucleuos field, then the photon produces the µ + µ − pair in a close-by nucleus. The direct production process through a virtual photon exchange, e − N → e − N γ * → e − N µ + µ − , is instead negligible, due to the much lower flux of virtual photons. The kinematics of produced muons is strongly peaked in the forward direction, with the typical magnitude of transverse momentum being p T,µ m µ . Being the pair-production on nuclei a coherent process, the yield of muons produced by the neutral component of the electro-magnetic shower through this reaction is almost independent on the target material (see [19]. Eq. 9).
In the low energy range, instead, the dominant production mechanism is the decay-in-flight of photo-produced π's and K's. The kinematics of emitted muons is isotropic in the hadron rest frame, and shows a forward peak in the laboratory frame whose width depends on the boost factor of the decaying parents. Since π's and K's photo-production on nuclei scales roughly as the atomic number, the corresponding muon yield shows a significant 1/Z material dependence, Z being the atomic number of the target. Figure 3 shows the kinetic energy distribution of muons produced by an 11 GeV e − beam impinging on a thick aluminum and water target, comparing the full yield to the contribution of photo-produced muons only. This result was obtained from a FLUKA-based [20,21] Montecarlo simulation, where we implemented the specific target geometry and composition foreseen in the BDX experiment (see Sec. V B). The overall muon yield per electron on target (EOT), in the energy range 0.5−11 GeV (2−11 GeV) is 7.65 · 10 −5 muons/EOT (0.91 · 10 −5 muons/EOT). For a primary e − beam current of about 10 µA, this corresponds to a muon rate of O(10 9 − 10 10 ) muons/second, depending on the considered energy range. This large flux demonstrates the potential of secondary muons at a multi-GeV electron beam dump facility, in particular when compared to the typical intensity of existing muon beams with similar energies (the Fermilab accelerator complex, for example, can deliver a muon beam of about 10 7 muons/second [22]). Figure 4 shows the kinematic distribution of produced muons, in terms of the emission angle and vertex longitudinal coordinate z (for the latter, z = 0 corresponds to the target front face.). Higher-energy muons are mostly (Colors online). The differential muons yield per EOT for 11 GeV e − beam a impinging on a thick aluminum and water target, as a function of the muon kinetic energy. The continuous black curve refers to all produced muons, while the dashed red curve refers only to pairproduced muons. The peak in the full distribution at E = 152 MeV is due to the kaon decay-at-rest process, K → µνµ. produced in the forward direction, in the first target radiation lengths, while for lower energies the angular distribution is wider, and particles are produced over the full target length.

B. S production
Muons produced in the thick target penetrate deeply in it and surrounding materials, losing energy mainly through ionization [23]. While traveling, they may produce a scalar particle through radiative emission on a nucleus.
In the simplest case where the target sizes are larger than the average muons range, and considering an uniform material, the total S production yield is given by where A and ρ, are, respectively, the target material atomic mass and mass density, E 0 is the primary electron-beam energy, N A is Avogadro's number, σ(E µ ) is the energy-dependent S production cross-section, and E min m S is the minimal muon kinetic energy required to produce a scalar with mass m S through radiative emission on a nucleus. Finally, T µ (E µ ) is the muons differen- tial track-length distribution in the target, defined as the integral over the target volume of the differential muons fluence Φ µ (E µ ), corresponding to the density of particle tracks in the volume [24]. Intuitively, the quantity T µ (E µ )dE µ represents the total path length in the thick target of muons with kinetic energy in the interval between E µ and E µ + dE µ . At the first order, neglecting multiple scattering and considering a constant stopping where n µ (E) is the differential yield of muons in the target. From previous considerations about muon production in the target, and given the shape of n µ (E), it follows that T (E µ ) depends on the target material as (ρZ) −1 in the low E µ region, and as ρ −1 log Z at higher muon energy. For a scalar in the mass range here considered most of the contribution to the integral in Eq. 12 comes from the low E µ region (E µ 1 GeV). Therefore, considering the Z 2 dependence of σ(E µ ) [6], it follows that the total scalar yield depends weekly on target material.
The kinematic distribution of scalar particles, including the production vertex, replicates the distribution of secondary muons in the beam dump, convolved with the differential production cross-section and distorted by the material-dependent mean free path for the µN → µN S process. In particular, the energy distribution of scalar particles (integrated over the full angular range) reads with dσ dE S given by Eq. 4. From the previous discussion about the shape of dσ dE S , it follows that for large m S val- . For lower mass values, instead, the broadening induced by dσ dE S is larger. Finally, we note that the shape of dN S dE S depends significantly on the target material. For example, the most energetic part of the spectrum, for E S E 0 , is due to the high-energy part of T (E µ ), and thus exhibits a Z log Z dependence.
Previous considerations were derived for a simplified geometry of a large and homogeneous target. In a more realistic case, the target may be non-homogeneous, and its size may be smaller than the range of muons, requiring a torough description of surrounding materials. The evaluation of the total S yield and of the corresponding kinematics thus requires to compute separately T (E µ ) in each geometry element, by tracking muons in the full experimental setup. For a realistic and accurate evaluation of real experimental setups, we performed this calculation numerically through an ad-hoc Montecarlo simulation, as described in Sec. IV.

C. Signal yield in the detector
After being produced, S propagates straight with a differential decay probability per unit path given by where is the energy-dependent S decay length. Particles from the S decay are emitted on a cone with typical aperture with respect to the S direction. The total signal yield is obtained combining the S angular and vertex distribution with the decay kinematics, and projecting the result over the geometrical acceptance of the detector. The latter can be roughly determined as the product of a longitudinal factor ε L depending on the shielding L sh and decay region L dec length and a transverse factor ε T related to the detector transverse area A T (see Fig. 2). For a single scalar particle produced with momentum p S and longitudinal vertex coordinate z S , the longitudinal factor reads The longitudinal acceptance ε L is obtained convolving P L with the scalar particles kinematic distribution. The transverse factor, instead, reads with θ rms S being the width of the scalar angular distribution, θ D the typical opening angle between the e + e − or γγ pair from S decay, z S the average S production vertex, and L D the average distance between the S decay point and the detector -in case of a BDX-like setup, with S decaying within the detector volume, L D =0.

IV. REACH EVALUATION
To evaluate the BDX reach and to determine the exclusion limit set by E137 null result for the dark scalar models, one has to determine, for a given combination of the model parameters (m S and g µ ), the expected number of signal events within the detector acceptance. We performed the calculation of the signal yield numerically, decoupling the evaluation of S particles production from the subsequent propagation and decay as follows.
First, we pre-computed the yield of muons produced by the primary electron beam through an optimized simulation using FLUKA 2011.2x.3. In the simulation, we implemented the description of each experimental setup (geometry and materials), including the thick target, the surrounding materials, and the subsequent shielding. We applied different biasing techniques, namely cross-section enhancement and leading-particle biasing for electromagnetic showers, to enhance forward muon production. The obtained result is a list of produced muons momenta, vertexes and statistical weights, together with the total yield per EOT.
Each muon was then individually tracked through a Geant4-based simulation [25], implementing the same geometry as implemented in FLUKA. We modified the Geant4 4.10.02 source code to include the new scalar particle S and the radiative S emission process (more details are given in App. B). When a scalar particle is produced, the corresponding four momentum and production vertex are saved to an output file, and the particle is not further tracked. The output of this second step is a list of produced S particles, together with the total normalization.
This result was finally used as input for a custom Montecarlo code that handles the S propagation and subsequent decay to an e + e − or γγ pair and computes the experimental acceptance of a detector placed downstream of the dump, ε det , including geometrical and detection threshold effects.
The overall normalization and the signal yield dependence on model parameters were handled as follows. The first computation step only considers S production, and the corresponding result depends on model parameters through the total number of produced S particles per EOT, N 0 S ∝ g 2 µ . The corresponding kinematics, instead, does not depend on g µ . Therefore, only the mass of the scalar, selected at the beginning of the simulation, was left as free parameter while the coupling g µ was fixed at a conventional value of g 0 µ ≡ 3.87 · 10 −4 . The detection acceptance computed in the second step, instead, depends critically on the g µ value, affecting the S decay length. Therefore, the evaluation of ε det has been performed as a function of this parameter. For a given combination of m S and g µ , the total signal yield in the detector per EOT thus reads We found this multi-step method more effective than performing a single Geant4 simulation, including both the muon production by the developing electromagnetic shower and the propagation and decay of S, since it saves a significant computation time in the evaluation of the sensitivity of a beam-dump experiment as a function of the S mass and coupling. Indeed, for a given experimental setup, the muon yield estimation is performed only once. Then, for a given mass value, the S production is evaluated, and only the computation of detector acceptance is repeated for different values of the scalar coupling.

V. RESULTS
In this section, we present the exclusion limit in the m S vs g µ parameters space obtained from the E137 experiment, and the expected sensitivity for the and BDX experiment, for both the lepton-specific and the muonspecific models. Results are summarized in Fig. 7, reporting also published limits from other experiments (see Figure caption for further details).

A. E137 exclusion limits
The E137 experiment searched for long-lived neutral objects produced in the electromagnetic shower initiated by 20 GeV electrons in the SLAC beam dump. Since muons penetrate deeply in the thick target and in the surrounding materials, a precise description of the experimental setup geometry and materials is required. The E137 beam dump was a 720 cm long aluminum tank completely filled with water [27,28], with an external diameter of 152 cm. The thickness of the aluminum vessel was 0.9525 cm, a part from the hemispherical front window, with a reduced thickness of 0.475 cm. The vessel radius was 75 cm. The underground hall hosting the dump was made by concrete walls, with thickness of about 90 cm. The distance between the dump end and the front wall was 160 cm. To protect the concrete wall from low-energy radiation escaping the dump, a 56 cm long iron block was added. The transverse size of the iron block, determined from the original engineering drawings, was comparable to that of the beam dump. Scalar particles produced in the thick target would have to penetrate 179 m of earth shielding and decay in the 204 m region downstream of the shield. We assumed a uniform dirt density of ρ dirt = 1.7 g/cm 3 , a value that seems to us reasonable considering what reported in Ref. [29] and given the average depth of soil of about 10 m for the terrain following the underground hall [15]. To evaluate the effect of this parameter on the final result, we repeated the E137 reach computation changing it by ±10%, obtaining a negligible variation of order ∼ 1%.
The E137 detector was an 8-radiation length electromagnetic calorimeter made by a sandwich of plastic scintillator paddles and iron (or aluminum) converters. Multi-wire proportional chambers provided an accurate angular resolution, essential to keep the cosmic background to a negligible level. A total charge of ∼ 30 C was dumped during the live-time of the experiment in two slightly different experimental setups: in the first run (accumulated charge 10 Coulomb), the detector had a transverse size of 2×3 m 2 , while in the second run this was 3×3 m 2 .
The original data analysis searched for axion-like particles decaying in e + e − pairs, requiring a deposited energy in the calorimeter larger than 1 GeV with a track pointing to the beam dump production vertex. The absence of any signal provided stringent limits on axions or photinos.
To derive the E137 exclusion limits on dark scalar production, we used the Montecarlo-based numerical approach described above. The experimental acceptance was evaluated separately for the two E137 runs and combined with proper weights to account for the different accumulated charges. In the calculation, we employed the same selection cuts used in the original analysis, requiring that the energy of the impinging particle is larger than 1 GeV and that the corresponding angle of impact on the detector, measured with respect to the primary beam axis, is smaller than 30 mrad. We found that both particles from S decay hit the detector in a nonnegligible fraction of events. In these cases, we applied previous selection cuts respectively considering the sum of the two energies to be greater than 1 GeV and the energy-averaged impinging angle to be less than 30 mrad. Based on the null observation reported by E137, we derived the exclusion contour considering a 95% C.L. upper limit of 3 events.

B. BDX sensitivity
BDX is a planned electron-beam dump experiment at JLab that will improve the E137 sensitivity by order of magnitudes by using the high intensity 11 GeV CEBAF beam [30], running for ∼1 year with currents up to 60 µA, thus collecting 10 22 EOT. BDX will make use of the Hall-A beam dump [31], with the detector placed 20 m downstream in a new, dedicated experimental Hall (see also Fig. 6). In this configuration, the experiment will produce ∼ 10 17 GeV-energy muons in the dump, a total yield much larger than what can be obtained in similar fixed-target efforts employing a primary muon beam [6].
The BDX detector consists of a 50 × 40 × 300 cm 3 electromagnetic calorimeter made by 800 CsI(Tl) scintillat- ing crystals, surrounded by two active veto layers made by plastic scintillator for cosmic backgrounds rejection. An iron and concrete shielding layer will be installed between the existing Hall-A beam dump vault and the detector hall to range out all other particles produced in the target (except neutrinos).
The primary physics goal of BDX is the search for light dark matter particles, produced by the primary electron beam in the target, by measuring their scattering on atomic electrons in the detector, resulting in high-energy electromagnetic shower [32]. In particular, the experiment was designed and optimized considering the physics case of an invisibly-decaying dark photon [33]. In case of a negative result, BDX is expected to improve current exclusion limits by about two order of magnitudes. The experiment is also sensitive to long-lived particles decaying visibly to an e + e − or γγ pair, if the decay occurs within the detector volume. In this work, we derived the expected BDX sensitivity to a scalar particle considering, conservatively, an expected background contribution of 5 events, corresponding to a 95% CL sensitivity of ∼ 6 signal events.

A. Anomalous magnetic moments
One of the most sensitive bounds on light scalars for lepton-specific or muon-specific models is from the measurement of muon anomalous magnetic moment a µ ≡ (g − 2) µ /2. Current measured value of a µ is larger than the SM prediction by [34] ∆a µ = a µ (EXP) − a µ (SM) = (268 ± 63 ± 43) × 10 −11 , (21) where the first and the second error bars indicate the experimental and theoretical systematical uncertainties, respectively. Introducing a new scalar that couples to muons can uplift the predicted value of a µ by and makes it agree with the current measured value. In the g µ -m S parameter space, we show the favored regions consistent with the measured value of a µ as the shaded green region with label a µ (±2σ). We also report the excluded regions that yield a 5σ-discrepancy with respect to the measured value as the shaded gray region (joint with other exclusions) with label a µ (5σ). On-going experimental and theoretical developments may reduce the experimental and theoretical systematical uncertainties by a factor of 4 and 2 respectively in near future [22,[36][37][38]. We indicated the favored and excluded regions for future a µ measurements as the hashed green and gray regions.

B. Cooling of SN 1987A
A core-collapsed supernova (SN) behaves like a protoneutron star and cools mainly through neutrino diffusion under the standard SN model. The measured SN 1987A neutrino burst flux agrees with SN model predictions [41] and hence can be used to constrain light dark scalars, which may provide an extra cooling channel. For the extra cooling to be effective, dark scalars should be abundantly produced in the supernova and not decay or trapped by nucleons and leptons along their exiting path. Reference [40] has considered an axion-like particle (ALP) with photon coupling −(g γγ /4)aFF . The dominant production and trapping mechanisms are through a Primakoff process γN → N S and an inverse-Primakoff process SN → N γ, respectively. Given the production and trapping mechanisms for ALPs are identical to those of muon-specific dark scalars, we ignored the difference in the CP quantum number and directly translated the resulting SN bounds on g γγ from [40] to those on g µ for muon-specific dark scalar through the relation The translated bound is shown as the gray shaded region with label "SN 1987A" in the right panel of Fig. 7 [42]. Note that the exclusion region is cut off at m S = 2m µ : a heavier scalar would be always trapped in the supernova core, without contributing to the supernova cooling [7]. A similar analysis can be applied to the lepton-specific  [34]. The gray regions are combined exclusions based on existing experiments or observations. For the lepton-specific model, the gray region includes current aµ measurement with 5σ exclusions, BaBar limits based on e + e− → µ + µ − S searches [35], and constraints from electron bremsstrahlung at Orsay and E137 [5]. For the muon-specific model, the gray region includes the same exclusions from aµ and BaBar as the ones for the lepton-specific model. Besides it also includes constraints from SN 1987A (see Sec. VI B for more details). The hashed regions in both panel indicate the projected sensitivity of future aµ measurements, assuming the current central value stays the same while the experimental and theoretical uncertainties will be improved by a factor of 4 and 2, respectively [22,[36][37][38]. The dotted contour shows expected sensitivity if BaBar or Belle could search for e + e− → τ + τ − S [5] or COMPASS could search for µN → µN S [39] with existing data. Note that we also indicate constraints from the secondary photons at E137 [40] as a dotted line given the uncertainties in the analysis as we discussed in Sec. VI C. Projected sensitivities for leptophilic dark scalars produced by future muon-beam experiments or secondary muons/photons at future proton-beam experiments can be found in [6,12].
dark scalars except that the electron-dark scalar interaction should be also included for the trapping. Nevertheless, the resulting bound is weaker than the constraints of other electron beam-dump experiments for the parameter space of interest. Hence it is not explicitly labeled in the joint exclusion region in the left panel of Fig. 7.

C. Other fixed-target and beam dump experiments
In case of the lepton-specific model, a constraint to g µ can be derived considering the emission of scalar particles by e + and e − in the beam-induced electromagnetic shower. Even if the corresponding radiative cross section is suppressed by the factor (m e /m µ ) 2 ∼ 2 · 10 −5 with respect to the muon case, the much larger track-length of e + and e − compensates for this, resulting in a non-negligible yield. In this case, since that the dump radiation length is much smaller than the target-detector distance, it is possible to assume that the scalar production happens always at the beginning of the target. This allows to exclusion limit for scalars from the sensitivity to a visibly-decaying dark photon [43,44], by following the procedure depicted in [6]. In the left plot of Fig. 7, results from E137 and BDX are reported.
In case of muon-specific model, there may be a significant contribution due to Primakoff process γN → N S by bremmstrahlung photons in the beam dump. Although g µ contribution to the g γγ is loop-suppressed, a large flux of photons produced by the primary beam may still yield a large number of signal events and consequently a competitive bound on g µ . In particular, in Ref. [40] the original E137 bound on the mass and coupling of ALPs was cast to a limit on g γγ . In the calculation, authors used the photons track-length distribution reported in the original E137 paper [15], and evaluated the experimental acceptance only including the angular effects due to the opening angle between the ALP decay products. We note the following critical issues in this procedure. First, it ignores the contribution to the detector geometrical acceptance from the finite production angle of ALPs in the beam dump. Secondly, we found the photons track length distribution reported in Ref. [15] to be inconsistent, in particular at low energy, with the result obtained from a FLUKA-based calculation that we performed (we cross-checked our result with the analytic predictions from showering theory [45,46], finding an excellent agreement). We still report the result obtained translating the aforementioned limit on g γγ to a limit on g µ , following the procedure depicted in Sec. VI B. We show the result as the grey-dotted line with label "E137(sec. γ)" on the right panel of Fig. 7, leaving a more detailed study for the future. Note that we explicitly checked the momentumsquared of the t−channel photon from MC simulations and found it to be close to the on-shell condition, thus allowing us to use Eq. 23.

VII. CONCLUSIONS
Muons play a special role in the Standard Model of particles and elementary interactions. Some discrepancies between the observations and the SM predictions may be an indication of new physics. In particular, the coupling of muon to a new scalar particle (S) could reconcile some of the tensions between data and model. In this paper we elaborated on the idea that an O(10) GeV electron beam hitting the dump is a copious source of muons and electron beam-dump experiments have the ideal set up to detect the subsequent S → e + e − and S → γγ decays. A combination of simulation codes (FLUKA and Geant4) has been used to generate, propagate, and project muonproduced scalars on E137 and BDX detectors. The numerical simulation realistically takes into account the different processes involved: muons production from the primary electron, resulting energy and angular spread, S production, including production vertex distribution, S propagation and decay products distribution within the detector volumes. Results obtained for the E137 setup extended the exclusion limits, consistently with the reported null results, to cover a broader area in g µ vs. m S parameter space. An even larger exclusion zone was obtained considering the planned BDX experiment. In conclusion, we demonstrated that electron beam-dump experiments (past and future) provide an enhanced sensitivity to new physics that may be specifically coupled not only to electrons or photons but also to muons. and electron mass. For the specific case of hydrogen atom, we adopted the following form-factor parametrization [49,50] G el For each chemical element and m S value, we computed σ TOT as a function of E µ on a grid of 2000 points, from 10 MeV to 20 GeV. This interval corresponds to the energy range of interest in this work. The value of the coupling was fixed to g 0 µ = 3.87 · 10 −4 . The numerical integration was performed by first adapting the Vegas grid with it 1 = 10 runs, each with N 1 = 100k iterations, then performing the actual calculation through it 2 = 10 runs, each with N 2 = 100k iterations. The target accuracy was set to 1%. For each point, we also generated 100 · 10 3 Montecarlo events, that we used to sample the T S , θ S distribution. Figure 8 shows the total cross-section for the µ + N → µ+N +S process, as a function of the impinging µ kinetic energy, for different values of the scalar mass, and for the materials that are more relevant in the simulation.
In order to confirm the numerical results, we compared them with those obtained from MadGraph5 aMC@NLO 2.6.3.2 [51] (MG5 for short in below), finding an excellent agreement for all materials and scalar masses [52]. As an example, Fig. 9 shows the comparison between the total cross section obtained from CalcHEP and MG5 for Iron, as a function of the impinging muon kinetic energy, for the two values of the scalar mass m S = 10 MeV and m S = 100 MeV.
Appendix B: Geant implementation of scalar radiative emission process We implemented the scalar radiative emission process by muons in Geant4 as follows. A new class G4Scalar, inheriting from G4ParticleDefinition, was introduced to describe the scalar particle S. The class contains a singleton G4ParticleDefinition instance, that effectively introduces the new particle in the simulation. The singleton is instantiated through a static method with a single parameter, the mass of the scalar particle. This allows to set m S dynamically at the simulation start, in order to perform the multi-step numerical evaluation of the event yield, as described in Sect. IV.
The production process was implemented through the new class G4MuonScalarProduction, inheriting from G4VDiscreteProcess and implementing the following mother class virtual methods: • GetMeanFreePath, that returns the mean free path associated to the S production by a muon with given energy, in a given material. • PostStepDoIt, that specifies the actions to perform when a new scalar is produced, instructing the Geant4 application to create a new scalar particle at the interaction vertex with given momentum, sampled from the differential cross-section, and to alter the impinging muon momentum to ensure total four-momentum conservation.
By making use of a class inheriting from G4VDiscreteProcess, and adding it to the list of allowed physics processes for muons in the main physics list, the Geant4 kernel automatically handles the S production coherently with the other muons interaction and propagation mechanisms.
The mean free path λ for a muon with kinetic energy T µ propagating in a certain chemical element is computed from the relation λ −1 = nσ TOT , with n being the number of atomic nuclei per unit volume. The value of n is obtained from the internal Geant4 database, while σ TOT is loaded by the G4MuonScalarProduction class at creation, from pre-computed data -more details are given in App. A. In case of compounds, λ is obtained as an average: with the index "i" running over the chemical elements forming the compound, and σ i TOT , n i being the corresponding cross-section and atomic number density .
To describe the emission of the new scalar particle, first the particle kinetic energy T S and polar angle θ S with respect to the impinging muon direction are sampled from a pre-computed distribution, for the corresponding E µ value and chemical element. The azimuthal angle φ S is generated assuming a flat distribution. This allows to define the corresponding four-momentum P S . The impinging muon four-momentum is thus altered to enforce four-momentum conservation, P µ = P µ − P S (this assumes that the target nucleus stays at rest during the interaction). In case of compounds, the chemical element through which the impinging muon interacts is randomly selected from those forming the element, associating to each a statistical weight w i = n i σ i TOT . Then, the same procedure described before is adopted.