Dark Scalars and Heavy Neutral Leptons at DarkQuest

The proposed DarkQuest beam dump experiment, a modest upgrade to the existing SeaQuest/SpinQuest experiment, has great potential for uncovering new physics within a dark sector. We explore both the near-term and long-term prospects for observing two distinct, highly-motivated hidden sector benchmark models: heavy neutral leptons and Higgs-mixed scalars. We comprehensively examine the particle production and detector acceptance at DarkQuest, including an updated treatment of meson production, and light scalar production through both bremsstrahlung and gluon-gluon fusion. In both benchmark models, DarkQuest will provide an opportunity to probe previously inaccessible interesting regions of parameter space on a fairly short timescale when compared to other proposed experiments.


Introduction
The hypothesis of a light, weakly coupled 'dark' or 'hidden' sector has received considerable attention in recent years. Though neutral under the Standard Model (SM) gauge group, dark sectors may exhibit rich dynamics, such as new forms of matter, new dark symmetries and forces, confinement, or spontaneous symmetry breaking, that could address some of the deficiencies of the SM. For example, the dark matter may be part of such a sector, communicating with the visible sector through a weakly coupled mediator, or the neutrino mass generation could be connected to new gauge singlet fermions within a dark sector.
A vibrant experimental program to search for light weakly coupled particles has emerged over the last decade and promises to be a fertile area of research for many years to come; for a recent summary of existing and planned efforts, see the community studies [1][2][3][4]. Among the critical components of this program, particularly in exploring GeV scale dark states, are proton beam fixed target experiments [5][6][7]. In these experiments, an intense proton beam impinges on a target, producing a torrent of SM particles alongside a smaller flux of relativistic dark sector particles. Due to their suppressed coupling to the SM, once produced these dark particles can travel macroscopic distances before decaying downstream into visible particles. Given a suitable detector apparatus, the visible decay products can then be identified, characterized, and discriminated from potential background sources, which provides a promising means to probe and discover new light weakly coupled states.
One particularly promising experiment is DarkQuest, a mild augmentation of the SeaQuest and SpinQuest experiments [8]. The proposed DarkQuest upgrade entails the addition of an electromagnetic calorimeter (ECAL) to the existing SeaQuest muon spectrometer, which will extend the physics capabilities of the experiment. These new capabilities will allow for Dark-Quest to produce a suite of sensitive searches for dark particles decaying to a wide variety of SM final states such as electrons, muons, charged hadrons, and photons [9][10][11][12][13][14][15]. The experiment's high luminosity coupled with its short baseline would allow for sensitivity to both fairly short-lived particles (cτ 1 m) and more weakly-coupled particles with fairly low production rates. Although a variety of other experimental proposals targeting dark sectors exist, DarkQuest is exceptional because most of the detector and infrastructure currently exists, is one of the few beam dump experiments with access to a high energy proton beam, would have an impressive range of sensitivity, and could provide novel results in comparatively short timescale.
In this work, we will study the potential sensitivity of DarkQuest to two highly motivated dark sector particles -dark scalars and heavy neutral leptons (HNLs). Dark scalars that mix through the Higgs portal provide one of the simplest extensions the SM and may be connected to a variety of puzzles such as dark matter [16], inflation [17], and naturalness [18]. Heavy neutral leptons (also called right-handed neutrinos or sterile neutrinos) are strongly motivated by the observation of neutrino masses [19][20][21][22][23][24] and GeV-scale HNLs may also play a role in the generation of the matter-antimatter asymmetry [25,26]. As we will demonstrate, DarkQuest has excellent prospects to explore substantial new regions of parameter space in these scenarios. Along with previous studies targeting a variety of dark sector models [9][10][11][12][13][14][15], our results lend further strong motivation for the DarkQuest ECAL upgrade, which will provide the basis for a rich and exciting experimental search program in the coming 5-10 years.
The paper is organized as follows. In Section 2 we provide an overview of DarkQuest along with a general discussion of the methodology used in our sensitivity estimates. In Section 3 we consider the prospects for HNLs searches at DarkQuest, while searches for dark scalars are covered in Section 4. We present a summary of our results in Section 5. In the Appendix, we provide the details about our calculation of dark scalar production.

The DarkQuest Experiment
The E906/E1039 SeaQuest/SpinQuest experiment is a proton fixed target beam dump spectrometer experiment on the neutrino-muon beam line of the Fermilab Accelerator Complex [8]. A schematic layout of the experiment is shown in Figure 1. A high-intensity beam of 120 GeV protons (center of mass energy √ s 15 GeV) is delivered to a thin nuclear target. The target is situated ∼ 1 m upstream of a 5 m long, closed-aperture, solid iron dipole focusing magnet ("FMAG"), which magnetically deflects soft SM radiation and also functions as a beam dump for the large majority of protons that do not interact in the target. This effectively allows only high energy muons, neutral kaons, and neutrinos to traverse the FMAG. The spectrometer consists of a high precision tracking system (St-1/2/3 tracking) and a muon identification system (absorber and St-4 muon ID). An additional 3 m long open-aperture magnet ("KMAG") is positioned at z = (9 − 12) m and delivers a transverse momentum impulse of ∆p KMAG T ∼ 0.4 GeV, enabling accurate momentum reconstruction of charged particles. In addition, in 2017 displaced vertex trigger hodoscopes were installed on both sides of the KMAG (see Figure 1), allowing for the detection of muons originating from the decays of exotic light long-lived particles after the dump. The experiment has been approved to collect ∼ 10 18 protons on target in the coming two years, until 2023.
On the horizon, there are plans to install a refurbished electromagnetic calorimeter (ECAL) from the PHENIX experiment [27] between St-3 and the absorber wall (see brown region in Figure 1). This will allow the upgraded experiment, DarkQuest, to search for a much broader set of dark sector displaced signatures, including electrons, charged pions and kaons, and photons. The DarkQuest experiment has a relatively compact geometry, making it well-suited to search for dark particles with O(10 cm − 1 m) lifetimes that are currently hidden to previous beam dump experiments with a much longer baseline.
Additional possible upgrades of the experiment ("LongQuest") have been also proposed [14]. This includes additional trackers and calorimeters after station 4 of the SeaQuest spectrometer. Figure 1. Layout of the DarkQuest experiment. The SeaQuest experiment has the same layout, except for the ECAL (dashed brown region located near z ∼ 19 m) [28].
The ultimate detectability of long lived dark particles at DarkQuest depends on several key factors. These include the production rate and kinematical properties of dark particles, their decay properties including branching ratios to final states containing charged particles and lifetime, the detector acceptance, and any potential SM background processes. In the remainder of this section we provide a brief discussion of these issues, which will motivate us to define two distinct run scenarios to be used later in our sensitivity projections for HNLs and dark scalars.

DarkQuest luminosity scenarios, Phase I and Phase II
At DarkQuest both HNLs and dark scalars can be produced in meson decays (e.g., K, D, and B mesons), while scalars can also be produced in the primary proton interactions through the proton bremsstrahlung and gluon fusion processes. Assuming every proton interacts in the dump, an estimate of the effective integrated luminosity at DarkQuest is given by 1 where N p is the total number of protons on target, λ int = 16.77 cm [29] is the nuclear interaction length in iron, ρ = 7.87 g cm −3 is the density of iron, and N A is the Avogadro's number. In the second equality, we assume the per nucleon cross-section is the total crosssection on iron times the mass number, A = 56. A related quantity often seen in the literature is the the total hadronic cross per nucleon, which in iron is given by σ pN ≡ (λ int ρ N A ) −1 12.6 mb. We will consider two benchmark luminosity scenarios in our projections below: a "Phase I" corresponding to N p = 10 18 (L ∼ 79 ab −1 of integrated luminosity) which can be achieved on the couple of years time scale, and a "Phase II" scenario corresponding to N p = 10 20 (L ∼ 7.9 zb −1 of integrated luminosity) which could potentially be collected over a longer time frame [30].

Meson production at DarkQuest
Given the considerable energy of the Main Injector protons and the substantial anticipated luminosity, mesons such as kaons, D-mesons, and B-mesons, as well as τ -leptons, are abundantly produced at DarkQuest. Much of hidden sector particle production at DarkQuest thus occurs through the decays of these SM states. Here we discuss our approach to modeling meson production at DarkQuest.
Kaons have an enormous production rate in primary proton collisions at DarkQuest, with an order one number of kaons produced per proton on target. However, since kaons are long lived and typically produced with boosts of order 10, their lab frame decay length is generally much longer than the characteristic hadronic interaction length, causing a significant Table 1. Number of mesons and τ s produced for N p = 10 18 . For kaons, we present the number that decay before one nuclear interaction length, λ int , where the asterisk merely serves to flag that these are not the total amount produced. For taus, we present both those produced directly from electroweak interactions (first entry) and those from D s decays (second entry). The values shown are the sum of the production of the two mesons (e.g., particle + anti-particle).
attenuation of the kaon flux as they traverse the dump. Taking this into account, the number of kaons that decay before the first interaction length can serve as a useful proxy for the opportunities to produce hidden sector particles, where λ K ≈ 20 cm is the kaon interaction length 2 , n K i ∼ 0.2 is the number of kaons produced per proton on target at DarkQuest for each of K + , K − , K 0 L , and K 0 S , and γ −1 K ∼ 0.1 is the mean inverse Lorentz boost. Both n K i and γ −1 K were estimated using PYTHIA 8 [31]. The values for N K ± , N K 0 L , and N K 0 S that decay before the first interaction length are shown in Table 1. As expected, the number of K 0 S is much larger than the number of K 0 L and K ± due to their much shorter lifetime.
For D-meson production, we follow an approach that is similar to the one used by the SHiP experiment at CERN [32]. We compute the pp → D 0 ,D 0 production cross section as a function of √ s, using PYTHIA 8 [31] with CTEQ6 LO parton distribution functions (PDFs) [33]. We rescale these cross sections to match the cross sections measured in the interval √ s = (20 − 40) GeV [34,35]. Using this rescaling, we estimate σ(D 0 ,D 0 ) ∼ 1 µb at √ s = 15 GeV. Using the fragmentation fractions for charm production, we obtain a charm production cross section σ cc = σ(D 0 ,D 0 )/f (c → D 0 ) ∼ 1.6 µb. To estimate the fragmentation fractions, we generate hard cc processes in PYTHIA 8 [31] at the DarkQuest energy and extract the ratios. As a cross check, we have also used PYTHIA 8 to estimate the B and D fragmentation fractions at SHiP and LHC energies, finding relatively good agreement with the values quoted in Ref. [36]. The number of charm mesons produced for N p = 10 18 is shown in Table 1 for D ± , D 0 andD 0 , and D ± s . We follow a similar procedure to compute the production rate of B-mesons. In Table 1, we report the number of mesons produced for N p = 10 18 . Due to 2m B + 2m p ∼ √ s, there is substantial uncertainty on σ bb at DarkQuest beam energies. In particular, Monte Carlo estimates with differing PDF choices can result in largely different values for the projected cross-section. This can be primarily understood from the high uncertainty at large momentum fraction. Unlike in the case of charm, we do not have empirical data to extrapolate from in a controlled manner. Through exploring a variety of PDF choices, we found roughly an order of magnitude discrepancy for the projected cross-sections σ(pp → bb) ∼ 0.5 − 5 pb. Given this range, we choose σ(pp → bb) = 1 pb throughout this work.
In addition to meson decays, τ ± decays can produce dark sector particles. At DarkQuest, the primary way of producing a τ lepton is through the decay of a D s meson with Br(D s → τ ± ν τ ) = (5.55 ± 0.24)% [37], which provides over an order of magnitude more τ s than the direct electroweak production (see Table 1, where the first entry represents the number of τ ± directly produced through electroweak processes).
We can compare the numbers in Table 1 to the numbers obtained for higher energy proton beams as, for example, the 400 GeV SPS proton beam. The number of kaons [38], D-mesons [36], and taus [39] produced per proton on target is suppressed only by roughly an order of magnitude at the Fermilab Main Injector. A much larger suppression applies to B-meson production [36], for which the Main Injector loses roughly three orders of magnitude. For this reason, we generally expect DarkQuest to achieve a similar reach for dark sector states produced from light meson or tau decays. Importantly, with the exception of D-mesons, most of these estimates consider only the particles produced in the incident protons primary interaction. Secondary interactions of hard particles and beam remnants within the beam dump can also produce additional kaons and taus, which could potentially enhance the flux of dark particles. The differential rates for these secondaries should be carefully evaluated in order for DarkQuest to most precisely state their sensitivity to a variety of models. In this sense, our estimate of the reach should be considered conservative.

Detector acceptance of DarkQuest
Next, we turn to the issue of the detector acceptance. Our considerations and approach to modeling the effect of the KMAG magnetic field and acceptance closely follows Ref. [11].
A Monte Carlo simulation is used to compute the total detection efficiency. In particular, we will consider signal events to be those in which the dark particle decays to final states containing two quasi-stable charged particles (i.e., electrons, muons, charged pions, and charged kaons) within a fiducial decay region at position z ∈ (z min , z max ), located downstream of the FMAG. The daughter charged particles are then required to intersect tracking station 3, assumed to be a 2 m × 2 m square centered about the beam line and located approximately 18.5 m downstream of the dump (see Figure 1). We also model the effect of the KMAG magnetic field on charged particles trajectories by an instantaneous transverse momentum impulse of ∆p T = 0.4 GeV × (∆z K /3m) applied in thex direction halfway through the particle's KMAG traverse, where ∆z K is the distance traveled by the daughter particles through the KMAG 3 .
The total detection efficiency is then estimated according to [11] where m, Γ, and p z are the mass, width, andẑ-component of the momentum of the dark particle, respectively. The sum in (2.3) is carried out over those events falling within the geometric acceptance as described above, and N MC represents the total number of simulated events. We will define two fiducial decay regions for our study that will be associated with our near future and long term run scenarios. As we will discuss in Secs. 3.3, 4.3, the detection efficiency for the two fiducial decay regions is relatively sizable, ranging from ∼ few ×10 −2 to ∼ 1, depending on the particular production and decay mode of the dark particle.
For our Phase I scenario, we require that the dark particle decays within the 5 m − 6 m region immediately downstream of the FMAG. The main advantages of this choice are that the charged daughter particles are tracked in Station I and their trajectories are bent by the KMAG magnetic field, making accurate momentum reconstruction feasible and greatly helping with particle identification, vertex reconstruction, and background rejection.
For our Phase II scenario we will consider the longer fiducial decay region of 7 m − 12 m. Given the higher luminosity in our Phase II scenario, we expect more background events, e.g., from K 0 L particles which pass through the FMAG and decay semileptonically. As discussed in Ref. [11], these backgrounds could be further mitigated with additional shielding in the 5 m − 7 m region, partially explaining the motivation of the 7 m − 12 m fiducial region. In addition, the 7 m − 12 m fiducial region would increase the geometric acceptance. While this choice allows for an appreciable enhancement of the overall signal rate and for additional suppression of backgrounds, it is not without additional challenges. For example, momentum reconstruction will be more challenging since the daughter particles would not pass through the first tracking station.
Our benchmark scenarios discussed here should be considered as preliminary, and a dedicated study of the potential backgrounds and signal region optimization is warranted. The DarkQuest collaboration is currently investigating the several sources of backgrounds, with a focus on the e + e − signature characteristic of dark photons [40]. While awaiting a definitive study from the collaboration, a crude estimate suggests that it will be possible to observe signals over the K 0 L decay backgrounds. For the signatures investigated in this paper, the dominant sources come from the production of K 0 will be produced in the beam dump during Phase I. Taking the kaon interaction length in iron to be ∼ 20 cm, we expect approximately ∼ 10 6 kaons to escape FMAG, and O(10 4 ) of which will decay in 5 m -6 m. Accounting for branching ratios and geometric acceptance, we find that, depending on the particular final state, O(100 − 1000) K 0 L will decay in the fiducial region with decay products detected by DarkQuest. Despite the substantial increase in luminosity, the situation during Phase II can be much improved over Phase I provided additional shielding is in place between 5 m -7 m. While approximately ∼ 10 19 K 0 L will be produced in Phase II, a similar estimate as given for Phase I suggests that depending on the specific final state, O(1 − 10) K 0 L will traverse 7 m of iron, decay in the 7 m -12 m fiducial region, and will lead to detectable decay products. Depending on the final state signature, additional handles can be utilized to further mitigate these backgrounds. In Secs. 3.4, 4.4, we will estimate how many of these K 0 L will result in background events for the several signatures. When discussing the DarkQuest reach for dark scalars and HNLs, we will require 10 signal events, but the true requirement against background may be more or less depending on the expected background population specific to the mass and decay paths.

Heavy Neutral Leptons
Heavy neutral leptons (HNL),N i , can interact with the SM neutrinos through the neutrino portal operator where H is the SM Higgs doublet andL i = (ν i , i ) T is the SM lepton doublet of flavor i. Because of these operators, after electroweak symmetry breaking, the HNLs will mix with the SM neutrinos. We will refer to the unhatted fields ν i and N i as the corresponding mass eigenstates of the light SM neutrinos and HNLs, respectively, and the relation between the flavor and mass bases is described by a mixing matrix, U . The phenomenology of HNLs largely follows from their induced couplings to electroweak bosons, which in the limit of small mixing angles are given by Additionally, we will assume that N is a Majorana particle throughout this work. Majorana HNLs are particularly motivated as they arise in the Type-I seesaw mechanism for neutrino mass generation. While the Type-I seesaw naively leads to mixing angles of parametric size ∼ m ν /m N , which is extremely small for GeV-scale HNLs, we note that there are schemes such as the inverse seesaw [41][42][43] and linear seesaw [44] where the mixing angles can be much larger. For the purposes of characterizing the DarkQuest sensitivity, we will take a phenomenological approach, as is commonly done in the literature, assuming the existence of a single HNL state, N , in the mass range of interest, which dominantly mixes with a particular neutrino flavor, i.e., dominant electron-, muon-, or tau-flavor mixing. In this case, the phenomenology is dictated by the HNL mass, m N , and mixing angle, denoted by U e , U µ , or U τ , respectively, for the three mixing scenarios. If these assumptions were relaxed, we expect the phenomenological implications relevant for DarkQuest are typically only slightly different than in a flavor-aligned case.

HNL production
As a consequence of the interactions in (3.2), HNLs can be copiously produced at DarkQuest through the decays of mesons and τ leptons. Meson and τ production at DarkQuest is   Table 1. For example, HNLs can be produced in the two body decays of charged pseudoscalar mesons, P → i N . In the regime m m N m P , the branching ratio is given by [45] Br where τ P , f P , and m P are the meson lifetime, decay constant, and mass, respectively, and the CKM matrix element, V αβ , is dictated by the valence quark content of P (e.g., V cd for D ± , etc.). The two body decay rates (3.3) scale as m 2 N as a consequence of the chirality flip, and are thus enhanced for heavier HNLs.
Three body decays of mesons to HNLs are also important and can even be the dominant production mechanism depending on the HNL mass. Although phase space suppressed, the three body meson decay rates do not suffer from the CKM or chirality flip suppressions characteristic of the two body decays in (3.3). HNLs can furthermore be produced through τ decays (e.g., two body decays involving hadronic resonances, or three body leptonic decays) and are subject to similar considerations.
For all meson and τ branching ratios, we use the expressions in Ref. [45]. The total number of HNLs produced at DarkQuest through different pathways is summarized in Figure 2, where we utilized a luminosity of 10 18 protons on target.

HNL decays
Once produced at DarkQuest, HNLs will decay through the weak interactions (3.2) to a variety of SM final states. Since their decays proceed through an off-shell heavy electroweak boson, GeV-scale HNLs are generically long lived and can easily traverse the beam dump at In each figure, we show the branching ratios into three SM neutrinos ννν (blue), e ± π ∓ or µ ± π ∓ (gold), νπ 0 (green), one neutrino and two charged leptons of any flavors (red), and one neutrino and two muons (dotted red). The thick black curve represents the sum of the branching ratios into two or more charged tracks.
DarkQuest before decaying. There is a rich variety of HNL decay modes, including a pseudoscalar meson and a lepton, a vector meson and a lepton, a lepton and two or more pions, or three leptons (including three neutrinos). We note that there is some disagreement in the literature about the corresponding rates. We have verified the results of Refs. [45,46], and utilize these expressions for the neutrino decays.
In Figure 3 we show the branching ratios of HNLs in the e−aligned, µ−aligned, and τ −aligned case (left, center, and right panel, respectively). For HNL masses below 1.5 GeV, we determine the total hadronic rate as the sum of exclusive meson decay rates, while above 1.5 GeV, we switch to using the inclusive N → qq rate, assuming exclusive rates are contained within this value. As we can observe from the figure, the branching ratio into the invisible ννν final state (in blue in the figure) is quite subdominant as long as the HNL has a mass above the pion mass. The other channels presented in the figure contain visible particles that are in principle observable by DarkQuest. The red dotted curve represents the decay into one neutrino and two muons. The corresponding branching ratio is also relatively suppressed, especially in the e−aligned, and τ −aligned scenarios. This is the only channel that can be easily identified now by the SeaQuest experiment, without the ECAL upgrade.
Provided the ECAL upgrade is installed, DarkQuest will have the capability to also search for a variety of HNL decays containing multiple charged particles in addition to muons. Among all visible channels, the π 0 ν channel is likely to be the most difficult one because of the challenging π 0 identification and large sources of backgrounds arising e.g., from the SM K 0 L → 3π 0 , K 0 S → π 0 π 0 processes, where some of the pions are missed or misidentified by the detector. For this reason, in the calculation of the DarkQuest reach on HNLs, we conservatively do not include this channel. The bold black line in Figure 3 shows the observable branching ratio used in this work, which is obtained by summing all branching ratios resulting in at least two charged particles.
In estimating the sensitivity below we will require 10 signal events, working under the assumption that backgrounds can be brought down to the level of a few events. The FMAG, i.e., the 5 m magnetized beam dump, serves to mitigate most of the backgrounds by sweeping away charged particles and largely blocking the most dangerous neutrals. Several potential sources remain and the ultimate size of these is the subject of current study [40]. One of the most relevant backgrounds comes from K 0 L particles that penetrate the dump and decay in the fiducial region. As we discussed in Sec. 2.3, we expect O(100 − 1000) of such K 0 L in Phase I and O(1 − 10) in Phase II. The decay K 0 L → π ± e ∓ ν will be background to the N → e + e − ν and N → e ± π ∓ signatures presented in Fig. 3. For the former, a pion rejection factor of order ∼ 1% will be sufficient to suppress the K 0 L → π ± e ∓ ν background to O(10) (< 1) events for Phase I (Phase II). This level of electron-pion discrimination should be feasible with the planned ECAL upgrade [27]. For the latter signal, the background could be suppressed through suitable kinematic cuts such as a cut on the m eπ invariant mass. However, a detailed study of these possibilities requires a careful modeling of K 0 L production in the FMAG, which is beyond our current scope. For signatures involving muons, the existing SeaQuest spectrometer already has the capability to distinguish muons, which pass through the absorber and are detected in the Muon-ID system (see Figure 1), from charged hadrons, which do not penetrate the absorber. As above, muonic backgrounds to the N → µ ± π ∓ signature can arise from decays such K 0 L → π ± µ ∓ ν, while the N → µ + µ − ν channel should have very small backgrounds.

Detector acceptance
We follow the procedure outlined in Sec. 2.3 to compute the geometric acceptance for HNLs at DarkQuest. To reduce the complexity for a clear presentation, we show in Figure 4 the normalized geometric efficiency in the large lifetime limit. To compute these curves, we consider the µ-aligned scenario and the large lifetime regime, i.e., we assume that the HNL decay length is much larger than the detector size so that the differential probability to decay is a constant with distance, and normalize to only the particles that decay within the fiducial region. This limit is relevant for small mixing angles. The different colored curves in Figure 4 correspond to several representative production and decay modes of the HNL. The lighter (darker) curves represent the acceptance for Phase I (5 m -6 m) (Phase II (7 m -12 m)). Overall, the acceptance is relatively large, ranging from a few % to ∼ 20% depending on the HNL production/decay mode, and is fairly constant with the HNL mass. As expected, the acceptance for Phase I is somewhat smaller than that for Phase II, since for Phase II the HNLs decay typically closer to tracking station 3.

The DarkQuest reach for HNLs
With our estimates for HNL production, decays, and experimental acceptance in hand, we can compute the total number of signal events in the SM final state i expected at DarkQuest Figure 4. Geometric acceptance as a function of the HNL mass normalized to the number of HNLs decaying within the fiducial decay region in the large lifetime limit (i.e., the HNL decay length is much larger than the detector size). We show separately the efficiency for HNLs that are produced and decay through several representative channels, including K → N, N → eeν (green), D → N , N → eeν (blue), D → K N, N → µπ (orange) and B → D N, N → µµν (red), and for two run scenarios: Phase II, 5m -6m (lighter darker), Phase II, 7m -12m, (darker color).
Here N N is the number of HNLs produced in a given production channel (see Section 3.1 and Figure 2), Br i is the branching ratio for N → i (see Section 3.2, and Figure 3), and eff i is the experimental efficiency to detect the final state i, computed using (2.3). A summary of the projected reach is shown in Figure 5 for µ-and τ -flavored HNLs decaying inclusively to final states containing two or more detected charged tracks. The solid black (dashed black) contour specifies the HNL mass -squared mixing angle parameters leading to 10 signal events according to (3.4) for the Phase I (Phase II) run scenario. We note that the projected reach for e-aligned HNLs is very similar to the µ-aligned reach shown in Figure 5. For this reason, we do not show the e-aligned scenario in the figure. We also show in the shaded gray regions the existing experimental or observational limits, including CHARM [47], PS191 [48], DELPHI [49], NuTeV [50], E949 [51], MicroBooNE [52], T2K [53], ATLAS [54], Belle [55], and Big Bang Nucleosynthesis (BBN) [56] 4 (see e.g., Ref. [3] for a thorough discussion of these limits). For comparison, we also display the projected sensitivities to HNLs from several proposed experiments, including NA62++ [57], FASER [58], CODEX-b [59], MATHUSLA [60] and SHiP [61]. For additional proposals to probe GeV-scale HNLs see e.g., Refs.  We conclude that DarkQuest Phase I can probe a significant region of currently unexplored parameter space for τ -aligned HNLs. For the Phase II scenario, DarkQuest will be able to extend the sensitivity by more than one order of magnitude in the squared mixing angle compared to Phase I, while also covering new regions of parameter space in the µ-aligned scenario which are presently unconstrained.

Dark Scalars
We now consider dark scalars interacting through the Higgs portal. A new singlet scalar can couple to the SM Higgs through two renormalizable portal couplings, − L ⊃ (AŜ + λŜ 2 )Ĥ †Ĥ . (4.1) The dark scalar may acquire a small coupling to SM fermions and gauge bosons through its mass mixing with the Higgs, which will occur if the A = 0 in (4.1) or if the dark scalar obtains a non-zero vacuum expectation value. Then, in the physical basis, the phenomenology at DarkQuest is governed by the dark scalar mass, m S , and the scalar-Higgs mixing angle, θ: Given the experimental constraints on the mixing angle for dark scalars at the GeV-scale, we will always be working in the regime θ 1. We will not study the phenomenological consequences of additional couplings between the scalar and the Higgs, such as the cubic interaction hSS. While such a coupling can lead to additional scalar production processes such as B → KSS, these are typically not as important at DarkQuest as processes involving singly produced scalars. Such coupling also leads to Higgs exotic decays of the type h → SS Figure 6. Number of scalars produced at DarkQuest for K → πS (green), B → KS (blue), proton bremsstrahlung (red), and gluon fusion (black), assuming 10 18 protons on target and a mixing angle equal to 1.
[67] that can be searched for at the LHC. We do not include the corresponding bounds in our summary plot in Fig. 9, since these bounds depend on the hSS coupling that is independent from the mixing angle θ. We now discuss in more detail the production of scalars, their decays, the experimental acceptance, and the DarkQuest reach.

Scalar production at DarkQuest
At DarkQuest scalars are produced through three main processes: meson decays, proton bremsstrahlung, and gluon-gluon fusion. The sensitivity of DarkQuest to scalars produced through B meson decays was already studied in Ref. [11]. In this work we will also examine the potential additional sensitivity from scalars produced through kaon decays, proton bremsstrahlung, and gluon-gluon fusion. Figure 6 shows the number of dark scalars produced through these three production channels as a function of the scalar mass, assuming 10 18 protons on target. Low mass scalars are dominantly produced in kaon decays. Above the m K − m π threshold and in the vicinity of m S ∼ 1 GeV, proton bremsstrahlung dominates, while heavier scalars can be produced through B-meson decays and gluon fusion.

Meson decays
We first consider scalar production through meson decays. We refer the reader to Sec. 2.1 and Table 1 for a summary of meson production at the DarkQuest. We first consider scalars produced through kaon decays, K → πS, which is especially relevant for lighter scalars. The partial decay width for K ± → π ± S is [17,[68][69][70][71] 5 Using these partial widths and (2.2), the number of scalars produced from kaon decays in a thick target can be estimated as [71] where n K ∼ 0.6 is the number of K ± and K 0 L produced per proton on target. Next, we consider scalars produced through B meson decays, which proceeds through b − s − S penguin transitions. The inclusive branching ratio for B → X s S can be written as [72][73][74] (see also [71,75] for exclusive B decays) where Φ ≈ 0.5 is a phase space factor. Using the measured inclusive rate for B → X c e ν [37], Since B-mesons decay promptly, we can estimate the number of scalars produced in their decays as N S = N p n B Br(B → X s S) ∼ 10 9 × θ 2 N p 10 18 (B meson decays), (4.6) where n B ∼ 10 −10 is the number of B mesons produced per proton on target at DarkQuest.

Proton bremsstrahlung
Next, we turn to scalars produced through proton bremsstrahlung, p + p → S + X. The cross section is obtained following the calculation in Ref. [76], which is based on the generalized Weizsacker-Williams method [77]; further details are provided in Appendix A.1. Specifically, scalar events are generated by sampling the differential cross section dσ brem /dz dp 2 T , where z ≡ p S /p p is the fraction of the proton beam momentum, p p , carried by the emitted scalar, with p S the scalar momentum, and p T is the scalar transverse momentum. The validity of the Weizsacker-Williams approach relies on the kinematic conditions p p , p S , p p − p S m p , |p T |. To satisfy these conditions for DarkQuest that uses 120 GeV protons, we follow Ref. [11] and restrict the phase space to the range z ∈ (0.1, 0.9) and p T < 1 GeV. We note that these conditions are slightly more restrictive than those used in Ref. [76], leading to an integrated cross section that is smaller by an order one factor. 5 Although the branching fractions are different, the partial widths are very similar, and the total width cancels out of the estimate (2.2) as long as λK γ −1 K cτi. In fact, K 0 S has λK ∼ γ −1 K cτK S suggesting its total width could also cancel out of the expression (up to an O(1) factor). However, K 0 S is not included in our analysis since the partial width Γ(K 0 S → π 0 S) Γ(K 0 L → π 0 S), so it can be neglected for that reason.
The total bremsstrahlung cross section is estimated to be where σ pp ≈ 40 mb is the total inelastic proton-proton cross section and the factor in parentheses gives the approximate integrated probability of scalar emission. The parameter g SN N is the zero momentum scalar nucleon coupling (for θ = 1) and F S (p 2 S ) is a time-like scalarnucleon form factor, which is discussed in detail in Appendix A.1. Including order one factors arising from phase space integration, we estimate the total number of scalars produced in proton bremsstrahlung to be N S ∼ 10 11 θ 2 N p 10 18 (Proton Bremsstrahlung). (4.8) Figure 6 shows the total number of scalars produced at DarkQuest as a function of the scalar mass. The large resonant enhancement near m S ∼ 1 GeV is a consequence of mixing with the narrow f 0 (980) scalar resonance, while the bremsstrahlung cross section drops steeply for m S 1 GeV due to the form factor suppression. It is likely that the zoo of heavy f 0 resonances would belay this high mass suppression, but we make no attempt to model that here. The uncertainty band is obtained by varying the lower integration limit for z between 0.05 and 0.2 as well as the scalar resonance masses and widths in the form factor F S (p 2 S ). We note that the rates for scalar production from bremsstrahlung have a rather mild dependence on the proton beam energy, and thus the production rate at higher energy facilities such as the CERN SPS (400 GeV protons) is very similar to that at DarkQuest.

Gluon fusion
The final process we consider is scalar production via gluon fusion. As in the case of the SM Higgs boson, this process proceeds at one loop through the heavy quark triangle diagrams. The full leading order cross section is discussed in Appendix A.2. We restrict our analysis to scalar masses above O(1 GeV) where the perturbative QCD computation is valid. In this mass range, the cross section is of order σ ggS ∼ 30 pb × θ 2 (m S /1GeV) −2 , and the number of scalars produced is therefore As in the case of the SM Higgs boson, we expect higher order corrections to enhance the rate by an order one factor, although we are not aware of an existing calculation in the literature that can be applied to such light scalars. While it would be interesting to study this question further, we will simply apply a K-factor equal to 1.5 in our estimate of the rate, which is similar to that of the SM Higgs boson. For our simulation, we use the HEFT model in MadGraph5 amc@nlo [78] to generate scalar events, which are then passed to PYTHIA 8 [31] for showering. While we find that gluon fusion is generally subdominant to ��� θ other production mechanisms (see the black curve in Figure 6), it can give some additional sensitivity in the 1-2 GeV scalar mass range, particularly for the Phase II scenario. For comparison, we find that the scalar production via gluon fusion is only about a factor of 2 larger at the higher energy CERN SPS.

Scalar decays
Through its mixing with the Higgs, the scalar will decay to SM final states. For example, the dark scalar can decay to charged leptons with a partial decay width, Γ S→ + − θ 2 m 2 m S /(8πv 2 ). Above the two pion threshold the scalar can also decay to hadronic final states. The theoretical description of such decays is complicated by strong interaction effects, leading to significant uncertainties in the predictions for masses of order 1 GeV. In our study we will use the results and prescriptions from the recent study in Ref. [71]. In particular, for relatively low scalar masses in the few hundred MeV range, the hadronic decays are well described using Chiral Perturbation Theory [79,80]. At higher masses, m S 2 GeV, the perturbative spectator model can be used to compute the decay rates to quarks and gluons [81]. In the intermediate regime of m S ∼ 1 − 2 GeV an analysis based on dispersion relations can be employed to estimate the partial decay widths for scalar decays to pairs of pions and kaons [71,80,[82][83][84]. Furthermore, Ref. [71] includes an additional contribution to the scalar decay width to account for other hadronic channels above the 4π threshold. Despite the formidable calculations involved in estimating the decays in these regimes, these are uncontrolled approximations and should be viewed with healthy skepticism [85]. The scalar branching ratios in the e + e − , µ + µ − , π + π − , and K + K − channels, as well as the scalar decay Figure 8. Geometric acceptance as a function of scalar mass normalized to the number of scalars decaying within the fiducial decay region in the infinite lifetime limit. We show separately the efficiency for scalars produced via proton bremsstrahlung (red), B decays (blue), and kaon decays (green), and for three run scenarios: Phase I, 5 m − 6 m (light shading), Phase II, 7 m − 12 m, (medium shading) and Phase II, 7 m − 12 m, without the KMAG (dark shading). The acceptance combines the e + e − , µ + µ − , π + π − , and K + K − final states weighted by their relative decay rates. length, are shown in Figure 7.
As with our HNL projections presented in Sec. 3.4, we will require 10 signal events in our dark scalar sensitivity estimates. The considerations leading to this assumption are similar to those outlined in Secs. 2.3 and 3.2. In particular, for the signatures arising from scalar decays to leptons, S → + − , there can be backgrounds from K 0 L that pass through the FMAG and decay via K 0 L → π ± ± ν, though we expect that detector level pion-lepton discrimination can be used to bring these backgrounds at the level of O(10) (< 1) events for Phase I (Phase II). For the hadronic scalar signatures such as S → π + π − , K + K − , there are backgrounds from the decays K 0 L → π − π + π 0 and K 0 L → π + π − . The corresponding background rates, particularly for the two pion decay, are further suppressed by the small branching ratios (BR(K 0 L → π + π − ) ∼ 2 × 10 −3 ), and we expect that kinematic information will be helpful in distinguishing the signal, though this remains to be studied in detail.

Detector acceptance
We follow the procedure discussed in Sec. 2.3 to account for the geometric acceptance of the experiment, with the total detector efficiency computed according to Eq. (2.3).
In Figure 8 we display the geometric acceptance as a function of scalar mass in the infinite lifetime limit, normalized to the number of scalars decaying within the fiducial decay region. This limit is of practical importance for much of the small θ parameter space. Several notable features can be observed in Figure 8. First, the overall efficiency is higher for dark scalars produced in proton bremsstrahlung compared to those from B and kaon decays. This is due to the larger typical Lorentz boosts of scalars originating in the former process, which inherit an order one fraction of the beam energy. Second, an increase in the efficiency is typically observed as m S increases beyond the dimuon threshold. Due to phase space suppression, heavier particles produced through scalar decays will typically be more collinear with the parent scalar, which leads to a higher overall acceptance. Furthermore, in the decays to electrons, the emitted particles are highly relativistic in the scalar rest frame and the fraction emitted towards the negative z direction can have a small lab frame longitudinal momentum. Such electrons can be swept out of the detector as they pass through the KMAG, explaining in the lower observed efficiency when the KMAG is present. Furthermore, we see that for heavy scalars produced via bremsstrahlung and B-meson decays, the efficiency tends to decrease as the the scalar mass increases beyond O(1 GeV) since in this regime the daughter particle p T inherited from the scalar mother increases approximately in proportion to m S and is generally larger than that imparted by the KMAG. Another trend observed in all production channels is the increased efficiency in Phase II (medium shading) over that in Phase I (lighter shading), which stems from the fact that for the Phase II scenario the scalars decay closer to tracking station 3.
Finally, we have displayed the efficiency for an alternate Phase II scenario in which the KMAG is removed and the charged daughters are not deflected. In this case, the daughter particles have a smaller characteristic transverse momentum, leading to a higher geometric acceptance as seen in Figure 8. However, it should also be emphasized that in this run scenario particle momenta measurement capability is likely to be significantly degraded. In fact, the magnetic field strength of the KMAG is tunable [86] and could impart a smaller p T kick than the 0.4 GeV used in this work. It would be interesting to study in detail its impact on the geometric acceptance and reconstruction capabilities.

DarkQuest sensitivity to dark scalars
Given the scalar production rates, decay branching ratios, lifetime, and experimental efficiency, we can now estimate the total number of signal events in the SM final state i according to the formula where N S is the number of scalars produced in a given production channel (see Eqs. (4.4, 4.6, 4.8, 4.9) for the number of scalars produced via K decay, B decay, bremsstrahlung, and gluon fusion, respectively). In Figure 9 we show the projected per-production-channel sensitivity of DarkQuest Phase I for scalars decaying inclusively to pairs of charged particles, specifically e + e − , µ + µ − , π + π − , and K + K − . Each contour indicates the scalar mass -mixing angle parameters predicting 10 signal events according to (4.10). We show three contours corresponding to distinct scalar production mechanisms, including kaon decays, B-meson decays, and proton bremsstrahlung. No sensitivity is obtained from the gluon fusion process alone in the Phase I run scenario. The gray shaded regions indicate parameter points that are ��� θ Figure 9. DarkQuest Phase I sensitivity to dark scalars corresponding to N p = 10 18 and 5 m -6 m decay region. The contours correspond to 10 signal events as obtained by adding the e + e − , µ + µ − , π + π − , K + K − channels, for dark scalars produced via K → πS (green), B → KS (blue), and proton bremsstrahlung (red). The gray shaded regions correspond to existing limits from past experiments; see text for further details.
excluded by past experiments, which will be discussed in more details below. We observe from Figure 9 that DarkQuest Phase I (5m -6m, N p = 10 18 ) will be able to explore a significant new region of parameter space, in particular for scalars produced through kaon decays and proton bremsstrahlung. Next, in Figure 10 we show the full DarkQuest sensitivity to scalars decaying inclusively to pairs of charged particles, now combining all S production channels, for both Phase I (solid, black) and Phase II (dashed, black) scenarios. In comparison to Ref. [11], which studied scalars produced only in B-decays, we find that the additional scalar production from kaon decays and proton bremsstrahlung can significantly expand the parameter space that can be probed by DarkQuest. 6 In the figure, we also show the current experimental bounds on dark scalar parameter space, including those from CHARM [71,87], LSND [88], E787/E949 [89,90], LHCb [91,92], and NA62 [93]. In addition, we also display sensitivity projections from several ongoing or proposed future experiments, including NA62 [3,94], SBND and ICARUS [95], Belle II [96] (see also Ref. [97]), FASER [98], CODEX-b [73], MATHUSLA [60] and SHiP [61]. See also e.g., Refs. [3,63,99,100] for further proposals to probe Higgs portal scalars in this mass range. 7 We observe that DarkQuest Phase I has the potential to cover a significant region of unexplored parameter space for scalar masses between about 200 MeV and 2 GeV. Phase II will probe angles as small as θ 5 × 10 −6 and as large as θ 10 −3 . 6 We have compared our projections with Ref. [11] for scalars produced via B decays and find good agreement. 7 We also note that a recent excess observed by the KOTO experiment can be explained in this scenario for scalar masses mS ∼ 150 MeV and mixing angles θ ∼ few × 10 −4 [101].

����� ��
�� �� �� �� ��� θ Figure 10. DarkQuest sensitivity to dark scalars. The contours correspond to 10 signal events as obtained by adding the e + e − , µ + µ − , π + π − , K + K − channels, for combined dark scalar production via K → πS, B → KS, proton bremsstrahlung and gluon fusion. We display both the DarkQuest Phase I sensitivity (solid, black) corresponding to N p = 10 18 and 5 m -6 m decay region, as well as the DarkQuest Phase II sensitivity (dashed, black) corresponding to N p = 10 20 and 7 m -12 m decay region. The gray shaded regions correspond to existing limits from past experiments. Also displayed are estimates from a variety of proposed experiments; see the text for further details and discussion.

Summary
We have investigated the sensitivity of the Fermilab DarkQuest experiment to two simple and well-motivated dark sector scenarios, heavy neutral leptons and Higgs-mixed scalars. The proposed DarkQuest ECAL upgrade will allow for sensitive searches to a variety of displaced final states containing charged particles and photons, which arise in the models considered here from the decay of long lived HNLs or scalars. We have carefully estimated the production and decay rates of these dark sector particles as well as the detector acceptance to derive projections under two benchmark run scenarios. During the Phase I scenario based on 10 18 protons on target and a 5m -6m fiducial decay region, DarkQuest will be able to explore significant new parameter space for τ -mixed HNLs and dark scalars in the mass range of a few hundred MeV -2 GeV. It is conceivable that this could be achieved on the 5 year time scale, putting DarkQuest on a competitive footing with other proposed experiments. Looking down the road, a potential Phase II scenario with 10 20 protons on target and a 7m-12m fiducial decay region would allow for improvements by more than one order of magnitude in terms of the interaction rates with SM particles (proportional to squared mixing angle). Our results build on past phenomenological studies [9][10][11][12][13][14][15] and provide further motivation for the DarkQuest ECAL upgrade. This upgrade can be realized with a relatively modest investment and will leverage the existing experimental infrastructure already in place to build an exciting dark sector physics program at Fermilab. We are not aware of any studies of the time-like scalar-nucleon form factor F S (p 2 S ) in the literature. In analogy with vector meson dominance model of the time-like electromagnetic form factor discussed in Ref. [103] (commonly used for dark photon production via proton bremsstrahlung [104,105]), we will assume that F S (p 2 S ) incorporates mixing with J P C = 0 ++ scalar resonances through a sum of Breit-Wigner components, where we include the three low-lying scalar resonances, φ ∈ {f 0 (500), f 0 (980), f 0 (1370)}. The decay constants f φ for each resonance are obtained by imposing the conditions F S (0) = 1 and F S (p 2 S ) ∼ 1/p 4 S as p 2 S → ∞ [106]. A central value is defined by taking the mean values of the masses, m φ = {475, 980, 1350} MeV, and widths, Γ φ = {550, 55, 350} MeV, leading to the decay constants f φ = {280, 1800, −990} MeV. To provide a naive estimate of the uncertainty, we vary the masses and widths of the resonances within their quoted uncertainty ranges [37], and take the envelope of the maximum and minimum values of |F S (p 2 S )|. The magnitude of the form factor is plotted in Figure 11.

A.2 Gluon Fusion
For scalars heavier than O(1 GeV), one can consider perturbative QCD production processes. In analogy with the SM Higgs boson, the dominant production channel is gluon fusion, gg → S. The production cross section can be written as σ ggS = θ 2 α 2 S (µ 2 R ) 1024 π v 2 q A 1/2 (τ q ) where τ q = m 2 S /4m 2 q , A 1/2 is a loop function (see e.g., [107]), A 1/2 (τ ) = 2 [τ + (τ − 1)f (τ )] τ −2 , (A. 6) with f (τ ) defined as Furthermore, L gg is the gluon parton luminosity function is the gluon PDF, and µ R (µ F ) is the renormalization (factorization) scale. To estimate the scale uncertainty in the cross section we fix µ F = µ R = µ and vary the scale between µ ∈ [ 2 3 m S , 4 3 m S ]. Our projections in the gluon fusion channel are made with the CT18NLO PDF set [108] and use the ManeParse package [109] for reading the PDF sets.