First neutrino interaction candidates at the LHC

FASER$\nu$ at the CERN Large Hadron Collider (LHC) is designed to directly detect collider neutrinos for the first time and study their cross sections at TeV energies, where no such measurements currently exist. In 2018, a pilot detector employing emulsion films was installed in the far-forward region of ATLAS, 480 m from the interaction point, and collected 12.2 fb$^{-1}$ of proton-proton collision data at a center-of-mass energy of 13 TeV. We describe the analysis of this pilot run data and the observation of the first neutrino interaction candidates at the LHC. This milestone paves the way for high-energy neutrino measurements at current and future colliders.


I. INTRODUCTION
There has been a longstanding interest in detecting neutrinos produced at colliders [1][2][3][4][5][6], but to date no collider neutrino has ever been directly detected. Protonproton (pp) collisions at a center-of-mass energy of 14 TeV during LHC Run-3, with an expected integrated luminosity of 150 fb −1 , will produce a high-intensity beam of O (10 12 ) neutrinos in the far-forward direction with mean interaction energy of about 1 TeV. FASERν [7] is designed to detect these neutrinos and study their properties. The detector was approved in December 2019, will be installed 480 m downstream of the ATLAS interaction point (IP) in 2021, and will take data starting in 2022. Deployment on the beam collision axis maximizes the flux of all three neutrino flavors, and allows FASERν to measure their interaction cross sections in the currently unexplored TeV energy range. For electron and tau neutrinos, these measurements will extend existing cross section measurements to significantly higher energies. For muon neutrinos, they will probe the gap between accelerator measurements (E ν < 360 GeV) [8] and IceCube data (E ν > 6.3 TeV) [9]. Following the FASERν approval, the SND@LHC experiment [10], designed to measure neutrinos at the LHC in a complementary rapidity region to FASERν, was approved in 2021.
In 2018, we performed a pilot run in the LHC tunnel to measure background and demonstrate neutrino detection at the LHC for the first time. Although the pilot detector lacked the ability to identify muons, given its depth of only 0.6λ int , much shorter than the 8λ int of the full FASERν detector, these data have aided reconstruction tool development and proven the feasibility of neutrino measurements in this experimental environment. Here we report the detection of neutrino interaction candidates in the pilot run data.

II. THE PILOT RUN IN LHC RUN-2
In 2018, we installed a 29 kg pilot detector in the TI18 tunnel, 480 m from the ATLAS IP, to measure neutrino interactions. With respect to the ATLAS IP, the TI18 tunnel is symmetric to TI12, where FASERν will be located in LHC Run-3. Previously, we reported charged particle flux measurements made with other emulsion detectors installed in the TI12 and TI18 tunnels in 2018 [11]. Here we focus on the pilot detector.
The pilot detector is divided into a 14 kg module with 101 1-mm-thick lead plates and a 15 kg module with 120 0.5-mm-thick tungsten plates, each containing the corresponding number of emulsion films [12], as shown in Fig. 1. These emulsion films and target plates were spare parts of the NA65/DsTau experiment [13]. Each module was vacuum-packed to preserve the alignment between films and placed in a 21-cm deep acrylic chamber. The transverse dimensions of the plates and films are 12.5 cm wide and 10 cm high. The beam collision axis in TI18 was mapped out by the CERN survey team with mm precision. The two modules were placed side by side with the collision axis passing between them. The estimated uncertainty of the detector position is ±1 cm in both dimensions transverse to the collision axis. An integrated luminosity of 12.2 fb −1 with an uncertainty of 2%, measured by the ATLAS experiment [14,15], was collected during 4 weeks of data taking from September to October with pp collisions at 13-TeV center-of-mass energy. The beam half-crossing-angle was about 150 µrad vertically upwards, which moves the collision axis at the FASER location upwards by ∼7 cm. The detector temperature was stable at 17.94 • C with a standard deviation of 0.07 • C [16]. Temperature stability is important to avoid displacement of the emulsion films and metallic plates and to ensure good alignment. The entire lead module and 15% of the tungsten module were used in the following analysis; the remaining spare films in the tungsten module had data quality problems.

III. SIMULATION
Neutrinos produced in the forward direction at the LHC originate from the decay of hadrons, mainly pions, kaons, and D mesons. Light hadron production is simulated using the EPOS-LHC [17], QGSJET-II-04 [18], Sibyll 2.3c [19,20], and DPMJET-III 2017.1 [21,22] simulation tools, as implemented in the CRMC [23] package, while heavy hadron production is simulated using Sibyll 2.3c, DPMJET-III 2017.1 and Pythia 8.2 [24,25] with the Monash tune [26]. Long-lived hadrons are then propagated through the forward LHC beam pipe and magnetic fields using a dedicated simulation [27] implemented as a Rivet module [28], using the geometry and beam optics for Run-2 as modeled by BDSIM [29]. We use 13 TeV collision energy and a beam half-crossingangle of 150 µrad vertically upwards. The hadrons are decayed at multiple locations along their trajectory according to decay branching fractions and kinematics pro- vided by Pythia 8.2, and the spectra of neutrinos passing through the pilot detector are tabulated. We then use Genie [30,31] with the configuration outlined in Ref. [7] to simulate neutrino interactions.
The dominant source of background to neutrino interactions in the pilot run is inelastic interactions of neutral hadrons produced in muon photonuclear interactions upstream of the detector. The flux and spectrum of muons have been estimated by the CERN sources, targets, and interactions group, which performed FLUKA simulations [32,33]. 10 8 pp collisions were simulated, and muons were propagated to the location when the beam collision axis leaves the concrete lining of the LHC tunnel (409 m from the IP) by the FLUKA simulation. The estimated muon flux as a function of energy at the 409 m position is shown in Fig. 2. The expected uncertainty on the FLUKA flux is of the order of 50%. The muons were further propagated through 67 m of rock to reach close to the pilot detector by a Geant4 simulation [34]. The expected muon fluxes at the pilot detector position are 9.4×10 3 µ − /cm 2 /fb −1 and 3.9×10 3 µ + /cm 2 /fb −1 for E µ > 100 GeV, and 1.5 × 10 4 µ − /cm 2 /fb −1 and 9.3 × 10 3 µ + /cm 2 /fb −1 for E µ > 10 GeV.
To simulate background hadron production and interactions, Geant4 simulations of muons passing through the last 8 m of rock before reaching the pilot detector were performed. Neutral hadrons produced in the last 2 m in front of the pilot detector are the most relevant, because those produced further upstream are absorbed in the rock before reaching the detector. The average rock density around CERN is measured to be about 2.5 g/cm 3 [35]. To reproduce this density, the rock was modeled as a mixture of 41% CaCO 3 and 59% SiO 2 . 10 9 negative muons and 10 9 positive muons were simulated. Kaons and neutrons are the relevant secondary neutral hadrons produced, with a small contribution from Λ baryons. Table I shows the production rates of neutral hadrons per incident muon for negative muons and positive muons. Muon-induced neutral hadrons have a steeply falling energy spectrum; a 10 GeV minimum en-TABLE I. The production rates of neutral hadrons per incident muon with an energy threshold of 10 GeV. The difference between µ − and µ + is mainly due to the difference in the energy spectra.

Negative Muons
Positive Muons ergy threshold is applied to the simulation, since lower energy hadrons cannot satisfy the vertex reconstruction criteria used in our analysis. Neutral hadron interactions with the pilot detector were simulated by Geant4 using the FTFP BERT and QGSP BERT physics lists, which correspond to different high-energy hadronic models [36].

IV. DATA ANALYSIS
The data analysis is based on the readout of the full emulsion films by the Hyper Track Selector (HTS) system [37] with a readout speed of 0.45 m 2 /hour/layer. The HTS identifies track segments ("microtracks") in the top and bottom emulsion layers of each film. A "basetrack" is formed by linking the two microtracks on a film. Each basetrack provides a 3D coordinate, 3D vector, and energy deposit (dE/dx) estimator.
Data processing is broken up into sub-volumes with a maximum size of 2 cm × 2 cm × 25 emulsion films. A preliminary alignment between each two consecutive films (position shifts and gap) is obtained using recorded tracks. To further improve the tracking resolution, an additional alignment calibration is applied by selecting tracks crossing many plates. Track reconstruction then links basetracks on different films by correlating their positions and angles. The track density in the data sample, some 10 5 /cm 2 in a small angular space of 10 mrad, is relatively high compared to other emulsion experiments. A dedicated tracking algorithm for high density environments [13] is therefore employed on top of the software framework developed for the OPERA experiment [38].
The majority of the tracks observed in the detector are expected to be background muons and related electromagnetic showers. These background charged particles were analyzed using 10 emulsion films in the lead module. The position resolution in this data set is 0.5 µm and the angular resolution is 0.2 mrad. The observed angular distribution is peaked in the direction of the ATLAS IP. The angular spread of the peak is 2.3 mrad horizontally and 1.1 mrad vertically. The spatial distribution was uniform within the detector volume. Track detection efficiency was determined from the single film efficiency  measured for tracks penetrating 10 plates. The estimated track detection efficiency for this flux measurement was (88 ± 5)%. After the efficiency correction, the charged particle flux within 10 mrad of the angular peak, which is dominated by energetic muons, is (1.7 ± 0.1) × 10 4 tracks/cm 2 /fb −1 normalized by luminosity. This result is consistent with the values previously reported by other detectors [7,16] and close to the FLUKA prediction of 2.5 × 10 4 tracks/cm 2 /fb −1 for E µ > 10 GeV.
For the neutrino analysis, using reconstructed tracks passing through at least 3 plates, vertex reconstruction was performed by searching for converging patterns of tracks with a minimum distance within 5 µm. Converging patterns with 5 or more tracks were then identified as vertices, rejecting the photon background. Collimation cuts were applied to these vertices to select high-energy interactions and suppress neutral hadron backgrounds: (1) the number of tracks with tan θ ≤ 0.1 with respect to the beam direction is required to be 5 or more, and (2) the number of tracks with tan θ > 0.1 with respect to the beam direction is required to be 4 or less. Vertices are categorized as charged or neutral based on the presence or absence, respectively, of charged parent tracks. A looser track selection is used for the charged parent track search with a higher track detection efficiency of 99.8 +0.1 −0.3 %. The background from charged vertices being reconstructed as neutral vertices is therefore negligible. The estimated selection efficiencies for neutrino signal and neutral hadron background vertices are shown in Table II. Signal classification is not performed in this analysis and interactions of all neutrino flavors are combined in the data.
The fiducial volume is defined by removing 7 films upstream, 5 films downstream, and 5 mm from the sides of the detector, corresponding to an 11 kg target mass. Within this volume, 18 neutral vertices passed the vertex selection criteria. Fig. 3 shows two selected neutral vertices in lead, with 11 and 9 associated charged particles, respectively.
The expected number of neutrino signal vertices after all selections is 3.3 +1.7 −0.9 , dominated by muon neutrino interactions. The uncertainty reflects only the range of predictions obtained from different Monte Carlo simulations. The expected numbers of neutral hadron background ver- tices are 11.0 (FTFP BERT) and 10.1 (QGSP BERT).
Since the difference of the two physics lists is not significant, FTFP BERT is used in the following analysis.
To validate the interaction features for the multivariate analysis described later, charged vertices (vertices with charged parent tracks attached), which simulation studies show also originate from muons, were analyzed. Although our muon flux measurement is close to the FLUKA prediction, no estimate of the uncertainty of the muon energy spectrum is available. The expected number of charged vertices satisfying the selection criteria is 115.4 (40.4 charged hadron interactions and 75.0 muon interactions), compared to the 78 charged vertices observed in the data.
Since the pilot detector lacked the ability to identify muons, which could allow a clean separation of neutrino charged-current and neutral hadron vertices, we introduced the following multivariate approach as a much less powerful, but necessary, alternative. A multivariate discriminant, based on a boosted decision tree (BDT) algorithm, has been developed to distinguish neutrino signal from neutral hadron background in the neutral vertex sample. The BDT was implemented using the Toolkit for Multivariate Data Analysis [39] and trained with Monte Carlo events passing the vertex selection criteria. To define input variables for the BDT analysis, we selected high-energy interactions and checked the momentum bal- ance of particles in the transverse plane. Using track angles with respect to the collision axis (θ), the following variables were defined: (1) the number of tracks with tan θ ≤ 0.1 with respect to the beam direction, (2) the number of tracks with 0.1 < tan θ ≤ 0.3 with respect to the beam direction, (3) the absolute value of the vector sum of transverse angles calculated considering all the tracks as unit vectors in the plane transverse to the beam direction (a sum ), (4) for each track in the event, we calculate the mean azimuthal angle between that track and all others, in the plane transverse to the beam direction, and then take the maximum value in the event (φ mean ), (5) for each track in the event, we calculate the ratio of the number of tracks with azimuthal opening angle ≤ 90 • and > 90 • in the plane transverse to the beam direction, and then take the maximum value in the event (r).
The expected distributions of the input variables for the neutrino signal and for the neutral hadron background, compared with the data, are shown in Fig. 4. The charged vertices mentioned above can be used to validate the modeling of the BDT input variables in simulated data. Fig. 5 shows that the BDT inputs for simulated charged hadron and muon interactions agree well with the charged vertex data.
The BDT estimator values for data and simulated neutral vertices are compared in Fig. 6. Here, the normalisation of the signal and background distributions is freely fitted to data, resulting in the best fit values of 6.1 and 11.9 events, respectively. The vertices shown in Fig. 3 correspond to the first and second largest BDT values. An excess of events over the background expectation is observed at high BDT estimator, which is in agreement with the background plus signal hypothesis. A hypothesis test using the RooStats tool implemented in the CERN ROOT framework [40] is carried out on the binned BDT estimator distribution. The backgroundonly hypothesis is disfavored with a statistical significance of 2.7σ. The expected significance is estimated with pseudo experiments with the signal expectation of 3.3 events to be 1.7σ.
A systematic uncertainty related to the shape of the BDT distribution for neutrino events was estimated by varying the generator used for neutrino production and redoing the analysis. This resulted in a small (<0.2 events) change in the fitted neutrino yield.
Systematic uncertainties on the shape of the background BDT distribution were also evaluated by varying the shape of the muon energy distribution, by varying the modeling of the photonuclear interactions in the rock that produce the background neutral hadrons from the incoming muons, and varying the physics lists for the hadron interactions. These effects can change the energy and type of the neutral hadron interacting in the detector and therefore can influence the shape of the BDT distributions. The muon distribution was scaled up and down by a factor (1+E/3 TeV) distorting the spectrum as a function of energy, and the analysis repeated. Fitting the data with the updated background BDT shapes changed the fitted neutrino yield by 0.1 events. In addition, the analysis was repeated using FLUKA to model the production of neutral hadrons instead of Geant4, this leads to a change in the fitted neutrino yield of 0.1 events. Also, the analysis was repeated using the physics list QGSP BERT to model the hadron interactions instead of FTFP BERT, this leads to a change in the fitted neutrino yield of 0.1 events.

V. CONCLUSIONS AND OUTLOOK
A search for neutrino interactions is presented based on a small emulsion detector installed at the LHC in 2018. We observe the first candidate vertices consistent with neutrino interactions at the LHC. A 2.7σ excess of neutrino-like signal above muon-induced backgrounds is measured. These results demonstrate FASERν's ability to detect neutrinos at the LHC and pave the way for future collider neutrino experiments.
We are currently preparing for data taking in LHC Run-3. With a deeper detector and lepton identification capability, FASERν will perform better than the pilot run detector. In addition, the FASER spectrometer will measure the muon flux, reducing uncertainties on background estimates. In the 2022-2024 run, we expect to collect ∼10,000 flavor-tagged charged-current neutrino interactions.