Sensitivity of the intensity frontier experiments for neutrino and scalar portals: analytic estimates

In recent years, a number of intensity-frontier experiments have been proposed to search for feebly interacting particles with a mass in the GeV range. We show analytically how the characteristic shape of the sensitivity regions of such experiments - upper and lower boundaries of the probed region, the maximal mass reach - depends on the parameters of the experiments, taking the SHiP and the MATHUSLA experiments as an example. We find a good agreement of our estimates with the results of the Monte Carlo simulations.

1 Introduction: searching for feebly coupled particles The construction of the Standard Model has culminated with the confirmation of one of its most important predictions -the discovery of the Higgs boson. The quest for new particles has not ended, however. The observed but unexplained phenomena in particle physics and cosmology (such as neutrino masses and oscillations, dark matter, baryon asymmetry of the Universe) indicate that other particles exist in the Universe. It is possible that these particles evaded detection so far because they are too heavy to be created at accelerators. Alternatively, some of the hypothetical particles can be sufficiently light (lighter than the Higgs or W boson), but interact very weakly with the Standard Model sector (we will use the term feeble interaction to distinguish this from the weak interaction of the Standard Model). In order to explore this latter possibility, the particle physics community is turning its attention to the so-called Intensity Frontier experiments. Such experiments aim to create high-intensity particle beams and use large detectors to search for rare interactions of feebly interacting hypothetical particles. New particles with masses much lighter than the electroweak scale may be directly responsible for some of the BSM phenomena, or can serve as mediators (or "portals"), coupling to states in the "hidden sectors" and at the same time interacting with the Standard Model particles. Such portals can be renormalizable (mass dimension ≤ 4) or be realized as higher-dimensional operators suppressed by the dimensional couplings Λ −n , with Λ being the new energy scale of the hidden sector. In the Standard Model there can only be three renormalizable portals: -a scalar portal that couples gauge singlet scalar to the H † H term constructed from a Higgs doublet field H a , a = 1, 2; -a neutrino portal that couples new gauge singlet fermion to the abLa H b where L a is the SU(2) lepton doublet and ab is completely antisymmetric tensor in two dimensions; -a vector portal that couples the field strength of a new U(1) field to the U (1) hyperfield field strength.
Let us denote a new particle by X. The interaction of X with the SM is controlled by the mixing angle θ X -a dimensionless parameter that specifies the mixing between X and the corresponding SM particle: the SM neutrinos for the neutrino portal, the Higgs boson for the scalar portal and the hyperfield for the vector portal. The searches for such particles are included in the scientific programs of many existing experiments [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Although the LHC is a flagship of the Energy Frontier exploration, its high luminosity (especially in Run 3 and beyond) means that huge numbers of heavy flavored mesons and vector bosons are created. This opens the possibility to supplement the High Luminosity phase of the LHC with Intensity Frontier experiments associated with the existing interaction points. Several such experiments have been proposed: CODEX-b [16], MATHUSLA [17,18], FASER [19], and AL3X [20].
Given that all these experiments can probe similar parameter spaces, it is important to be able to assess their scientific reach in a consistent way, under clearly specified identical assumptions. While Monte Carlo simulations of both production and decay, complemented with background studies, will eventually offer the ultimate sensitivity curve for each of the experiments, it is crucial first have a sufficiently simple and fully controlled analytic estimate. Such an estimate shows the main factors that influence the sensitivity, which permits scanning over different models and evaluate different experimental designs in a quick and effective way.
As it turns out, to a large extent, the ratio between the sensitivities of the experiments to a great extent does not depend on the specific model of new physics, but is mainly determined by the geometry and the collision energies of the experiments, which allow for a comparison in a largely model-independent way. To illustrate this point in this paper, we compare the potentials of two proposed experiments: the LHC-based MATHUSLA experiment [17,18,[21][22][23] and a proton fixed target experiment using the proton beam of the Super Proton Synchrotron (SPS) at CERN -SHiP [24][25][26]. We analyze their sensitivity to the neutrino [27][28][29][30][31][32] and scalar [33][34][35][36][37] portals. For particle masses M X m Bc 1 the main production channels are decays of heavy flavored mesons and W bosons [22,25] (see also Appendix B.1 for a brief overview). We concentrate on the mass range M X m K , since the domain of lower masses for the HNL and Higgs-like scalar is expected to be probed by the currently running NA62 experiment [14,38].
The sensitivity of the experiments is determined by the number of events that one expects to detect for a set of given parameters. In realistic experiments such events should be disentangled from the "background" signals.
For SHiP, detailed simulations have shown that the number of background events is expected to be very low, so that the experiment is "background free" [24,[39][40][41]. For MATHUSLA, the background is also expected to be low [18,21], although no simulation studies of background have been performed. Even in the most favorable case of N bg 1 one needs on averageN events = 2.3 expected signal events to observe at least one event with the probability higher than 90%. 2 However, due to the lack of spectrometer, mass reconstruction and particle identification at MATHUSLA, the meaning of the discovery of 2.3 events in the two experiments is very different as there is no way to associate the signal with a model in MATHUSLA and further consolidate the discover.
For both experiments considered here the production point ("target") is separated from the detector decay volume (of length l det ) by some macroscopic distance l target-det (see Appendix D). For such experiments the sensitivity curve has a typical "cigar-like shape" in the plane "mass vs. interaction strength". The upper boundary is determined by the condition that the decay length l decay of a produced particle becomes comparable to the distance to the decay volume, l decay ∼ l target-det . The lower boundary of the sensitivity region is determined by the parameters at which decays become too rare (i.e. l decay l det ). The intersection of these bounds defines the maximal mass which can be probed at the experiment.
The number of decay events in the decay volume factorizes into 1) where N prod,M is the number of particles X that are produced from a mother particle M and P decay,M is the decay probability. For N prod,M we have Here, N M is the number of mother particles produced at the experiment; in the case of mesons N M = N meson = 2Nq q ×f meson , where f meson is the fragmentation fraction of quark q into a given meson, and N M = N W in the case of the W bosons. BR M →X is the total branching ratio of decay of the mother particle into X (see Appendix B.1). Finally, decay is the decay acceptance -the fraction of particles X whose trajectory intersects the decay volume, so that they could decay inside it.
The probability of decay into a state that can be detected is given by 3 where the branching ratio BR vis is the fraction of all decays producing final states that can be registered. Finally, det ≤ 1 is the detection efficiency -a fraction of all decays inside the decay volume for which the decay products could be detected. In the absence of detector simulations we optimistically assume a detector efficiency of MATHUSLA of det = 1. The decay length l decay in Eq. (1.3) is defined as where τ X is the lifetime of the particle X (see Appendix B.2), β X is its velocity and γ X is the γ factor (which depends on the mother particle that produces X).
The production branching ratio and the lifetime behave with the mixing angle as At the lower bound of the sensitivity the decay probability behaves as P decay ≈ l det /l decay , and as a consequence of (1.5) the number of events scales as At the upper bound P decay ≈ e −l target-det /l decay , and where C is some numerical factor (that depends on properties of X). Larger γ factor suppresses the exponents in the expression for the decay probability (1.3). From (1.6), (1.7) we see that this affects the upper and lower bounds of the cigar-like sensitivity plots in the opposite ways. For the lower bound, an experiment with the smaller average γ factor is sensitive to small coupling constants. For sufficiently large couplings, larger γ factor ensures that particles do not decay before reaching the detector, thus increasing the sensitivity to the upper range of the sensitivity curve.
The paper is organized as follows. In Sections 2-4 we discuss the lower and upper boundaries of the sensitivity region, the maximal mass that can be probed and experimental parameters that affect them. In Sec. 5 we discuss the total amount and energy distribution of charm-and beauty mesons at both SHiP and MATHUSLA experiments, as well as the contribution from the W bosons. In Sec. 6 we summarize and discuss our results. Finally, in Sec. 8 we make conclusions. Appendices A-H provide details of computations and relevant supplementary information.
2 Lower boundary of the sensitivity region: main factors As we will see later (see Section 6), the production from the W bosons does not give a contribution to the lower bound of the sensitivity curve for neither of the two experiments, and for neither of the two models discussed. So, in this Section we will consider only the production from the mesons.
Let us first estimate the lower boundary of the sensitivity region, where l decay l det . For the number of events (1.1) we have where X ≡ prod × decay × BR vis is the overall efficiency and τ X is the lifetime of the particle X (see the discussion below Eq. (1.4)). The particles are assumed to be relativistic (we will see below when this assumption is justified), so that β X ≈ 1. We estimate the average γ factor γ X using the average γ factor of the mother meson as where E rest X is the average energy of X in the rest frame of the meson (see Appendix C).
Since at the lower bound N events ∝ θ 4 X (see Eq. (1.6)), for the ratio of the mixing angles at the lower bound, we have where we assumed that the same meson is the main production channel at both the SHiP and MATHUSLA experiments for the given mass M X of the new particle, so the branching ratio BR meson→X from Eq. (2.1) disappears. Therefore, to make a comparison between the experiments we only need to know the total number of mesons, their average γ factor, the decay volume length and the overall efficiency.

Upper boundary of the sensitivity curve
For sufficiently large interaction strength (i.e., the mixing angles) particles decay before reaching the decay volume. This determines the upper bound of the sensitivity curve, that we call θ 2 X,upper . A useful quantity to consider is a mixing angle for which the amount of decays inside the decay volume is maximal, θ X,max . It can be found using the asymptotic behavior for the number of events N events from the estimations (1.6), (1.7). In the domain l decay l target-det for a fixed mass M X it follows that N events monotonically grows as θ 4 X with the increase of θ X , while in the domain l decay l target-det it falls exponentially. Therefore, the position of the maximum θ X,max can be estimated by requiring that the decay length is equal to the distance from the target to the decay volume: Using θ X,max , we can estimate the value of θ X,upper assuming that all the particles X have the same energy equal to the average energy E X . If we neglect the second exponent in the expression for the decay probability (1.3), then the formula for the number of events (1.1) becomes We can estimate the exponent in (3.2) as l target-det /l decay ≈ θ 2 X /θ 2 X,max , see Eq. (3.1). So imposing the condition N events 1 in Eq. (3.2) with the logarithmic precision we get An example of the dependence of the number of events on θ 2 X for the fixed mass M X , together with the estimation of the θ X for the maximal number of events given by (3.1) and the upper bound predicted by (3.3), is shown in Fig. 1.
However, it is not sufficient to use only the average energy E X to estimate the position of the upper boundary. Indeed, the decrease of cτ X with the growth of θ 2 X can be compensated by the increase of the energy E X and, therefore, of the γ-factor. As a result the particles with E X > E X can reach the detector even if the mixing angle θ X is larger than the estimate (3.3). Consider the energy distribution of the X particles, Taking into account this distribution, the formula for the decay probability (1.3) at the upper bound should be modified as where the function π(y), defined via determines a "window" of energies in which the shape of f X (E X ) distribution (rather than the averange number of particles) contributes to the overall probability. π(y) is shown in Fig. 2. For small energies (small y) π(y) is exponentially small, while for large energies (large y) π(y) is inversely proportional to energy and decreases slowly. Therefore, a sufficiently long "tail" of high-energy mesons can contribute to the integral in (3.5), but this range cannot be estimated without knowledge of the distribution function f X .

Maximal mass probed
The maximal mass probed by the experiment is defined as the mass at which the lower sensitivity bound meets the upper sensitivity bound. It can be estimated from the condition that the decay length, calculated at the lower bound θ X,lower (see Sec. 2) is equal to the distance from the target to the decay volume of the given experiment: l decay (M X,max , θ 2 X,lower (M X,max )) l target−det . , where the term α in the exponent approximates the behaviour of the lifetime with the mass, and the term 1 comes from the γ factor.
Using the condition (4.1), the maximal mass probed can be estimated as and the ratio of the maximal mass probedes at SHiP and MATHUSLA experiments is For the scalar we have α ≈ 2, while for the HNLs it is α ≈ 5, see Appendix B.2. The estimate of the maximal mass probed (4.2) is applicable only if the result does not exceed the kinematic threshold; for the production from B mesons for the HNLs it is m Bc −m l or m B −m l depending on whether amount of produced B c mesons is large enough to be relevant for the production (see the discussion in Sec. 5.1), and for the scalars it is m B − m π .

Number and momentum distribution of mesons and W 's at SHiP and MATHUSLA
In this Section, we discuss the number and distribution of charm and beauty mesons and of W bosons at SHiP and MATHUSLA experiments. As we have seen, to estimate the lower boundary we need only the number of parent particles and their average γ factors (see Eqs. (2.1), (4.3)). On the other hand, for the estimation of the upper boundary we also need the energy distribution of the mesons and W (see Sec. 3).

B and D mesons
The main production channel of the HNLs in the mass range M N m Ds is the two-body leptonic decay of D s mesons. For the masses m Ds M N m Bc the main contribution comes from the decays of B mesons, see, e.g., [42]. 4 For masses M N 3 GeV the main HNL production channel depends on the value of the fragmentation fraction of B c mesons, f Bc : if f Bc 10 −4 then it is the two-body decay of the B c meson, while for smaller values it is the two-body decay of the B + meson [43]. For the scalars the production from D mesons is negligible as compared to the B +/0 mesons decays even for masses m K M S m D . The B c mesons are not relevant for their production (see, e.g., [44,45]). The branching ratios of the production of the HNLs and the scalars used for our estimations are given in Appendix B.1.
For the LHC energies the fragmentation fraction f Bc was measured at the LHCb [46] and found to be f Bc ≈ (2.6 ± 1.3) × 10 −3 . Earlier measurements at the Tevatron give a similar value f Bc ≈ (2 ± 1) × 10 −3 [47][48][49], which is in good agreement with [46]. Therefore, at the LHC the B c decay is the main production channel for heavy HNLs. However, at the energies of the SHiP experiment, √ s 30 GeV, currently there is no experimental data on f Bc . Additionally, the theoretical predictions of f Bc (see, e.g., [50][51][52]) disagree with the LHC and Tevatron measurements at least by an order of magnitude, which also makes them untrustable at SHiP's energies. As a result, the value of f Bc at SHiP experiment is unknown. In order to estimate the effect of this uncertainty, we perform our analysis of the sensitivity of the SHiP experiment for two extreme cases: (i) SHiP's f Bc at the same level as at the LHC, and (ii) f Bc = 0, "no B c mesons".
Let us now discuss the available data. For the SHiP experiment, the amounts of produced charmed and beauty mesons (except the B c mesons) were obtained in detailed PYTHIA simulations; the corresponding numbers can be found in [53] and are reproduced in Table 1. We estimate the spectrum of the B c mesons from the spectrum of the B + mesons by rescaling the energy E Bc = (m Bc /m B )E B for the events with B + mesons. For MATHUSLA experiment, the situation is different: there is no available data with detailed simulations that give us the relevant properties of the mesons, so we discuss them below.

Mesons at MATHUSLA
In order to estimate the number of mesons and their γ factors for the MATHUSLA experiment, one needs to know their p T distribution at ATLAS/CMS in the MATH-USLA pseudorapidity range 0.9 < η < 1.6 (see Appendix D). The relevant distributions were measured for B + mesons by the CMS collaboration [54] (13 TeV) with the p T cut p B T 10 GeV, and for D + /D 0 mesons by the ATLAS collaboration [55] (7 TeV) for p D T 3.5 GeV. We show the spectra obtained in these papers in Fig. 4. The low p T mesons, unaccounted for these studies, are the most relevant for the MATHUSLA sensitivity estimate because of two reasons. Firstly, the p T spectrum of the hadrons produced in pp collisions has a maximum at p T ∼ few GeV (see, e.g., experimental papers [56,57], theoretical paper [58] and references therein), and therefore we expect that most of the D or B mesons have p T s below the LHC cuts. Secondly, low p T mesons produce decay products with the smallest γ factor, and therefore with the shortest decay length (1.4) and the largest probability to decay inside the decay volume (here we consider the case l decay l target-det ). Therefore, by shifting the position of the peak to smaller p T s, we increase the number of mesons and decrease their average γ factor, and both of these effects enhance the number of   Figure 3. The prediction of FONLL for the ratio of the cross-sections of the production of the ψ(2S) and J/ψ mesons from B mesons decays as a function of their p T , compared to experimental data from CDF at the Tevatron and CMS and LHCb at the LHC. The figure is reproduced from [58].
events at the lower bound (2.1). Therefore an accurate prediction of the distribution dσ/dp T in the domain of low p T s is very important. In order to evaluate the distribution of heavy flavored mesons at low p T s (and also to estimate the D meson production cross-section at √ s = 13 TeV) we use FONLL [58][59][60][61]. The predictions of FONLL have been calibrated against the accelerator data and were found to be in very good agreement, see e.g. [54,55,[62][63][64]. In particular, by comparing the simulations of FONLL with the measurements at the Tevatron and LHC of the production of the B + s, we see that FONLL predicts the low p T distribution accurately, see Fig. 3. We show the central values of the FONLL predictions down to p T = 0, confronted with the measurements of the CMS [54] and ATLAS [55] collaborations in Fig. 4. As expected, the distributions have maxima, after which they fall. We see, however, that the central predictions of FONLL for the differential cross-sections typically lie below the uncertainty range of the experimental cross-section, which results in a somewhat lower total cross-sections. Namely, integrating the central predictions over the experimentally measured p T s, we have σ D,FONLL /σ D,exp ≈ 0.4 and σ B,FONLL /σ B,exp ≈ 0.7. However, as is demonstrated in the same papers [54,55], the agreement between the FONLL predictions and the experimental data is much better if one uses the upper bound of the FONLL predictions defined by the theoretical uncertainties.
Using the results of the FONLL simulations, we find the amounts of low p T mesons traveling in the MATHUSLA direction:    This justify our statement that most of the B and D mesons have the p T below the cuts in the currently available experimental papers [54,55]. FONLL does not provide the distributions of the D s and the B c mesons. We approximate their distributions by the D + and B + distributions. In the case of the B c mesons we justify this approximation by comparing the distributions provided by BCVEGPY 2.0 package [65] (that simulates the distribution of the B c mesons and was tested at the LHC energies) for the B c meson with that of FONLL for the B + meson. We conclude that the p T and η distributions of B c and B + have similar shapes.
The relevant parameters -the total number of mesons, the average γ factor of the mesons that are produced in the direction of the decay volume of the experiments and the geometric acceptances geom,meson for the mesons -are given in Table 1.

W bosons
The production channel from the decays of W bosons is only relevant for the MATH-USLA experiment since the center of mass energy at SHiP experiment is not enough to produce on-shell W bosons. The total W boson production cross-section at the LHC energies √ s = 13 TeV was measured in [66] as σ W →N +l ≈ 20.5 nb. The corresponding number of W bosons produced during the high luminosity phase of the LHC is The p T distribution of the W bosons at the LHC in the pseudorapidity range |η| < 2.5 and for energies √ s = 7 − 8 TeV was measured by the ATLAS and CMS collaborations [67,68]. The p T spectrum for √ s = 8 TeV based on the data [68] is shown in Fig. 5. Their results show that most of the vector bosons are produced with low p T s. However, owing to the lack of precise knowledge of the W rapidity distribution, these spectra do not give us the magnitude of the W 's average momentum p W , needed to estimate the decay acceptance and the average momentum of HNLs.
In order to obtain p W we have simulated the process p + p → W ± in Mad-Graph5 [69]. In the leading order we have obtained σ W →ν+l ≈ 15.7 nb, which is in reasonable agreement with the prediction [66]. The resulting momentum distribution of W bosons is shown in Fig. 6 (left). A remark is in order here: at the leading order MadGraph5 does not predict the p T distribution of W s, since it is 2 → 1 process; therefore, all of the W bosons in the simulations fly along the beam line, and the magnitude of their momentum is given by the longitudinal momentum p L . The realistic p T spectrum can only be obtained after implementation of the parton showering. However, based on the above-mentioned measurements [67,68], the typical p T 's of W bosons are significantly smaller than their typical p L and therefore we chose to neglect the p T momentum of the W bosons in what follows.
Having the W boson distribution dN W /dp W , we can obtain the decay,W and the average HNL momentum p X by calculating the distribution of the particles in the 10 -3 energy E X and the angle θ X between the direction of motion of the X and the beam: Here d 2 BR W →X /dθ X dE X is the differential production branching ratio, and P (θ X ) is a projector which takes the unit value if θ X lies inside MATHUSLA's polar angle range and zero otherwise. Let us compare the amounts of the X particles produced from the W bosons and from B mesons and flying in the direction of the decay volume. We have where we used the amount of B mesons at the LHC and the decay acceptance from the Table 2, the number of W s at the LHC (5.2) and the branching ratios of the scalar and HNL production from Appendix B.1. Therefore we conclude that for the scalars the production from the W s is not relevant, while for the HNLs careful estimation is needed.
In the case of the HNL the differential branching ratio in the Eq. (5.3) is The energy and angular distributions of the HNLs from the W bosons at MATH-USLA are almost independent of the HNL mass in the mass range of interest, M N m W . It is an expected result because the kinematic in this limit should not depend on small HNL masses. The energy distribution for M N = 1 GeV is shown in Fig. 6. The decay acceptance was found to be decay,W 2%, while the average momentum of the produced HNLs is p N ≈ 62 GeV.
The shape of the energy spectrum of the HNLs can be qualitatively understood in the following way. For a given value of the angle θ N of the HNL, the energy distribution has a maximum at E N,max (θ N ) = m W /2 sin(θ N ), 5 which corresponds to HNLs produced from the W mesons with some momentum p W,max (θ). As a consequence, the biggest amount of HNLs flying in direction θ N has energy close to their maximum, see the dashed line at the right panel of Fig. 6. The total energy spectrum is a superposition of different angles and has a peak at E N,peak ≈ 58 GeV which corresponds to the maximal angle possible at MATHUSLA, θ max MAT ≈ 44 • . From the other side, the maximal energy possible for HNLs at MATHUSLA is defined by the minimal angle θ min MAT ≈ 22 • , which explains why the spectrum tends to zero near the energy E N,max ≈ 106 GeV.

.1 Efficiencies
Using the results of Sec. 5, we have almost all ingredients needed to estimate the lower bound, the upper bound, the maximal mass probed and the total sensitivity curve. The only questions remaining are the following. The first one is the relation between the mesons spectra and the X particles spectra. The second one is the value of the overall efficiency = decay × det × BR vis , (6.1) where the quantities decay , det , BR vis are the decay acceptance, detection efficiency and the visible branching correspondingly; they are defined by Eqs. (1.2), (1.3).
We approximate the spectra of the X particles originating from the mesons and flying to the decay volume by the distributions of the mesons flying in the direction of the decay volume. To take into account the kinematics of the meson decays, we use the relation (2.2) between γ factors of the X particle and the meson in the expressions (1.3), (3.5) for the decay probability.
Let us discuss the efficiencies. For the HNLs at SHiP experiment, we used the values of decay and det provided by detailed FairSHiP simulations [70]. The results of the SHiP collaboration on the sensitivity to the scalars are not currently available, and for the product of decay · det we used the value for the HNL averaged over its mass, decay · det ≈ 0.2.
For the MATHUSLA experiment there currently is no such detailed analysis of the efficiencies and background. In [17,21] it is claimed that all the SM background can be rejected with high efficiency, but detailed simulations are needed for the justification of this statement. Here we optimistically use det = 1. For the decay acceptance of the particles produced from the mesons we use the geometric acceptance of the mesons at MATHUSLA, which we obtained using FONLL. 6 For the decay acceptance of the HNLs produced in the decays of the W bosons we used the value decay,W ≈ 0.02 obtained in Sec. 5.3. All the parameters above, together with geometrical properties of the experiments are summarized in Table 2. We estimate l det and l target-det using an assumption that the angular distribution of the X particles in the angular range of the decay volume is isotropic, see Appendix D for details.
The last parameter needed parameter is the visible decay branching fraction. Following [21,43], for the visible decay branching fractions for both MATHUSLA and SHiP experiments we include only the decay channels of the X particle that contain at least two charged tracks. Our estimation of BR vis is described in Appendix B.2.3. The plots of the visible branching ratios for the HNLs and for the scalars are shown in Fig. 7.

Lower bound
Let us first compare the relevant parameters of the experiments summarized in Tables 1, 2. One sees that the effective number of D mesons is approximately two orders of magnitude larger at SHiP, 7 the effective numbers of B mesons are comparable be- Qualitatively, for the particles produced in the decays of the B mesons (HNLs with masses M N > m D and scalars with masses M S > m K ) MATHUSLA can probe mixing angles a factor 5 smaller than SHiP due to the smaller γ factor of the B mesons and the larger overall efficiency. For the HNLs in the mass range m K M N m D the smallness of γ factor of the D mesons at MATHUSLA and the suppression of the number of events at SHiP by the overall efficiency cannot compensate the difference of two orders of magnitude in the effective numbers of the D mesons, and therefore the SHiP reaches a sensitivity which is about half an order of magnitude lower in U 2 . We note again that the result (6.2) was obtained under the optimistic condition det = 1 for MATHUSLA; after using a realistic efficiency the lower bound of the sensitivity at MATHUSLA will be changed by a factor 1/ √ det , which will affect the ratio (6.2).

Upper bound
We show the dependence of the number of events at θ 2 X = θ 2 max as a function of the mass for the HNLs mixing with ν e and the scalars in Fig. 8. We see that by the maximal number of events the SHiP experiment is much better than the MATHUSLA experiment, which is explained by the shorter length to the decay volume and higher value of the average gamma factor.
With the energy distributions of the mesons and the W bosons obtained in Sec. 5, let us now estimate their effect on the upper bound of the sensitivity. To do this, we N events (θ max events ) Figure 8. The dependence of the number of events at SHiP and MATHUSLA evaluated at U 2 = θ 2 max for the HNLs mixing with ν e (left) and for scalars (right). Dashed lines denote the values for U 2 max for which the sensitivity of SHiP and MATHUSLA intersects the domain that has been closed by previous experiments (see, e.g., [25]).
introduce the width of the upper bound defined by We take the HNLs as an example, commenting later on the difference with the scalar. We will be interested in the HNLs with M N 2 GeV (for smaller masses θ 2 upper lies deep inside the region excluded by the previous searches, see e.g. [42]). The HNLs in question are produced from the decays of B mesons and W bosons. Our procedure of the estimation of the upper bound width is based on (3.5). As we already mentioned at the beginning of this Section, in the case of the production from B mesons we approximate the spectra of the HNLs by the spectra of the B mesons (so that the HNLs fly in the same direction as the B mesons) and take into the account the relation (2.2) between the B meson and the energies of HNLs. In the case of the production from the W bosons, we use the energy spectrum of the HNLs from Fig. 6. We approximate the shapes of the high-energy tails of these spectra by simple analytic functions. For the B mesons at SHiP, the fit is an exponential function, for the B mesons at MATHUSLA the fit is a power law function, while for the HNLs from the W bosons the fit is a linear function, see Appendix E.1. Using the fits, we calculate the upper bound θ 2 X,upper using the steepest descent method for the evaluation of the integral (3.5). The derivation of θ 2 X,upper is given in Appendix E.2. Using θ 2 X,upper , we present the upper bound width (6.3) in Fig. 9. We also show there the prediction of the estimations of the upper bound width which assume that all of the produced particles have the same energy, see Eq. (3.3).
We see that for the particles from B mesons at SHiP and for the HNLs from the W bosons at MATHUSLA the broadening of the width due to the distribution is small, while for the particles from B mesons the distribution contributes significantly. This is a direct consequence of the behavior of the shape of the high-energy tails of the distributions. Namely, for the B mesons at SHiP, the number of high-energy HNLs is exponentially suppressed. For the HNLs originating from the W bosons the tail falls linearly, and naively the upper bound would be significantly improved. However, the distribution becomes zero not very far from p N , and the effect of the contribution is insignificant. Only for the B mesons at MATHUSLA the tail causes significant improvement of the width of the upper bound.
Finally, let us comment on the difference between the shapes of the width between the HNL and scalar cases. The lifetime τ S is changed with the mass slower than τ N , see the discussion in Sec. B.2. In addition, Br B→S behaves with the mass monotonically, while for HNLs new production channels appear at different masses. Therefore the upper bound of the sensitivity region for the scalars changes less steeply and more smoothly with their mass, see Fig. 9 (right panel).
The comparison of the upper bound of the sensitivity for the HNLs originating from W bosons and B mesons is shown in Fig. 10. Our method of obtaining the sensitivity is summarized in Appendix F. We see that the W s determine the upper bound. The reason for this is that the HNLs from W s have sufficiently larger average momentum, which compensates the production suppression (see Eq. (5.4)).

Maximal mass probed
The smaller γ factor of the mesons at MATHUSLA adversely affects the upper bound of the sensitivity curve and thus the maximal mass probed. In particular, for the HNLs mixing with ν e/µ , the maximal mass probed ratio SHiP experiment based on the definition above exceeds the kinematic threshold, and therefore the result (4.3) is not valid. However, for the scalars the maximal mass for the MATHUSLA experiment is smaller than the kinematic threshold, which is still a consequence of smaller γ factor.

Comparison with simulations
Next we compare our sensitivity estimates with the results of the SHiP [43,71,72] and MATHUSLA [21] collaborations. Our method of obtaining the sensitivity curves is summarized in Appendix F.

HNLs
The results for the HNLs are shown in Fig. 11. To facilitate the cross-check of our results, we also provide simple analytic estimates of the lower boundary for several HNL masses (see Appendix G). Small discrepancies between the simple estimation of the lower bound and numeric result are caused by the difference in the values of 1/ p meson and 1/p meson , which actually defines the lower bound.
For the sensitivity of the SHiP experiment, there is good agreement of the sensitivity curves, with a slight difference in the maximal mass probed. We think that this is due to the difference in the average γ factors used in our estimation and those obtained in Monte Carlo simulations by the SHiP collaboration. Indeed, using the SHiP simulations results available in [70], we have found that for the masses  Figure 11. Comparison of the sensitivities to the HNLs mixing with the electron flavor obtained in this paper (solid lines) with the results of SHiP [43] and MATHUSLA [21] (dotted lines). For the MATHUSLA experiment, the contributions from both B, D mesons and from W bosons are shown separately. For the SHiP experiment, we consider the case of maximally possible contribution of B c mesons, given by the fragmentation fraction f Bc = 2.6 · 10 −3 measured at LHC energies √ s = 13 TeV [46].
M N m Bc the ratio of average γ factors is γ analytic N / γ simulations N 0.8, which seems to explain the difference.
For the sensitivity of MATHUSLA to the HNLs from W s there is good agreement for the entire mass range probed. For the sensitivity to the HNLs from B and D mesons, the situation is somewhat different. In the mass range M N m Ds , where the main production channel is the decay of the B mesons, there is reasonable agreement with our estimate. The discrepancy can be caused by higher average energy of the HNLs in the simulations, which simultaneously lifts up the lower and upper bounds of the sensitivity. The reason for the difference at masses M N < 1 GeV is not known, see a discussion in Appendix H.

Scalars
The comparison of our sensitivity estimates with the results of the SHiP and MATH-USLA experiments is presented in Fig. 12. We also show the results of a simple analytic estimate of the lower bound for particular masses from Appendix G. Our model of the scalar phenomenology described in Appendix B.1.2, B.2.2 differs from the models used in [21]. Namely, for obtaining the scalar production ratio we have summed over exclusive channels B → X s/d + S [45], 8 while in [21] the estimation of 8 There is 30% level uncertainty in the total production rate because of the B → K * 0 (700) + S channel. The meson K * 0 (700) is not observed experimentally and it could be either a di-quark or a tetra-quark state (see, e.g. [74] and references therein). We did an estimation assuming that K * 0 (700) is the di-quark state, in the other case this production channel is absent.  [71,72] and MATHUSLA [21]. Note that in order to facilitate comparison with collaboration results we have used different scalar production and decay models in left and right panels: for comparison with MATH-USLA results we took the model from [21], while for SHiP we used the model from [73], see the text for details. Orange points, based on analytical estimates of the lower boundary, allow for simple cross-check of our results, see Appendix G for details.
the production from B based on the free quark model was used, which causes differences in the magnitude of the branching ratio and kinematic production threshold. We note that the free quark model breaks down for large scalar masses M S 3 GeV since the QCD enters the non-perturbative regime, and therefore it gives meaningless predictions for the production rate. As for the scalar decay width, because of theoretical uncertainty for the mass range 2m π M S 2m D there is no agreement in the literature how to describe the scalar decays in this domain, see [75]. Our decay width differs significantly from the decay width used in [21]. In the case of SHiP [72], the decay widths used in the estimates are the same, but for the production channel was used the results from [73]. Therefore, in order to make a meaningful comparison we used the scalar production branching ratio and the decay width from [21] for the MATHUSLA experiment and the production branching ratio from [73] for the SHiP experiment.
The sensitivity curves are in good agreement. For MATHUSLA, a small difference in the position of the maximal mass probed can be explained by different energy distributions of the scalars used in our estimate and in [21] and in [72].

Conclusions
In this work, we investigated the sensitivity of Intensity Frontier experiments to two models of new super-weakly interacting physics: heavy neutral leptons and dark   Figure 13. Comparison of the sensitivity of SHiP and MATHUSLA for the HNL. The production fraction of B c mesons at SHiP energies √ s ≈ 28 GeV is not known, and the largest possible contribution is based on the production fraction measured at the LHC, f (b → B c ) = 2.6 × 10 −3 . In the case of the SHiP experiment we used the overall efficiency calibrated against the Monte Carlo simulations [43] and also selected only those channels where at least two charged tracks from the HNL decay appear. In the case of the MATHUSLA experiment we optimistically used det = 1 for the detection efficiency.
scalars. We explored analytically the characteristic features of the experiment's sensitivity regions: upper and lower boundaries and the maximal mass of new particles that can be probed. Our analytic analysis allows identifying the parameters responsible for the positions of the main "features" of these curves and to cross-check/validate the results of the Monte Carlo simulations. We analyse a number of experimental factors that contribute to the sensitivity estimates: (i) the number of heavy flavour mesons traveling in the direction of the detector; (ii) their average momentum and the high-energy tail of the momentum distribution; (iii) geometry of the experiment; (iv) efficiency. We use SHiP and MATHUSLA as examples of the fixed target and LHC-accompanying Intensity Frontier experiments, respectively. Our analytic estimates agree well with the Monte Carlo-based sensitivities provided by the SHiP [43] and MATHUSLA [21] collaborations under similar assumptions about the overall efficiencies of the experiments.
Our main results are as follows. Our estimates of the sensitivities of the SHiP and MATHUSLA experiments to the HNLs are shown in Figs. 13 and to the dark scalars in Fig. 14.
Qualitatively both experiments can probe similar ranges of parameters. The SHiP has higher average γ factors of the mesons ( γ ship meson / γ mat meson O(10)) and, as a result, significantly higher upper boundary of the sensitivity region than MATH-USLA (as the upper boundary is exponentially sensitive to the γ factor). As the consequence, the SHiP can probe higher masses for both HNLs and scalars than MATHUSLA. However, the W boson decays at the LHC would produce some highly boosted HNLs traveling to the MATHUSLA decay volume, partly mitigating this difference.
The SHiP experiment is able to probe lower mixing angles for HNLs with M N m Ds owing to the larger number of D mesons. MATHUSLA can probe lower mixing angles for the HNLs with M N m Ds and the scalars for all masses, owing to the larger number of the B +/0 mesons at the LHC (as charmed mesons contribute negligibly to the scalar production).
Uncertainties. According to the theoretical predictions the dσ/dp T distribution of B mesons at the LHC has a maximum at p T ∼ GeV, see Fig. 4. The region of low p T is complicated for the theoretical predictions because of limitations of the applicability of the perturbative QCD. At the same time, these cross-sections have not been measured by neither the ATLAS, nor the CMS collaborations in the required kinematic range. The increase in the overall amount of low-momentum mesons shifts leftwards the position of the peak of the dσ/dp T distribution, thus decreasing their average momentum. Both factors lead to the increase of the number of events at the lower boundary. Therefore the uncertainty in the position of the lower boundary of the sensitivity region depends on both of these numbers such that the uncertainty in the position of the peak enters into the sensitivity estimate squared.
Another uncertainty comes from the background estimates. For the SHiP experiment, comprehensive background studies have proven that the yield of background events passing the online and offline event selections is negligible [24,26]. For MATH-USLA such an analysis is not available at the time of writing. The Standard Model background at MATHUSLA is non-zero (due to neutrinos from LHC and atmosphere, cosmic rays, muons, etc), however, it is claimed to be rejected with high efficiency based on the topology of the events [17,21]. It is not known how much this rejection affects the detection efficiency, det . In this work, we conservatively assumed det = 1 for MATHUSLA, while for SHiP it was taken from the actual Monte Carlo simulations [43]. More detailed analysis of the MATHUSLA background should be performed, which could influence the sensitivities.
In case of the SHiP experiment, the main uncertainty for HNLs is the unknown production fraction of the B c mesons at √ s ∼ 30 GeV. It changes the position of the lower bound and consequently the maximal mass probed in a significant way, see Fig. 13.
Comparison with other works. We have compared our sensitivity estimates with the results of the Monte Carlo simulations presented by the SHiP collaboration [43,71,72] and with the estimates of the MATHUSLA physics case paper [21] (Figs. [11][12]. For the HNLs, the estimates are in good agreement with the results of the SHiP collaboration. In the case of MATHUSLA, there is a difference for HNLs with mass smaller than 1 GeV. It can be attributed to different branching for the HNL production used in our estimates and in the Monte Carlo simulations of [21], see discussion in Appendix H. For the scalars, our estimates are in good agreement with the results from the SHiP and MATHUSLA collaboration. Small discrepancies between the sensitivities at the upper bound can be explained mainly by the difference in the meson energy spectrum used in our estimation and obtained in the Monte Carlo simulations.  Table 3. The fragmentation fractions of the at the energies of the LHC [46,79,80] and of the SHiP experiment [81]. where g is the coupling constant and L int are interaction terms that play no role in our analysis. After the spontaneous symmetry breaking the cubic term in (A.1) gives rise to the Higgs-like interaction of the scalar S with all massive particles with their mass times a small mixing parameter where g is the coupling in (A.1); v is the Higgs VEV; m H is the Higgs mass; sum in (A.2) goes over all massive fermions (leptons and quarks); W ± µ is the W boson and · · · denote other interaction terms, not relevant for this work. The details of the phenomenology of the scalar portal are provided in [45] (see also [37,[76][77][78]). The computation of hadronic decay width of S is subject to large uncertainties at masses M S ∼ few GeV, where neither chiral perturbation theory not perturbative QCD can provide reliable results (see a discussion in [75]).
We have also considered the neutrino portal where one adds to the Standard Model new gauge-singlet fermion -heavy neutral lepton N -that couples to the abLa H b where L a is the SU(2) lepton doublet and ab is absolutely antisymmetric tensor in 2 dimensions. Phenomenologically, HNL is massive Majorana particle that possesses "neutrino-like" interactions with W and Z bosons (the interaction with the Higgs boson does not play a role in our analysis and will be ignored). The interaction strength is suppressed as compared to that of ordinary neutrinos by a flavour-dependent factors (known as mixing angles) U α 1 (α = {e, µ, τ }).

B Production and detection of portal particles B.1 Production in proton-proton collisions
The number of mesons is determined by the number of produced qq pairs and fragmentation fractions f meson , that can be extracted from the experimental data [46,79,80]. We summarize the fragmentation fractions that we use for MATHUSLA in the Table 3. For the SHiP experiment, all fragmentation fractions except for B c meson are known to be close to the MATHUSLA's ones [81]. The B c meson fragmentation fraction at the energy of the SHiP experiment is unknown. In our estimations, we take it the same as for the MATHUSLA experiment. The production of the HNL in the decay of charmed and beauty mesons have been considered in [82,83], see [42] for the recent review and summary of the results. The branching ratios multiplied by fragmentation fractions of D and B for the most relevant channels and the values of the fragmentation fractions from the Table 3 are presented at the Fig. 15. We see that that for the HNL mass range m N 3.5 GeV the main production channel is B c meson decay B c → N + l. This is a quite surprising fact, taking into account that B c fragmentation fraction is of order 10 −3 . To understand this result let us compare HNL production from B c with production from the two-body B + decay. The decay widths for both cases are given by where we take m N m , and K is a kinematic suppression. Neglecting them, for the ratio for the numbers of HNLs produced by B c and B + we obtain We see that the small fragmentation fraction of B c meson is compensated by the ratio of the CKM matrix elements and meson decay constants. HNLs can also be produced in the decays of the W bosons, W → N + l. The corresponding branching ratio is where we have neglected the HNL and the lepton masses.

B.1.2 Scalar production
The main difference in the phenomenology of the Higgs-like scalar S in comparison to HNLs is that the interaction of S with fermions is proportional to their mass, see Sec. A. Therefore, its production at the mass range M S > M K is dominated by the decay of the B + , B 0 , while the contribution from D mesons is negligible [44,45]. The main production process is the 2-body decay where X q is a hadron that contains the quark q. The branching ratios for these states will be discussed in the forthcoming paper [45]. Here we only state the main points. For B → X s + S we choose two lightest resonances X s for each given spin and parity.
Exceptions are pseudo-scalar and tensor mesons (there is only one known meson that has these properties, see [84]). We have found that each heavier meson from the "family" gives a smaller contribution to the branching ratio than the lighter one. For B → X d + S we take only one meson X d = π since this channel has the largest kinematic threshold m B − m π . We summarize the list of the final states below: • Spin 0, odd parity: X q = π, K; • Spin 0, even parity: X q = K * 0 (700), K * 0 (1430); • Spin 1, odd parity: K 1 (1270), K 1 (1400); • Spin 1, even parity: K * (892), K * (1410); • Spin 2, even parity: K * 2 (1430). The main source of the uncertainty is unknown quark squad of the K * 0 (700) meson: it can be either a di-quark or a tetra-quark (see e.g. [74]). In the second case, the K * 0 (700) contribution to the scalar production is unknown, which causes an uncertainty up to 30%. We consider it as the di-quark state.
The dependence of the branching ratios of the process (B.4) on the scalar mass is shown in Fig. 16.
We estimate the production from of the scalars the W bosons by the decay W → S + f +f , where the summation over all the SM fermions species f = l, q is taken. We obtained BR W →S /θ 2 4 · 10 −3 .

B.2 Main decay channels B.2.1 HNL
The HNL has 3-body leptonic decays and different semileptonic modes. Following the paper [42], we estimate the decay width of HNL into hadronic states as a sum of decay widths of specific channels for the HNL with a mass lower as 1 GeV and use the decay into the quarks with QCD corrections for larger masses. In the latter mass region, the HNL decay width can be approximated by the formula where N eff = N ν + N + N c N q is the effective number of channels, N ν = 1 takes into account decays into neutrino, N and N q is a number of kinematically allowed lepton or quark generations and N c = 3 is a number of the QCD colors. The dependence of the proper lifetime cτ N on the HNL mass at U 2 = 1 is given at the left panel of Fig. 17.

B.2.2 Scalar
The decay width of the scalar particle has large uncertainty in the scalar mass region 0.5 GeV < M S < 2 GeV because of resonant nature of S → 2π decay [77]. At higher scalar mass the decay width is fixed by the perturbative QCD calculations [85]. We omit the problem of pion resonance in this work using continuous interpolation between the sum of the decay channel at low mass and perturbative QCD at high ones.
For scalar mass region above 2 GeV one can naively estimate S decay width as Γ S ∝ f θ 2 y 2 f M S . This estimation does not take into account three effects: 1. For the decay into quarks parameter y q depends on scalar mass as y q ≡ M q (M 2 S )/v, where M q (M 2 S ) is quark running mass, which gives logarithmic correction; 2. The decay into gluons has different M S dependence, Γ S ∝ θ 2 M 3 S /v 2 , and dominates in the region 2 GeV < M S < 3.5 GeV [44]; 3. In the region M S near 3.5 GeV new decay channels appear (into τ and c quark), and the kinematical factor is important.
Taking them into account, for the mass domain 3.5 GeV < M S < 5 GeV, near the threshold of production from B mesons, we made a fit to the total Γ S and found that its behavior is Γ S ∝ M 2 S . The dependence of the proper lifetime cτ S on the scalar mass at θ = 1 is given at the right panel of Fig. 17.

B.2.3 Visible branching ratio
We define the "visible" decay channels as those that contain at least two charged particles α in the final state. The corresponding decays are where Y is arbitrary state, F is uncharged state that decays to n charged particles andỸ is a state with at least 2 − n charged particles (assuming n < 2). Using this definition, the decay N → 3ν is identified as invisible decay, the decay N →μν µ e -as visible decay, while the decay N → ην -as the visible decay if η meson decays into two charged particles, and as invisible decay otherwise. To take into account only visible decays of F , we include the factor BR F →vis to the partial decay width Γ X→FỸ . We take the values of BR F →vis from [84]. For HNL/scalar masses M > 2 GeV we describe hadronic decays as having quarks and gluons in the final states. In this case we assume that any such decay will contain at least 2 charged tracks and therefore the whole hardonic width is visible.

C Relation between momentum of HNL and meson momentum
The energy of the particle X at lab frame, E X , is related to energy of X at meson's rest frame, E rest X , and meson energy E meson at lab frame as where θ is the angle between the direction of motion of meson in lab frame and the direction of motion of the particle X in the meson's rest frame. At meson frame the angle distribution is isotropic, so for the average energy we obtain where γ meson ≡ E meson /m meson .

D Geometry of the experiments D.1 SHiP
The SHiP experiment [26] is a fixed-target experiment using the proton beam of the Super Proton Synchrotron (SPS) at CERN. SPS can deliver N p.o.t. = 2 × 10 20 protons with the energy 400 GeV over a 5 year term. The SHiP will be searching for new physics in the largely unexplored domain of very weakly interacting particles with masses below O(10) GeV and cτ exceeding tens of meters. The overview of the experiment is as follows.
The proton beam hits a target [24,25]. The target will be followed by a 5 m hadron stopper, intended to stop all π ± and K mesons before they decay, and by a system of shielding magnets called active muon shield, constructed to sweep muons away from the fiducial decay volume. The whole active muon shield system is 34 m long.
The decay volume is a long pyramidal frustum vacuum camera with the length downstream of the primary target respectively. The SHIP spectrometer downstream of the decay volume consists of a four-station tracker, timing detector, and an electromagnetic calorimeter and muon detector for particle identification. The detectors are seen from the interaction point at an angle θ ≈ 3.2 · 10 −3 rad. Figure 18. MATHUSLA experiment geometry. Adapted from [18] MATHUSLA (MAssive Timing Hodoscope for Ultra Stable neutraL pArticles) is a proposed experiment [17,22] that consists of a 20 m × 200 m × 200 m surface detector, installed above ATLAS or CMS detectors (see Fig. 18). The long-lived particles, created at the LHC collisions, travel 100+ meters of rock and decay within a large decay volume (8×10 5 m 3 ). Multi-layer tracker at the roof of the detector would catch charged tracks, originating from the particle decays. The ground between the ATLAS/CMS and MATHUSLA detector would serve as a passive shield, significantly reducing the Standard Model background (with the exception of neutrinos, muons and K 0 L created near the surface). Assuming isotropic angular distribution of a given particle traveling to MATHUSLA, the average distance that it should travel to reach the MATHUSLA decay volume is equal tō

D.2 MATHUSLA
where L ground = 100 m. The average distance a particle travels inside the decay volume,L det , is given byl  Table 4.

E Analytic estimation of the upper bound: details
In this Section we estimate the ratio between θ max and θ upper -the quantities defined in Section 3.  Table 4. Parameters of MATHUSLA experiment [21]. For the definition of angles θ 1,2 see Fig. 18, and ∆φ is the azimuthal size of MATHUSLA

E.1 Fits of the spectra
The high-energy tail of the B mesons distribution function at SHiP is well described by the exponential distribution, see the left Fig. 19: The distribution of the high energy B mesons at the LHC for energies E B 300 GeV can be approximated by the power law function, see the right panel in Fig. 19: ≈f 0 E −α ,f 0 ≈ 1.6 × 10 4 GeV α−1 and α 4.6 (E.2) Finally, the distribution of the HNLs originating from the W bosons can be approximated by the expression

E.2 Upper bound estimation
We start from the number of the events N events (M X , θ 2 X ) =Ñ prod (M X , θ 2 X ) ×P decay (E.4)

F Details of the sensitivity curve drawing
We draw the sensitivity curve for HNLs and scalars by using the formula N events (θ X , M X ) = meson N events,meson + N events,W , (F.1) where the numbers of the particles from the mesons and the W bosons is estimated as N events,meson = N qq × f q→meson Br meson→X × meson dp meson f pmeson ×P decay (p meson ), N events,W = N W,LHC × Br W →X × W × dp X f p X ,W ×P decay (p X ) (F.3) Here N qq is the total number of the qq pairs that are produced in pp collisions at the high luminosity LHC or in p − target collisions at SHiP, f pmeson is the momentum distribution of the mesons that fly to the decay volume of the experiment (see Sec. 5.1), and is the overall efficiency (see Sec. 6.1). In the expression for the decay probabilityP decay the γ factor of the X particle is related to the meson momentum by the relation (C.2). N W,LHC is the number of the W bosons produced at the high luminosity LHC, and f p X ,W is the momentum distribution function of the X particles (see Sec. 5.3).

G Analytic estimation of the lower bound for particular masses
Here we make an analytic estimation of the lower bound for particular masses using the formula (2.1). The relevant parameters and the value θ 2 X,lower for the SHiP and MATHUSLA experiments are given in Tables 5, 6.  Table 5. Table of parameters used in simple analytic estimation of the lower bound (2.1) for the SHiP experiment for particular masses of the HNLs with the pure electron mixing and the scalars. The columns are as follows: the type of the particle X and the mother particle, the branching ratio of the production of X at θ 2 X = 1, average γ factor, proper decay length cτ X at θ 2 X = 1, overall efficiency (6.1), the mixing angle at the lower bound estimated as N events,lower = 2.3, where N events,lower is given by (2.1).

H HNLs at MATHUSLA for small mass
For M N m Ds , where the sensitivity curve is determined by the production channel from D mesons, in the range 1 GeV M N M D the curves are in agreement, while for smaller masses we have a discrepancy of order of a factor 3 (which corresponds to a difference in a factor 10 between the number of events) and different shapes of the sensitivity curve. There are a few possible reasons for this: a) difference in numbers of the D mesons that are produced in the direction of the MATHUSLA's decay volume, b) extra amount of the small mass HNLs that are produced from the mesons that does not fly in the decay volume of MATHUSLA (they were not taken into account in our estimate), c) different HNL phenomenology (production and decay) used in comparison.
Since the positions of the lower bounds in the mass range 1 GeV M N m Ds are in good agreement, we conclude that case a) with different amounts of D mesons is less probable, as it should shift the lower bound for the whole mass range M N m D . Next, in order to estimate the amount of very light HNLs produced from all D mesons, we perform a MadGraph 5 simulation of a process pp →ce + ν e s, which corresponds to the main process for a production of light HNLs M N 1 GeV -D → N e + K. We computed the ratio χ = σc e + νes, mat σc e + νes, tot / σ cc, mat σ cc, tot , (H.1) where the first fraction is the amount of the HNLs that fly in the decay volume of MATHUSLA, while the second one is the amount of cc pairs that fly in the same direction. We found χ ≈ 1.7, which is not enough to explain the discrepancy. For c), we compared the decay widths of the HNLs used in our estimate and in [21] (see [86]), and found that they are in good agreement, see Fig. 20. We have not . The ratio of the decay widths of the HNL at U 2 e = U 2 µ = U 2 τ = 1 used in our work and in the simulations of the MATHUSLA collaboration [21] using the results of [86].
found the information about the production branching ratios used in [21]. As we see, the cases a) and b) are not enough to explain a factor 10 in the number of events, and we assume that the main reason for the discrepancy is different production branching ratios.