Probing light dark scalars with future experiments

We investigate a dark sector containing a pair of light non-degenerate scalar particles, with masses in the MeV-GeV range, coupled to the visible sector through heavier mediators. The heaviest dark state is long-lived, and its decays offer new testable signals. We analyze the prospects for detection with the proposed beam-dump facility SHiP, and the proposed LHC experiments FASER and MATHUSLA. Moreover, we consider bounds from the beam-dump experiment CHARM and from colliders (LEP, LHC and BaBar). We present our results both in terms of an effective field theory, where the heavy mediators have been integrated out, and of a simplified model containing a vector boson mediator, which can be heavy $\gtrsim\mathcal{O}(1)$ TeV, or light $\mathcal{O}(10)$ GeV. We show that future experiments can test large portions of the parameter space currently unexplored, and that they are complementary to future High-Luminosity LHC searches.


Contents
1 Introduction and framework Although elusive, it is a concrete possibility that dark sectors exist in nature. Their physics is interesting for various reasons: they may appear in connection to dark matter (DM) [1], they could produce new signals, such emerging jets at colliders (see e.g. hidden valley models [2][3][4]), and they might be connected to the solutions of some open questions of particle physics, for instance the hierarchy problem (e.g. in Twin Higgs [5] and relaxion [6] models). Several new experiments have been proposed in the recent years to target such sectors [7][8][9][10][11][12][13][14][15][16], and investigate possible extensions of the Standard Model (SM) of particle physics. Particular attention have been devoted to light dark states, with masses in the MeV to GeV range. These scenarios can be explored with experiments both at the intensity frontier [17] and at the energy frontier, like the Large Hadron Collider (LHC), with a nice interplay between the two experimental programs.
Dark sector physics can be broadly characterized in terms of the particle content in the dark sector itself, and the mediators that connect them to the SM. The interactions between the dark states and the SM can thus be written in terms of renormalizable or non-renormalizable portal operators. The latter option is appropriate when the mediators are heavier than the energy scale involved in the process under consideration, like the production of dark sector particles in an experiment. In this case, the mediators can be integrated-out, and the physics is described in terms of contact non-renormalizable operators. The advantage of this description in terms of an effective field theory (EFT) is its model-independence: it captures at the same time different microscopic ways in which the dark sector communicates with the SM. When the mediators are light they can be produced on-shell in the experiments considered. A small drawback is that the description, now usually in terms of renormalizable operators, is more model dependent. On the positive side, processes in which the mediator is produced on-shell can now be used to study the scenario, a possibility that was precluded in the EFT description.
In this paper we will consider both cases in the context of a dark sector containing non-degenerate scalar states. As we are going to see, this framework has a rich phenomenology, since the decays of the heavier state open up new experimental venues for detection. A similar situation with non-degenerate fermion states has been considered in ref. [18]. The prototypical model that we have in mind is given by a complex scalar singlet (1.1) whose real and imaginary components are mass splitted. The mass splitting can be generated dynamically, and although in the following its origin will not be important, it is useful to consider a concrete example 1 . Take for instance the potential V = −µ 2 H |H| 2 + λ H |H| 4 + µ 2 φ |φ| 2 + λ φ |φ| 4 + λ Hφ |H| 2 |φ| 2 + Bφ 2 + h.c. (1.2) In the B = 0 limit the U (1) symmetry of the potential forces φ 1 and φ 2 to be degenerate (unless the U (1) is spontaneously broken, a situation that we will not consider in this paper). Once B is turned on, it is easy to see that a mass splitting is obtained. To quantify the mass difference we introduce the parameter where m 2 and m 1 are the masses of the heavier and lighter states, respectively. Generically it is technically natural to have a small δ, since in the δ = 0 limit a symmetry forcing the two components of the complex singlet to be degenerate is recovered. The phenomenology of scalar states interacting with the SM via the Higgs portal coupling λ Hφ have been studied thoroughly in the literature (see for instance [19] for a recent review). In the following we will thus assume that the Higgsportal coupling λ Hφ in eq. (1.2) is negligible, with the dark sector-SM interactions mediated by some new particle. As an example, we can imagine a vector boson Z interacting via the Lagrangian where f L and f R are the left-handed and right-handed SM fermions, and the dark current is given by Another possibility is to introduce vector-like fermions Ψ f L,R interacting with the SM via Our focus here is on mediators with masses above few tens of GeV. The phenomenology of a light (below the GeV) mediator is different, and it has been studied elsewhere, see e.g. [9,11,14,[20][21][22][23]. Assuming the mediators to be heavy allows to give a unified treatment of the phenomenology at low energy. Indeed, integrating out the vector in eq. (1.4) or the fermions in eq. (1.6) we obtain effective operators of the form [24]: where the dots represent other operators that are generated. 2 Notice that we write the operators in a SU (3) c × SU (2) L × U (1) Y invariant manner. In terms of vector and axial currents they are defined as where the sum is extended over all SM fermions. The Wilson coefficients are simply (1.9) 2 In the case of heavy fermions, gauge invariance allows for Yukawa interactions of the type gΨ f L HΨ f R in addition to the operators in eq. (1.6). These interactions would generate additional dimension-6 operators of the form φ † φf L Hf R which may contribute to the phenomenology of the model. In the following we will suppose the Wilson coefficient of such operators to be suppressed. This can be achieved for instance requiring Minimal Flavor Violation, in which case the Wilson coefficient is proportional to the fermion mass and is typically suppressed.
Analogously, the vector and axial couplings of the Z vector with the SM fermions are , . (1.10) For definitiveness, in this paper we will consider only the case c f L = c f R for the EFT operators, and g f L = g f R for the couplings of the Z boson. This means that the SM fermions interact with the dark scalars only through vector currents, with the exception of the neutrinos, which couple through the usual V-A current. We leave a detailed analysis of the phenomenology of axial interactions to future work. Notice also that this choice allows us to neglect SM radiative corrections that could otherwise be important. In fact, it has been shown that, for vector currents, the renormalization group evolution of the effective operators from the scale Λ to the scale of low energy experiments is small [25][26][27][28][29][30][31][32][33][34][35][36]. Therefore, we will neglect such effect.
The operators introduced above control the production of the dark states from the scattering of SM particles, and the decay of φ 2 into φ 1 and SM states. The latter process is suppressed in the limit of a small mass splitting δ, and/or feeble interactions among the SM and the dark sector. Interestingly, this opens the possibility that φ 2 is long-lived, and travels macroscopic distances before decaying. Such scenario can be tested in experiments where the dark states are produced in high-intensity or high-energy facilities, and then show up in a far-placed detector through the decays of φ 2 . The aim of this paper is to explore the sensitivity of such experimental facilities for the scenario presented above. Concretely, we will consider fixed target experiments, and proposals for detectors to be placed near the LHC interactions points. For the first class of experiments, we will compute the constraints from searches of the CHARM experiment [37], and the projected sensitivity for the proposed SHiP facility [7,8]. We will then focus on the LHC experiments FASER [9,10] and MATH-USLA [11,12]. We will also compare the sensitivity regions of these experiments with the collider bounds from searches at LEP, BaBar, and LHC (namely LHCb, ATLAS and CMS).
In our work we assume that φ 1 is stable or very long-lived, such that once produced, it decays far outside the detectors that we are considering. Although not essential, it is interesting to conceive the possibility that it is a DM candidate. We will briefly comment about this option later.
Before entering into the details of the calculations, let us explain the structure of the paper. In section 2.1 and 2.2 we describe how the the constraints and future sensitivities are computed in fixed target experiments and LHC experiments. In section 2.3 we present collider bounds which do not rely on the inelastic nature of the dark sector. More specifically we will discuss constraints from searches of missing energy events at LEP, LHC and BaBar, as well as limits from the invisible decays of heavy quarkonia states. In section 3 we present our results in the context of the EFT  branching ratio of φ 2 decays into φ 1 and SM states. In both panels we take c f L = c f R = 1, δ = 10 and Λ = 1 TeV. In the right panel the threshold at m 1 0.2 GeV denotes the transition of the treatment of φ 2 hadronic decays from chiral perturbation theory to perturbative QCD. See appendix A for details. operators in eq. (1.7). As mentioned before, we will consider mediators with masses above few tens of GeV. If light enough, these particles can be produced on-shell at LHC. Motivated by this consideration, in section 4 we will move to the case of the Z model in eq. (1.4). This will allow also a more precise comparison with current LHC searches. We will consider the case of a heavy Z , with a mass above few TeV, and the benchmark case of a dark photon with a mass of 40 GeV. In section 5 we present cosmological and astrophysical constraints. Part of this discussion applies only when φ 1 is the dominant component of DM. Finally, we conclude in section 6. We also add two appendices. In appendix A we collect the decay widths used in the EFT analysis, explaining in detail how they are obtained. In appendix B we instead present the decay widths of the Z boson, used for the computations in section 4.

Overview of experimental bounds and forecasts
In the following two sections we examine current and future experiments that can test the interactions in eq. (1.7) via φ 2 decays, i.e. in which the inelastic nature of the dark sector is instrumental in probing the parameter space. We classify them according to their center-of-mass energy: √ s 28 GeV (section 2.1) and √ s = 14 TeV (section 2.2). The physical processes of interest can be summarized as If the decay of the heaviest state φ 2 occurs inside the decay volume of the experiment, a signal event may be detected. φ 1 φ 2 pairs can be produced directly from parton collisions, and in the decays of the mesons generated by the proton interactions 3 .

Experiments at √ s 28 GeV: CHARM and SHiP
The CHARM [37] and the proposed SHiP [7,8] experiments exploit a 400 GeV proton beam impinging on a fixed target, i.e. √ s 28 GeV. For these experimental facilities, we compute the fluxes of dark scalars generated by meson decays. We do not consider the production of φ 1 φ 2 pairs from parton level processes for reasons that will be explained in sec. 2.2. We include in our analysis light neutral pseudo-scalar mesons (π 0 , η and η ), neutral vector mesons (ρ and ω) and the charmonium J/ψ and the bottomonium Υ(1S) resonances.
The light pseudo-scalar and vector mesons are more abundantly produced by the proton interactions than the heavier ones. Nevertheless, the latter tend to have larger branching ratios into dark states, that may compensate for the smaller production cross sections. We will come back to this point later. Furthermore, heavy mesons are obviously relevant for heavy enough dark states, whose production from light mesons is kinematically forbidden. We show in the left panel of fig. 1 the branching ratios for the decays of the mesons into φ 1 φ 2 states, for the following benchmark case: Λ = 1 TeV, δ = 10. In this plot, and in the rest of the paper for the discussion of the EFT operators in eq. (1.7), we adopt a democratic vector coupling to all the SM fermions, specifically c f L = c f R = 1 for all the SM fermions (of course c f R = 0 for the neutrinos). Detailed formulas for the decay widths can be found in appendix A. For vector interactions, the dominant decay mode of the light pseudo-scalar mesons into dark sector particles is via the three-body process M → φ 1 φ 2 γ. Notice also that the branching ratio for ρ → φ 1 φ 2 is suppressed with respect to the one for the analogous processes involving the other vector mesons (ω, J/ψ and Υ). This is because the effective coupling of the ρ meson to the dark current is suppressed, for the choice of couplings above (see eq. (A.9)). As mentioned before, J/ψ and Υ have larger branching fraction into the dark states than the lighter mesons.
Our goal is to compute the number of signal events in a detector. This can be done using: where N φ 2 M is the number of φ 2 particles produced in the decays of the meson M, f dec M is the fraction of events in which φ 2 decays inside the detector, and f rec M is the efficiency for the reconstruction of the signal event in the detector. We are going to discuss how to compute each term.
Production from meson decays. The number of dark particles produced in the decays of the mesons M = {π 0 , η, η , ρ, ω, J/ψ, Υ} can be computed according to where N POT is the number of protons on target collected (for CHARM) or expected (for SHiP) by the experiment, N M the average number of mesons produced per proton interaction, and BR(M → φ 1 φ 2 + X) is the branching ratio of the meson M into the dark states. The number of protons on target for the CHARM and SHiP experiments [7,8,38] are summarized in table 2. The average number of mesons produced per proton interaction, N M , is computed in the following way. We simulate pp collisions using two different softwares: EPOS-LHC [39] (as found in the CRMC package [40]) for light mesons, and PYTHIA8 [41] for the J/ψ and Υ mesons. Since EPOS-LHC has been tuned to correctly reproduce mesons multiplicities and energy distributions of LHC data, it is particularly useful for our purposes. For the heavy mesons we use PYTHIA8, which allows to switch on only the charmonium or bottomonium production processes 4 . We find the total production cross-sections σ J/ψ 250 nb and σb b 1.3 nb, in reasonable agreement with the literature [42][43][44]. From the results of our simulations, and using a total proton-proton cross section σ pp 40 mb [45], we compute the average number of mesons produced per proton interaction, reported in table 1. For light mesons, we have checked that similar multiplicities are obtained using PYTHIA8. A comparison between the output of PYTHIA8 and experimental data [46] for the production of π 0 and η mesons can be found in [47]     Finally, the branching ratios BR(M → φ 1 φ 2 + X) are obtained using the equations in appendix A. They are shown in the left panel of fig. 1 for the same benchmark case discussed before. In fig. 2 (left panel) we show the expected number of φ 2 particles at SHiP. The production via mesons decays is dominated by the ω contribution for small dark scalar masses, and by the J/ψ and Υ decays for larger ones. Notice that the contributions from J/ψ and Υ exceed those from the light mesons, except for the ω. As anticipated, this is due to their larger branching ratio into dark states.
Let us also mention that the relative contributions are different when the mediator of the interaction between the dark and visible sectors is light so that it can be produced on-shell, see e.g. [14,20].
Dark state decays. The fraction of φ 2 states which decay inside the detector can be computed as: The first factor is a geometrical cut and represents the fraction of φ 2 particles whose trajectories intersect the detector. The second term takes into account the probability that those particles decay inside the detector. It is computed in terms of the distance between the interaction point and the detector (L), the detector size (L dec ) and the of hadrons in the target [20]. For the scenario that they are considering, they find that cascades affect the signal rate by 15-40% for π 0 decays, while the impact is negligible for the case of η mesons.
decay length of φ 2 in the laborarory (LAB) frame, given by Here τ φ 2 is the decay time of φ 2 , while β φ 2 and γ φ 2 are respectively its speed in units of speed of light (c), and its Lorentz factor, in the LAB frame. The decay time τ φ 2 is obtained using the equations of appendix A. The additional quantities are obtained as follows. We consider the samples of mesons obtained with EPOS-LHC and PYTHIA8. For each event, we simulate the M → φ 1 φ 2 + X decay in the rest frame of the meson, and then compute the 4-momentum of φ 2 in the LAB frame performing the appropriate Lorentz boost. From this, one can compute the β φ 2 and γ φ 2 factors. Then, we select only those events whose trajectories intersect the detector, which allows us to compute f geom M . For this purpose, we approximate the detectors as cylinders, and we impose a cut on θ φ 2 , the angle between the direction of flight of φ 2 and the beam axis. In the case of SHiP the area of the cylinder is taken to be A dec = 5 × 10 m 2 [7,8]. The maximum opening angle θ φ 2 is thus tan (θ dec ) = A dec /π/(L + L dec ). For CHARM we follow ref. [49] to take into account the fact that the detector is off-axis. We summarize all the relevant quantities in table 2. The second term in eq. (2.4) is computed averaging ( · ) the probabilities over all the events in the direction of the detector.
Reconstruction of the events in the detector. To mimic event selection cuts of experimental analysis, we impose a requirement on the energy of the visible particles produced in φ 2 decays. More specifically we adopt the following simplified strategy. We focus on the dominant decays φ 2 → φ 1 + V (where V is ρ or ω) and φ 2 → φ 1 + ll (where l is a charged lepton), see fig. 1. 6 From the sample of φ 2 events produced in the decay of a given meson M , we compute the energy of the meson V or of one of the charged leptons l produced in the decay. Then, we select the events where the energy of these visible particles is larger than a certain threshold E cut to be specified below. We compute the fraction of events which satisfy this energy requirement, f rec,i M , for each of the φ 2 decay channel i under consideration. Finally, we define the efficiency of this selection cut as where BR(φ 2 → φ 1 + X i ) is the branching ration for the decay of φ 2 in the channel i. In the equation above, we have also introduced an efficiency i for the reconstruction of the visible states in the decay. For CHARM, we recast a search of heavy neutrinos decaying into light neutrinos and an electrons or muons pair [38]. We impose a cut E cut = 2 GeV and we adopt the efficiencies for the the electron and muon channels in table 2. For SHiP, instead, we simply assume E cut = 2 GeV and a 100% efficiency Table 2.
Specifications used for each experiment to compute the total number of φ 2 signals in the detector according to eq. (2.2): N POT is the number of proton on target (for CHARM and SHiP), L is the total luminosity of the experiment (for FASER and MATHUSLA), θ φ 2 is the angle between the direction of flight of φ 2 and the beam axis, L is the distance between the interaction point and the detector, L dec is the detector length and is the efficiency for detection of the φ 2 decay products. For each experiment we also specify the production mechanism for φ 2 that we have considered (mes = meson decays, pp = direct parton production). (*) In the case of the CHARM experiment the efficiencies must be multiplied by a factor 0.095 to take into account that the detector covers only 9.5% of the circular corona defined by the angular cut.
in all channels. The branching ratio of φ 2 in the relevant channels can be computed using the formulas in appendix A.

Experiments at √ s = 14 TeV: FASER and Mathusla
We consider two proposed experiments that aim to leverage the proton collisions occurring at √ s = 14 TeV at LHC: FASER [9,10] and MATHUSLA [11,12]. The FASER detector is planned to be placed downstream of the ATLAS experiment along the beam axis. MATHUSLA instead will be located on the surface, near the CMS or ATLAS interaction point. We summarize the FASER geometry in table 2 (we have considered the design of FASER-2, which will have a cylindrical detector of radius of 1 m [10]). See instead eq. (2.8) below for a discussion of the configuration of the MATHUSLA detector. In the following we consider the production of dark particles in meson decays or directly from parton collisions.
Dark particles production in mesons decays has already been discussed in sec. 2.1. We follow the same procedure also for the experiments at √ s = 14 TeV. The average number of mesons produced in the collisions has been computed using EPOS-LHC and PYTHIA8 and is summarized in table 1. For the production of the J/Ψ and Υ mesons, we have introduced a correction factor to the output of the PYTHIA8 simulations, to match the production cross-sections measured by the LHCb  We take c f L = c f R = 1 in both panels. Green (red) lines are for the production from meson decays (direct parton production) only. The black lines include both production processes. Gray regions are excluded by the constraints discussed in sec. 2.3.
collaboration [50,51]. It is 1.6 and 0.25 respectively for the J/Ψ and Υ mesons. A small modification of eq. (2.3) is needed to compute the total number of mesons produced at the experiment. The equation now reads where L is the integrated luminosity collected by the experiment and σ tot pp the total proton-proton cross section at √ s = 14 TeV. Using EPOS-LHC we obtain σ pp 110.69 mb, in agreement with the experimental measurement [52] (which has been performed at √ s = 13 TeV). The number of signal events in the detector is computed using eqs. (2.2), (2.4), (2.5) and (2.6), and following the procedure of sec. 2.1.
Dark particles production from parton collisions is computed as follows. We implement the effective operators in eq. (1.7) in Feynrules [53] and use MadGraph5_aMC@NLO [54] to simulate pp → φ 1 φ 2 events at √ s = 14 TeV. From this sample, we can compute the number of signal events in the detector using the same procedure already outlined in sec. 2.1. More in detail, eqs. (2.4) and (2.5) can be applied directly also in this case. As for the number of φ 2 produced in the parton collisions, it is obtained from: where the σ pp→φ 1 φ 2 cross section is computed with MadGraph5_aMC@NLO. We must, however, ensure that the EFT in eq. (1.7) is used inside its domain of validity. Any EFT is a valid description of a more fundamental theory only for processes occurring at energy scales smaller than the cut-off of the EFT, M cut . The latter can be written as M cut = g * Λ, where g * is a combination of couplings of the UV theory, for example g * = √ g φ g V f for the Z model in sec. 1 (we remind that we are taking c f R = c f L = 1 for all the fermions). As an example of a weakly coupled theory we take g * = 1. Following [55,56], we then impose √ŝ ≤ Λ, with √ŝ the center of mass energy of the partonic event. In other words, from our sample of events simulated with MadGraph5_aMC@NLO, we select only those satisfying this requirement. Of course, this cut is not needed when we consider the Z model in sec. 4. In that case, apart from this difference, we implement the operators in eq. (1.4) in Feynrules and we compute the signal events following the same steps explained above.
Another important concern is the uncertainty in the parton distribution functions, and even the use of perturbative QCD to describe the production process. This point can be understood as follows (see ref. [9] for a discussion). For small opening angles θ φ 2 , (i.e. if the detector is small in the transverse direction or placed at a large distance from the interaction point) only φ 2 particles with low transverse momenta reach the detector. This implies that one is interested in events with small momentum transfer, and the parton distribution functions has to be evaluated at low factorization scales Q and momentum fractions x, in a regime where they suffer from large uncertainties, or even the description of the hadrons in terms of partons in perturbative QCD breaks down. Inspecting table 2 we see that CHARM, SHiP and FASER have very small angular acceptances, and we have checked that for these experiments only events with relatively small transverse momentum ( 2 GeV) are involved. For this reason, we do not include direct parton production in the computation of the signal events for those experiments.
The case of MATHUSLA is different: given its location, only particles with relatively large transverse momentum reach the detector on the surface. Therefore, we consider φ 2 production both via meson decays and parton collisions. We show the relative importance of the various contributions in the right panel of fig. 2. As it can be seen, for MATHUSLA parton production dominates. The geometry of the detector is the most recent one in ref. [12]. The detector is delimited by: where the coordinate system is centered at the LHC interaction point, the z axis is along the beam direction, and x denotes the vertical to to the surface.
We conclude mentioning that, for both FASER and MATHUSLA, the cut on the energy of the φ 2 decay products discussed in sec. 2.1 is E cut = 2 GeV, and we assume an efficiency of reconstruction of 100% for all channels.

Collider searches
In this section we discuss constraints from accelerators which do not rely on the inelastic nature of the dark sector. More specifically we will discuss bounds from searches of missing energy events at LEP, LHC and BaBar, as well as limits from the invisible decays of heavy quarkonia states.
LEP can be used to test our scenario in different ways: (i) using the excellent measurement of the Z decay width to constrain exotic decay modes, and (ii) from searches of mono-photon events with large missing energy. We start by considering the first possibility. The EFT of eq. (1.7) induces the 4-body decay Z → φ 1 φ 2f f and the invisible decay Z → φ 1 φ 2 via a fermion loop. The corresponding decay widths scale as Λ −4 and are thus suppressed at large Λ. We have used MadGraph5_aMC@NLO to compute the 4-body process. Focusing on values of Λ where the EFT is valid, say M Z (see discussion in sec. 2.1), we found that the decay width of the exotic process is always at least three orders of magnitude smaller than the current uncertainty on the Z width. Therefore, we can simply neglect this constraint. Analogously, also the Z radiative decay into φ 1 φ 2 gives a very weak bound [57].
A more interesting limit is obtained from searches of mono-photon events. We consider the analysis of ref. [58], where bounds from these searches at LEP II have been used to constrain the interaction of fermionic DM with the SM, in the context of an EFT approach. We recast these results for our effective operator in eq. (1.7) with the following simplified method. In ref. [58] a bound Λ > 500 GeV has been obtained for the vector operatorχγ µ χēγ µ e/Λ 2 coupling electrons to a DM candidate χ lighter than (60 ÷ 70) GeV. To estimate the bound in our case, we simply factorize the production cross section as σ(e + e − →χχγ) σ(e + e − →χχ) R(e → eγ). We take the function R(r → eγ) to be universal, and ignore for simplicity the difference in the angular distributions between the fermionic and scalar cases. From the scaling σ(e + e − →χχ) ∼ 1/Λ 4 and the ratio of the cross sections of the processes e + e − →χχ and e + e − → φ 1 φ 2 , we obtain the constrain Λ 350 GeV. This is reported with a solid gray line in figs. 3 and 4. Notice however that, as already mentioned in sec. 2.2, one should ensure that the EFT is applied inside its domain of validity. Following the previous section and taking M cut = Λ, we conclude that the constraint derived above for the EFT can not be consistently applied for Λ smaller than the center of mass energy at LEP, i.e. E cm 200 GeV. This lower limit is depicted with a dashed gray line in figs. 3,4. To test small values of Λ, outside the validity of the EFT, it is necessary to specify the microscopic origin of the EFT. For instance, ref. [58] has explored models where the EFT operators arise from the exchange of an s-channel mediator. We will come back to this case in sec. 4. Another issue that must be taken into account is that, for the mono-photon search to apply, we need to make sure that φ 2 decays outside the detector. To estimate the region of the parameter space where this happens we simulate e + e − → γφ 1 φ 2 events with a center-of-mass energy of 200 GeV, and compute the average proper decay length of φ 2 in the LAB frame: cτ φ 2 γ φ 2 β φ 2 . Requiring a proper decay length larger than the size the detector, which we take to be 5 m, we obtain an upper limit m 1 0.06 GeV (2.9 GeV) for δ = 10 (0.1). These cuts on the mono-photon bounds are visible in figs. 3,4. We shall now consider LHC searches. The problem of the validity of the EFT arises also in the interpretation of missing energy signatures at the LHC. Still, conservative but consistent constraints can be obtained including in the analysis of LHC data only those signal events with a center of mass energy below M cut [55]. This procedure has been applied in [56] to recast mono-jet searches in the framework of the EFT of a DM singlet fermion, finding that for g * = 1 no bound could be set, while some region of the parameter space was probed assuming the large coupling g * = 4π, representative of a strongly coupled UV completion. A similar calculation could be repeated for our scenario with current LHC data. However, for the analysis of the LHC searches, we prefer to resort to a simple UV completion in sec. 4, which will allow a more pertinent investigation of current collider bounds.
Let us now move to colliders working at lower energies. We consider searches performed by BaBar in e + e − collisions at √ s 10 GeV, near the resonances Υ(2S), Υ(3S) and Υ(4S) [59]. Invisible decays of heavy quarkonium states is a handle to probe light dark particles [60,61]. The current upper limit on the invisible decay of the resonance Υ(1s) has been obtained by the BaBar collaboration, and reads BR(Υ(1s) → inv) < 3 × 10 −4 [62]. Using the decay width of Υ(1s) into φ 1 φ 2 computed in appendix A, this leads to the constraint shown in figs. 3,4. As for the case of mono-photon searches at LEP, we need to ensure that φ 2 decays outside the detector. We use a procedure similar to the one described above, imposing for φ 2 a proper decay length in the LAB frame larger than 3 m. The effect of this requirement can be seen in the vertical cuts in the Υ(1S) → inv bounds in figs. 3,4. A less stringent bound is obtained from the upper limit on the invisible decay of the J/ψ [63]. Also searches for mono-photon events and large missing energy in e + e − collisions could put bounds on the Wilson coefficients of the model. In our scenario, the corresponding bound turns out to be less stringent than the one derived above from the Υ(1s) invisible decay, in the kinematical range where the latter applies, i.e. m 1 + m 2 < M Υ [57].

Sensitivities on the EFT operators
We are now in a position to discuss how the experiments described in sec. 2 can be used to exclude or probe the parameter space of the EFT defined in eq. (1.7). Let us remind that we are taking the following choice of vector democratic couplings: c f L = c f R = 1 for all the SM fermions. The exclusion region from CHARM is obtained recasting the search in ref. [38] where no events have been observed. A limit at 90% CL can be obtained imposing that the number of signal events is N sig < 2.3. The same criterion is used to derive the sensitivity reach of the SHiP, FASER and MATHUSLA experiments. This assumes that backgrounds can be reduced at a negligible level, which is expected to be the case according to present studies [8,10,12]. Moreover, our sensitivity contours are only marginally affected by an imperfect reconstruction of the signal (i.e. assuming an efficiency i < 1), or a moderate number of background events. This is due to the steep dependence of the number of signal events on Λ. In fact, the number of φ 2 particles produced in meson decays or via direct parton collisions scales as Λ −4 . Moreover, the fraction of those events which decay inside the detector depends exponentially on the φ 2 lifetime, see eq. (2.4), which in turn grows like Λ 4.
Let us now present our results in the (m 1 , Λ) plane for the two representative cases of a large, δ = 10, and small, δ = 0.1, mass splitting among the dark states. First we show our sensitivity contours for MATHUSLA in fig. 3, including separately only events from meson decays (green line), direct parton production (red line), and then combining the two processes (black line). As expected from fig. 2 (right panel), parton production is generically more relevant. In particular it allows to extend the reach at masses of the dark states that can not be explored with meson production for kinematical reasons (simply because m 1 + m 2 is larger than the mass of the mesons). However, there are regions of the parameter space where the two processes give comparable sensitivities, or the production via meson decays is more relevant. The reason for that is the term in eq. (2.4), which accounts for the probability that the φ 2 states decay inside the detector. As discussed below eq. (2.4), it depends on the distance travelled by the φ 2 particle, and therefore its Lorentz factor γ φ 2 . Production from meson decays or via direct parton collisions lead to different distributions for γ φ 2 , centred around larger values for the latter process. Therefore the term in eq. (2.4) can be very different in the two cases. Notice also, from the right panel of fig. 2, that the number of dark pairs produced by heavy meson decays is larger than those from the lighter ones.
One can notice in fig. 3 that for a given mass m 1 , MATHUSLA is able to probe a limited range of Λ. For small values of Λ the φ 2 particles decay before reaching the detector, while in the opposite case their production is suppressed, and/or they decay at large distances. Moreover, for the direct parton production, small values of Λ are not included inside our sensitivity reach, see the right panel of fig. 3. This is because the cut introduced in sec. 2.2 to ensure the validity of the EFT approach, removes most of the signal events in this region.
In fig. 4 we compare the sensitivities of the different experiments that we have analyzed. The constraints discussed in sec. 2.3 are shown in gray. The left panel is for δ = 10. As evident, the CHARM exclusion region extends well above the bounds from LEP. FASER will be able to further improve these limits. SHiP and MATHUSLA can probe a much larger region of parameter space, reaching up to Λ (5 ÷ 8) TeV for m 2 ∼ (1÷3) GeV. In the case of SHiP (and less pronounced for FASER) a feature is visible around m 2 0.6 GeV, which corresponds to the threshold at which the main production channel for φ 2 moves from ω decays to J/Ψ decays. We now turn to the right panel, with δ = 0.1. Qualitatively, the results are similar to the previous case, but smaller values of Λ can be probed. The effect is particularly evident in the case of SHiP and, to a lesser extent, for CHARM. Consider a fixed value of Λ, for instance Λ = 2 TeV. With this small mass splitting, φ 2 has a large decay length, which causes most of the decays to happen after the detector, and diminishes the number of events. The decay length can be reduced increasing the mass of the scalars, but if they are too heavy they can not be produced in meson decays. Again, the features appearing in the contour regions of SHiP, FASER and CHARM, signal the transitions from ω, to J/Ψ and Υ decays as the main production mechanism for φ 2 . A difference with respect to the left panel is the range of φ 1 masses which can be probed. In this case, with a quite compressed mass spectrum, larger dark scalar masses are needed to lead a detectable signal through the process φ 2 → φ 1 + X i . For m 1 10 −2 GeV even the electron channel (X i = e + e − ) is closed, and φ 2 can only decay into neutrinos. To give an idea of the lifetimes that can probed with these experiments, we compute cτ φ 2 in the rest frame of the φ 2 particle for several choices of the dark particles masses and Λ. For δ = 10, we obtain cτ φ 2 = 5.6 × 10 4 m for (m 1 , Λ) = (2 × 10 −3 GeV, 100 GeV) and cτ φ 2 = 1.4 m for (m 1 , Λ) = (0.3 GeV, 6 TeV). Notice that it is trivial to rescale these results for other values of Λ, since the decay length simply scales as cτ φ 2 ∝ Λ 4 . Turning to δ = 0.1, we find cτ φ 2 = 2.5 × 10 5 m for (m 1 , Λ) = (0.1 GeV, 100 GeV) and cτ φ 2 = 4.2 m for (m 1 , Λ) = (10 GeV, 3 TeV). As evident, future experiments will be able to probe decay lengths that span many orders of magnitude. We remind that in the computation of the experimental sensitivities one should property take into account the relativistic γ φ 2 factor to determine the decay length of φ 2 in the LAB frame. This has been done as explained in sec. 2.1.
In the bottom panel of fig. 4, we the show a different slice of the parameter space. We fix the mass of the heaviest dark scalar M 2 = 1 GeV, and we present our sensitivities as a function of the mass splitting δ. The contours saturate at large values of δ simply because φ 1 can be considered effectively massless for large enough δ. Notice also that there is a region of the parameter space at small δ and Λ which is not covered by MATHUSLA. In this corner of the parameter space, the production of dark states via meson decays do not lead to detectable signals because the φ 2 decay products are too soft to satisfy our cut on their energy in sec. 2.2. The situation is different for FASER, CHARM and SHiP, since these detectors are placed along the beam axis, and they are typically crossed by φ 2 particles with larger energy, for which it is easier to produce decay products that satisfy the energy cut. For events produced directly from parton collisions, the cut on the center of mass energy of the partonic event, √ŝ < Λ (see sec. 2.2), becomes increasingly more stringent at small Λ, and selects events with softer φ 2 's. Then, these low energy events are typically removed by the requirement on the energy of the φ 2 decay products.
Let us close this section mentioning some other relevant experimental results. New limits on dark particles have been obtained by the MiniBooNE collaboration from a search performed with an 8 GeV proton beam dump [64]. We have recasted these results for our scenario with a simplified method to implement the experimental analysis cuts. We have found that the region of the parameter space excluded by this search (for the same benchmark scenarios in fig. 4) is already tested by the CHARM experiment. A more detailed analysis would be necessary to outline precisely the exclusion limits. Other complementary constraints could be derived from searches at LSND [65]. Searches of missing energy in the electron fixed target experiment NA64 [66] leads to weak bounds in our scenario, see [57]. We have also estimated the constraints from the electron beam dump E137 [67]. We have found that it probes a region of parameter space already excluded by CHARM and the collider bounds in fig. 4.

Sensitivities on the Z model
We now focus on a simplified model that provides a partial UV completion for the EFT defined in eq. (1.7). Specifically, we are considering the Z model introduced in sec. 1. Our aim is to properly study the sensitivities of the various experiments under consideration, in the case where the mediator of the interaction among the dark states and the SM can be produced on-shell. In this situation, as mentioned in sec. 2.3, the use of a simplified model allows us an in-depth comparison with the bounds from high-energy accelerators, in particular with those from the ATLAS and CMS experiments at the LHC. We consider two benchmark cases: a heavy Z (with a mass above 1 TeV) coupled to the SM fermions via vector currents, and a light dark photon (with a mass of 40 GeV) interacting with the SM via kinetic mixing (see for instance [68]).

A heavy Z
The simplified model that we are considering is specified by the lagrangian in eq. (1.4). At low energies, these interactions lead to the EFT in eq. (1.7), with the Wilson coefficients: where M Z is the mass of the Z boson. Since in this section we are considering M Z 1 TeV, for the CHARM, SHiP and FASER experiments we can apply the analysis on the EFT operators explained in sec. 3 and use eq. (4.1) 7 .
For the MATHUSLA experiment the situation is different, since the Z boson can be produced on shell at the LHC, a situation which is not captured by the  EFT. Therefore, we implement the interactions of eq. (1.4) in Feynrules [53] and use MadGraph5_aMC@NLO [54] to simulate pp → Z → φ 1 φ 2 events at √ s = 14 TeV. The analysis then proceeds as explained in sec. 2.2.
We focus on vectorial couplings, defined in eq. (1.10), and we consider two benchmark scenarios. In the first case, the Z is hadrophilic and its couplings to the SM quarks and leptons are respectively g Vq = 0.25, g V l = 0. In the second scenario an interaction with the leptons is turned on: g Vq = 0.1, g V l = 0.01. In both cases the interaction is flavor blind (same coupling for all the flavors), and g φ = 1. These choices are motivated by the possibility to confront with the LHC searches presented in ref. [69]. For the hadrophilic scenario, the strongest constraint is from searches of a new resonance using dijets. In the other case, the best limit is obtained through searches of missing energy and photons or jets (E miss T + X in [69]). For these analysis, and for the range of masses of the dark scalars that we are considering, the dark states can be effectively considered massless, and the mass splitting plays no role. The constraints in [69] are presented for a simplified model including a Dirac dark fermion, besides the Z boson. Instead, we are considering dark scalars. We can recast these bounds for our scenario using the cross-sections of the processes p p → Z → j j, p p → Z → j φ 1 φ 2 , and the analogous ones for the Dirac dark fermion 8 . With this approximate procedure, we find the exclusion limits presented as dotted lines in fig. 5. They correspond to a shift of 5% and 9% with respect to the bounds for the Dirac fermion, respectively for the first and second scenarios. Following a similar method, and assuming that future constraints on the signal cross-sections will improve as the square root of the luminosity, we estimate the reach of future searches at the High-Luminosity LHC (HL-LHC), with an integrated luminosity of 3 ab −1 . The results are shown in fig. 5 with dashed lines. Care must be taken in the interpretation of the bounds from E miss T + X. The problem is the same already discussed in sec. 2.3 for LEP and Babar: one should ensure that the φ 2 lifetime is long enough such that it decays outside the detector. To estimate the region of the parameter space where this happens we follow the same procedure outlined for LEP, requiring the proper decay length of φ 2 in the LAB frame to be larger than the detector radius, which we take equal to 10 m. This explains the vertical cuts in the excluded regions in the lower panels of fig. 5. Other searches (like prompt or displaced decays) may be relevant in the region of the parameter space where φ 2 decays inside the detector. Although interesting, they lie beyond the scope of the paper. Now we compare these searches at the LHC with the sensitivities of the CHARM, SHiP, FASER and MATHUSLA experiments. Our 90% CL sensitivity limits are presented in fig. 5. The top (bottom) panel refers to the first (second) benchmark scenario, and the left (right) panel is for δ = 10 (δ = 0.1). We find that the CHARM and FASER experiments are able to test Z masses smaller those shown in the plot, and therefore already excluded by the current bounds from LHC. The same is true for SHiP, except for the case in the bottom right panel of fig. 5. For this choice of the parameters, SHiP will be able to improve present constraints and will be competitive with future HL-LHC limits. For the same couplings, but a smaller mass splitting (bottom right panel of fig. 5), φ 2 only decays into SM leptons and φ 1 , since hadronic decays only occurs for masses such that the dark scalars can not be produced in meson decays. These leptonic decays are suppressed for our choice of couplings (g Vq = 0.1, g V l = 0.01), and it turns out that the SHiP sensitivity falls below the current LHC constraints.
The situation is drastically different for MATHUSLA. For all the four cases in fig. 5, MATHUSLA will be able to surpass current LHC limits, and even probe regions of the parameter space outside the reach of the HL-LHC. This is especially the case for the second benchmark scenario, where the couplings of the Z boson to quarks are reduced. Already in fig. 4 for the case of the EFT operators, one can notice that MATHUSLA is able to significantly extend the reach of the other experiments under consideration at large Λ. The examples studied here show that MATHUSLA is very useful to explore dark sectors, and they underline its complementary to other LHC experiments, namely CMS and ATLAS.
Let us conclude mentioning the constraints from electroweak precision measurements (EWPM). Using the results in ref. [70] one can verify that in our scenario only the following electroweak oblique parameters are generated: W , Y , X and V . The future reach of the HL-LHC on the W and Y parameters have been computed in [71] analyzing the reach on Drell-Yan processes. For a coupling g V l = 0.01 we obtain that the region that can be probed is M Z 325 GeV, not competitive with current and future direct searches.

A light dark photon
We now analyze the case of a relatively light vector boson. Specifically we focus on a dark photon, which has been extensively been studied as a portal through which the dark and the visible sector communicate. This particle interacts with the SM via the kinetic mixing operator [72][73][74]: where B µν and F µν are the field strengths of the hypercharge and the new vector boson respectively. The adimensional parameter controls the kinetic mixing, and c w is the cosine of the Weinberg angle. After a suitable field redefinition which diagonalize the kinetic and mass terms of the vector bosons, the dark photon acquires a coupling with the SM fermions. At leading order in 1 and in M Z /M Z (M Z is the mass of the SM Z boson) these couplings are [68]: where e Q f the electric charge of fermion f . As mentioned before, the dark photon has been widely studied in the literature, and relevant constraints for this scenario have been derived. We now discuss the most stringent ones focusing on the benchmark case of a dark photon mass of M Z = 40 GeV. A bound from EWPM has been computed in ref. [68]. For M Z = 40 GeV it  2.2 × 10 −2 . Additional limits come from searches for a light resonance decaying into leptons at LHCb [75] and CMS [76]. They lead to a similar bound: 2 3 × 10 −6 for M Z = 40 GeV. However, these analysis can not be applied directly to our case, since they assume that the dark photon decays into SM states only. We should instead take into account its decays into dark sector particles. This can be done computing the branching fraction for decays into leptons for the two cases where the dark photon couples to the dark sector or only to SM particles. The bound can be rescaled using the ratio of these quantities. Explicit formulas are presented in appendix B. More specifically we write: 3 × 10 −6 (4.4) and for g φ = 1 we obtain 3.3 × 10 −2 , comparable to the bound obtained from EWPM. Constraints from searches of of mono-photon events with large missing energy at LEP have been derived in [77], using the analysis of [58]. For our case the limit is comparable but less stringent than the previous ones. Future EWPM are expected to improve the sensitivity on by about a factor of 2 [68]. We use ref. [78] to derive future bounds from LHCb (see their figure 3.4.1). All these constraints and sensitivities are shown in fig. 6. There we also present the sensitivities for CHARM, SHiP, FASER and MATHUSLA, for the two representative values of the mass splitting among the dark scalars adopted in the previous sections: δ = 10 (left panel) and δ = 0.1 (right panel). These contours are derived with the same procedure outlined in sec. 4.1. For δ = 10, the current exclusion region from CHARM is comparable to the limit from EWPM, ∼ (2 ÷ 3) × 10 −2 , but for a rather limited range of φ 1 masses around m 1 ∼ 0.05 GeV. The reach of the other experiments extends well below the current excluded region, probing kinetic mixing down to 5 × 10 −5 in the case of MATHUSLA. For δ = 0.1, smaller values of can be tested, analogously to the situation presented for the EFT in sec. 3. Still large regions of the parameter space can be probed with SHiP and MATHUSLA. Once again, these results demonstrate that these future experiments are useful to search for the dark sector under consideration. A similar complementarity of bounds has been explored in the case of inelastic fermion dark matter interacting with a dark photon in [22,79,80].
We conclude this section pointing out that the dark scalar responsible for the mass of the dark photon is expected to have a mass around M Z . For simplicity we assume that this state is sufficiently decoupled and/or weakly coupled to the dark and visible sector to be ignored.

Bounds from astrophysics and cosmology
In this section we present astrophysical and cosmological bounds that can be relevant for the scenario considered in this paper.

DM abundance
Although not essential, we are assuming throughout our work that φ 1 is a stable particle. This is an attractive possibility since φ 1 could then play the role of DM. Let us briefly discuss few scenarios which could determine its cosmological abundance. Assuming that DM annihilations in the early Universe are dominated by the processes induced by the effective operator in eq. (1.7), one can compute the DM relic density after its thermal freeze-out. It turns out that, under this assumption, DM is overabundant in most of the parameter space that we are studying [57,81]. This is because, for large enough Λ, the annihilation cross-section is too small to efficiently deplete the dark matter density.
On the other hand, the dark sector might be richer of what we are assuming, and include other states into which DM can annihilate. Still, even decoupled (or very weakly coupled) dark sectors are subject to cosmological constraints. For instance, the abundance of the additional dark states should not excessively alter the expansion of the Universe during the Big Bang Nucleosynthesis (BBN) [82,83]. Moreover, cosmic microwave background (CMB) observations strongly constrain extra energy deposition into the primordial plasma due to annihilations or decays of dark states into SM particles, e.g. [84][85][86][87][88][89]. Several possibilities exist to satisfy the cosmological bounds. Examples of those are models where DM annihilates into a stable dark sector particle (which has a small cosmological abundance) [90], or where the processes determining the DM relic abundance involve dark states almost degenerate in mass with the DM, and they are suppressed at later times [91][92][93][94].
Alternatively, scenarios with an excessive DM abundance might be brought in agreement with observations if a period of entropy dilution arises after the DM freezeout, diluting the DM abundance, see e.g. [95] for light DM interacting through an heavy mediator.
Finally, let us remind that the predictions for the DM the relic abundance can be dramatically altered if the thermal history of the Universe before BBN was different from what is usually assumed. For instance, this happens if the reheating temperature of the Universe was small enough to prevent the DM thermalization with the SM plasma [96,97]. In this case, the abundance is set by the freeze-in mechanism rather than the freeze-out.
Since it is not fundamental for our discussion, in this work we do not focus on a specific model where the correct DM abundance can be obtained. It would be interesting to explore such possibility in a future work.

BBN and CMB constraints
DM with O(MeV) masses in thermal equilibrium with the SM plasma at the epoch of the neutrino decoupling (T D ≈ 2.3 MeV) are constrained by BBN and CMB observations. The reason is that these particles modify the expansion of the Universe and change the ratio between neutrino and photon temperatures, altering the production of light elements during BBN, and shifting the value of the effective number of neutrinos (N eff ), which is inferred quite precisely from CMB measurements. Current bounds for complex scalar DM exclude masses below few MeV, the precise value depending on the cosmological dataset considered, and the relative size of the DM annihilation cross-section into neutrinos and electrons or photons [98][99][100][101]. For Λ few (100) GeV and a negligible mass splitting δ, the dark scalars are decoupled from the SM bath at T D . Therefore BBN bounds do not apply in these regions of the parameter space. Notice also that the lifetime of φ 2 is shorter than one second in the regions of the parameter space probed by the CHARM, SHiP, MATHUSLA and FASER experiments, and shown in figs. 4, 5 and 6. This implies that φ 2 decays before the BBN epoch avoiding, for example, potential bounds related to the photodissociation of light elements from its decays.
As discussed before, CMB measurements also constrain DM annihilations into SM particles at early times. In our scenario φ 2 decays before the recombination era, so the relevant annihilation processes involves only pair annihilations of φ 1 into two (induced at loop level) or four SM particles. The bounds on these processes are very weak. It is also worth mentioning that φ 1 φ 2 (co)-annihilations are p-wave suppressed at low velocities. For the same reasons, indirect dark matter searches give weak limits.
As mentioned before, additional ingredients should be added to our framework in order to obtain the correct DM abundance. Depending on the concrete model, the constraints discussed in this section might apply.

Supernova bounds
Observations of supernova explosions are a precious tool to test light dark sectors weakly coupled to the SM. In the dense and hot supernova environment, with temperatures of the order of few tens of MeV, dark particles with masses up to O(0.1) GeV can be produced. If these particles interact weakly enough with the supernova material, they can efficiently escape, providing therefore an additional mechanism to cool the supernova. On the other hand, if the dark particles interacts too much, they are trapped inside the core of the supernova, and the cooling rate becomes negligible with respect to the one due to the neutrino emission. Therefore, there is an optimal range of couplings between the dark particles and the SM that can be tested with this argument, for which the dark states are abundantly produced, and they efficiently escape from the supernova. The observed cooling time of the supernova 1987A has been used to set bounds on the pair production of light DM [57,[102][103][104][105][106][107][108][109][110][111][112]. Ref. [57] considered a complex scalar coupled to electrons through an heavy Z , leading to the effective vector operator in eq. (1.7). They found that supernova limits extend up to masses of 0.2 GeV, testing Λ 1TeV at this mass scale, where the excluded region closes. Our model differs from the one in [57] in several respects. We are focusing on a democratic coupling of the dark scalars to all the SM fermions. This implies that additional processes should be consider for the production of the dark particles and their interaction with the medium 9 . Moreover, in our scenario the dark scalars present a mass splitting. This is crucial to determine whether these particles are trapped or not inside the supernova. Once produced, φ 2 decays into φ 1 and SM particles, and if the process is fast enough, only a population of φ 1 remains. The lightest scalar φ 1 can interact with the ambient particles scattering into φ 2 . However, if the mass splitting is large, this process is strongly suppressed. Therefore the bounds at large couplings (i.e. small Λ) are completely altered with respect to the degenerate case. This can be appreciated in ref. [110] for the case of inelastic fermionic dark matter. In conclusion, a dedicated analysis is needed to investigate the impact of supernova constraints in our scenario. We leave this for future work. Meantime, we shall comment that these constraints are expected to test masses up to m 1 + m 2 O(0.1) GeV. We have shown that the experiments discussed in sec. 2 attain their strongest sensitivities at larger masses. Therefore, that regions of the parameter space are left untouched by supernova constraints.

Conclusions
Let us summarize the main results of our work. We have focused on a dark sector containing a pair of non-degenerate dark scalars with masses in the MeV-GeV range. The lightest one is stable, at least on the timescales relevant for our analysis, and it could play the role of the DM. We have entertained the possibility that the dark sector communicates with the SM fermions via the effective operators in eq. (1.7). We have studied the sensitivities to this scenario of two fixed-target experiments, CHARM and the SHiP proposal, and two proposed LHC experiments, FASER and MATHUSLA. In these accelerator experiments dark scalars are produced from the collision of SM particles. If the interactions in eq. (1.7) are weak and/or the mass splitting between the scalars is small, the heaviest dark particle travels macroscopic distances before decaying and leading to a signal in a far placed detector. We have compared the projected sensitivities of these experiments with constraints from other probes, in particular searches at LEP, LHC and BaBar.
The operators in eq. (1.7) could be generated for instance by the exchange of a Z boson, coupled both to the visible and the dark sectors, see (1.4). Other possibilities exist, as mentioned in sec. 1. Motivated by the fact that the mediator can be produced on-shell at the LHC, if light enough, we have investigated this Z model. We have considered both the case of a heavy Z with a mass above the TeV and of a relatively light Z (concretely, a dark photon with a mass of 40 GeV).
Our main results are shown in figs. 4, 5 and 6. As evident, the future experiments that we have studied can significantly improve current constraints. In particular, we have shown that the sensitivity reach of the MATHUSLA experiment can surpass those from future searches at the HL-LHC. This have been demonstrated for the case of the Z model, for which an appropriate comparison with LHC constraints have been performed. The phenomenology explored here is different from that arising in models with light mediators ( 1 GeV) since, in the experiments that we have considered, the latter can be directly produced on-shell in decays of mesons, see e.g. [9,11,14,[20][21][22][23]. Therefore our analysis complement other studies for the search of dark sectors. "LINDARK," the research grant "The Dark Universe: A Synergic Multimessenger Approach No. 2017X7X85" funded by MIUR, and the project "Theoretical Astroparticle Physics (TAsP)" funded by the INFN.

A Useful equations for computations in the EFT
We collect in this appendix the expressions that we have used to compute the production and decay of dark particles. For concreteness, we will write the effective operators in terms of the vector and axial couplings defined in eq. (1.9). To compute the interactions of light mesons with a pair of dark particles we follow ref. [56,113], including the dark current in the Chiral Perturbation Theory Lagrangian. We also compute the Wess-Zumino-Term following ref. [114] to consider interactions involving photons. The relevant interactions are given by MeV is the pion decay constant, F αβ is the photon field strength and J φ is the dark current defined in eq. (1.5). The first sum is taken over the pseudoscalar mesons P = {π 0 , η 1 , η 8 } 10 , while the second sum is taken over the conjugate pairs P = {π + , K + , K 0 } andP = π − , K − ,K 0 , with coefficients and Notice that the pseudoscalar mesons η 8 and η 1 appearing in eq. (A.1) are not the mass eigenstates. To rotate to the mass basis we will follow Ref. [115] and write where [115] f The couplingsc P refer to the case where the dark current couples to an axial fermion current (see eq. 1.8) while we focus on vector interactions. They are reported here only for sake of completeness. For the vector mesons we use instead the matrix elements that can be found in the literature. More specifically for the heavy J/Ψ and Υ we follow ref. [116] and write where M is the vector meson mass. All other matrix elements vanish. Introducing the Wilson coefficients of eq. (1.9), we define f V = c Vq f q V , where q is quark flavor correspondent to the meson M. Since f V is connected to the quarkonium wave function and it is difficult to compute it from first principles, we will express all our results in terms of the (known) decay widths into electrons, as explained later. For the light vector mesons (ρ and ω) we instead use the results of ref. [117]. Considering the quark composition of the mesons wave functions, we write f V = q c Vq f q V , with Notice that for the democratic choice of the Wilson coefficients in fig. 1 (c V f = 1) the effective coupling of the ρ meson is suppressed with respect to the analogous one for the ω meson. We are now in the position to give explicit formulas for the decay widths.

A.1 Production of dark particles via mesons decays
Under the assumption c f R = c f L the relevant production channels are (i) V → φ 1 φ 2 with V = {ρ, ω, J/Ψ, Υ} the vector mesons, and (ii) P → φ 1 φ 2 γ with P = {π 0 , η, η } the pseudoscalar mesons. Denoting by M the meson mass and defining A(x, y, z) = 1 + (x − y) 2 z 2 − 2 x + y z , (A.10) the decay widths are given by In the second equation e is the electric charge, the couplings c P are defined in eq. (A.7), and the total decay width is obtained integrating from (m 1 + m 2 ) 2 and M 2 . As for the first equation, in the case of the ρ and ω mesons we use eq. (A.9) to express f V in terms of the Wilson coefficients. In the case of the heavy quarkonia vectors J/Ψ and Υ we instead use Γ(V → e + e − ) = f q 2 V e 4 Q 2 q 12 π M (A. 12) to express f q V in terms of the known decay width into an electron-positron pair. The electric charge of the relevant quark is denoted by Q q , respectively 2/3 and -1/3 for the J/Ψ and Υ mesons.

A.2 Decays of φ 2
We collect now the equations used to compute the decay widths of φ 2 . We define the auxiliary functions The total decay width for the φ 2 → φ 1f f decay is obtained integrating s between (m 1 + m f ) 2 and (m 2 − m f ) 2 , while for the φ 2 → φ 1 γP and φ 2 → φ 1 PP decays the integration extends between (M + m 1 ) 2 and m 2 2 , and between (m 1 + M ) 2 and (m 2 − M ) 2 , respectively. The couplings c P and c P are defined in eq. (A.7) and (A.4), while e is the electric charge. For the decay φ 2 → φ 1f f we show both the contributions generated by c V f and c A f . For all fermions apart the neutrinos we always consider c A f = 0, leaving only the vector contribution. For the neutrinos, on the contrary, we are considering a V-A current, i.e. c Aν = −c Vν . Notice that the decays φ 2 → φ 1 PP are not shown in fig. 1 because the effective couplings involved in these decays vanish for the choice of the Wilson coefficients c V f = 1 adopted in that plot.
When the momentum transfer in the decay is sufficiently large, chiral perturbation theory breaks down. Therefore we adopt the following prescription: for M 2 − M 1 > 2 GeV the hadronic decays of φ 2 are computed using perturbative QCD (we consider the processes φ 2 → φ 1 ff , described by the third eq. A.14 multiplied by N c = 3), while in the opposite regime we use the description in terms of mesons presented above. Notice however that for momentum transfers 0.5 − 2 GeV both approaches are not appropriate. In this window the hadronic decays are difficult to model. An approach in terms of form factors have been pursued in [81,118].