FORESEE: FORward Experiment SEnsitivity Estimator for the LHC and future hadron colliders

We introduce a numerical package FORward Experiment SEnsitivity Estimator, or FORESEE, that can be used to simulate the expected sensitivity reach of experiments placed in the far-forward direction from the proton-proton interaction point. The simulations can be performed for $14$ TeV collision energy characteristic for the LHC, as well as for larger energies: $27$ and $100$ TeV. In the package, a comprehensive list of validated forward spectra of various SM species is also provided. The capabilities of FORESEE are illustrated for the popular dark photon and dark Higgs boson models, as well as for the search for light up-philic scalars. For the dark photon portal, we also comment on the complementarity between such searches and dark matter direct detection bounds. Additionally, for the first time, we discuss the prospects for the LLP searches in the proposed future hadron colliders: High-Energy LHC (HE-LHC), Super proton-proton Collider (SppC), and Future Circular Collider (FCC-hh).


I. INTRODUCTION
The Large Hadron Collider (LHC) has proven to be a powerful tool for studying both new physics and the Standard Model (SM), with its most remarkable achievement related to the discovery of the Higgs boson [1,2]. While further such outstanding signatures of new particles are much awaited, the data collected in the LHC have also already been used to constrain many beyond the Standard Model (BSM) scenarios in a tremendous number of phenomenological analyses. This became feasible due to the development of convenient modeling tools that were made available to a wider community.
This relentless quest for new physics will not only be continued in the upcoming LHC Run 3, but it will also be further extended to better encompass growing interest in searches for new light and long-lived particles (LLPs), cf. Refs [3,4] for recent reviews. While the traditional BSM physics programs target new particles with masses typically above O(10 GeV) up to a few TeV scale, lighter new physics species could become more accessible with modified experimental strategies [5]. In particular, this observation has lead to the establishment of a new direction in the BSM searches at the LHC in its far-forward region, as originally proposed [6][7][8][9] for the FASER experiment [10][11][12]. Further such opportunities could be explored in the high-luminosity LHC (HL-LHC) era with an expanded research agenda of the Forward Physics Facility (FPF) [13].
In order to fully exploit the relevant physics potential, the experimental efforts should be supplemented with a comprehensive program of theoretical and phe- * felixk@slac.stanford.edu † strojanowski@camk.edu.pl nomenological studies. To facilitate this, we introduce a numerical package, namely the FORward Experiment SEnsitivity Estimator, or FORESEE, which could be used to obtain the expected sensitivity reach for BSM models in various far-forward experiments. 1 The package allows one to perform quick but accurate simulations for selected popular BSM simplified models. This can be done for user-defined experimental geometries and basic cuts applied to the visible signal. The package also provides a set of useful numerical data, including e.g. the far-forward spectra of light mesons, that can easily be accessed and employed in separate studies estimating the new physics sensitivity reach in other BSM scenarios. We illustrate below the capabilities of FORESEE for the popular dark photon and dark Higgs boson models, as well as for the model with a hadrophilic dark scalar with the dominant couplings to the up quarks [14,15]. In addition, for the first time, we study the possibility to perform such far-forward searches in the future 27 TeV and 100 TeV hadron colliders.
This study is organized as follows. In Sec. II, we introduce the FORESEE package. We briefly present the validated far-forward spectra of hadrons and other SM particles that can be used for BSM modeling. We then discuss the capabilities of the automated numerical tool that we provide in studying the sensitivity reach of far-forward experiments. In Sec. III, the FORESEE is employed to analyze specific examples of the far-forward LLP searches both in the LHC and future hadron colliders. We conclude in Sec. IV. Distribution of π 0 (left) and B 0 (right) mesons in the forward hemisphere in the (θ,p) plane, where θ and p are the meson's angle with respect to the beam axis and momentum, respectively. We use 20 bins per decade and present the spectra obtained for 14 TeV pp collision energy. The π 0 spectrum is obtained with EPOS-LHC [16], while the B-meson one with Pythia [17] with the Monash tune [18]. The diagonal black dashed lines highlight the characteristic transverse momentum scale pT ∼ ΛQCD ∼ 250 MeV for pions and pT ∼ mB for B mesons. The angular acceptances for the FASER and FASER 2 experiments to take data during LHC Run 3 and in the HL-LHC era, respectively, are indicated by the vertical black dotted lines.

II. FORESEE PACKAGE
A. Forward spectra of SM particles Light new physics particles can be efficiently produced in the far-forward region of the LHC due to at least several different production mechanisms, cf. Refs [11,19,20] for further discussion. A prominent role among them is played by rare decays of light mesons that are abundantly created in pp collisions. Notably, one expects e.g. about N π 0 ∼ 4 × 10 17 neutral pions to be produced during Run 3. In particular, high-energy mesons typically travel far-forward at angles θ M ∼ p T,M /p M 1 with respect to the beam axis. Here p T,M ∼ m M ∼ O(100 MeV − a few GeV) is the meson's characteristic transverse momentum of order the meson mass, and p M ∼ TeV is the meson's total momentum. Importantly, such meson decays are typically the dominant LLP production channels as long as they are kinematically available.
Properly predicting the meson spectra is then essential for modeling new physics effects in far-forward experiments. This can be done using dedicated hadronic interaction models, designed to describe inelastic collisions at both particle colliders and cosmic ray experiments. Those models have greatly improved in recent years, partially due to the dedicated forward physics program at the LHC [21]. In the left panel of Fig. 1, we present exemplary such spectra of neutral pions and B mesons in the (θ M , p M ) plane. These were obtained, respectively, with the EPOS-LHC [16] event generator implemented in the CRMC simulation package [22], and by using Pythia 8 [17,23] with the Monash tune [18]. As expected, the spectra are concentrated along the p T,M ∼ m M lines.
In the FORESEE package, these spectra are then used to generate forward LLP flux due to the relevant meson decays. Since in some applications it might be beneficial to access only these spectra and employ them in independent simulations, we provide them in a simple format in separate text files for different mesons and three pp collision energies: 14 TeV, 27 TeV, and 100 TeV. The relevant files can be found in the package's files/hadrons directory. The meson spectrum is provided in tables with the consecutive columns corresponding to the bin position in (log 10 θ M , log 10 (p M /GeV)) and the weights of each bin (3rd column) in units of pb/bin. Here θ M in radians is the angle with respect to the beam axis. To obtain the number of mesons per bin, the bin weights should be multiplied by the relevant integrated luminosity (e.g., 150 fb −1 for LHC Run 3, 3000 fb −1 for HL-LHC). The number in each of the filenames is the meson ID in the MC particle numbering scheme [24]. The production cross section is given for the forward hemisphere only. We similarly treat selected other SM species listed below.
We provide the spectra for the following light mesons: π 0 , π ± , K ± , K L , K S , K * 0 , K * ± , η, η , ω, ρ and φ, for the baryons: n, p and Λ and their anti-particles, as well as for photons γ. The user can access the results obtained for three different MC generators that should also be referred to when using the spectra, namely EPOSLHC [16], QGSJET II-04 [25], and SIBYLL 2.3c [26,27] as implemented in CRMC [22]. For heavier SM species, we employ SIBYLL and Pythia 8 [17,23] with the Monash tune [18] for charmed hadrons (D 0 , D + , D s , Λ c , and c,c quarks), while we use Pythia 8 for bottom hadrons (B 0 , B + , B s , B c , Λ b and b,b quarks), heavy gauge bosons W and Z, and for the SM Higgs boson h. In addition, we also provide the spectra of J/ψ, ψ(2S), and Υ(1, 2, 3S) mesons for 14 TeV LHC, see Ref. [28] for further discussion and validation of these spectra. In Table I, we summarize the SM spectra available in FORESEE and the relevant simulation tools used to obtain them.

B. Sensitivity Reach Estimator
On top of providing the users with the meson spectra files, the FORESEE package can also be used as a standalone simulation tool to study the sensitivity reach of far-forward experiments at the LHC and in future hadron colliders. The simulations can be performed for all the three aforementioned pp collision energies. Below, we briefly describe the consecutive steps of such simulations. Further instructions on how to use the code are provided in the package in tutorial jupyter notebooks.
In the first step of the simulation, the user must define the model by specifying i) the LLP production rates, ii) their lifetimes, and iii) the LLP decay branching fractions. The last information is optional and, if it is not defined, it is assumed that the experiment is sensitive to all the LLP decays happening inside the detector.
The FORESEE package supports three different production modes: SM particle decays, mixing with SM particles, and direct production. In most models, the main production modes of LLPs are the decays of SM particles. For 2-body decays and unpolarized beams, the kinematics are fully specified by the masses of the parent and daughter particles, in which case one only needs to provide the total decay branching fraction. There is also the option for 3-body decays, p 0 → p 1 p 2 p 3 , with p 3 being the LLP. Here ones needs to provide the differential branching fraction dBR/(dq 2 d cos ϑ), where q 2 = (p 2 +p 3 ) 2 and ϑ is the angle between p 3 in the rest frame of p 2 + p 3 , and the direction of p 2 + p 3 in the rest frame of p 0 . When modeling the production of LLPs in decays of long-lived mesons (charged pions, charged/neutral kaons), we additionally require the mesons to decay before they hit the beampipe at radius r = 5 cm, the neutral particle absorber TAN at a distance of z = 140 m from the pp Interaction Point (IP), or the inner triplet quadrupole absorber TAS at z = 20m. Charged particles are always required to decay before the first magnets at z = 20 m.
In addition, LLPs can also be produced via their mixing with the SM mesons. Examples are the mixing of axion-like particles with the SM pseudoscalar mesons or dark photons with the SM vector bosons, which were discussed in Ref. [11] and Ref. [29], respectively. Here we assume that the LLP production rate is related to the SM particle production rate via σ(LLP ) = κ 2 × σ(SM ), where κ describes the mixing which must be provided by the user.
Finally, LLPs can also be produced directly. Examples for this include dark photon production via Bremsstrahlung or Drell-Yan production (see Ref. [6] and Ref. [29] for details on these production modes). In this case, the user must provide the full two-dimensional LLP spectra for different LLP masses in the same format as previously used for the SM particles.
The LLP lifetime, cτ , and decay branching fractions need to be provided by the user in form of a table. The user can choose either a one-dimensional or a twodimensional parameterization. In the first case, the lifetime is given for the specific value of the coupling constant g * as a function of the varying LLP mass m. For different values of the coupling constant g, one evaluates the lifetime as cτ (m, g) = cτ (m, g * )g 2 * /g 2 . In the two-dimensional parameterization, cτ (m, g) is provided directly as a function of both m and g. In both cases, the code then interpolates the lifetime throughout the parameter space. The same options are available for the branching fractions.
After specifying the model, FORESEE generates the twodimensional LLP spectra in terms of (θ LLP , p LLP ). For the next step of the analysis, the user can specify the i) distance L between the IP and experiment, ii) acceptance in terms of the LLPs momentum and position, iii) luminosity, iv) production channels to consider, and v) allowed LLP decay channels. FORESEE then counts the number of signal events that pass the selection criteria.
In particular, the package allows the user to define the detector position and geometry, as well as the analysis cuts, by restricting the allowed position and momentum 3-vectors of the LLP. The reach plots are automatically generated for a userdefined number of BSM signal events. In the plots, current bounds are also shown, as well as the expected reach for selected proposed future experiments. For the detector under study, the default number of signal events, which defines the sensitivity line to be drawn, is N ev = 3. This appropriately estimates the reach in searches with negligible SM backgrounds. It should be stressed, however, that FORESEE only estimates the BSM signal rate. Depending on the detector location and design, there could be a large number of background events that could affect the searches for new physics. For example, in detectors with a material-filled decay volume, scattering signatures from neutrinos, photons, or neutral hadrons could mimic LLP decays if they are not carefully rejected in the analysis. Alternatively, locations not well shielded from the LHC beam and infrastructure could suffer from significant backgrounds from the beam debris, which would have to be taken into account in the modeling. In order to allow the user to include such effects when estimating the BSM sensitivity reach of a given experiment, we leave freedom to define the aforementioned number of signal events N ev to be presented in the plots.
Further developments of the package are planned for the future that will add more popular LLP models, the relevant production and decay modes, as well as will expand related modeling capabilities.

III. PHYSICS EXAMPLES
In this section, we illustrate the capabilities of FORESEE by presenting the expected sensitivity reach plots for selected specific BSM scenarios and the LLP decay signature. We first focus on the FASER and FASER 2 experiments at the LHC and then move to discuss the prospects for LLP searches in the far-forward region of the future hadronic colliders.
When presenting our results, we assume that background in the searches for decays of high-energy LLPs, typically with E LLP 100 GeV, can be reduced to negligible levels. The possibility of such a background reduc-tion has been studied in detail for the FASER experiment to take data during LHC Run 3 [10,12]. In the case of other detectors under study, we expect that a similar suppression will be possible, as they would be shielded from the pp collision point by at least O(100 m) of the rock and other elements of the collider infrastructure.

A. Models
For illustration, we will present the results for the models predicting the existence of sub-GeV dark photons A , dark Higgs bosons φ, and the up-philic scalars S. The corresponding interaction Lagrangians are given by: wheref f represents the sum over all the SM fermions and h is the SM Higgs boson. Above, we also define the relevant coupling constants: the kinetic mixing parameter for A , the mixing angle θ and trilinear coupling λ for the dark Higgs boson, and the g u coupling for the upphilic scalar. Further details of the far-forward modeling for these scenarios can be found in Ref. [6] for A s and in Ref. [7] for the dark Higgs boson, while the case of the up-philic scalar is described below. These example models are characterized by different dominant production or decay modes of LLPs. In particular, light dark photons with m A O(100 MeV) are dominantly produced in π 0 decays, the up-philic scalars come from η decays, and the dark Higgs bosons are mostly generated in rare decays of B mesons. Notably, while they typically lead to LLP decays into two opposite-charge SM particles, the light up-philic scalars can also have the leading decay channel into two photons, cf. Refs [9,30] for other examples of far-forward searches for LLPs with the dominant di-photon decay channel.

B. Dark photon
In the left panel of Fig. 2, we show the expected sensitivity reach of the FASER and FASER 2 experiments in the search for the vanilla dark photons, cf. Eq. (1). The assumed detector geometries and the collider luminosities are given in Table II. The total sensitivity reach lines obtained with FORESEE and presented in the plot, reproduce the results obtained previously by the FASER collaboration [11].
In addition, we present in the plot the expected reach of FASER 2 assuming that only the simplest A → e + e − decay channel is considered in the analysis. As expected, this has no impact on the reach for the dark photon mass below the di-muon threshold, while it moderately reduces the sensitivity for m A 500 MeV. Here, possible A decays into hadronic final states start to play the dominant role and they would have to be taken into account in the analysis. The relevant hadronic branching fraction is dominated by the decays into two charged pions, A → π + π − , while in narrow ranges of the A mass around the ω or φ resonances decays with neutral pions or kaons in the final state become also important [31].
The sensitivity reach can also be degraded if the detector is not positioned along the beam collision axis, cf. Refs [11], also Ref. [32] for the recent discussion about the SND@LHC experiment. However, this effect is not substantial as long as the displacement is not too large and the beam axis crosses the detector volume. In order to illustrate this, in the left panel of Fig. 2 we also present the corresponding sensitivity reach for the dark photon model of a FASER 2-like detector with the radius R = 1 m assuming that its center has been shifted by 1 m off the beam axis. The geometrical acceptance of the detector can be varied freely in FORESEE.
One of the most important motivations for searching for LLPs is the role they can play as mediators between the SM and dark matter (DM) particles χ. In particular, this allows for obtaining the correct value of the relic density, Ω total χ h 2 0.12 [48], via the popular thermal production mechanism [49][50][51]. Since such scenarios often predict both the mediator and DM masses to be in the sub-GeV range, this could naturally lead to complementarity between searches for light mediator particles and those targeting DM direct detection (DD) signatures [52].
In the left panel of Fig. 2, we show the example thermal relic contour obtained for the scenario with the A mediator and light complex scalar DM described by the following Lagrangian where D µ = ∂ µ − ig D A µ and we denote α D = g 2 D /(4π). For the line presented in the plot, we fix the dark coupling constant at a value α D = 0.1 and we assume that the mass ratio between the two dark species is always equal to m χ /m A = 0.6 for varying dark photon mass.
This choice of the mass ratio guarantees that the dark photon decays visibly into the SM species and can then be searched for in FASER and FASER 2 (see Ref. [53] for the discussion about invisibly decaying forward-going dark photons). On the other hand, we require χ to be lighter than the dark photon. In this case, the DM relic density is set by annihilations into the SM fermions, χχ → A * → ff , and it is not suppressed by otherwise (if kinematically allowed) very efficient annihilations in the dark sector, χχ → A A . The particular value of the benchmark mass ratio used in Fig. 2 has been chosen to minimize the possible impact on Ω χ h 2 of the forbidden DM annihilations χχ → A A [54], the dark photon initial state radiation χχ → A A * → A [55], and the resonant annihilations into the SM particles, χχ → A → , with m A 2m χ , cf. Refs [56][57][58] for recent discussion. Last but not least, we note that the annihilation cross section for scalar DM is p-wave suppressed, which allows one to avoid stringent bounds from the Cosmic Microwave Background (CMB) radiation [48,59].
As can be seen in the plot, a small FASER detector during Run 3 will already be able to probe an important part of the allowed region of the parameter space which corresponds to the correct value of the DM relic density. Notably, in this scenario, the FASER reach corresponds to the most cosmologically relevant region of the parameter space, in which χ DM is not overproduced in the early Universe. Instead, for the region below the thermal relic line, the correct DM abundance cannot be reproduced in this simple model without introducing further modifications in the allowed DM interactions or in the standard cosmological scenario.
The role of far-forward search for A in FASER at the LHC can be further emphasized by studying current constraints on such DM particles from DD searches in the underground detectors. In the left panel of Fig. 2, we present with a light-gray shaded region current such dominant bounds obtained by the Xenon10 [60] and Xenon1T [61,62] detectors. In the plot, we combine bounds from the DM-electron scatterings following Ref. [63], as well from the searches for DM-nucleus interactions that employ the Migdal effect [64]. When presenting these bounds, we also assume that Ω χ h 2 Ω total χ in the entire relevant region of the parameter space.
After taking into account the DM DD searches, there remains only a narrow allowed region in the parameter space of the A /χ model under study, in which we do not expect DM to be thermally overproduced, i.e. we require Ω thermal   (m A , ) and (mS, gu) planes, respectively, obtained for the FASER (brown solid lines) and FASER 2 (red) detectors. In the left panel, we also show for the dark photon model the expected FASER 2 sensitivity to only A → e + e − decay channel (magenta) and to all the possible visible decays but happening in a displaced detector by a 1 m distance off the beam axis (orange). In both panels, we present previous bounds with dark gray-shaded regions and future sensitivity (colorful dashed lines) for selected other searches, as indicated in the plots. In the left panel, we also show with light-gray color current constraints on the dark photon parameter space from dark matter direct detection searches. In this case, we assume the complex scalar χ DM model with the fixed mass ratio mχ/m A = 0.6 and the dark coupling constant αD = 0.1. See the text for more details and references for the bounds.
We stress that the DM DD bounds could be evaded e.g. for the inelastic DM scenario with the mass splitting between the two dark sector particles sufficient to suppress DM DD rates. Importantly, the presence of such a mass splitting between the two dark species would not affect the FASER sensitivity reach, unless it is large enough so that decays of the excited state can also contribute to the visible signal in the detector, cf. Refs [19,29] for further similar discussions.

C. Up-philic scalar
We next consider the search for the light scalars coupled dominantly to the up quarks, cf. Eq. (3). While typical LLP studies focus on the Yukawa-like interactions of BSM scalars that resemble the SM Higgs boson, coupling light scalars to only selected SM mass eigenstates leads to distinct phenomenology. In particular, by restricting the couplings to only up quarks, one can avoid stringent constraints from B meson decays and other flavor bounds, while simultaneously keeping the production rate of such scalars in pp collisions at a high level [14,15].
The dominant production modes of the up-philic scalars are in rare decays of η and η mesons, η, η → π 0 S, while kaon and B meson decays are suppressed in this case, differently from the light dark Higgs model. The scalars S decay preferentially into hadronic final states, if kinematically allowed, with the dominant such decay channel into a pion pair. Instead, for lower scalar masses, the lifetime of S is driven by loop-induced decays into photons and it, therefore, becomes significantly larger.
The drop in the lifetime of S for m S growing above the pion threshold leads to specific features in the reach plot presented in the right panel of Fig. 2. For m S < 2m π 0 , the expected sensitivity is limited by the large S lifetime that allows LLPs to easily overshoot the detector. As a result, only values of the coupling constant above certain threshold, g u 10 −5 (10 −6 ), can be probed in FASER (FASER 2). In addition, the relevant region of the parameter space has already been excluded by past searches in CHARM [66], E137 [35], and MAMI [67] experiments. For completeness, following Ref. [15], we also show in the plot the constraints from the energy-loss rate in the core of SN1987A and from the number of effective degrees of freedom that could affect the Big Bang Nucleosynthesis (BBN).
On the other hand, for heavier S with m S > 2m π 0 , we can take advantage of large LHC energies and boost factors to improve the previous bounds on more short-lived up-philic scalars. The corresponding expected sensitivity of both FASER and FASER 2 is, however, concentrated in the limited region of the parameter space. It corresponds to the dark scalar masses lying in the range 2 m π < m S < m η − m π . The far-forward searches at the LHC will allow one to probe the currently not excluded values of the dark coupling constant in between the lower CHARM bound and the upper constraints from the KLOE experiment [68]. Notably, the latter search was sensitive to S → π + π − decays and, therefore, it did not constrain dark scalar masses below the corresponding threshold, even though possible decays to lighter neutral pions would keep the S lifetime relatively low in the re- maining small mass window, 2 m π 0 < m S < 2 m π ± . Instead, when estimating the FASER and FASER 2 reach, we take into account both decays into charged and neutral pions, with the latter leading to a striking 4γ signature in the detector. In the figure, we also show for comparison the expected future sensitivity of the proposed SHiP and REDTOP [69] experiments following Ref. [15].

D. Future Colliders
While we have so far focused on the LHC, the idea to search for LLP decay signatures in the far-forward region of the hadron colliders could also be employed in future such experimental facilities beyond the HL-LHC era. It is then useful to briefly analyze the relevant detection prospects in connection to the proposed such colliders, see, e.g., Ref. [70,71] for recent reviews. In particular, we will discuss the sensitivity reach of the detector similar in size and position to FASER 2, which could operate at the 27 TeV pp collision energy in the era of the High-Energy LHC (HE-LHC) [72]. We also present the sensitivity reach for 100 TeV energy. In this case, we choose for illustration a possible forward detector location based on the proposed concepts of the proton-proton Future Circular Collider (FCC-hh) [73] and the Super Proton-Proton Collider (SppC) [74]. We give the relevant parameters of the detectors in Table II.
While for the HE-LHC we employ the same location as for FASER and FASER 2, the 100 TeV colliders require shifting the detector to a larger distance L from the IP. This is due to a smaller curvature of the beam pipe characteristic for larger colliders. In particular, to allow for sufficient shielding from the beam-induced background, we require the minimal distance between the far-forward detector and the (already curved) beam tunnel to be of order 5 m. We then expect that the minimal corresponding values of L are equal to 1.1 km for FCC-hh and 865 m for SppC, respectively. In the following, for simplicity, we assume L = 1 km for the detector operating at 100 TeV pp collision energy.
The relevant sensitivity reach lines for FASER 2 during HL-LHC, HE-LHC and FCC-hh are shown in Fig. 3 for the dark photon and dark Higgs models, cf. Eqs. (1) and (2), respectively. As shown in the left panel, increasing the pp collision energy and the boost factor of the produced LLPs will improve the expected sensitivity reach in the dark photon model. This is especially true for larger values of both and m A that corresponds to the smaller A lifetime. It is worth stressing, though, that for the vanilla dark photon model this improvement is not significant. In fact, if fully exploited, the far-forward searches in FASER and FASER 2, along with other experimental probes, could almost saturate the relevant detection prospects already during the LHC and HL-LHC eras.
On the other hand, the detectors operating in both the HE-LHC and FCC-hh/SppC could have a more important impact on constraining the allowed parameter space of the dark Higgs model. We illustrate this in the right panel of Fig. 3. In the plot, we also show with a gray-shaded region current bounds. The dominant ones come from the searches for dark scalars in kaon decays K → πφ by NA62 [75] and BNL-E949 [76], B-meson decays B → Kφ by LHCb [77,78], and long-lived particle searches at CHARM [79], LSND [80] and Micro-BooNE [81]. We also show the limits from astrophysical observations of SN 1987A and the BBN bounds, following Ref. [79].
For a light scalar φ, which can be efficiently produced in rare B meson decays, the search performed at 100 TeV energies will improve the sensitivity reach by an additional order of magnitude in the φ − h mixing angle θ in comparison with the FASER 2 experiment operating during HL-LHC. For the lowest values of θ, the dominant production mode of φ in this region can be due to the B meson decays into pairs of dark scalars via an off-shell Higgs boson, b → sh * → sφφ. This production mode remains relevant if the corresponding SM Higgs branching fraction is larger than O(0.1%), as dictated by the trilinear coupling λ in Eq. (2). This, however, only affects a narrow region in the covered parameter space, while the majority of the observed reach in this dark scalar mass range is not sensitive to λ.
The situation changes for dark scalar masses larger than the B meson mass. Here, FASER 2 in the HL-LHC era will have almost no sensitivity reach, besides a very narrow region in the parameter space, in which m φ m h /2 [82]. In this case, the dark scalars are pairproduced in decays of on-shell SM Higgs bosons, h → φφ. This production mode becomes even more important for similar such detectors operating at future colliders. The pp collision energy then becomes large enough so that even the SM Higgs bosons with m h 125 GeV are sufficiently boosted to become relevant for far-forward searches. In particular, in the right panel of Fig. 3, we assume BR(h → φφ) = 5% when obtaining the reach lines. In this case, the sensitivity reach could be extended to the values of m φ spanning the entire mass range in between the B meson mass and m φ m h /2. This will also correspond to very low values of the mixing angle, θ ∼ 10 −9 . Again, the impact of direct φ production in decays of the SM Higgs bosons will become much suppressed if the corresponding branching fraction falls below 0.1%.

IV. CONCLUSION
The upcoming LHC Run 3 will soon begin the next chapter in the continued efforts towards discovering signatures of new physics. Besides extending the existing experimental strategies, it will also initiate a new direction in BSM searches at the LHC in its far-forward region. In this case, when properly placed, even relatively small-size detectors can offer excellent discovery prospects. Notably, this has been recently demonstrated with the first observation of neutrino candidate events at the LHC [83].
The most commonly discussed strategy to detect new forward-going light and long-lived particles produced abundantly at the pp interaction point is via their displaced decays into SM species. These will be searched for in the distant FASER experiment. In the more distant future, a similar research agenda could be promoted to an even higher level by detectors taking data during the HL-LHC era that will be placed in the Forward Physics Facility.
This experimental program could much benefit from further phenomenological and theoretical analyses of new physics models that could be probed this way. In this study, we have introduced a numerical package FORward Experiment REach Estimator, or FORESEE, that allows one to estimate the sensitivity reach of the far-forward experiments operating at 14, 27, or 100 TeV pp collision energies in their searches for highly-displaced LLP decays. In addition, in the package, the validated far-forward spectra of light mesons and many other SM species have been made available to facilitate their use in separate studies, cf. Sec. II A for the discussion about the simulation details and the relevant references.
We have first illustrated the capabilities of the package by presenting the expected reach of the FASER and FASER 2 experiments in the dark photon and up-philic scalar models. In the first case, we have also discussed the possible complementarity between the FASER search for visibly decaying dark photons and DM direct detection experiments, assuming a simple model with complex scalar DM and the A mediator. We stress that such a complementarity could best be seen for only a limited mass range of both dark species, m A /2 m χ m A . The cosmologically viable such scenarios are further constrained by requiring that the DM is not thermally overproduced in the early Universe. Interestingly, however, these conditions can be satisfied in the simple benchmark case precisely in the region of the parameter space of the model that will be covered by FASER during Run 3.
Last but not least, for the first time, we have also analyzed the discovery potential of possible far-forward experiments to operate in the future hadron colliders: HE-LHC, SppC, and FCC-hh. To this end, we have focused on the popular dark photon and dark Higgs boson models. We show that, while in the former case the capabilities of the far-forward searches can basically be saturated already during the HL-LHC phase, the latter BSM scenario could benefit from employing the increased pp collision energy. This is, however, mainly due to possible decays of the on-shell Higgs bosons. The precise impact of employing the future colliders on the dark Higgs sensitivity reach will, therefore, depend on earlier measurements of the relevant invisible Higgs branching fraction.
Astrophysics Science and Technology Centre, carried out within the International Research Agendas programme of the Foundation for Polish Science financed by the European Union under the European Regional Development Fund. S.T. is also supported in part by the Polish Ministry of Science and Higher Education through its scholarship for young and outstanding scientists (decision no 1190/E-78/STYP/14/2019).