# Diphoton Signals of Muon-philic Scalars at DarkQuest

Nikita Blinov<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Victoria, Victoria, BC V8P 5C2, Canada  Department of Physics
and Astronomy, York University, Toronto, Ontario, M3J 1P3, Canada  Stefania
Gori<EMAIL_ADDRESS>Nick Hamer<EMAIL_ADDRESS>Department of Physics,
University of California Santa Cruz, 1156 High St., Santa Cruz, CA 95064, USA
and Santa Cruz Institute for Particle Physics, 1156 High St., Santa Cruz, CA
95064, USA

###### Abstract

We analyze the unique capability of the DarkQuest proton beam-dump experiment
at Fermilab to discover new light resonances decaying into photons. As an
example model, we focus on muon-philic scalar particles that decay to photons.
This is one of the few minimal models that can address the $(g-2)_{\mu}$
anomaly at low mass. These scalars can be copiously produced by meson decays
and muon bremsstrahlung. We point out that thanks to DarkQuest’s compact
geometry, muons can propagate through the dump and efficiently produce dark
scalars near the end of the dump. This mechanism enables DarkQuest to be
sensitive to both long-lived and prompt scalars. At the same time, di-photon
signatures are generically not background free, and we discuss in detail the
different sources of background and strategies to mitigate them. We find that
the backgrounds can be sufficiently reduced for DarkQuest to test currently-
viable $(g-2)_{\mu}$ parameter space.

##  I Introduction

The existence of a dark sector, neutral under the Standard Model (SM) gauge
symmetries, can offer an attractive explanation for Dark Matter (DM) and
presents a compelling avenue for new physics exploration. This could be
relevant for fundamental questions such as neutrino masses, the hierarchy
problem, and the Universe’s matter-antimatter asymmetry.

A vibrant set of efforts to look for sub-GeV dark sector particles - dark
photons, scalars, and a variety of DM models - is now underway Gori:2022vri ;
Cooley:2022ufh ; Coloma:2022dng . An especially interesting direction of
investigation has focused on the use of high-energy proton beams in fixed-
target setups to produce and detect long-lived particles Gardner:2015wea ;
Berlin:2018pwi ; Ahdida:2023okr . In particular, it has been shown that the
DarkQuest experiment, the upgrade of SpinQuest at Fermilab, will be able to
probe a plethora of dark sector models Apyan:2022tsd . So far, signatures with
electrons and muons have been mainly investigated, with studies focusing on DM
models, and on addressing the $(g-2)_{\mu}$ anomaly.

In fact, the measurement of the anomalous magnetic moment of the muon
Muong-2:2006rrc ; Muong-2:2021ojo ; Muong-2:2023cdq ,
$a_{\mu}\equiv(g-2)_{\mu}/2$, remains in significant tension with Standard
Model predictions. The difference between the Run I measurement and the
calculation from the $(g-2)_{\mu}$ theory initiative is Aoyama:2020ynm :

$a_{\mu}({\mathrm{Exp}})-a_{\mu}({\mathrm{Theory}})=(251\pm 59)\times
10^{-11}.$ (1)

While recent lattice QCD results may point toward a potential problem with
this SM prediction Borsanyi:2020mff ; Blum:2023qou ; Bazavov:2023has , new
particles beyond the SM (BSM) still provide a compelling solution to this
anomaly that is worth exploring further.

A large variety of new particle solutions have been proposed, including new
vector bosons, and CP-even and CP-odd scalar particles. Because their effects
on $(g-2)_{\mu}$ are through loops, this anomaly does not point to a very
specific new physics scale and the viable mass range is very broad, spanning
from $\sim$ MeV to several TeV Capdevilla:2021kcf . Therefore many experiments
running at different energies are required to comprehensively probe this
parameter space.

The minimal requirement for a new particle to alleviate the $(g-2)_{\mu}$
problem is a non-negligible interaction with muons. Realistic models, however
often contain couplings to other leptons, quarks or photons. These other
interactions can provide a powerful way of testing these models; for example
the minimal dark photon $(A^{\prime})$ explanation of $(g-2)_{\mu}$ has been
excluded in both visible NA482:2015wmo ; BaBar:2020jma  and invisible decay
channels BaBar:2017tiz , mainly by virtue of the $A^{\prime}$ coupling to
electrons and mesons. The searches utilized to probe these models, however, do
not directly probe the interaction with muons which would be desirable to make
robust claims about any BSM explanation of $(g-2)_{\mu}$ Chen:2017awl .

In this paper we consider searching directly for such a coupling using
secondary particle beams at the proposed proton beam dump DarkQuest experiment
Apyan:2022tsd . At this experiment, a high energy proton impinging on a target
produces a large flux of secondary mesons. If these mesons decay into muons
(like $\pi^{\pm}$ and $K^{\pm}$ do), they can produce muon-coupled force
carriers in two ways: the new particle can be produced in the meson decays
themselves, or via bremsstrahlung as the muons propagate through the beam
dump. Finally, some mesons can also decay into photons which can also produce
these particles in photonuclear reactions via a muon-loop-induced coupling.

We illustrate the sensitivity of DarkQuest in the context of the muon-philic
scalar model whose only tree-level interactions are with muons.111The
sensitivity of DarkQuest to a variety of dark sector scenarios has been
studied in Refs. Gardner:2015wea ; Berlin:2018pwi ; Berlin:2018tvf ;
Dobrich:2019dxc ; Tsai:2019buq ; Batell:2020vqn ; Blinov:2021say ;
Forbes:2022bvo ; Ariga:2023fjg .  Loop-suppressed couplings to photons are
naturally expected and give rise to scalar decays when its mass is below the
dimuon threshold. Thus the model’s main signature at a beam dump like
DarkQuest is a displaced decay into photons when the scalar is lighter than
$2m_{\mu}$. Ref. Forbes:2022bvo  studied the complimentary case where the
scalar is kinematically allowed to decay to muon pairs promptly. Our work
compliments the existing literature on proposed fixed-target probes of muon-
coupled particles at NA64, Fermilab, BDX, NA62, and SHiP Gninenko:2014pea ;
Chen:2017awl ; Berlin:2018pwi ; Chen:2018vkr ; Kahn:2018cqs ;
Marsicano:2018vin ; Krnjaic:2019rsv ; Rella:2022len ; Cesarotti:2023udo .
Compared to previous studies of similar experimental setups and models (e.g.,
Berlin:2018pwi  focusing on a leptophilic scalar model at DarkQuest and
Rella:2022len  focusing on NA64 and SHiP), we improve the modelling of known
production channels (muon bremsstrahlung) and identify additional relevant
channels (meson decays, $K,\pi\to\mu\nu S$222these decays have also been
studied in the context of kaon factories Krnjaic:2019rsv ), which turn out to
dominate the sensitivity projections. We also provide the first more in-depth
study of di-photon backgrounds at DarkQuest.

This work is organized as follows. In Section II we briefly describe the
DarkQuest experiment to establish key experimental parameters that determine
its sensitivity (for more details see Ref. Apyan:2022tsd ). We also estimate
the fluxes of secondary particles that are relevant for signal and background
production. We then specify the interactions of the muon-philic scalar model
in Section III and highlight the parameter space that can explain the
$(g-2)_{\mu}$ anomaly. Section IV details the modelling and simulation of
signal events in three channels: muon bremsstrahlung, meson decays, and
photon-induced processes. The compact nature of the DarkQuest experiment
enables unique sensitivity to short BSM particle lifetimes at the cost of
introducing several backgrounds in the case of photon signatures, which we
discuss in Section V. There we identify processes that can mimic our signal
and describe mitigation strategies. Finally in Section VI we present our
sensitivity projections, showing that DarkQuest can decisively probe parts of
the still-allowed $(g-2)_{\mu}$ parameter space.

We collect more details of our analysis in the appendices. Muon-philic scalar
decays and certain production channels rely on the loop-induced scalar-photon
coupling which is computed in Appendix A (while the corresponding on-shell
amplitude/form-factor is well known, we also provide a corrected off-shell
amplitude). Appendix B discusses a possible ultraviolet completion of the
muon-philic scalar model based on vector-like leptons, which can be
constrained with high energy colliders. Appendix C describes the propagation
of muons through the DarkQuest dump which is a key aspect of our
bremsstrahlung and muon-induced background calculations. In Appendix D we
present the full expression for the three-body kaon decay partial width which
dominates the signal reach.

##  II The DarkQuest Experiment

###  II.1 The Experimental Setup

Figure 1: Side view of the DarkQuest experiment including proposed the ECal,
target, and tracking detector upgrades, in red. A dark sector signature is
illustrated in blue. Adapted from Apyan:2022tsd .

DarkQuest is a dark sector-focused upgrade of SpinQuest, a proton fixed-target
spectrometer experiment on the neutrino-muon beamline of the Fermilab
Accelerator Complex. This beamline supplies a high-intensity beam of 120 GeV
protons from the Main Injector. Details of the set up are described in Refs.
SeaQuest:2017kjt  and Apyan:2022tsd . Here we summarize the key features
relevant for our analysis.

The setup of the DarkQuest experiment is shown in Fig. 1. A 5 m long
magnetized iron block (“FMAG”) is placed 1 m downstream from a thin nuclear
target, serving as a focusing magnet and as a beam dump. Its magnetic field of
$\sim 1.8$ T imparts a transverse kick of $\Delta p_{T}\sim 2.6$ GeV to
ultrarelativistic charge-one particles traversing the entire dump on average
(the actual $p_{T}$ for individual particles is sensitive to multiple
scattering and energy loss fluctuations in addition to the magnetic field, see
Appendix C). The FMAG stops most of the primary beam and secondary particles
produced by proton-iron interactions, other than neutrinos and muons; despite
this, the severely attenuated fluxes of protons and kaons can still be
relevant for background processes.

The detector extends for $\sim 20$ m after the dump and is composed of a
series of tracking stations and an open-aperture magnet (“KMAG”) that further
sweeps away charged SM radiation thanks to a transverse kick of $\Delta
p_{T}\sim 0.4$ GeV. The above components are already present in the SpinQuest
experiment. The DarkQuest upgrade will install: an electromagnetic calorimeter
(“ECal”) at $\sim 19$ m between station 3 and the muon station, additional
tracking layers to improve particle identification and hodoscopes for better
triggering. We will also consider modifications of this set up that include
more shielding before or after the FMAG. This is one strategy that can be used
to suppress some (but not all) backgrounds as we discuss in Section V.

The DarkQuest experiment will collect $\sim 10^{18}$ protons on target (POT)
in a two year long run. Depending on the status of the Fermilab accelerator
complex Gori:2024zbs , longer DarkQuest runs will be possible in the future.
In our analysis, we will discuss the realistic scenario of $10^{18}$ POT and
also a futuristic possibility in which DarkQuest will accumulate $10^{20}$
POT.

###  II.2 Production of SM Particles

The high intensity proton beam produces many SM particles at the beginning of
the dump, including mesons, baryons, muons, and photons. These secondary
particles open different production channels for the muon-philic scalar, and
also produce backgrounds. In Table 1, we summarize the average counts of the
relevant mesons and hadrons produced in the first proton interaction length
$\lambda_{\rm{int}}=16.77$ cm (approximately $1-e^{-1}\approx 0.63$ of the
primary proton beam interacts within this region). We estimated these yields
using Pythia 8.306 Bierlich:2022pfr .

Short-lived mesons like $\pi^{0}$ decay promptly after production. $\pi^{0}$
along with $\eta^{(\prime)}$ decays generate a flux of photons at the
beginning of the FMAG. As we will discuss, decays of the long-lived mesons
$\pi^{\pm}$ and $K^{\pm}$ contribute to signal production in two ways: first,
through rare three body decays, and second, by generating a flux of muons that
can undergo bremsstrahlung.333Muons are also produced in the Drell-Yan
process, but the corresponding cross-section is smaller than the total $pp$
cross-section by a factor of $\sim 10^{7}$, so this mechanism is subdominant
to production from meson decays. We focus on production and decay of these
mesons in the first meson interaction length for simplicity. This allows us to
neglect energy losses of hadrons that decay or interact deep in the dump. We
expect this to be a good $\mathcal{O}(1)$ approximation for the signal. This
is because while including a longer decay region for $\pi^{\pm}$, $K^{\pm}$
would increase the muon and scalar yields linearly with the decay region size,
the meson intensity falls exponentially with depth. Therefore decays in the
first interaction length account for the majority of the (scalar or muon) flux
that is relevant for the signal and background estimates. We estimate the muon
and photon yield, and the number of $\pi^{\pm}$ and $K^{\pm}$ decays in the
first interaction length in Table 2. The latter counts are obtained via

$N_{M,\;\mathrm{dec}}\simeq
N_{p}n_{M}\Gamma_{M}\langle\gamma_{M}^{-1}\rangle\lambda_{M}\,,$ (2)

where $N_{p}$ is the number of POT, $n_{M}$ is the meson yield per POT in the
first interaction length (from Table 1), $\Gamma_{M}$ and $\lambda_{M}$ are
the meson decay rate and interaction length, respectively. The average over
inverse boosts $\langle\gamma_{M}^{-1}\rangle$ is performed using Monte Carlo
samples from Pythia.

For certain backgrounds the attenuated proton beam or secondary $K_{L}$ and
$\Lambda$ fluxes _are_ an important source, and we will comment on them in
Section V.

$\pi^{0}$ | $\eta$ | $\eta^{\prime}$ | $\pi^{\pm}$ | $K^{\pm}$ | $K_{L}$ | $K_{S}$ | $\Lambda$
---|---|---|---|---|---|---|---
$1.5$ | $0.19$ | $0.021$ | $3.2$ | $0.25$ | $0.11$ | $0.21$ | $0.082$
Table 1: Average number of mesons or baryons produced in the first interaction length of the FMAG per 120 GeV proton on target. These counts are estimated using Pythia.  $\mu^{\pm}$ | $\gamma$ | $\text{prompt }\pi^{\pm}\text{ dec.}$ | $\text{prompt }K^{\pm}\text{ dec.}$
---|---|---|---
$1.0\times 10^{16}$ | $3.8\times 10^{18}$ | $8.6\times 10^{15}$ | $1.9\times 10^{15}$
Table 2: Average counts for particles that are produced or decay near the
front of FMAG for $10^{18}$ POT. The first two columns contain the number of
muons and photons produced in the first proton interaction length of the FMAG,
while the last two columns contain the number of prompt $\pi^{\pm}$ and
$K^{\pm}$ decays as estimated via Eq. 2. Here “prompt” refers to decays within
the first meson interaction length.

###  II.3 Di-photon mass resolution

As we will discuss in Section V, di-photon signatures are not background-free
at DarkQuest. The most important backgrounds is $\pi^{0}$ and $\eta$
production at the end of the dump produced either from the attenuated proton
beam or from the secondary muon beam undergoing deep inelastic scattering. One
experimental handle that can disentangle these events from the signal are
selections on diphoton invariant mass. Having a good mass resolution will
improve the discrimination between signal and background. The resolution can
be estimated using the energy and position resolution of the ECAL
Apyan:2022tsd . In particular, we take a position resolution of 3 cm and an
energy resolution of $1\%+7\%/\sqrt{E(\rm{GeV})}$.

By simulating $\pi^{0}$, $\eta^{(\prime)}$ production and their decays at the
back of the dump we find that the mass resolution varies between $10\%$ and
$30\%$ for typical meson energies, with a better resolution for the
$\eta^{(\prime)}$, if compared to the $\pi^{0}$.444We thank Yongbin Feng for
useful discussions regarding these calculations. We will find that even taking
the optimistic $10\%$ figure, selections of diphoton invariant mass can reduce
background rates by a factor of a few, but not completely eliminate them.
Other mitigation strategies will still be required (see Section VI).

##  III The Scalar Singlet Model

###  III.1 The Lagrangian of the Muon-philic Scalar

We focus on minimal scalar model that at the renormalizable level contains
only the interaction with muons

$\mathscr{L}\supset\frac{1}{2}(\partial
S)^{2}-\frac{1}{2}m_{S}^{2}S^{2}-g_{S}S\bar{\mu}\mu.$ (3)

We emphasize that this scenario does not have a coupling to electrons at tree
level; in the parlance of Ref. Chen:2017awl  this is “Model B”. The effective
theory also contains interactions with photons at dimension five

$\mathscr{L}\supset-\frac{1}{4}g_{S\gamma\gamma}SF_{\mu\nu}F^{\mu\nu}.$ (4)

The coupling $g_{S\gamma\gamma}$ receives contributions from ultraviolet (UV)
physics and from loops of muons:

$g_{S\gamma\gamma}=g_{S\gamma\gamma}^{(\mathrm{UV})}+\frac{\alpha g_{S}}{2\pi
m_{\mu}}f_{1/2}\left(\frac{4m_{\mu}^{2}}{m_{S}^{2}},0\right),$ (5)

where in writing the loop function $f_{1/2}$ we took $S$ and the photons to be
on-shell. We will neglect the UV contribution in what follows in order to be
consistent with previous studies. The usual explanation given for this is that
one naively expects $g_{S\gamma\gamma}^{(\mathrm{UV})}$ to be smaller by a
factor of $m_{\mu}/M$ where $M$ is the UV scale, see Appendix A. However, in
realistic UV completions, the $S$ coupling to muons (i.e., $g_{S}$) and the
heavy new physics scale, $M$, are independent, so the UV and infrared (IR)
contributions to $S\gamma\gamma$ can potentially compete for some choices of
parameters, as discussed in Appendix B. This photon interaction gives rise to
$S$ decays into photons with rate:

$\Gamma_{S}=\frac{\alpha^{2}g_{S}^{2}m_{S}^{3}}{256\pi^{3}m_{\mu}^{2}}\left|f_{1/2}\left(\frac{4m_{\mu}^{2}}{m_{S}^{2}},0\right)\right|^{2},$
(6)

where the loop function $f_{1/2}\to 4/3$ as $m_{S}/m_{\mu}\ll 1$ (the full
result is given in Appendix A). The muon loop also generates couplings of off-
shell photons to $S$ which are relevant for, e.g., the photon fusion and
Primakoff $S$ production. The $S\gamma\gamma^{*}$ amplitude needed for this
calculation is given in Appendix A.

Clearly the model in Eq. 3 is not electroweak gauge-invariant. From the point
of view of UV completions, the inclusion of couplings to electrons and taus is
natural in minimal models with Higgs portal interactions or extra Higgs
doublets above the electroweak scale Batell:2016ove . However, these models
have recently been excluded as solutions of $(g-2)_{\mu}$ BaBar:2020jma . The
_flavor-specific_ interaction in Eq. 3, while more minimal in the infrared,
requires more engineering in the UV. Ref. Batell:2017kty  delineated the
conditions under which couplings to specific SM fermion mass eigenstates (the
muon in our case) are obtained in a technically-natural and experimentally-
viable way. As a first step, the flavor-specific coupling to muons can be
obtained from the electroweak-invariant effective operator

$-\frac{S}{M}H^{\dagger}L\mathbf{c}_{S}\ell_{R}^{c}+\mathrm{h.c.},$ (7)

where $H$ ($L$) is the Higgs (lepton) $SU(2)_{L}$ doublet, $\ell_{R}^{c}$ is
the $SU(2)_{L}$ singlet lepton (in Weyl notation), $M$ is a UV scale, and
$\mathbf{c}_{S}$ is a matrix in flavor space. In order for this interaction to
reduce to Eq. 3, $\mathbf{c}_{S}$ has to obey two conditions: 1) it must be
diagonal in the lepton mass basis (i.e., $\mathbf{c}_{S}$ is simultaneously
diagonalizable with the lepton Yukawa matrix) and 2)
$\mathbf{c}_{S}=\operatorname{diag}(0,c_{S},0)$ in the charged lepton mass
basis. This structure generates the desired coupling to muons with
$g_{S}=c_{S}v/(\sqrt{2}M)$, where $v=246$ GeV, and it minimizes flavor-
changing neutral currents. In general, $c_{S}$ is complex, but its imaginary
part shifts $(g-2)_{\mu}$ in the wrong direction Batell:2017kty . We therefore
focus on purely real $c_{S}$. The spurion analysis of Ref. Batell:2017kty
reveals that the intricate structure of $\mathbf{c}_{S}$ is radiatively
stable.555The parity $S\to-S$ ensures that corrections to $c_{S}$ are always
proportional to $c_{S}$, while any off-diagonal elements can only arise from
the non-diagonal neutrino mass matrix.

The electroweak-invariant interaction in Eq. 7 is still non-renormalizable and
therefore requires a UV completion. This is easily achieved by adding a pair
of vector-like leptons with mass $M$ Batell:2017kty , in which case
$\mathbf{c}_{S}$ is proportional to the Yukawa couplings of these new
particles with the SM Higgs and with $S$. We explore such a UV completion in
Appendix B, which allows us to connect the low energy muon and photon
couplings to UV parameters, and to consider naturalness and direct
experimental bounds that must be satisfied. This scenario still does not
explain the flavor alignment of $\mathbf{c}_{S}$ and the SM lepton Yukawa
matrices.

###  III.2 The New Physics Contribution to $(g-2)_{\mu}$

The model in Eq. 3 is one of the few available minimal explanations of the
$(g-2)_{\mu}$ anomaly with new physics at or below the several GeV scale
Capdevilla:2021kcf . The contribution of $S$ to $(g-2)_{\mu}$ is
Batell:2016ove

$\Delta
a_{\mu}^{S}=\frac{g_{S}^{2}}{8\pi^{2}}\int_{0}^{1}dx\frac{(1+x)(1-x)^{2}}{(1-x)^{2}+x(m_{S}/m_{\mu})^{2}}.$
(8)

The anomaly in Eq. 1 is explained if $g_{S}\sim 4\times 10^{-4}$ for
$m_{S}\lesssim 2m_{\mu}$.

##  IV Signal Production and Simulations

The dark scalar can be produced in several processes enabled by secondary
muons, mesons and photons present in the beam dump. The energy and transverse
momentum spectra of these particles are shown in Fig. 2. We find that muon
bremsstrahlung and meson decays dominate the dark scalar yield in interesting
regions of parameter space, while the rate of photon-initiated reactions is
loop-suppressed. We discuss these processes in more detail below.

Figure 2: Counts of secondary particles per 120 GeV proton on target as a
function of energy (left panel) and transverse momentum (right panel). The
$\pi^{\pm}$, $K^{\pm}$ ($\gamma$, $\mu^{\pm}$) spectra sum to the average
yields in the first interaction length of the dump given in Table 1 (Table 2,
divided by $10^{18}$). Note that the muon spectrum includes the probability
for the parent particles (mostly $\pi^{\pm}$ and $K^{\pm}$) to decay in the
first interaction length (this spectrum is also multiplied by 100 for
convenient visualization).

###  IV.1 Muon Bremsstrahlung

The primary proton collisions with the dump lead to a large flux of secondary
muons. These muons are generated from meson decays ($\pi^{\pm}$ and $K^{\pm}$
in particular), and from Drell-Yan production. These muons can then interact
within the dump, radiating a scalar. The scalar bremsstrahlung cross-section
has been computed exactly in Tsai:1986tx ; Tsai:1989vw ; Liu:2016mqv , but
most recent works use the Weiszacker-Williams (WW) approximation, see, e.g.,
Ref. Chen:2018vkr ; Rella:2022len . This cross-section enters the signal
prediction calculation in two important ways: first, it normalizes the overall
rate and, second, its differential is used to generate $S$ kinematics. We have
compared three ways of evaluating the overall cross-section: (1) using the WW
approximation; (2) numerically computing it in MadGraph Alwall:2014hca ; (3)
integrating the exact results of Ref. Liu:2016mqv  using VEGAS Lepage:1977sw ;
Lepage:2020tgj . In principle the latter two methods should give very similar
results, but in practice we find that the custom VEGAS implementation yields
results that are more numerically stable.666The difficulty of simulating
light-particle bremsstrahlung using off-the-shelf Monte Carlo tools has been
noted in, e.g., Celentano:2020vtu ; Forbes:2022bvo . The result of this
calculation is illustrated in Fig. 3 for several values of $m_{S}$ \- we find
excellent agreement with results of Ref. Marsicano:2018vin  which used a
different numerical method, which is also based on the integration of the full
matrix element squared. For the main results in this work we use the VEGAS
calculation of the total cross-section, noting that the WW approximation
predicts scalar yields that are a factor of 2-5 larger for $m_{S}\lesssim
2m_{\mu}$ (such discrepancies were also observed in Refs. Liu:2016mqv ;
Rella:2022len ). We will use MadGraph to generate final state scalar and muon
kinematics.

Figure 3: Scalar bremsstrahlung ($\mu N\to\mu NS$) cross-section, in iron as a
function of muon energy for several choices of $m_{S}$. Note that to obtain a
realistic scalar yield at DarkQuest for this process one must convolve these
cross-sections against the muon spectrum and account for muon propagation
effects.

Most of the muons are produced near the front of the dump, so they can
propagate a significant distance before undergoing bremsstrahlung. During this
propagation they experience energy losses, multiple Coulomb scattering, and
their trajectory is bent in the magnetic field. There have been two different
approaches to modelling these effects. First, input muon spectrum can be used
to predict the signal yield by convolving it with the differential cross-
section Berlin:2018pwi ; Rella:2022len

$\displaystyle N_{\text{signal}}\simeq n_{\text{atom}}\int
dE_{\mu}^{0}\frac{dN_{\mu}}{dE_{\mu}^{0}}\int_{E_{\mu}^{\text{min}}}^{E_{\mu}^{0}}dE_{\mu}\frac{\sigma_{\text{Brem.}}(E_{\mu})}{|dE_{\mu}/dz_{\mu}|}\int_{z_{\text{min}}-z_{\mu}}^{z_{\text{max}}-z_{\mu}}dz\frac{1}{N_{\text{MC}}}\sum\limits_{\text{events
$\in$ geom.}}\frac{e^{-z/(\gamma_{S}c\tau_{S})}}{\gamma_{S}c\tau_{S}}.$ (9)

where $z_{\rm{min}}\simeq 5$ m, $z_{\rm{max}}\simeq 19$ m and $N_{\rm{MC}}$ is
the total number of simulated events. The sum is performed over only those
events that remain within the geometry of the DarkQuest spectrometer. Here
$E_{\mu}^{0}$ is the initial energy of the muon at production, and $E_{\mu}$
is the muon energy after traversing some finite distance, $z_{\mu}$, in the
FMAG before radiating the scalar; $\gamma_{S}$ and $\tau_{S}$ are the boost
and the lifetime of the scalar, respectively; $dN_{\mu}/dE_{\mu}^{0}$ is the
total number of muons produced per initial energy bin within the first meson
interaction length (output of Pythia); $n_{\rm{atom}}=8.5\times 10^{22}$/cm3
is the number density of target iron nuclei and
$\sigma_{\text{Brem.}}(E_{\mu})$ is the muon bremsstrahlung cross-section, as
a function of the muon energy (see Fig. 3). Note that in this approach muon
propagation is simply encoded in the $dE_{\mu}$ integration, and the effect of
a magnetic field is applied “by hand” as a shift in $p_{T}$ of each muon.
Multiple scattering is not included.

Alternatively, each muon can be tracked through the dump using its equations
of motion, until it undergoes a bremsstrahlung event. At this point $S$
kinematics is sampled from the differential cross-section. This approach makes
it easier to include the effects of multiple Coulomb scattering, energy loss
fluctuations, and magnetic fields. In Ref. Marsicano:2018vin  this method was
implemented within GEANT4 GEANT4:2002zbu ; Allison:2006ve ; Allison:2016lfl .
Here we follow a similar approach, but instead model the propagation in a
custom Python program described in Appendix C.

Finally, we note that muon bremsstrahlung in a thick target has the peculiar
feature that scalar particles can be produced anywhere in the dump with
similar probability, including very close to the end ($z_{\mu}\sim
z_{\rm{min}}$). This can be easily seen from Eq. 9 in the limit of large muon
energy (such that one can neglect energy losses) and short scalar lifetime
($\gamma_{S}c\tau_{S}\ll z_{\mathrm{min}}$). In this case the probability of
the scalar to be produced in the dump and decay beyond it is approximately

$\begin{split}n_{\mathrm{atom}}\sigma_{\text{Brem.}}(E_{\mu})\int_{0}^{z_{\mathrm{min}}}dz_{\mu}e^{-(z_{\mathrm{min}}-z_{\mu})/(\gamma_{S}c\tau_{S})}\\\
\approx
n_{\mathrm{atom}}\sigma_{\text{Brem.}}(E_{\mu})\gamma_{S}c\tau_{S},\end{split}$
(10)

where we assumed that the cross-section and energy loss rates are constant.
This shows that the signal yield will be dominated by the scalar production in
the last $S$ decay length of the shield. In contrast, if the scalar was
produced near the front of the dump, the rate would be proportional to
$e^{-z_{\mathrm{min}}/(\gamma_{S}c\tau_{S})}$. Therefore the projected
sensitivity of this channel extends to large couplings, which is qualitatively
different to many other long-lived-particle (LLP) searches which feature a
very thick shield. For example, the proton beam-dump experiment CHARM
CHARM:1983ayi ; CHARM:1985nku ; Gninenko:2012eq  had a shield thickness of
$\sim 480$ m so that the LLP signal in the large coupling/short lifetime
regime is exponentially suppressed despite the higher beam energy of $400$
GeV. Note that potential signals from secondary muons (which can propagate
over such distances) have not been studied, but the muons would likely
experience significant energy loss, making it difficult for signals to pass
detector selections.

###  IV.2 Meson Decays

Dark scalars are also produced from the charged current meson decays
$\pi^{\pm}\to\mu^{\pm}\nu S$ and $K^{\pm}\to\mu^{\pm}\nu S$. Using the
leading-order chiral perturbation theory Lagrangian

$\mathscr{L}_{\chi\mathrm{PT}}\supset\sqrt{2}G_{F}f_{P}V_{q_{1}q_{2}}(\partial_{\alpha}P^{-})(\bar{\mu}\gamma^{\alpha}P_{L}\nu_{\mu})+\mathrm{h.c.},$
(11)

we find the partial width for meson $P=\pi,K$

$\displaystyle\Gamma(P^{\pm}\to\mu^{\pm}\nu S)$
$\displaystyle\approx\frac{g_{S}^{2}f_{P}^{2}G_{F}^{2}|V_{q_{1}q_{2}}|^{2}m_{P}^{3}}{768\pi^{3}}$
(12a)
$\displaystyle\approx\frac{g_{S}^{2}m_{P}^{2}}{96\pi^{2}m_{\mu}^{2}}\Gamma(P^{\pm}\to\mu^{\pm}\nu),$
(12b)

in the $m_{P}\gg m_{S},m_{\mu}$ limit (the full expression is given in
Appendix D). Here $f_{P}$ is the meson decay constant
($f_{\pi}=130\;\mathrm{MeV}$ and $f_{K}=156\;\mathrm{MeV}$ in our
conventions); $V_{q_{1}q_{2}}$ is the relevant CKM matrix element:
$|V_{ud}|\approx 0.9737$ ($|V_{us}|\approx 0.2246$) for $\pi$ ($K$) decay.
These are the most important meson production modes for the scalar thanks to
their sizable branching ratio into a muon and a neutrino.

An additional production mode is $K_{L}\to\pi^{\pm}\mu^{\mp}\nu_{\mu}S$. The
contribution from this process, however, is generically smaller than from the
meson three body decays discussed above.

###  IV.3 Photon-induced Processes

The scalar-photon interaction in Eq. 4 enables several production mechanisms
that are analogous to axion-like particle production in beam dumps. These
include the Primakoff process $\gamma+N\to S+N$ and photon fusion
$\gamma^{*}\gamma^{*}\to S$, where $\gamma^{*}$ are virtual photons
corresponding to the electromagnetic fields of a beam proton and target
nucleus. Primakoff production was shown to dominate over photon fusion in
thick targets due to the large secondary photon flux Blinov:2021say . The
scalar-photon interaction needed to evaluate the Primakoff cross-section is
described in Appendix A. A potential difference of the present model with
respect to axion-like particles is that the photon coupling is generated by
muons, leading to a non-trivial energy dependence of that coupling (because
the muon mass is not much larger than other experimental energy scales). We
have estimated the yield of dark scalars from this process, finding that it is
1–2 orders of magnitude smaller than for bremsstrahlung in the mass range of
interest.777Moreover, the corresponding signal would not have a on-shell muon
in the final state which may be an important event selection designed to
reduce diphoton backgrounds as we discuss in the next two sections.

##  V Backgrounds

The compact geometry of DarkQuest enables searches for particles with short
decay lengths (especially if they are radiated by muons propagating though the
dump). This unique feature leads to a complementary reach compared to other
experiments with a larger dump and smaller angular acceptance. The
comparatively short shield also means that potential signals must contend with
backgrounds. There is a variety of SM processes that can mimic the production
and decay of a BSM long-lived particle into photons. The three types of
potentially problematic SM events are: 1) production of genuine SM long-lived
particles, such as $K_{L}$; 2) neutral mesons ($\pi^{0}$, $\eta^{(\prime)}$)
produced at the very back of the dump from the attenuated proton beam; and 3)
production of neutral mesons ($\pi^{0}$, $\eta^{(\prime)}$) via Deep-Inelastic
Scattering (DIS) of muons towards the end of the dump. We will show that the
events of the third kind are the most dangerous, but that they can be
mitigated with several event selections, while preserving experimental
coverage of the most interesting $(g-2)_{\mu}$ band at $m_{S}\lesssim
2m_{\mu}$. These backgrounds and potential mitigation strategies are
summarized in Table 3 – they will be discussed in more detail below and in
Section VI.

| Raw counts | ECAL hit |  $\gamma$ sep. |  $\mu$ hit | $m_{\gamma\gamma}$ | Extra Shielding
---|---|---|---|---|---|---
$p$-induced $K_{L}$ | $6\times 10^{4}$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$
$p$-induced $\pi^{0},\;\eta^{(\prime)}$ | $2\times 10^{4}$ | $\checkmark$ |  | $\checkmark$ | $\checkmark$ | $\checkmark$
$\mu$-DIS-induced $\pi^{0},\;\eta^{(\prime)}$ | $10^{6}$ | $\checkmark$ |  |  | $\checkmark$ |
Table 3: Summary of SM background events that can mimic the signal
$S\to\gamma\gamma$ and experimental techniques that can reduce or eliminate
them. The first two rows correspond to mesons produced by the primary proton
beam. The last row corresponds to mesons produced deep in the dump by a
secondary muon beam via deep-inelastic-scattering. The columns to the right of
the double-line divider indicate whether the listed background events can be
significantly reduced using the corresponding experimental handle. “ECAL hit”
requires at least one photon with a sufficient energy in the ECAL; “$\gamma$
sep.” corresponds to exactly two well-separated photons in the ECAL; “$\mu$
hit” is the requirement of detecting a muon; “$m_{\gamma\gamma}$” is a
selection based on the diphoton invariant mass potentially measurable with the
ECAL; “Extra shielding” corresponds to adding additional material behind the
FMAG. The impact of these selections on the background yield is quantified in
Section V and Section VI.

###  V.1 SM Long-Lived Particles

Standard Model LLPs, such as $K_{L}$ and $\Lambda$, which decay into neutral
pions can produce signal-like photons. We will focus on $K_{L}$ decays in the
fiducial region $z>5$ m and $z<19$ m (before the ECal). $\Lambda$’s have a
shorter attenuation length compared to mesons, so their contribution to
background ends up being two orders of magnitude smaller. The number of LLP
background events can be estimated as

$\begin{split}&n_{K_{L}}N_{p}\mathrm{BR}(K_{L}\to
3\pi^{0})e^{-z_{\mathrm{shield}}/\lambda_{K\;\mathrm{int}}}\\\
&\times\left[e^{-z_{\mathrm{shield}}/(\gamma
c\tau_{K_{L}})}-e^{-z_{\mathrm{max}}/(\gamma c\tau_{K_{L}})}\right]\\\
&\approx 6\times 10^{4},\end{split}$ (13)

where $n_{K_{L}}$ is the number of $K_{L}$ produced per POT from Table 1,
$\mathrm{BR}(K_{L}\to 3\pi^{0})\approx 0.2$ is the branching fraction of the
most important $K_{L}$ decay channel that gives photons in the final state,
$z_{\mathrm{shield}}=5$ m is the length of shielding,
$\lambda_{K\;\mathrm{int}}\approx 20\;\mathrm{cm}$ is the kaon interaction
length in iron. The factor on the second line is the probability of the kaons
to decay after the shielding and before the ECAL. The number of events in the
last line is computed for $N_{p}=10^{18}$ POT. Each of these events can lead
to a $\leq 6$ photons signal in our detector. Note that most $K_{L}$’s are
produced near the front of the dump as the incoming proton beam is
exponentially attenuated.888$K_{L}$’s produced deeper in the dump experience
less attenuation before they exit which somewhat off-sets the attenuation of
the beam that produced them, compared to $K_{L}$’s near the front. However,
because the kaon interaction length is longer than the proton interaction
length, the latter $K_{L}$ flux still dominates.

There are two approaches to deal with this background. First, one can use
additional shielding in front or behind the FMAG (“Extra shielding” in Table
3). Using Pythia to estimate the median meson production rates and boosts, we
find that an additional 2 m of (iron) shielding would bring the $K_{L}$
background rate below 10 events for $10^{18}$ POT. Similar shielding
requirements were discussed in Refs. Berlin:2018pwi ; Blinov:2021say  in the
context of axion-like particle searches. It is important to note that if the
only goal is to close the $(g-2)_{\mu}$ window, we can significantly relax the
shielding requirement. This is because we expect $\mathcal{O}(10^{3}-10^{6})$
signal events in the currently-allowed $(g-2)_{\mu}$ band at $m_{S}\lesssim
2m_{\mu}$ (see Fig. 8). The amount of extra shielding needed to bring the
background below $10^{3}$ events is only $90$ cm. This can be further reduced
by using a different material, such as tungsten.

The second approach to reducing SM LLP backgrounds is to select events with
exactly two photons in the ECAL and a single muon hit (“$\gamma$ sep.” and
“$\mu$ hit” in Table 3). Both the kaon decay and muon bremsstrahlung $S$
production mechanisms can in principle pass these cuts, while the backgrounds
either have too many photons, do not have a muon in the final state or both.
The leading $K_{L}$ decay that leads to background events is $K_{L}\to
3\pi^{0}\to 6\gamma$, which would be completely eliminated by such cuts,
without requiring any additional shielding. In Section VI, we will show that
enough signal events pass these selections to retain sensitivity to the open
$(g-2)_{\mu}$ region at $m_{S}\lesssim 2m_{\mu}$.

###  V.2 Secondary Mesons from the Back of the Dump

Another type of potentially problematic SM events is the production of
$\pi^{0}$ and $\eta^{(\prime)}$ from the attenuated proton beam in the last
radiation length of the dump. The prompt photons from their decays would have
a large probability of escaping the dump and potentially ending up in the
ECAL. The number of such events is approximately

$\begin{split}&(n_{\pi^{0}}+n_{\eta}+n_{\eta^{\prime}})N_{p}e^{-z_{\rm{shield}}/\lambda_{p\;\rm{int}}}\frac{\lambda_{\rm{rl}}}{\lambda_{p\;\rm{int}}}\\\
&~{}~{}~{}~{}\approx 2\times 10^{4},\end{split}$ (14)

where $\lambda_{p\;\rm{int}}$ = 16.77 cm is the nuclear interaction length for
protons, $n_{M}$ are given in Table 1 and $\lambda_{\rm{rl}}$= 1.757 cm is the
radiation length in iron.999The attenuated proton beam looses $\sim 10$ GeV of
energy as it propagates through the FMAG, but the meson yields at this lower
energy are approximately the same as for in original 120 GeV beam. The last
factor, $\lambda_{\rm{rl}}/\lambda_{p\;\rm{int}}$, is the probability for
protons to interact in the final radiation length of the dump. The background
events in the last line of LABEL:eq:LLPBack can be eliminated in a few
different ways. Thanks to the exponential suppression in LABEL:eq:LLPBack, an
additional 1.8 m (70 cm) of iron shielding would suffice to reduce this
background to below 1 ($10^{3}$) event for $N_{p}=10^{18}$. Alternatively, as
for SM LLPs, event selections requiring two photon hits in the ECAL and a muon
hit would also help thanks to a requirement on the di-photon invariant mass
and distributions, see Section VI. We will investigate how such selections
affect the signal yield in the next section.

###  V.3 Mesons from $\mu$ DIS

Finally, we consider backgrounds generated by muons undergoing DIS at the end
of the dump and producing $\pi^{0}$ and other mesons that can decay to
photons. The reason why this reaction chain is potentially important is
because muons are not exponentially attenuated by the dump, unlike the
original proton beam (c.f., the secondary meson background in Section V.2). In
order to estimate the rate, we first produce a distribution of muons at the
end of the dump as follows. We simulate the primary $pp$ collisions in Pythia
that produce $\pi^{\pm}$ which are then decayed to muons in the first
interaction length. The muons are then propagated to the back of the dump as
described in Appendix C. We find that the number of muons at the back of the 5
m dump with sufficient energy to initiate DIS is
$n_{\mu}^{(\mathrm{back})}\sim 5\times 10^{-6}$ per $pp$ interaction. Note
that this number is not sensitive to the dump thickness because the muons are
not strongly attenuated. This is because the muons lose only about $\sim
10\;\mathrm{MeV}/\mathrm{cm}$, so even an additional 2 m of iron does not
drastically affect their spectrum at the back of the dump.

The propagated sample of muons is then fed into Pythia to simulate DIS, $\mu
p\to\mu+X$, which is used to estimate the meson-count weighted cross-sections,
$\langle\sigma(\mu p\to\mu+X)n_{M}\rangle$ where $M=K_{L},\;\pi^{0}$. We find
that

$\begin{split}\langle\sigma(\mu p\to\mu+X)n_{\pi^{0}}\rangle&=2\times
10^{-5}\;\mathrm{mb},\\\ \langle\sigma(\mu p\to\mu+X)n_{K_{L}}\rangle&=5\times
10^{-7}\;\mathrm{mb}.\end{split}$ (15)

We can now put all of these ingredients together. First, the expected number
of $\mu$-DIS-induced $\pi^{0}$ produced in the final radiation length of the
dump is

$N_{p}n_{\mu}^{(\mathrm{back})}An_{A}\lambda_{\rm rl}\langle\sigma
n_{\pi^{0}}\rangle\sim 10^{6}$ (16)

where $\langle\sigma n_{\pi^{0}}\rangle$ is given in Eq. 15. The target
(nucleon) density in iron is $An_{A}\approx 4.7\times
10^{24}\;\mathrm{cm}^{-3}$.

The number of kaons that are produced in a similar manner in the last kaon
interaction length, and that decay beyond the dump is

$\begin{split}&N_{p}n_{\mu}^{(\mathrm{back})}An_{A}\lambda_{K\;\mathrm{int}}\langle\sigma
n_{K_{L}}\rangle\mathrm{Br}(K_{L}\to 3\pi^{0})\\\
&\times\left[1-e^{-(z_{\mathrm{max}}-z_{\mathrm{dump}})/(\gamma
c\tau_{K_{L}})}\right]\sim 10^{4}\end{split}$ (17)

where the quantity the brackets is the probability of $K_{L}$ produced at the
back of the dump to decay between $z=z_{\mathrm{dump}}$ and
$z_{\mathrm{max}}=19\;\mathrm{m}$ (about 0.14 for the mean kaon energy
produced in the DIS reaction).

While the raw counts for the backgrounds in Eqs. 16 and 17 appear problematic,
it turns out that simple selections on photon energy and angle distributions
can make these backgrounds significantly smaller. Nevertheless, as we will
discuss in the next section, it will be challenging to completely suppress
this source of background.

Figure 4: Scalar yield as a function of the mass $m_{S}$ for $10^{18}$ POT
from different production mechanisms without imposing any event selections.

Figure 5: Histograms of signal photon energy and angle from $S\to\gamma\gamma$
with $S$ produced in rare kaon decays (left panel) and muon bremsstrahlung
(right panel). $\theta_{\gamma}$ is the angle between the photon and the beam
axis. In both panels $m_{S}=100$ MeV and the bin normalization is such that
they sum to 1. The vertical dashed white line corresponds to an energy cut of
1 GeV and the horizontal line corresponds to photons being roughly in
geometric acceptance of the detector for scalars decaying just after the dump.
Figure 6: Acceptance probability of signal photons as a function of scalar
mass for different $S$ production channels, decay positions and event
selections. Dashed curves indicate a further $5.5$ cm separation cut was
applied to the photons hitting the detector.

##  VI Results and Discussion

Figure 7: Projected sensitivity of the DarkQuest experiment to the muon-philic
scalar model in the $m_{S}-g_{S}$ plane for $m_{S}\leq 2m_{\mu}$. In the left
panel we show contributions of different production channels to the
sensitivity for $10^{18}$ POT (the dashed line is for $10^{20}$ POT). The
signal photons are required to have $E_{\gamma}>1\;\mathrm{GeV}$ and to be
separated by more than 5.5 cm in the transverse direction when they enter the
ECAL. In the right panel we show how the sensitivity changes with a lower cut
on $E_{\gamma}$, and how the reach is improved with $10^{20}$ POT. In both
plots the green band is the parameter space that can address the $(g-2)_{\mu}$
anomaly; the gray regions are excluded either by $(g-2)_{\mu}$ or by the
recast of E137 Marsicano:2018vin . We also show projections for the proposed
resonant search at SpinQuest Forbes:2022bvo  and at NA62 Krnjaic:2019rsv .
Figure 8: Number of signal events for $10^{18}$ POT, after requiring two
photons with $E_{\gamma}>1\;\mathrm{GeV}$ each, separated by $\geq 5.5$ cm.

(a)  (b)

(c)  (d)

Figure 9: Normalized histograms of kinematic variables for a parameter point
at high scalar mass, inside the $(g-2)_{\mu}$-favored region
($m_{S}=0.14\;\mathrm{GeV},g_{S}\sim 10^{-3}$). The bremsstrahlung and the
three body kaon decay $K\to\mu\nu S$ events are selected by requiring each
photon in the pair have $E_{\gamma}>1\;\mathrm{GeV}$ and $\geq 5.5$ cm
separation in the ECal. The $K_{L}$ background events are selected by grouping
the 6 photons into 5.5 cm square calorimeter towers and requiring that two
distinct towers have a total $E_{\gamma}>1\;\mathrm{GeV}$. The variables which
are being plotted are then totaled for that tower.  Figure 10: Comparison of
the (background-free) experimental reach for different signal event selections
for the bremsstrahlung production mechanism. All projections include a minimum
photon energy cut of $1$ GeV; the green dashed (orange dotted) line also
requires transverse photon separation of at least $5.5$ cm (a muon hit with
$E_{\mu}>1\;\mathrm{GeV}$ in the muon station). The pink dot-dashed line
combines all of these selections. Even for this more restrictive set of cuts
designed to eliminate the SM LLP background, DarkQuest maintains sensitivity
to the $(g-2)_{\mu}$ band at $m_{S}\lesssim 2m_{\mu}$.

In this section we put together our signal and background simulations to
project the sensitivity of future DarkQuest searches for muon-philic scalars
in the $m_{S}\lesssim 2m_{\mu}$ regime. First, we compare the raw $S$ yields
from different production channels without making any selections in Fig. 4. We
observe that production from meson decays dominates the $S$ yield throughout
our mass range. While rare $\pi^{0}$ decays produce the largest $S$ flux for
$m_{S}\lesssim m_{\pi}-m_{\mu}$, we will see that even minimal event
selections make this channel less important than the other two. Kaon decays
produce more $S$ that $\mu$-bremsstrahlung despite the initial flux of $\mu$
being larger – see Table 2. This occurs for two related reasons: first, the
spectrum-averaged muon interaction probability in the dump is smaller than the
$K^{\pm}$ branching fraction for the $m_{S}$ of interest; second, while the
bremsstrahlung cross-section grows with initial muon energy, the probability
of those muons to be produced in pion decays near the front of the dump is
smaller, because of their higher boost. Thus the interaction probabilities are
anti-correlated with the energy-differential muon flux, which is not captured
in Table 2.

Examples of distributions of the energy and the angle between the photon and
the beam axis are shown in Fig. 5 for the two dominant production mechanisms.
The probability of these photons ending up in the detector acceptance depends
on the precise location of the $S$ decay and any selections imposed on the
events. In Fig. 6 we show the effects of several simple cuts on energy and
photon separation for a few fixed decay positions. We consider minimum photon
energy cuts of 0.5 and 1 GeV as representative values, but looser selections
are also possible Apyan:2022tsd . Requiring photon separation will enable the
selection of events with exactly two photons which may be a useful handle for
rejecting backgrounds such as $K_{L}\to 3\pi^{0}$. The choice of 5.5 cm
minimum transverse separation corresponds to the transverse granularity of the
ECAL towers used by DarkQuest PHENIX:2003fvo ; Apyan:2022tsd . In Fig. 6 we
see that photon separation strongly reduces sensitivity to light scalars. As
we discuss next this is not an issue since the novel parameter space that
DarkQuest can test lies at larger masses. In Fig. 7 we show background-free
sensitivities for easy comparison with previous literature. The reach curves
correspond to 10 signal event contours in the $m_{S}-g_{S}$ plane for phase 1
($N_{p}\approx 10^{18}$) and phase 2 ($N_{p}\sim 10^{20}$). One remarkable
feature of these projections is the sensitivity to relatively large couplings
where the lab-frame lifetime of $S$ is expected to be significantly shorter
than the FMAG depth of $5$ m. It arises due to the secondary muons’ ability to
interact deep in the dump, leading to $S$ production much closer to $z=5$ m as
discussed in Section IV. From the figure, we can see that already phase 1
could probe entirely the $(g-2)_{\mu}$ region at higher masses. Phase 2 will
also have access to lower masses.

The projections in Fig. 7 are valid only if the backgrounds discussed in
Section V are reduced or eliminated. A simple, brute-force way to achieve this
for the SM LLP backgrounds and diphotons from meson produced by the proton
attenuated beam in the back of the dump is to include substantial extra
shielding. Such a modification of the DarkQuest setup may not be practically
feasible so it is worth considering other mitigation strategies.

We note that while DarkQuest is sensitive to very small couplings, it does not
surpass the (recasted) limits from E137 Marsicano:2018vin  in the long-
lifetime regime for $N_{p}=10^{18}$ POT. In fact, the biggest impact of
DarkQuest will be to probe the parameter space with $m_{S}\lesssim 2m_{\mu}$
that can explain $(g-2)_{\mu}$; in this regime the couplings are substantial
$g_{S}\lesssim 10^{-3}$, leading to signal rates of
$\mathcal{O}(10^{3}-10^{6})$ as we show in Fig. 8. This suggests that the
analysis can tolerate additional cuts that lower signal acceptance while still
testing $(g-2)_{\mu}$. We delineate a few background reduction strategies
below.

$N_{\textrm{towers}}$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$
---|---|---|---|---|---|---|---
$E^{\textrm{tower}}_{\textrm{cut}}=1.0\mathrm{GeV}$ | $55\%$ | $8\%$ | $8\%$ | $9\%$ | $10\%$ | $7\%$ | $3\%$
Table 4: Distribution of the number of distinct calorimeter towers whose total
energy deposition is above $E^{\textrm{towers}}_{\textrm{cut}}$ for photons
produced from the $K_{L}\to 3\pi^{0}$ background. Similar numbers are obtained
using the 0.5 GeV cut on the photons.

Requiring exactly two calorimeter towers with an energy deposition of at least
1 GeV eliminates a large fraction of the $K_{L}\to 3\pi^{0}$ background (see
Table 4 with a $92\%$ reduction). The remaining events have photon energy and
$p_{T}$ distributions that are drastically different from those of the signal
as we show in Fig. 9. This suggests additional selections that can further
reduce this background. For example by increasing the minimum photon energy
cut to 2.8 GeV one can reduce the $K_{L}$ backgrounds by an order of
magnitude. This selection can be further optimized by using different cuts for
the lower and higher energy photons in the ECAL. For example, a lower energy
photon cut of $E>4$ GeV together with a higher energy photon cut of $E>10$ GeV
would reduce the background by a factor of 50 while reducing signal by a
factor of 3.

Diphoton production in the last radiation length of the dump from the
attenuated proton beam (Section V.2) and displaced decays of $K_{L}$ (Section
V.1) do not produce associated muons, while the dominant $S$ production
channels do. Thus, selecting events with a muon can be an effective way of
eliminating these backgrounds. In Fig. 10 we show how the background-free
sensitivity changes with the addition of various cuts, including the muon
selection. While requiring a muon in the acceptance of the detector incurs a
significant penalty on the signal rate, DarkQuest remains sensitive to the
parameter space that can explain $(g-2)_{\mu}$. The impact of the muon
selection requirement was estimated by propagating each muon after it
underwent bremsstrahlung in the FMAG through the remainder of the dump,
through KMAG and to the muon station at $\sim 21$ m. The muon hit was counted
if it had an energy over 1 GeV and it intersected the station 4 detectors
SeaQuest:2017kjt  (labelled “muon station” in Fig. 1).

Requiring a muon hit does not eliminate the muon DIS-induced neutral meson
production at the back of the dump (see Section V.3). However, we find that
simple geometric acceptance and photon energy selections reduce the estimated
$\pi^{0}$ yield in Eq. 16 by a factor of $>200$. Similar conclusions hold for
DIS-induced $K_{L}$ backgrounds: requiring exactly two photons with energy
above a GeV in the geometric acceptance reduces this by a factor $>1000$, even
before cutting on the final state muon or photon separation.

We have also explored a diphoton invariant mass cut to reduce backgrounds
coming from $\pi^{0}$ decays. While excluding events with
$m_{\gamma\gamma}\approx m_{\pi}$ does decrease these by a factor of at most
$\sim 2$, the poor diphoton mass resolution discussed in Section II.3 means
that the signal yield is also negatively affected since DarkQuest’s most
important sensitivity lies in the window $m_{\mu}\lesssim m_{S}\lesssim
2m_{\mu}$, i.e., near the pion mass.

##  VII Conclusions

The proposed DarkQuest experiment offers a powerful test of scalars coupled to
muons. In this paper, we studied di-photon signatures at DarkQuest and showed
that DarkQuest will be able to probe new regions of parameter space that
address the $(g-2)_{\mu}$ anomaly already with $10^{18}$ POT, for $m_{S}\leq
2m_{\mu}$. The large number of secondary mesons and muons produced in the
proton-dump collisions open multiple production channels for new light states
whose displaced decays can be observed with the downstream detectors. We
identified new production channels for scalars coupled to muons and background
processes relevant for the displaced diphoton signatures at DarkQuest. We
carefully modeled production and detection of signal and background events and
outlined several strategies for achieving the projected sensitivities in Fig.
7. Other searches (prompt $S\to\mu\mu)$ using SpinQuest will be able to probe
higher masses Forbes:2022bvo .

## Acknowledgments

We thank Cristina Mantilla Suarez, Christian Herwig, Yongbin Feng, Nhan Tran,
Dean Robinson, Shirley Li, Wendy Taylor and Claudia Cornella for valuable
discussions. This research was enabled in part by support provided by the BC
DRI Group, Compute Ontario and the Digital Research Alliance of Canada
(alliancecan.ca). NB was supported in part by Natural Sciences and Engineering
Research Council of Canada (NSERC). The research of SG and NH is supported in
part by the U.S. Department of Energy grant number DE-SC0023093 and by the
U.S. Department of Energy grant number DE-SC0010107. Part of this work was
performed at the Aspen Center for Physics, which is supported by National
Science Foundation grant PHY-1607611.

##  Appendix A Scalar-Photon Coupling

In this section we calculate the form factor associated with the amplitude
$S\to\gamma\gamma$ which gives rise to the photon coupling in Eq. 4. This
amplitude provides the main $S$ decay channel (for $m_{S}<2m_{\mu}$), and an
$S$ production mechanism through secondary photons in the beam dump. The
former involves two on-shell photons, while for the latter at least one photon
is off-shell. In general, this amplitude has both UV and IR contributions
(from heavy matter that generates the defining $S\bar{\mu}\mu$ interaction and
from the muon loop, respectively). The one loop amplitude with a single on-
shell photon, generated by a generic fermion with mass $M$ and coupling
$g_{f}$ to $S$ is

$\mathcal{M}(S(p_{S})\to\gamma(p_{1})\gamma^{*}(p_{2}))\equiv\frac{g_{f}\alpha
Q^{2}}{2\pi
M}f_{1/2}\left(\frac{4M^{2}}{m_{S}^{2}},\frac{p_{2}^{2}}{4M^{2}}\right)(p_{1}^{\mu}p_{2}^{\nu}-p_{1}\cdot
p_{2}\eta^{\mu\nu}),$ (18)

where the form factor is

$f_{1/2}(\tau,\rho)=\frac{\tau\left((L_{1}^{2}-L_{2}^{2})(\rho\tau+\tau-1)+4\rho\tau\left(L_{1}\sqrt{\frac{1-\tau}{\tau^{2}}}\tau-1\right)-4L_{2}\sqrt{(\rho-1)\rho}\tau+4\right)}{2(\rho\tau-1)^{2}},$
(19)

with

$\tau=\frac{4M^{2}}{p_{S}^{2}},~{}p_{1}^{2}=0,~{}\rho=\frac{p_{2}^{2}}{4M^{2}},$
(20)

and

$L_{1}=\ln\left(2\sqrt{\frac{1-\tau}{\tau^{2}}}-\frac{2}{\tau}+1\right),\;\;\;L_{2}=\ln\left(-2\rho+2\sqrt{(\rho-1)\rho}+1\right).$
(21)

This form factor reduces to the familiar expressions found for the fermionic
contribution to $h\to\gamma\gamma$ Gunion:1989we , when both photons are on-
shell. For muons $g_{f}=g_{S}$ and $M=m_{\mu}$, but we can also use this
result to estimate the contribution of heavy fermions $M^{2}\gg
m_{S}^{2},\;p_{2}^{2}$ to the diphoton coupling
($g_{S\gamma\gamma}^{(\mathrm{UV})}$ in Eq. 5):

$\lim_{M\to\infty}f_{1/2}\left(\frac{4M^{2}}{m_{S}^{2}},\frac{p_{2}^{2}}{4M^{2}}\right)=\frac{4}{3},$
(22)

so that

$g_{S\gamma\gamma}^{(\mathrm{UV})}\sim\frac{2g_{f}\alpha Q^{2}}{3\pi M}.$ (23)

Note that generally $g_{f}\neq g_{S}$ (see next section) so the UV
contribution to $S\gamma\gamma$ is not necessarily negligible compared to the
muon one as we discuss in the next section.

##  Appendix B An Ultraviolet Completion

In this section we explore a simple ultraviolet completion for the effective
operator in Eq. 7, $-\frac{S}{M}H^{\dagger}L\mathbf{c}_{S}\ell_{R}^{c},$ that
generates the muon-philic interaction. This will allow to discuss how the low
energy couplings are related to the underlying UV parameters, and to mention
additional constraints on realistic scenarios. A similar scenario was
considered in the Supplementary Materials of Ref. Krnjaic:2019rsv .

The operator in Eq. 7 can be generated at tree-level by a heavy vector-like
lepton (VLL) pair. The simplest possibility is a VLL pair of $SU(2)_{L}$
singlets $\psi$ and $\psi^{c}$, which have the same quantum numbers as the
right-handed leptons $\ell^{c}_{R}$ in the SM (a completion with $SU(2)_{L}$
doublets is discussed in Ref. Batell:2017kty ). The allowed renormalizable
terms (in two-component notation) are

$\mathscr{L}\supset-M\psi\psi^{c}-H^{\dagger}L^{T}y_{\psi}\psi^{c}-y_{S}^{\prime}S\psi\psi^{c}-S\psi
y_{S}^{T}\ell_{R}^{c}+\mathrm{h.c.},$ (24)

where $L$ and $\ell_{R}^{c}$ are SM leptons (with flavor index suppressed),
and the transposes are in flavor space. Since we consider a single generation
of VLL, $M$ and $y_{S}^{\prime}$ are numbers, while $y_{\psi}$ and $y_{S}$ are
vectors (with each entry corresponding to a SM fermion generation). In the
last two terms, we assumed that $\psi$ and $\psi^{c}$ have a Yukawa
interaction with a scalar mediator $S$.101010 The mass-mixing term
$\mu_{\ell\psi}\ell_{R}^{c}\psi$ is allowed by gauge charges, but it can be
rotated away by performing an orthogonal rotation on
$(\ell_{R}^{c}\;\;\psi^{c})$; this rotation results in redefinitions of the
(unknown) parameters $M$, $y_{S}^{\prime}$, $y_{S}$, $y_{\psi}$ and
$y_{\ell}$, which can then be relabelled. Integrating out $\psi$, and
$\psi^{c}$ yields the effective operator

$\mathcal{L}\supset\frac{S}{M}H^{\dagger}L^{T}y_{\psi}y_{S}^{T}\ell_{R}^{c}+\mathrm{h.c.}\to-\frac{v}{\sqrt{2}M}SL_{L}\mathbf{c}_{S}\ell_{R}^{c}+\mathrm{h.c.},$
(25)

where in the second step we transformed to the SM-lepton mass basis, which
gives $\mathbf{c}_{S}\equiv-V_{L}^{T}y_{\psi}y^{T}_{S}V_{R}$ ($V_{L,R}$ are
the left- and right-handed lepton unitary rotations that diagonalize their
mass matrix). A muon-only interaction therefore requires that
$\mathbf{c}_{S}\propto\operatorname{diag}(0,1,0)$ as described in the
introduction, necessitating a non-trivial flavor model. If this is realized
then

$g_{S}=\frac{y_{S}y_{\psi}v}{\sqrt{2}M}.$ (26)

There are several theoretical constraints on $g_{S}$ that follow from
naturalness of the $S$ mass (i.e., the desire for quantum corrections to
$m_{S}$ to not exceed $m_{S}$ itself). These can be derived either directly
from the UV completion above or from the EFT interaction as in Ref.
Batell:2017kty . The two most relevant quantum corrections to $m_{S}$ arise
from 1) a two loop contribution from a Higgs-left handed lepton-right handed
lepton loop; 2) from the one loop radiatively-generated $|H|^{2}S^{2}$
interaction. Requiring that the two corrections do not induce a mass shift
larger than $m_{S}$ leads to a theoretical upper bound on $g_{S}$:

$g_{S}\lesssim\frac{4\pi
v}{\sqrt{2}M}\min\left(\frac{\sqrt{2}m_{S}}{v},4\pi\frac{m_{S}}{M}\right),$
(27)

where the first (second) term comes from the two (one) loop correction. The
one loop effect dominates for $M\lesssim 4\pi v/\sqrt{2}\simeq 2$ TeV. At the
same time, $M$ cannot be too small in order to satisfy Eq. 27 in the region of
interest (see Batell:2017kty ).

We can also estimate the UV contribution to the $S\gamma\gamma$ coupling, by
using Eq. 23. At one loop, only the $y_{S}^{\prime}$ coupling from Eq. 24
contributes, yielding

$g_{S\gamma\gamma}^{(\mathrm{UV})}\sim\frac{2y_{S}^{\prime}\alpha Q^{2}}{3\pi
M}.$ (28)

We see that $g_{S\gamma\gamma}^{(\mathrm{UV})}$ is not necessarily correlated
with $g_{S}$, and it can therefore be comparable to
$g_{S\gamma\gamma}^{(\mathrm{IR})}$ (generated by a muon loop) for
$y_{S}^{\prime}\gtrsim g_{S}(M/m_{\mu})$.111111This, however, would require a
large hierarchy between dimensionless parameters.

Finally let us briefly discuss high energy collider constraints on VLL. Recent
theoretical studies of these models in the context of the LHC and future
colliders include Refs. Dermisek:2014qca ; Kumar:2015tna ; Bhattiprolu:2019vdu
. The production rate (and the resulting bounds) depend strongly on the VLL
$SU(2)_{L}$ representation Kumar:2015tna , with singlets usually being the
least constrained. In the singlet case, the main production channel is
$pp\to\gamma,Z^{*}\to\psi^{+}\psi^{-}$. The heavy leptons decay to SM
particles via the interactions of Eq. 24. Our model differs slightly from the
minimal VLL scenarios studied in these works, as it features new production
and decay mechanisms involving $S$. For example, in addition to the “standard”
decays $\psi\to W\nu,\;Z\ell$, the VLL can decay $\psi\to S\ell$. Even more
interestingly, interference of SM and $S$-mediated channels can significantly
weaken existing constraints Egana-Ugrinovic:2018roi . The standard and
$S$-mediated effects are parametrized by independent couplings, making it
difficult to directly translate existing bounds on VLL models. It would
therefore be interesting to carry out a dedicated study for the VLL+$S$ model
at the LHC.

With these important caveats in mind, we note some constraints on the minimal
VLL model (without $S$) that may be relevant for the full scenario: LEP ruled
out many kinds of charged particles with mass below $\sim 100$ GeV ADLO_SUSY ;
L3:2001xsz . The best current LHC limits on the singlet VLL scenario are from
ATLAS and CMS ATLAS:2015qoy ; CMS:2022nty , with $M$ excluded in the ranges
$114-176$ GeV ($\psi\to Z+e/\mu$) or $125-150$ GeV ($\psi\to Z+\tau$).

##  Appendix C Muon Transport

In this Appendix, we describe our simulation of muon transport in the dump. We
developed our own propagation code instead of using off-the shelf tools such
as GEANT in order to have a simple way to incorporate rare processes relevant
for signal and background events (e.g., $S$ bremsstrahlung and muon deep
inelastic scattering). Simulations of these processes are performed using
other codes (e.g., Madgraph and Pythia), so having minimal muon code enables a
simple interface between them.

Muons produced near the front of the beam dump can propagate a significant
distance before undergoing a dark bremsstrahlung reaction. The energy and
angular distributions of the resulting muon-philic scalars are sensitive to
the details of muon transport through the dump. Muon propagation is also
relevant for predicting the rates of some of the backgrounds, such as DIS-
induced meson production at the very back of the dump (see Section V). Below
we describe our treatment of muon production and propagation in the magnetic
field, including energy loss and multiple scattering.

We generate an initial muon sample using Pythia to model $pp$ collisions with
incoming $120$ GeV proton beam, assuming the production takes place in the
first nucleon interaction length. The primary proton flux is exponentially
diminished beyond the first interaction length, and the surviving protons lose
energy; neglecting these protons is therefore a conservative simplifying
assumption.

The muons come dominantly from charged pion decays. We require that the pions
decay within the first pion interaction length, $\lambda_{\pi\text{ int}}$.
The fraction of decayed pions grows linearly with distance (in the long life-
time regime), while their flux is exponentially suppressed by
$\exp(-z/\lambda_{\pi\text{ int}})$. As a result the total muon flux is well
approximated by the first interaction length alone. This assumption allows us
to neglect pion transport in the dump. The energy and transverse momentum
distributions of pions and muons obtained this way are shown in Fig. 2.

The secondary muons can undergo bremsstrahlung in the dump. Interactions close
to the end of the dump are particularly interesting since they enable
sensitivity to shorter scalar lifetimes. We propagate the muons through the
FMAG in small steps, including the following effects:

  * a.

The 1.8 T magnetic field of the FMAG;

  * b.

Ionization energy loss;

  * c.

Multiple Coulomb scattering.

The magnetic field is included by solving the relativistic equation of motion
in each step. Energy loss and multiple scattering are stochastic processes, so
they are applied at the end of each step. We describe our treatment of these
processes below.

While $dE/dx$ is frequently used to estimate the energy losses of particles
propagating through a medium, it is important to note that $dE/dx$ is an
_average_ loss, appropriate for describing very large ensembles of particles.
This average is sensitive to rare, high loss collisions. This does not mean
that it is an accurate representation of the typical losses of individual
particles Workman:2022ynf . Indeed, the definition of $dE/dx$ includes large
losses due to rare high-energy transfer events, skewing the mean loss to
larger values. For example, the most likely loss for a 5 GeV muon traversing 1
cm of iron is $10.6\;\mathrm{MeV}$, while $dE/dx$ is $14.6\;\mathrm{MeV}$
Workman:2022ynf . In practice, this means that the losses of individual
particles follow a distribution where the peak of the distribution (the most
probable loss) is different from the average. In our case, energy loss
fluctuations can be modeled by the Landau distribution (see Refs.
Talman:1978tc ; Bak:1987cz ; Bichsel:1988if ; mcparland2016medical  for
helpful reviews):

$f(\Delta)=\frac{1}{2\pi i\xi}\int_{c-i\infty}^{c+i\infty}ds\exp{(\lambda
s+s\ln s)},\;\;\;\lambda=\frac{\Delta-\Delta_{p}}{\xi}-0.22278,$ (29)

where $c$ is a positive number, $\Delta$ is the energy loss, $\xi$ and
$\Delta_{p}$ are given by Workman:2022ynf :

$\Delta_{p}=\xi\left[\ln\frac{2m_{e}\beta^{2}\gamma^{2}}{I}+\ln\frac{\xi}{I}+0.2-\delta\right],\;\;\xi=0.153538\;\mathrm{MeV}\;\mathrm{cm^{2}}\;\mathrm{mol}^{-1}\frac{Z}{A}\frac{t}{\beta^{2}},$
(30)

where $t$ is the thickness in g/cm2 of the material traversed, $\beta$ the
particle speed, and $A$ and $Z$ are the atomic mass and number, respectively.
The material dependent parameters $I$ and $\delta$ are the mean excitation
energy and density correction which are evaluated using the tabulated data of
Ref. Workman:2022ynf . The Landau distribution is an appropriate description
of energy loss fluctuations when the losses are small compared to the maximum
kinematically allowed loss and the initial kinetic energy
($\xi/E_{\mathrm{max}}\ll 1$ with $E_{\mathrm{max}}$ given in Workman:2022ynf
), but large compared to typical atomic binding energies ($\xi/I\gg 1$)
mcparland2016medical . Both conditions are easily satisfied for muon energies
of $\sim$ GeV and above, and for propagation distances of a few cm. Eq. 29 is
heavy-tailed with an undefined mean, implying that particles are allowed to
lose arbitrarily large amounts of energy. This is of course unphysical since
there is a maximum energy loss dictated by kinematics. Thus, the distribution
should be truncated, i.e., $\lambda\leq\lambda_{\mathrm{max}}$ for some
$\lambda_{\mathrm{max}}$. We implement such a truncation by requiring that the
mean energy loss of the truncated distribution is equal to $dE/dx$. Sampling
from this truncated distribution is achieved using a modified version of the
landaupy package. We checked that our implementation of energy losses
reproduces those of GEANT for muons with initial energies from 1 to 100 GeV
propagating through iron targets of 1 to 500 cm, covering the parameter ranges
relevant for DarkQuest.

Multiple scattering is modeled by applying a transverse momentum kick drawn
from a Gaussian approximation to the Moliere distribution of scattering
angles. The variance $\sigma^{2}$ of this distribution is given in Sec. 7 of
Ref. Lynch:1990sq :

$\sigma^{2}=\frac{\chi_{c}^{2}}{1+F}\left[\frac{1+v}{v}\log(1+v)-1\right],\;\;\chi_{c}^{2}=\frac{4\pi
n_{A}t\alpha^{2}Z(Z+1)}{(p\beta)^{2}},\;\;v=\frac{\Omega}{2(1-F)},\;\;F=0.98$
(31)

where $n_{A}$ is the atomic number density, $p$ is the particle momentum and
$\Omega$ is the mean number of scatters

$\Omega=\frac{7800~{}t~{}(Z+1)Z^{1/3}}{\beta^{2}A\left[1+3.34\left(\frac{Z\alpha}{\beta}\right)^{2}\right]}$
(32)

with $\beta$, $t$, $Z$, $A$ defined above. The Gaussian approximation provides
an adequate description of multiple Coulomb scattering when $\Omega\gg 1$;
this is the case for $t\gtrsim 10^{-3}X_{0}$, $\beta=1$ and arbitrary $Z$,
where $X_{0}$ is the radiation length Lynch:1990sq .

We chose a 1 cm step size to ensure the validity of the various approximations
for the treatments of ionization and multiple scattering. Specifically, 1 cm
is large enough such that the number of Coulomb collisions is large and the
Gaussian approximation is valid, but small enough that the losses do not
significantly affect the muon trajectory in a single step. A larger step size
is also desirable for computational speed. We have verified that this
algorithm reproduces the GEANT4-transported muon spectrum used in Ref.
Forbes:2022bvo .121212We thank Cristina Mantilla Suarez and Christian Herwig
for extensive discussions regarding this check. It also reproduces the spectra
of monoenergetic muons propagated through $5$ m of iron (this cross-check
eliminates the convolution of propagation effects with the primary muon
spectrum). Note that while Forbes:2022bvo  focused on a dimuon resonance
search with a monoenergetic muon beam, here we take into account the full muon
spectrum from beam-target interactions and its propagation through the dump.

The muons are thus transported until they either leave the dump or reach an
energy below a GeV (in our reach estimates we always require the signal
photons to each have an energy of at least 500 MeV, which implies a minimum
muon energy of $\mathcal{O}(1)$ GeV). The magnetic field and energy loss are
by far the most important effects: The $p_{T}$ imparted by the $B$ field is
$\sim 2.6\;\mathrm{GeV}$ for relativistic muons that traverse the entire FMAG;
in practice the $p_{T}$ fluctuates around this value due to multiple
scattering and energy losses; it can also be substantially smaller if the muon
ranges out inside the dump. The most probable energy loss ($\Delta_{p}$ in Eq.
30) is $\sim 7\;\mathrm{GeV}$ after traversing the FMAG.

The phase-space trajectories obtained using the above procedure are used to
sample the momenta and interaction vertices at the interaction time. The
probability distribution of the muon path length $\ell$ through the FMAG at
interaction is given by (in the small coupling limit)

$\frac{dP}{d\ell}\approx n_{A}\sigma_{\rm{Brem.}}(E_{\mu}(\ell)),$ (33)

where $n_{A}$ is the number density of target nuclei,
$\sigma_{\mathrm{Brem.}}$ is the $\mu N\to\mu N+S$ cross-section, and
$E_{\mu}(\ell)$ is obtained from the muon trajectory. This distribution is
sampled using the inverse transform method. In Fig. 11 we show the energy and
transverse momentum spectra of the initial muons, the muons at their
interaction positions and the muons that traverse the entire FMAG without
interacting. We clearly see the effects of transport, as the muons typically
have a lower energy and larger $p_{T}$ the more material they traverse. In our
simulation, we force every generated muon to interact within the FMAG. In this
approach, each event must be weighted by the probability of the muon do so:

$n_{A}\int_{\mathrm{FMAG}}\sigma_{\rm{Brem.}}(E_{\mu}(\ell))d\ell.$ (34)

The geometric acceptance-weighted sum of these probabilities (normalized to
the total number of simulated $pp$ collisions) gives the total $S$ yield (see
Eq. 9).

Figure 11: Muon spectra at creation (blue lines), at the $\mu N\to\mu N+S$
interaction (orange lines), and at FMAG exit (gray lines). The initial spectra
are generated using Pythia (see Fig. 2) and transported as described in the
text. The transport includes the effects of energy loss, multiple scattering,
and the magnetic field which are evident in the left and right panels,
respectively. Note that blue and orange histograms sum to $0.0028$ per POT,
which is smaller than the total muon yield of $0.01$ per POT in Table 2 and
the corresponding histogram in Fig. 2; this is because here we impose
$E_{\mu}>1\;\mathrm{GeV}$. The gray histogram contains a smaller number of
muons still as not all muons are able to penetrate the entire FMAG.

##  Appendix D Rare Meson Decays

The full result for the three-body decay of the kaon to produce a scalar is
given by,131313 The Mathematica version of this expression can be found here.

$\displaystyle\Gamma\left(K\rightarrow\mu\nu
S\right)=\frac{1}{(2\pi)^{3}}\frac{(G_{f}f_{K}V_{\textrm{us}})^{2}}{4m_{K}}g_{S}^{2}\quad\Bigg{\\{}\Bigg{.}$
$\displaystyle\frac{1}{48m_{K}^{2}}\left(2m_{K}^{4}+26m_{K}^{2}m_{S}^{2}+2m_{S}^{4}-m_{\mu}^{2}\left(115m_{K}^{2}+61m_{S}^{2}\right)+119m_{\mu}^{4}\right)\lambda^{1/2}\left(m_{K}^{2},m_{\mu}^{2},m_{S}^{2}\right)-$
$\displaystyle\frac{1}{8m_{K}^{2}m_{\mu}^{2}}m_{S}^{2}\left(m_{K}^{2}-m_{\mu}^{2}\right)\left(2m_{K}^{2}m_{\mu}^{2}+m_{K}^{2}m_{S}^{2}-18m_{\mu}^{4}+3m_{\mu}^{2}m_{S}^{2}\right)g\left(m_{\mu}^{2}/m_{S}^{2}\right)\tan^{-1}\left[\frac{g\left(m_{\mu}^{2}/m_{S}^{2}\right)\,\lambda^{1/2}\left(m_{K}^{2},m_{\mu}^{2},m_{S}^{2}\right)}{m_{K}^{2}-m_{S}^{2}+3m_{\mu}^{2}}\right]-$
$\displaystyle\frac{1}{8m_{\mu}^{2}}m_{K}^{2}\left(m_{S}^{2}-m_{\mu}^{2}\right)^{2}\log\left(\frac{m_{K}^{2}\left(m_{S}^{2}+m_{\mu}^{2}\right)+\left(m_{S}^{2}-m_{\mu}^{2}\right)\left[-m_{S}^{2}+m_{\mu}^{2}+\lambda^{1/2}\left(m_{K}^{2},m_{\mu}^{2},m_{S}^{2}\right)\right]}{m_{K}^{2}\left[m_{K}^{2}-m_{S}^{2}-m_{\mu}^{2}-\lambda^{1/2}\left(m_{K}^{2},m_{\mu}^{2},m_{S}^{2}\right)\right]}\right)+$
$\displaystyle\frac{1}{8m_{K}^{2}m_{\mu}^{2}}\Big{[}\Big{.}\qquad
m_{K}^{4}\left(m_{S}^{4}-4m_{\mu}^{2}m_{S}^{2}+8m_{\mu}^{4}\right)-\,m_{K}^{2}m_{\mu}^{2}\left(2m_{S}^{4}-24m_{S}^{2}m_{\mu}^{2}+16m_{\mu}^{4}\right)+\,m_{\mu}^{4}\left(3m_{S}^{4}-36m_{\mu}^{2}m_{S}^{2}+8m_{\mu}^{4}\right)\quad\Big{.}\Big{]}\times$
$\displaystyle\log\left(\frac{2m_{S}m_{\mu}}{m_{K}^{2}-m_{S}^{2}-m_{\mu}^{2}-\lambda^{1/2}\left(m_{K}^{2},m_{\mu}^{2},m_{S}^{2}\right)}\right)\Bigg{.}\Bigg{\\}},$
(35)

where

$\displaystyle g\left(m_{\mu}^{2}/m_{S}^{2}\right)$
$\displaystyle\equiv\sqrt{4m_{\mu}^{2}/m_{S}^{2}-1},$ (36)
$\displaystyle\lambda^{1/2}\left(m_{K}^{2},m_{\mu}^{2},m_{S}^{2}\right)$
$\displaystyle\equiv\sqrt{m_{K}^{4}+m_{S}^{4}+m_{\mu}^{4}-2m_{K}^{2}m_{S}^{2}-2m_{K}^{2}m_{\mu}^{2}-2m_{S}^{2}m_{\mu}^{2}},$
(37)

where $f_{K}=156\;\mathrm{MeV}$ is the kaon decay constant. A similar
expression holds for the corresponding pion three-body decay.

## References

  * (1) S. Gori et al., Dark Sector Physics at High-Intensity Experiments, arXiv:2209.04671.
  * (2) J. Cooley et al., Report of the Topical Group on Particle Dark Matter for Snowmass 2021, arXiv:2209.07426.
  * (3) P. Coloma, L. W. Koerner, I. M. Shoemaker, and J. Yu, Neutrino Frontier Topical Group Report (NF03): Physics Beyond the Standard Model, arXiv:2209.10362.
  * (4) S. Gardner, R. J. Holt, and A. S. Tadepalli, New Prospects in Fixed Target Searches for Dark Forces with the SeaQuest Experiment at Fermilab, Phys. Rev. D 93 (2016), no. 11 115015, [arXiv:1509.00050].
  * (5) A. Berlin, S. Gori, P. Schuster, and N. Toro, Dark Sectors at the Fermilab SeaQuest Experiment, Phys. Rev. D 98 (2018), no. 3 035011, [arXiv:1804.00661].
  * (6) C. Ahdida et al., Post-LS3 Experimental Options in ECN3, arXiv:2310.17726.
  * (7) A. Apyan et al., DarkQuest: A dark sector upgrade to SpinQuest at the 120 GeV Fermilab Main Injector, in 2022 Snowmass Summer Study, 3, 2022.  arXiv:2203.08322.
  * (8) Muon g-2 Collaboration, G. W. Bennett et al., Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL, Phys. Rev. D 73 (2006) 072003, [hep-ex/0602035].
  * (9) Muon g-2 Collaboration, B. Abi et al., Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm, Phys. Rev. Lett. 126 (2021), no. 14 141801, [arXiv:2104.03281].
  * (10) Muon g-2 Collaboration, D. P. Aguillard et al., Measurement of the Positive Muon Anomalous Magnetic Moment to 0.20 ppm, Phys. Rev. Lett. 131 (2023), no. 16 161802, [arXiv:2308.06230].
  * (11) T. Aoyama et al., The anomalous magnetic moment of the muon in the Standard Model, Phys. Rept. 887 (2020) 1–166, [arXiv:2006.04822].
  * (12) S. Borsanyi et al., Leading hadronic contribution to the muon magnetic moment from lattice QCD, Nature 593 (2021), no. 7857 51–55, [arXiv:2002.12347].
  * (13) T. Blum et al., An update of Euclidean windows of the hadronic vacuum polarization, arXiv:2301.08696.
  * (14) A. Bazavov et al., Light-quark connected intermediate-window contributions to the muon $g-2$ hadronic vacuum polarization from lattice QCD, arXiv:2301.08274.
  * (15) R. Capdevilla, D. Curtin, Y. Kahn, and G. Krnjaic, Systematically testing singlet models for (g $-$ 2)?, JHEP 04 (2022) 129, [arXiv:2112.08377].
  * (16) NA48/2 Collaboration, J. R. Batley et al., Search for the dark photon in $\pi^{0}$ decays, Phys. Lett. B 746 (2015) 178–185, [arXiv:1504.00607].
  * (17) BaBar Collaboration, J. P. Lees et al., Search for a Dark Leptophilic Scalar in $e^{+}e^{-}$ Collisions, Phys. Rev. Lett. 125 (2020), no. 18 181801, [arXiv:2005.01885].
  * (18) BaBar Collaboration, J. P. Lees et al., Search for Invisible Decays of a Dark Photon Produced in ${e}^{+}{e}^{-}$ Collisions at BaBar, Phys. Rev. Lett. 119 (2017), no. 13 131804, [arXiv:1702.03327].
  * (19) C.-Y. Chen, M. Pospelov, and Y.-M. Zhong, Muon Beam Experiments to Probe the Dark Sector, Phys. Rev. D 95 (2017), no. 11 115005, [arXiv:1701.07437].
  * (20) A. Berlin, N. Blinov, S. Gori, P. Schuster, and N. Toro, Cosmology and Accelerator Tests of Strongly Interacting Dark Matter, Phys. Rev. D 97 (2018), no. 5 055033, [arXiv:1801.05805].
  * (21) B. Döbrich, J. Jaeckel, and T. Spadaro, Light in the beam dump - ALP production from decay photons in proton beam-dumps, JHEP 05 (2019) 213, [arXiv:1904.02091]. [Erratum: JHEP 10, 046 (2020)].
  * (22) Y.-D. Tsai, P. deNiverville, and M. X. Liu, Dark Photon and Muon $g-2$ Inspired Inelastic Dark Matter Models at the High-Energy Intensity Frontier, Phys. Rev. Lett. 126 (2021), no. 18 181801, [arXiv:1908.07525].
  * (23) B. Batell, J. A. Evans, S. Gori, and M. Rai, Dark Scalars and Heavy Neutral Leptons at DarkQuest, JHEP 05 (2021) 049, [arXiv:2008.08108].
  * (24) N. Blinov, E. Kowalczyk, and M. Wynne, Axion-like particle searches at DarkQuest, JHEP 02 (2022) 036, [arXiv:2112.09814].
  * (25) D. Forbes, C. Herwig, Y. Kahn, G. Krnjaic, C. Mantilla Suarez, N. Tran, and A. Whitbeck, New searches for muonphilic particles at proton beam dump spectrometers, Phys. Rev. D 107 (2023), no. 11 116026, [arXiv:2212.00033].
  * (26) A. Ariga, R. Balkin, I. Galon, E. Kajomovitz, and Y. Soreq, Hunting muonic forces at emulsion detectors, Phys. Rev. D 109 (2024), no. 3 035003, [arXiv:2305.03102].
  * (27) S. N. Gninenko, N. V. Krasnikov, and V. A. Matveev, Muon g-2 and searches for a new leptophobic sub-GeV dark boson in a missing-energy experiment at CERN, Phys. Rev. D 91 (2015) 095015, [arXiv:1412.1400].
  * (28) C.-Y. Chen, J. Kozaczuk, and Y.-M. Zhong, Exploring leptophilic dark matter with NA64-$\mu$, JHEP 10 (2018) 154, [arXiv:1807.03790].
  * (29) Y. Kahn, G. Krnjaic, N. Tran, and A. Whitbeck, M3: a new muon missing momentum experiment to probe (g $-$ 2)μ and dark matter at Fermilab, JHEP 09 (2018) 153, [arXiv:1804.03144].
  * (30) L. Marsicano, M. Battaglieri, A. Celentano, R. De Vita, and Y.-M. Zhong, Probing Leptophilic Dark Sectors at Electron Beam-Dump Facilities, Phys. Rev. D 98 (2018), no. 11 115022, [arXiv:1812.03829].
  * (31) G. Krnjaic, G. Marques-Tavares, D. Redigolo, and K. Tobioka, Probing Muonphilic Force Carriers and Dark Matter at Kaon Factories, Phys. Rev. Lett. 124 (2020), no. 4 041802, [arXiv:1902.07715].
  * (32) C. Rella, B. Döbrich, and T.-T. Yu, Searching for muonphilic dark sectors with proton beams, Phys. Rev. D 106 (2022), no. 3 035023, [arXiv:2205.09870].
  * (33) C. Cesarotti, Y. Kahn, G. Krnjaic, D. Rocha, and J. Spitz, New $\mu$ Forces From $\nu_{\mu}$ Sources, arXiv:2311.10829.
  * (34) SeaQuest Collaboration, C. A. Aidala et al., The SeaQuest Spectrometer at Fermilab, Nucl. Instrum. Meth. A 930 (2019) 49–63, [arXiv:1706.09990].
  * (35) S. Gori et al., ACE Science Workshop Report, arXiv:2403.02422.
  * (36) C. Bierlich et al., A comprehensive guide to the physics and usage of PYTHIA 8.3, arXiv:2203.11601.
  * (37) B. Batell, N. Lange, D. McKeen, M. Pospelov, and A. Ritz, Muon anomalous magnetic moment through the leptonic Higgs portal, Phys. Rev. D 95 (2017), no. 7 075003, [arXiv:1606.04943].
  * (38) B. Batell, A. Freitas, A. Ismail, and D. Mckeen, Flavor-specific scalar mediators, Phys. Rev. D 98 (2018), no. 5 055026, [arXiv:1712.10022].
  * (39) Y.-S. Tsai, AXION BREMSSTRAHLUNG BY AN ELECTRON BEAM, Phys. Rev. D 34 (1986) 1326.
  * (40) Y.-S. Tsai, PRODUCTION OF NEUTRAL BOSONS BY AN ELECTRON BEAM, Phys. Rev. D 40 (1989) 760.
  * (41) Y.-S. Liu, D. McKeen, and G. A. Miller, Validity of the Weizsäcker-Williams approximation and the analysis of beam dump experiments: Production of a new scalar boson, Phys. Rev. D 95 (2017), no. 3 036010, [arXiv:1609.06781].
  * (42) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, [arXiv:1405.0301].
  * (43) G. P. Lepage, A New Algorithm for Adaptive Multidimensional Integration, J. Comput. Phys. 27 (1978) 192.
  * (44) G. P. Lepage, Adaptive multidimensional integration: VEGAS enhanced, J. Comput. Phys. 439 (2021) 110386, [arXiv:2009.05112].
  * (45) A. Celentano, L. Darmé, L. Marsicano, and E. Nardi, New production channels for light dark matter in hadronic showers, Phys. Rev. D 102 (2020), no. 7 075026, [arXiv:2006.09419].
  * (46) GEANT4 Collaboration, S. Agostinelli et al., GEANT4–a simulation toolkit, Nucl. Instrum. Meth. A 506 (2003) 250–303.
  * (47) J. Allison et al., Geant4 developments and applications, IEEE Trans. Nucl. Sci. 53 (2006) 270.
  * (48) J. Allison et al., Recent developments in Geant4, Nucl. Instrum. Meth. A 835 (2016) 186–225.
  * (49) CHARM Collaboration, F. Bergsma et al., A Search for Decays of Heavy Neutrinos, Phys. Lett. B 128 (1983) 361.
  * (50) CHARM Collaboration, F. Bergsma et al., A Search for Decays of Heavy Neutrinos in the Mass Range 0.5-GeV to 2.8-GeV, Phys. Lett. B 166 (1986) 473–478.
  * (51) S. N. Gninenko, Constraints on sub-GeV hidden sector gauge bosons from a search for heavy neutrino decays, Phys. Lett. B 713 (2012) 244–248, [arXiv:1204.3583].
  * (52) PHENIX Collaboration, L. Aphecetche et al., PHENIX calorimeter, Nucl. Instrum. Meth. A 499 (2003) 521–536.
  * (53) J. F. Gunion, H. E. Haber, G. L. Kane, and S. Dawson, The Higgs Hunter’s Guide, vol. 80.  2000\.
  * (54) R. Dermisek, J. P. Hall, E. Lunghi, and S. Shin, Limits on Vectorlike Leptons from Searches for Anomalous Production of Multi-Lepton Events, JHEP 12 (2014) 013, [arXiv:1408.3123].
  * (55) N. Kumar and S. P. Martin, Vectorlike Leptons at the Large Hadron Collider, Phys. Rev. D 92 (2015), no. 11 115018, [arXiv:1510.03456].
  * (56) P. N. Bhattiprolu and S. P. Martin, Prospects for vectorlike leptons at future proton-proton colliders, Phys. Rev. D 100 (2019), no. 1 015033, [arXiv:1905.00498].
  * (57) D. Egana-Ugrinovic, M. Low, and J. T. Ruderman, Charged Fermions Below 100 GeV, JHEP 05 (2018) 012, [arXiv:1801.05432].
  * (58) LEPSUSYWG, Combined lep chargino results, up to 208 gev for large m0, 2001.  https://lepsusy.web.cern.ch/lepsusy/www/inos_moriond01/charginos_pub.html.
  * (59) L3 Collaboration, P. Achard et al., Search for heavy neutral and charged leptons in $e^{+}e^{-}$ annihilation at LEP, Phys. Lett. B 517 (2001) 75–85, [hep-ex/0107015].
  * (60) ATLAS Collaboration, G. Aad et al., Search for heavy lepton resonances decaying to a $Z$ boson and a lepton in $pp$ collisions at $\sqrt{s}=8$ TeV with the ATLAS detector, JHEP 09 (2015) 108, [arXiv:1506.01291].
  * (61) CMS Collaboration, A. Tumasyan et al., Inclusive nonresonant multilepton probes of new phenomena at $\sqrt{s}$=13 TeV, Phys. Rev. D 105 (2022), no. 11 112007, [arXiv:2202.08676].
  * (62) Particle Data Group Collaboration, R. L. Workman and Others, Review of Particle Physics, PTEP 2022 (2022) 083C01.
  * (63) R. Talman, On the Statistics of Particle Identification Using Ionization, Nucl. Instrum. Meth. 159 (1979) 189.
  * (64) J. F. Bak, A. Burenkov, J. B. B. Petersen, E. Uggerhoj, S. P. Moller, and P. Siffert, Large Departures From Landau Distributions for High-energy Particles Traversing Thin Si and Ge Targets, Nucl. Phys. B 288 (1987) 681–716.
  * (65) H. Bichsel, Straggling in Thin Silicon Detectors, Rev. Mod. Phys. 60 (1988) 663–699.
  * (66) B. J. McParland, Medical radiation dosimetry.  Springer, 2016.
  * (67) G. R. Lynch and O. I. Dahl, Approximations to multiple Coulomb scattering, Nucl. Instrum. Meth. B 58 (1991) 6–10.