aainstitutetext: Institut für Astroteilchen Physik, Karlsruher Institut für
Technologie (KIT), Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-
Leopoldshafen, Germanybbinstitutetext: Instituut-Lorentz, Leiden University,
Niels Bohrweg 2, 2333 CA Leiden, The Netherlandsccinstitutetext: Departamento
de Física Teórica and Instituto de Física Teórica UAM/CSIC, Universidad
Autónoma de Madrid, Cantoblanco, 28049, Madrid, Spainddinstitutetext: IFPU,
Institute for Fundamental Physics of the Universe, via Beirut 2, I-34014
Trieste, Italyeeinstitutetext: SISSA, via Bonomea 265, I-34132 Trieste,
Italyffinstitutetext: INFN, Sezione di Trieste, SISSA, Via Bonomea 265, 34136,
Trieste, Italy

# Sensitivities to feebly interacting particles: public and unified
calculations

Maksym Ovchynnikov c Jean-Loup Tastet b Oleksii Mikulenko d,e,f Kyrylo
Bondarenko<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The idea that new physics could take the form of feebly interacting particles
(FIPs) — particles with a mass below the electroweak scale, but which may have
evaded detection due to their tiny couplings or very long lifetime — has
gained a lot of traction in the last decade, and numerous experiments have
been proposed to search for such particles. It is important, and now very
timely, to consistently compare the potential of these experiments for
exploring the parameter space of various well-motivated FIPs. The present
paper addresses this pressing issue by presenting an open-source tool to
estimate the sensitivity of many experiments — located at Fermilab or at the
CERN’s SPS, LHC, and FCC-hh — to various models of FIPs in a unified way: the
Mathematica-based code SensCalc.

##  1 Introduction

The well-known shortcomings of the Standard Model suggest to us the existence
of new physics “Beyond the Standard Model” (BSM), that is generally expected
to involve new particles. There is currently no clear theoretical guidance,
nor experimental hints, about the mass of the hypothetical new particles,
which could range from sub-eV all the way up to the Planck scale. Particles
with a mass below the electroweak scale are of particular interest
experimentally, since they may be numerously produced at accelerators. Past
experiments have already excluded the largest values of the couplings for such
particles; hence they are called feebly interacting particles, or FIPs for
short. FIPs may be searched for both at the main detectors of colliders
(ATLAS, CMS and LHCb at the LHC, or their equivalents at future colliders such
as the FCC-hh) which are located very close to the collision point, or at so-
called lifetime-frontier experiments, which re-use existing facilities or
infrastructure and place a displaced decay volume near an interaction point or
target. Lifetime-frontier experiments may be broadly split into two classes
Beacham:2019nyx : collider-based, which consist of a decay volume placed near
the interaction points of ATLAS, CMS, and LHCb, and extracted-beam
experiments, which use an extracted beam line hitting a target.

During the last few years, many lifetime-frontier experiments have been
proposed. Among extracted-beam experiments, we can list SHiP Alekhin:2015byh ;
Aberle:2839677 , SHADOWS Alviggi:2839484 , and HIKE CortinaGil:2839661  at the
SPS, and DUNE DUNE:2020ypp ; DUNE:2020fgq  and DarkQuest Batell:2020vqn  at
Fermilab. The proposed LHC-based experiments include MATHUSLA MATHUSLA:2019qpy
and FACET Cerci:2021nlb , associated with CMS; FASER FASER:2018bac , SND@LHC
SHiP:2020sos  (together with their upgrades, AdvSND and FASER2) and ANUBIS
Bauer:2019vqk , close to the ATLAS interaction point; Codex-b Aielli:2019ivi
near LHCb; and AL3X Dercks:2018wum  at ALICE. Furthermore, lifetime-frontier
experiments will likely remain part of the physics program of future
colliders, such as the FCC-hh Boyarsky:2022epg .

In order to evaluate the potential of those experiments to search for generic
FIPs, the PBC initiative has proposed Beacham:2019nyx  a few benchmark models.
They include dark photons, millicharged particles, dark scalars, heavy neutral
leptons, and axion-like particles coupled to various SM particles.

While some of the experiments from the above list are already running, many
are still at the status of proposals. Their design is not finalized yet, and
is still undergoing optimization. Their sensitivity can be optimized by
focusing on two key aspects: increasing the rate of events that contain FIPs,
and reducing the Standard Model (SM) backgrounds. Studying the background
requires knowing the detailed specifications of the experimental setup,
background-reducing systems, and surrounding infrastructure. As a result, full
simulations are required, which accurately trace each event starting from the
initial proton collision and ending with the interactions of the background
particles with the detector material. Most of the experimental proposals claim
to achieve zero background level. In contrast, the evaluation of the FIP event
rate is comparatively less affected by these complexities. This is the case,
in particular, when the FIPs are produced at the collision point. They would
then propagate through the infrastructure without being affected (due to their
tiny interaction strength), and decay or scatter inside the decay volume with
some tiny probability. If the reaction products reach the detector and satisfy
some simple kinematic cuts, they could typically be detected with $\approx 1$
efficiency. Therefore, the sensitivity111When talking about “sensitivity”, it
is important to point out the distinction between “exclusion” sensitivity
(rejecting the New Physics hypothesis in the absence of signal) and
“discovery” sensitivity (rejecting the Standard Model in favor of New Physics
if a signal is observed). While the former is not very sensitive to the exact
background expectation as long as it is $\lesssim 1$, the latter strongly
depends on it. Throughout this paper, we mean “exclusion sensitivity” whenever
we use the word “sensitivity” unqualified. However, we advise the reader to
keep this distinction in mind when comparing the physics potential of various
experiments: indeed, two experiments with the same exclusion sensitivity could
in principle have significantly different discovery sensitivities. of a given
experiment to FIPs is determined mainly by 1) the distribution of FIPs at the
facility housing the experiment and 2) the geometry of the experiment itself.

Despite the relative simplicity of estimating the sensitivity to FIPs, there
exists a caveat that can make their comparison challenging: the sensitivity
estimates performed by the various experimental collaborations are not
publicly accessible. This is crucial for three reasons. First, there is often
no unique description of the production and decay of a given FIP in the
literature. This is related to either theoretical uncertainties in the
description of the FIP phenomenology, or different conventions in the
definition of the model. As a result, different collaborations can end up
using different FIP descriptions; sometimes, even the definition of the FIP
coupling is different (see Appendix A). Secondly, due to the rapid pace of
change as the experiment’s design is being optimized, there may exist a
mismatch between, on the one hand, the experimental setup and/or the
assumptions used and, on the other hand, the reported sensitivity, even within
a same document (see Fig. 4 and the corresponding discussion). Indeed, in
order to update the sensitivity while the setup is undergoing optimization,
collaborations would need to re-launch full-scale simulations, which require a
lot of time, computational resources and person-power. Finally, these
calculations are black-box: they do not provide a qualitative understanding of
the sensitivity. This problem becomes especially important when comparing the
sensitivities of various experiments to understand which one is better suited
to probe a given region of the FIP parameter space.

To address these issues, a public tool that can calculate the sensitivity of
various experiments to FIPs — in a unified and transparent way — is required.
Several publicly available packages can already perform such sensitivity
calculations Kling:2021fwx ; Jerhot:2022chi . However, they are limited to a
specific type of facilities: either beam dump experiments or colliders. This
paper presents the Mathematica Mathematica  code SensCalc SensCalc-Zenodo ,
that can evaluate the sensitivity of the various experiments proposed at LBNE,
SPS, LHC, and FCC-hh to various FIPs.222Available at
https://doi.org/10.5281/zenodo.7957785 . The code is based on a semi-analytic
approach developed in Ref. Bondarenko:2019yob  and further improved and cross-
checked in Refs. Boiarska:2021yho ; Boyarsky:2022epg ; Ovchynnikov:2022its
(see also Coloma:2023adi ; Batell:2023mdn ), where the number of events is
approximated by the integral of several quantities: the FIP angle-energy
distribution, its decay probability, the geometric acceptance of FIPs, and the
acceptance of its decay products. Most of these quantities can be accurately
computed analytically, which is especially attractive as it improves the
transparency of the computations.

The present paper is organized as follows. In Sec. 2, we discuss the semi-
analytic method that we use for calculating the sensitivity, together with its
validation and limitations. In Sec. 3, we provide a brief description of
SensCalc, specifying the list of currently implemented experiments and models
of FIPs. We also compare it with other publicly available packages for
computing the sensitivity, as well as with SensMC SensMC-GitHub , a simplified
Monte-Carlo simulation that we have specifically developed to validate
SensCalc. Finally, we conclude in Sec. 4.

##  2 Semi-analytic approach to calculate sensitivities

###  2.1 Method

This work concentrates on FIPs produced at the collision point (or close to
it). In this case, the production is unaffected by the surrounding
infrastructure. We calculate the number of events involving a decaying FIP
using the following expression:

$N_{\text{ev}}=\sum_{i}N^{(i)}_{\text{prod}}\int dEd\theta dz\
f^{(i)}(\theta,E)\cdot\epsilon_{\text{az}}(\theta,z)\cdot\frac{dP_{\text{dec}}}{dz}\cdot\epsilon_{\text{dec}}(m,\theta,E,z)\cdot\epsilon_{\text{rec}}$
(1)

Figure 1: Illustration of the impact of different contributions to the number
of events (1), taking as an example a beam dump experiment with a detector
located downstream of the decay volume. Consider a FIP decaying at coordinates
$(\theta,z)$, where $\theta$ is the polar angle relative to the beamline, $z$
is the longitudinal displacement from the target, and the azimuthal angle
$\phi$ has been omitted from the diagram. The differential probability for a
FIP with energy $E$ to decay there is $f(\theta,E)dP_{\text{dec}}/dz$. The
azimuthal coordinate $\phi$ of the decaying FIP (whose trajectory is shown by
the red arrow) must be within the decay volume, which restricts the available
decay positions to the blue dashed line. These limitations are included in the
azimuthal acceptance $\epsilon_{\text{az}}$. Next, at least two of the FIP
decay products (the green arrows) have to point to the detector; this is
accounted for in $\epsilon_{\text{dec}}$. Depending on the setup and the FIP,
this requirement may significantly limit the decay volume’s “useful” angular
coverage. In particular, for 2-body decays into stable particles, the decay
products can only point to the detector if the decayed FIP also points to the
detector. Only the narrow angular domain that the detector covers contributes
to the number of events.

The quantities entering Eq. (1) are the following (see also Fig. 1):

  * –

$N_{\text{prod}}^{(i)}$ is the total number of FIPs produced by the process
$i$, e.g., decays of mesons, direct production by proton-target collisions,
etc. (see fig. 2).

  * –

$z$, $\theta$ and $E$ are, respectively, the position along the beam axis, the
polar angle, and the energy of the FIP.

  * –

$f^{(i)}(\theta,E)$ is the differential distribution of FIPs in polar angle
and energy, for FIPs produced through the process $i$.

  * –

$\epsilon_{\text{az}}(\theta,z)$ is the azimuthal acceptance:

$\epsilon_{\text{az}}=\frac{\Delta\phi_{\text{decay volume}}(\theta,z)}{2\pi}$
(2)

where $\Delta\phi$ is the fraction of azimuthal coverage for which FIPs
decaying at $(z,\theta)$ are inside the decay volume.

  * –

$\frac{dP_{\text{dec}}}{dz}$ is the differential decay probability:

$\frac{dP_{\text{dec}}}{dz}=\frac{\exp[-r(z,\theta)/l_{\text{dec}}]}{l_{\text{dec}}}\frac{dr(z,\theta)}{dz},$
(3)

with $r=z/\cos(\theta)$ being the modulus of the displacement of the FIP decay
position from its production point, and
$l_{\text{dec}}=c\tau\sqrt{\gamma^{2}-1}$ is the FIP decay length in the lab
frame.

  * –

$\epsilon_{\text{dec}}(m,\theta,E,z)$ is the decay products acceptance, i.e.
among those FIPs that are within the azimuthal acceptance, the fraction of
FIPs that have at least two decay products that point to the detector and that
may be reconstructed. Schematically,

$\epsilon_{\text{dec}}=\text{Br}_{\text{vis}}(m)\times\epsilon_{\text{dec}}^{\text{(geom)}}\times\epsilon_{\text{dec}}^{\text{(other
cuts)}}$ (4)

Here, $\text{Br}_{\text{vis}}$ denotes the branching ratio of the FIP decays
into final states that are detectable; depending on the presence of a
calorimeter (EM and/or hadronic), $\text{Br}_{\text{vis}}$ may encompass only
those states featuring at least two charged particles, or it may also include
some neutral states such as photons and $K^{0}_{L}$.
$\epsilon_{\text{dec}}^{\text{(geom)}}$ denotes the fraction of decay products
that point to the end of the detector, and
$\epsilon_{\text{dec}}^{\text{(other cuts)}}$ is the fraction of these decay
products that additionally satisfy the remaining cuts (e.g., the energy cut,
etc.).

  * –

$\epsilon_{\text{rec}}$ is the reconstruction efficiency, i.e. among the FIP
decays that are within $\epsilon_{\text{az}}$ and $\epsilon_{\text{dec}}$, the
fraction of them that the detector can successfully reconstruct.

Figure 2: Examples of production processes for various FIPs: (a) proton
bremsstrahlung (for the dark photon $V$), (b) coherent scattering off nuclei
(for the ALP $a$ coupling to photons), (c) decays of $B$ mesons into a FIP and
another meson $h$ (for HNLs $N$).

Most quantities entering Eq. (1) can be accurately estimated analytically and
cross-checked separately, which makes the approach (1) very transparent.
Namely, the azimuthal acceptance is completely determined by the geometry of
the decay volume, which is typically very simple. Once $\epsilon_{\text{az}}$
is computed, a simple way to cross-check it is to verify that the integral

$\mathcal{V}=2\pi\int d\theta
drr^{2}(z,\theta)\sin(\theta)\epsilon_{\text{az}}=2\pi\int d\theta
dz\frac{z^{2}}{\cos^{3}(\theta)}\sin(\theta)\epsilon_{\text{az}}$ (5)

matches the total volume of the decay volume.

Depending on the production channel, evaluating the FIP distribution function
$f^{(i)}(\theta,E)$ may require some external input. This is the case, for
instance, for FIPs produced directly in inelastic proton collisions — where
one needs to simulate $f^{(i)}(\theta,E)$ using e.g. PYTHIA 8 to account for
showering and hadronization — or for FIPs produced in the interactions of
secondary particles, either in their decays or scattering with the material
(see Fig. 2 for examples) — where one needs to know the distribution of
secondaries $f_{\text{secondary}}(\theta,E)$. Nevertheless, in the latter
case, once $f_{\text{secondary}}(\theta,E)$ has been computed, the
distribution of FIPs can then be derived analytically without the need for
external tools.

$\epsilon_{\text{dec}}$ may in principle be estimated qualitatively by
comparing the opening angle $\Delta\theta_{\text{dec}}$ between the decay
products with the angle $\Delta\theta_{\text{det}}$ covered by the detector as
seen from the production point. In the simplest case of a 2-body decay into
two massless particles, the opening angle is $\Delta\theta_{\text{dec}}\simeq
2\arcsin(\gamma^{-1})$, where $\gamma$ is the boost factor of the FIP. If
$\Delta\theta_{\text{dec}}\gtrsim\Delta\theta_{\text{det}}$, then
$\epsilon_{\text{dec}}\approx 0$, otherwise $\epsilon_{\text{dec}}\approx 1$.
Because the detector angle is smallest at the beginning of the decay volume
while the opening angle decreases as $E_{\text{FIP}}^{-1}$,
$\epsilon_{\text{dec}}$ effectively imposes a cut from below on the FIP energy
and the displacement of its decay position from the beginning of the decay
volume. If the decay volume is simultaneously the detector (as in the case of,
e.g., neutrino detectors), $\epsilon_{\text{dec}}^{\text{(geom)}}\equiv 1$.

In order to more accurately estimate $\epsilon_{\text{dec}}$ — by accounting
for the experiment geometry, the presence of a dipole magnet, different FIP
decay topologies (such as multi-body decays or decays into unstable
particles), and various other selections imposed on the decay products — a
_separate_ simulation can be performed (see details about the simulation in
Sec. 3).

Finally, the computation of $\epsilon_{\text{rec}}$ would require running the
full simulation, including the detector response. As such, it goes beyond the
scope of the present semi-analytic approach. However, we believe that it is
possible to perform an adequate pre-selection with the help of
$\epsilon_{\text{dec}}^{\text{(other cuts)}}$ — for instance, by requiring a
minimum energy or $p_{\mathrm{T}}$ above which the particles are detected with
high efficiency (see, e.g., Aberle:2839677 ) — such that, conditioned on this
pre-selection, $\epsilon_{\text{rec}}\sim\mathcal{O}(1)$.

Last but not least, this semi-analytic method allows for a simple analysis of
the number of events in the limit of long lifetimes
$c\tau\langle\sqrt{\gamma^{2}-1}\rangle\gg l_{\text{experiment}}$, where
$l_{\text{experiment}}$ is the length scale of the experiment. The number of
events then reduces to a simple expression:

$N_{\text{ev}}\approx\sum_{i}N_{\text{prod}}^{(i)}\cdot\frac{(z_{\text{max}}-z_{\text{min}})}{c\tau}\times\epsilon,$
(6)

where $\epsilon$ is the total acceptance:

$\epsilon=\frac{1}{z_{\text{max}}-z_{\text{min}}}\int d\theta dEdz\
\epsilon_{\text{az}}\cdot\frac{\epsilon_{\text{dec}}}{\cos(\theta)\sqrt{\gamma^{2}-1}}$
(7)

It may be decomposed as

$\epsilon=\langle\epsilon_{\text{FIP}}\rangle\times\langle\epsilon_{\text{decay}}\rangle\times\langle(\gamma^{2}-1)^{-1/2}\rangle,$
(8)

where $\langle\epsilon_{\text{FIP}}\rangle$ is the mean probability for the
FIP to intersect the decay volume, $\langle\epsilon_{\text{decay}}\rangle$ is
the mean probability for the decay products to meet the decay products
acceptance criteria, and $\langle(\gamma^{2}-1)^{-1/2}\rangle$ is the mean
inverse $p/m$ among the FIPs meeting the azimuthal and decay acceptance
criteria. This representation is particularly useful when discussing the
impact of the geometry on the event rate and when comparing the potential of
various experimental setups Bondarenko:2023fex .

The semi-analytic approach presented here is also well suited for estimating
the sensitivity to FIP scatterings, which is the main signature in models of
light dark matter. In this case, the differential decay probability should be
replaced with the scattering probability

$\frac{dP_{\text{scatt}}}{d\theta
dEdz}=n_{\text{detector}}\frac{d^{2}\sigma_{\text{scatt}}}{d\theta dE},$ (9)

where $n_{\text{detector}}$ is the number density of target particles inside
the detector, and $d^{2}\sigma_{\text{scatt}}/d\theta dE$ is the differential
cross-section for the scattering of FIPs off the target particles.

###  2.2 Validation and limitations

The semi-analytic approach presented above has been used to estimate the
sensitivities of various experiments at the SPS Boiarska:2021yho , LHC
Boyarsky:2022epg ; Ovchynnikov:2022its , and FCC-hh Boyarsky:2022epg . The
experimental setups considered cover various options: on-axis and off-axis
placements, different decay volume shapes, and different detector orientations
relative to the beamline. These estimates, carried out using our semi-
analytical method, have been found to be in good agreement with the estimates
available in the literature, including simulations-based ones. In particular,
Fig. 3 shows the comparison of the sensitivity of the SHiP experiment to heavy
neutral leptons (HNLs) and dark photons obtained using Eq. (1) with the
sensitivity obtained by the SHiP collaboration using the FairShip simulation.
In the case of dark photons, the slight differences in the sensitivity can be
explained by the different elastic proton form factors used to describe the
production probability. In the case of HNLs, the discrepancy at the upper
bound follows from the monochromatic approximation of the HNL energy spectrum
used when computing the sensitivity shown in the SHiP paper SHiP:2018xqw .

Figure 3: Comparison of the 90% CL sensitivity of the SHiP experiment to heavy neutral leptons (left panel) and dark photons (right panel), obtained using Eq. (1) within the framework of SensCalc and derived using the FairShip simulations SHiP:2020vbd ; SHiP:2018xqw . The old ECN4 configuration of SHiP has been considered here. Experiment | SHADOWS | MATHUSLA@CMS
---|---|---
$(x,y,z)_{\text{min}}$, m | (-1,0,14) | (0,60,68)
Fid. dim, $\text{m}^{3}$ | $2.5\times 2.5\times 20$ | $100\times 25\times 100$
Det. dim., $\text{m}^{3}$ | $2.5\times 2.5\times 12$ | $100\times 5\times 100$
Detector plane | $xy$ | $xz$
Requirement for decay products |  Point to the end of detector Oppositely charged, or neutral No other cuts |  Point to the end of detector Oppositely charged No other cuts
$B$ distribution | CERN-SHiP-NOTE-2015-009  | Kling:2021fwx
Scalar production | Exclusive production, Boiarska:2019jym
Scalar decays | Following Boiarska:2019jym
Table 1: Description of the experimental setups and of the scalar
phenomenology used to obtain the sensitivity shown in Fig. 4. The rows
indicate respectively: the closest distance from the collision point to the
decay volume (the $z$ axis being along the beamline), the decay volume
dimensions, the detector dimensions, the orientation of detector layers, the
decay products acceptance criteria, the distribution of $B$ mesons used to
calculate the flux of scalars, the scalar production branching ratios, and the
description of the scalar lifetime and decays. The description of the
experiments has been taken from Refs. Alviggi:2839484  (SHADOWS) and
MATHUSLA:2022sze  (MATHUSLA@CMS). For the description of the scalar
production, we followed the PBC recommendations Beacham:2019nyx .

If the assumptions are well-controlled, the semi-analytic approach can agree
very well with simulations. Fig. 4 compares the sensitivity of SHADOWS and
MATHUSLA to dark scalars as computed via Eq. (1) and calculated independently
by SensMC, a simple weight-based Monte-Carlo that we have implemented as
described in Appendix C (see Table 1 for the detailed description of the setup
and of the scalar phenomenology used to compute the sensitivity). In these
calculations, we did not impose any cuts on the decay products apart from the
geometric requirement $\epsilon_{\text{decay}}^{\text{geom}}$; therefore the
sensitivities shown are optimistic. The agreement between the two approaches
is very good for all masses.

Figure 4: Comparison of the predictions of SensCalc (solid blue) with the
SensMC Monte-Carlo code used for validation (dashed black; described in App.
C), for the sensitivity of the experiments located off-axis. SHADOWS (left,
with the setup described in Ref. Alviggi:2839484 ) and MATHUSLA
MATHUSLA:2022sze  (right) are considered for the comparison. The description
of the experiments has been taken from the collaboration papers. To simplify
the comparison and because SensMC cannot simulate it, the effect of the dipole
magnet has not been included. The numbers of events produced by the two
approaches (both with and without including $\epsilon_{\text{dec}}$) mostly
agree within 20%. Discrepancies are caused mostly by different treatments of
the scalar decay and small differences in the numbers describing the
production and decay phenomenology of the scalar. The solid red lines show the
sensitivities reported in the collaboration documents Alviggi:2839484 ;
MATHUSLA:2022sze . The discrepancy between these calculations and our
estimates is discussed in the main text.

As a further demonstration of the importance of having open-access sensitivity
calculations with clear and controllable assumptions and inputs, we have also
included in Fig. 4 the sensitivities reported in the respective collaboration
papers: the SHADOWS LoI Alviggi:2839484 , and the MATHUSLA EoI
MATHUSLA:2022sze . These sensitivities differ greatly from those we obtained
for two main reasons. First, both collaborations use a different description
of the scalar production, based on the inclusive estimate where the decay of a
$B$ meson into a scalar is described as the decay of its constituent $b$
quark. Second, the assumptions about the experimental setups that have been
used to compute the sensitivity differ from what is actually described in the
documents. In the case of SHADOWS, Ref. Alviggi:2839484  used not the setup
described within that same work (and summarized in Table 1), but a more
optimistic setup located closer to the target and the beamline.333From private
communications with the representatives of the SHADOWS and MATHUSLA
collaborations. In the case of MATHUSLA, the acceptance of the decay products
was assumed to be $1$ in MATHUSLA:2022sze . These differences can
significantly affect the reported sensitivity.

Finally, the predictions of our method agree with other publicly available
packages — FORESEE and ALPINIST, as will be discussed in more detail in Sec.
3.2.

The simplicity of our semi-analytic method incurs a number of limitations.
First, SensCalc cannot provide the full event record associated with each FIP
decay or interaction, i.e. the set of all initial, intermediate and final-
state particles, including their full kinematics. Instead, it averages over
all events that pass the selection. Therefore, it does not allow studying the
reconstruction of the FIP parameters, such as its mass, for which the detailed
event information is essential. Second, this approach assumes that the
surrounding infrastructure does not influence the production of the FIPs.
While this is often true in the case of FIPs produced at the collision point
or close to it, the situation is different for non-prompt production; examples
of which include FIPs originating from the decays of long-lived $K^{\pm/0}$
mesons or from neutrino up-scatterings (the neutrino dipole portal
Ovchynnikov:2022rqj ; Ballett:2019bgd ), as well as the conversion of photons
into axion-like particles (ALPs) in the magnetic field at the LHC
Kling:2022ehv .

##  3 SensCalc

###  3.1 Description

Facility | List of experiments
---|---
SPS | SHiP Aberle:2839677 , NA62${}_{\text{dump}}$ CortinaGil:2839661 , HIKE${}_{\text{dump}}$ CortinaGil:2839661 , SHADOWS Alviggi:2839484
Fermilab (dump) | DUNE, DUNE-PRISM DUNE:2021tad , DarkQuest Batell:2020vqn
LHC |  FASER/FASER2/FASER$\nu$/FASER$\nu 2$ FASER:2019aik ; FASER:2020gpr ; Feng:2022inv  SND@LHC/advSND SHiP:2020sos ; Feng:2022inv  FACET Cerci:2021nlb , MATHUSLA MATHUSLA:2022sze , Codex-b Aielli:2019ivi , ANUBIS Bauer:2019vqk
FCC-hh | Analogs of the LHC-based experiments Boyarsky:2022epg
Table 2: List of the experiments whose geometry is currently implemented in
SensCalc, along with, for each experiment, a reference containing a
description of the setup used.

The code SensCalc consists of a few Mathematica notebooks that compute the
number of events for various FIPs (see Table 3 for the list of the currently
available models). Four notebooks have to be run sequentially: Acceptances.nb,
FIP distribution.nb, FIP sensitivity.nb, and Plots.nb, see Fig. 5.

Figure 5: Sketch of the modular structure of SensCalc. The notebook
Acceptances.nb produces the list of acceptances $\epsilon_{\text{az}}$ and
$\epsilon_{\text{dec}}$ entering Eq. (1) for the selected experiment. The
notebook FIP distribution.nb computes the distribution of FIPs $f(m,\theta,E)$
at the facility housing the experiment. The notebook FIP sensitivity.nb uses
as input the outputs of the two previous notebooks to calculate the tabulated
number of events, and then calculates the sensitivity in the mass-coupling
plane as a function of the remaining parameters such as the minimal number of
events and any additional model-specific parameters. Finally, Plots.nb
produces the sensitivity plots from the output of the previous notebook.

Figure 6: Visualizations of the geometries of the SHiP (left) and MATHUSLA
(right) experiments, as implemented in SensCalc (in the notebook
Acceptance.nb). The blue domain corresponds to the decay volume, while the red
domain shows the detector. The descriptions of the two geometries have been
taken from the SHiP LoI Aberle:2839677  and Ref. MATHUSLA:2022sze .

#### In the first notebook, Acceptances.nb,

the user specifies the experimental setup — the geometry and dimensions of the
decay volume and detector, as well as some details about the detector such as
the presence of an ECAL and dipole magnet, see Fig. 6. The list of the
experiments currently implemented in SensCalc is provided in Table 2. The user
can easily implement new experiments, or modify one of the already implemented
setups, which may be useful when optimizing an experiment. Some past
experiments are also included: CHARM CHARM:1983ayi  and BEBC BEBCWA66:1986err
at the SPS. In this notebook, the user must also provide all the relevant
quantities such as the number of protons on target (or the integrated
luminosity for LHC- and FCC-hh-based experiments), the target material, and
the production cross-sections for secondary particles (mesons and heavy
bosons). For the implemented experiments, these parameters are already listed
in the notebook.

Once the setup is fixed, the notebook evaluates the angular coverage of the
experiment and $\epsilon_{\text{dec}}$ for various FIPs. Concretely, it first
defines the grid of the FIP masses $m$, FIP energies $E$, and its decay
coordinates within the decay volume: the polar angle $\theta$ and the
longitudinal displacement from the target along the beam axis $z$. Using these
coordinates, the notebook then evaluates $\epsilon_{\text{az}}(\theta,z)$ and
the list of azimuthal angles $\phi$ for which the FIP is inside the decay
volume.

Having the grid $(E,\theta,z,\phi)$, the notebook then simulates the FIP
decays using its dominant decay channels and calculates the decay acceptance
$\epsilon_{\text{dec}}(m,\theta,E,z)$ by averaging over these decays and
$\phi$. The averaging over $\phi$ is already possible at this stage since the
other quantities that determine the number of events (1) do not depend on the
azimuthal angle. Namely, the differential decay probability
$dP_{\text{dec}}/dz$ only depends on $z$ and $\theta$, while the FIP
distribution function is typically isotropic in $\phi$.

The decay channels implemented for each FIP are listed in Table 4. For 3-body
decays, the distribution of the decay products is generated taking into
account both the phase space and the matrix element of the process. If the FIP
decay products are short-lived, the routine decays them until only metastable
particles are left. By default, those are
$\gamma,e,\mu,K^{0}_{L},\pi^{\pm},K^{\pm}$. The decays of particles with many
modes (such as $D,B,\tau$) are approximated by some representative decays; for
example, for $\tau$ this is a 3-body decay into one charged particle and two
neutrinos, while for $D,B$ these are multi-hadronic decays.

Let us now discuss the computation of $\epsilon_{\text{dec}}$ in more details.
The main acceptance criterion is the requirement that the trajectories of at
least two decay products with zero total electric charge are within the
acceptance of the detector until its final plane. Decays into pure neutral
final states (i.e., photons or $K^{0}_{L}$) are also included if a calorimeter
is present. If the detector includes a magnetic spectrometer, the components
of the charged particles’ coordinates and momenta are shifted by a kick right
after the magnet in order to approximate the effect of the magnetic field. In
addition to this geometric requirement, $\epsilon_{\text{dec}}$ may also
include various kinematic cuts. The implemented cuts include the energy cut,
the transverse momentum cut, the transverse impact parameter cut, and the
spatial separation cut for neutral particles in the calorimeter. By complete
analogy, the user may impose further kinematic cuts. Although the cuts are
applied at the Monte-Carlo truth level, i.e. they are implemented without
considering reconstruction effects such as the finite resolution of 4-momenta
measurements, they can already give us some understanding of the effects that
a realistic event reconstruction would have on the signal yield. Such
reconstruction effects could in principle be approximated by e.g. applying
some smearing to the kinematics variables of the decay products, according to
the detector resolution. Note that the acceptance criterion includes partially
reconstructible states, i.e. final states for which the FIP invariant mass
cannot be reconstructed from the detected decay products.

The output of the first notebook is a table with the following columns:

$\\{m,\theta,E,z,\epsilon_{\text{az}},\epsilon_{\text{dec}}\\}$ (10)

Model | Ref. | Production channels
---|---|---
BC1 | Ilten:2018crw ; SHiP:2020vbd  |  Decays of $\pi,\eta,\eta^{\prime}$, mixing with $\rho^{0}$ Proton bremsstrahlung, Drell-Yan process
BC4, BC5 | Boiarska:2019jym ; Winkler:2018qyg  |  2-/3-body decays of $B$, decay $h\to SS$
BC6-8 |  Bondarenko:2018ptm  |  2-/3-body decays of $B,D,W$
BC9 | Jerhot:2022chi ; Dobrich:2019dxc  |  Coherent production: Primakov process, $pZ$ scattering Decays of $\pi^{0},\eta$
BC10 | Beacham:2019nyx  |  $B$ decay
BC11 | Aloni:2018vki ; Chakraborty:2021wda ; Jerhot:2022chi  |  Decays of $B$, mixing with $\pi^{0}/\eta/\eta^{\prime}$ Deep-inelastic production
Table 3: FIP production channels in the various models implemented in SensCalc
with: the benchmark names according to PBC Beacham:2019nyx , the reference
used to describe the production channels, and the list of the production
channels implemented in SensCalc. The models are dark photons (BC1), dark
scalars with Higgs mixing (BC4) and also with the quartic coupling (BC5),
Heavy neutral leptons with arbitrary mixing patterns (including the limiting
cases of the pure mixing with $\nu_{e}$, $\nu_{\mu}$, or $\nu_{\tau}$
(BC6–BC8)), and ALPs coupling to photons (BC9), fermions (BC10) and gluons
(BC11).

#### The second notebook, FIP distribution.nb,

computes the angle-energy distribution of the FIPs produced by various
facilities and mechanisms. The list of implemented production channels, along
with the relevant references used to describe the production, can be found in
Table 3. Many production mechanisms require knowing the distributions of the
parent particles at the given facility, such as mesons, heavy SM bosons, and
photons — including those produced in secondary interactions. We provide them
as tabulated distributions in polar angle and energy, which we generate
following the literature or just using available distributions from existing
studies (see also Appendix B for a description of how we have generated the
distributions of parent particles). Users may easily replace the included
distributions with their own differential flux. With the distribution of
parent particles at hand, we then derive the distribution of FIPs. If the FIPs
are produced in decays, we compute their phase space in the rest frame of the
parent particle and then boost it to the lab frame. In the case of 3-body
decays, the phase space takes into account the matrix element of the process.
For FIPs produced via elastic scattering, we adopt the differential cross-
section of the process from existing studies, and then convolve it with the
distribution of the parent particles. Should the need arise, new production
channels may be added by the user, following the above examples.

Such a derivation of the FIP distribution is not possible, however, in the
case of FIPs that are produced inelastically in proton-proton collisions (such
as via the Drell-Yan process for dark photons or deep-inelastic production of
ALPs through the gluon coupling), which require an external simulation. In
this case, we use MadGraph5_aMC@NLO (v3.4.2) Alwall:2014hca  with a model
implemented in FeynRules Alloul:2013bka  and exported to the UFO format
Degrande:2011ua . To account for showering and hadronization, the events
simulated in MadGraph are further processed by PYTHIA 8 Sjostrand:2014zea ;
see also Appendix B.1 for details. The UFO files and the tabulated FIP
distributions are provided alongside SensCalc.

The output of the second notebook is a tabulated distribution of the form

$\\{m,\theta,E,f^{(i)}\\},$ (11)

where the last column is the value of the FIP distribution function for the
given $(m,\theta,E)$ and the production mechanism $i$. Some examples of
computed distribution functions are shown in Fig. 7.

Model | Ref. | Decay channels | Decay channels ($\epsilon_{\text{dec}}$)
---|---|---|---
BC1 | Ilten:2018crw ; SHiP:2020vbd  |  $ee,\mu\mu,\tau\tau$ $\pi\pi,3\pi,4\pi,KK,m\lesssim 2\text{ GeV}$ $q\bar{q},m\gtrsim 2\text{ GeV}$ |  $ee,\mu\mu$ $\pi\pi,4\pi,m_{V}<2\text{ GeV}$ $q\bar{q},m_{V}>2\text{ GeV}$
BC4, BC5 | Boiarska:2019jym ; Winkler:2018qyg  |  $ee,\mu\mu,\tau\tau$, $\pi\pi,KK,4\pi,DD,BB$ |  $ee,\mu\mu,\tau\tau$ $\pi\pi,KK,DD,\tau\tau,BB$
BC6-8 | Bondarenko:2018ptm  |  $3\nu,ll\nu$ $\text{meson}+l/\nu,m\lesssim 1\text{ GeV}$ $\nu q\bar{q},lq\bar{q}^{\prime},m\gtrsim 1\text{ GeV}$ |  $ll\nu$ $lq\bar{q}^{\prime}/\nu q\bar{q},m\geq 1.5\text{ GeV}$ $\text{meson}+l/\nu,m<1.5\text{ GeV}$
BC9 | Jerhot:2022chi  | $\gamma\gamma$ | $\gamma\gamma$
BC10 | Beacham:2019nyx  | $ee,\mu\mu,\tau\tau$ | $ee,\mu\mu,\tau\tau$
BC11 | Jerhot:2022chi ; Aloni:2018vki  |  $\gamma\gamma$ $\gamma\pi\pi,\eta\pi\pi,3\pi,m<1.5\text{ GeV}$ $GG,m>1.5\text{ GeV}$ |  $\gamma\gamma,GG$ $\gamma\pi\pi,\eta\pi\pi$
Table 4: Decay channels of the FIPs implemented in SensCalc. From left to
right: the model name according to PBC Beacham:2019nyx  (see also the caption
of Table 3), the reference used to describe the decays, the decay channels
used to calculate the lifetime $\tau_{\text{FIP}}$ and the branching ratio of
visible decays, and the decay channels used to calculate the decay acceptance
$\epsilon_{\text{dec}}$.

Figure 7: Examples of angle-energy distributions $f^{(i)}(\theta,E)$ for ALPs
coupled to photons (left) and dark scalars with a non-zero quartic coupling
(right), produced by the notebook FIP distribution.nb.

Let us highlight an important point. Since the FIP distributions are
determined mainly by the kinematics of the collisions, they can be considered
identical for the different experiments housed at a same facility, assuming
that the colliding particles are the same.444This is not the case for non-
prompt production of FIPs, which goes beyond the scope of the present
discussion. For collider experiments, we typically deal with proton-proton
collisions, and this notebook only needs to be run once to obtain the
distributions. In the case of beam dump experiments, some differences may
arise as a result of different target/beam dump compositions. When the FIP is
produced via the decays of secondaries, this only affects the overall scaling
of the secondaries production cross section, which depends on the atomic
number $A$: $\sigma_{\text{prod,second}}\propto A^{0.29}$ Carvalho:2003pza .
Therefore, as in the collider case, the notebook only needs to be run once.
If, however, the FIP is produced in scattering processes, then different
targets may affect not only the normalization but also the shape of the
distribution. To take this into account, we generate the fluxes for a few
common types of targets.

#### The notebooks <FIP> sensitivity.nb

(with <FIP> replaced by the actual FIP) evaluate the sensitivity of the chosen
experiment to the corresponding FIP. This is done via computing a tabulated
number of events. There is a dedicated notebook for each FIP. First, the
notebook imports the acceptance data computed by Acceptances.nb, the
distributions produced by FIP distribution.nb, as well as the relevant
quantities defining the FIP phenomenology, such as the production branching
ratios, lifetimes, and branching ratios of the decays into visible states at
the given experiment. It then maps them to logarithmic scale and interpolates
them to obtain the functions entering Eq. (1).

Depending on the FIP, there may exist uncertainties in the description of its
production and decay, which may significantly affect the event rate. This is
the case, e.g., for dark scalars, where one may describe their production
inclusively or exclusively; and for dark photons, for which the description of
the proton bremsstrahlung channel depends on the maximal allowed $p_{T}$ and
on the minimal energy allowed to be transferred to the dark photon. The user
has the freedom to tune these parameters.

In addition, there may exist model-specific parameters that must be selected
before performing the computation. For instance, in the case of HNLs, this is
their nature (Dirac or Majorana) and their mixing pattern
$U_{e}^{2}:U_{\mu}^{2}:U_{\tau}^{2}$.

Figure 8: Examples of the output produced by the notebook FIP sensitivity.nb.
Left panel: differential number of events with respect to the FIP’s energy for
various production channels. Right panel: density plot of the total number of
events as a function of the FIP mass and coupling. As an example, dark photons
at FACET are considered. No cuts on the decay products other than the
geometric acceptance have been applied.

During the computation, this notebook produces intermediate results that may
be useful for the sensitivity analysis. This includes the differential number
of events with respect to $\theta,E$, or $z$, as well as the number of events
as a function of the mass and coupling, see Fig. 8. Last but not least, the
notebook also shows the behavior of the overall acceptances $\epsilon$, cf.
Eq. (7).

Once the tabulated number of events has been produced, the notebook computes
the sensitivities. To this end, the user needs to select the critical number
of events as well as some model-specific parameters. For example, for dark
scalars, one needs to specify the value of the branching ratio $\text{Br}(h\to
SS)$, which is non-zero in the presence of the quartic coupling
$\mathcal{L}\propto hSS$ (see Appendix A for detail). Because the critical
number of events can be freely specified, the user can compute both
“exclusion” sensitivity limits — corresponding to $2.3$ expected events at
$90\%$ CL — or “discovery” sensitivity limits by (externally) providing the
critical $N_{\text{ev}}$ corresponding to the desired significance level and
background expectation.

#### Finally, the notebook Plots.nb

plots the sensitivities obtained in the previous notebook. It scans over
available sensitivity files, imports those needed by the user, and finally
produces the figures (see e.g. Fig. 9).

Figure 9: Example of a sensitivity plot produced by the notebook Plots.nb, for
the model of dark scalars.

The user interaction with the various notebooks, such as choosing the
experiment, selecting the cuts and the particular FIP model, is organized via
dialog windows. This makes running the notebooks straightforward for FIPs and
experiments that are already implemented.

To successfully run the notebooks, the user needs to install two dependencies:
FeynCalc Shtabovenko:2020gxv , which is a Mathematica package for the symbolic
evaluation of Feynman diagrams, and a C compiler that is recognized by
Mathematica.

The performance of the code has been tested on various machines and operating
systems. For instance, on a Windows laptop with 16 GB of RAM, 8 CPU cores, and
Mathematica 12.1, the typical time required to compute the sensitivity from
scratch is $\mathcal{O}(1\text{ hour})$ — depending on the FIP type and on the
mass-coupling grid density. This time is reduced if the FIP distribution has
already been pre-generated.

SensCalc still offers significant potential for further improvement. Of
particular interest would be the possibility to compute the sensitivity to
additional FIP models, including those for which the main signature is
scatterings with the detector material. Another well-motivated extension would
be to support ALPs with an arbitrary coupling pattern.

A particularly interesting improvement would be to implement an approximate
description of the hadronic decays of heavy FIPs ($m_{\text{FIP}}\gg 1\text{
GeV}$). Currently, SensCalc describes them perturbatively – via decays into
quarks and gluons. Their subsequent showering and hadronization in-flight
would require external tools such as PYTHIA 8, and hence goes beyond the scope
of the present semi-analytic approach.555Nevertheless, the accuracy of the
current method remains surprisingly good even in this regime. This can be
understood as follows: because the momentum flow of the hadronized decay
products is determined by the kinematics of the incoming jets, the geometric
part of $\epsilon_{\text{dec}}$ (Eq. (4)) still describes the decay
adequately. This qualitative argument is in agreement with Fig. 3: the
dominant decays of the HNLs and dark photons with mass $m_{\text{FIP}}\gtrsim
1\text{ GeV}$ include a quark-antiquark pair, and nevertheless, the semi-
analytic approach agrees well with simulations. However, it may be possible to
implement this feature approximately, for instance by simulating the phase
space of FIP decays for several masses in PYTHIA 8, then selecting typical
sets of decay products, and finally using their pre-computed phase space when
evaluating the decay products acceptance.

Finally, the implementations of the various experiments should be updated
according to their latest specifications, which may differ from those listed
in currently available documents. This may be done by contacting the
representatives of the collaborations.

We are planning to add the above features in a future code update.

###  3.2 Comparison with similar software packages

At the moment of releasing SensCalc, there are two publicly available codes
for computing the sensitivity of lifetime-frontier experiments to decaying
FIPs: FORESEE Kling:2021fwx  and ALPINIST Jerhot:2022chi .

FORESEE is a Python-based code developed to evaluate the sensitivities of the
far-forward experiments at the LHC and FCC-hh. The currently implemented
models of FIPs include dark scalars, dark photons, ALPs coupling to $W$
bosons, millicharged particles, and up-philic scalars. The package includes
the tabulated distributions of various SM particles, including photons,
mesons, and heavy bosons. Apart from the tabulated number of events as a
function of the FIP mass and coupling, it can additionally produce detailed
event records in the HepMC format, which may then be passed to e.g. a detector
simulation software. By default, FORESEE does not calculate the acceptance of
the decay products; instead it only requires the FIP to decay inside the decay
volume, although the user may impose various cuts.

ALPINIST computes the sensitivity of extracted-beam experiments — including
those at the SPS, Fermilab, and some past experiments — to ALPs couplings to
various SM particles. Its modules use Mathematica, ROOT, and Python. The
prominent feature of the code is that it can handle generic ALPs with
simultaneous couplings to $W$ bosons, gluons, and the $U_{Y}(1)$ field. Unlike
FORESEE, the computation also incorporates the reconstruction of the decay
products inside a detector, at the price of a longer computation time to
obtain the tabulated number of events. Only fully reconstructible states are
considered. The output of ALPINIST consists of data files with the mass-
coupling dependence of the number of events for various production and decay
modes.

Figure 10: Comparison of the sensitivity of SHiP to ALPs coupling to photons
as computed by SencCalc (blue line) and ALPINIST Jerhot:2022chi  (red line).
The definition of the coupling, as well as the fractions of mesons leading to
the production of ALPs via photo-conversion, are taken from Ref.
Jerhot:2022chi . The use of different setups for SHiP in the sensitivity
calculations may explain the discrepancy at the upper bound: we have used the
currently considered ECN3 configuration from Ref. Aberle:2839677 , for which
the decay volume is located 12 m closer to the target, while Ref.
Jerhot:2022chi  considers the old configuration from Ref. Ahdida:2654870 .

The predictions of SensCalc agree well with the results of ALPINIST (see Fig.
10) and FORESEE (the comparison between the semi-analytic approach and FORESEE
is discussed in Ref. Ovchynnikov:2022its ).

Unlike these two software packages, SensCalc is not restricted to a particular
facility. In addition, among the implemented FIP models, it considers for the
first time HNLs with arbitrary mixing patterns. The main limitation of
SensCalc compared to FORESEE is that it cannot generate detailed event
records, while compared to ALPINIST, it is that it does not (currently)
consider generic ALPs and does not perform a detailed event reconstruction.

##  4 Conclusion

Feebly interacting particles (FIPs) are present in a broad class of new-
physics scenarios that attempt to resolve the known problems of the Standard
Model. Their search at various facilities and experiments collectively forms
the lifetime frontier of particle physics. During the last decade, many
lifetime-frontier experiments have been proposed, that differ in the housing
facility, geometric location, and detector technology. With a few exceptions,
most of these experiments are not approved yet, and their design is not
finalized. Their sensitivities to FIPs are computed by the collaborations
themselves, using internal tools which are not publicly accessible. This makes
it difficult to control the inputs to the computations, such as the model of
the production and decay. It is therefore crucial to have a publicly available
tool for computing the sensitivity of those experiments to various FIPs in a
uniform, fast and well-controlled way.

The present paper addresses this issue by presenting SensCalc — a Mathematica-
based code for evaluating the sensitivity of various experiments to decaying,
long-lived FIPs, based on a semi-analytic approach developed in a number of
previous studies (see Sec. 2.1) and cross-checked against various state-of-
the-art packages (see Sec. 2.2).

SensCalc already supports a broad range of models and experiments (see Sec.
3.1), with more to be added in future versions. Models currently implemented
include dark photons, dark scalars, heavy neutral leptons with various mixing
patterns, and axion-like particles coupled to different SM particles. Numerous
experiments have been implemented, located at any of the following facilities:
the SPS, LBNE, LHC, and FCC-hh. The user retains full control over every
aspect of the sensitivity calculation: from the geometry of the experiment and
the distribution of the FIP’s parent particles to the branching ratios of the
FIP production/decay modes and the requirements on the decay products. Besides
contributing to the transparency and trustworthiness of the results, this also
allows users to easily modify the underlying assumptions as needed, or to add
their own models and experiments to SensCalc.

By publicly providing a transparent, semi-analytic method to consistently
compute the expected signal at various lifetime-frontier experiments, SensCalc
can help address the discrepancies that currently exist in the literature
between the descriptions of FIPs and acceptances employed by different
collaborations. This is a timely and necessary contribution to the field of
FIP searches, as many experiments are currently undergoing active development
and optimization, while funding bodies and hosting facilities must decide
which projects to prioritize. SensCalc can help with the former by providing
fast (re-)calculation of the expected signal as the experiment’s design
evolves, and with the latter by ensuring a fair and consistent comparison of
the expected signals between the proposed experiments, with well-controlled
assumptions thanks to a uniform and well-validated implementation of the
official PBC benchmarks. This could be particularly relevant in the context of
the ECN3 hall upgrade at the CERN SPS, in which a number of experiments are
currently being considered for inclusion, namely HIKE, SHiP, and SHADOWS.

### Acknowledgments

We thank Alexey Boyarsky and Oleg Ruchayskiy, who supervised the authors on
the present topic in the past, and helped develop the foundations of the
approach described in this paper. We thank Felix Kahlhoefer and Jan Jerhot for
helpful discussions on the phenomenology of ALPs and the ALPINIST code, and
Felix Kling for discussions on FORESEE. We also thank Thomas Schwetz, Nashwan
Sabti, Vsevolod Syvolap, Felix Kahlhoefer and Inar Timiryasov for reading the
manuscript at different stages of its writing. Finally, we thank the users of
the mathematica.stackexchange.com website, who greatly helped us optimize some
elements of the code. MO received support from the European Union’s Horizon
2020 research and innovation program under the Marie Sklodowska-Curie grant
agreement No. 860881-HIDDeN. OM is supported by the NWO Physics Vrij Programme
“The Hidden Universe of Weakly Interacting Particles” with project number
680.92.18.03 (NWO Vrije Programma), which is (partly) financed by the Dutch
Research Council (NWO). KB is partly funded by the INFN PD51 INDARK grant. JLT
acknowledges partial financial support by the Spanish Research Agency (Agencia
Estatal de Investigación) through the grant IFT Centro de Excelencia Severo
Ochoa No CEX2020-001007-S, by the grant PID2019-108892RB-I00 funded by
MCIN/AEI/ 10.13039/501100011033, by the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No
860881-HIDDeN, and by the grant Juan de la Cierva FJC2021-047666-I funded by
MCIN/AEI/10.13039/501100011033 and by the European Union
“NextGenerationEU”/PRTR.

### Conflict of Interest Statement

The authors of the present manuscript are also members of the SHiP
collaboration, which represents one of the experimental proposals currently
competing for funding and access to facilities, notably as part of the ongoing
Physics Beyond Colliders study and in the context of the upcoming upgrade of
the ECN3 hall at CERN. The present manuscript solely reflects the authors’
views, and not those of the SHiP collaboration.

##  Appendix A Uncertainties in the description of FIPs

###  A.1 Discrepancies in the literature

The description of the FIP production and decay, and sometimes even the
definition of the FIP couplings, may vary among the sensitivity estimates
performed by the different collaborations. One example is the dark scalar $S$.
Following the PBC report Beacham:2019nyx , the SHiP collaboration uses the
exclusive description of the production of $S$, while other collaborations
adopt instead the inclusive description (see the discussion in Ref.
Boiarska:2019jcw ). In the domain $m_{S}\gtrsim 2-3\text{ GeV}$, where the
inclusive approach breaks down, the difference in the number of produced
scalars between these two descriptions may be a factor of 20 or more. Another
problem arises from the theoretical uncertainty on the hadronic decay width,
which may be as large as a factor of 100 Winkler:2018qyg ; Boiarska:2019jcw
(see also a recent discussion in Ref. Gorbunov:2023lga ). While SHiP and
SHADOWS assume the decay width computed in Ref. Winkler:2018qyg , the FASER
collaboration Feng:2017vli  uses the decay width from Ref. Bezrukov:2009yw .
Depending on the calculation used, the sensitivity may therefore differ
significantly.

Another example is with ALPs $a$ coupling to gluons. The PBC report defines an
interaction of the form $\mathcal{L}\propto
ag_{a}G^{\mu\nu,a}\tilde{G}_{\mu\nu}^{a}$, where $G^{\mu\nu}$ is the gluon
field strength and $g_{a}$ is a fixed dimensionful coupling. Theoretical works
often Aloni:2018vki ; Chakraborty:2021wda  adopt a different definition,
$\mathcal{L}\propto ag_{s}^{2}g_{a}G^{\mu\nu,a}\tilde{G}_{\mu\nu}^{a}$, where
$g_{s}=g_{s}(m_{a})$ is the QCD coupling. The latter definition is used by
ALPINIST Jerhot:2022chi  for computing the sensitivity of beam dump
experiments to ALPs (and their results are used by the SHiP, HIKE, and SHADOWS
collaborations in Ref. Antel:2023hkf ). Furthermore, while some collaborations
Aielli:2019ivi  include the production of ALPs through gluon fusion, others do
not (this is the case in particular of ALPINIST Jerhot:2022chi ).

Another problem arises with ALPs that couple to fermions. The PBC
Beacham:2019nyx  recommends including only the decays into leptons in the
total width — even though it may be dominated by hadronic decays in the mass
range $m_{a}\gtrsim 2m_{\pi}$ — while some collaborations also include
hadronic channels Aielli:2019ivi .

Such mismatches between the assumptions used to compute different
sensitivities are particularly problematic when said sensitivities are shown
in the same plot — such as e.g. in the FIPs 2022 proceedings Antel:2023hkf  —
without emphasizing that the underlying assumptions differ.

###  A.2 Definition of the FIP couplings used in SensCalc

The effective Lagrangians of the models implemented in SensCalc are:

  * –

BC1 (dark photons):

$\mathcal{L}_{\text{int}}=-\epsilon eV_{\mu}J^{\mu}_{\text{EM}}$ (12)

where $V_{\mu}$ is the dark photon field, $J_{\text{EM}}^{\mu}$ is the EM
current, and $e=\sqrt{4\pi\alpha_{\text{EM}}}$ is the EM coupling.

  * –

BC4 and BC5 (dark scalars):

$\mathcal{L}_{\text{eff}}\supset m_{h}^{2}\theta hS+\frac{\alpha}{2}hS^{2},$
(13)

where $\theta$ is the mixing angle and $\alpha$ is the quartic coupling. By
default, the sensitivity is evaluated assuming a constant branching ratio
$\text{Br}(h\to SS)\propto\alpha^{2}$.

  * –

BC6, BC7, BC8 (HNLs):

$\mathcal{L}_{\text{int}}=\sum_{\alpha=e,\mu,\tau}U_{\alpha}\bar{N}\left(\frac{g}{\sqrt{2}}\gamma^{\mu}P_{L}l_{\alpha}W_{\mu}+\frac{g}{2\cos(\theta_{W})}\gamma^{\mu}P_{L}\nu_{\alpha}Z_{\mu}\right)+\text{h.c.},$
(14)

where $N$ is the HNL, $U_{\alpha}$ the mixing angle, $g$ the weak coupling,
and $l_{\alpha},\nu_{\alpha},W,Z$ the SM fields. The HNL may be either a Dirac
or a Majorana particle.

  * –

BC9 (ALPs coupling to photons):

$\mathcal{L}_{\text{int}}=\frac{g_{a}}{4}aF_{\mu\nu}\tilde{F}^{\mu\nu},$ (15)

where $a$ is the ALP field, $g_{a}$ is a dimensionful coupling, and
$F_{\mu\nu},\tilde{F}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}F^{\alpha\beta}$
are the EM field strength and its dual.

  * –

BC10 (ALPs coupling to fermions):

$\mathcal{L}_{\text{int}}=\frac{g_{Y}}{2v_{H}}(\partial_{\mu}a)\sum_{\alpha}\bar{f}\gamma^{\mu}\gamma_{5}f,$
(16)

where $g_{Y}$ is a dimensionless coupling, $v_{H}\approx 246\text{ GeV}$ is
the Higgs VEV, and $f$ are SM fermions.

  * –

BC11 (ALPs coupling to gluons):

$\mathcal{L}_{\text{int}}=g_{a}g_{s}^{2}aG_{\mu\nu}^{a}\tilde{G}^{\mu\nu,a},$
(17)

where $g_{s}$ is the strong coupling constant, $a$ is the ALP field, $g_{a}$
is a dimensionful constant, $G_{\mu\nu}^{a}$ is the gluon field strength, and
$\tilde{G}_{\mu\nu}^{a}=\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}G^{\alpha\beta,a}$
is its dual field strength. Everywhere except for the production of ALPs from
DIS, we follow the definition of $g_{s}$ from Ref. Jerhot:2022chi . In the DIS
case, we employ the running of $g_{s}$ associated with the default PDF set in
MadGraph.

##  Appendix B Inputs used for generating the FIP distributions

###  B.1 DIS processes

DIS production channels are relevant for several FIPs considered in SensCalc:
dark photons and ALPs with the gluon coupling. We simulate the production of
these particles in MadGraph5_aMC@NLO, interfaced with PYTHIA 8 for showering
and hadronization. The hard processes that we simulate are the lowest-order
and next-to-lowest-order processes for quark and gluon fusion:

$q+\bar{q}\to V,\quad q+\bar{q}\to V+j,\quad G+G\to a,\quad G+G\to a+j$ (18)

As for the scale of the process, we choose the invariant mass of the quark-
antiquark pair (dynamical_scale_choice = 4). Although SensCalc already
includes the tabulated angle-energy distributions of the FIPs produced by DIS,
it also includes the UFO files for the models of dark photons and ALPs,
allowing the user to re-generate the distributions under different assumptions
if needed.

Particle | LBNE | SPS | LHC | FCC-hh
---|---|---|---|---
$\pi^{0}/\eta/\eta^{\prime}/\gamma$ | Dobrich:2019dxc  | Dobrich:2019dxc  | Pierog:2013ria  | Pierog:2013ria
$B,D$ | Jerhot:2022chi  | CERN-SHiP-NOTE-2015-009  | Kling:2021fwx  | Kling:2021fwx
$W,h,Z$ | – | – | Kling:2021fwx  | Kling:2021fwx
Table 5: List of the references used to generate, or directly take, the
distributions of secondary particles that may produce FIPs.

The DIS production suffers from significant theoretical uncertainties. First,
the choice of scale becomes important for light FIPs with masses
$m_{\text{FIP}}\simeq 1-2\text{ GeV}$, where the uncertainties in the
production cross-section may become $\mathcal{O}(1)$. Second, the minimal
parton energy fraction required to produce a FIP is
$x_{\text{min}}=m_{\text{FIP}}^{2}/s_{\text{pp}}$. For experiments like the
LHC/FCC-hh and GeV-scale FIPs, $x_{\text{min}}$ can be as tiny as $10^{-8}$;
this domain is only explored experimentally and is therefore subject to
theoretical uncertainties (see Ref. Berlin:2018jbm ). This becomes especially
problematic in the case of the FCC-hh. Because of this, we do not consider the
DIS production channel for the FCC-hh-based experiments.

###  B.2 Production by secondary particles

Another important FIP production mechanism is through secondary particles —
either in their decays or scatterings. We handle this case by either
generating the distributions of secondary particles using approaches from the
literature, or directly using pre-calculated distributions. The list of
references is provided in Table 5.

##  Appendix C SensMC: a simplified Monte-Carlo used for validation

As an additional cross-check of SensCalc, we have implemented SensMC SensMC-
GitHub , a small, customizable weight-based Monte-Carlo simulation, as an
alternative way of numerically integrating Eq. (1) for FIPs produced in meson
decays. It makes extensive use of importance sampling in order to handle the
(typically tiny) branching ratios of mesons to FIPs and the (possibly very
displaced) decay vertex of the FIP. SensMC is written in the Julia programming
language julialang  in order to combine performance and readability, and it is
released alongside SensCalc in the same repository SensCalc-Zenodo , as well
as on GitHub.666The GitHub repository can be found at
https://github.com/JLTastet/SensMC .

SensMC numerically estimates Eq. (1) using Monte-Carlo integration with
importance sampling, by randomly generating a large number of weighted samples
whose expectation values are $N_{\mathrm{ev}}$, and finally averaging them.
The value of each random sample is computed as follows:

  1. 1.

A meson species is randomly sampled based on the proportion of produced mesons
of this species, with the event weight initially set to the total number of
mesons produced across all species. The meson momentum is then randomly
sampled from a precomputed spectrum (either a list for the spectrums from
FairShip SHiP:2018xqw  or a grid for those from FORESEE Kling:2021fwx ). To
account for potential variations in the atomic weight of the target, that
would affect the overall normalization of the spectrums, the event is
optionally reweighted using the formula $w_{A}=w_{\mathrm{Mo}}(A/96)^{0.29}$
Carvalho:2003pza , with $A$ denoting the atomic weight of the target and
assuming that the spectrums were initially computed for a molybdenum target
(as is the case for the FairShip spectrums).

  2. 2.

The FIP production channel is randomly selected with a probability
proportional to its branching ratio, and the event is reweighted by the total
branching ratio to FIPs of the particular meson species. Upon the decay of the
parent meson, the momenta of its decay products, including the FIP, are
uniformly sampled in phase space. The present simulation currently does not
take into account the matrix elements because it cannot compute them all,
however the logic needed to handle them is already present, allowing the user
to implement their own matrix elements if needed.

  3. 3.

The FIP’s decay vertex is then selected randomly along its trajectory by
either a) sampling the proper lifetime from an exponential distribution and
calculating the corresponding distance in the lab frame or b) employing
importance sampling, which restricts the position of the decay vertex to a
shell covering the full decay volume, and then reweights the event by the
ratio of the true decay distribution to the importance distribution. The FIP
decay mode is selected similarly to its production mode, with a sampling
probability proportional (and in most cases equal) to its branching ratio; and
the event is reweighted by the total branching ratio of the implemented
channels. The momenta of the FIP decay products are uniformly sampled in phase
space in the current version (but matrix elements could in principle be taken
into account, just like for the FIP production).

  4. 4.

Following a similar procedure, any unstable Standard Model particles are
recursively decayed until only metastable particles (that live long enough to
be detected) remain, assuming the branching ratios listed in the particletools
Python package. The acceptance condition is then evaluated on the set of final
metastable particles produced in the FIP decay. The event weight is recorded,
along with whether the event is accepted or not.

Because each event is initially weighted by the total number of mesons, all
event weights must finally be divided by the number of generated events. The
sum of weights then provides a numerical estimate of the total number of
physical events (with the FIP decay within the ”shell” in case importance
sampling is used), while the sum of event weights multiplied by their
corresponding (binary) acceptances gives the total number of accepted events,
and is independent of the specific importance distribution as long as it fully
covers the decay volume.

The sensitivity curve is computed iteratively, starting from a coarse grid in
$(\log(m),\allowbreak\log(\theta))$ that covers the region where the
experiment is susceptible to be sensitive. The expected number of accepted
events is computed at each grid point. The multi-dimensional bisection method
(MDBM) MDBM  is then used to iteratively refine the grid in the vicinity of
the iso-contour corresponding to (for example) 2.3 accepted events (for an
exclusion sensitivity at the 90% confidence level), effectively bisecting it
without the need to evaluate a dense grid, which would be computationally
costly. The final curve is then obtained from bilinear interpolation of the
sparse grid values.

## References

  * (1) J. Beacham et al., Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report, J. Phys. G 47 (2020), no. 1 010501, [arXiv:1901.09966].
  * (2) S. Alekhin et al., A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case, Rept. Prog. Phys. 79 (2016), no. 12 124201, [arXiv:1504.04855].
  * (3) SHiP Collaboration, O. Aberle et al., BDF/SHiP at the ECN3 high-intensity beam facility, tech. rep., CERN, Geneva, 2022.
  * (4) SHADOWS Collaboration, M. Alviggi et al., SHADOWS Letter of Intent, tech. rep., CERN, Geneva, 2022.
  * (5) HIKE Collaboration, E. Cortina Gil et al., HIKE, High Intensity Kaon Experiments at the CERN SPS: Letter of Intent, tech. rep., CERN, Geneva, 2022.  Letter of Intent submitted to CERN SPSC. Address all correspondence to<EMAIL_ADDRESS>
  * (6) DUNE Collaboration, B. Abi et al., Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume II: DUNE Physics, arXiv:2002.03005.
  * (7) DUNE Collaboration, B. Abi et al., Prospects for beyond the Standard Model physics searches at the Deep Underground Neutrino Experiment, Eur. Phys. J. C 81 (2021), no. 4 322, [arXiv:2008.12769].
  * (8) B. Batell, J. A. Evans, S. Gori, and M. Rai, Dark Scalars and Heavy Neutral Leptons at DarkQuest, JHEP 05 (2021) 049, [arXiv:2008.08108].
  * (9) MATHUSLA Collaboration, H. Lubatti et al., Explore the lifetime frontier with MATHUSLA, JINST 15 (2020), no. 06 C06026, [arXiv:1901.04040].
  * (10) S. Cerci et al., FACET: A new long-lived particle detector in the very forward region of the CMS experiment, JHEP 2022 (2022), no. 06 110, [arXiv:2201.00019].
  * (11) FASER Collaboration, A. Ariga et al., Technical Proposal for FASER: ForwArd Search ExpeRiment at the LHC, arXiv:1812.09139.
  * (12) SHiP Collaboration, C. Ahdida et al., SND@LHC, arXiv:2002.08722.
  * (13) M. Bauer, O. Brandt, L. Lee, and C. Ohm, ANUBIS: Proposal to search for long-lived neutral particles in CERN service shafts, arXiv:1909.13022.
  * (14) G. Aielli et al., Expression of interest for the CODEX-b detector, Eur. Phys. J. C 80 (2020), no. 12 1177, [arXiv:1911.00481].
  * (15) D. Dercks, H. K. Dreiner, M. Hirsch, and Z. S. Wang, Long-Lived Fermions at AL3X, Phys. Rev. D 99 (2019), no. 5 055020, [arXiv:1811.01995].
  * (16) A. Boyarsky, O. Mikulenko, M. Ovchynnikov, and L. Shchutska, Exploring the potential of FCC-hh to search for particles from B mesons, JHEP 01 (2023) 042, [arXiv:2204.01622].
  * (17) F. Kling and S. Trojanowski, Forward experiment sensitivity estimator for the LHC and future hadron colliders, Phys. Rev. D 104 (2021), no. 3 035012, [arXiv:2105.07077].
  * (18) J. Jerhot, B. Döbrich, F. Ertas, F. Kahlhoefer, and T. Spadaro, ALPINIST: Axion-Like Particles In Numerous Interactions Simulated and Tabulated, JHEP 07 (2022) 094, [arXiv:2201.05170].
  * (19) Wolfram Research, Inc., “Mathematica, Version 13.2.”  Champaign, IL, 2022.
  * (20) M. Ovchynnikov, SensCalc, https://doi.org/10.5281/zenodo.7957785, 2023\.
  * (21) K. Bondarenko, A. Boyarsky, M. Ovchynnikov, and O. Ruchayskiy, Sensitivity of the intensity frontier experiments for neutrino and scalar portals: analytic estimates, JHEP 08 (2019) 061, [arXiv:1902.06240].
  * (22) I. Boiarska, A. Boyarsky, O. Mikulenko, and M. Ovchynnikov, Constraints from the CHARM experiment on heavy neutral leptons with tau mixing, Phys. Rev. D 104 (2021), no. 9 095019, [arXiv:2107.14685].
  * (23) M. Ovchynnikov, V. Kryshtal, and K. Bondarenko, Sensitivity of the FACET experiment to Heavy Neutral Leptons and Dark Scalars, JHEP 02 (2023) 056, [arXiv:2209.14870].
  * (24) P. Coloma, J. López-Pavón, L. Molina-Bueno, and S. Urrea, New Physics searches using ProtoDUNE and the CERN SPS accelerator, arXiv:2304.06765.
  * (25) B. Batell, W. Huang, and K. J. Kelly, Keeping it Simple: Simplified Frameworks for Long-Lived Particles at Neutrino Facilities, arXiv:2304.11189.
  * (26) J.-L. Tastet, SensMC, https://github.com/JLTastet/SensMC, 2023.
  * (27) K. Bondarenko, A. Boyarsky, O. Mikulenko, R. Jacobsson, and M. Ovchynnikov, Towards the optimal beam dump experiment to search for feebly interacting particles, arXiv:2304.02511.
  * (28) SHiP Collaboration, C. Ahdida et al., Sensitivity of the SHiP experiment to Heavy Neutral Leptons, JHEP 04 (2019) 077, [arXiv:1811.00930].
  * (29) SHiP Collaboration, C. Ahdida et al., Sensitivity of the SHiP experiment to dark photons decaying to a pair of charged particles, Eur. Phys. J. C 81 (2021), no. 5 451, [arXiv:2011.05115].
  * (30) SHiP Collaboration, H. Dijkstra and T. Ruf, Heavy Flavour Cascade Production in a Beam Dump, Dec, 2015.  CERN-SHiP-NOTE-2015-009.
  * (31) I. Boiarska, K. Bondarenko, A. Boyarsky, V. Gorkavenko, M. Ovchynnikov, and A. Sokolenko, Phenomenology of GeV-scale scalar portal, JHEP 11 (2019) 162, [arXiv:1904.10447].
  * (32) MATHUSLA Collaboration, C. Alpigiani et al., Recent Progress and Next Steps for the MATHUSLA LLP Detector, in 2022 Snowmass Summer Study, 3, 2022.  arXiv:2203.08126.
  * (33) M. Ovchynnikov, T. Schwetz, and J.-Y. Zhu, Dipole portal and neutrinophilic scalars at DUNE revisited: The importance of the high-energy neutrino tail, Phys. Rev. D 107 (2023), no. 5 055029, [arXiv:2210.13141].
  * (34) P. Ballett, T. Boschi, and S. Pascoli, Heavy Neutral Leptons from low-scale seesaws at the DUNE Near Detector, JHEP 03 (2020) 111, [arXiv:1905.00284].
  * (35) F. Kling and P. Quílez, ALP searches at the LHC: FASER as a light-shining-through-walls experiment, Phys. Rev. D 106 (2022), no. 5 055036, [arXiv:2204.03599].
  * (36) DUNE Collaboration, V. Hewes et al., Deep Underground Neutrino Experiment (DUNE) Near Detector Conceptual Design Report, Instruments 5 (2021), no. 4 31, [arXiv:2103.13910].
  * (37) FASER Collaboration, A. Ariga et al., FASER: ForwArd Search ExpeRiment at the LHC, arXiv:1901.04468.
  * (38) FASER Collaboration, H. Abreu et al., Technical Proposal: FASERnu, arXiv:2001.03073.
  * (39) J. L. Feng et al., The Forward Physics Facility at the High-Luminosity LHC, J. Phys. G 50 (2023), no. 3 030501, [arXiv:2203.05090].
  * (40) CHARM Collaboration, F. Bergsma et al., A Search for Decays of Heavy Neutrinos, Phys. Lett. B 128 (1983) 361.
  * (41) BEBC WA66 Collaboration, H. Grassler et al., Prompt Neutrino Production in 400-GeV Proton Copper Interactions, Nucl. Phys. B 273 (1986) 253–274.
  * (42) P. Ilten, Y. Soreq, M. Williams, and W. Xue, Serendipity in dark photon searches, JHEP 06 (2018) 004, [arXiv:1801.04847].
  * (43) M. W. Winkler, Decay and detection of a light scalar boson mixing with the Higgs boson, Phys. Rev. D 99 (2019), no. 1 015018, [arXiv:1809.01876].
  * (44) K. Bondarenko, A. Boyarsky, D. Gorbunov, and O. Ruchayskiy, Phenomenology of GeV-scale Heavy Neutral Leptons, JHEP 11 (2018) 032, [arXiv:1805.08567].
  * (45) B. Döbrich, J. Jaeckel, and T. Spadaro, Light in the beam dump - ALP production from decay photons in proton beam-dumps, JHEP 05 (2019) 213, [arXiv:1904.02091]. [Erratum: JHEP 10, 046 (2020)].
  * (46) D. Aloni, Y. Soreq, and M. Williams, Coupling QCD-Scale Axionlike Particles to Gluons, Phys. Rev. Lett. 123 (2019), no. 3 031803, [arXiv:1811.03474].
  * (47) S. Chakraborty, M. Kraus, V. Loladze, T. Okui, and K. Tobioka, Heavy QCD axion in b→s transition: Enhanced limits and projections, Phys. Rev. D 104 (2021), no. 5 055036, [arXiv:2102.04474].
  * (48) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, [arXiv:1405.0301].
  * (49) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, FeynRules 2.0 - A complete toolbox for tree-level phenomenology, Comput. Phys. Commun. 185 (2014) 2250–2300, [arXiv:1310.1921].
  * (50) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer, and T. Reiter, UFO - The Universal FeynRules Output, Comput. Phys. Commun. 183 (2012) 1201–1214, [arXiv:1108.2040].
  * (51) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, An introduction to PYTHIA 8.2, Comput. Phys. Commun. 191 (2015) 159–177, [arXiv:1410.3012].
  * (52) J. Carvalho, Compilation of cross sections for proton nucleus interactions at the HERA energy, Nucl. Phys. A 725 (2003) 269–275.
  * (53) V. Shtabovenko, R. Mertig, and F. Orellana, FeynCalc 9.3: New features and improvements, Comput. Phys. Commun. 256 (2020) 107478, [arXiv:2001.04407].
  * (54) SHiP Collaboration, C. Ahdida, R. Albanese, A. Alexandrov, et al., SHiP Experiment - Progress Report, tech. rep., CERN, Geneva, 2019.
  * (55) I. Boiarska, K. Bondarenko, A. Boyarsky, S. Eijima, M. Ovchynnikov, O. Ruchayskiy, and I. Timiryasov, Probing baryon asymmetry of the Universe at LHC and SHiP, arXiv:1902.04535.
  * (56) D. Gorbunov, E. Kriukova, and O. Teryaev, Scalar decay into pions via Higgs portal, arXiv:2303.12847.
  * (57) J. L. Feng, I. Galon, F. Kling, and S. Trojanowski, Dark Higgs bosons at the ForwArd Search ExpeRiment, Phys. Rev. D 97 (2018), no. 5 055034, [arXiv:1710.09387].
  * (58) F. Bezrukov and D. Gorbunov, Light inflaton Hunter’s Guide, JHEP 05 (2010) 010, [arXiv:0912.0390].
  * (59) C. Antel et al., Feebly Interacting Particles: FIPs 2022 workshop report, arXiv:2305.01715.
  * (60) T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider, Phys. Rev. C 92 (2015), no. 3 034906, [arXiv:1306.0121].
  * (61) A. Berlin and F. Kling, Inelastic Dark Matter at the LHC Lifetime Frontier: ATLAS, CMS, LHCb, CODEX-b, FASER, and MATHUSLA, Phys. Rev. D 99 (2019), no. 1 015021, [arXiv:1810.01879].
  * (62) J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review 59 (2017), no. 1 65–98, [https://doi.org/10.1137/141000671].
  * (63) D. Bachrathy and G. Stépán, Bisection method in higher dimensions and the efficiency number, Periodica Polytechnica Mechanical Engineering 56 (2012), no. 2 81–86.