text:
\section{Introduction}

The quest for beyond the Standard Model (BSM) physics is undoubtedly among the biggest challenges in contemporary physics. It has typically been driven by the lack of answers to some fundamental questions, such as understanding the unification of the forces of nature or the mechanism of baryogenesis, as well as by hints of new physics appearing in anomalous experimental observations. The latter become especially intriguing when the known Standard Model (SM) background (BG) can be reduced to negligible levels so that even the observation of a few events can be attributed to BSM effects.

The flavor-changing rare decays of kaons are exceptionally promising examples of such discovery channels given their highly suppressed SM rates. In particular, the KOTO collaboration has recently reported a possible observation of three anomalous events in the search for $K_L\to\pi^0\nu\bar{\nu}$~\cite{KOTOanomaly}, in which no SM signal was expected at the current level of the experimental sensitivity. If taken as a signal, this observation would indicate that the relevant branching fraction exceeds the SM prediction by about two orders of magnitude~\cite{Kitahara:2019lws}. It has been quickly noticed that the discrepancy could be resolved by introducing a new light long-lived particle (LLP) with mass of order $100~\text{MeV}$, light enough to be produced in kaon decays. By properly adjusting the lifetime of the LLP, one can explain the neutral kaon anomaly while avoiding stringent bounds on such models from searches for charged kaon decays.

In this study, we illustrate how such proposed solutions to the KOTO anomaly could be thoroughly probed by the recently approved ForwArd Search ExpeRiment, or FASER~\cite{Feng:2017uoz, Ariga:2018zuc, Ariga:2018pin}. The FASER experiment will search for LLPs abundantly produced in the far-forward direction of the Large Hadron Collider (LHC) during its upcoming Run 3. Several features of its setup make it an ideal probe for light LLPs produced in kaon decays:
i) Given a large expected boost factor of light particles traveling towards FASER, and the location of the detector, the experiment's sensitivity is optimal for the precise combination of the lifetime and mass that are required to fit the anomaly. In particular, for high-energy LLPs with $E\sim 1~\text{TeV}$ that are produced at the ATLAS IP, the sensitivity reach of FASER is maximized for
\begin{equation}
\frac{\tau_{\textrm{LLP}}}{m_{\textrm{LLP}}} \sim \frac{0.1~\textrm{ns}}{100~\text{MeV}},\hspace{1cm}\textrm{(best FASER reach)}.
\end{equation}
ii) The FASER sensitivity additionally benefits from the $\mathcal{O}(150~\text{m})$ long part of the LHC beam pipe, in which forward-going neutral kaons can travel and decay before being absorbed, as well as the presence of additional LLP production modes at the LHC.
iii) The FASER detector has the capability to not only detect LLP decays into charged particles but also into photon pairs, due to its dedicated pre-shower detector. As a result, during Run 3 of the LHC, FASER can turn the three currently observed anomalous events at KOTO into up to $\mathcal{O}(10^4)$ spectacular high-energy LLP decay events with expected negligible BG.

The rest of this paper is organized as follows. In \cref{sec:KOTO} and \cref{sec:FASER} we provide more details about the KOTO anomaly and the FASER experiment, respectively. The FASER sensitivity reach in several selected BSM scenarios is presented in \cref{sec:results}. In \cref{sec:ALPs}, we further illustrate the possible interplay between FASER and searches for LLPs in kaon factories. We conclude in \cref{sec:conclusion}. Several useful expressions for BSM meson decay branching fractions are listed in \cref{app:BRs}.

\section{\label{sec:KOTO} KOTO experiment and the anomaly}

The KOTO experiment~\cite{Yamanaka:2012yma} has been designed to study decays of kaons into neutral pions and a neutrino/anti-neutrino pair, $K_L\rightarrow\pi^0\nu\bar{\nu}$. The kaons are produced in collisions of $30~\text{GeV}$ protons from the Japan Proton Accelerator Research Complex~\cite{Nagamiya:2012tma} (J-PARC) main ring accelerator with the target made out of gold. A high-intensity beam of neutral kaons produced at an angle of $16^\circ$ with respect to the proton beam line is created with the use of dedicated collimators. At a distance of about $24~\text{m}$ away from the target, the $\sim\!3$ m long vacuum chamber of the KOTO detector begins. Here, kaon decays are identified by detecting photons produced in prompt neutral pion decays.

The data collected by the KOTO collaboration in 2015 allowed them to obtain the leading bound on the branching fraction of the aforementioned decay process. At the $90\%$ confidence level (CL) it reads~\cite{Ahn:2018mvc}
\begin{equation}
\mathcal{B}_{\textrm{KOTO,bound}}(K_L\rightarrow\pi^0\nu\bar{\nu}) \leq 3\times 10^{-9}.
\label{eq:K0Lneutbound}
\end{equation}
Notably, \cref{eq:K0Lneutbound} is consistent with the SM prediction which remains about two orders of magnitude lower, $\mathcal{B}_\textrm{SM}(K_L\rightarrow\pi^0\nu\bar{\nu})= (3.4\pm 0.6)\times 10^{-11}$~\cite{Buras:2015qea}. A similar bound on the two-body decay $K^0\rightarrow\pi^0\,X$ has also been obtained, where $X$ is a stable or long-lived BSM bosonic particle,
\begin{equation}
\mathcal{B}_{\textrm{KOTO,bound}}(K_L\rightarrow\pi^0\,X) \lesssim 2.4\times 10^{-9}.
\end{equation}

In the subsequent analysis of the data from 2016-18~\cite{KOTOanomaly}, however, a total of four candidate events were found, only one of which had the properties of BG. If taken as a signal, the remaining three anomalous events point to the branching fraction which significantly exceeds the SM prediction. At the $95\%$ CL it is given by~\cite{Kitahara:2019lws}
\begin{equation}
\mathcal{B}_{\textrm{KOTO,anom.}}(K_L\rightarrow\pi^0\nu\bar{\nu}) = (2.1^{+4.1}_{-1.7})\times 10^{-9}.
\label{eq:K0Lneutanomaly}
\end{equation}

Importantly, in the SM the aforementioned neutral kaon decay mode is also related by the isospin symmetry to the value of the branching fraction of a charged kaon decay into a pion $\pi^+$ and a neutrino/anti-neutrino pair, which proceeds via the same parton level process $s\to d\nu\bar{\nu}$. The relevant so-called Grossman-Nir (GN) bound~\cite{Grossman:1997sk} reads
\begin{equation}
\mathcal{B}(K_L\rightarrow\pi^0\nu\bar{\nu}) < 4.3\,\mathcal{B}(K^+\rightarrow\pi^+\nu\bar{\nu}).
\label{eq:GNbound}
\end{equation}
The current upper limit on the $K^+$ decay branching fraction is obtained from the results of the E949 experiment, which observed several charged kaon decay events~\cite{Artamonov:2008qb,Artamonov:2009sz}, and it is given by $\mathcal{B}_{\textrm{E949,bound}}(K^+\rightarrow\pi^+\nu\bar{\nu})<3.35\times 10^{-10}$ at $90\%$ CL. This is consistent with the SM expectation $\mathcal{B}_\textrm{SM}(K^+\rightarrow\pi^+\nu\bar{\nu}) = (8.4\pm 1)\times 10^{-11}$~\cite{Buras:2015qea}. An improved analysis has also recently been preliminarily presented by the NA62 collaboration with the corresponding branching fraction measured to be at the $68\%$ CL~\cite{NA62bound},
\begin{equation}
\mathcal{B}_{\textrm{NA62}}(K^+\rightarrow\pi^+\nu\bar{\nu}) = 0.47^{+0.72}_{-0.47}\times 10^{-10},
\label{eq:NA62bound}
\end{equation}
and the upper bound at the $95\%$ CL given by $\mathcal{B}_{\textrm{NA62,bound}}(K^+\rightarrow\pi^+\nu\bar{\nu})<2.44\times 10^{-10}$. The measurements of the charged kaon decay branching fraction, together with the GN bound \cref{eq:GNbound}, lead to $(2-3)\sigma$ tension with the anomalous observation by the KOTO collaboration, cf. \cref{eq:K0Lneutanomaly}, depending on how the possible impact of new physics is taken into account~\cite{Kitahara:2019lws}. A number of BSM scenarios have been proposed to explain this discrepancy~\cite{Kitahara:2019lws, Fabbrichesi:2019bmo, Egana-Ugrinovic:2019wzj, Li:2019fhz,Dev:2019hho, Liu:2020qgx,Jho:2020jsa, Cline:2020mdt,He:2020jzn, Ziegler:2020ize, Liao:2020boe, He:2020jly, Gori:2020xvq, Hostert:2020gou, Datta:2020auq, Dutta:2020scq,Altmannshofer:2020pjb,Liu:2020ser, Archer-Smith:2020hqq}.

Among these scenarios, a prominent role is played by models predicting the existence of a new LLP with the mass $m_X\sim 100~\text{MeV}$, which can be produced in rare kaon decays, $K_L\rightarrow\pi^0\,X$~\cite{Kitahara:2019lws}. In particular, if such a light BSM particle has a lifetime of order $\tau_X\sim 0.1$~ns, it can be effectively stable and invisible in the search performed by the KOTO collaboration. At the same time, $X$ will typically decay inside the E949 and NA62 detectors, and it does not contribute to the measured branching fraction of $K^+_L\rightarrow\pi^+\nu\bar{\nu}$. This leads to an apparent violation of the GN bound when the results of these experiments are compared with each other, which is further referred to as the \textsl{lifetime solution} to the anomaly.

A very long-lived $X$ can also explain the anomaly while avoiding the bounds from past beam-dump searches. In addition, if its mass is close to the pion mass, it can escape detection in E949 and NA62 experiments by hiding in the SM BG~\cite{Fuyuto:2014cya}. Another distinct possibility is to allow $K_L$ decays into neutral dark states, $K_L\rightarrow X_1\,X_2$. Subsequent $X_i$ decays that produce photons can be detected in the KOTO electromagnetic calorimeter and resemble a neutral pion signature~\cite{Hostert:2020gou}. In this case, charged kaon three-body decays can be suppressed or even kinematically forbidden.

Below, we analyze the prospect of probing such selected BSM scenarios in the FASER experiment at the LHC.

\section{\label{sec:FASER}FASER and FASER~2 experiments}

The FASER experiment~\cite{Feng:2017uoz,Ariga:2018zuc,Ariga:2018pin} has been proposed to search for LLPs~\cite{Beacham:2019nyx, Alimena:2019zri} produced in the forward region of the LHC~\cite{Feng:2017vli, Kling:2018wct, Feng:2018noy, Ariga:2018uku, Berlin:2018jbm, Ariga:2019ufm, Jodlowski:2019ycu}, as well as to study interactions of high-energy neutrinos~\cite{Abreu:2019yak, Abreu:2020ddv, Kling:2020iar}. It utilizes the fact that light and high-energy particles produced in $pp$ collisions at the LHC and e.g. subsequent decays of light mesons $M$, will typically travel in the forward direction, as dictated by their estimated angular spread around the beam collision axis, $\theta_X\sim m_M/E_X\ll 1$. As a result, even a small detector placed along the beam collision axis can search for displaced decays of LLPs, provided that it is shielded to suppress the SM BG.
The FASER detector will operate during the LHC Run 3 in the former service tunnel TI12 located at a distance $L=480~\text{m}$ away from the ATLAS IP. Below, we also show the sensitivity reach for the same detector with continued operation during the High-Luminosity LHC phase (HL-LHC). We refer to this version of the experiment as FASER~HL. We also present the results obtained for a possible larger version of the detector to operate during HL-LHC, dubbed FASER~2. In particular, we assume a cylindrical decay volume with length $\Delta$, radius $R$, and an integrated luminosity $\mathcal{L}$ for each of the three aforementioned cases:
\begin{equation} \begin{aligned}
\textrm{\textbf{FASER:}}
&\ \Delta = 1.5~\text{m} ,
\ R=10~\text{cm},
\ \mathcal{L} = 150~\text{fb}^{-1},\nonumber\\
\textrm{\textbf{FASER~HL:}}
&\ \Delta = 1.5~\text{m},
\ R=10~\text{cm},
\ \mathcal{L} = 3~\text{ab}^{-1},\nonumber\\
\textrm{\textbf{FASER~2:}}
&\ \Delta = 5~\text{m},
\ R=1~\text{m}, \phantom{0c}
\ \mathcal{L} = 3~\text{ab}^{-1}.\nonumber
\end{aligned} \end{equation}
Our main focus in this study is to highlight possible complementarity between FASER searches for LLPs and kaon factories, with a particular emphasis on the recently observed KOTO anomaly. It is then useful to briefly summarize the main advantages of FASER that make it an ideal tool to study related BSM scenarios predicting a new unstable light species:
\begin{itemize}
\item FASER has the capability to study di-photon final states in LLP decays, on top of more often considered electron-positron pairs.
\item Forward-going kaons produced in $pp$ collisions at the ATLAS IP can travel about $140~\text{m}$, and hence decay with sizable probability, before they are absorbed. Additional production modes can further improve detection prospects.
\item FASER's sensitivity is optimal for typical LLP mass and lifetime proposed to explain the KOTO anomaly, cf. discussion in \cref{sec:KOTO}.
\end{itemize}

The signatures that we focus on consist of LLP $X$ decays into mainly $\gamma\gamma$ or $e^+e^-$ pairs that can be resolved in the detector. For sufficiently large energy, which we take to be $E_{X}\gtrsim 100~\text{GeV}$, the search at the FASER location can be considered BG-free~\cite{Ariga:2018zuc,Ariga:2018pin}. In the following, we assume a $100\%$ detection efficiency. Notably, the preliminary efficiency studies show that it typically has a minor impact on the sensitivity reach plots~\cite{Ariga:2018zuc}.

\paragraph{FASER capabilities to study the $\gamma\gamma$ final state} As far as the $\gamma\gamma$ final state is concerned, a very good separation efficiency between the photons is achieved thanks to a dedicated pre-shower detector placed in front of the FASER calorimeter, i.e. $2~\text{m}$ after the end of the FASER decay vessel. A preliminary relevant discussion can be found in Ref.~\cite{Feng:2018noy} in the context of FASER search for axion-like particles with the dominant di-photon coupling. As shown there, taking into account a finite separation resolution between the two $\gamma$s of order $\delta\sim 200$~nm, which is achievable in FASER, has a negligible impact on the sensitivity reach. Instead, the substantial effect was expected only if the resolution was worse than about $1$~mm.

Additionally, LLP decays into photons could also be proved in the single photon channel. Importantly, single high-energy BG photons with $E_\gamma\gtrsim 100~\text{GeV}$ at FASER location are typically associated with time-coincident muons activating veto layers. Hence, a large number of very energetic photon pairs produced in LLP decays even closer to the pre-shower detector, which could be misidentified as a single EM shower, would still be indicative of new physics.

\paragraph{Forward going kaons at the LHC} It is also important to stress that forward-going kaons produced at the ATLAS IP will not be immediately absorbed. Instead, they can travel a long distance through the LHC beam pipe, which allows them to decay with a sizable probability. In particular, for neutral kaons produced within the angular coverage of FASER, $\theta_K < \theta_{\textrm{FASER}}\approx 0.2~\text{mrad}$, the closest element of the LHC infrastructure that they hit is the TAN neutral particle absorber placed at a distance $140~\text{m}$ away from the ATLAS IP. We take this into account in our modeling along with the TAS absorber placed $20~\text{m}$ away from the IP that affects neutral kaons produced with larger $\theta_K$. The TAS also marks the end of the region where charged kaons are not deflected away by strong LHC magnets.

Below, we also illustrate how in specific models additional production modes of LLPs, for example in heavy meson decay, can further improve FASER sensitivity.

\paragraph{Lifetime regime} Last but not least, it is worth highlighting that in order for the LLP with $m_X\sim 100~\text{MeV}$ and energy $E_X\sim \text{TeV}$ to reach FASER, the required lifetime is of order $\tau_X\sim L/(\gamma_X\,c) \sim 0.1$~ns. This precisely corresponds to the sweet spot between the KOTO and NA62 searches discussed in \cref{sec:KOTO}. Therefore, as we will see below, it is not a surprise that FASER can effectively exclude many such explanations of the anomaly even with the first $10~\text{fb}^{-1}$ of data. Instead, in the case of discovery, FASER can confirm the anomaly with large statistics of related LLP decay events. This can reach up to $\mathcal{O}(10^4)$ events in FASER during LHC Run 3, while it could grow even larger for FASER~HL and FASER~2.

\section{\label{sec:results}FASER sensitivity reach}

As discussed in \cref{sec:KOTO}, if the anomalous neutral kaon decay events observed in the KOTO experiment are confirmed, this would require a BSM explanation. In the following, we illustrate how such models can be independently probed in FASER. To this end, we focus on several scenarios predicting LLPs.

\subsection{\label{sec:genericX} Generic LLP and the lifetime solution}

\begin{figure*}[tb]
\centering
\includegraphics[width=0.49\textwidth]{GenericX.pdf}
\hfill
\includegraphics[width=0.49\textwidth]{SpecificX.pdf}
\caption{Results for the \textsl{lifetime solution} to the KOTO anomaly (white regions in the plots) presented in the $(m_X,c\tau_X)$ plane, where $m_X$ and $\tau_X$ correspond to the LLP mass and lifetime, respectively. The current bounds are shown as gray-shaded regions.
\textbf{Left:} Reach for the generic X scenario in which a new light and unstable particle is produced specifically in the $K_L$ and $K^+$ decays, cf. discussion in \cref{sec:genericX}. From top to bottom, dark and light red solid lines correspond to the upper limit in $c\tau_X$ of the sensitivity reach of the FASER~HL and FASER experiments, respectively. The experiments can probe scenarios relevant for the entire white region in the plot below the sensitivity lines. Also shown, are the lines with a fixed number of events $N_{\textrm{ev}}=1000$ expected for the FASER~2 (orange dotted line) and SHiP (blue dashed line) experiments. The number of events in both cases grow for lower values of $c\tau_X$ in the white region of the plot. At the bottom of the plot, the dark red solid line marks the region of the parameter space which is excluded by the untagged $K_L \to \pi^0 X $ branching fraction.
\textbf{Right:} Event contours (from top to bottom: $10$ (solid), $10^2$ (dash-dotted), $10^3$ (dashed), $10^4$ (dotted), $10^5$ (solid) and $10^6$ (dash-dotted)) at FASER at the LHC Run 3 for the dark Higgs boson model with non-universal couplings to the SM leptons and quarks, as discussed in \cref{sec:2HDM}. The green-shaded band corresponds to the $2\sigma$ fit to the $(g-2)_\mu$ anomaly.
}
\label{fig:genericandspecificX}
\end{figure*}

The simplest BSM scenario proposed to explain the anomaly and to avoid stringent bounds on rare decays of charged kaons, employs a new unstable $X$ particle produced in kaon decays, $K_L\rightarrow \pi^0\,X$, and $K^+\rightarrow \pi^+\,X$. In particular, if the lifetime of $X$ and its mass are of order $\tau_X\sim 0.1$~ns and $m_X\sim 100~\text{MeV}$, respectively, then $X$ typically decays within the volume of the E949 and NA62 detectors and does not contribute to the measured value of $\mathcal{B}(K^+\rightarrow\pi^+\nu\bar{\nu})$. On the other hand, the probability for $X$ to escape the KOTO decay volume without decaying remains larger~\cite{Kitahara:2019lws}.

In the left panel of \cref{fig:genericandspecificX}, we show the FASER and FASER~HL sensitivity reach contours in a generic scenario, in which $X$ is produced only in kaon  decays. We assume that the real branching fraction of a two-body neutral kaon decay is of order $\mathcal{B}(K^0\rightarrow\pi^0 X)\sim 10^{-8}-10^{-6}$. Its precise value is chosen depending on $m_X$ and $\tau_X$, as well as the KOTO detector efficiency~\cite{Ahn:2018mvc}, as detailed in Ref.~\cite{Kitahara:2019lws,Liu:2020qgx}. This is done by requiring that the measured branching fraction fits the KOTO anomaly, cf. \cref{eq:K0Lneutanomaly}, after correcting for finite escape probability from the detector. In the plot, it is also assumed that the charged kaon decay branching fraction saturates the GN bound, cf. \cref{eq:GNbound}. The gray-shaded region is excluded by the constraints from the aforementioned searches for $K^+\rightarrow\pi^+\nu\bar{\nu}$ at E949~\cite{Artamonov:2009sz} and NA62~\cite{NA62bound}, the untagged $K_L \to \pi^0 X $ branching fraction~\cite{Tanabashi:2018oca}, beam-dump searches at CHARM~\cite{Bergsma:1985qz} and NuCal~\cite{Blumlein:1990ay,Blumlein:1991xh} experiments that we implement following Refs.~\cite{Liu:2020qgx,Egana-Ugrinovic:2019wzj}, and the measurement of the branching fraction for $K_L \rightarrow \pi^0 \gamma\gamma $ by KTEV~\cite{Abouzaid:2008xm}, assuming that $X$ decays dominantly into a $\gamma\gamma$ final state.\footnote{In particular, the NuCal bounds shown in \cref{fig:genericandspecificX} have been obtained following a ``conservative'' approach discussed in Ref.~\cite{Egana-Ugrinovic:2019wzj} and by taking into account $X$ production in decays of kaons (left panel) or all the relevant mesons (right panel). We note that the bounds from hadronic beam dumps are subject to uncertainties in modeling of the high-energy tail of the meson spectrum in the target.}

The anomalous KOTO events can be fitted within the remaining portion of the parameter space shown in \cref{fig:genericandspecificX}. Notably, the reconstructed transverse momentum of the pion can more easily resemble the observed such distribution by KOTO when $m_X\lesssim 180~\text{MeV}$~\cite{Kitahara:2019lws}, or even $m_X\lesssim m_\pi$~\cite{Liao:2020boe}.

In order to obtain the FASER sensitivity reach, we generate the meson spectra using \texttt{EPOS-LHC}~\cite{Pierog:2013ria} as implemented in the simulation package \texttt{CRMC}~\cite{CRMC} and subsequently decay the mesons using the branching fractions obtained as discussed above. As can be seen, already the FASER experiment operating during the LHC Run 3 can already cover the entire region of the parameter space corresponding to the anomaly. The relevant expected number of events can be as large as $\mathcal{O}(10^2)$ for FASER, and grow by an additional factor of $20$ and a few hundred for FASER~HL and FASER~2, respectively. Notably, thanks to the use of the dedicated pre-shower detector, FASER can probe this scenario even for leptophobic $X$ which decays dominantly into a di-photon pair, as assumed in the bounds shown in the left panel of \cref{fig:genericandspecificX}.

The lifetime solution to the KOTO anomaly can also be explored by other experiments dedicated to LLP searches, although this can be limited by the lack of di-photon signal detection capabilities. We illustrate this in the left panel of \cref{fig:genericandspecificX} for the proposed SHiP experiment~\cite{Alekhin:2015byh} to operate in a similar time frame to FASER~HL and FASER~2. The sensitivity of SHiP is analyzed by conservatively requiring that the kaons decay within one nuclear interaction length in the SHiP target made out of molybdenum, $\lambda = 15.25~\text{cm}$. We show the relevant contour line corresponding to a fixed number of expected events, $N_{\textrm{ev}}=1000$, similarly to the FASER~2 line also shown in the plot. The number of events grows larger for lower values of the LLP lifetime within the allowed region in the $(m_X,c\tau_X)$ plane. Interestingly, although kaon absorption in the target at distances greater than $\lambda$ could limit SHiP sensitivity, a large number of protons on target (POT), $N_{\textrm{POT}} = 2\times 10^{20}$, as well as the size of the decay volume, $\Delta \simeq 50~\text{m}$, sufficiently compensate for this effect. The number of expected events in SHiP could be further enhanced once more detailed modeling of the kaon propagation in the target is performed.

\subsection{\label{sec:2HDM}Specific example employing 2HDM}

\begin{table*}
\centering
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|c||c|c||c|c|c|c|c|}
\hline
\hline
Experiment & \multicolumn{2}{c||}{Benchmark} & \multicolumn{4}{c|}{$X$ prod. in meson decays and decay in FASER} & Proton brem \\
\cline{2-8}
\cline{2-8}
&$m_X$ [MeV] & $c\tau_X [\textrm{cm}]$ & $K_L\to\pi^0 X$ & $K^+\to\pi^+ X$ & $\eta,\eta^\prime,K_S\to\pi^0 X$ & $b\to s X$ & $pp\to X+\ldots$\\
\hline
\hline
 FASER &
 \multirow{2}{*}{50} & \multirow{2}{*}{1.5} &
 150& 50& 650& 2000&
 135\\
\cline{4-8}
 FASER~2 &
  &  &
 80k & 30k & 250k & 11M &
 16k\\
\hline
\hline
\end{tabular}
\caption{The number of $X$ decay events in FASER (FASER~2) and the integrated luminosity of $\mathcal{L} = 150~\text{fb}^{-1}$ ($\mathcal{L} = 3~\text{ab}^{-1}$) for various production modes of the LLP. The number of events for FASER~HL are $20$ times larger than for FASER during LHC Run 3. The results correspond to the model discussed in \cref{sec:2HDM} and to the benchmark values of the model parameters given in the text and indicated in the table. The scenario can explain both the KOTO and $(g-2)_\mu$ anomalies.
\label{tab:productionmodes}}
\end{table*}

Once a specific model of the LLP that corresponds to the aforementioned lifetime solution is considered, typically more production and decay modes appear for $X$ that should be taken into account. This has an impact on present constraints on such $X$, but can also significantly increase the expected number of events in the future detectors. We illustrate this below for FASER employing the model with a leptophilic dark scalar discussed in Ref.~\cite{Liu:2020qgx}. Interestingly, such a scalar with $40~\text{MeV}\lesssim m_X\lesssim 70~\text{MeV}$ and $\tau_X\sim 0.1$~ns could simultaneously explain the recently observed KOTO events and the measurement of the anomalous magnetic moment of the muon, $(g-2)_\mu$~\cite{Bennett:2006fi}.

The model of our interest can effectively be described by three parameters corresponding to the LLP mass, as well as its coupling constants to leptons, $\epsilon_\ell$, and quarks, $\epsilon_q$. The relevant Lagrangian reads
\begin{equation} \begin{aligned}
\mathcal{L}  \supset & - m_X^2 X^2
+ \sum \epsilon_q \frac{m_q}{v} \, X \bar{q} q
+ \sum \epsilon_\ell \frac{m_\ell}{v} \, X \bar{\ell} \ell \\
&  + \epsilon_W \frac{m_W^2}{v}\, X W^+_\mu W^{\mu -},
\label{eq:L2HDM}
\end{aligned} \end{equation}
where the SM Higgs boson vev is equal to $v\simeq 246~\text{GeV}$ and we set $\epsilon_W=\epsilon_q$. The couplings $\epsilon_\ell$ and $\epsilon_q$ can arise as the mixing angles in the scalar sector of the type-X two Higgs doublet model (2HDM). In particular, the model allows one to disjointly treat interactions of $X$ with the SM leptons and quarks.

Assuming a hierarchy between the couplings, $\epsilon_\ell\gg\epsilon_q$, allows one to reduce the lifetime of $X$, which is governed by the leptonic coupling $\epsilon_\ell$. Thanks to this, stringent bounds from hadronic beam-dump searches can be avoided and the kaon decay branching fraction into $X$ can be kept low, as it depends on $\epsilon_q\sim 10^{-3}-10^{-2}$.  In particular, for $\epsilon_\ell\sim \mathcal{O}(1)$ the most important bounds on long-lived $X$ come from electron beam-dump experiments E137~\cite{Bjorken:1988as} and Orsay~\cite{Davier:1989wz}. However, they do not exclude scenarios with $\tau_X\sim 0.1$~ns.

The dominant decay mode for $X$ in the considered mass regime is into the $e^+e^-$ final state, although a loop-induced decay into the $\gamma\gamma$ pair can take values of up to $\mathcal{O}(10\%)$ of the decay branching fraction of $X$. The additional bounds from the KTEV search for $K_L\rightarrow\pi^0\,e^+e^-$ decays exclude a part of the region of the parameter space of the model with $m_X\gtrsim 140~\text{MeV}$~\cite{AlaviHarati:2003mr}, although they are not relevant for lighter $X$. We employ the decay width of $X$ following Ref.~\cite{Liu:2020qgx}.

A relative smallness of $\epsilon_q$ is consistent with only a few events currently observed in KOTO. However, as far as FASER is concerned, additional production modes of $X$ can become more important than rare kaon decays, therefore increasing the expected number of signal events. In particular, when obtaining the FASER sensitivity reach, we take into account the following production channels:
\begin{description}
\item [Meson decays] In our modeling, we include rare two-body decays of charged and neutral long-lived kaons that propagate in the LHC beam pipe, as discussed in \cref{sec:FASER}. We also analyze the impact of short-lived kaon decays, $K_S \to \pi^0 X$. Notably, the relevant BSM branching fraction of $K_S$ is suppressed with respect to $K_L$ by an additional factor of order $\mathcal{O}(10^{-4})$ due to a larger value of the total decay width of $K_S$. On the other hand, short-lived kaons can more easily decay prior to hitting any element of the LHC forward infrastructure.

On top of kaons, a number of other mesons can be abundantly produced in the forward direction of the LHC. In particular, we study prompt rare two-body decays $\eta\to\pi^0 X$ and $\eta^\prime\to\eta X$, as well as inclusive decays of $B$ mesons into final states containing strange hadronic states, $b\to s X$, where the light meson spectra were obtained using \texttt{EPOS-LHC}, while the $b$-quark spectrum was obtained using \texttt{Pythia~8}~\cite{Sjostrand:2006za, Sjostrand:2014zea}. We implement the kaon and $B$-meson branching fractions following Ref.~\cite{Feng:2017vli}, while for $\eta$ mesons we employ the results from Ref.~\cite{Egana-Ugrinovic:2019wzj}. We give the relevant expressions for the branching fractions in \cref{app:BRs}.

The dominant contribution to the $K$ and $B$ decay widths into $X$ corresponds to a top-$W$ loop. Unlike for kaon decays, the relevant branching fractions of $B$ mesons do not suffer from a strong CKM suppression, $V^2_{tb}\gg V_{td}^2$, which makes this channel the dominant production mode despite the suppressed $b$-quark production rate and the broader $p_T$ spectrum. Instead, the decays of $\eta$ mesons are suppressed for small Yukawa-like couplings of $X$ to the first-generation quarks and for a loop-induced coupling to gluons. Similarly to $K_S$, however, prompt decays of $\eta$ mesons make them a non-negligible production mode of $X$.

\item[Scalar bremsstrahlung] Another production mode of light and high-energy scalar BSM particles produced in the forward direction of the LHC is through their bremsstrahlung in proton-proton collisions. We study this process following the discussion in Ref.~\cite{Boiarska:2019jym} and find that the relevant contribution to the event rate in FASER is typically of order a few $\%$ or smaller.

\item [Lepton-induced production] We have also analyzed a number of production modes that depend on the lepton coupling, $\epsilon_\ell\sim \mathcal{O}(1)$. These include loop-induced Primakoff production from high-energy photons hitting the TAN, cf. Ref.~\cite{Feng:2018noy}, as well as various secondary processes employing electrons and muons pair-produced in the absorber material or traversing the rock shielding of FASER. However, since the difference in $\epsilon_\ell$ and $\epsilon_q$ is typically not sufficient to compensate for a large Yukawa suppression of the corresponding coupling constants, as well as due to additional suppression factors relevant for these production modes, we have found that they play a subdominant role for high-energy scalars of our interest.
\end{description}

In the right panel of \cref{fig:genericandspecificX}, we show the contours with the different number of events in FASER, $N=10,10^2\ldots$, that correspond to the total integrated luminosity of $\mathcal{L}=150~\text{fb}^{-1}$ for the LHC Run 3. In the white region of the plot, the KOTO anomaly can be explained without violating current bounds shown as the gray-shaded region, following Ref.~\cite{Liu:2020qgx}. In addition, we have added the constraint from the NuCal experiment following Ref.~\cite{Egana-Ugrinovic:2019wzj}. The green-shaded band in the plot indicates the region in the parameter space of the model in which the $(g-2)_\mu$ anomaly can be resolved. As can be seen, FASER can detect up to $\mathcal{O}(10^4)$ LLP decays in the region of the parameter space corresponding to both the anomalies. In this case, even a few weeks of the operation of the experiment would be enough to test such scenarios. Once more data are accumulated, the entire allowed region of the parameter space will be covered with at least a few hundred expected events.

The breakdown of the number of the expected events in FASER and FASER~2 for different production modes of $X$ is shown in \cref{tab:productionmodes} for the benchmark scenario with $m_X=50~\text{MeV}$, $\epsilon_q = 0.016$ and $\epsilon_\ell = 1.22$. The values of the parameters of the model have been chosen so that both the KOTO and $(g-2)_\mu$ anomalies can be fitted to their central values~\cite{Liu:2020qgx}. The dominant production mode in FASER, in this case, is due to rare $B$ meson decays, although decays of lighter mesons give a contribution of order $30\%$ in total. This is dominated by rare decays of $\eta$ mesons. Instead, the kaon decays give only a few $\%$ contribution.

\subsection{\label{sec:longlifeX}A very long-lived dark scalar}

If a new LLP is produced in rare kaon decays with the mass close to the pion mass, $m_X\sim m_\pi$, and the relevant branching fraction fitting the anomaly, cf. \cref{eq:K0Lneutanomaly}, the stringent NA62 and E949 bounds on $K^+$ decays can be avoided even for a very long-lived or effectively stable $X$. This is due to increased BG in the respective searches for $K^+$ decays (see the discussion in Ref.~\cite{Kitahara:2019lws} and references therein). In addition, such scenarios can also escape constraints from beam-dump searches provided that $\tau_X$ is large enough and the LLPs typically overshoot the detector.

An interesting example of such a scenario that can fit the KOTO anomaly is a dark Higgs boson $X$ with a universal mixing angle with the SM species, cf. \cref{eq:L2HDM} with $\epsilon_\ell = \epsilon_q = \epsilon_W$ and Ref.~\cite{Egana-Ugrinovic:2019wzj} for further discussion. In this scenario, the lifetime of the LLP is typically large, $c\tau_X\sim 100~\text{km}$. As a result, probing this model goes beyond the capabilities of FASER and FASER~2 experiments, which are designed to focus on more short-lived BSM species.

On the other hand, in the long-lifetime regime relevant here, the sensitivity reach in searches for displaced decays of $X$ can be improved by increasing the size of the detector and its angular coverage. We illustrate this in the left panel of \cref{fig:darkHiggsandX1X2}, where we show the expected reach of proposed future Codex-b~\cite{Gligorov:2017nwh}, MATHUSLA~\cite{Curtin:2018mvb} and SHiP~\cite{Alekhin:2015byh} experiments in their searches for dark Higgs boson decays into electron-positron pairs, $X\to e^+e^-$, following Ref.~\cite{Beacham:2019nyx}.

For comparison, we also present the expected sensitivity of the enlarged FASER~2 detector assuming the total integrated luminosity of $\mathcal{L} = 3~\text{ab}^{-1}$ relevant for HL-LHC. As can be seen, probing the region in the parameter space of the model that corresponds to the KOTO anomaly would require a quite substantial increase in the detector radius, $R\gtrsim 5~\text{m}$, or its length, $\Delta\gtrsim 50~\text{m}$, with respect to the design discussed in \cref{sec:FASER}. On the other hand, for lower values of $\tau_X$, even a much smaller and properly placed detector can have very good detection prospects, as we show for other scenarios discussed in this study.

\begin{figure*}[tb]
\centering
\includegraphics[width=0.49\textwidth]{DarkHiggs.pdf}
\hfill
\includegraphics[width=0.49\textwidth]{Pospelov.pdf}
\caption{The results for the models discussed in \cref{sec:longlifeX} and \cref{sec:X1X2}. In both panels, the $2\sigma$ regions of the parameter space in which the KOTO anomaly can be explained are shown as the green-shaded bands. \textbf{Left:} Reach for the dark Higgs boson model discussed in \cref{sec:longlifeX}. The sensitivity contour for FASER~2 is shown with the orange solid line, while purple and yellow solid lines correspond to larger versions on FASER~2 with the radius $R=5~\text{m}$ or length $\Delta=50~\text{m}$, respectively. The current bounds are shown as gray-shaded regions following Ref.~\cite{Egana-Ugrinovic:2019wzj}. The blue dashed lines correspond to the expected sensitivity of other proposed future experiments (from top to bottom): CODEX-b, SHiP and MATHUSLA.
\textbf{Right:} Event contours (from outside to inside: 3 (solid red line), 10 (dash-dotted), 100 (dashed), 1000 (dotted)) at FASER at LHC Run 3 for the model with pure dark sector decays of kaons, $K_L\to X_1X_2$, discussed in \cref{sec:X1X2}. The current bounds are shown as gray-shaded regions following Ref.~\cite{Hostert:2020gou}.
}
\label{fig:darkHiggsandX1X2}
\end{figure*}

\subsection{\label{sec:X1X2} $K_L$ decays into dark sector particles}

A different approach to the BSM explanation of the KOTO anomaly has been proposed in Refs~\cite{Gori:2020xvq,Hostert:2020gou}; this approach employs two-body decays of neutral kaons into dark species, $K_L\to \psi_1 \psi_2$. The relevant three-body decay processes of charged kaons, e.g. $K^+\to \pi^+ \psi_1 \psi_2$, can then be kinematically suppressed or even forbidden depending on the masses of the $\psi_i$s. The anomaly can be fitted if at least one of the dark species is unstable and can mimic the di-photon decay signature of neutral pions inside the KOTO detector.

In the following, we focus on the model in which a non-diagonal coupling of kaons to two dark fermionic states $\psi_i$ arises due to a vector portal and a new gauge field $X_\mu$ that mixes with the SM $Z$ boson~\cite{Babu:1997st, Davoudiasl:2012ag, Dror:2018wfl}
\begin{equation}
\mathcal{L}\supset \frac{1}{2}m_X^2 X_\mu X^\mu - \epsilon_Zm_Z^2X_\mu Z^\mu,
\label{eq:LX1X2}
\end{equation}
where $\epsilon_Z\sim 10^{-3}-10^{-2}$ denotes the relevant mixing parameter. After the electroweak symmetry breaking, a new $Z^\prime$ gauge boson acquires a mass $m^2_{Z^\prime} = m_X^2 - \epsilon_Z^2 m_Z^2$. The interactions of $Z^\prime$ with the dark fermions can then be described by
\begin{equation}
\mathcal{L}\supset g_XZ^\prime_\mu\,(c_V\bar{\psi}_2\gamma^\mu\psi_1 + c_A\bar{\psi}_2\gamma^5\gamma^\mu\psi_1) + h.c.
\end{equation}
In particular, an effective operator which couples $\psi_i$s to $s$ and $d$ quarks is induced at a loop level with $W$ boson exchange~\cite{Dror:2018wfl}. In the following, we set  $g_X=1$, $c_V=0$, $c_A=1$, $m_{Z^\prime} = 10~\text{GeV}$, and $m_{2} = 11\,m_{1}$, where $m_i\equiv m_{\psi_i}$.

The $\psi_2$ decay width and branching fractions for $\psi_1\psi_2$-pair production in various meson decays can be found in Ref.~\cite{Hostert:2020gou}. In particular, as discussed therein, the assumption about a relatively large mass splitting between the dark species, $m_{1}\ll m_{2}-m_{\pi^0}$, allows one to obtain about an $80\%$ branching fraction of two-body decays $\psi_2\to\psi_1\pi^0$. The subsequent prompt decays of high-energy neutral pions with $E_{\pi^0}\gtrsim 100~\text{GeV}$ generate a di-photon signature inside the FASER detector.

In the right panel of \cref{fig:darkHiggsandX1X2}, we show the contours with the expected number of events in FASER in the parameter space spanned by $m_{2}$ and $G_X/G_F$ parameters, where $G_F$ is the Fermi coupling constant and $G_X = (\epsilon_Z g g_X)/(2\sqrt{2}c_Wm_{Z^\prime}^2)$ with $g$ corresponding to the weak coupling constant. The region of the parameter space in which a $2\sigma$ fit can be obtained to the KOTO anomaly is indicated by a green-shaded band in the plot. The current bounds on the model from the CHARM experiment, the previous KOTO search for $K_L\to \pi^0\nu\bar{\nu}$, cf. \cref{eq:K0Lneutbound}, atomic parity violation (APV)~\cite{Dror:2018wfl}, and the Belle-II search for $B\to K\nu\bar{\nu}$~\cite{Grygier:2017tzo} are shown as the gray-shaded region following Ref.~\cite{Hostert:2020gou}.\footnote{As discussed in Ref.~\cite{Hostert:2020gou}, a more detailed simulation of the CHARM and NuCal experiments is needed to refine the lower part of the gray-shaded region in \cref{fig:darkHiggsandX1X2}. This, however, is expected to have a small impact on the upper part of the plot corresponding to the KOTO anomaly.}

The region of the parameter space of the model in which the KOTO anomaly can be explained by rare decays $K_L\to \psi_1\psi_2$ corresponds to $350~\text{MeV}\lesssim m_{2}\lesssim 450~\text{MeV}$. Here, one requires that $\psi_2$ decays within the KOTO detector volume so that it can mimic neutral pion decays. However, such decays that are too fast lead to larger production rates that are already excluded by previous KOTO studies. As a result, the anomalous events can be best fitted for $\tau_{\psi_2}\sim (0.01-0.1)$~ns and such scenarios can be well tested in FASER. In particular, as shown in \cref{fig:darkHiggsandX1X2}, one expects $\mathcal{O}(100)$ of high-energy visible $\psi_2$ decay events in the detector for LHC Run 3.

\section{\label{sec:ALPs} Beyond the anomaly -- ALP coupling to SU(2)$\mathbf{_W}$ gauge bosons}

The complementarity between FASER and kaon factories extends beyond possible explanations to the currently observed anomalous events in the KOTO experiment. In particular, more general BSM scenarios leading to two-body kaon decays, $K\to\pi^0 X$ and $K^+\to\pi^+ X$, can be constrained in searches for detector-stable $X$ acting as neutrino impostors, as well as in studies focused on displaced decays $X\to\gamma\gamma$ leading to $4\gamma$ or $\pi^+\!+2\gamma$ signatures observed in the detector. Some examples of such studies have recently been discussed in Ref.~\cite{Gori:2020xvq} for models with axion-like particles (ALPs), $X\equiv a$, which couple dominantly to gluons or $SU(2)_W$ gauge bosons. The sensitivity reach of FASER in the former scenario has been studied in Ref.~\cite{Ariga:2018uku}.

In the following, we will focus on the latter model in which the ALP couples to the SM field strength tensor $W_{\mu\nu}^a$ of the $SU(2)$ group
\begin{equation}
\mathcal{L}\supset -\frac{1}{2}m_a^2 a^2 - \frac{g_{aWW}}{4}a W_{\mu\nu}^a \widetilde{W}_{\mu\nu}^a.
\label{eq:LALP}
\end{equation}
The coupling to $W$ bosons in \cref{eq:LALP} gives rise to both kaon decays, $s\!\to\! d\, a$, and $B$-meson decays, $b\!\to\! s \,a$, via the usual loop diagrams. We implement the relevant branching fractions following Ref.~\cite{Izaguirre:2016dfi,Aloni:2018vki,Gori:2020xvq} as discussed in \cref{app:BRs}.

After the EWSB, additional couplings of $a$ to $\gamma\gamma$, $Z\gamma$, and $ZZ$ arise with the relative strength dictated by the weak mixing angle. As a result, the ALP can also be produced through the Primakoff process, $\gamma N\to a\gamma$, employing the di-photon coupling $g_{a\gamma\gamma} = g_{aWW}s_W^2$, cf. Ref.~\cite{Feng:2018noy} for relevant discussion for FASER, providing a subleading contribution to the ALP event rate. In contrast, $a$ production in rare $Z$ decays is more isotropic and can be neglected. For the low-mass ALPs of our interest, the dominant decay mode is into two photons, $a \to \gamma\gamma $~\cite{Gavela:2019wzg}.

We present the FASER, FASER~HL, and FASER~2 sensitivity reach lines for this model in \cref{fig:ALPs}. Additional dashed contours in the plot correspond to the expected future sensitivity in searches in KOTO and NA62 for the aforementioned two-body BSM kaon decays, $K^{(+)}\to \pi^{0(+)}a$. Depending on whether the ALP escapes the detector without decaying or decays inside to two photons, the searches focus on $2\gamma$ or $4\gamma$ signatures in KOTO, and $\pi^++0\gamma$ or $\pi^++2\gamma$ signal in NA62. We show these searches following Ref.~\cite{Gori:2020xvq}. The current bounds on the model are taken from Refs.~\cite{Dolan:2017osp,Gori:2020xvq}. The most important constraints in the region of the parameter space of our interest come from the beam-dump experiment E137 and past searches in kaon factories.

As can be seen, FASER will cover the entire region of the parameter space corresponding to KOTO and NA62 searches for detector-stable $a$. It will also independently probe scenarios relevant for unstable ALP decaying inside these detectors. In addition, FASER reach extends towards lower values of the $g_{aWW}$ coupling and larger mass of the ALP. The larger values of the coupling constant will, instead, be probed by e.g. Belle-II searches for the $3\gamma$ signature from the ALP production and subsequent prompt decay inside the detector, $e^+e^-\to\gamma+(a\to2\gamma)$. We show the relevant reach for $50~\text{fb}^{-1}$ in \cref{fig:ALPs} following Ref.~\cite{Dolan:2017osp}. Similarly, future ATLAS and CMS searches for $Z\to\gamma+(a\to2\gamma)$ signature can also cover an important part of the parameter space of the model corresponding to $g_{aWW}\gtrsim 10^{-5}\,\textrm{GeV}^{-1}$, as shown in the plot based on Ref.~\cite{Bauer:2017ris}.

Both FASER~HL and FASER~2 experiments will be able to further improve the detection prospects, as well as constrain scenarios with $m_a$ up to $1~\text{GeV}$ or even close to the limit $m_a\lesssim m_B-m_K$, respectively. In a similar time frame, additional bounds on the model can come from e.g. the SHiP experiment (not shown in the plot).

\begin{figure}[tb]
\centering
\includegraphics[width=0.49\textwidth]{ALP-W.pdf}
\caption{Sensitivity reach of the FASER (red solid line), FASER~HL (brown) and FASER~2 (orange) experiments in the model with axion-like particles coupled to $SU(2)_W$ discussed in \cref{sec:ALPs}. The currently excluded gray-shaded region and the dashed blue lines corresponding to the future sensitivity reach in the KOTO and NA62 searches for promptly decaying and detector-stable ALPs are taken from Ref.~\cite{Dolan:2017osp,Gori:2020xvq}. The expected reach of the Belle-II search for the $3\gamma$ signature is shown following Ref.~\cite{Dolan:2017osp}. The future ATLAS/CMS reach in the $Z\to 3\gamma$ signature from production and prompt decays of the ALP to a di-photon pair is shown following Ref.~\cite{Bauer:2017ris}.
}
\label{fig:ALPs}
\end{figure}

\section{\label{sec:conclusion} Conclusion and Outlook}

A long-awaited definite experimental discovery of new physics effects should first manifest itself as a clear deviation from the SM predictions and excess over expected BG. In this study, we have focused on recently observed and intriguing anomalous neutral kaon decay events in the KOTO experiment. The current KOTO observation corresponds to a striking excess over the expected SM branching fraction of the decay process $K_L\to\pi^0\nu\bar{\nu}$, while it awaits further dedicated analyses to be performed by the KOTO collaboration and other experiments before it can be fully confirmed. In particular, we have analyzed the prospects of independently probing the relevant BSM scenarios proposed to explain the anomaly in the upcoming FASER experiment at the LHC~\cite{Feng:2017uoz,Ariga:2018zuc,Ariga:2018pin}.

We have shown that FASER could start probing such models, which predict light and unstable new particles, immediately after the beginning of its operation. The $pp$ collisions at the LHC produce large numbers of kaons that typically hit elements of the infrastructure before decaying unless they go down the beam pipe i.e. towards FASER. This can lead to a large flux of forward-going LLPs produced in rare kaon decays and in other production modes. The FASER baseline compatible with the predicted lifetime $c\tau_X$ of such LLPs, and its excellent di-photon detection capabilities, will then result in up to $\mathcal{O}(10^4)$ high-energy LLP decay events observed during LHC Run 3.

Our results are presented for selected distinct BSM scenarios corresponding to the \textsl{lifetime solution} to the KOTO anomaly, as well as for some other models that either predict very long-lived LLPs with a mass close to the pion mass or employ purely dark sector decays of neutral kaons. Interestingly, even in the less promising case of large $\tau_X$, FASER can also be indirectly sensitive to KOTO-related models that predict a larger set of LLPs and simultaneously address a number of outstanding experimental and theoretical issues. An example of such a scenario with below $\text{GeV}$-scale heavy neutral leptons within the reach of FASER 2~\cite{Kling:2018wct,Ariga:2018uku} has recently been discussed in Ref.~\cite{Cline:2020mdt}. Although in this study we focus on collider searches for new physics, it is important to note that the KOTO anomaly can also have profound cosmological consequences that would have to be thoroughly investigated, especially if the observed excess persists, cf. Ref.~\cite{Altmannshofer:2020pjb} for a recent discussion.

Last but not least, FASER complementarity to BSM searches in kaon factories extends beyond the current anomaly. We have illustrated this for a simplified model with ALPs coupled to $SU(2)_W$ gauge bosons that dominantly decay to $\gamma\gamma$ pairs. FASER will probe regions of the parameter space relevant for several distinct searches for LLPs decaying either inside or outside the KOTO and NA62 detectors. It will, therefore, provide an independent probe of such scenarios. In a longer term, improved bounds can be obtained by FASER~HL or FASER~2 detectors operating during the HL-LHC phase.

The upcoming run of the LHC will push the limits of intensity-frontier exploration of new physics to the next level. The FASER experiment will play a vital role in this endeavor toward potentially groundbreaking discoveries.

\acknowledgments
We would like to thank Jonathan Feng, Sam Homiller and Ben Lillard for useful discussions and comments on the manuscript. We are also grateful to the authors and maintainers of many open-source software packages, including
\texttt{CRMC}~\cite{CRMC},
\texttt{EPOS}~\cite{Pierog:2013ria}
\texttt{Jupyter} notebooks~\cite{soton403913},
\texttt{Matplotlib}~\cite{Hunter:2007},
\texttt{NumPy}~\cite{numpy},
\texttt{pylhe}~\cite{lukas_2018_1217032},
\texttt{Pythia~8}~\cite{Sjostrand:2014zea},
and \texttt{scikit-hep}~\cite{Rodrigues:2019nct}.
FK is supported by the Department of Energy under Grant DE-AC02-76SF00515. ST is supported by the Lancaster-Manchester-Sheffield Consortium for Fundamental Physics under STFC grant ST/P000800/1. ST is partially supported 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).

redpajama_set_name: RedPajamaArXiv


text:
\section{Late infall}

There are compelling reasons to believe that the dominant component of
the dark matter of the Universe is non-baryonic collisionless particles,
such as axions, weakly interacting massive particles (WIMPs) and massive
neutrinos. The word ``collisionless'' signifies that the particles are
so weakly interacting that they have moved purely under the influence of
gravity since their decoupling at a very early time (of order $10^{-4}$
sec for axions, of order 1 sec for neutrinos and WIMPs).  In the limit
where the primordial velocity dispersion of the particles is neglected,
they all lie on the same 3-dimensional (D) `sheet' in 6-D phase-space.
Their phase-space evolution must obey Liouville's theorem.  This implies
that the sheet cannot tear and hence that it satisfies certain topological
constraints.

Because their phase-space sheet cannot tear, collisionless dark matter
(CDM) particles must be present everywhere in space, including specifically
intergalactic space.  The space density may be reduced by stretching of the
sheet, but it cannot vanish.  Moreover, the average space density is
recovered as soon as the average is taken over distances larger than
the distance CDM may have moved locally away from perfect Hubble flow.
In a region that is sparsely populated with galaxies, this distance
is much smaller than the distance between galaxies.  The implication is that
isolated galaxies are surrounded by unseen CDM and hence, because of gravity,
CDM keeps falling continuously onto such galaxies from all directions.
If the galaxy joins other galaxies to form a cluster, infall onto the galaxy
gets shut off because of lack of material, but infall onto the cluster
continues assuming that the cluster is itself isolated.  In an open
universe $(\Omega < 1)$, the infall process eventually turns off because
the universe becomes very dilute.  However, even if our universe is open,
we are far from having reached the turn-off time.

\section{Velocity peaks}

Consider the infall of CDM onto an isolated galaxy.  Let us at first
neglect the velocity dispersion of the infalling particles.  In practice
it is sufficient that their velocity dispersion is much smaller than
the rotation velocity of the galaxy.  Consider the time evolution of all
CDM particles that are about to fall onto the galaxy for the first time
in their history at time $t$.  For the sake of definiteness, we may consider
all particles which have zero radial velocity $(\dot{r} =0)$ for the first
time then.  Such particles are said to be at their `first turnaround';
they were receding from the galaxy as part of the general Hubble flow
before $t$ and will be falling onto the galaxy just after $t$.  They form
a closed surface, enclosing the galaxy, called the turnaround 'sphere' at
time $t$.  The present turnaround sphere of the Milky Way galaxy has a
radius of order 2 Mpc.  The turnaround sphere at time $t$ falls through
the central parts of the galaxy at a time of order $2t$.  Particles falling
through the galactic disk (assuming the galaxy is a spiral) get scattered
through an angle $\Delta \theta \sim 10^{-3}$ by the gravitational fields
of various inhomogeneities such as molecular clouds, globular clusters,
and stars \cite{is}. However, most particles carry too much angular momentum
to reach the luminous parts of the galaxy and are scattered much less.
Because the scattering is small, the particles on the turnaround sphere
at time $t$, after falling through the galaxy, form a new sphere which
reaches its maximum radius $R'$ at some time $t'$.  The radius $R'$ at
the second turnaround is smaller than the radius $R$ at the first
turnaround because the galaxy has grown by infall in the meantime.  The
sphere continues oscillating in this way although it gets progressively
fuzzier because of scattering off inhomogeneities in the galaxy.

\begin{figure}[top]
\psfig{figure=infall.eps,height=4.5in}
\vskip 2cm
\caption{\small{Infall of a turnaround sphere.  The sphere has net
angular momentum about the vertical axis.  It crosses itself between
frames b) and c).  After frame e) the sphere has completed the process
of turning itself inside out.  The cusps in frames d) and e) are at
the intersection of a ring caustic with the plane of the figure.}}
\label{fig1}
\end{figure}

Fig.~1 describes the fall of one turnaround sphere through the galaxy
as a succession of time frames.  No particular symmetry is assumed.  It
is assumed, for the sake of definiteness, that the particles on the
turnaround sphere carry net angular momentum about the vertical axis.
The particles which are near the top in frame a) carry little angular
momentum.  They fall through the center of the galaxy and end up near the
bottom of the turnaround sphere in frame f).  Similarly the particles
which are near the bottom in frame a) end up near the top in frame f).
The particles which are near the equator carry the most angular momentum.
They form a ring whose size decreases to some minimum value and then
increases again.  In the process of falling through the galaxy the
turnaround sphere turns itself inside out.  Note also that the sphere
crosses at least twice each point which is inside the sphere both at
the initial time of frame a) and at the final time of frame f).

Fig.~1 shows the motion of just one sphere in a continous flow of such
spheres.  Moreover, there are many flows in and out of the galaxy going
on at the same time.  To each flow in and out of the galaxy is associated
a pair of peaks in the velocity distribution of CDM particles at every
physical point in the galactic halo \cite{is}.  One peak is due to
particles falling onto the galaxy for the first time, one peak is due to
particles falling out of the galaxy for the first time, one peak is due
to particles falling onto the galaxy for the second time, and so on.  In
particular this is true on Earth.  Igor Tkachev, Yun Wang and I \cite{twi}
obtained estimates of the velocity magnitudes and the local density
fractions associated with these peaks using the self-similar infall model
of galactic halo formation.  We generalized the existing version of the
model \cite{ssi} to take account of the angular momentum of the dark matter
particles.  We find that the first twenty peaks contribute each between
one to four percent of the local density.  This implies a large
non-thermal component to the galactic CDM distribution.

\section{Caustic rings}

There is a caustic ring associated with each flow in and out of the
galaxy \cite{ring}.  A caustic is a place in physical space where the
density is large because the phase-space sheet `folds back' there.  At
the caustic, the space density diverges in the limit of zero velocity
dispersion.  In Fig.~1, the caustic ring intersects the plane of the
figure at the location of the cusps in frames d) and e).  One can prove
mathematically \cite{tbp} that the density diverges at these points in
the limit of zero velocity dispersion of the infalling particles.  In
reality the infalling particles have some velocity dispersion.  However
only when this velocity dispersion is as large as 30 km/s do the caustic
rings in our galaxy get washed out.  It is generally believed that the
velocity dispersion of infalling CDM particles is at most 10 km/s.

For an arbitrary angular momentum distribution on the turnaround sphere,
the ring is a closed loop of arbitrary shape.  However, if the angular
momentum distribution is dominated by a smooth component that carries
\emph{net} angular momentum, the ring resembles a circle.  If there is
no angular momentum at all, the ring reduces to a point at the galactic
center.

The caustic ring is located in physical space near where the particles
with the largest amount of angular momentum are at their distance
of closest approach to the galactic center.  In the particular case where
the turnaround sphere is initially rotating with a definite angular
velocity as if it were a rigid body, the density distribution near the
caustic is:
\begin{equation}
d(a;\rho,z) \simeq {dM \over d\Omega dt}~ {2 \over v}~
{1 \over \sqrt{(r^2-a^2)^2 + 4 a^2 z^2}}
\end{equation}
where $(\rho,z,\theta)$ are cylindrical coordinates, $a$ is the caustic
ring radius assumed to be much smaller than $R$, $r= \sqrt{\rho^2 + z^2},
~{dM \over d\Omega dt}$ is the rate at which mass falls in per unit time
and unit solid angle, and $v$ is the velocity magnitude of the particles
at the caustic.  Near the caustic, the density diverges as the inverse
distance to the ring:
\begin{equation}
d(a;\rho,z) \simeq {dM\over d\Omega dt}~ {1\over v~ a~ \sigma}
\label{(3)}
\end{equation}
with $\sigma = \sqrt{(\rho -a)^2 + z^2}$.

The self-similar infall model mentioned earlier \cite{twi} predicts the
values of successive radii to be \cite{ring}:
\begin{equation}
\{a_n: n=1,2,3...\} \simeq (39,19.5,13,10,8 ...) {\rm kpc}
\left({j_{max} \over 0.25}\right) \left({0.7 \over h}\right)
\left({v_{rot} \over 220 {\rm km/s}}\right)
\end{equation}
for $\epsilon=0.3$.  $\epsilon$ is a parameter which is predicted by
theories of large scale structure formation to be in the range 0.2 to
0.35, $j_{max}$ is the maximum value of the angular momentum of the
particles on the turnaround sphere in the dimensionless units
defined in Ref.[2], $h$ is the Hubble rate in units of
100 km/sec.Mpc and $v_{rot}$ is the rotation velocity of the galaxy.
For $\epsilon = 0.2,~ a_1 \simeq 36 {\rm kpc}~\left({j_{max}\over 0.25}
\right) \left({0.7\over h}\right) \left({v_{rot}\over 220
{\rm km/s}}\right)$, but the ratios $a_n/a_1$ are almost the same as
in the $\epsilon = 0.3$ case.  Thus in the range of interest $0.2 \leq
\epsilon \leq 0.35$, the ratios of ring radii are nearly $\epsilon$
independent and follow the approximate law: $a_n \sim 1/n$.  The
self-similar model also predicts for each ring the prefactor
${dM \over d\Omega dt}~ {2 \over v}$ that appears in Eq.(1).  For
example, for $\epsilon = 0.2$,
$\{{dM_n \over d\Omega dt}~{2 \over v_n}:n=1,2 ...\} =
(26,~11,~7,~5,~4,~...)10^{-2}~v_{rot}^2/4 \pi G$.

The amount of angular momentum is related to the effective core radius
of the galactic halo \cite{twi}.  For galaxies like our own, the average
of the dimensionless angular momentum distribution $\bar j \simeq 0.2$.
For a given $j$ distribution, $\bar j$ and $j_{max}$ are related by
some numerical factor.  If the turnaround sphere is initially rigidly
rotating, $j_{max} = {4 \over \pi} \bar j$.

If the caustic rings lie close to the galactic plane, they may manifest
themselves as bumps in the rotation curve.  Galactic rotation curves
often do have bumps.  Of special interest here are those which occur
at radii larger than the disc radius because they cannot readily
be attributed to inhomogeneities in the luminous matter distribution.
Consider the rotation curve of NGC 3198 \cite{vA},  one of the best
measured and often cited as providing compelling evidence for the
existence of dark halos.  It appears to have bumps near 28, 13.5 and
9 kpc, assuming $h = 0.75$.  Although the statistical significance of
these bumps is not great, let's assume for the moment that they are real
effects.  Note then that their existence is inconsistent with the
assumption that the dark halo is a perfect isothermal sphere.  On the
other hand, the radii at which they occur are in close agreement with
the ratios predicted by the self-similar model assuming that the bumps
are caused by the gravitational fields of the first three caustic rings
of NGC 3198.  Since $v_{rot} = 150$ km/sec, we find that $j_{max} = 0.28$
in this case if $\epsilon = 0.3$.  The uncertainty in $h$ drops out.  A
fit of the infall model to our own galaxy \cite{twi} produced
$\bar j \simeq 0.2$ for $\epsilon = 0.2$ to $0.3$.  If the turnaround
sphere is initially rigidly rotating, one has
$j_{max} = {4\over \pi}\bar j$.  Thus the values of $\bar j$ for our
own halo and that of NGC 3198 are found to be similar.

Let's discuss the implications of the self-similar model for our own halo.
As an example, we use the model parameters $\epsilon = 0.28$,
$j_{max} = 0.25$ and $h=0.7$.  Table 1 shows the caustic ring radii
$a_n$, the local velocities $\vec v_n$ and the local densities $d^{\pm}_n$,
at the position of the Sun, associated with the first 20 pairs of flows.
$\hat z$ is the direction perpendicular to the galactic plane, $\hat r$ is
the radial direction in the galactic plane and $\hat \phi$ is the direction
of galactic rotation.  The velocities are given in the rest frame of the
Galaxy.  $d^{\pm}_n$ is the density of {\it each} of the two $n$th flows.
We at 8.5 kpc are between the 4th and 5th ring.  The local densities of the
4th and 5th flows are large because of our proximity to the corresponding
rings.

\begin{table}[t]
\caption{Caustic ring radii, local velocities and local densities of the
first 20 pairs of flows in the self-similar infall model with
$\epsilon = 0.28, j_{max} = 0.25$ and $h = 0.7$~.}
\vspace{0.2cm}
\begin{center}
\footnotesize
\begin{tabular}{|r|c|c|c|c|c|c|}
\hline
$n$ & $a_n$ & $v_n$ & $v_{n\phi}$ & $v_{nz}$ & $v_{nr}$ & $d_n^{\pm}$\\
& (kpc) & (km/s) & (km/s) & (km/s) & (km/s) & ($10^{-26}$ gr/cm$^3$)\\
\hline
1 & 38. & 620 & 140 & $\pm$605 & ~~0~~ & 0.4 \\
2 & 19. & 565 & 255 & $\pm$505 & ~~0~~ & 1.0 \\
3 & 13. & 530 & 350 & $\pm$390 & ~~0~~ & 2.0 \\
4 & 9.7 & 500 & 440 & $\pm$240 & ~~0~~ & 6.3 \\
5 & 7.8 & 480 & 440 & ~~0~~ & $\pm$190 & 9.2 \\
6 & 6.5 & 460 & 355 & ~~0~~ & $\pm$295 & 2.9 \\
7 & 5.6 & 445 & 290 & ~~0~~ & $\pm$330 & 1.9 \\
8 & 4.9 & 430 & 250 & ~~0~~ & $\pm$350 & 1.4 \\
9 & 4.4 & 415 & 215 & ~~0~~ & $\pm$355 & 1.1 \\
10 & 4.0 & 400 & 190 & ~~0~~& $\pm$355 & 1.0 \\
11 & 3.6 & 390 & 170 & ~~0~~& $\pm$355 & 0.9 \\
12 & 3.3 & 380 & 150 & ~~0~~& $\pm$350 & 0.8 \\
13 & 3.1 & 370 & 135 & ~~0~~& $\pm$345 & 0.7 \\
14 & 2.9 & 360 & 120 & ~~0~~& $\pm$340 & 0.6 \\
15 & 2.7 & 350 & 110 & ~~0~~& $\pm$330 & 0.6 \\
16 & 2.5 & 340 & 100 & ~~0~~& $\pm$325 & 0.55 \\
17 & 2.4 & 330 & ~90 & ~~0~~& $\pm$320 & 0.50 \\
18 & 2.2 & 320 & ~85 & ~~0~~& $\pm$310 & 0.50 \\
19 & 2.1 & 315 & ~80 & ~~0~~& $\pm$305 & 0.45 \\
20 & 2.0 & 310 & ~75 & ~~0~~& $\pm$300 & 0.45 \\
\hline
\end{tabular}
\end{center}
\end{table}

There is evidence for the 6th through 13th caustic rings in that there are
sudden rises in the inner rotation curve \cite{Cl} of our galaxy at radii
very near those listed in the Table.  Similarly there is some evidence for
the 2d and 3d caustic rings in the averaged outer rotation curve \cite{Tre}.
This will be discussed in detail in a future paper\cite{tbp}.

This study was motivated by the axion \cite{ad} and WIMP \cite{sl} dark
matter searches.  These experiments may some day measure the energy spectrum
of CDM particles on Earth.  I would like to stress here the relevance of the
model for the annual modulation of the signal in direct searches for WIMP dark
matter.  The $n$ = 2, 3 ... 8 pairs of flows have velocities exceeding the
220 km/s velocity of the Sun in the direction $\hat \phi$ of the Sun's
motion.  A simple calculation \cite{tbp,Krau} shows that they dominate the
annual modulation because of their large contributions to the local halo
density.  Thus the annual modulation of the WIMP signal has the opposite
sign in this model from that predicted by the isothermal sphere model.

The model is also relevant to calculations \cite{Berg} of the sky-map of
gamma radiation from WIMP annihilation in the galactic halo.  The
annihilation rate is proportional to the square of the density.  If the
diffuse gamma ray component seen by EGRET \cite{Dix} is attributed to
annihilations in the halo away from the disk, there ought to be hot spots
in the disk at the location of the rings.

\section*{Acknowledgments}
I am grateful to K. Freese, B. Fuchs and L. Krauss for stimulating
comments.  This work is supported in part by the U.S. Department of
Energy under contract DE-FG05-86ER40272 .

\section*{References}