Sensitivity of the SHiP experiment to light dark matter

Dark matter is a well-established theoretical addition to the Standard Model supported by many observations in modern astrophysics and cosmology. In this context, the existence of weakly interacting massive particles represents an appealing solution to the observed thermal relic in the Universe. Indeed, a large experimental campaign is ongoing for the detection of such particles in the sub-GeV mass range. Adopting the benchmark scenario for light dark matter particles produced in the decay of a dark photon, with αD = 0.1 and mA′ = 3mχ, we study the potential of the SHiP experiment to detect such elusive particles through its Scattering and Neutrino detector (SND). In its 5-years run, corresponding to 2 · 1020 protons on target from the CERN SPS, we find that SHiP will improve the current limits in the mass range for the dark matter from about 1 MeV to 300 MeV. In particular, we show that SHiP will probe the thermal target for Majorana candidates in most of this mass window and even reach the Pseudo-Dirac thermal relic.


Introduction
One of the main challenges in particle physics today is figuring out the microscopic identity and the cosmological origin of dark matter (DM). The theoretical landscape is broad and it spans over many orders of magnitude in the mass/coupling parameter space. A compelling idea to explore is DM as a thermal relic of the early universe. The canonical example of this scenario is the Weakly Interacting Massive Particle (WIMP), a particle in the GeV-TeV mass range interacting with the visible sector via weak-sized interactions. Searches for WIMPs are in full swing [1,2]. However, the interesting parameter space goes beyond what has been explored in the past decade: thermal DM can be as heavy as 100 TeV or as light as a few keV. Recently, a lot of attention has been directed towards light DM (LDM) in the keV-GeV mass range [3].
Direct detection has traditionally employed the Migdal Effect [4] using liquid Argon [5,6] or liquid Xenon [7][8][9][10], while a novel strategy based on silicon devices has allowed to design new experiments optimised for sub-GeV DM, as SENSEI [11]. Since current DM direct detection experiments searching for elastic nuclear recoils rapidly lose sensitivity once the candidate mass drops below a few GeV [1,12], experiments at the intensity frontier represent an alternative yet appealing route and play an important role in this quest [3]. Fixed target experiments represent the prototype for such searches, although other collider experiments might be relevant in the same parameter space, as showed by the mono-photon searches at BaBar [13] and Belle II [14].
Here we present the sensitivity of the SHiP scattering and neutrino detector (SND), to LDM. We focus on the hypothesis that the DM couples to the SM through the exchange of a massive vector mediator, dubbed in the literature dark photon, and we have considered the cleanest signature given by the LDM-electron scattering. The scattering with nuclei, both coherent and deep inelastic scattering, although plagued by a larger neutrino background, might be an alternative detection strategy and will be the subject of a forthcoming dedicated analysis.
In a proton beam dump experiment signal yields are largely reduced as the interaction with the dark photon A is probed twice, if compared to electron fixed target experiments which make use of search strategies based on missing energy, such as NA64 [27], or missing momentum, such as the LDMX proposal [28]. Indeed, the LDM detection is achieved through its scattering within the downstream detector. Hence, the expected LDM yield scales as 4 α D ( being the interaction strength of the dark photon to SM particles and α D the LDM-A coupling), where a factor 2 comes from production and another 2 α D is due to detection. This has to be compared to the 2 scaling of typical missing energy/momentum experiments, which prove however to be not sensitive to LDM coupling constant α D . Due to their higher penetrating power and enhancements from meson decay reactions and/or strong interactions, proton beams are characterised by dark photon production rates larger than the ones achievable in electron beams of comparable intensity, which in part compensate for the detection suppression factor.
The potential to directly probe the dark sector mediator coupling α D , together with a wider sensitivity which encompasses other viable dark matter models, shows to a large extent the complementarity between the two approaches. This is even a more pressing aspect in the light of a possible discovery, as in general the observation of an excess alone is not sufficient to claim a discovery of a Dark matter particle. Indeed, intensity frontier probes do not depend on whether the particle χ produced through prompt DP decay is DM or not, as the only necessary ingredient is its stability concerning the target-detector distance. The observed excess might have an instrumental origin rather than a genuine New Physics (NP) effect. This applies also to the constraints that the SHiP experiment can place. With this regard, invaluable contribution could come from complementary DM observations from a cosmic source to unequivocally probe its thermal origin. In addition, since the SHiP experiment has a direct sensitivity to LDM interactions, we anticipate the possibility to use the time of flight measurement to uncontroversially discriminate massive NP particles from the SM neutrino background.
The paper is organised as follows: in section 2 we give a brief presentation of the model focusing on the main motivations. After introducing the SHiP experiment in section 3, we discuss the relevant production and detection mechanisms, in section 4. The detailed analysis of the neutrino background is the topic of section 5. We then pass to the discussion of the signal reviewing the main processes taken into account in our simulation. Finally, we show the main results on the sensitivity reach of the SHiP experiment in section 6 and we give our conclusions in section 7.

JHEP04(2021)199 2 Vector portal
Thermal freeze-out can naturally explain the origin of the DM relic density for a sub-GeV particle provided the interaction with the visible sector is mediated by a new light force carrier [2,29]. Here, we will consider as benchmark model the dark photon (DP) [30] vector portal where the DP A µ is the gauge boson of a new dark gauge group U(1) D kinetically mixed with the photon, and a scalar χ charged under U(1) D that serves as a DM candidate. Then, the low-energy effective Lagrangian describing the DM reads where: where is the DP-photon kinetic mixing parameter and m A is the mass of the DP, while: coupling and m χ is the mass of the dark matter particle. The region of the parameter space relevant for χ searches at beamdump facilities corresponds to m A > 2m χ and g D e which implies BR(A → χχ † ) ∼ 1. In case χ is DM, precise measurements of the temperature anisotropies of the cosmic microwave background (CMB) radiation significantly constrain the parameter space. In particular, they rule out Dirac fermions with mass below 10 GeV as a thermal DM candidate and more in general every DM candidate that acquires its relic abundance via s-wave annihilation into SM particles [31,32]. Hence, a complex scalar dark matter candidate χ is safe from such constraints as well as a Majorana or Pseudo-Dirac fermion. Tighter bounds come instead from the Planck measurement of the effective number of neutrino species N eff [32] and rule out the minimal DP model considered here if the complex scalar is lighter than 9 MeV [33].
In order to show the region of parameter space relevant for thermal freeze-out, we will present the SHiP sensitivity in the (m χ , Y ) plane where Y is defined as: . (2.4) In the assumption m A > 2m χ , the parameter Y is linked to the DM annihilation cross section via the formula [34]: where v is the relative velocity between the colliding DM particles.

The SHiP experiment
The Search for Hidden Particles (SHiP) experiment has been proposed as a general-purpose experiment [35] at the CERN Super-Proton-Synchrotron (SPS), addressed to explore the   high-intensity frontier for NP searches, thus complementing the LHC research program [35]. It is particularly targeted at the observation of long-lived weakly interacting particles of mass below 10 GeV/c 2 , foreseen in many Standard Model (SM) extensions. The use of a beam-dump facility [36] will result in a copious flux of charmed hadrons, from which not only hidden sector particles originate [37], but also tau neutrinos and their corresponding anti-particles. Therefore, being also a neutrino factory, SHiP will perform a wide neutrino physics program, as well as a first direct observation of the tau anti-neutrino, which represents the last particle to be directly observed to complete the SM framework. The SHiP Scattering Neutrino Detector (SND) is an apparatus designed for LDM particles searches, since it exploits an optimised combination of a dense target and high-granularity scattering detector, being it based on nuclear emulsion technology. In figure 1 a sketch of the experimental facility as currently implemented in the official simulation framework of the experiment FairShip [38] is shown. A brief overview of the simulated processes within FairShip and corresponding simulation software is reported in table 1.

JHEP04(2021)199
A 400 GeV/c proton beam will be delivered onto a thick heavy-metal hybrid target, specifically designed to maximise the charm production yield and thus hidden sector particles and tau neutrino yields. Over five years of SPS operations, a total of 2×10 20 protons on target (p.o.t.) collisions will be delivered, where each proton spill will have nominally 4×10 13 p.o.t.. A hadron stopper follows the beam-dump target, with the goal to absorb the SM particles produced in the beam interaction. In addition, a series of sweeping magnets, referred to as Muon Shield [39], act as a deflecting device for the residual muons, further cleaning the flux of particles from leftover backgrounds to hidden sector particles and neutrino searches.
The SND, shown in figure 2 in the setup adopted for this study, is located downstream of the muon sweeper. Placed in a magnetised region of 1.2 Tesla in the horizontal direction and perpendicular to the beam axis [40], it consists of a (90×75×321) cm 3 high-granularity tracking device which exploits the Emulsion Cloud Chamber (ECC) technique developed by the OPERA experiment [41], which was successfully used for tau neutrino detection [42,43]. Each elementary unit of the modular detector, called brick (figure 3), consists of alternating 56 lead plates of 1 mm thickness, passive material to increase the interaction probability, and 57 nuclear emulsion films of 0.3 mm thickness, acting as tracking detector with micrometric accuracy. It is worth noting that the brick also functions as a high-granularity sampling calorimeter with more than five active layers for every radiation length X 0 over a total thickness of 10 X 0 [44]. The ECC technology is also particularly efficient in the e/π 0 separation. The Compact Emulsion Spectrometer (CES), made up of a sequence of emulsion films and air gaps, is attached immediately downstream of the brick for electric charge measurement for particles not reaching the spectrometer. Despite the magnetic field, electron charge measurement is not possible due to early showering happening within the brick and the consequent information loss. The resulting weight of each ECC brick is approximately 8.3 kg, adding up to ∼ 8 tons for the whole SND. The bricks are then assembled to shape 19 walls of ∼ 50 units each, alternated with planes of electronic detector, called Target Tracker (TT), planes. For the time being, we consider the SciFi detector [45] as a feasible TT technological option. The TT additionally provide the time stamp of the event and help in linking the emulsion tracks to those reconstructed in the spectrometer and the muon system downstream of the SND. These features make the SND perfectly tailored for neutrino physics using all three flavours, as well as detection of light dark matter particles scattering off of electrons and nuclei of the SND.
An approximately 50 m long vacuum decay vessel is positioned downstream of the SND. The proposed facility is completed with a Hidden Sector Detector (HSD), equipped with calorimeters and muon detectors for the identification of long-lived Beyond Standard Model (BSM) particles.

Light dark matter production and detection
At a proton beam dump, DP can be abundantly produced in several channels: 1. Light meson decay: proton collisions on a heavy target result in a copious flux of outgoing mesons. Hence, DP may be produced in radiative decays of neutral mesons,

JHEP04(2021)199
Simulation Software  whereas a final state photon converts into a DP. The production cross-section is proportional to 2 and the relevant contributions come from the lightest mesons, because of decay modes with photons with relatively high branching ratio: π 0 , η, ω [15].

2.
Proton bremsstrahlung: being a charged particle surrounded by its own electromagnetic field, the proton radiates low-frequency and/or quasi-collinear photons with high probability when it undergoes a scattering process. Vector states like DP mediators can then be generated via radiative process p A → p A A [50] in proton interactions with the target nuclei. 3. Direct perturbative QCD production: it corresponds to the dominant production mechanism for higher masses (m A 1 GeV). At the lowest order in the strong interaction, DP are produced through the quark-antiquark annihilation process qq → A [15]. At higher orders, one can also have the associated production with a jet, according to the quark-gluon scattering process qg → qA , and with multiple jets.

JHEP04(2021)199
In addition, secondary leptons produced in the dump can contribute to the flux of photons, and thereby of DPs, by different re-scattering processes occurring within the target. Such lepton-induced processes are usually sub-dominant at a proton beam dump. However, they are not completely negligible, as nicely shown in a the dedicated analysis [51], and might be relevant in scenarios in which the New Physics does not couple with coloured particles. We do not include them in this work. Therefore, our estimates should be considered conservative in this regard. The minimal DP model can be probed by the SHiP experiment through the direct detection of LDM elastic scattering process off of the electrons and nuclei of the SND (figure 4) For the majority of the events χ e − → χ e − , the scattered electron is sufficiently energetic to generate an electromagnetic shower within the brick. Given that the ECC device acts as a high-granularity sampling calorimeter, it is thus possible to reconstruct the electron and measure its energy. Furthermore any activity in the proximity of the primary vertex can be spotted down to 100 MeV or below, thanks to the micrometrical position resolution of the nuclear emulsion device and the high sampling rate. These features translate into capability to accurately identify and tackle background events to LDM searches, as further described in section 5. As a consequence, LDM scattering events can be distinguished from a large neutrino-induced background.
An estimate of the order of magnitude of the expected yield of LDM interactions at SHiP can be determined as follows. The number of LDM-electron scattering events in the SND detector is given by the standard formula where N e − is the numbers of scattering centres, i.e. the electrons in the detector, φ is the flux of incident LDM particles and A SND represents the transverse area in (x, y) of the JHEP04(2021)199 SND. The elastic LDM-electron scattering cross section is roughly given by The flux φ mainly depends on the specific value of the DP mass which in turn determines the relative importance of the different production mechanisms. For example, for m A m π , LDM production in the beam dump is dominated by pion decays. In this case and under the assumption that all the primary proton impinging on the target will eventually interact in the beam dump, φ can be written as is the total number of p.o.t. delivered in the five years of data-taking; λ π 0 denotes the multiplicity of π 0 s per p.o.t.; A geo embeds the geometrical acceptance of the SND to LDM interaction vertices, corresponding to an angular coverage |(θ x , θ y )| ≤ (12, 10) mrad from the proton beam dump. If considering an average value of λ π 0 ∼ 6 as provided by the simulation of prompt proton-nucleon collisions with Pythia 1 [47], a geometrical acceptance A geo ∼ 30% and if assuming a coupling close to the current experimental constraints ∼ 5 × 10 −5 for a 10 MeV-DP and α D ∼ 0.1, the expected number of LDM candidates foreseen in SHiP is ∼ 1.3 × 10 4 .
We used MadDump [52] as the principal tool for the simulation of signal events. Its general philosophy and all the technical details are outlined in ref. [52]. We generate the event samples at the particle level and apply the selection criteria on the recoil electrons without taking into account other detector effects besides the geometrical acceptance. This is consistent with what has been done in the estimate of the background event rate. Since the target length is way larger than the proton interaction length in the material, we assume all of them to interact within the beam dump. In the following, we give further details for each production mechanism.

Meson decay
The relevant parameter space within the reach of the SHiP SND corresponds to m A > 2m χ and g D e. Indeed, in this scenario, the DP decays almost entirely into DM after travelling a very short distance, maximising the DM flux reaching the SHiP SND. The

JHEP04(2021)199
decay rate for light mesons decaying into dark photons is then dominated by the formation of an on-shell dark photon which decays promptly into dark matter, BR(A → χχ) 1. The production process is then well described by an effective Lagrangian with mesons as dynamical degrees of freedom leading to interaction vertices like XγA (X = π 0 , η, see figure 5) and ωπ 0 A . The corresponding branching ratios scale with 2 and are given by: An interested reader can find useful insights about the formulas above in [16,53,54]. The full simulation process is performed in three steps: i. production of the input meson fluxes originating from the incoming protons impinging and interacting within the target (beam dump); ii. generation of DM fluxes from the BSM meson decays in the relevant DM mass range; iii. generation of the corresponding DM − e − scattering events within the detector acceptance and the selection criteria.
MadDump provides a unified framework to handle the last two steps, in which all the new physics content is involved. The main source of uncertainties comes from the meson fluxes. Indeed, the description of the proton-nucleus interactions is highly nontrivial and experimental data are available only for a limited collection of energies and nuclear targets. Phenomenological and data-driven parametrisations for the distributions of the light mesons have been proposed in the literature [55]. An alternative strategy is provided by Monte Carlo programs like Pythia [47]. Recently, Pythia(8) results have been compared with existing experimental data for the inclusive production of π 0 and η mesons in pp collisions [56]. A fairly good agreement has been found for mesons with high momentum and within middle-high rapidity range (where the Feynman variable 0.025 < x F < 0.3), which represent the bulk of our events in acceptance.
Furthermore, secondary interactions of hadrons in the beam-dump target may affect the particle multiplicities, which in turn may increase the LDM yields and impact the sensitivity reach of the experiment. It is thus important to estimate the so-called cascade effects [57]. As the main input for the lightest mesons (π 0 , η) we use the fluxes generated with GEANT4(v10.3.2) within the FairShip software framework, which takes into account the secondary interactions adapting what has been used in ref. [58]. We also consider samples of mesons from primary proton-nucleon interactions generated with Pythia, as a reference to assess the impact of the cascade. For the ω, we rely on the Pythia samples only.
In tables 2, 3 and 4, we report a selection of results for π 0 and η, comparing the FairShip and Pythia samples. An important parameter in the FairShip simulation is the energy cutoff E cut applied to the particles produced at each interaction vertex: particles with energy less than E cut are removed from the list of those considered for new interactions within the JHEP04(2021)199

FairShip
Pythia . Therefore, we observe that secondary interactions occurring within the beam-dump target greatly increase the particle multiplicities and, in turn, lead to a rise of the DM yield by the same amount. However, this does not translate directly into an enhancement of the signal yield in the SND. Indeed, in order to produce a detectable scattering event one should take into account • the geometrical acceptance, • the path travelled within the detector, • the cross section for the scattering process.
We consider separately the effect due to the geometrical acceptance, defining an effective number of mesons per p.o.t. N eff X /p.o.t. as the average number of mesons of species X per p.o.t. which produce a DM particle impinging on the detector surface. For different m A values, we report in tables 3 and 4 our estimate of N eff π 0 /p.o.t. and N eff η /p.o.t. as estimated with Pythia and FairShip. The comparison shows that the increase due to the cascade is around 50 − 70%. The explanation is that the secondary particles mainly populate the soft part of the spectrum, as it is clearly shown in the left panels of figure 6 and figure 7 which have to be compared with the corresponding right panels describing the spectrum from prompt yields. Moreover, the cross section for the elastic LDM-e − scattering grows with the energy of the incoming dark-matter particle before saturating to a constant behaviour [59]. Hence, we expect that scattering events initiated by LDM particles produced in secondary interactions, being softer, will be less probable. This is clearly demonstrated by the last two columns in tables 3 and 4 in which we report the final signal yields N s (corresponding to the benchmark point α D = 0.1 and m A = 3m χ and = 10 −4 ) due to FairShip and Pythia samples respectively. From the comparison, we see that the impact of the secondary interactions is reduced to that given by the geometrical acceptance only. In conclusion, our finding is that for π 0 , the cascade modestly affects (∼ 15 − 40%) the signal event yields within the detector volume, while for η it is negligible.

Proton Bremsstrahlung
In the mass range 500 MeV m A 1 GeV, the production of A is dominated by the proton bremsstrahlung mechanism. The photon emission cross section is indeed enlarged in the collinear direction so that a sizeable fraction of A is produced within the geometrical acceptance for an on-axis detector as the SND (∼ 20%). In this limit, the process can be described by a generalisation of the Fermi-Williams-Weizsäcker method [60][61][62], based on the assumption that the p − N scattering is dominated by the exchange in the 1 −− channel. We extend MadDump include the bremsstrahlung from the primary protons. Following refs. [50,63], we parametrise the four-momentum vector of the emitted A as p A = (E A , p T cos(φ), p T sin(φ), zP ), with E A = zP +(p 2 T +m 2 A )/(2zP ), where P is the momentum of the incident proton, z is the fraction of the proton momentum carried by the outgoing A , p T is the momentum perpendicular to the beam momentum and φ is the azimuthal angle. We generate unweighted A events according to the differential production rate where s = 2m p (E p − E A ), s = 2m p E p and the photon splitting function is In the above formula, F 1,p is the time-like proton formfactor, as provided by the parameterisation in ref. [64]. It effectively incorporates off-shell mixing with vector mesons such as ρ and ω, corresponding to a resonance effect around m A ∼ 770 MeV [65]. In ref. [63], the authors compare the description of the peak region JHEP04(2021)199 by adopting the time-like proton form factors and by adding by hand the vector mixing within an on-shell computation, finding small deviations in the peak region. Assessing the size of this uncertainty is beyond the scope of this work.
The next steps of the simulation, namely the decay A → χχ and the χ − e − scattering in the SND, are handled by standard MadDump functions. The whole process has been integrated into the new MadDump mode bremsstrahlung-interaction.
The normalisation of the flux of the original A is given by the integral of the differential production rate eq. (4.6) in the validity range of the equivalent photon approximation, given by the kinematical conditions Following refs. [50,63,66], we adopt z ∈ [0.1, 0.9]. For a relatively high energy experiment such as SHiP, the minimum DP energy E P corresponds then to ∼ 40 GeV and we can set its maximum transverse momentum p T to 4 GeV, i.e. an order of magnitude less. We expect electron bremsstrahlung to be sub-dominant as discussed for example in [66,67]. As for the cascade effects, extra dark photons may arise from the bremsstrahlung of secondary charged hadrons. Similarly to what happens in the case of mesons, the picture is complicated by the impact of the geometrical acceptance and the convolution with the scattering cross section. For the case the proton undergoes a chain of elastic proton-nucleon collisions, so that it retains all of its initial energy, we can make a rough estimate by means of the following simplified calculation. Let p el be the probability that the incoming proton undergoes an elastic scattering interaction with a nucleon in the target and p brem the probability of a dark photon produced in the proton bremsstrahlung. Under the assumption that p brem does not depend much on the number of previous elastic collisions, the probability that a dark photon is produced in this chain is At the energy of SHiP, p el 0.25 so that we estimate a mild increment of ∼ 30%. In the following, we keep the conservative estimate based only upon the bremsstrahlung of the primary protons.

QCD prompt production
QCD parton processes become relevant for m A 1 GeV, at the edge where perturbation theory starts to become reliable. Indeed, at scales 1 GeV the strong coupling α s is O(1) and the description of the hadrons in terms of constituent partons is spoiled by the confinement. In the attempt of estimating the relative importance of this production mechanism, we have tried to push the perturbative computation down to m A ∼ 300 MeV. The main tree-level diagrams correspond to the Drell-Yan-like production and the associated production with QCD radiation, figure 8. The latter allows for smaller m A values since the characteristic scale of the process, given by the transverse momentum of the emitted parton, can be kept to be higher than the Λ QCD scale. A minimum cut on the JHEP04(2021)199 p T of the QCD radiation is physically required to tame infrared singularities. The cross section diverges logarithmically up to scales of order O(Λ QCD ), when perturbation theory eventually breaks down. The transverse momentum cut-off is a severe requirement for an on-axis set-up as SHiP, due to its small angular acceptance. We find that even for relatively small values of the cut-off, p T ∼ 800 MeV, the production rate is not sufficient to produce a significant yield of LDM within the geometrical acceptance. Therefore, we focus only on Drell-Yan-like production. We rely on MadGraph5(v.2.66) [68], which is integrated in MadDump as it is based on the former package, for the generation of the events, and we use the NNPDF2.3LO [69,70] set as our choice of the proton parton distribution function (PDF). In the normalisation of the number of produced LDM particles, we effectively take into account nuclear effects in the following way where A = A Molybdenum = 96; the nuclear rescaling as A 0.71 is taken from ref. [71] and σ pp→all = 40 mb [72]. In this case, the characteristic scale of the process coincides with m A . As mentioned before, we cannot use scales 1 GeV, where both the strong coupling and PDF are illdefined from the perturbative point of view. To push our projection into the sub-GeV range, we adopt the following prescription: we set both the re-normalisation scale µ R (at which the strong coupling constant is evaluated) and the factorisation scale µ F (at which the PDF is evaluated) to a fixed value chosen to be µ R = µ F = 1.5 GeV, the lowest scale variation point associated to open charm production.

Background estimate
Neutrinos emerging from the beam-dump target and interacting in the SND are the relevant background source to the detection of LDM elastic scattering, whenever the topology at the primary vertex consists of a single outgoing charged track, an electron. The expected background yield for five years of data-taking has been estimated by means of the GENIE [49] Monte Carlo software, supplied with the spectrum of neutrinos produced at the beam dump as simulated with Pythia v6.4.28 within FairShip and including secondary production, for the generation of the following neutrino interactions in the whole kinematic phase space:

JHEP04(2021)199
• Elastic scattering (EL) of ν e (ν e ), ν µ (ν µ ) off the electrons of the SND, which is a source of irreducible background as it shares the same topology of LDM elastic interactions: • • Quasi-elastic scattering (QE) of ν e (ν e ), with the primary proton undetected because it is below the energy threshold: Charged current interactions of ν (ν ) with = µ, τ do not represent a concern because they are easily discernible from LDM events by reconstructing the charged lepton produced in the final state. Electron decay modes of the τ lepton are a negligible background source, since an early decay of the parent track (∼ 1% occurrence) leading to an undetected τ would occur with less than a per-mill probability. In addition, we do not consider ν τ (ν τ ) elastic scattering processes as background, due to the suppression resulting from the combination of smaller flux φ ντ (∼ 1 order of magnitude smaller than φ νe and ∼ 2 orders of magnitude smaller than φ νµ ) and cross section. The whole ν spectrum is made to interact within the SND and the surrounding materials. Moreover, for this study we assume the detection efficiency to be unitary [73].
The simulated sample of neutrinos undergoes a two-steps selection procedure, in order to be tagged as residual background.
First, only interactions occurring within geometrical acceptance and associated with a single charged final state track, an electron, are selected: ν vertices are further considered in the analysis only if located inside the SND volume, whereas all the out-coming charged tracks are inspected in order to assess their visibility in the nuclear emulsion medium. The visibility threshold depends crucially on the exploited tracking device technology; for this study we assume 170 MeV/c for the protons, 100 MeV/c for the other charged particles including the electrons. These are derived as benchmark values from the OPERA experiment, where charged-particle reconstruction is possible only if two consecutive straight track segments, before and after a lead plate, are found to be in agreement [74]. A further handle considered here for signal against background discrimination is the presence of neutral particles, e.g. photons or π 0 s, nearby the interaction vertex, since it is not foreseen in any LDM elastic scattering event.

JHEP04(2021)199
The second step of the event identification procedure consists of a kinematic selection in the energy E e and polar angle with respect to the incoming neutrino/LDM direction θ e of the scattered electron. For the elastic case, these quantities are constrained by the kinematic relation E e θ 2 e ≤ 2 m e , valid in the regime E in m e , m χ , where E in is the energy of the incident neutrino/LDM particle. In order to choose the energy and angle ranges for the selection, an optimisation procedure is performed, aiming at maximising the following significance: where S denotes the signal yield, while B i are the individual contributions to the background yield B per interaction category and neutrino flavour, each of them weighted by a factor κ i accounting for the systematic uncertainty. We have focused on the relevant systematics, arising from the uncertainty on the neutrino cross sections (assumed flavour-independent, κ i ) and on the incoming neutrino flux produced at the beam dump (interaction-independent,κ ), so that we have assumed κ i = κ 2 i +κ 2 . As for the former, we assume the following: 5% for DIS [75], 18% for RES [76], 8% for QE [77], while we neglect the uncertainty on the EL cross section that is well known within the SM [78]. As for the uncertainty on the incoming neutrino flux, this will be well constrained by an independent measurement of the abundant CC-DIS interactions occurring within the SHiP detector (expected ∼ 10 6 for ν e, µ ). Since the corresponding cross section is lepton-universal and known within ∼ 5% accuracy down to E ν of 2.5 GeV [75], we assume it to be the driving systematic uncertainty on the neutrino flux. While SHiP is capable of disentangling ν µ fromν µ interactions by measuring the charge of the primary muon, thus providing a different estimate for ν µ andν µ fluxes, with regard to electron species it will measure a combination of the lepton and anti-lepton initiated events. As the relative abundance of ν e andν e produced at the beam dump can be assessed, the individual fluxes can be estimated accordingly. For neutrino energies below 2.5 GeV we double the uncertainty on the flux assuming them to be at 10%.
Since the signal yield S depends on the mass hypothesis placed on the LDM candidate and thus on the DP, we adopt the most-general assumption of maximising the experimental sensitivity with respect to the broadest possible range of masses. Therefore, S is given as the average of the signal yields for three DP mass hypotheses: 50 MeV, 250 MeV and 500 MeV.
The selection optimisation strategy is based on a grid-search method and proceeds as follows: • an energy window [E min , E max ] is identified, according to the signal events distributions; • in the given energy range, the significance Σ values are determined in uniform angular intervals of 5 mrad spread; • the selection ranges, corresponding to the maximum Σ, are chosen for the analysis.  As shown in figure 9, signal events are mostly concentrated at energies below 10 GeV. Two energy windows have thus been considered: [1,5] GeV and [1 , 10] GeV, where the lower cut is placed as a minimum requirement for the recoil electron to produce a detectable electromagnetic shower within the ECC brick. The motivation to consider an additional tighter energy range resides in the opportunity to further suppress the high energetic components of the neutrino background, as illustrated in figure 10 which shows the relevant EL and QE contributions. DIS and RES processes are not shown since they exhibit signatures with higher multiplicities of charged tracks at the primary vertex.
The results of the optimisation are reported in figure 11, showing indeed a preference for the tighter energy window E e ∈ [1,5] GeV and an angular range θ e ∈ [10,30] mrad.
The corresponding background yield estimate is reported in table 5.
Neutrino elastic scattering processes, involving either electronic and muonic species, represent the dominant background source and are to some degree irreducible, since they share the same topology as the signal.
With regard to quasi-elastic ν e andν e interactions, a small but non-negligible contribution is observed. The process ν e n → e − p mimics the signal when the proton at the primary vertex is not identified, because of the 170 MeV/c threshold. Improvements in the proton identification efficiency with dedicated techniques, including Machine Learning clustering algorithms, will be the subject of future studies. When considering anti-neutrinos, events asν e p → e + n are topologically irreducible since we assume for the present study the neutron to be undetectable within the SND. This effect compensates the larger (by a factor of ∼ 3) neutrino flux, thus making the two contributions comparable.
In the case of resonant neutrino scattering, the outgoing electron is often accompanied by a further charged track, which helps discriminating between background and signal.  Nevertheless, some topologically irreducible interactions are present as well: where the K 0 L/S is considered undetectable within the SND for this study. Future improvements lie in the employment of combined information of ECC and TT, coming from the linking of the emulsion tracks with those reconstructed in the electronic tracking system. Moreover, some final states with the pattern e + (n) γ contribute, when the emitted photon is too soft to be identified via the reconstruction of the electron-positron pairs from its conversion.
The contribution from neutrino deep inelastic scattering processes is, on the contrary, negligible, as a consequence of the high rejection power observed on these event categories, which exhibit a topology with a high multiplicity of charged tracks.    Figure 11. Grid-search optimisation of the significance Σ as a function of the angular cut for a fixed energy window. The left axis represents the lower cut value for θ e whereas the upper axis is the higher one. The plots in the two panels share the same normalisation. The best selection corresponds to the tighter energy window (a) and the angular range [10,30] mrad.
In the eventuality of an observed excess in the number of events, SHiP may collect data in a bunched beam mode, exploiting the time of flight measurement to separate massive particles like LDM from neutrinos.

Sensitivity
Once the significance of eq. (5.1) is maximised, the optimal energy and angle ranges are employed to determine the yields of signal and background, following a cut-and-count procedure, per each fixed value of the mediator mass m A . The 90% confidence level (C.L.) exclusion limits on the coupling at fixed mass m A are then retrieved by adopting JHEP04(2021)199 a single-tail Poissonian statistics. Statistical and systematic uncertainties are combined as reported after eq. (5.1).
In figure 12, we report our projection for the SHiP SND exclusion limit at 90% C.L. in the m χ − Y plane of the dark-photon model. As stated above, we consider the benchmark scenario α D = 0.1 and m χ = m A /3. In figure 13, we separate the contributions given by the different production mechanisms. In the low mass range m χ 150 MeV, the main contribution comes from the decay of the lightest mesons. π decays dominate the A yield up to masses close to the m π → γA kinematic threshold. When approaching this threshold, the decay rate rapidly closes due to the steep suppression given by the phase space factor and with further increasing m χ mass the η → γA starts to dominate. The contribution of the ω is subdominant in the whole available mass range, which justifies a posteriori the fact that we do not include in our analysis A production from decays of heavier meson like the η .
We find that the contribution due to pQCD is very small in the mass region explored. By varying the factorisation scale in the range 800 MeV < µ F < 3 GeV, we estimate the uncertainty associated with missing higher orders to be about 15% on the signal yield within acceptance. We believe that this is an underestimation of the uncertainty as at next-to-leading order the process starts to receive radiative corrections proportional to the strong coupling constant at a scale close to Λ QCD , and new production channels open. While we do not expect that this will lead to a sizeable impact on the sensitivity, neglecting it leads anyway to a conservative estimate of the signal; hence, we have not considered the contribution of pQCD in our final result.
In the mass range 1 MeV < m χ < 300 MeV, the SHiP upper limit fairly improves the current strongest experimental limits (BaBar [13], Na64 [27]), even by more than an order of magnitude in the central region (5 MeV < m χ < 100 MeV). In this range and for the benchmark point under investigation, SHiP will cover the still unexplored parameter space corresponding to the solution of the relic density given by a scalar LDM. In the range 3 MeV < m χ < 300 MeV, SHiP will reach the thermal target for a Majorana candidate. Furthermore, it will exceed the thermal target for a Pseudo-Dirac candidate for masses around 10 MeV < m χ < 40 MeV.
We notice that for m χ 5 MeV the SHiP line saturates. In this region, the dark matter mass starts to become negligible and the selection requirements affect similarly the signal and the background. The rise in the signal production rate due to a lower mass is then balanced by a smaller fraction of events passing the kinematics selection, leading to the observed flat sensitivity in the small mass range. The distinctive peak at m χ 257 GeV corresponds to the ρ−ω resonant region, which is effectively taken into account by the timelike proton form factors used in the modelling of the proton bremsstrahlung mechanism.
In figure 14, the comparison between the SHiP sensitivity reach and that of other concurrent experiments clearly shows strengths and the complementarity offered by the proposed experimental scenario. Indeed the SHiP experiment will place constraints in unexplored regions of parameters space by exploiting a high intensity proton beam dump at 400 GeV and a micrometrical resolution tracking capability with the ECC. Thus, it offers a diverse approach to this NP search with respect to other experimental scenarios including direct searches and electron beam-line technologies.

JHEP04(2021)199
Relic Density SHiP Figure 12. SHiP SND exclusion limit at 90% CL relative to a A decaying into χχ pairs for the benchmark point α D = 0.1 and m A = 3 m χ . The current strongest experimental limits are also shown (BaBar [13], NA64 [27]), together with the three thermal relic lines corresponding to the scalar and the Majorana [3], and the Pseudo-Dirac DM [79] hypothesis.

Conclusions
Light dark matter particles χ with masses in the sub-GeV region represent an appealing scenario for the explanation of the observed thermal relic density in the Universe. In this work, we have studied the potential offered by the SHiP SND to reveal LDM which couple to SM particles via a new gauge force mediated by a vector boson, A . We have assumed the simplest DP model, with coupling g D to χ and A kinetically mixed with the SM photon with mixing parameter . We have focused on the relevant scenario for the SHiP SND: m A > 2m χ and g D e. Our main result is that for DM masses in [1,300] MeV the SHiP experiment will reach an unexplored region of the parameter space. For the benchmark point considered, the sensitivity of the SHiP SND is even below the thermal relic line corresponding to a Majorana DM candidate in the mass window [3,300] MeV and it will reach the target for a Pseudo-Dirac candidate within [15,30] MeV. Our analysis is based on a robust simulation framework for both the signal and the background which includes the relevant physical processes propagated within the detector. In particular, interactions of secondary particles in the beam-dump target have been taken into account in the neutrino background modelling, assuming unitary detection efficiency. As for the signal, we have consistently adopted 100% detection efficiency within the selection requirements and we have studied the impact of the cascade effect on meson multiplicities and the resulting dark matter production yield. We have observed that the impact of the cascade is quite modest and affects mainly the low masses.
In this work, we have focused on the elastic χ e → χ e signature detectable within the SHiP Scattering and Neutrino Detector. Other signatures, as the elastic scattering with nuclei, may lead to an improvement of the sensitivity. We leave their study to forthcoming works. In our case, the main background sources arise from elastic ν /ν -electron and quasi-elastic ν e /ν e scattering. We have considered the region E e ∈ [1,5] GeV and θ e ∈ [10,30] mrad, where e is the recoil electron. We have found that about 230 neutrino events survive the selection requirements, for 2 × 10 20 p.o.t. corresponding to 5 years of datataking.
We conclude by mentioning that, should an excess of events be observed, a time of flight measurement of particles scattering within the SND might represent a smoking gun to discriminate LDM from neutrino events, thus leading to an inarguable discovery.