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

TL;DR

This paper provides analytical estimates of the sensitivity regions for intensity-frontier experiments searching for feebly interacting particles, illustrating how experimental parameters influence the maximal mass reach and boundaries, with validation against simulations.

Contribution

It introduces analytical methods to estimate sensitivity regions of experiments like SHiP and MATHUSLA, enhancing understanding of their capabilities without extensive simulations.

Findings

Analytical estimates match Monte Carlo simulation results.

Sensitivity boundaries depend on experimental parameters.

Maximal mass reach can be analytically predicted.

Abstract

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.

Figures

Tables

Table 1. Table 1: Parameters of the SHiP and MATHUSLA experiments: the total number of all charmed/beauty hadrons; the average γ 𝛾 \gamma factor of mesons flying in the direction of the decay volumes of the experiment; the geometric acceptances for these hadrons. We take B c subscript 𝐵 𝑐 B_{c} meson distribution to be proportional to that of B + superscript 𝐵 B^{+} mesons, scaled by f B c subscript 𝑓 subscript 𝐵 𝑐 f_{B_{c}} . As a result, B c subscript 𝐵 𝑐 B_{c} gamma factor is the same as for B + superscript 𝐵 B^{+} mesons for SHiP and scaled by m B / m B c subscript 𝑚 𝐵 subscript 𝑚 subscript 𝐵 𝑐 m_{B}/m_{B_{c}} for MATHUSLA, see discussion in Section  5.1 and  5.2 For SHiP we assumed 5 years of operation ( 2 × 10 20 2 superscript 10 20 2\times 10^{20} protons on target) and for MATHUSLA we took the luminosity of the HL phase, ℒ h = 3000 ​  fb − 1 subscript ℒ ℎ 3000 superscript  fb 1 \mathcal{L}_{h}=3000\text{ fb}^{-1} . Predictions are based on the FairSHiP simulations (SHiP) and on the FONLL simulations (MATHUSLA). See text for details.

Experiment	$N_D$	$N_B$	$⟨{\gamma }_D⟩$	$⟨{\gamma }_B⟩$	${ϵ}_{\text{geom},D}$	${ϵ}_{\text{geom},B}$
MATHUSLA	$4.4\times {10}^{16}$	$3\times {10}^{15}$	2.6	2.3	$1.3\times {10}^{-2}$	$1.8\times {10}^{-2}$
SHiP	$1.6\times {10}^{18}$	$1.1\times {10}^{14}$	19.2	16.6	$-$	$-$

Table 2. Table 2: Parameters of the SHiP and MATHUSLA experiments: the average length from the interaction point to the decay volume ⟨ l target-det ⟩ delimited-⟨⟩ subscript 𝑙 target-det \langle l_{\text{target-det}}\rangle , the average length of the decay volume ⟨ l det ⟩ delimited-⟨⟩ subscript 𝑙 det \langle l_{\text{det}}\rangle (see Appendix  D for details), values of the overall efficiencies ( 25 ) averaged over the probed mass range of X 𝑋 X for the particles X 𝑋 X produced from D 𝐷 D and B 𝐵 B mesons, the effective number of the D 𝐷 D and B 𝐵 B mesons and W 𝑊 W bosons defined by N M , eff = N M × ϵ ¯ X , M subscript 𝑁 𝑀 eff subscript 𝑁 𝑀 subscript ¯ italic-ϵ 𝑋 𝑀 N_{M,\text{eff}}=N_{M}\times\bar{\epsilon}_{X,M} .

Exp.	$⟨l_{\text{target-det}}⟩$	$⟨l_{\text{det}}⟩$	${\overline{ϵ}}_{X,D}$	${\overline{ϵ}}_{X,B}$	${\overline{ϵ}}_{X,W}$	$N_{D,\text{eff}}$	$N_{B,\text{eff}}$	$N_{W,\text{eff}}$
MAT	$192\text{ m}$	$38\text{ m}$	$0.013$	$0.018$	0.02	$5.7\cdot {10}^{14}$	$5.4\cdot {10}^{13}$	$1.2\cdot {10}^{10}$
SHiP	$50\text{ m}$	$50\text{ m}$	$0.09$	$0.12$	–	$1.4\cdot {10}^{17}$	$1.3\cdot {10}^{13}$	–

Table 3. Table 3: The fragmentation fractions for heavy mesons at the LHC energies  Aaij:2011jp ; Aaij:2017kea ; Kling:2018wct and of the SHiP experiment  Graverini:2133817 ; Bondarenko:2018ptm . For SHiP the contribution of flavoured baryons or quarkonia states can be neglected, see  Bondarenko:2018ptm . For the LHC energies, the remaining 20-25% come of all heavy flavour quarks hadronize into baryons, mostly Λ b subscript Λ 𝑏 \Lambda_{b} states  Aaij:2011jp .

Meson M	$B^{+}$	$B^0$	$B_s$	$B_c^{+}$	$D^{+}$	$D^0$	$D_s^{+}$
MATHUSLA	0.324	0.324	0.088	$2.6\cdot {10}^{-3}$	0.225	0.553	0.105
SHiP	0.417	0.418	0.09	?	0.207	0.632	0.088

Table 4. Table 4: Parameters of MATHUSLA experiment  Curtin:2018mvb . For the definition of angles θ 1 , 2 subscript 𝜃 1 2 \theta_{1,2} see Fig.  17 , and Δ ​ ϕ Δ italic-ϕ \Delta\phi is the azimuthal size of MATHUSLA

Parameter	${\theta }_1$	${\theta }_2$	${\eta }_1$	${\eta }_2$	${\overline{l}}_{\text{target-det}}$, m	${\overline{l}}_{\text{det}}$, m	$\Delta \varphi$
Value	${44.3}^{\circ }$	${22.9}^{\circ }$	0.9	1.6	192.5	38.5	$\pi /2$

Table 5. Table 5: Table of parameters used in simple analytic estimation of the lower bound ( 8 ) 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 𝑋 X and the mother particle, the branching ratio of the production of X 𝑋 X at θ X 2 = 1 superscript subscript 𝜃 𝑋 2 1 \theta_{X}^{2}=1 , average γ 𝛾 \gamma factor, proper decay length c ​ τ X 𝑐 subscript 𝜏 𝑋 c\tau_{X} at θ X 2 = 1 superscript subscript 𝜃 𝑋 2 1 \theta_{X}^{2}=1 , overall efficiency ( 25 ), the mixing angle at the lower bound estimated as N events , lower = 2.3 subscript 𝑁 events lower 2.3 N_{\text{\rm events},\text{lower}}=2.3 , where N events , lower subscript 𝑁 events lower N_{\text{\rm events},\text{lower}} is given by ( 8 ).

$X/M$	$M_X,\text{ GeV}$	$f_M\times {\text{BR}}_{M\to X}$	$⟨{\gamma }_X⟩$	$c{\tau }_X,\text{ m}$	$ϵ$	$U_{e,\text{lower}}^2$
$N,{\nu }_e/D$	$0.5$	$4.5\cdot {10}^{-2}$	$30.1$	$1.1\cdot {10}^{-2}$	$6.5\cdot {10}^{-2}$	$1.9\cdot {10}^{-9}$
$N,{\nu }_e/B$	$3$	$4\cdot {10}^{-4}$	$20.2$	$2\cdot {10}^{-6}$	$0.13$	$1.7\cdot {10}^{-8}$
$S/B$	$0.5$	$88.2$	$30.1$	$4.3\cdot {10}^{-9}$	$0.14$	$2.5\cdot {10}^{-11}$
$S/B$	$2.5$	$21.4$	$20.2$	$1.\cdot {10}^{-9}$	$0.2$	$4.8\cdot {10}^{-12}$

Table 6. Table 6: Table of parameters used in simple analytic estimation ( 8 ) for the MATHUSLA experiment for particular masses of the HNLs and scalars. We use the description of the scalar phenomenology from  Curtin:2018mvb . The columns are as follows: the type of the particle X 𝑋 X and the mother particle, the branching ratio of the production of X 𝑋 X at θ X 2 = 1 superscript subscript 𝜃 𝑋 2 1 \theta_{X}^{2}=1 , average γ 𝛾 \gamma factor, proper decay length c ​ τ X 𝑐 subscript 𝜏 𝑋 c\tau_{X} at θ X 2 = 1 superscript subscript 𝜃 𝑋 2 1 \theta_{X}^{2}=1 , overall efficiency ( 25 ), the mixing angle at the lower bound estimated as N events , lower = 4 subscript 𝑁 events lower 4 N_{\text{\rm events},\text{lower}}=4 , where N events , lower subscript 𝑁 events lower N_{\text{\rm events},\text{lower}} is given by ( 8 ).

X/M	$M_X,\text{ GeV}$	$f_M\times {\text{BR}}_{M\to X}$	$⟨{\gamma }_X⟩$	$c{\tau }_X,\text{ m}$	$ϵ$	${\theta }_{\text{lower}}^2$
$N,{\nu }_e/D$	$0.5$	$4.5\cdot {10}^{-2}$	$2.4$	$1.1\cdot {10}^{-2}$	$6.5\cdot {10}^{-2}$	$1.2\cdot {10}^{-8}$
$N,{\nu }_e/B$	$2.5$	$2.6\cdot {10}^{-3}$	$2.2$	$5.4\cdot {10}^{-6}$	$1.8\cdot {10}^{-2}$	$2.4\cdot {10}^{-9}$
$N,{\nu }_e/W$	$1$	$0.11$	$62$	$6.5\cdot {10}^{-2}$	$1.9\cdot {10}^{-2}$	$1.7\cdot {10}^{-6}$
$N,{\nu }_e/W$	$2.5$	$0.11$	$24.8$	$5.4\cdot {10}^{-6}$	$1.9\cdot {10}^{-2}$	$1.1\cdot {10}^{-7}$
$S/B$	$0.5$	$5.3$	$9$	$4.3\cdot {10}^{-9}$	$1.8\cdot {10}^{-2}$	$4.2\cdot {10}^{-12}$
$S/B$	$2.5$	$2.7$	$2$	$1.\cdot {10}^{-9}$	$0.1$	$4.9\cdot {10}^{-13}$

Equations

$$N_{\text{events}}=\sum _M{N}_{\text{prod},M}\times P_{\text{decay},M},$$

$$N_{\text{prod},M}\approx N_{\text{M}}\times {\mathrm{BR}}_{M\to X}\times {ϵ}_{\text{decay},M}$$

$$P_{\text{decay},M}=\left[\exp \left(-\frac{l_{\text{target-det}}}{l_{\text{decay}}}\right)-\exp \left(-\frac{l_{\text{target-det}}+l_{\text{det}}}{l_{\text{decay}}}\right)\right]\times {ϵ}_{\det }\times {\text{BR}}_{\mathrm{vis}},$$

$$l_{\text{decay}}=c{\tau }_X{\beta }_X{\gamma }_X,$$

$${\mathrm{BR}}_{\text{meson}\to X}\propto {\theta }_X^2,{\tau }_X\propto {\theta }_X^{-2}$$

$$N_{\text{events},\text{lower}}\propto {\theta }_X^4/{\gamma }_X$$

$$N_{\text{events},\text{upper}}\propto {\theta }_X^2{e}^{-C{\theta }_X^2/{\gamma }_X},$$

$$N_{\text{events},\text{lower}}\approx N_{\text{meson}}\times {\mathrm{BR}}_{M\to X}\times \frac{⟨l_{\det }⟩}{c{\tau }_X⟨{\gamma }_X⟩}\times {ϵ}_X,$$

$${\gamma }_X\approx {\gamma }_{\text{meson}}\frac{⟨E_X^{\text{rest}}⟩}{M_X},$$

$$\boxed{\frac{({\theta }_{X,\text{lower}}^{\text{SHiP}}{)}^2}{({\theta }_{X,\text{lower}}^{\text{MAT}}{)}^2}=\sqrt{\frac{N_{\text{events}}^{\text{textsc}mat}}{N_{\text{events}}^{\text{textsc}ship}}}\simeq \sqrt{\frac{N_{\text{meson}}^{\text{textsc}mat}}{N_{\text{meson}}^{\text{textsc}ship}}\times \frac{l_{\det }^{\text{textsc}mat}}{l_{\det }^{\text{textsc}ship}}\times \frac{⟨{\gamma }_{\text{meson}}^{\text{textsc}ship}⟩}{⟨{\gamma }_{\text{meson}}^{\text{textsc}mat}⟩}\times \frac{{ϵ}_{\text{textsc}mat}}{{ϵ}_{\text{textsc}ship}}}},$$

$$l_{{\text{\rm decay}}}\bigl{(}M_{X},\theta^{2}_{\text{max}}\bigr{)}\simeq\left\{\begin{aligned} 1.5l_{\text{target-det}},&\quad\text{if}\quad l_{\text{target-det}}\simeq l_{\text{det}}\\ 0.5l_{\text{target-det}},&\quad\text{if}\quad l_{\text{target-det}}\gg l_{\text{det}}\\ \end{aligned}\right.$$

$$N_{\text{events}}\simeq N_{\text{prod}}\times {ϵ}_{\det }\times {\mathrm{BR}}_{\text{vis}}\times e^{-l_{\text{target-det}}/l_{\text{decay}}}$$

$$\boxed{\theta^{2}_{\text{upper}}\simeq\theta^{2}_{\text{max}}\times\log\Bigl{[}N_{\text{prod}}(\theta^{2}_{\text{max}})\,\epsilon_{\det}\mathop{\text{BR}}\nolimits_{\text{vis}}\Bigr{]}}.$$

$$f_X(E_X)=\frac{1}{N_X}\frac{d{N}_X}{d{E}_X}.$$

$$P_{\text{decay}}={ϵ}_{\det }\times {\mathrm{BR}}_{\text{vis}}\times \underset{0}{\overset{\infty }{\int }}d{E}_X{f}_X(E_X)\times \pi \left(\frac{{\tau }_X{p}_X}{l_{\text{target-det}}{M}_X}\right),$$

$$\pi (y)\equiv \exp \left(-\frac{1}{y}\right)-\exp \left(-\frac{l_{\text{target-det}}+l_{\text{det}}}{l_{\text{target-det}}}\frac{1}{y}\right),$$

$$l_{\text{decay}}(M_{X,\text{max}},{\theta }_{\text{lower}}^2(M_{X,\max }))\simeq l_{\mathrm{target}-\det }.$$

$$M_{X,\max }\propto {\left(\frac{⟨E_X⟩}{|{\theta }_{\text{lower}}{|}^2{l}_{\mathrm{target}-\det }}\right)}^{\frac{1}{\alpha +1}},$$

$$\boxed{\frac{M_{X,\max }^{\text{textsc}ship}}{M_{X,\max }^{\text{textsc}mat}}\simeq {\left(\frac{⟨E_X{⟩}^{\text{textsc}ship}}{⟨E_X{⟩}^{\text{textsc}mat}}\times \frac{|{\theta }_{X,\text{lower}}^{\text{textsc}mat}{|}^2}{|{\theta }_{X,\text{lower}}^{\text{textsc}ship}{|}^2}\times \frac{l_{\text{target-det}}^{\text{textsc}mat}}{l_{\text{target-det}}^{\text{textsc}ship}}\right)}^{\frac{1}{\alpha +1}}}.$$

$$\frac{N_D{|}_{p_T<3.5\text{ GeV}}}{N_D{|}_{p_T>3.5\text{ GeV}}}=3.8,\frac{N_B{|}_{p_T<10\text{ GeV}}}{N_B{|}_{p_T>10\text{ GeV}}}=5.7$$

$$N_{W,\text{total}}\approx 6\cdot 10^{11}$$

$$\frac{d^2{N}_X^W}{d{E}_{X}d\cos ({\theta }_X)}=\int d{p}_W\frac{d{N}_W}{d{p}_W}\times \frac{d^2{\mathrm{BR}}_{W\to X}}{d{\theta }_{X}d{E}_X}\times P({\theta }_X)$$

$$N_{\text{prod},W}/N_{\text{prod},B}\approx \frac{N_W}{N_B}\times \frac{{\mathrm{BR}}_{W\to X}}{{\mathrm{BR}}_{B\to X}}\times \frac{{ϵ}_{\text{decay},W}}{{ϵ}_{\text{decay},B}}\approx \{\begin{array}{l}10^{-3}{ϵ}_{\text{decay},W},\text{scalars}\\ 10{ϵ}_{\text{decay},W},\text{HNLs},\end{array}$$

$$\frac{d^2{\mathrm{BR}}_{W\to N}}{d{\theta }_{N}d{E}_N}=\frac{1}{{\Gamma }_W}\frac{|{\mathcal{M}}_{W\to e+N}{|}^2}{8\pi }\frac{p_N}{E_W}\delta (M_N^2+m_W^2-2E_N{E}_W+2|{\mathbf{p}}_N||{\mathbf{p}}_W|\cos ({\theta }_N))$$

$$ϵ={ϵ}_{\text{decay}}\times {ϵ}_{\det }\times {\mathrm{BR}}_{\text{vis}},$$

$$\frac{U^{2}_{\text{lower},\textsc{ship}}}{U^{2}_{\text{lower},{\textsc{mat}}}}\bigg{|}_{M_{N}\lesssim m_{D}}\approx\frac{1}{5},\quad\frac{U^{2}_{\text{lower},\textsc{ship}}}{U^{2}_{\text{lower},{\textsc{mat}}}}\bigg{|}_{M_{N}\gtrsim m_{D}}\simeq\frac{\theta^{2}_{\text{lower},\textsc{ship}}}{\theta^{2}_{\text{lower},{\textsc{mat}}}}\bigg{|}_{M_{S}\gtrsim m_{K}}\approx 5$$

$$R={\theta }_{\text{upper}}^2/{\theta }_{\text{max}}^2$$

$$M_{\max }^{N,\text{textsc}ship}/M_{\max }^{N,\text{textsc}mat}\approx 1.3,$$

$${\mathcal{L}}_{\text{scalar}}={\mathcal{L}}_{\mathrm{SM}}+\frac{1}{2}({\partial }_{\mu }S{)}^2-\frac{M_S^2}{2}{S}^2+gS{H}^{†}H+{\mathcal{L}}_{\mathrm{int}}$$

$$\mathcal{L}_{S,{\rm int}}=\theta S\biggl{[}\sum_{f}m_{f}\bar{f}f+M_{W}W^{+}_{\mu}W^{-}_{\mu}+\dots\biggr{]}\quad\theta\equiv\frac{gv}{m_{H}}\ll 1$$

Full text

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, see e.g. Beacham:2019nyx for an overview. 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 $\leq 4$) or be realized as higher-dimensional operators suppressed by the dimensional couplings $\Lambda^{-n}$, with $\Lambda$ 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^{\dagger}H$ term constructed from a Higgs doublet field $H_{a}$, $a=1,2$;

• –

a neutrino portal that couples new gauge singlet fermion to the $\epsilon_{ab}\bar{L}_{a}H_{b}$ where $L_{a}$ is the SU(2) lepton doublet and $\epsilon_{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) hypercharge field strength.

Let us denote a new particle by $X$. The interaction of $X$ with the SM is controlled by the mixing angle $\theta_{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 Hyun:2010an ; Lees:2012ra ; Adams:2013qkq ; Lees:2014xha ; Aaij:2014aba ; Aad:2015xaa ; Khachatryan:2015gha ; TheBelle:2015mwa ; Banerjee:2016tad ; Lees:2017lec ; Mermod:2017ceo ; Izmaylov:2017lkv ; Aaij:2017rft ; Dobrich:2017yoq ; Sirunyan:2018mtv . Although the LHC is a flagship of the Energy Frontier exploration, its high luminosity (especially in the Run 3 and beyond) means that huge numbers of heavy flavored mesons and vector bosons are created. This opens the possibility of supplementing 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 Gligorov:2017nwh , MATHUSLA Chou:2016lxi ; Curtin:2018mvb , FASER Feng:2017uoz ; Ariga:2019ufm , and AL3X Gligorov:2018vkc . 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.

Detailed Monte Carlo (MC) simulations of both production and decays, complemented with background studies and detector simulations, offer ultimate sensitivity curves for each of the experiments.

Such simulations are however difficult to reproduce and modify. The modifications are nevertheless routinely needed because

1. (a)

Geometrical configurations of most experiments are not fully fixed yet and it is important to explore changes of the science reach with the modification of experimental designs; 2. (b)

Production or decays of GeV-mass feeble interacting particles involving quarks and mesons often requires the description outside of the validity range of both perturbative QCD and low-energy meson physics and is, therefore, subject to large uncertainties. This is the case for example for both scalar and neutrino portals (see e.g. Bezrukov:2009yw ; Bondarenko:2018ptm ; Monin:2018lee ; Bezrukov:2018yvd ; Winkler:2018qyg ; Boiarska:2019jym as well as the discussion in Section 7). In particular,

• –

Different groups use different prescription for scalar production Krnjaic:2015mbs ; Curtin:2018mvb ; Winkler:2018qyg ; Beacham:2019nyx ; Boiarska:2019jym

• –

the decay width and hadronic branching fractions for scalars with masses from $\sim 0.5$ GeV to few GeV are subject to large uncertainties, see Monin:2018lee ; Bezrukov:2018yvd ;

• –

multi-hadronic HNL decays are not accounted for by any of the existing simulation tools. Yet they account for the largest part of the HNLs with masses around few GeV Bondarenko:2018ptm ; Dib:2019ztn . 3. (c)

Monte Carlo simulations are done for a limited set of model parameters and it is difficult to explore the overall parameter space and/or modify the sensitivity estimates for extended models (see e.g. the discussion and approach in SHiP:2018xqw )

With this in mind we gathered in one place a sufficiently simple and fully controlled (semi)analytic estimates. Such estimates emphasize the main factors that influence the sensitivity:

(i) dependence on the model (parameters, physical assumptions);

(ii) dependence on the geometry of the experiment;

(iii) factors, related to the beam energy, etc.

We present the final number as a convolution of these factors, which allows to modify any of them at will. As a result one can efficiently compare between several experimental designs; to identify the main factors that influence the sensitivity reach of a particular experiment/model; to reuse existing Monte Carlo sensitivities by separating them into the experimental efficiencies and physical input (model, production/decay phenomenology) with the subsequent modification of one of these factors; to scan over the parameter space of different models as compared to those used in the MC simulations.

It turns out that the ratio between the sensitivities of the experiments to a great extent does not depend on the specific model of new physics, and is determined mainly by the geometry and collision energies of the experiments, which allow a comparison of the sensitivities in a largely model-independent way. To illustrate this point, we compare the potentials of two proposed experiments: the LHC-based MATHUSLA experiment Curtin:2018mvb ; Chou:2016lxi ; Curtin:2017izq ; Evans:2017lvd ; Helo:2018qej and a proton fixed target experiment using the proton beam of the Super Proton Synchrotron (SPS) at CERN – SHiP Anelli:2015pba ; Alekhin:2015byh ; SHIP:2018yqc . We analyze their sensitivity to the neutrino Minkowski:1977sc ; Yanagida:1979as ; Glashow:1979nm ; GellMann:1980vs ; Mohapatra:1979ia ; Mohapatra:1980yp and scalar Silveira:1985rk ; McDonald:1993ex ; Binoth:1996au ; Burgess:2000yq ; Schabinger:2005ei ; Strassler:2006im ; Patt:2006fw ; OConnell:2006rsp portals. For particle masses $M_{X}\lesssim m_{B_{c}}$111By $m_{\dots}$ we denote the masses of lightest flavour mesons, for example, kaons ($m_{K}$), $D^{+}$ ($m_{D}$), $B^{+}$ ($m_{B}$), etc. the main production channels are decays of heavy flavored mesons and $W$ bosons Alekhin:2015byh ; Evans:2017lvd (see also Appendix B.1 for a brief overview). We concentrate on the mass range $M_{X}\gtrsim 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 Dobrich:2017yoq ; Drewes:2018gkc .

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” Anelli:2015pba ; SHiP:2015gkj ; Baranov:2017chy ; Akmete:2017bpl . For MATHUSLA, the background is also expected to be low Chou:2016lxi ; Curtin:2018mvb , although no simulation studies of background have been performed. Even in the most favorable case of $N_{\rm bg}\ll 1$ one needs on average $\bar{N}_{\text{\rm events}}=2.3$ expected signal events to observe at least one event with the probability higher than $90\%$.222To obtain $95\%$ confidence limit one should assume $\bar{N}_{\text{\rm events}}=3$, as the Poisson probability to see at least one event, while expecting 3 “on average” is $0.9502$. 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 discovery.

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_{\text{target-det}}$ (see Appendix D). For such experiments the sensitivity curve has a typical “cigar-like shape” in the plane “mass vs. interaction strength”, see Fig. 1.

The number of decay events in the decay volume factorizes into

[TABLE]

where $N_{\text{\rm prod},M}$ is the number of particles $X$ that are produced from a mother particle $M$ and $P_{{\text{\rm decay}},M}$ is the decay probability. For $N_{\text{\rm prod},M}$ we have

[TABLE]

Here, $N_{M}$ is the number of parent particles produced at the experiment; in the case of mesons $N_{M}=N_{\text{meson}}=2N_{\bar{q}q}\times f_{\text{meson}}$, where $f_{\text{meson}}$ is the fragmentation fraction of a quark $q$ into a given hadron, and $N_{M}=N_{W}$ in the case of the $W$ bosons. $\mathop{\text{BR}}\nolimits_{M\to X}$ is the total branching ratio of decay of the parent particle into $X$ (see Appendix B.1). Finally, $\epsilon_{\text{\rm 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 by333Here we ignored that particles travel slightly different distances depending on their off-axis angle. Eq. (3) also neglects the energy distribution of the produced particles, assuming that all of them travel with the same average energy. This is essential for proper determining of the upper boundary and we will return to this in Section 3.

[TABLE]

where the branching ratio $\mathop{\text{BR}}\nolimits_{\rm vis}$ is the fraction of all decays producing final states that can be registered. Finally, $\epsilon_{\det}\leq 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 $\epsilon_{\det}=1$. The decay length $l_{\text{\rm decay}}$ in Eq. (3) is defined as

[TABLE]

where $\tau_{X}$ is the lifetime of the particle $X$ (see Appendix B.2), $\beta_{X}$ is its velocity and $\gamma_{X}$ is the $\gamma$ factor (which depends on the mother particle that produces $X$).

The production branching ratio and the lifetime behave with the mixing angle as

[TABLE]

At the lower bound of the sensitivity the decay probability behaves as $P_{{\text{\rm decay}}}\propto l_{\det}/l_{{\text{\rm decay}}}$, and as a consequence of (5) the number of events scales as

[TABLE]

At the upper bound $P_{{\text{\rm decay}}}\approx e^{-l_{\text{target-det}}/l_{{\text{\rm decay}}}}$, and

[TABLE]

where $C$ is some numerical factor (that depends on properties of $X$).

Larger $\gamma$ factor suppresses the exponents in the expression for the decay probability (3). From (6), (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 $\gamma$ factor is sensitive to small coupling constants. For sufficiently large couplings, larger $\gamma$ 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, while in Sec. 7 we compare our approach with results of official simulations. 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_{\text{\rm decay}}\gg l_{\det},l_{\text{target-det}}$. For the number of events (1) we have

[TABLE]

where $\epsilon_{X}\equiv\epsilon_{\text{prod}}\times\epsilon_{\text{decay}}\times\mathop{\text{BR}}\nolimits_{\text{vis}}$ is the overall efficiency and $\tau_{X}$ is the lifetime of the particle $X$ (see the discussion below Eq. (4)). The particles are assumed to be relativistic (we will see below when this assumption is justified), so that $\beta_{X}\approx 1$. We estimate the $\gamma$ factor $\gamma_{X}$ from that of the parent meson:

[TABLE]

The average formula (9) does not take into account the distribution of HNLs (scalars) in the meson rest frame – some of the new particles fly in the direction of the parent meson and have $\gamma_{X}$ larger than (9), while the other fly in the opposite direction. We show below that this does not play a role for the lower boundary of the sensitivity curve while the upper boundary is exponentially sensitive to the high $\gamma$-factor tail of the distribution and therefore cannot be determined from Eq. (9). For the experiments like FASER this difference plays an essential role, see Boiarska:FASER .

Since at the lower bound $N_{\text{\rm events}}\propto\theta_{X}^{4}$ (see Eq. (6)), for the ratio of the mixing angles at the lower bound, we have

[TABLE]

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 $\mathop{\text{BR}}\nolimits_{\text{meson}\to X}$ from Eq. (8) disappears. Therefore, to make a comparison between the experiments we only need to know the total number of mesons, their average $\gamma$ factor, the decay volume length and the overall efficiency.

3 Upper boundary of the sensitivity curve

If particles have sufficiently large interaction strength (i.e., the mixing angles), they decay before reaching the decay volume. This determines the upper bound of the sensitivity curve, that we call $\theta_{X,\text{upper}}^{2}$.

A useful quantity to consider is a mixing angle for which the amount of decays inside the decay volume is maximal, $\theta_{X,\text{max}}$. It can be found using the asymptotic behavior for the number of events $N_{\text{\rm events}}$ from the estimations (6), (7). In the domain $l_{{\text{\rm decay}}}\gg l_{\text{target-det}}$ for a fixed mass $M_{X}$ it follows that $N_{\text{\rm events}}$ monotonically grows as $\theta_{X}^{4}$ with the increase of $\theta_{X}$, while in the domain $l_{{\text{\rm decay}}}\ll l_{\text{target-det}}$ it falls exponentially. The position of the maximum $\theta_{\text{max}}$ can be found from

[TABLE]

Using $\theta_{\text{max}}$, we can estimate the value of $\theta_{\text{upper}}$ assuming that all the particles $X$ have the same (average) energy $\langle E_{X}\rangle$. If we neglect the second exponent in the expression for the decay probability (3), then the formula for the number of events (1) becomes

[TABLE]

We can estimate the exponent in (12) as $l_{\text{target-det}}/l_{{\text{\rm decay}}}\approx\theta_{X}^{2}/\theta^{2}_{\text{max}}$, see Eq. (11). So imposing the condition $N_{\text{\rm events}}\simeq 1$ in Eq. (12) with the logarithmic precision we get

[TABLE]

An example of the dependence of the number of events on $\theta_{X}^{2}$ for the fixed mass $M_{X}$, together with the estimation of the $\theta_{X}$ for the maximal number of events given by (11) and the upper bound predicted by (13), is shown in Fig. 2.

Of course, it is not sufficient to use only the average energy $\langle E_{X}\rangle$ to estimate the position of the upper boundary. Indeed, the decrease of $c\tau_{X}$ with the growth of $\theta_{X}^{2}$ can be compensated by the increase of the energy $E_{X}$ and, therefore, of the $\gamma$-factor. As a result the particles with $E_{X}>\langle E_{X}\rangle$ can reach the detector even if the mixing angle $\theta_{X}$ is larger than the estimate (13).

The expression (13) helps to estimate how the sensitivity curve depends on the parameters of the experiment and on various assumptions. In particular, we can now estimate how large is a mistake from using $\langle E_{X}^{\rm rest}\rangle$ in Eq. (9) rather than the actual $E_{X}$ distribution. In order to do that we replaced $\langle E_{X}^{\rm rest}\rangle\to m_{\rm meson}$ – the maximal energy of the particle $X$ in the meson’s rest frame. This substitution increases the $\gamma_{X}$ by a factor of $2$. The estimates (11)–(13) show that $\theta_{\text{max}}^{2}$ and as a result $\theta_{\text{upper}}^{2}$ will shift by the same factor of $2$. This number indicates an upper bound on the possible error, introduced by the approximate treatment.

Next, we turn to the exact treatment. To this end we consider the energy distribution of the $X$ particles,

[TABLE]

Taking into account this distribution, the formula for the decay probability (3) at the upper bound should be modified as

[TABLE]

where an argument of $\pi$ function is $l_{{\text{\rm decay}}}/l_{\text{target-det}}$ and we used the expression for the decay length (4). The function $\pi(y)$, defined via

[TABLE]

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. $\pi(y)$ is shown in Fig. 3. For small energies (small $y$) $\pi(y)$ is exponentially small, while for large energies (large $y$) $\pi(y)$ is inversely proportional to energy and decreases slowly. Therefore, a sufficiently long “tail” of high-energy mesons can contribute to the integral in (15), but this range cannot be estimated without knowledge of the distribution function $f_{X}$. We will discuss $f_{X}$ for mesons and $W$ bosons in Sec. 5.

4 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 $\theta_{\text{lower}}$ (see Sec. 2), is equal to the distance from the target to the decay volume of the given experiment:

[TABLE]

The decay length (4) depends on the mass as $l_{{\text{\rm decay}}}\propto M_{X}^{-\alpha-1}$, where the term $\alpha$ in the exponent approximates the behaviour of the lifetime with the mass, and the term $1$ comes from the $\gamma$ factor.

Using the condition (17), the maximal mass probed can be estimated as

[TABLE]

which results in the following ratio of the maximal mass probed at the SHiP and MATHUSLA experiments:

[TABLE]

For Higgs-like scalars we have $\alpha\approx 2$, while for HNLs it is $\alpha\approx 5$, see Appendix B.2.

The estimate of the maximal mass probed (18) is applicable only if the result does not exceed the kinematic threshold; for the production from $B$ mesons for the HNLs it is $m_{B_{c}}-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_{\pi}$.

5 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 $\gamma$ factors (see Eqs. (8), (19)). 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).

5.1 $B$ and $D$ mesons

The main production channel of HNLs in the mass range $M_{N}\lesssim m_{D_{s}}$ is the two-body leptonic decay of $D_{s}$ mesons. For masses $m_{D_{s}}\lesssim M_{N}\lesssim m_{B_{c}}$ the main contribution comes from decays of $B$ mesons, see, e.g., Bondarenko:2018ptm .444This statement is true for HNLs with dominant mixing with $\nu_{e/\mu}$. For dominant mixing with $\nu_{\tau}$ the main production channel is from $\tau$ leptons for $M_{N}\lesssim m_{\tau}$ and from $B$ mesons for larger masses Bondarenko:2018ptm . For masses $M_{N}\gtrsim 3\text{ GeV}$ the main HNL production channel is determined by the value of the fragmentation fraction of $B_{c}$ mesons, $f_{B_{c}}$: in the case $f_{B_{c}}\gtrsim 10^{-4}$ 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 SHiP:2018xqw . For the scalars the production from $D$ mesons is negligible as compared to the $B^{+/0}$ mesons decays even for masses $m_{K}\lesssim M_{S}\lesssim m_{D}$. The $B_{c}$ mesons are not relevant for their production (see, e.g., Bezrukov:2009yw ; Boiarska:2019jym ). 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_{B_{c}}$ was measured at the LHCb Aaij:2017kea and found to be $f_{B_{c}}\approx(2.6\pm 1.3)\times 10^{-3}$. Earlier measurements at the Tevatron give a similar value $f_{B_{c}}\approx(2\pm 1)\times 10^{-3}$ Abe:1998wi ; Abe:1998fb ; Cheung:1999ir , which is in good agreement with Aaij:2017kea . Therefore, at the LHC the $B_{c}$ decay is the main production channel for heavy HNLs. However, at the energies of the SHiP experiment, $\sqrt{s}\simeq 30\text{ GeV}$, currently there is no experimental data on $f_{B_{c}}$. Additionally, the theoretical predictions of $f_{B_{c}}$ (see, e.g., Berezhnoy:2004gc ; Kolodziej:1997um ; Kolodziej:1995nv ) 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_{B_{c}}$ 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_{B_{c}}$ at the same level as at the LHC, and (ii) $f_{B_{c}}=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 CERN-SHiP-NOTE-2015-009 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_{B_{c}}=(m_{B_{c}}/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.

5.2 Mesons at MATHUSLA

In order to estimate the number of mesons and their $\gamma$ factors for the MATHUSLA experiment, one needs to know their $p_{T}$ distribution at ATLAS/CMS in the MATHUSLA pseudorapidity range $0.9<\eta<1.6$ (see Appendix D). The relevant distributions were measured for $B^{+}$ mesons by the CMS collaboration Khachatryan:2016csy (13 TeV) with the $p_{T}$ cut $p_{T}^{B}\geqslant 10\text{ GeV}$, and for $D^{+}/D^{0}$ mesons by the ATLAS collaboration Aad:2015zix (7 TeV) for $p_{T}^{D}\geqslant 3.5\text{ 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}\sim{}$few GeV (see, e.g., experimental papers Aaij:2013noa ; Khachatryan:2010us , theoretical paper Cacciari:2012ny 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 $\gamma$ factor, and therefore with the shortest decay length (4) and the largest probability to decay inside the decay volume (here we consider the case $l_{\text{\rm decay}}\gg l_{\text{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 $\gamma$ factor, and both of these effects enhance the number of events at the lower bound (8). Therefore an accurate prediction of the distribution $d\sigma/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}$ and also to estimate $D$ meson production cross-section at $\sqrt{s}=13$ TeV we use FONLL (Fixed Order + Next-to-Leading Logarithms) – a model for calculating the single inclusive heavy quark production cross section which convolutes perturbative cross section with non-perturbative fragmentation function, see Cacciari:1998it ; Cacciari:2001td ; Cacciari:2012ny ; Cacciari:2015fta for details.

Predictions of FONLL have been calibrated against the accelerator data and were found to be in very good agreement, see e.g. Aad:2015zix ; Khachatryan:2016csy ; Aad:2015gna ; Aaij:2015bpa ; Aaij:2013lla . In particular, comparison of the FONLL simulations of the production of the $B^{+}$s with the measurements at the Tevatron and at the LHC showed that FONLL predicts the low $p_{T}$ distribution accurately. We show the central values of the FONLL predictions down to $p_{T}=0$, confronted with the measurements of the CMS Khachatryan:2016csy and ATLAS Aad:2015zix 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 $\sigma_{D,\text{FONLL}}/\sigma_{D,\text{exp}}\approx 0.4$ and $\sigma_{B,\text{FONLL}}/\sigma_{B,\text{exp}}\approx 0.7$. However, as is demonstrated in the same papers Khachatryan:2016csy ; Aad:2015zix , 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:

[TABLE]

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 Aad:2015zix ; Khachatryan:2016csy .

FONLL does not provide the distributions of the $D_{s}$ and the $B_{c}$ mesons. We approximate their distributions by those 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 Chang:2003cq (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 $\eta$ distributions of $B_{c}$ and $B^{+}$ have similar shapes.

The relevant parameters — the total number of mesons, the average $\gamma$ factor of the mesons that are produced in the direction of the decay volume of the experiments and the geometric acceptances $\epsilon_{\text{geom},\text{meson}}$ for the mesons — are given in Table 1.

5.3 $W$ bosons

The production channel from the decays of $W$ bosons is only relevant for the MATHUSLA 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 $\sqrt{s}=13\text{ TeV}$ was measured in Aad:2016naf as $\sigma_{W\to N+l}\approx 20.5\text{ nb}$. The corresponding number of $W$ bosons produced during the high luminosity phase of the LHC is

[TABLE]

The $p_{T}$ distribution of the $W$ bosons at the LHC in the pseudorapidity range $|\eta|<2.5$ and for energies $\sqrt{s}=7-8\text{ TeV}$ was measured by the ATLAS and CMS collaborations Aad:2011fp ; Khachatryan:2016nbe . Their results show that most of the vector bosons are produced with low $p_{T}$ (below 10 GeV or so). However, these results do not give us the magnitude of the $W$’s average momentum $\langle p_{W}\rangle$, needed to estimate the decay acceptance and the average momentum of HNLs.

In order to obtain $\langle p_{W}\rangle$ we have simulated the process $p+p\to W^{\pm}$ in MadGraph5 Alwall:2014hca . In the leading order we have obtained $\sigma_{W\to\nu+l}\approx 15.7\text{ nb}$, which is in reasonable agreement with the prediction Aad:2016naf . The resulting momentum distribution of $W$ bosons is shown in Fig. 5 (left). A remark is in order here: at the leading order MadGraph5 does not predict the $p_{T}$ distribution of $W$s, since the production process is $2\to 1$ process and the colliding partons have $p_{T}=0$; 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 Aad:2011fp ; Khachatryan:2016nbe , 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 $\epsilon_{{\text{\rm decay}},W}$ and the average HNL momentum $\langle p_{X}\rangle$ by calculating the distribution of the particles in the energy $E_{X}$ and the angle $\theta_{X}$ between the direction of motion of the $X$ and the beam:

[TABLE]

Here $d^{2}\mathop{\text{BR}}\nolimits_{W\to X}/d\theta_{X}dE_{X}$ is the differential production branching ratio, and $P(\theta_{X})$ is a projector which takes the unit value if $\theta_{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

[TABLE]

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 (21) and the branching ratios of the scalar and HNL production from Appendix B.1. Therefore we conclude that for scalars the production from the $W$s is not relevant, while for HNLs careful estimation is needed.

In the case of HNL, the differential branching ratio in the Eq. (22) is

[TABLE]

The energy and angular distributions of the HNLs from the $W$ bosons at MATHUSLA are almost independent of the HNL mass in the mass range of interest, $M_{N}\ll 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\text{ GeV}$ is shown in Fig. 5. The decay acceptance was found to be $\epsilon_{{\text{\rm decay}},W}\simeq 2\%$, while the average momentum of the produced HNLs is $\langle p_{N}\rangle\approx 62\text{ 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 $\theta_{N}$ of the HNL, the energy distribution has a maximum at $E_{N,\text{max}}(\theta_{N})=m_{W}/2\sin(\theta_{N})$,555This formula is valid for 2-body decay into massless particles. which corresponds to HNLs produced from the $W$ bosons with some momentum $p_{W,\max}(\theta)$. As a consequence, the largest amount of HNLs flying in the direction $\theta_{N}$ has an energy close to their maximum, see the dashed line at the right panel of Fig. 5. The total energy spectrum is a superposition of different angles and has a peak at $E_{N,\text{peak}}\approx 58\text{ GeV}$ corresponding to the maximal angle possible at MATHUSLA, $\theta^{\max}_{\text{MAT}}\approx 44^{\circ}$. From the other side, the maximal energy possible for HNLs at MATHUSLA is defined by the minimal angle $\theta^{\text{min}}_{\text{MAT}}\approx 22^{\circ}$, which explains why the spectrum tends to zero near the energy $E_{N,\text{max}}\approx 106\text{ GeV}$.

6 Calculation of sensitivities

6.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

[TABLE]

where the quantities $\epsilon_{{\text{\rm decay}}},\epsilon_{\det},\mathop{\text{BR}}\nolimits_{\text{vis}}$ are the decay acceptance, detection efficiency and the visible branching correspondingly; they are defined by Eqs. (2), (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 (9) between $\gamma$ factors of the $X$ particle and the meson in the expressions (3), (15) for the decay probability.

Let us discuss the efficiencies. For the HNLs at SHiP experiment, we used the values of $\epsilon_{{\text{\rm decay}}}$ and $\epsilon_{\det}$ provided by detailed FairSHiP simulations zenodo . The results of the SHiP collaboration on the sensitivity to the scalars are not currently available, and for the product of $\epsilon_{{\text{\rm decay}}}\cdot\epsilon_{\det}$ we used the value for the HNL averaged over its mass, $\overline{\epsilon_{{\text{\rm decay}}}\cdot\epsilon_{\det}}\approx 0.2$.

For the MATHUSLA experiment there currently is no such detailed analysis of the efficiencies and background. In Curtin:2017izq ; Curtin:2018mvb 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 $\epsilon_{\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.666For the geometric acceptance as MATHUSLA we use the definition $\epsilon_{\text{geom}}=N_{\text{meson}}^{\eta\in\eta_{{\textsc{mat}}}}/N_{\text{meson}}\times\Delta\varphi/(2\pi)$, where $\eta_{{\textsc{mat}}}\in(0.9;1.6)$ and $\Delta\varphi=\pi/2$ are correspondingly pseudorapidity range of the MATHUSLA experiment and azimuthal size, see Sec. D. For the decay acceptance of the HNLs produced in the decays of the $W$ bosons we used the value $\epsilon_{{\text{\rm decay}},W}\approx 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 $\langle l_{\det}\rangle$ and $\langle l_{\text{target-det}}\rangle$ 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 needed parameter is the visible decay branching fraction. Following SHiP:2018xqw ; Curtin:2018mvb , 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 $\mathop{\text{BR}}\nolimits_{\text{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. 6.

6.2 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,777By the effective number of the mesons we call the production of the number of mesons multiplied by the overall efficiency, $N_{\text{meson}}\cdot\epsilon_{X,\text{meson}}$. the effective numbers of $B$ mesons are comparable between the experiments, and the average momenta (and therefore the $\gamma$ factors) of the mesons produced in the direction of the decay volume are $\simeq 7-8$ times smaller at MATHUSLA. The latter is caused by (i) different beam configurations (colliding beams for MATHUSLA, fixed target for SHiP) (ii) their different geometric orientation relative to the proton beam direction (the decay volume of the SHiP experiment is located in the forward direction, while the one of MATHUSLA’s is about $20^{\circ}$ off-axis.)

Using the numbers from the Tables 1, 2, for the ratio of the mixing angles at the lower bound (10) we have

[TABLE]

Qualitatively, for 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 $\simeq 5$ smaller than SHiP due to the smaller $\gamma$ factor of the $B$ mesons and larger effective number of $B$ mesons (i.e. the total number of $B$ mesons times the overall efficiency (25)). For the HNLs in the mass range $m_{K}\lesssim M_{N}\lesssim m_{D}$ the smallness of $\gamma$ 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 (26) was obtained under the optimistic condition $\epsilon_{\det}=1$ for MATHUSLA; after using a realistic efficiency the lower bound of the sensitivity at MATHUSLA will be changed by a factor $1/\sqrt{\epsilon_{\det}}$, which will affect the ratio (26).

6.3 Upper bound

We show the dependence of the number of events at $\theta^{2}_{X}=\theta^{2}_{\text{max}}$ as a function of the mass for the HNLs mixing with $\nu_{e}$ and the scalars in Fig. 7. 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 introduce the width of the upper bound defined by

[TABLE]

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}\gtrsim 2$ GeV (for smaller masses $\theta^{2}_{\text{upper}}$ lies deep inside the region excluded by the previous searches, see e.g. Bondarenko:2018ptm ). 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 (15). 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 (9) 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. 5. 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 $\theta^{2}_{\text{upper}}$ using the steepest descent method for the evaluation of the integral (15). The derivation of $\theta^{2}_{\text{upper}}$ is given in Appendix E.2.

Using $\theta^{2}_{\text{upper}}$, we present the upper bound width (27) in Fig. 8. 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. (13).

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 $\langle p_{N}\rangle$, 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 $\tau_{S}$ is changed with the mass slower than $\tau_{N}$, see the discussion in Sec. B.2. In addition, $\text{Br}_{B\to 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. 8 (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. 9. 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. (23)).

6.4 Maximal mass probed

The smaller $\gamma$ 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 $\nu_{e/\mu}$, the maximal mass probed ratio (19) becomes

[TABLE]

which agrees well with the sensitivity plot from Fig. 12. For the other cases — the HNLs mixing with $\nu_{\tau}$ and the scalars — the estimation of the maximal mass for the SHiP experiment based on the definition above exceeds the kinematic threshold, and therefore the result (19) 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 $\gamma$ factor.

7 Comparison with simulations

Next we compare our sensitivity estimates with the results of the SHiP and MATHUSLA collaborations. Different groups used different phenomenology for HNL and, especially, for scalars. Therefore, we will use different prescriptions for production and/or decay in different sections below, in order to facilitate the comparison of our approach with the Monte Carlo results of other groups. Our current view of the HNL phenomenology is summarized in Bondarenko:2018ptm and for scalar in Boiarska:2019jym . Our method of obtaining the sensitivity curves is summarized in Appendix F.

7.1 HNLs

The results for the HNLs are shown in Fig. 10. 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/\langle p_{\text{meson}}\rangle$ and $\langle 1/p_{\text{meson}}\rangle$, 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 $\gamma$ factors used in our estimation and those obtained in Monte Carlo simulations by the SHiP collaboration. Indeed, using the SHiP simulations results available in SHiP:2018xqw ; zenodo , we have found that for the masses $M_{N}\simeq m_{B_{c}}$ the ratio of average $\gamma$ factors is $\langle\gamma_{N}^{\text{analytic}}\rangle/\langle\gamma_{N}^{\text{simulations}}\rangle\simeq 0.8$, which seems to explain the difference.

For the sensitivity of MATHUSLA Curtin:2018mvb to the HNLs produced in $W$ decays 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}\gtrsim m_{D_{s}}$, 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.

7.2 Scalars

The comparison of our sensitivity estimates with the results of the SHiP and MATHUSLA experiments is presented in Fig. B.3. We also show the results of a simple analytic estimate of the lower bound for particular masses from Appendix G. For the comparison with the sensitivity provided by the MATHUSLA collaboration we used the model of scalar production and decay given in Curtin:2018mvb , while comparing with the results of the SHiP collaboration – from Beacham:2019nyx . A description of the models is given in Appendix B.3.

The sensitivity curves are in good agreement. Small differences in the position of the maximal mass probed can be explained by different energy distributions of the scalars used in our estimate and in Curtin:2018mvb and in Beacham:2019nyx .

8 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 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 SHiP:2018xqw and MATHUSLA Curtin:2018mvb 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. 12 and to the dark scalars in Fig. 13.

Qualitatively both experiments can probe similar ranges of parameters. The SHiP has higher average $\gamma$ factors of the mesons ($\langle\gamma_{\text{meson}}^{\textsc{ship}}\rangle/\langle\gamma_{\text{meson}}^{{\textsc{mat}}}\rangle\simeq\mathcal{O}(10)$) and, as a result, significantly higher upper boundary of the sensitivity region than MATHUSLA (as the upper boundary is exponentially sensitive to the $\gamma$ factor). As the consequence, the SHiP can probe higher masses for both HNLs and scalars than MATHUSLA (except of HNLs with dominant mixing with tau flavor). 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}\lesssim m_{D_{s}}$ owing to the larger number of $D$ mesons. MATHUSLA can probe lower mixing angles for the HNLs with $M_{N}\gtrsim m_{D_{s}}$ 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\sigma/dp_{T}$ distribution of $B$ mesons at the LHC has a maximum at $p_{T}\sim\mathrm{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\sigma/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 Anelli:2015pba ; SHIP:2018yqc . For MATHUSLA 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 Curtin:2017izq ; Curtin:2018mvb . It is not known how much this rejection affects the detection efficiency, $\epsilon_{\det}$. In this work, we conservatively assumed $\epsilon_{\det}=1$ for MATHUSLA, while for SHiP it was taken from the actual Monte Carlo simulations SHiP:2018xqw . 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 $\sqrt{s}\approx 28$ GeV. It changes the position of the lower bound and consequently the maximal mass probed in a significant way, see Fig. 12.

Comparison with other works.

We have compared our sensitivity estimates with the results of the Monte Carlo simulations presented by the SHiP collaboration SHiP:2018xqw ; Beacham:2019nyx ; SHiP:2019 and with the estimates of the MATHUSLA physics case paper Curtin:2018mvb (Figs. 10–11). 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 Curtin:2018mvb , 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.

Acknowledgements.

We thank D. Curtin, J. Evans, R. Jacobsson and W. Valkenburg for fruitful discussions and comments on the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (GA 694896) and from the Netherlands Science Foundation (NWO/OCW).

Appendix A Portals

New particles with masses much lighter than the electroweak scale can couple to the Standard Model fields via renormalizable interactions with small dimensionless coupling constants (sometimes called “portals” as they can mediate interactions between the Standard Model and “hidden sectors”). In this work, we considered two renormalizable portals: scalar (or Higgs) portal and neutrino portal.

The scalar portal couples a gauge-singlet scalar $S$ to the gauge invariant combination $H^{\dagger}H$ made of the Higgs doublet:

[TABLE]

where $g$ is the coupling constant and $\mathcal{L}_{\rm int}$ are interaction terms that play no role in our analysis. After the spontaneous symmetry breaking the cubic term in (29) gives rise to the Higgs-like interaction of the scalar $S$ with all massive particles with their mass times a small mixing parameter

[TABLE]

where $g$ is the coupling in (29); $v$ is the Higgs VEV; $m_{H}$ is the Higgs mass; sum in (30) goes over all massive fermions (leptons and quarks); $W^{\pm}_{\mu}$ is the $W$ boson and $\cdots$ denote other interaction terms, not relevant for this work. The details of the phenomenology of the scalar portal are provided in Boiarska:2019jym (see also Patt:2006fw ; Batell:2009di ; Clarke:2013aya ; Schmidt-Hoberg:2013hba ). The computation of hadronic decay width of $S$ is subject to large uncertainties at masses $M_{S}\sim{}$few GeV, where neither chiral perturbation theory not perturbative QCD can provide reliable results (see a discussion in Monin:2018lee ).

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 $\epsilon_{ab}\bar{L}_{a}H_{b}$ where $L_{a}$ is the SU(2) lepton doublet and $\epsilon_{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_{\alpha}\ll 1$ ($\alpha=\{e,\mu,\tau\}$).

Appendix 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 $q\bar{q}$ pairs and fragmentation fractions $f_{\text{meson}}$, that can be extracted from the experimental data Aaij:2011jp ; Aaij:2017kea ; Kling:2018wct . 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 Graverini:2133817 . 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.

B.1.1 HNL production

The production of the HNL in the decay of charmed and beauty mesons have been considered in Gorbunov:2007ak ; Atre:2009rg , see Bondarenko:2018ptm 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. 14. We see that for the HNL mass range $m_{N}\gtrsim 3.5$ GeV the main production channel is $B_{c}$ meson decay $B_{c}\to 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

[TABLE]

where we take $m_{N}\gg m_{\ell}$, and $K$ is a kinematic suppression. Neglecting them, for the ratio for the numbers of HNLs produced by $B_{c}$ and $B^{+}$ we obtain

[TABLE]

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\to N+l$. The corresponding branching ratio is

[TABLE]

where we have neglected the HNL and the lepton masses.

B.1.2 Quarkonia and heavy flavour baryons

Quarkonia states (especially $\Upsilon$ meson) can produce HNLs reaching $10$ GeV in mass above the beauty meson threshold. The contribution of quarkonia decays to the production were found negligible at SHiP energies, see Bondarenko:2018ptm .

The LHC experiments have measured $\Upsilon$ production at both ATLAS and CMS Aad:2012dlq ; Chatrchyan:2013yna . In the rapidity range $|y|<2$, relevant for MATHUSLA, the cross-section is given by Hu:2017pat

[TABLE]

(as this is an order of magnitude estimate, we combine production of $1S$, $2S$ and $3S$ bottomonium states and neglected both statistical and systematic uncertainties of the cross-section measurement). Using $\mathop{\text{BR}}\nolimits(\Upsilon\to\mu^{+}\mu^{-})\simeq 2.4\times 10^{-2}$ Tanabashi:2018oca we find that during the high luminosity phase one can expect $N_{\Upsilon}\sim 10^{12}$. Large fraction of this mesons are traveling into the direction of the fiducial decay volume of MATHUSLA, as their distribution is sufficiently flat in the $|y|<2$ rapidity range. This number should be multiplied by the branching ration $\mathop{\text{BR}}\nolimits(\Upsilon\to N\nu)$, estimated in Bondarenko:2018ptm to be at the level $\mathop{\text{BR}}\nolimits(\Upsilon\to N\nu)\sim 10^{-5}U^{2}$, so that overall one expects in MATHUSLA detector about $10^{7}U^{2}$ HNLs from $\Upsilon$ decays.

This number should be compared with those, produced from $W$-bosons (as we are above the $B$-meson threshold): $N_{W}\times\mathop{\text{BR}}\nolimits(W\to N+l)\times\epsilon_{N}$, where $N_{W}$ is given by (21), $\epsilon_{N}\approx 0.02$ is the geometric acceptance for the HNLs produced from $W$ and flying into the MATHUSLA fiducial volume and the branching fraction is given by (33). The resulting number is $\sim 6\times 10^{8}U^{2}$ – exceeding the number of HNLs from $\Upsilon$-mesons by about 2 orders of magnitude.

As Table 3 demonstrates, about $25\%$ of $b$-quarks at the LHC hadronize into the $\Lambda_{b}^{0}$ baryons. These baryons produce HNLs in the 3-body semi-leptonic decay $\Lambda_{b}^{0}\to B+\ell+N$ where $B$ is a baryon. The mass of the $\Lambda_{b}^{0}$ is $m_{\Lambda_{b}^{0}}\simeq 5.62\,\mathrm{GeV}$. The decays $\Lambda_{b}^{0}\to p+\ell^{-}+N$ are suppressed by the CKM matrix element $V_{bu}$, while the decays $\Lambda_{b}^{0}\to\Lambda_{c}^{+}+\ell^{-}+N$ can only produce HNLs with $M_{N}<m_{\Lambda_{b}^{0}}-m_{\Lambda_{c}^{+}}\simeq 3.35\,\mathrm{GeV}$. HNLs of this mass are produced from more copious $B$-mesons and therefore $\Lambda$ baryons can be neglected.

The contribution of heavy flavour baryon decays to the production were found negligible at SHiP energies, see Bondarenko:2018ptm .

B.1.3 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 Bezrukov:2009yw ; Boiarska:2019jym . The main production process is the 2-body decay

[TABLE]

where $X_{q}$ is a hadron that contains the quark $q$. The branching ratios for these states are discussed in Boiarska:2019jym . Here we only state the main points. For $B\to 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 Tanabashi:2018oca ). We have found that each heavier meson from the “family” gives a smaller contribution to the branching ratio than the lighter one. For $B\to X_{d}+S$ we take only one meson $X_{d}=\pi$ since this channel has the largest kinematic threshold $m_{B}-m_{\pi}$. We summarize the list of the final states below:

• –

Spin 0, odd parity: $X_{q}=\pi,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. Cheng:2013fba ). 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 (35) on the scalar mass is shown in Fig. 15.

We estimate the production of the scalars from the $W$ bosons by the decay $W\to S+f+\bar{f^{\prime}}$, where the summation over all the SM fermions species $f=l,q$ is taken. We obtained $\mathop{\text{BR}}\nolimits_{W\to S}/\theta^{2}\simeq 4\cdot 10^{-3}$.

We mention in passing that the production of scalars from $\Upsilon$ (due to $b\to s$ transition) is not playing essential role, as the mass difference $m_{\Upsilon}-m_{B}<m_{B}$ and therefore one should compare the number of scalars produced from the bottomonium decays with the number of scalars from $B$-meson decays. The latter of $B$-mesons is several orders of magnitude higher (see Table 1). In addition to that the branching ratio of $\Upsilon\to B+S$ is much smaller than $B\to K+S$ because the width of $\Upsilon$ is dominated by electromagnetic decays.

The $b\to s$ transitions also generate decays $\Lambda_{b}^{0}\to\Lambda^{0}+S$. However, the mass of thus produced scalar, $M_{S}<m_{\Lambda_{b}^{0}}-m_{\Lambda^{0}}\simeq 4.5\,\mathrm{GeV}$ and thus is subdominant to the production from $B$ mesons.

B.2 Main decay channels

B.2.1 HNL

The HNL has 3-body leptonic decays and different semileptonic modes. Following the paper Bondarenko:2018ptm , 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 decay width of HNLs mixing with the flavor $\alpha$ can be approximated by the formula

[TABLE]

where $g_{\text{eff}}$ is a dimensionless factor that depends on the mass of HNL and changes from 1 to $\sim 10$, see e.g. Bondarenko:2018ptm for details.

The dependence of the proper lifetime $c\tau_{N}$ on the HNL mass at $U^{2}=1$ is given at the left panel of Fig. 16.

B.2.2 Scalar

The decay width of the scalar particle has large uncertainty in the scalar mass region $0.5\text{ GeV}<M_{S}<2\text{ GeV}$ because of resonant nature of $S\to 2\pi$ decay, see Monin:2018lee for the recent overview. At higher masses the decay width is determined by perturbative QCD calculations Spira:1997dg . We omit the problem of pion resonance in this work using continuous interpolation between the sum of the decay channel at low masses and perturbative QCD at high ones.

For scalar mass region above 2 GeV one can naively estimate $S$ decay width as $\Gamma_{S}\propto\sum_{f}\theta^{2}y_{f}^{2}M_{S}$. This estimation does not take into account three effects:

For the decay into quarks parameter $y_{q}$ depends on scalar mass as $y_{q}\equiv M_{q}(M_{S}^{2})/v$, where $M_{q}(M_{S}^{2})$ is quark running mass, which gives logarithmic correction; 2. 2.

The decay into gluons has different $M_{S}$ dependence, $\Gamma_{S}\propto\theta^{2}M_{S}^{3}/v^{2}$, and dominates in the region $2\text{ GeV}<M_{S}<3.5\text{ GeV}$ Bezrukov:2009yw ; 3. 3.

In the region $M_{S}$ near $3.5\text{ GeV}$ new decay channels appear (into $\tau$ and $c$ quark), and the kinematical factor is important.

Taking them into account, for the mass domain $3.5\text{ GeV}<M_{S}<5\text{ GeV}$, near the threshold of production from $B$ mesons, we made a fit to the total $\Gamma_{S}$ and found that its behavior is $\Gamma_{S}\propto M_{S}^{2}$.

The dependence of the proper lifetime $c\tau_{S}$ on the scalar mass for $\theta^{2}=1$ is shown in Fig. 16 (right panel).

B.2.3 Visible branching ratio

We define the “visible” decay channels as those that contain at least two charged particles $\alpha$ in the final state. The corresponding decays are

[TABLE]

where $Y$ is arbitrary state, $F$ is uncharged state that decays to $n$ charged particles and $\tilde{Y}$ is a state with at least $2-n$ charged particles (assuming $n<2$). Using this definition, the decay $N\to 3\nu$ is identified as invisible decay, the decay $N\to\bar{\mu}\nu_{\mu}e$ — as visible decay, while the decay $N\to\eta\nu$ — as the visible decay if $\eta$ meson decays into two charged particles, and as invisible decay otherwise. To take into account only visible decays of $F$, we include the factor $\mathop{\text{BR}}\nolimits_{F\to\text{vis}}$ to the partial decay width $\Gamma_{X\to F\tilde{Y}}$. We take the values of $\mathop{\text{BR}}\nolimits_{F\to\text{vis}}$ from Tanabashi:2018oca . 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.

B.3 Comparison with scalar models used by SHiP and MATHUSLA collaborations

Our model of the scalar production and decay described above differs from those used in Curtin:2018mvb ; Beacham:2019nyx for estimating the sensitivity. Namely, for scalar production in Ref. Boiarska:2019jym has summed over main exclusive channels $B\to X_{s/d}+S$.888There is $30\%$ level uncertainty in the total production rate because of the $B\to 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. Cheng:2013fba 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. In Ref. Curtin:2018mvb the production from $B$ mesons is estimated using the free quark model, while Beacham:2019nyx considers only the channel $B\to KS$. This causes differences in the magnitude of the branching ratio and kinematic production thresholds. In particular, we note that the free quark model breaks down for large scalar masses $M_{S}\simeq 3\text{ GeV}$ since the QCD enters the non-perturbative regime, and therefore it gives meaningless predictions for the production rate of heavy scalars.

As for the scalar decay width, because of theoretical uncertainty for the mass range $2m_{\pi}\lesssim M_{S}\lesssim 2m_{D}$ there is no agreement in the literature how to describe the scalar decays in this domain, see Monin:2018lee . Our decay width differs significantly from the decay width used in Curtin:2018mvb ; Beacham:2019nyx .

Appendix 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_{X}^{\text{rest}}$, and meson energy $E_{\text{meson}}$ at lab frame as

[TABLE]

where $\theta$ 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

[TABLE]

where $\gamma_{\text{meson}}\equiv E_{\text{meson}}/m_{\text{meson}}$.

Appendix D Geometry of the experiments

D.1 SHiP

The SHiP experiment SHIP:2018yqc is a fixed-target experiment using the proton beam of the Super Proton Synchrotron (SPS) at CERN. SPS can deliver $N_{\text{p.o.t.}}=2\times 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 $\mathcal{O}(10)$ GeV and $c\tau$ exceeding tens of meters. The overview of the experiment is as follows.

The proton beam hits a target Alekhin:2015byh ; Anelli:2015pba . The target will be followed by a $5\leavevmode\nobreak\ \text{m}$ hadron stopper, intended to stop all $\pi^{\pm}$ 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\leavevmode\nobreak\ \text{m}$ long.

The decay volume is a long pyramidal frustum vacuum chamber with the length

[TABLE]

and the cross-section $5\text{ m}\times 10\text{ m}$. It begins at

[TABLE]

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 $\theta\approx 25\text{ mrad}$.

D.2 MATHUSLA

MATHUSLA (MAssive Timing Hodoscope for Ultra Stable neutraL pArticles) is a proposed experiment Evans:2017lvd ; Curtin:2017izq that consists of a $20\text{ m}\times 200\text{ m}\times 200\text{ m}$ surface detector, installed above ATLAS or CMS detectors (see Fig. 17). The long-lived particles, created at the LHC collisions, travel 100+ meters of rock and decay within a large decay volume ($8\times 10^{5}\text{ 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_{L}^{0}$ 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 to

[TABLE]

where $L_{\text{ground}}=100\text{ m}$. The average distance a particle travels inside the decay volume, $\bar{L}_{\text{det}}$, is given by

[TABLE]

Geometrical parameters of MATHUSLA experiment are summarized in Table 4.

Appendix E Analytic estimation of the upper bound: details

In this Section we estimate the ratio between $\theta_{\text{max}}$ and $\theta_{\text{upper}}$ – the quantities defined in Section 3.

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. 18:

[TABLE]

The distribution of the high energy $B$ mesons at the LHC for energies $E_{B}\lesssim 300\text{ GeV}$ can be approximated by the power law function, see the right panel in Fig. 18:

[TABLE]

Finally, the distribution of the HNLs originating from the $W$ bosons can be approximated by the expression

[TABLE]

E.2 Upper bound estimation

We start from the number of the events

[TABLE]

Here for simplicity we defined quantities $\tilde{P}_{{\text{\rm decay}}}$ and $\tilde{N}_{{\text{\rm decay}}}$ defined by

[TABLE]

The decay probability (15) can be rewritten in the form

[TABLE]

where $g(E)=l_{\rm target-det}\Gamma_{X}M_{X}/E-\log\left(\frac{dN_{X}}{dE}\right)$ (for clarity we assumed $\frac{dN_{X}}{dE}$ to be dimensionless). The integral (15) can be evaluated as

[TABLE]

where $E_{\mathrm{peak}}$ is the minimum of $g(E)=0$, here we used the steepest descent approximation.

For the exponential spectrum (44) the peak energy and the probability are, correspondingly,

[TABLE]

while for the power law spectrum $f_{E}=f_{0}E^{-\alpha}$ they are

[TABLE]

Expressing then $\Gamma_{X}\propto U^{2}$ and $\theta^{2}_{X}\equiv\theta^{2}_{\text{max}}(M_{X})\times R$, for the upper bound given by the particles produced from $B$ mesons one immediately arrives to

[TABLE]

Similarly we can estimate the upper bound for the $W$ bosons.

We estimate $\theta_{\text{max}}$ for the production from $B$ mesons using the HNLs momenta given by Eq. (9) with $\langle p_{B}\rangle=12\text{ GeV}$ for MATHUSLA and $\langle p_{B}\rangle=88\text{ GeV}$ for SHiP (see Table 2), while for the production for $W$ we use $\langle p_{N}\rangle\approx 62\text{ GeV}$, see Sec. 5.3.

Appendix F Details of the sensitivity curve drawing

We draw the sensitivity curve for HNLs and scalars requiring

[TABLE]

where numbers of decay events of particles produced from $B,D$ mesons and $W$ bosons are estimated as

[TABLE]

[TABLE]

Here, $N_{q\bar{q}}$ is the total number of the $q\bar{q}$ pairs that are produced in $pp$ collisions at the high luminosity LHC or in $p-\text{target}$ collisions at SHiP, $f_{p_{\text{meson}}}$ is the momentum distribution of the mesons that fly to the decay volume of the experiment (see Sec. 5.1), and $\epsilon$ is the overall efficiency (see Sec. 6.1). In the expression for the decay probability $\tilde{P}_{{\text{\rm decay}}}$ the $\gamma$ factor of the $X$ particle is related to the meson momentum by the relation (39). $N_{W,\text{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).

Appendix 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 (8). The relevant parameters and the value $\theta^{2}_{\text{lower}}$ for the SHiP and MATHUSLA experiments are given in Tables 5, 6.

Appendix H HNLs at MATHUSLA for small mass

For HNLs with $M_{N}\lesssim m_{D_{s}}$, where the sensitivity curve is determined by the production from $D$ mesons our sensitivity curve reproduces that of the MATHUSLA collaboration Curtin:2018mvb in the range up to 1 GeV (Fig. 10). For smaller masses our estimates of the lower boundary differ by a factor $\sim 3$ (which would corresponds to the order-of-magnitude difference between the number of decay events). Moreover, the shapes of the sensitivity curves also differ.

Below we list several possible reasons for this discrepancy:

1. a)

different estimate of the number of the parent $D$ mesons produced 2. b)

HNLs that are produced from the mesons that do not fly into the fiducial decay volume of MATHUSLA (they were not taken into account in our estimate) 3. c)

different HNL phenomenology (production and decay) used in comparison. 4. d)

Production from $K$-mesons was not taken into account in our estimates

Since the positions of the lower bounds in the mass range $1\,\mathrm{GeV}\lesssim M_{N}\lesssim m_{D_{s}}$ 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}\lesssim m_{D}$.

In order to estimate the amount of light HNLs produced from $D$ mesons, we performed a MadGraph 5 simulation of a process $pp\to\bar{c}e^{+}\nu_{e}s$, whose kinematics corresponds to the main process for a production of light HNLs $M_{N}\lesssim 0.5\,\mathrm{GeV}$ — $D\to Ne^{+}K$ Bondarenko:2018ptm . We computed the ratio

[TABLE]

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 $c\bar{c}$ pairs that fly in the same direction. We found $\chi\approx 1.7$, which is not enough to explain the discrepancy.

For c), we compared the decay widths of the HNLs used in our analysis with those, used in Curtin:2018mvb (based on Helo:2010cw ). We found them to be different by $20-40\%$ (for $U_{e}:U_{\mu}:U_{\tau}=1$) with the decay width from Helo:2010cw being smaller than that from Bondarenko:2018ptm . The difference can reach up to 80% at small masses (below $\mathcal{O}(500)$ MeV).

Finally, production from $K$ mesons would not explain why the discrepancy starts close to 1 GeV, much higher than production threshold from kaons.

We did not find the information about the HNL production ratios used in Curtin:2018mvb . As we see, the cases b) – d) 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 adopted in  Curtin:2018mvb .