Extending the reach of FASER, MATHUSLA and SHiP towards smaller lifetimes using secondary production

Many existing or proposed intensity-frontier search experiments look for decay signatures of light long-lived particles (LLPs), highly displaced from the interaction point, in a distant detector that is well-shielded from SM background. This approach is, however, limited to new particles with decay lengths similar or larger than the baseline of those experiment. In this study, we discuss how this basic constraint can be overcome in BSM models that go beyond the simplest scenarios. If more than one light new particle is present in the model, an additional secondary production of LLPs may take place right in front of the detector, opening this way a new lifetime regime to be probed. We illustrate the prospects of such searches in the future experiments FASER, MATHUSLA and SHiP, for representative models, emphasizing possible connections to dark matter or anomalous magnetic moment of muon. We also analyze additional advantages from employing dedicated neutrino detectors placed in front of the main decay volume.

Motivated by a successful history of discoveries of new elementary particles, it has long been the dominant paradigm in experimental searches to look for heavier and heavier particles that could manifest their existence in increasingly more powerful colliders. This approach led to a well-established experimental program that is now being continued at the Large Hadron Collider (LHC). In parallel with these persevering efforts, however, there has also been a growing interest in recent years in exploring scenarios that might have been overlooked in the past.
Among them, a prime focus has been directed on light and very weakly-interacting beyond the Standard Model (BSM) particles that could have escaped detection in previous searches due to the lack of sufficient luminosity. The corresponding efforts are often referred to as intensity frontier searches for light long-lived particles (LLPs). This captures the fact that the relevant detection prospects depend upon the ability to study very rare events that should also be discriminated from a priori overwhelming Standard Model (SM) background (BG).
A variety of atypical experimental signatures have been proposed to search for LLPs [1,2]. In particular, very clean searches for LLPs have been considered to employ their displaced decays in a distant detector that is physically separated and shielded from the primary interaction point (IP). This allows one to drastically reduce BG, often to negligible levels. However, an obvious limitation of this strategy is the lifetime of BSM particles that can be probed this way. If a LLP is too short-lived, it primarily decays before reaching a detector's fiducial volume in which case prospects for discovery are greatly reduced. This is evident in a variety of simplified models that have been proposed as benchmark scenarios for intensity frontier studies [1]. In particular, in models employing a single coupling of an LLP to the SM particles, the shape of current exclusion bounds and future sensitivity lines is often driven by the tension between increasing production rates and diminishing lifetime of the LLP. This can be overcome in less simplified models employing e.g. non-universal couplings of the LLP to SM hadrons and leptons, or by the use of compressed mass spectrum to increase the lifetime of decaying LLP. By tuning different types of couplings or a mass splitting between the particles present in the model, one can keep relatively high production rate of LLPs, while simultaneously increasing their decay length.
On the other hand, in more complete models, if an LLP is present, it is often accompanied by other light BSM species. A notable examples can be light dark matter (DM) with a mediator particle of similar mass that couples to the SM and yields correct thermal DM relic density [3,4], or various realizations of the twin Higgs model [5], that have been advocated for in the context of neutral naturalness and often predict the existence of several light BSM particles in the mirror sector.
In such scenarios, the aforementioned effective decoupling of production and decay of the LLPs can appear even more naturally, without requiring a tuning between the model parameters or introducing exotic non-universal couplings. In particular, if one of the LLPs can be effectively produced in interactions of the other, that happen in front of the detector or inside it, this can lead to secondary production of LLPs at a position much closer to a decay vessel. We illustrate this schematically in Fig. 1.
As a result, the range of the LLP lifetimes, τ LLP , that can be probed is extended towards smaller values. In addition, for a range of intermediate values of τ LLP , both primary and secondary production give comparable con-tributions to the total number of events in the detector. We illustrate this in Fig. 2 for selected intensity frontier experiments and benchmark points (BPs). In the figure, we show the expected number of events associated with both the primary production at the interaction point (IP) and secondary production mechanisms as a function of τ LLP .
In this study, we analyze the prospects for such searches in some of the future experiments that have been proposed to look for the LLPs and could naturally employ the aforementioned secondary production mechanism to increase their discovery potential. Among them, we study the FASER experiment [6][7][8][9] that will be taking data during LHC Run 3 and can then be extended towards the High Luminosity LHC (HL-LHC) phase. We also discuss the prospects of two other experiments with a possible timeline similar, or identical, to HL-LHC, namely the proposed MATHUSLA [10][11][12] and SHiP [13][14][15] detectors.
In the case of the proposed SHiP detector, we additionally discuss two other possible signatures of models with more than a single LLP that employ a dedicated neutrino subdetector. Interactions of LLPs can be discovered there based on their scatterings off electrons, as well as such scatterings followed by a prompt decay in the subdetector that leads to two simultaneous electromagnetic (EM) showers. We analyze how a combination of various signatures can deliver important information, on top of the one from the standard search for LLP decays in the decay vessel, that could help to resolve the nature of the LLP. This paper is organized as follows. In Sec. II, we introduce simplified BSM models of interest to us. In Sec. III, we discuss basic aspects of the LHC experiments under study, while more details of our analysis are given in Sec. IV. The results for the secondary LLP production in scatterings off nuclei are presented in Sec. V, while additional signatures for scatterings off electrons in the SHiP neutrino subdetector are discussed in Sec. VI. We conclude in Sec. VII. More technical aspects of the study are described in appendices. In Appendix A we discuss how the sensitivity reach from the secondary production can be affected by varying cuts imposed in the analysis. In Appendices B and C, respectively, we present the expressions relevant for the primary production and decay of LLPs, as well as the scattering cross sections corresponding to the secondary production. Appendices D and E are devoted to a more detailed discussion of the detector designs used in modeling, as well as relevant experimental cuts.

II. MODELS
In order to illustrate the impact of secondary production mechanism on sensitivity reach in intensity frontier searches, we consider some simplified BSM models that can lead to production of LLPs in scatterings in front  Figure 1: Schematic illustration of primary and secondary production of the light long-lived particle LLP 2 , shown with the dashed blue line, that subsequently decays inside the detector. In case of secondary production, intermediate role of other light species, LLP 1 , shown with the dotted red line, allows for producing LLP 2 much closer to the detector. In models studied below, the (LLP 1 , LLP 2 ) pair can be identified with: (χ, A ) for dark bremsstrahlung (see Sec. II A), (χ 1 , χ 2 ) in case of inelastic dark matter (see Sec. II B) and (S, A ) for the model with dark photon and secluded dark Higgs boson (see Sec. II C). and SHiP (right) coming from primary (red) and secondary (blue) production of a light long-lived particle. For the purpose of illustration, the plots were obtained for benchmark points indicated in the captions that correspond to the models described in Sec. II. The currently excluded ranges of lifetimes for these BPs are shown as grey-shaded regions.
of the detector and their subsequent decays within its fiducial volume.
To this end, we focus on popular scenarios with a light dark vector particle appearing in minimal extensions of the SM with a new U (1) d symmetry group. We assume that the resulting gauge boson, dubbed dark photon or A , couples to the SM solely via a kinetic mixing term ( /2) F µν F µν , where F µν and F µν are the field strength tensors of the SM photon and dark photon, respectively. The parameter can be naturally small when induced at a loop level due to the presence of new heavy charged particles [16]. After a field redefinition to remove the nondiagonal term in the field strength tensors, SM fermions acquire milli-charges under the U (1) D group with the corresponding interaction mediated by the dark photon (for recent reviews see [17][18][19]). The relevant Lagrangian terms for the dark photon with mass m A then read where the sum in the second term spans over SM fermions f with electromagnetic charges q f . In the following, we will focus on scenarios with m A ∼ MeV−GeV. In this mass range, dark photons are promising targets for intensity frontier searches and one of just a few renormalizable BSM portals to study at a simplified level. In particular, the values of the kinetic mixing parameter in the allowed region of the parameter space of the dark photon model spanned by just two parameters, m A and , can be as large as ∼ 10 −3 and lie within reach of many current and future experiments [1,2].
Once produced, dark photons can decay back into SM particles with the relevant decay length given in Eq. (B6). If dark photons are produced in secondary production processes right in front of the detector, this can allow one to probe boosted A s with the decay length d A of order meters Additional LLPs with mass of a similar order often arise in such models, e.g. as dark sector particles comprising DM, or in connection to the dark Higgs mechanism that can generate non-zero mass of A . In such scenarios, dark photons can either decay visibly or decay predominantly into dark sector particles. Below, we briefly discuss several simple scenarios with A accompanied by additional LLPs.

A. Dark bremsstrahlung
One of the most important motivations to search for new light subweakly coupled particles is the role they can play in cosmology and astrophysics, acting as mediators between the SM and DM (for a recent review see e.g. [20]). In particular, light dark photon serves as an important example of such a portal that can yield correct relic density of thermally produced DM due to a generalized WIMP miracle [3,4] or in secluded WIMP DM scenario [21], as well as due to non-zero temperature effects in the forbidden DM regime [22] (see also Ref. [23] for recent study in this direction). In addition, depending on the hierarchy between DM and dark photon masses, this scenario can belong to a more general framework of self-interacting DM and can help to resolve problems in understanding DM distribution towards center of galaxies and galaxy clusters [24] (for review see [25]).
We focus on the model with fermionic DM χ coupling to SM via a dark photon portal described by Eq. (1) and the following additional terms in the Lagrangian where D µ = ∂ µ − i g D A µ , m χ is the DM mass and g D is the dark coupling constant that governs DM interactions with A . 1 As a result, the parameter space of this model is spanned by four parameters, m A , m χ , and α D = g 2 D /(4π). In the case of a massive dark photon, that we focus on, dark fermion χ remains electrically neutral after gauge field transformations are applied [26].
In the following, we set α D = 0.1, a value that lies within perturbativity limits but at the same time is large enough so that α D -dependent secondary production of LLPs can become sizeable. For illustrative purposes, we also assume a fixed mass ratio m χ : m A = 0.6 : 1. This corresponds to a particularly interesting mass regime in which m χ < m A < 2 m χ . Here, dark photons decay visibly into SM particles with possible striking experimental signatures in distant detectors. At the same time, DM relic density is set by a freeze-out of direct χχ annihilations into SM particles via intermediate A [27], with only a small impact of very efficient annihilation into two dark photons relevant for the forbidden DM regime. Instead, in case of a larger mass ratio, m χ > m A , characteristic for secluded DM, one would typically obtain tiny DM relic density for a chosen value of α D . Further comments about DM relic density in this model are given in Sec. V. Once the aforementioned assumptions about the mass ratio and the value of α D are taken into account, there remain only two free parameters of the model: m A and .
While DM particles χ are stable, spectacular signatures of this model can come from dark photon decays inside the detector. Importantly, on top of A s produced at the primary IP, further dark photons can come from dark bremsstrahlung processes, χ T → χ T A . Here, χ scatters off electron or proton/nucleus target T in the material in front of the detector and radiates off the dark photon (see e.g. Refs. [28][29][30]). The relevant Feynman diagram is shown in the left panel of Fig. 3. Alternatively, one can study scattering signatures of DM particles off electrons, χ e − → χ e − , that can lead to excess in the number of high-energy EM showers in the detector with no significant nuclear recoil over expected neutrinoinduced BG [31][32][33][34][35][36].

B. Inelastic dark matter
Generalizing the aforementioned case of Dirac fermion DM, we also consider scenario with a pair of Weyl (dark) fermions χ L and χ R with opposite chirality that share the same charge under the dark gauge group U (1) D , see e.g. Refs. [37,38] for recent discussion. When U (1) D symmetry is spontaneously broken, Majorana mass terms can be generated on top of the Dirac mass. The relevant Lagrangian can be written as L ⊃ Ψ T C M Ψ + h.c. , with Ψ T = (χ L , χ c R ) and the (real) mass matrix where the Majorana masses are denoted by m L and m R and the Dirac mass by M χ . By rotating to the mass basis, two dark fermion states appear with masses The effective Lagrangian then reads where we focus on the case with g ii g 12 and the dark photon coupling to the SM as shown in Eq. (1). In fact, the diagonal and non-diagonal couplings between dark fermions depend on their mixing angle θ. They are proportional to cos 2θ and sin 2θ, respectively, with the mixing angle defined as tan 2θ = 2 M χ /(m L + m R ). Figure 3: Feynman diagrams relevant for secondary production processes discussed in the text. Left: dark bremsstrahlung of A in DM scattering, Center: upscattering of lighter χ 1 fermion into heavier one χ 2 in inelastic DM model, Right: scattering of secluded dark Higgs boson with A production.
As can be seen, in the pseudo-Dirac limit with M χ m L , m R , diagonal couplings are suppressed and two dark fermion states are characterized by a small mass splitting. This corresponds to a well-known scenario first discussed in the context of DAMA anomaly in DM direct detection searches [39], which has recently received renewed attention due to possible interesting signatures in intensity frontier searches and other experiments from displaced decays of a heavier fermion, χ 2 → χ 1 e + e − , see e.g. Refs. [40][41][42][43][44][45][46][47][48].
On the other hand, similar suppression of diagonal couplings can be achieved even for larger Majorana masses by requiring m L ≈ −m R . This then allows one to consider a larger mass splitting between χ 1 and χ 2 , while the mass eigenstates read m χ1,2 |m L ∓ M χ |. In particular, for m L −m R M χ , one obtains very low mass of χ 1 , i.e. m χ1 ≈ 0, while the χ 2 mass can remain significantly larger, m χ2 2 M χ . On the other hand, for χ 1 to be cold DM candidate, we require m χ1 to be not suppressed too much. This can be achieved by e.g.
assuming m L −m R 2 M χ leading to m χ2 3 M χ 3 m χ1 . In this case, for increasing mass of dark fermions, other decay channels might open up on top of the leading one to an electron-positron pair. These include e.g. χ 2 → χ 1 µ + µ − , as well as decays with hadronic particles in the final state.
While in the pseudo-Dirac limit χ 2 can become longlived due to suppressed mass spectrum in its 3-body decays, the lifetime of heavier dark fermion becomes smaller for increasing ∆ χ . As a result, χ 2 s often struggle to reach a distant detector before decaying. This opens up a region in the parameter space of the model in which correct DM relic density can be obtained, while current bounds from beam-dump experiments are weakened [38,42].
Due to suppression of diagonal couplings, lighter dark fermions, when scattering off the electron or proton target, preferably produce heavier state, χ 1 T → χ 2 T , if kinematically allowed. This is illustrated in the central panel of Fig. 3. If such upscattering occurs in front of the detector, subsequent decays of χ 2 might lead to a spectacular signature inside the fiducial volume of the detector. For large mass splitting between both dark fermions, an approximate decay length of boosted χ 2 reads while a full expression for the relevant decay width is given in Eq. (B7).
When presenting the results in Sec. V, we follow a simple mass scaling mentioned above with m χ2 ∼ 3 m χ1 and take both masses in the MeV − GeV range, which has been chosen for illustrative purposes. As the upscattering cross section decreases with growing dark photon mass, we additionally focus on the case when m A saturates the minimal value required for on-shell A to decay invisibly into χ 1 χ 2 pairs, i.e., we assume m χ1 : m χ2 : m A ∼ 1 : 3 : 4. In addition, similarly to Sec. II A, we assume α D = 0.1, where we define α D = g 2 12 /(4π). Hence, there remain only two free parameters of the inelastic DM (iDM) model that we vary freely when obtaining sensitivity reach plots: m χ1 and .

C. Dark photon with secluded dark Higgs boson
A natural way to introduce a non-zero dark photon mass is to employ a dark Higgs mechanism in which m A is driven by a vacuum expectation value (vev), v S , of a new SM-singlet complex scalar field S added to the model. In comparison to the SM Higgs boson vev v h , the vev of the new field is assumed to be small, v S v h , as expected for low-mass dark gauge boson. This new dark scalar, dubbed dark Higgs boson, can have non-negligible couplings to the SM fermions that either arise thanks to the mixing between S and the SM Higgs boson H, or at a loop level with the exchange of intermediate A s. The relevant Lagrangian terms read [31,49] where D µ = ∂ µ − i g D A µ , while dark photon coupling to the SM is given by Eq. (1). The phenomenology of new BSM light scalars in connection to intensity frontier searches has been extensively studied in the literature, see, e.g., Refs. [1,15,[50][51][52] and references therein. Assuming small mixing, λ SH λ S , and by solving relevant tadpole equations, one can rewrite the dark scalar mass in terms of v S which reads m 2 S . At the same time, the dark photon mass induced by the vev of S is given by m 2 i.e., both dark scalar and dark vector masses are naturally of similar order. In particular, keeping a small value of λ SH allows one to suppress a contribution to the dark Higgs boson mass from the vev of the SM Higgs, λ SH v 2 h . In the following, we require the mixing term to be very small, λ SH ∼ (m 2 S /v 2 h ) 10 −6 . This results in highly suppressed values of a mixing angle between the dark and SM Higgs bosons that typically lies below the reach of current and future searches. It then plays a negligible role in our phenomenological analysis, while we will comment on it when discussing current bounds on this scenario.
In fact, the dominant couplings of such a secluded S to the SM fermions arise via the dark photon portal and an unsuppressed coupling between S and A that appears after U (1) D symmetry breaking, L ⊃ g D m A S A µ A µ . This can lead to efficient co-production of light scalars in any process leading to A production, where S can be bremmed off the vector leg. As a result, a flux of dark scalars going towards the detector can be obtained along with dark photons produced at the primary IP.
Importantly, unlike to dark photons which can decay promptly depending on the value of the kinetic mixing parameter , dark Higgs bosons in such a scenario are typically very long-lived if m S < m A . This is because their dominant decays into a pair of SM fermions, e.g. S → e + e − , proceeds at a loop level through a triangle diagram with intermediate vector states [31,43,49]. A typical lifetime of S then reads For larger values of m S , the di-muon final state becomes possible, as well as 3-or even 2-body decays with one or two on-shell dark photons in the final state that reduce the lifetime of S. These, however, turn out to be irrelevant for the mass range of our interest and for the assumed mass ratio in our benchmark scenario m S = (3/4) m A , as discussed in Sec. V. We also fix the dark coupling constant to be α D = 0.1, which leads to only two free parameters of the model varied to obtain the sensitivity reach plots: m A and .
As evident from Eq. (10), dark scalars in the considered scenario are effectively stable at collider scales. Thus, once a flux of dark Higgs bosons is produced, they will only rarely decay before reaching the detector. Instead, when traveling, they can scatter off nuclei and electrons producing secondary dark photons, S T → A T , as illustrated in the right panel of Fig. 3. The dark photons can then decay inside the detector. If such a secondary A production can take place in the vicinity of the decay vessel, it allows one to probe small dark photon lifetimes.

III. EXPERIMENTS
In order to illustrate the impact of secondary production of LLPs on the reach, we study the expected sensitivity for three representative proposed experiments. Since this effect becomes more evident for highly boosted particles, we focus first on the LHC experiments, namely FASER and MATHUSLA. We also study the possible impact of secondary production of LLPs on the reach plots relevant for the SHiP experiment. 2 As discussed in more detail in Sec. IV and Appendix C, the dominant contribution to the secondary production rate comes from Z 2 -enhanced coherent scatterings off nuclei. They are characterized by a low-momentum transfer to the nuclear target. This is different from the neutrinoinduced neutral hadron BG, which, in order for the recoiled hadron to be energetic, require larger momentum transfer.
As a result, for the signal of interest to us, no significant recoil energy is expected and veto layers in front of the decay vessel often remain inactive, unless the scatterings process takes place right in front of them. We then always require in our analysis that p recoil < 1 GeV, while the typical recoil momentum is even smaller, p recoil ∼ O(100 MeV). In addition, we exclude all scattering processes happening in the last three hadronic interaction lengths, 3 λ had,int , of the material lying in the most immediate neighborhood of the veto layers. While the last condition follows our generally conservative approach, we also discuss in Appendix A how the sensitivity reach can be changed by ameliorating, or strengthening, this cut.
The geometry for each experiment is illustrated in Fig. 4. Below, we briefly describe basic details of the experiments that are relevant for our discussion. More detailed description of simplified detector designs employed in our study is given in Appendix D, and physics cuts applied to signal events in each experiment are discussed in Appendix E.
A. FASER a. Basic detector design The FASER experiment has been proposed [6] to search for LLPs that can be abundantly produced in the forward direction of the LHC and subsequently decay in a distant detector [6,8,51,56,57]. Following a preparation of the detailed detector design [7,9], the experiment was approved to take data during LHC Run 3. In addition, a larger version of the detector has been considered [8], dubbed FASER 2, that could collect data during the HL-LHC era. In the following, we will present the results for both these experiments, as well as for a small version of the FASER detector left for the entire HL-LHC era, that we further denote by FASER (HL-LHC).
The FASER detector consists of a cylindrically-shaped decay vessel with length ∆ and radius R, followed by a spectrometer and calorimeter. At the front of the detector, there is a veto layer with a primary role to reject muon-induced BG. On top of the main FASER detector described above, recently an additional detector component has been proposed [58], dubbed FASERν, to be placed in front of the main detector with a primary role to study interactions of high-energy neutrinos. This is an emulsion detector covering the transverse size of FASER with tungsten layers interleaved with emulsion films. The total length of tungsten in FASERν is 1 m. Both in front of FASERν, as well as after the emulsion detector, there are scintillator layers acting as veto for muon-induced BG.
As far as secondary production of LLPs is concerned, FASER can employ all the material in front of the detector that could serve as a target. In particular, these include the aforementioned concrete shielding and rock, but also FASERν subdetector. In the following, we will take into account secondary production processes happening in all these places. To this end, we will assume that both FASER (Run 3) and FASER (HL-LHC) are equipped with the FASERν subdetector, while FASER 2 does not contain it due to its larger size that needs to be fitted in the TI12 tunnel. Further details about detector design are shown in Figs. 4a and 4b for FASER and FASER 2, respectively.
b. LLP decay signature in the decay vessel The LLPs produced at the ATLAS IP or in secondary processes in the material in front of FASER, can decay inside the fiducial volume of the detector. The expected BG in searches for such decays can be suppressed to negligible levels by employing basic signal characteristics. First, a high-energy cut on visible energy in the detector, E visible > 100 GeV, can be applied to reject all soft BG particles with only a small impact on the expected number of signal events. In addition, pointing and timing information can be used to further associate two charged tracks detected from LLP decay with pp collisions happening at the ATLAS IP. As a result, FASER can operate in an essentially BG free environment [7]. In the following, we also assume 100% detection efficiency for easier comparison with other experiments.

B. MATHUSLA
a. Basic detector design Similarly to FASER, the MATHUSLA experiment has been proposed [10] to take advantage of possible abundant production of long-lived BSM particles in pp collisions at the LHC that could have escaped detection in existing searches [12]. It employs a large-scale detector to be placed above the CMS IP to take data during the HL-LHC era. Here, we follow the proposed MATHUSLA100 design [11] with additional excavation performed to position the detector partially underground, more close to the CMS IP, in order to maximize the physics reach [59].
The main part of MATHUSLA is an air-filled decay volume of the size of 20 m × 100 m × 100 m which is followed by a tracking system that spans over the entire area of the detector and is placed on top of the decay volume. In addition, scintillator layers are placed at the bottom and on the sides of the decay volume to veto charged particles entering the detector. The precise geometry of the decay volume that we use reads where x corresponds to an upward direction, while z is the direction along the LHC beam pipe and the origin of the coordinate system is placed in the CMS IP. When modeling secondary production of LLPs, we take into account the rock separating MATHUSLA from the CMS cavern. We illustrate this design in Fig. 4c. b. LLP decay signature in the decay vessel In our analysis, we take into account all the LLP decay events producing two charged SM tracks that happen inside the decay volume of MATHUSLA. We assume 100% efficiency of detection and the daughter track separation for  energies p daughter > 1 GeV. Position and timing information about the tracks is used to identify the vertex in the decay volume. This, along with the direction reconstruction of the coming LLP, allows one to greatly reduce various sources of BG including cosmic-rays, muons from pp collisions at the CMS IP and atmospheric-neutrinoinduced BG. Since the direction of the LLP is only mildly changed by the recoil in the secondary production processes, as we discuss in Appendix C, we will assume in the following that the search of our interest can be performed with zero BG.
C. SHiP a. Basic detector design The proposed SHiP detector [13] is a fixed-target experiment that will use the SPS beam of 400 GeV protons incident on a target material made out of titanium-zirconium doped molybdenum alloy and tungsten. A nominal value of N POT = 2 × 10 20 protons on target (POT) allows one to potentially produce a large number of LLPs that could subsequently decay in a distant decay vessel.
In between the target and decay vessel, there is a place for a hadronic stopper and an active muon shield with an essential role to reduce muon-induced BG to negli-gible levels. The active muon shield is followed by the Scattering and Neutrino Detector (SND) which has been designed to study the interactions of SM neutrinos and light DM particles. After the SND, a 50 m long decay vessel begins. Decays of LLPs into charged SM particles are detected by observing the resulting tracks in the Decay Spectrometer (DS).
Since the initial release of the technical proposal [14], SHiP detector design underwent revision primarily in order to reduce the cost and weight of the active muon shield while maintaining its assumed high-quality performance. In our analysis, we follow the recent technical update [60], while we note that further possible modifications to the design might require updating the results in the future. We simplify a complicated design of the planned SHiP detector. However, when doing so, we keep its key components that could play an essential role for secondary production of LLPs. A schematic drawing of the SHiP detector design is shown in Fig. 4d.
The most important part of the detector with regards to secondary production of LLPs is the SND and its surrounding magnet. Scattering processes occuring there can lead to LLPs produced only several meters in front of the decay vessel. The SND consists of an emulsion detector which is followed by the SND muon system. Since the latter can partially act as a front veto for the decay vessel, we exclude from our analysis all scattering events happening in the material lying in its close neighborhood, within 3 λ had,int . Notably, on top of absorption of soft particles with p < 1 GeV in the emulsion detector in front of the SND muon system, charged soft particles produced in the emulsion detector can also be swept away by the SND magnet and then never reach the SND muon system. This then prevents the events from being rejected as combinatorial BG.
The emulsion detector in the SND is also equipped with electronic tracking layers that can time stamp the events. These could detect even soft activity, corresponding to energy O(100 MeV), in case of secondary production taking place inside the emulsion detector, even though the SND muon system will not be activated. In the following, we assume that events will not be rejected based solely on this soft activity in the emulsion detector. On the other hand, even if such rejection is present, additional secondary production processes can take place in the surrounding magnets and in the muon shield that will not be vetoed. These events will then always contribute to the secondary production rate.
b. LLP decay signature in the decay vessel A detailed reconstruction of signal events in SHiP with two charged and energetic tracks from LLP decays within the decay volume, as well as BG discrimination, employs a number of observables including e.g. the momentum of detected tracks and their impact factor with respect to the target at the IP. In the following, we apply a simplified acceptance procedure that is primarily based on the momentum of the visible tracks coming from the LLP decay. In particular, we require each visible track to have p 1 GeV.
In fact, similarly to MATHUSLA, we find that the precise value of this low momentum threshold does not play an important role, at least as far as secondary production of LLPs is concerned. In this regime, even setting the lower limit for the visible energy at the level of E vis > 10 GeV leads to very similar results, as discussed in Appendix A.
When obtaining the sensitivity reach plots, we assume that SHiP can operate in zero BG environment and can detect signal events with 100% efficiency. See e.g. Ref. [60] for further discussion about BG in SHiP and Ref. [14] for more realistic studies of the efficiency.

IV. DETAILS OF MODELING
As discussed above and illustrated in Fig. 1, models with more than a single LLP can effectively lead to production of BSM species both in initial pp and pN interactions at the LHC or in the target material, as well as in scattering processes taking place more closely to the decay vessel. We refer to the former as primary production, while to the latter as secondary.
Once produced, LLPs can travel some distance and, eventually, decay in the detector leading to visible signal of two oppositely-charged tracks. The probability of this to happen in the decay vessel depends on the decay length of a boosted LLP, as well as on the geometrical acceptance of the detector.
Below, we briefly summarize our analysis, while further details, including the expressions for relevant branching ratios, decay widths and scattering cross sections are given in Appendices B and C.

A. Primary production of the LLPs
Light new physics particles can be produced in highenergy hadronic interactions in a number of processes. The dominant production channels depend on the BSM model of interest and on the mass of the LLP.
a. Dark photons Given our focus on MeV − GeV mass range dark photons, the dominant production processes relevant for our simulations are: Meson decays: Light dark photons can most efficiently be produced in rare decays of mesons, if kinematically available. We simulate meson distributions produced in pp collisions at the LHC and pN collisions with the molybdenum target at SHiP with the Monte-Carlo (MC) event generator EPOS-LHC [61], as implemented in the CRMC simulation package [62]. In our simulations, we take into account possible rare BSM decays of pions, η and η mesons, as well as vector mesons ρ and ω. We focus on the dominant decay channels π 0 , η, η → γA and ρ, ω → πA .
Proton bremsstrahlung of A : Dark photons can also be efficiently produced due to bremsstrahlung in coherent proton scatterings. This is especially relevant for A heavier than the threshold for production in rare pion and η meson decays.
We model the bremsstrahlung of A following the Fermi-Weizsacker-Williams (FWW) approximation and taking into account additional proton form factor to allow for off-shell mixing with vector mesons ρ and ω, see e.g. Refs. [6,33,63,64] for a recent discussion. The mixing with vector mesons leads to an increased production rate of dark photons in proton bremsstrahlung for m A ∼ 775 MeV.
In case of MATHUSLA, focus on high-p T -regime of the LHC makes it more challenging to directly apply the FWW formalism. Instead, here we approximate the relevant contribution to the dark photon spectrum by rescaling the spectrum of vector meson ρ by the appropriate mixing angle θ V as discussed in Appendix D of Ref. [47]. We note, however, that this contribution plays a subdominant role with respect to the one from meson decays in the region of parameter space of the inelastic DM model probed by MATHUSLA.
Hard processes: For even heavier dark photons, with masses m A 1.5 GeV, hard scattering contribution from Drell-Yan A production can start to play a dominant role. We take this into account, although it only concerns a small part of the parameter space of the inelastic DM model under study that is relevant for SHiP.
Other (subdominant) processes: Dark photons could also be produced in various other processes, e.g. in subsequent showers, but these have been found to be subdominant and, therefore, we neglect them in the following. As discussed below, secondary production becomes the most prominent in regions of the parameter space characterized by relatively short lifetime where large boost factor is required for the LLP to reach the detector before decaying. It is then sufficient for us to focus on the high-energy part of the LLP spectrum that is dominantly associated with the initial proton interactions in the target.
b. Single dark matter particle For the benchmark scenario described in Sec. II A, a flux of dark matter particles going towards the detector comes primarily from 3-body decays of light pseudoscalar mesons via off-shell dark photon, i.e. π 0 , η → γ A * → γ χχ [30]. The contribution from other processes like proton bremsstrahlung could become important for heavier χ, but it is subdominant in the mass-range of our interest.
c. Inelastic dark matter The benchmark model with inelastic DM discussed in Sec. II B is characterized by the dark photon mass exceeding the masses of two dark fermions, m A > m χ1 + m χ2 . In this case, the dominant 2-body decays of on-shell dark photons into χ 1 χ 2 pair become possible. The parent dark photon flux is governed by one of the production processes discussed above, depending on the A mass.
In addition, since χ 2 s are not stable, they decay into χ 1 and, typically, an electron-positron pair. These decays allow one to detect χ 2 s, when they happen inside the decay vessel. On the other hand, when χ 2 s decay before reaching the decay vessel, the resulting lighter fermions additionally contribute to the flux of χ 1 relevant for secondary production discussed below. We take such displaced χ 2 decays into account in our simulations.
d. Secluded dark Higgs boson As discussed in Sec. II C, secluded dark Higgs boson that couples to the SM predominantly via dark photon, can be efficiently co-produced in any of the A production mechanisms described above. In particular, for the mass range of our interest, the most important such production comes from 3-body meson decays π 0 , η → γ A * → γA S with intermediate dark photon, as well as from 2-body decays of vector mesons such as ρ → S A .

B. Secondary production of the LLPs
On top of the LLPs produced in the vicinity of the primary proton IP, in models of our interest it is also possible for them to appear in secondary production processes in material closer to the decay vessel. The relevant Feynman diagrams are shown in Fig. 3 in which an incoming particle LLP 1 , that corresponds to χ, χ 1 or S in the models described in Secs. II A to II C, respectively, produces an outgoing species LLP 2 with possibly much smaller lifetime, where LLP 2 = A or χ 2 in the considered scenarios.
These scatterings can occur both on electrons and nuclei with the latter giving the dominant contribution, especially in the regime of coherent scatterings with Z 2 enhancement. In addition, the coherent scatterings are associated with a small nuclear target recoil, p recoil ∼ O(100 MeV), that typically do not cause any veto activation, as well as only mildly affect the momentum of LLPs. While we perform MC simulations to account for the latter effect, we note that a very good approximation of the sensitivity reach can be obtained when working in the simplified, collinear regime with p LLP1 ≈ p LLP2 . A more detailed discussion of the relevant scattering cross sections is given in Appendix C.

C. Event rate
When LLPs produced in one of either primary or secondary production processes decay inside the decay vessel, this can lead to visible signal in the detector. The number of expected signal events depends on the relevant production rates, as well as on decay in volume probability that takes into account the acceptance factor A. The latter depends on the geometry of the detector, as well as on the efficiency to generate and detect visible charged tracks satisfying experimental cuts.
In this study, we perform full numerical MC simulations, which takes into account the interaction kinematics and experimental geometry, to obtain the sensitivity reach plots presented in Sec. V. Although in a strict sense not identical, it is illustrative to consider the probabilities P prim. for a LLP 2 produced at the IP and P sec. for an LLP 1 with subsequent interaction producing a LLP 2 to lead to the signal in the detector.
In case of primary production, the decay in volume probability of the LLP with momentum p can be written as where L is the distance to the beginning of the decay vessel and ∆ is the length of the vessel. The decay length of the LLP in ultrarelativistic regime reads d = cτ βγ cτ E/m where τ is the LLP lifetime, m corresponds to its mass, while the LLP energy is given by E. The acceptance factors A relevant for each of the experiments are included in numerical simulations.
For the LLPs coming from secondary production, the decay in volume probability needs to be convoluted with the scattering rate, which can differ for the various materials in which secondary production occurs. For a single incoming LLP 1 , the relevant probability to produce a LLP 2 which then decays inside the detector reads Here the integration limits x min and x max correspond to the distance to the decay vessel. They are dictated by the geometry of the detector and surrounding material, as well as by veto-related requirements discussed above. The interaction length is given by − where σ is the scattering cross section per nucleus relevant for secondary production mechanism, ρ is the density of material, m T is the mass of the target nucleus and in Eq. (13) we assumed that int.
(x max − x min ) which is always the case for scenarios of our interest. Contributions associated with various detector components are then summed up to obtain the total event rate.
As manifest in Eq. (13), secondary production processes are typically subdominant with respect to primary ones due to additional suppression for small scattering cross section. However, since secondary production can take place much closer to the detector, x min L, it allows one to probe short lifetime regime where contribution from primary production is already highly suppressed by the exponential factor in Eq. (12), exp (−L/d) 1.

V. RESULTS FOR SCATTERING OFF NUCLEI
In order to illustrate the interplay between the primary and secondary production mechanisms, we have studied the sensitivity reach for selected LLP models in the FASER, MATHUSLA and SHiP experiments. The respective results are shown in Fig. 5 for the model with A produced from dark bremsstrahlung, in Fig. 6 for the model with inelastic DM, and in Fig. 7 for the model with the secluded dark Higgs boson, that are described in Secs. II A to II C, respectively.

A. Current bounds from past and existing experiments
Before analyzing the reach of future experiments, we first discuss current bounds on the benchmark scenarios under study. These are associated with both searches for dark photons and for other LLPs present in the models. Below, we briefly summarize these constraints, while the most stringent of them are also shown in Figs. 5 to 7 as grey-shaded regions. A more comprehensive discussion can be found e.g. in Ref. [1] and references therein.
a. Dark photon with visible decays Current constraints on visibly decaying dark photons are relevant for models with dark bremsstrahlung of A and with secluded dark Higgs boson. The upper limit on the kinetic mixing parameter, 10 −3 , shown in Figs. 5 and 7, come from null searches for narrow resonances in e + e − → γ A → γ + − events and rare pion decays, π 0 → γ e + e − , reported by the BaBar [65] and NA48/2 [66] collaborations, respectively. On the other hand, bounds on lower values of have been derived from reinterpretation of old results of a number of fixed-target experiments, including E141 [67] and NuCal [64], as well as the limits obtained by current searches in e.g. the NA64 [68] detector. Depending on the value of m A , these bounds can effectively exclude part of the dark photon parameter space with 10 −7 −10 −4 . b. Dark bremsstrahlung The bounds on the model discussed in Sec. II A follow the above discussion relevant for visibly decaying dark photon. An additional constraint shown in Fig. 5 corresponds to the search for light DM particles scattering off electrons in LSND that follow Refs [29,31,69].
We note that in the benchmark scenario considered in our study, the region of the parameter space relevant to secondary production of dark photon in dark bremsstrahlung process corresponds to too low thermal DM relic density of χ. In this scenario, χ would be one of the DM components that constitutes only a fraction of the total DM abundance. This inefficiency of purely thermal production of χ could, in principle, be compensated in non-standard cosmological scenarios or by adding additional non-thermal component to its relic density in further extensions of the model. However, it is important to note that the model with Dirac fermion DM ef- ficiently annihilating via vector portal into SM particles is tightly constrained by the Planck CMB data [70,71]. We indicate this in Fig. 5 by showing the lower bound on as a function of m A that is relevant for scenario with the relic abundance of χ coming only from the standard freeze-out mechanism.
A notable exception to this obstacle could be a scenario in which total χ relic density is generated due to an initial asymmetry between χ andχ [72,73]. In this case, χ could even correspond to the entire DM relic density provided that the symmetric component is sufficiently suppressed [74]. On the other hand, it would then be a subject to additional constraints from DM direct detection (DD) searches, as discussed e.g. in Ref. [29]. These will partially exclude the region of the parameter space relevant for larger values of and small m A , which is shown as allowed in Fig. 5. In addition, some larger part of this region would then be covered by future DM DD searches and would provide a complementary method for discovery, focusing on the DM particle χ, with respect to the intensity frontier searches of A discussed in this study. We do not show corresponding bounds and future sensitivities in Fig. 5 as their presence depends on additional assumptions about cosmological scenario. c. Inelastic DM In the case of the inelastic DM model, and for sufficiently long-lived χ 2 , the most relevant bounds come from null searches for invisibly decaying A reported by the BaBar [75] and NA64 [76] collaborations. These correspond to the dark photon decaying promptly into a χ 1 χ 2 pair with χ 2 decaying outside the detector. However, for a large mass splitting between the dark fermions, ∆ χ ∼ 2, the lifetime of χ 2 is typically too small for the heavier fermion to decay outside the detector, as can be deduced from Eq. (8). In this case, the invisible bounds could only apply to a small region of the parameter space relevant for m χ1 ∼ O(10 MeV) and, generally, remain less important than other constraints discussed below. In addition, constraints on visibly decaying dark photons also do not apply to this scenario. This is because three-body decays of χ 2 into χ 1 and e + e − pair would not lead to a narrow resonance peak in the electron-positron spectrum seen in the detector, in opposite to e.g. A decays into e + e − pairs with no missing energy.
The dominant bounds are then associated with measurements of the anomalous magnetic moment of electron [77,78] and muon [79,80] that constrain large deviations of these quantities from the SM predictions, cf. Refs [81,82] for recent discussion, as well as from electroweak precision observables at LEP and LHC [83,84]. In addition, we also recast bounds on χ 1 particles scattering off electrons in the E137 experiment [32,85], see also e.g. Ref. [86] for recent relevant discussion, assuming E e − ,rec. 1 GeV threshold for the recoil energy. Notably, the bound from lighter fermion scattering in the LSND, which is discussed above for the model with a single DM particle, is suppressed in case of inelastic DM due to a large mass splitting between χ 1 and χ 2 , as well as a very small self-coupling of the lighter fermion.
On the other hand, heavier fermions produced at the IP could travel to the detector and leave a decay signature. This requires the χ 2 decay width to be suppressed enough, so that it does not decay too early.
The corresponding bounds from decays in the E137 [85] and LSND [87] detectors are shown in Fig. 6, following Ref. [38]. Importantly, they constrain the region of the parameter space relevant for values of the kinetic mixing parameter in the range 10 −6 10 −4 , while some regions of the parameter space with larger , but lying below the upper bound 10 −3 , is left unconstrained. Last, but not least, we note that in the model with inelastic DM the aforementioned stringent bounds from CMB can be avoided as annihilations of Majorana fermions, χ 1 χ 1 → SM, are p-wave suppressed. Similarly, there are no such bounds from coannihilation processes between χ 1 and χ 2 , as heavier χ 2 particles decay prior to a time of recombination.
When presenting the reach plots, we have assumed fixed mass ratios between all the three dark species, m χ1 : m χ2 : m A = 1 : 2.9 : 4 that were chosen for illustrative purposes. We note, however, that small changes in these benchmark ratios would have mild impact on the reach plots, as long as the dark photon decays predominantly into a χ 1 χ 2 pair. On the other hand, the relic density of χ 1 DM could change significantly, even by a few orders of magnitude, depending on how close are the assumed ratios to satisfy the resonance condition, m χ1 + m χ2 m A (see e.g. Ref. [88] for recent related discussion). The line with the correct value of Ω χ1 h 2 could then correspond to a wide range of , from values even lower than the region corresponding to the (g − 2) µ anomaly, up to large ones that are already excluded. For illustrative purposes, we present in Fig. 6 the line corresponding to the scenario with the resonance, following Ref. [38]. We take into account hadronic annihilation final states relevant for sufficiently large m χ1 , cf Ref. [40].
d. Dark photon and secluded dark Higgs boson The most important bounds on the model with the secluded dark Higgs boson come from the aforementioned searches for visibly decaying dark photons. Additional searches for S produced in rare decays of kaons, heavier mesons or the SM Higgs boson, that are based on the mixing be-tween S and H, do not constrain our parameter space of interest. This is because this mixing is highly suppressed for secluded dark Higgs boson.
The rate of light meson decays into S could, however, be enhanced thanks to their coupling to dark photons. This allows one to study corresponding constraints on such a scenario from subsequent delayed decays of S into e + e − pairs inside the detector, or from S scattering off electrons, S e − → A e − followed by a prompt A decay.
The constraints based on scattering of S can be derived similarly to the aforementioned bounds on fermionic DM from the LSND experiment. One important difference in this case is, though, the subsequent decay of A into e + e − that could additionally contribute to the EM signal in the detector. Given low energies relevant for LSND, typically E LLP O(100 MeV), these decays are prompt and, therefore, could be reconstructed as a singleelectron signal together with the recoiled electron. Assuming that this is the case and that the entire EM energy lies within the range relevant for LSND cuts [89], 18 MeV (E e − (recoil) + E e − + E e + ) E S 50 MeV or similarly in between 60 and 200 MeV [90], we have recasted the relevant bounds from Refs [29,43] and show them in Fig. 7 as the region excluded by the LSND.
In fact, the bounds that we present for LSND should be considered as an approximate indication of the already excluded region of the parameter space. A more detailed analysis would take into account that the opening angle between electrons and positrons in the final state could exceed 12 • , resulting in distinguishable lepton tracks which no longer mimic neutrino-electron elastic scatterings [91]. Such events with two or even three lepton tracks would, however, also be subject to constraints from the LSND, although deriving these bounds would require access to the experimental data and would then go beyond the scope of this study; see also Ref. [42] for relevant discussion.
On the other hand, the constraints based on delayed S decays into e + e − inside the detector are less severe. This is due to a very long lifetime of S, as indicated in Eq. (10). In addition, these bounds are also sensitive to a small mixing between S and H that could affect the exact value of τ S without having an impact on the main results presented below. For this reason, we only qualitatively discuss several such possible constraints, while we do not show them in Fig. 7.
In particular, as discussed in Refs [43,92], relevant constraints from electron beam-dump experiments are typically subdominant with respect to the ones that employ protons impinged on the target material. As an example, decays of S inside the LSND or MiniBoone detectors, that are interpreted as single electron events satisfying relevant cuts, can constrain an additional slice of the allowed parameter space, on top of standard dark photon searches. For negligible S-H mixing, this bound is, however, generally subdominant with respect to the aforementioned constraint from S scattering in the LSND. The same is true for possible bounds from S decays in the old high-energy proton beam-dump experiments, e.g. BEBC [93] or CHARM [94].
Importantly, due to the long lifetime of S, additional constraints can arise from Big Bang Nucleosynthesis (BBN) (see [95,96] for recent reviews). They are associated with possible co-production of S in scattering processes producing light dark photons in the early Universe. The dominant contributions come from annihilation processes qq, l + l − → A * → A S and could lead to a large number of long-lived dark Higgs bosons produced, depending on the value of and α D . On the other hand, a small but not negligible contribution from S-H mixing could allow one to always keep τ S below the value relevant for BBN, i.e. τ S 0.1 − 1 sec, while maintaining a secluded regime in terms of collider searches.

B. Sensitivity reach of future experiments
The sensitivity reach of all the experiments and models under study is shown in the left panels of Figs. 5 to 7. The reach corresponding the secondary production of LLPs is additionally marked with colorful shaded regions to distinguish it from the one relevant for the primary production. As can be seen, both types of production mechanisms can simultaneously cover some parts of the parameter space of the models, in which case the respective numbers of events add up to the total visible signal in the detector.
On top of this, some regions of the parameter space can only be covered by one type of the production process. This illustrates the complementarity between both mechanisms that allows one to probe additional scenarios that would otherwise seem to lie beyond the reach of the given experiment. It is important to note that a distinction that we focus on, corresponds to only the production mechanisms, while the signature in the detector is the same in both cases. In practice, this means that neglecting the impact of the secondary production could affect the interpretation of the results of the experiments.
Notably, the secondary production is most relevant for regions of parameter space with large couplings between the LLPs and the SM, which brightens the prospects for potential simultaneous discovery of such LLPs in many experiments in the near future. In order to properly interpret such a co-discovery, it would be essential to model secondary production of LLPs.
Focusing on the secondary production, we show in the right panels of Figs. 5 to 7 number of events contours for the experiments under study. As can be seen, depending on the model, we can typically expect up to O(10 3 ) events of this type, but this number can grow to even > 10 6 in certain scenarios, especially for the SHiP experiment. Below, we comment on the relevance of the secondary production mechanism for each of the experiments under study.

FASER (Run 3, HL-LHC) and FASER 2
The secondary production mechanism can improve the sensitivity of the FASER experiment both during LHC Run 3, as well as for the HL-LHC era. For all three considered benchmark models, the secondary production extends the FASER reach towards smaller LLP lifetimes, or larger values of . In the case of the inelastic DM benchmark model, even the FASER detector operating during LHC run 3 can probe a large region of currently unconstrained parameter space.
Note that the reach of FASER extends further to larger coupling in comparison to the reach of FASER 2, even though the latter detector has larger size. This is due to the additional dense material in the tungsten-based neutrino detector, FASERν, placed in front of the FASER decay volume. If such a block of tungsten, or a similarly dense material, would be placed in front of FASER 2, it would also have a positive impact on the physics reach provided that muon-induced BG can be successfully rejected by the front veto. This might be a useful observation for the detector design, especially in case some hints of new physics appear in the near future that correspond to the models with smaller LLP lifetimes than the FASER experiment can typically cover.
Interestingly, the regions of parameter space that are accessible at FASER through secondary production mechanisms can be related to outstanding problems in particle physics and anomalies: i) In case of the inelastic DM model, FASER probes the currently unexcluded region of the parameter space that yields the correct value of the DM relic density of χ 1 , as discussed in [38]. ii) In addition, the secondary production opens up the possibility to probe an unconstrained region of the parameter space of the inelastic DM model corresponding to the observed discrepancy between the measurements and SM predictions of the muon anomalous magnetic moment [79,80]. iii) Furthermore, the probed region of the dark photon parameter space can correspond to the existing anomaly observed in decays of the excited state of Beryllium-8 with a mass of m A = 17 MeV [97], as discussed e.g. in Ref. [1]. However, realistic BSM models that could accommodate for this anomaly [98,99] require going beyond the simplest dark photon case and modifying its couplings to nucleons. This would also affect the reach of both FASER detectors.
Notably, these particularly motivated regions of parameter space can already be probed by the FASER detector operating during the LHC Run 3, with additional positive role of FASERν, while this would not be possible if only the primary production at the IP were taken into account.

MATHUSLA
The primary target of the MATHUSLA detector, which is to be placed off the LHC beam collision axis, is new physics particles produced in rare decays of heavier SM species e.g. B mesons or the SM Higgs boson, or in hard pp scatterings that are the most relevant for larger LLP masses. It is then not a surprise that MATH-USLA has no reach in the currently allowed region of the parameter space of the models discussed above with below-GeV-mass A decaying in the detector.
Interestingly, however, in the model with inelastic DM illustrated in Fig. 6, MATHUSLA can have nonnegligible reach both for the primary and secondary production of LLPs. This is due to a larger angular spread in a DM flux produced with one additional decay of relatively heavier dark photon, A → χ 1 χ 2 , when compared to the A flux from e.g. pion decays. On top of this, the primary production benefits from a relatively large lifetime of χ 2 .
As a result, two distinct regions of the parameter space can be seen corresponding to both types of production mechanism with the secondary one covering the region with the correct value of DM relic density and, partially, the one corresponding to the aforementioned (g − 2) µ anomaly.

SHiP
The secondary production in case of the SHiP experiment has also a positive impact on the sensitivity to BSM searches. It can lead to not only an improved discovery prospects, but also to potentially a large number of such events. This is particularly relevant for the model with inelastic DM, as shown in the right panel of Fig. 6. Taking into account the secondary production could then also become crucial for the model parameter reconstruction, in case of discovery.
The sensitivity reach of SHiP for both models focusing on the dark photon decaying in the detector can also be improved, cf Figs. 5 and 7. Importantly, when obtaining the results, we have adapted a quite conservative approach to exclude all the secondary processes happening in the SHiP detector more close than about 2.5 m away from the decay vessel or so, as detailed in Sec. III and Appendices D and E. On the other hand, if these cuts can be weakened, the sensitivity of SHiP with respect to the secondary production will be further enhanced.
Interestingly, on top of the interplay between the secondary and primary production of LLPs, in case of the SHiP experiment, additional signatures of the models with more than a single LLP can arise by analyzing LLP scatterings off electrons. We discuss them in more detail in the section below.

VI. SCATTERING OFF ELECTRONS IN SHIP SND
We have so far been focusing on two possible production mechanisms of LLPs, but only on one standard LLP decay signature in the detector with two high-energy charged tracks emerging from the vertex in the decay vessel. On the other hand, models with more than a single LLP also offer additional types of experimental signatures that can be studied contingent on a specific detector design.
In particular, the presence of electronic tracker layers in the emulsion detector in the SHiP's SND can be markedly advantageous for BSM searches. Depending on their time resolution, they could time stamp the scattering events in the SND by detecting the corresponding recoil products. Such events could then be analyzed either separately or in conjunction with a subsequent LLP decay. Interactions of LLPs with the light electron target typically generate large recoil energy of e − , as discussed in Appendix C. This leads to a detectable EM cascade in the SND with no hadronic counterpart, for which the expected BG is greatly reduced.
While the actual capabilities of the SHiP experiment with respect to such signatures will depend on the final detector design, below we briefly discuss two such possible search strategies employing scatterings off electrons in the SND with possible secondary production of LLPs. They are also schematically shown in Fig. 8.
(Only) Electron scattering signature: The search for light dark sector particles scattering off electrons is one of the main aims of the SND [14]. The expected number of relevant BG events has been reported to be about 800 for 2 × 10 20 POT [60]. These events are mostly associated with quasi-elastic electron neutrino and anti-neutrino charged current (CC) interactions with nuclei. While it might be possible to further reduce this background, for example by using more sophisticated tracking algorithms to identify soft protons in ν e n → e − p interactions, there remain truly irreducible backgrounds related to elastic scattering off electrons νe → νe, or quasi-elastic scattering of anti-neutrinos,ν e p → e + n, with neutrons in the final state that easily avoid detection.
For our models of interest, a similar pure scattering signature is possible, especially if subsequent decay of the LLP from secondary production is delayed and happens outside the detector. In a simplified geometry employed in our study, this corresponds to decay of the LLP at a distance at least about 71 m away from the end of the SND, where we have assumed that the SHiP Decay Spectrometer has the total length of about 15 m. When obtaining the result, we follow Refs. [14,36] and apply cuts on the recoil energy, 1 GeV < E e < 20 GeV, and on the recoil electron angle, 10 mrad < θ e < 20 mrad.
"Double" signature inside the SND: Yet another spectacular signature of models with more than a single LLP can be associated with a simultaneous generation of two resolvable collinear EM showers with no hadronic recoil counterpart, that both can be seen inside the emulsion detector. These can be associated with the secondary production process followed by a prompt LLP decay. Importantly, neutrino-induced combinatorical BG to such process can be greatly suppressed.
While an optimized description of this search will require a separate study, here we discuss one such simple strategy in which both EM showers should satisfy the aforementioned cuts relevant for the pure scattering signature. In addition, we require both showers to be initiated at positions in the emulsion detector that are not too close to each other. As a benchmark gap between these positions, we choose a distance of 10 cm. This typically allows for each of the showers to be detected by different tracking layers, as these are placed every about 8 cm. Notably, a distance of 10 cm corresponds to more than 15 radiation lengths in lead. This significantly reduces any possible overlap between the showers in the emulsion films.
On the other hand, we note that requiring such a large gap between the two showers might be a too conservative approach. This is due to generally excellent spatial resolution of the emulsion detectors that could be combined with the signal detected by the tracker layers. We then also illustrate below how weakening of this cut might allow one to study even broader class of BSM scenarios.
While the signature based solely on scatterings in the SND for models with a single dark matter particle has already been discussed in the literature, see e.g. Refs [31][32][33][34][35][36], here we present exemplary results of such searches for the model with inelastic DM, as well as for the one with a dark photon and secluded dark Higgs boson. The relevant results corresponding to benchmark scenarios introduced in Secs. II B and II C are shown in Figs. 9a and 9b, respectively. In these plots, a region above dotdashed green lines corresponds to at least three electron scattering events happening in the SND. Hence, these lines indicate an absolute lower limit on as a function of the relevant mass that could be probed by any signature based on scatterings off electrons in the SND in the absence of BG. This can be compared with the reach from the standard decay in volume signature discussed in Sec. V. The corresponding results for both primary and secondary production of LLPs are shown as solid red and solid blue lines in the figure, respectively.
The results corresponding to the aforementioned pure scattering search are shown for the inelastic DM model in Fig. 9a as dotted purple lines with the fixed number of expected electron-scattering signal events, N = 3, 30, 300, 3000. Notably, the number of events in this case can reach up to O(1000) in the allowed region of the parameter space. This can be compared with order 100 events required to exceed the expected level of BG by 5σ where the error has been roughly estimated as √ N BG for the aforementioned number of BG events. These events are associated with χ 1 upscattering to χ 2 taking place in the SND, while χ 2 surviving for long enough so that  Both plots have been obtained for benchmark scenarios defined in Sec. II. The green dash-dotted lines correspond to N ev = 3 electron scattering events in the SND, while solid red and blue lines to the reach from the standard decay in volume signature and LLPs coming from primary or secondary production, respectively. In the left panel, the purple dotted (dashed golden) lines were obtained for the pure scattering (double) signature discussed in the text and correspond to different number of events as indicated in the plot. In the right panel, the dashed (dash-dotted, solid) golden lines represent the reach of the double signature with N ev = 3 events.
it will decay outside the decay vessel and DS. For sufficiently long-lived heavier fermion, an additional contribution to the total number of signal events would come from χ 2 downscattering to χ 1 in the SND.
Complementary search can be performed by employing the double signature of both upscattering and decay happening inside the SND with at least 10 cm gap in between. In this case, the maximum distance between the secondary production point and decay is the length of the SND, ∆ SND ∼ 1 m. This allows one to probe a regime of a smaller lifetime. As a result, the corresponding dashed golden lines with different number of events N ev = 3, 30, . . . cover a distinct region in the parameter space of the model, with larger values of , than in the pure scattering signature case discussed above.
As seen in Fig. 9a, a combination of various search strategies employing the SND and the decay spectrometer, either separately or in a joint signature, could further shed light on the nature and lifetime of the LLP.
A somewhat different scenario is illustrated in Fig. 9b for the model with secluded dark Higgs boson. Here, the LLP produced in secondary production processes in the SND is a very short-lived dark photon that promptly decays back to an electron-positron pair. As a result, no events with A being able to travel outside the detector are expected. Actually, even the 10 cm gap in between the two showers in the SND might be too large as dark photons will typically decay very fast. We then show several lines corresponding to 3 events with the double signature for a selection of minimal distances between the scattering and decay vertices ranging from 1 cm to 10 cm. As expected, the smaller is the allowed distance between the beginning of the two showers, the better is the reach, at least assuming zero BG. On the other hand, for too small such a distance, both showers could no longer be effectively resolved. In this case, we would again enter the regime of effectively pure scattering signature with a much larger expected BG. This would significantly affect the sensitivity reach. We leave for future studies a determination of the effective gap for which the showers can be distinguished based on an interplay between the signals measured by the tracker layers and in the emulsion films.
Last, but not least, we note that, if a sufficient number of events is observed with simultaneous EM showers that both can be resolved, more information useful for the model parameter reconstruction could be obtained from analyzing the distribution of the gaps between the vertices in such events.

VII. CONCLUSIONS
One of the most promising ways to find new light longlived particles is to search for their highly displaced decays in distant detectors separated from the primary interaction point. However, these potentially very clean searches are limited by the lifetime of LLPs that need to travel the entire distance to the detector. On the other hand, in models with more than a single LLP, often additional production mechanisms naturally appear due to efficient couplings between the BSM species. These can lead to the secondary production of LLPs in the close proximity of the decay vessel. As a result, smaller than usual LLP lifetimes can be probed, while the search still benefits from the shielding from the SM background, similarly to the case of LLPs produced at distant IP.
In this study, we have illustrated such a possibility for several theoretically motivated models extending the SM by additional U (1) d group and the corresponding gauge boson, namely the dark photon A . Particular models of our interest include: dark photon being a mediator between the SM and DM, with either a single Dirac fermion DM particle or a pair of Majorana dark fermions, as well as A gaining a non-zero mass via the dark Higgs mechanism and the vev of the secluded dark Higgs boson.
While these models serve as an example, the secondary production can play an important role in many types of BSM scenarios with light new particles. In particular, in a more realistic setup, A could both play a role of DM mediator and benefit from the presence of an additional scalar particle in the model. This would lead to combined signatures from different scenarios discussed in our study, cf Refs [38,43] for recent studies in this di-rection. In addition, more than a single LLP can arise e.g. in the context of the Twin Higgs scenario [5], supersymmetry (see e.g. Ref. [100]) or models with a linear dilaton action with possibly light Kaluza-Klein (KK) gravitons and scalars [101]. Such models have also been proposed to accommodate for past and present experimental anomalies, cf. e.g. Ref. [102] for a recent such study connected to the MiniBooNE excess and LSND anomaly. Remarkably, even in the presence of just a single LLP, the model can still be subject to effective secondary production e.g. if the BSM species couple to the SM neutrinos, cf. Ref. [103] for such study for heavy neutral leptons and the IceCube detector.
Importantly, the secondary production is the most important for the regions of the parameter space of the models that are often the most appealing due to prospects for a simultaneous discovery in many different experiments. By taking into account this production mechanism, intensity frontier searches can offer an independent way of verifying possible future hints of new physics that would otherwise seem to be beyond their reach.
This can also be true for the existing discrepancies between the theoretical predictions and measurements. In particular, we show how the secondary production of LLPs in front of the FASER or MATHUSLA detectors can allow them to probe the (g − 2) µ anomaly in the unconstrained region of the parameter space of the inelastic DM model. Notably, FASER capabilities with respect to this scenario are enhanced by a dedicated neutrino subdetector, FASERν, placed in front of the decay vessel.
Similar such Scattering and Neutrino Detector (SND) can also substantially add to discovery prospects of SHiP. On top of already carefully studied signature based on LLP scatterings in the subdetector, we also propose to search for simultaneous production and decay of the LLP in the SND with two separate and collinear electromagnetic showers that are coincident in time. In case of discovery, employing a combined information about all the signal events detected both in the SND and the decay vessel, could shed more light on the nature of the LLP.
Last but not least, while we analyze the impact of the secondary production for only selected intensity frontier experiments, namely FASER, MATHUSLA and SHiP, the same production mechanism can be relevant for other similar experiments e.g. CODEX-b [53,54], NA62 [55], SeaQuest [44] and further discussed i.a. in Ref. [1]. We encourage the experimental collaborations to take this possibility into account in detailed modeling, or even at the level of detector design, as it can crucially extend the sensitivity reach to new physics to promising scenarios with connection to current, or possible future, hints of new physics.

ACKNOWLEDGMENTS
We thank Kevin Kelly for providing the MG model file used in [29] and Luc Darmé for useful discussions.
ST would like to express a special thanks to the Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence PRISMA+ (Project ID 39083149) and the GGI Institute for Theoretical Physics for their hospitality and support. KJ and LR are supported in part by the National Science Centre (NCN) research grant No. 2015/18/A/ST2/00748. LR is also supported by the project "AstroCeNT: Particle Astrophysics Science and TechnologyCentre" carried out within the Interna-tional Research Agendas programme of the Foundation for Polish Science financed by the European Union under the European Regional Development Fund. FK is supported by U.S. Department of Energy Grant DE-AC02-76SF00515. F.K. performed part of this work at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. ST is supported by the Lancaster-Manchester-Sheffield Consortium for Fundamental Physics under STFC grant ST/P000800/1. The results presented in Sec. V correspond to a particular set of cuts that we impose when estimating the number of signal events, as detailed in Sec. III and Appendices D and E. In particular, for the energy thresholds we follow standard values relevant for each experiment, while as far as secondary production of LLPs is concerned in front of veto layers, we require the recoil momentum to be low, p recoil < 1 GeV, and reject all the events happening within the last three hadronic interaction lengths (3 λ had,int ) of the layers. This corresponds to about 1 m of standard rock, 0.5 m of lead and 30 cm of tungsten [104]. Here, we discuss how the result can be affected when one sharpens these cuts. For illustrative purposes, we focus on the model with inelastic DM discussed in Sec. II B.
It is instructive to first analyze the example energy distribution of the events seen in the detectors under study when they are associated with the secondary production of χ 2 in front of the decay vessel. This is shown in the top panels of Fig. 10 for the FASER 2, MATHUSLA and SHiP experiments, and for the benchmark scenarios with m A = 50 or 250 MeV and m χ1 : m χ2 : m A = 1 : 2.9 : 4, as well as for = 10 −3 or 10 −4 . As expected, the energy of the events grows with decreasing lifetime of χ 2 , i.e., with growing and LLP masses. This is due to larger boost factors required for χ 2 to reach the decay vessel before decaying. For this reason, one expects that the region of the parameter space of the model that is the most sensitive to increasing the energy threshold will correspond to lower values of and lighter LLP particles. Notably, these are typically scenarios for which the primary production of LLPs dominate over the secondary one.
In the bottom panels of Fig. 10, we show the distribution of the distance between the position of χ 1 upscattering to χ 2 , and the decay vessel. As can be seen, for smaller lifetime of χ 2 , the impact of the last part of the material right in front of the decay vessel becomes more important. This is illustrated as more steeply growing histograms for increasing and the LLP masses. Hence, one can expect that excluding from the analysis larger parts of the material in front of the detector, would most straightforwardly affect this region of the parameter space of the model.
In the opposite limit, for small and light LLPs, the distribution of the upscattering events becomes flat in the distance to the detector. In this regime, the lifetime of χ 2 is large enough for it to travel long distances. As a result, in this regime, the primary production of χ 2 dominates over the secondary one. A somewhat unexpected behavior of the solid and dashed blue lines in the figure relevant for the SHiP experiment, that corresponds to such long-lifetime regime and describe diminishing number of events from the region in the most immediate vicinity of the decay vessel, is connected to the geometry of the detector and the presence of the SND with smaller average density than surrounding magnets and the muon shielding.
We illustrate all these effects in Fig. 11 where we present the sensitivity reach plots for the considered experiments. In this figure, the solid black lines correspond to the default results discussed in Sec. V, while other lines were obtained by varying the length of material in front of the decay vessel acting as an effective veto (blue lines), and varying the energy thresholds (red lines) .
In our default analysis setup we require the secondary production to occur at least 3 λ had,int before the front vetos. However, this cut is rather conservative and could be relaxed. We require the target nucleus to have recoil momenta below p recoil < 1 GeV, resulting in velocities v/c < 0.1. At these velocities the nucleus will experience an energy loss dE/dx 10 MeV cm 2 /g. Therefore the nucleus will stop after traveling a distance ∆x = p recoil /(dE/dx) 100 g/cm 2 which is smaller than λ int,had for all considered targets. Changing the effective length of the front veto has a large effect on the sensitivity reach for LLPs with short lifetimes. This is illustrated by the blue lines in Fig. 11, where we exclude the last λ had,int (5 λ had,int ) of material, resulting in an improved (reduced) reach at large couplings .
The impact of varying the energy thresholds affects mostly the reach in lower values of the kinetic mixing parameter. Here, even soft LLPs can reach the detector before decaying due to generally larger lifetime τ χ2 . Therefore, increasing the energy threshold can have an impact on the number of signal events. However, for FASER 2 and SHiP, this mostly affects the region of the parameter space where sensitivity is, either way, driven by the contribution due to the primary production. This is not the case for the MATHUSLA detector where both primary and secondary production regimes are clearly separated. On the other hand, the MATHUSLA detector does not have the energy measurement,  Dark Photon Production: If the dark photon is light, it is mainly produced via the decay of a pseudoscalar meson P = π 0 , η, η into a into a photon and a dark photon. The corresponding branching fraction is given by Here we have used the Källén function λ(a, b, c) = (a − b − c) 2 − 4bc, where m P is the pseudoscalar meson mass and B P →γγ is the branching fraction of the di-photon decay channel of P . Additionally, the dark photon can be produced in decays of a vector meson V = ρ, ω into a pseudoscalar meson P and a dark photon (such as ω → π 0 A ), or via the decay of a pseudoscalar meson P into a vector meson V and a dark photon (such as η → ρ 0 A ). The corresponding branching fractions are given by [44]

(B2)
Fermionic Dark Matter: For heavy dark fermions, m 1 + m 2 > m A , the differential branching fraction of a pseduscalar meson P decaying into γχ 1 χ 2 via intermediate off-shell dark photon A * is given by [28,30]  (c) SHiP Here, s is the squared four-momentum of A * , while the angle θ corresponds to the momentum of χ 1 in the rest frame of A * that is measured with respect to the boost direction of the off-shell dark photon. The case with a single dark fermion corresponds to S = 1/2 and m 1 = m 2 ≡ m χ , while for distinct dark fermions S = 1.
For light dark fermions with m 1 + m 2 < m A , the dominant production mode of χ 1 χ 2 pairs is due to prompt decays of on-shell dark photons, A → χ 1 χ 2 , with the relevant branching ratio of B(A → χ 1 χ 2 ) 1 dictated by the hierarchy between the coupling constants, α D α 2 .
Secluded dark Higgs boson: The branching fraction of a pseudoscalar meson decay P → γA * → γ S A reads [43] where s and θ are defined similarly to Eq. (B3) but with χ replaced by S. Additionally, the branching fraction relevant for vector meson, e.g. ρ, into a dark Higgs boson and a dark photon, ρ → SA , is equal to [31,43] (B5)

Decays of LLPs
Dark Photon Decay: The decay width of an on-shell dark photon into SM particles is where m e is the electron mass and B e = B(A → e + e − ) is the branching fraction of a decay into an electron-positron pair. For m A below the di-muon threshold, B e = 1 which is typically the case in our analysis. For heavier dark photons, decays into µ + µ − , as well as hadronic final states start to play important role [105].
Decay of χ 2 in inelastic Dark Matter Models: The relevant differential branching fraction for decay into an electron-positron pair, χ 2 → χ 1 e + e − , can be found e.g. in Ref. [37]. In case of heavier dark fermions and a larger mass splitting between them, decays into other SM particles become kinematically allowed. We take this into account by including the branching fraction of an off-shell dark photon decaying into an electron-positron pair, B e (m A * = m ee ), evaluated at the invariant mass of the electron pair m ee . The decay width is then given by where Γ A is given in Eq. (B6) and (B9) In Fig. 12 we show sample distributions of the electron pair energy as well as the angle of electron pair momentum the with respect to the direction of χ 2 . As can be seen in the left panel, in the boosted regime, where E χ2 m χ2 , the distribution of the energy fraction of χ 2 that goes into the visible final state, x = E e + e − /E χ2 , only mildly depends on E χ2 . In particular, for a larger mass splitting between the dark fermions characteristic for scenarios of our interest, the visible energy often corresponds to more than a half of the energy of χ 2 , with a kinematical upper limit E e + e − 1 − (m χ1 /m χ2 ) 2 E χ2 . The relevant value for the assumed mass ratio between both dark fermions is given by E e + e − 0.9 E χ2 . The upper limit on E e + e − becomes more stringent in case of lower energies, E χ2 ∼ m χ2 , due to a non-zero mass of χ 1 in the final state.
In the right panel, we show the angle between electron pair momentum relative to the direction of decaying χ 2 . Notably, the visible charged tracks are typically collimated along the direction of χ 2 . This is especially true for energies E χ2 100 GeV i.e. above the energy threshold of FASER. On the other hand, for energies E χ2 ∼ 10 GeV, that are more relevant for the secondary production in the SHiP and MATHUSLA experiments, the deflection angle of order few tens of mrads leads to the impact parameter with respect to the IP of order meters, given the distance between the IP and the decay vessel exceeding 50 m. This might partially limit the actual reach of both experiments in this model once a more detailed experimental analysis is performed that considers a pointing requirement. This effect is even more important for some scenarios with lower energy O(GeV), corresponding to the regions of the parameter space dominated by the primary production of LLPs.
Appendix C: Scattering cross sections for secondary production of LLPs In this section, we will give the expressions for the scattering cross sections relevant for secondary production of new light particles, LLP 1 + T → LLP 2 + T , where the target can be electron or nuclei, T = e, N . For the models discussed in Secs. II A to II C, LLP 1 = χ, χ 1 or S, while LLP 2 = A , χ 2 or A , respectively. Notably, all the considered models are characterized by somewhat different kinematics when both the secondary production and subsequent decay of LLP 2 in the detector are taken into account: • 2 → 3 scattering followed by a 2-body decay (relevant for scenario with dark bremsstrahlung, cf Sec. II A), • 2 → 2 scattering followed by a 3-body decay (inelastic DM, cf Sec. II B), • 2 → 2 scattering followed by a 2-body decay (A with secluded dark Higgs boson, cf Sec. II C).
In the following, we present the expressions for the scattering cross sections relevant to our analysis.

2 → 2 scatterings
In case of models with inelastic DM and secluded dark Higgs boson, secondary production proceeds via 2 → 2 scattering processes illustrated in the center and right panels of Fig. 3, respectively. The relevant differential cross section in the lab frame reads where E T and m T is the energy and mass of the recoiling target, while m 1 ≡ m LLP1 is the mass of the initial state scattered LLP. The integration limits are: where E 1 and E 2 are energies of the initial and final state LLPs, respectively, and m 2 ≡ m LLP2 .

a. Scatterings off Electrons
Let us first consider the case on of scattering off electrons.
Inelastic DM: Following Ref. [41], the squared matrix element for upscattering χ 1 e → χ 2 e is given by where the amplitude is defined as with δm χ = m χ2 − m χ1 .
Secluded dark Higgs boson: Following Ref. [92], the squared matrix element for the process Se → A e is where E 1 = E S and E 2 = E A .

b. Scatterings on Nuclei
The coupling between the dark photon and protons is dependent on nuclear form factors. This introduces a nontrivial dependence on the momentum transfer squared, q 2 = −Q 2 < 0, as dictated by the nucleon electromagnetic current, J µ =ū(p 4 ) [F 1 γ µ − (σ µν q ν /2 m p )] u(p 2 ), where F 1 (q 2 ) and F 2 (q 2 ) are the Dirac and Pauli form factors, respectively. These are usually expressed through the Sachs form factors, G E = F 1 − τ F 2 and G M = F 1 + F 2 , where we defined τ = Q 2 /(4 m 2 p ) > 0. It is customary to express G E employing the dipole approximation, which is particularly useful in the regime of low momentum exchange relevant for our analysis, which is characterized by Q 2 1 GeV, where µ p = 2.79 is the proton magnetic moment and in the last step we have followed the usual approximation that a simple scaling G M µ p G E holds for small, but even non-zero, values of Q 2 . For our purposes, the most convenient parametrization consists of form factors which take the form As follows from Eq. (C7), for a small momentum transfer, τ 1, one obtains G 2 G 1 and the term in the cross section proportional to G 1 can typically be neglected. The term proportional to G 2 , instead, plays the dominant role.
In addition, for sufficiently small momentum transfer, the internal structure of the nuclei is not probed and the dominant contribution to the cross section comes from coherent scatterings off entire nuclei. In the following, we implement approximate nuclear form factors that lead to a relatively smooth transition between the incoherent and coherent regimes [106][107][108] where G 2,el and G 2,inel are the form factors corresponding to the coherent and incoherent regimes, respectively. They are given by where the relevant atomic form factors have parameters a = 111 Z −1/3 /m e and a = 773 Z −2/3 /m e for coherent and incoherent form factors, respectively, while the coherent nucleus form factor is characterized by d = 0.164 GeV 2 A −2/3 with Z (A) being the atomic (mass) number of the nucleus.
We note that our G 2 form factor differs slightly from the relevant expression in Ref. [108] due to a spurious square factor in the nuclear part of the form factor appearing in that reference. We also note that, given our focus on Q 2 1 GeV, we do not have to consider the regime of deep inelastic scattering (DIS).
Inelastic DM: For the upscattering process χ 1 N → χ 2 N , the squared matrix element is given by where M 0 is given by Eq. (C4) and M 1 reads Secluded dark Higgs Boson: Finally, for the process SN → A N the squared matrix element is given by As discussed above, the scattering cross section is dominated by the contributions proportional to G 2 , effectively replaced by G 2,tot given in Eq. (C8), while the terms proportional to G 1 can be neglected in case of low momentum exchange. c. Energy and angular distribution in the secondary production Importantly, as illustrated in Fig. 13 for the case of inelastic DM, the energies and direction of the LLP is not much affected by the recoil, i.e. p χ2 ≈ p χ1 . In particular, the outgoing LLP tends to inherit most of the energy of the incoming high-energy particle, as shown in the left panel.
On the other hand, even a small deflection angle θ χ2 can play a role in the physics analysis. In particular, for energies E χ2 ∼ 100 GeV relevant for FASER, a typical value of θ χ2 ∼ 3 mrad deduced from the right panel of Fig. 13 leads to a displacement that exceeds the 10 cm radius of the detector, if the upscattering happens more than about 30 m away from the decay vessel. In this case, it is particularly important to properly model the deflection in the secondary production during the simulation. Instead, for larger energies or LLPs with a smaller lifetime, that are produced in the last few meters in front of the decay vessel, as well as for large-size FASER 2 detector, small values of the deflection angle play less important role in modeling.
Similarly for SHiP, assuming that the interaction takes place close to the decay vessel at a distance about 50 m away from the IP, and that E χ2 ∼ 10 GeV, a typical deflection angle, θ χ2 30 mrad, leads to an impact parameter with respect to the target of order 150 cm. This value is smaller than 250 cm relevant for the loose selection cuts, as mentioned in Appendix E, and it would be further suppressed for larger energies of χ 1 . On the other hand, a small but non-zero value of the deflection angle might render it impossible to employ the tight selection cuts typical for e.g. dark photon searches. This would be further complicated by the presence of an additional missing energy in the final state, as discussed in Appendix B 2.
In addition, if the secondary production takes place away from the decay vessel, than the displacement O(m) could result in LLP missing the vessel. This, however, mainly concerns scenarios with relatively long-lived LLPs for which the dominant production mechanism is, either way, the primary production.

2 → 3 scattering
In case of the dark photon mediated single dark matter model, the secondary production of the dark photon proceeds through the dark bremsstrahlung processes illustrated in the Feynman diagram in left panel of Fig. 3. In this work, the elastic 2 → 3 scattering of the dark matter off the nucleus is modeled numerically using MadGraph 5 [109]. We have created a simplified model containing the DM particle χ, a nucleus N and the dark photon A , where the effective coupling of the dark photon to the nucleus takes into account the elastic form factor G 2 (t) defined in Eq. (C8). We simulate the interaction separately for all considered materials and for various energies of the incoming dark matter particle. The simulated events are then convoluted with the incoming DM flux, where we sample over the interaction location and weight by the decay in volume probability.
Here we summarize main features of the simplified benchmark detector designs for the experiments considered in our study. These are also illustrated in Fig. 4. Bold case is used for these pieces of material, for which we model secondary production of LLPs in scatterings of other light particles. Any possible secondary production processes happening in other parts of the detectors are not taken into account.
More detailed experimental setups have been discussed for FASER/FASER 2 detectors in Refs. [8,9,58], for MATHUSLA in Ref. [11,12], while in case of SHiP in Refs. [14,60]. In the following, we employ simplified geometries that encompass the most important aspects of detector designs under study and, therefore, remain sufficient for the purpose of our phenomenological analysis.
These benchmark geometries might be further updated in the future, especially for experiments to be performed at a time of HL-LHC era, in which case sensitivity plots presented in this study would have to be updated. We note, however, that the results presented in Sec. V are typically not sensitive to possible (mild) changes in detector designs. The most important are potential modifications of detector design in the most immediate neighborhood of the decay vessel. Depending on the nature of these changes, they could either improve or worsen reach in the regions of parameter space of BSM models corresponding to LLPs with the shortest lifetime. Even such modifications, however, would often not have siginificant impact on the reach, as we briefly discuss in Appendix A.
FASER (both Run3 and HL-LHC): (beginning with most distant places away from the detector) • We do not model in details all the LHC infrastructure that can be found in between ATLAS IP and concrete shielding close to the place when the arc section of the LHC tunnel begins. While such detailed modeling would go beyond this study, we also note that the most relevant sections of the LHC tunnel for LLP production considered here are either ATLAS IP (primary production) or the rock close to FASER (secondary production).
• The first piece of material that we include in our simulations is a region with 10 m of concrete shielding that starts about 100 m in front of the tunnel TI12, where FASER detector is located. We model this shielding as material with Z = 8.6, A = 17 and ρ = 2.3 g/cm 3 .
• Concrete shielding is followed by 89 m of the standard rock, modeled as material with Z = 11, A = 22 and ρ = 2.65 g/cm 3 [104].
• The remaining section of the rock with the length of 1 m, as well as a 30 cm long section beginning right after the rock that contains i.a. the first front veto layers, are treated separately. In particular, any secondary production processes happening here are rejected in the analysis as they might lead to veto activation. It is useful to note that 1 m of the standard PDG rock corresponds to roughly three hadronic interaction lengths, 3 λ had,int .
• Located behind the first front veto is FASER's neutrino detector, dubbed FASERν [58]. It contains 1 m of tungsten out of which the first 70 cm of tungsten are taken into account when modeling secondary production. This is modeled as material with Z = 74, A = 183.84, ρ = 19.3 g/cm 3 . Any scatterings that occur in the remaining 30 cm of tungsten, as well as within further 20 cm of FASER/FASERν equipment containing i.a. second section with front veto layers, are not taken into account when obtaining sensitivity plots. We note that FASERν detector consists of alternating layers of tungsten (with the total length of about 1 m) and emulsion films. This complicated geometry is simplified for our phenomenological study in which we model FASERν as a single tungsten block with an average position with respect to other elements of the infrastructure around. Precise details of this modeling have minor impact on the final results.
• After FASERν and the second front veto, the FASER decay vessel begins at a distance L 480 m away from the ATLAS IP. The decay vessel has a cylindrical shape with radius R = 10 cm and length ∆ = 1.5 m.

FASER 2 (HL-LHC):
• FASER 2 is characterized by a similar geometry to FASER, but without FASERν neutrino detector in front of the decay vessel. As a result, after the concrete shielding (10 m), rock (89 m+1 m) and the first front veto layers (30 cm), the decay vessel begins at a distance L 480 m away from the ATLAS IP. The decay vessel is a cylinder with radius R = 1 m and length ∆ = 5 m.
MATHUSLA (HL-LHC): (beginning with most distant places away from MATHUSLA) • The size of the CMS detector and its cavern is taken into account in the analysis in a simplified manner. In particular, we do not treat secondary production processes happening in the elements of the infrastructure there as it would be challenging to model them accurately. On the other hand, this region lies close to the CMS IP and the relevant production rate for LLPs that can reach MATHUSLA is dominated by the primary production processes. The CMS cavern has roughly a cylindrical shape with 13 m radius and the length of 60 m [110]. By taking into account planned position of the MATHUSLA detector and the size of the CMS cavern, one can see that about 17 m to 30 m of the distance between the CMS IP and the bottom part of the MATHUSLA detector (details below) will be occupied by the cavern.
• The rest of a distance between the CMS IP and MATHUSLA is modeled as the standard rock with Z = 11, A = 22 and ρ = 2.65 g/cm 3 [104] that extends up until the last 1 m long section of the rock in front of MATHUSLA. The last 1 m long section of the rock is, instead, treated separately, similarly to the aforementioned case of FASER/FASER 2. This is to make sure that secondary production processes that we take into account are shielded from veto layers by at least 3 λ had,int .
• MATHUSLA decay vessel begins right after the rock. The total distance between the CMS IP and MATHUSLA, as well as the size of the detector, is dictated by geometry described by Eq. (11).
SHiP: (beginning with most distant places away from the SHiP decay vessel) • The target and hadron stopper with additional shielding correspond to about 13.2 m of infrastructure elements that we neglect when modeling secondary production. This is justified by the fact that the dominant contribution to the signal yield from this region corresponds to primary production of LLPs initiated by protons hitting the target.
• After this front section, the remaining part of the active muon shield consists of 6 iron blocks each 5 m long with 10 cm empty gaps in between, where iron is modeled as material with Z = 26, A = 55.84 and ρ = 7.874 g/cm 3 .
• This is followed by an additional section with the length of about 1 m which contains other elements of the infrastructure that are not modeled in details in our analysis.
• The SND detector begins at a distance about 44.7 m away from the IP. In its internal 80 cm × 80 cm part around the collision axis, it corresponds to the emulsion detector and the following downstream tracker. The emulsion detector is surrounded by a magnet with a coil made out of copper (alternative design with aluminium is also considered), which is further surrounded by an iron yoke. For the purpose of modeling, we assume that the region of the SND surrounding the emulsion detector is mostly filled with material with properties similar to iron.
The emulsion detector and downstream tracker: The total length of the emulsion detector is 3 m and it consist of alternating layers of 1 mm thick lead plates and emulsion films that are additionally interleaved with electronic detectors. We model this as a single block of lead with the total length of about 1 m that corresponds to 19 walls that each has 57 lead plates, as discussed in section 4.2 of Ref. [60]. The second half of this block of lead (i.e. 50 cm of lead) is excluded from our analysis of secondary LLP production as it corresponds to the last 3 λ had,int in front of the downstream tracker and the following SND muon system acting as veto. In our simplified treatment of detector geometry, we assume that the block of lead occupies 1 m long section in the middle of a 3 m long emulsion detector. Hence, out of the emulsion detector with the total lenght of 3 m, only secondary production happening in the front 50 cm long part of the lead block is taken into account which is positioned in between 1 m and 1.5 m away from the beginning of the SND detector. It is modeled as material with Z = 82, A = 207.2 and ρ = 11.35 g/cm 3 . After the emulsion detector, the remaining 3 m of the internal 80 cm × 80 cm part of the SND detector is occupied by the downstream tracker. We do not study secondary production processes that occur here. Importantly, for the purpose of studying additional LLP signatures relevant for scatterings off electrons in the SND, that are discussed in Sec. VI, we model the emulsion detector as a full block of lead with 1 m length. Magnets surrounding the emulsion detector and the downstream tracker: As discussed above, we assume that the magnets are made out of material equivalent to iron with Z = 26, A = 55.845 and ρ = 7.874 g/cm 3 . They have a total length of about 6 m. We again exclude a section of the magnet towards the end with 0.5 m length, as it corresponds to the last 3 λ had,int before the SND muon system acting as veto. On the other hand, the first 5.5 m long section of the magnet is taken into account when modeling secondary production.
• The SND downstream tracker and surrounding magnet is followed by the 2 m long muon identification system of the SND. It partially acts as the front veto before the decay vessel begins, so we do not take into account any possible secondary production processes happening in the close vicinity or inside the SND muon system.
• The muon identification system is followed by the decay vessel that has length of 50 m. We take into account its irregular shape [60] when simulating the events.
When obtaining the sensitivity reach plots presented in Sec. V, we impose additional cuts that can limit the number of signal events, but, most importantly, allow one to discriminate between signal and background events more easily. We summarize these cuts briefly below for both two-body decays with all LLP energy going into visible charged tracks, as well as for three-body decays with missing energy in the final state. The latter is relevant for searches for inelastic DM. In the analysis, we also assume a 100% detection efficiency after the relevant cuts are imposed.
On top of these experimental cuts, we always require the recoil momentum of the target in secondary production processes not to exceed p recoil = 1 GeV. Similarly, no secondary production events that occur within the last three hadronic interaction lengths, 3 λ had,int , in front of veto are taken into account. We note that detailed detector simulations could allow one to alleviate these constraints. In particular, similar challenges have already been considered in the context of LLP searches in high-energy electron and muon beam-dump experiments, e.g. NA64-e [111] and NA64-µ [112] employing ∼ 100 GeV e/µ beams.

FASER/FASER 2:
We require the LLP decays to happen inside the decay volume and deposit at least 100 GeV of energy in two final state visible charged tracks. Due to large energy threshold, no additional direction reconstruction cut is imposed since charged particles produced in the decay are highly collimated along the LLP direction in the lab frame.

MATHUSLA:
We follow Ref. [12] and, for LLP decays happening in the decay volume, we implement lower threshold on three-momenta of both charged particles in the final state, | p e | > 1 GeV, that is relevant for e + /e − charged tracks. On top of this cut, in order to suppress various possible sources of BG, additional timing and pointing information can be used to correlate the event with pp collisions at the CMS IP.
As for the latter, due to the small momentum exchange with the target in scattering processes, the direction of LLPs entering the decay volume is not much affected by the recoil and this effect can be neglected in the analysis (see also discussion in Appendix C). The direction reconstruction of decaying LLPs can also be affected in the presence of missing energy in the final state, especially for low-energy signal events. As a result of such misreconstruction, isotropic BG induced by interactions of atmospheric neutrinos might become non-negligible as it corresponds to few tens of events per year [10]. This issue has been studied in Ref. [113] for LLPs produced in Drell-Yan processes in pp collisions at the CMS IP, as well as produced in invisible decays of the SM Higgs boson. In both cases, the sensitivity reach has been found to decrease by a factor of 3 for new particles with mass m LLP ∼ 10 GeV, once neutrino-induced BG is taken into account. On the other hand, as noted in Ref. [113], this effect could become even less sizable when dedicated background simulations are performed. In this study, we will then neglect possible negative impact of direction misreconstruction in searches for inelastic DM.
SHiP: We follow Ref. [60] and require three-momenta of both visible charged particles produced in LLP decays in the decay vessel to be large enough, | p daughter | > 1 GeV. Similarly to the case of MATHUSLA, we note that direction reconstruction is not much affected by the recoil in scattering processes relevant for secondary production of LLPs.
In case of 3-body decays with missing energy, loose selection efficiency could be imposed with impact parameter with respect to target not exceeding 250 cm [60]. While this condition allows one to more easily accept partially reconstructed signal events, it could also make BG rejection more challenging. Dedicated BG simulations that go beyond this analysis are required to thoroughly investigate this issue which we left for future studies. Instead, in our phenomenological analysis, we neglect this possible obstacle and assume that all such events can be properly reconstructed and discriminated from BG. This also allows for a better comparison to the reach of MATHUSLA experiment for which we make a similar assumption.
Experimental cuts applied to additional searches based on LLP scatterings off electrons in the SND are described in Sec. VI for both considered signatures.