# First Direct Observation of Collider Neutrinos with FASER at the LHC
FASER Collaboration

Henso Abreu  0000-0002-1599-2896 Department of Physics and Astronomy,
Technion—Israel Institute of Technology, Haifa 32000, Israel  John Anders
0000-0002-1846-0262 CERN, CH-1211 Geneva 23, Switzerland  Claire Antel
0000-0001-9683-0890 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Akitaka Ariga
0000-0002-6832-2466 Albert Einstein Center for Fundamental Physics, Laboratory
for High Energy Physics, University of Bern, Sidlerstrasse 5, CH-3012 Bern,
Switzerland  Department of Physics, Chiba University, 1-33 Yayoi-cho Inage-ku,
263-8522 Chiba, Japan  Tomoko Ariga  0000-0001-9880-3562 Kyushu University,
Nishi-ku, 819-0395 Fukuoka, Japan  Jeremy Atkinson  0009-0003-3287-2196 Albert
Einstein Center for Fundamental Physics, Laboratory for High Energy Physics,
University of Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland  Florian U.
Bernlochner  0000-0001-8153-2719 Universität Bonn, Regina-Pacis-Weg 3, D-53113
Bonn, Germany  Tobias Blesgen  0009-0007-2013-1031 Universität Bonn, Regina-
Pacis-Weg 3, D-53113 Bonn, Germany  Tobias Boeckh  0009-0000-7721-2114
Universität Bonn, Regina-Pacis-Weg 3, D-53113 Bonn, Germany  Jamie Boyd
0000-0001-7360-0726 CERN, CH-1211 Geneva 23, Switzerland  Lydia Brenner
0000-0001-5350-7081 Nikhef National Institute for Subatomic Physics, Science
Park 105, 1098 XG Amsterdam, Netherlands  Franck Cadoux  Département de
Physique Nucléaire et Corpusculaire, University of Geneva, CH-1211 Geneva 4,
Switzerland  David W. Casper  0000-0002-7618-1683 Department of Physics and
Astronomy, University of California, Irvine, CA 92697-4575, USA  Charlotte
Cavanagh  0009-0001-1146-5247 University of Liverpool, Liverpool L69 3BX,
United Kingdom  Xin Chen  0000-0003-4027-3305 Department of Physics, Tsinghua
University, Beijing, China  Andrea Coccaro  0000-0003-2368-4559 INFN Sezione
di Genova, Via Dodecaneso, 33–16146, Genova, Italy  Ansh Desai
0000-0002-5447-8304 University of Oregon, Eugene, OR 97403, USA  Sergey
Dmitrievsky  0000-0003-4247-8697 Affiliated with an international laboratory
covered by a cooperation agreement with CERN.  Monica D’Onofrio
0000-0003-2408-5099 University of Liverpool, Liverpool L69 3BX, United Kingdom
Yannick Favre  Département de Physique Nucléaire et Corpusculaire, University
of Geneva, CH-1211 Geneva 4, Switzerland  Deion Fellers  0000-0002-0731-9562
University of Oregon, Eugene, OR 97403, USA  Jonathan L. Feng
0000-0002-7713-2138 Department of Physics and Astronomy, University of
California, Irvine, CA 92697-4575, USA  Carlo Alberto Fenoglio
0009-0007-7567-8763 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Didier Ferrere
0000-0002-5687-9240 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Stephen Gibson
0000-0002-1236-9249 Royal Holloway, University of London, Egham, TW20 0EX,
United Kingdom  Sergio Gonzalez-Sevilla  0000-0003-4458-9403 Département de
Physique Nucléaire et Corpusculaire, University of Geneva, CH-1211 Geneva 4,
Switzerland  Yuri Gornushkin  0000-0003-3524-4032 Affiliated with an
international laboratory covered by a cooperation agreement with CERN.  Carl
Gwilliam  0000-0002-9401-5304 University of Liverpool, Liverpool L69 3BX,
United Kingdom  Daiki Hayakawa  0000-0003-4253-4484 Department of Physics,
Chiba University, 1-33 Yayoi-cho Inage-ku, 263-8522 Chiba, Japan  Shih-Chieh
Hsu  0000-0001-6214-8500 Department of Physics, University of Washington, PO
Box 351560, Seattle, WA 98195-1460, USA  Zhen Hu  0000-0001-8209-4343
Department of Physics, Tsinghua University, Beijing, China  Giuseppe Iacobucci
0000-0001-9965-5442 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Tomohiro Inada
0000-0002-6923-9314 Department of Physics, Tsinghua University, Beijing, China
Sune Jakobsen  0000-0002-6564-040X CERN, CH-1211 Geneva 23, Switzerland  Hans
Joos  0000-0003-4313-4255 CERN, CH-1211 Geneva 23, Switzerland  II.
Physikalisches Institut, Universität Göttingen, Göttingen, Germany  Enrique
Kajomovitz  0000-0002-8464-1790 Department of Physics and Astronomy,
Technion—Israel Institute of Technology, Haifa 32000, Israel  Hiroaki Kawahara
0009-0007-5657-9954 Kyushu University, Nishi-ku, 819-0395 Fukuoka, Japan  Alex
Keyken  Royal Holloway, University of London, Egham, TW20 0EX, United Kingdom
Felix Kling  0000-0002-3100-6144 Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany  Daniela Köck  0000-0002-9090-5502
University of Oregon, Eugene, OR 97403, USA  Umut Kose  0000-0001-5380-9354
CERN, CH-1211 Geneva 23, Switzerland  Rafaella Kotitsa  0000-0002-7886-2685
CERN, CH-1211 Geneva 23, Switzerland  Susanne Kuehn  0000-0001-5270-0920 CERN,
CH-1211 Geneva 23, Switzerland  Helena Lefebvre  0000-0002-7394-2408 Royal
Holloway, University of London, Egham, TW20 0EX, United Kingdom  Lorne
Levinson  0000-0003-4679-0485 Department of Particle Physics and Astrophysics,
Weizmann Institute of Science, Rehovot 76100, Israel  Ke Li
0000-0002-2545-0329 Department of Physics, University of Washington, PO Box
351560, Seattle, WA 98195-1460, USA  Jinfeng Liu  Department of Physics,
Tsinghua University, Beijing, China  Jack MacDonald  0000-0002-3150-3124
Institut für Physik, Universität Mainz, Mainz, Germany  Chiara Magliocca
0009-0009-4927-9253 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Fulvio Martinelli
0000-0003-4221-5862 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Josh McFayden
0000-0001-9273-2564 Department of Physics & Astronomy, University of Sussex,
Sussex House, Falmer, Brighton, BN1 9RH, United Kingdom  Matteo Milanesio
0000-0001-8778-9638 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Dimitar Mladenov
0000-0002-0692-0495 CERN, CH-1211 Geneva 23, Switzerland  Théo Moretti
0000-0001-7065-1923 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Magdalena Munker
0000-0003-2775-3291 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland  Mitsuhiro Nakamura
Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan  Toshiyuki
Nakano  Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan
Marzio Nessi  0000-0001-7316-0118 Département de Physique Nucléaire et
Corpusculaire, University of Geneva, CH-1211 Geneva 4, Switzerland  CERN,
CH-1211 Geneva 23, Switzerland  Friedemann Neuhaus  0000-0002-3819-2453
Institut für Physik, Universität Mainz, Mainz, Germany  Laurie Nevay
0000-0001-7225-9327 CERN, CH-1211 Geneva 23, Switzerland  Royal Holloway,
University of London, Egham, TW20 0EX, United Kingdom  Hidetoshi Otono
0000-0003-0760-5988 Kyushu University, Nishi-ku, 819-0395 Fukuoka, Japan  Hao
Pang  0000-0002-1946-1769 Department of Physics, Tsinghua University, Beijing,
China  Lorenzo Paolozzi  0000-0002-9281-1972 Département de Physique Nucléaire
et Corpusculaire, University of Geneva, CH-1211 Geneva 4, Switzerland  CERN,
CH-1211 Geneva 23, Switzerland  Brian Petersen  0000-0002-7380-6123 CERN,
CH-1211 Geneva 23, Switzerland  Francesco Pietropaolo  CERN, CH-1211 Geneva
23, Switzerland  Markus Prim  0000-0002-1407-7450 Universität Bonn, Regina-
Pacis-Weg 3, D-53113 Bonn, Germany  Michaela Queitsch-Maitland
0000-0003-4643-515X University of Manchester, School of Physics and Astronomy,
Schuster Building, Oxford Rd, Manchester M13 9PL, United Kingdom  Filippo
Resnati  0000-0003-1434-6435 CERN, CH-1211 Geneva 23, Switzerland  Hiroki
Rokujo  Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan  Elisa
Ruiz-Choliz  0000-0002-2417-7121 Institut für Physik, Universität Mainz,
Mainz, Germany  Jorge Sabater-Iglesias  0000-0003-2328-1952 Département de
Physique Nucléaire et Corpusculaire, University of Geneva, CH-1211 Geneva 4,
Switzerland  Osamu Sato  0000-0002-6307-7019 Nagoya University, Furo-cho,
Chikusa-ku, Nagoya 464-8602, Japan  Paola Scampoli  0000-0001-7500-2535 Albert
Einstein Center for Fundamental Physics, Laboratory for High Energy Physics,
University of Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland  Dipartimento
di Fisica “Ettore Pancini”, Università di Napoli Federico II, Complesso
Universitario di Monte S. Angelo, I-80126 Napoli, Italy  Kristof Schmieden
0000-0003-1978-4928 Institut für Physik, Universität Mainz, Mainz, Germany
Matthias Schott  0000-0002-4235-7265 Institut für Physik, Universität Mainz,
Mainz, Germany  Anna Sfyrla  0000-0002-3003-9905 Département de Physique
Nucléaire et Corpusculaire, University of Geneva, CH-1211 Geneva 4,
Switzerland  Savannah Shively  0000-0002-4691-3767 Department of Physics and
Astronomy, University of California, Irvine, CA 92697-4575, USA  Yosuke Takubo
0000-0002-3143-8510 Institute of Particle and Nuclear Studies, KEK, Oho 1-1,
Tsukuba, Ibaraki 305-0801, Japan  Noshin Tarannum  0000-0002-3246-2686
Département de Physique Nucléaire et Corpusculaire, University of Geneva,
CH-1211 Geneva 4, Switzerland  Ondrej Theiner  0000-0002-6558-7311 Département
de Physique Nucléaire et Corpusculaire, University of Geneva, CH-1211 Geneva
4, Switzerland  Eric Torrence  0000-0003-2911-8910 University of Oregon,
Eugene, OR 97403, USA  Serhan Tufanli  CERN, CH-1211 Geneva 23, Switzerland
Svetlana Vasina  0000-0003-2775-5721 Affiliated with an international
laboratory covered by a cooperation agreement with CERN.  Benedikt Vormwald
0000-0003-2607-7287 CERN, CH-1211 Geneva 23, Switzerland  Di Wang
0000-0002-0050-612X Department of Physics, Tsinghua University, Beijing, China
Eli Welch  0000-0001-6336-2912 Department of Physics and Astronomy, University
of California, Irvine, CA 92697-4575, USA  Stefano Zambito
0000-0002-4499-2545 Département de Physique Nucléaire et Corpusculaire,
University of Geneva, CH-1211 Geneva 4, Switzerland

(March 24, 2023)

###### Abstract

We report the first direct observation of neutrino interactions at a particle
collider experiment. Neutrino candidate events are identified in a
$13.6\text{\,}\mathrm{TeV}$ center-of-mass energy $pp$ collision data set of
$35.4\text{\,}{\mathrm{fb}}^{-1}$using the active electronic components of the
FASER detector at the Large Hadron Collider. The candidates are required to
have a track propagating through the entire length of the FASER detector and
be consistent with a muon neutrino charged-current interaction. We infer
$153^{+12}_{-13}$ neutrino interactions with a significance of $16$ standard
deviations above the background-only hypothesis. These events are consistent
with the characteristics expected from neutrino interactions in terms of
secondary particle production and spatial distribution, and they imply the
observation of both neutrinos and anti-neutrinos with an incident neutrino
energy of significantly above 200 GeV.

††preprint: CERN-EP-2023-056

© 2023 CERN for the benefit of the FASER Collaboration. Reproduction of this
article or parts of it is allowed as specified in the CC-BY-4.0 license.

#### Introduction

Since their discovery at a nuclear reactor in 1956 Cowan:1992xc , neutrinos
have been detected from a variety of sources: fixed target experiments
Danby:1962nd , cosmic ray interactions in the atmosphere Achar:1965ova ;
Reines:1965qk ; Fukuda:1998mi , the Sun Davis:1968cp ; Cleveland:1998nv , the
Earth Araki:2005qa , supernovae Hirata:1987hu ; Bionta:1987qt , and
astrophysical bodies outside our galaxy Aartsen:2014gkd . With each new source
has come new insights, with important implications for many fields, from
particle physics to geophysics to astrophysics and cosmology.

Until now, however, no neutrino produced at a particle collider has ever been
directly detected. Colliders copiously produce both neutrinos and anti-
neutrinos of all flavors, and they do so in a range of very high energies
where neutrino interactions have not yet been observed. Nevertheless, collider
neutrinos have escaped detection, because they interact extremely weakly, and
the highest energy neutrinos, which have the largest probability of
interacting, are dominantly produced in the forward region, parallel to the
beamline DeRujula:1984pg ; DeRujula:1984ns ; Vannucci:253670 ; DeRujula:1992sn
; Park:2011gh ; Feng:2017uoz ; FASER:2019dxq , where collider detectors
typically have uninstrumented regions to allow the entry and exit of the
colliding particle beams. In 2021, the FASER Collaboration identified the
first collider neutrino candidates FASER:2021mtu  using a
$29\text{\,}\mathrm{kg}$ pilot detector, highlighting the potential of
discovering collider neutrinos in the forward region of the Large Hadron
Collider (LHC) collisions. In addition to FASER, the SND@LHC experiment is
expected to observe and study neutrinos produced in the LHC collisions
Beni:2019gxv ; SHiP:2020sos ; SNDLHC:2022ihg  and recently reported
preliminary findings BLAH . The observation of collider neutrino interactions
will have wide ranging implications for the study of neutrino properties, QCD,
astroparticle physics, and searches for physics beyond the Standard Model
Feng:2022inv .

This letter reports the first direct observation of neutrinos produced at a
particle collider by analyzing $35.4\text{\,}{\mathrm{fb}}^{-1}$of proton-
proton ($pp$) collision data from Run 3 of the LHC at a center-of-mass energy
of $13.6\text{\,}\mathrm{TeV}$. Neutrinos of all flavors are produced in the
decays of light and heavy hadrons as a high-intensity beam along the collision
axis. In this paper we focus on the charged-current (CC) interactions of
$\nu_{\mu}$ and $\overline{\nu}_{\mu}$; in the following, charge conjugation
and natural units are implied throughout. The chosen analysis strategy is
designed to be independent of the simulation of the detector response and
therefore does not measure the neutrino interaction cross section, but rather
the significance of the observed number of neutrino candidate events over the
non-neutrino background. In addition to being the first collider neutrinos
ever observed, the neutrinos detected here are expected to be the most
energetic ever detected from a human source, with energies in the unexplored
range $360\text{\,}\mathrm{GeV}$-$6.3\text{\,}\mathrm{TeV}$ between fixed
target measurements Zyla:2020zbs  and astroparticle data Aartsen:2017kpd .

#### The FASER Detector

Figure 1:  Schematic side view of the FASER detector with a muon neutrino
undergoing a CC interaction in the emulsion-tungsten target.

FASER FASER:2022hcn ; Feng:2017uoz ; FASER:2018ceo ; FASER:2018bac , the
ForwArd Search ExpeRiment, is an apparatus dedicated to searching for light,
extremely weakly-interacting particles and studying neutrinos. A detailed
description can be found in Ref. FASER:2022hcn . The experiment is located in
the TI12 tunnel, which connects the Super Proton Synchrotron (SPS) and LHC
tunnels, approximately $480\text{\,}\mathrm{m}$ downstream of the ATLAS
interaction point (IP) and aligned with the collision axis line-of-sight
(LOS). Charged particles produced in the forward direction at the ATLAS IP are
deflected by LHC magnets, and FASER is also shielded from the ATLAS IP by
about $100\text{\,}\mathrm{m}$ of rock and concrete. FASER’s location
therefore ensures that a high-intensity beam of neutrinos traverses the
detector, while backgrounds are highly suppressed.

The FASER detector is partially immersed in a magnetic field and consists of a
passive tungsten-emulsion neutrino detector (FASER$\nu$), two scintillator-
based veto systems, additional scintillators for triggering, a tracking
spectrometer, a pre-shower scintillator station, and an electromagnetic
calorimeter. For the current analysis, the most important components are the
veto systems, the tracking spectrometer, and the tungsten target of the
FASER$\nu$ detector. A schematic of the FASER detector is depicted in Figure
1. The trigger and data acquisition system of FASER was designed to achieve
high efficiency and reliability FASER:2021cpr . Neutrino candidate events are
triggered by scintillator signals that exceed a preset threshold below that of
a single minimum-ionizing particle (MIP), resulting in a typical trigger rate
of 0.5-$1.3\text{\,}\mathrm{kHz}$.

FASER$\nu$ consists of 730 layers of $1.1\text{\,}\mathrm{mm}$-thick tungsten
plates interleaved with emulsion films. With a width of
$25\text{\,}\mathrm{cm}$ and a height of $30\text{\,}\mathrm{cm}$, it has a
total mass of 1.1 metric tons. Although the emulsion films provide excellent
position and angular resolution to identify CC neutrino interactions, they are
not used as the extraction, scanning and analysis is time intensive. Instead,
the FASER$\nu$ detector is used as a target for CC neutrino interactions, and
we rely on the active electronic detector components of FASER to identify
suitable muon neutrino candidates Arakawa:2022rmp .

The FASER scintillator stations are instrumental to identify suitable neutrino
candidates and veto charged particles originating from the interaction point
or from secondary interactions. The first veto system (FASER$\nu$ scintillator
station) is located in front of the FASER$\nu$ emulsion detector. It is
constructed from two modules of
$$30\text{\,}\mathrm{cm}$\times$35\text{\,}\mathrm{cm}$$,
$2\text{\,}\mathrm{cm}$-thick plastic scintillators, which are read out with
photomultiplier tubes (PMTs). The second veto system (veto scintillator
station) is located after the FASER$\nu$ emulsion detector and in front of the
first magnet. It is built from three planes of the same plastic scintillators,
arranged with a $10\text{\,}\mathrm{cm}$-thick lead block placed between the
first and second planes. The lead acts as an additional target for neutrino
interactions and to absorb or convert high-energy photons from muon
bremsstrahlung.

The tracking system consists of the interface tracking station (IFT) and the
three tracking spectrometer stations FASER:2021ljd . Each tracking station is
composed of three planes with eight ATLAS semiconductor tracker (SCT) barrel
modules Abdesselam:2006wt  per plane, arranged as two columns of four modules.
Each SCT module consists of a double-layer of single-sided silicon microstrips
with a $40\text{\,}\mathrm{mrad}$ stereo angle and an
$80\text{\,}\mathrm{\SIUnitSymbolMicro m}$ strip pitch. To identify muons from
CC interactions, only the tracking spectrometer stations are used, whereas the
IFT’s location after the tungsten-emulsion detector makes it ideal to study
remnants and secondary particles of CC deep inelastic scattering neutrino
interactions. A muon candidate traversing the full length of the spectrometer
produces 18 silicon hits. Adjacent silicon hits in the tracking stations are
combined into clusters. Between the three tracking spectrometer stations are
two $1\text{\,}\mathrm{m}$-long dipole magnets with magnetic field of 0.57 T,
with a similar $1.5\text{\,}\mathrm{m}$-long magnet in front of the
spectrometer. The magnets have an aperture of 200 mm diameter, which defines
the active transverse area of the detector, and bend charged particles in the
vertical plane. In addition, signals from the timing scintillator station,
located between the first and second magnet and in front of the first tracking
station of the spectrometer, are used. The scintillator stations in
combination with the tracking system are capable of reliably identifying
incoming charged particles passing through the full length of FASER with
inefficiencies smaller than $10^{-7}$, depending on the momentum and other
requirements in the selection.

#### Data Set and Simulated Samples

For this analysis we use data from runs with stable beam conditions collected
between July and November 2022, corresponding to a total luminosity of
$(35.4\pm 0.8)$ $\text{\,}{\mathrm{fb}}^{-1}$ ATLAS:2022hro ; ATLAS:2016fhk
after data quality selection. A detailed description of the analysis is
following; additional details are contained in Appendices A-E.

To study the detector response to neutrino interactions, we simulate
$4.3\times 10^{4}$ neutrino events corresponding to an integrated luminosity
of approximately $600\text{\,}{\mathrm{fb}}^{-1}$. The interaction with the
tungsten-emulsion detector is simulated using the GENIE event generator
Genie2010 ; Genie2015 . The neutrino energy spectra and relative flavor
composition are based on Ref. Kling:2021gos . To estimate the number of
expected neutrino events, we adjust several of the assumptions of Ref.
Kling:2021gos : we correct the center-of-mass energy, beam crossing angle, and
LOS alignment, and we use the average of the neutrino flux from the predicted
light and heavy hadron production of DPMJET Roesler:2000he ; Fedynitch:2015kcn
and SIBYLL Riehn:2019jet . The difference between the two individual
predictions and the average is 27% and is assigned as an uncertainty. All
interactions of particles traversing the FASER detector are simulated using
GEANT4 GEANT4:2002zbu .

The main background to neutrino signatures originates from high-momentum
muons. We use the energy and angular spectrum predicted by the FLUKA generator
Fluka2005 ; Fluka2015 , which includes a detailed description of the LHC
machine elements and infrastructure, to simulate a sample of $2\times 10^{6}$
muons for background studies. Two additional sources of backgrounds are
relevant: neutral hadrons produced by muon interactions in the concrete in
front of the FASER detector and geometric backgrounds from charged particles
missing the FASER$\nu$ scintillator.

We use simulated samples to study the neutral hadron backgrounds. The
contamination from geometric background events is studied using sidebands and
extrapolated into the signal region using simulations. The backgrounds from
cosmic rays and LHC beam background have been studied using events occurring
when there are no collisions, and are found to be negligible.

#### Selection and Background Rejection

We focus on identifying $\nu_{\mu}$ and $\overline{\nu}_{\mu}$ CC interactions
produced in the tungsten-emulsion detector. Such interactions will produce a
high-momentum $\mu$ that can be reconstructed in the three stations of the
FASER tracking spectrometer. In addition, we expect increased activity in the
veto and timing scintillator stations and in the IFT tracking station from
secondary particles produced in the CC interaction. To avoid unconscious bias,
a blind analysis was carried out where the event selection, background
estimations, and systematic uncertainties were fixed prior to looking at data
in the signal-enhanced region.

We select events triggered by any of the scintillators downstream of
FASER$\nu$. To discard signals from beam backgrounds and cosmic muons, we
further require a timing stamp consistent with a colliding bunch crossing
identifier. We use the FASER$\nu$ scintillator to identify backgrounds from
muons or other charged particles entering the FASER detector and reject events
that deposit a charge of more than $40\text{\,}\mathrm{pC}$ in the PMTs. Such
a charge deposition would be consistent with the presence of one or several
MIPs. We only look for CC interactions that produce a muon that traverses the
entire length of the FASER detector. The signals in the scintillators
downstream of the lead wall in the veto system, and in the calorimeter, are
therefore required to be compatible with those of a MIP. With the three
tracker stations we reconstruct events with exactly one track and require more
than 11 silicon hits in the tracking stations. The reconstructed tracks are
required to have a reasonable track fit quality, and we require the
reconstructed track momentum to fulfill
$p_{\mu}>$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$. To reject charged
particles, whose trajectory geometrically missed the FASER$\nu$ scintillator
station, we extrapolate the reconstructed track from the spectrometer back to
the IFT and FASER$\nu$ scintillator. The track’s extrapolation to the IFT must
lie within $95\text{\,}\mathrm{mm}$ of the detector’s central axis, and its
extrapolation to the FASER$\nu$ scintillator must be at a distance of
$r_{\mathrm{veto}\,\nu}<$ $120\text{\,}\mathrm{mm}$ from the FASER$\nu$
scintillator center.

#### Neutral Hadron and Geometric Backgrounds

To estimate the number of neutral hadrons that reach FASER, we simulate
$2.1\text{\times}{10}^{9}\text{\,}\ $$\mu$ events based on the FLUKA energy
spectrum, and use GEANT4 to propagate through the last $8\text{\,}\mathrm{m}$
of rock in front of FASER. From this sample we determine the number of neutral
hadrons with a momentum larger than
$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ that reach the detector. The
selection efficiency is evaluated with an additional sample of neutral kaons
and neutrons with momenta larger than
$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ in front of the FASER$\nu$
emulsion detector. Most simulated hadrons are absorbed in the tungsten or do
not produce a charged track with sufficient momentum to pass the signal
selection and only a small fraction of the simulated hadrons pass all
selection steps. From this we estimate the total neutral hadron background to
be $n_{\mathrm{had}}=0.11\pm 0.06$, with the uncertainty denoting the
statistical error. Further simulation studies show that in most cases the
parent muon enters the detector along with the neutral hadron. Such events
would be rejected by the FASER$\nu$ veto scintillator. The estimate assumes
that all neutral hadron events are not already vetoed by the accompanying
muon, and is therefore a conservative estimate of this background
contribution.

To estimate the geometric background contribution, we count the number of
background events $n_{\mathrm{geo}}$ using a sideband and apply a scaling to
the signal region of $f_{\mathrm{geo}}$, which is extracted from simulated
samples. The sideband is defined to enhance the contribution of muons that
miss the FASER$\nu$ scintillator station, but may be able to produce a track
in the spectrometer, which passes the selection by scattering in the tungsten
and/or bending in the magnetic field. We modify the event selection outlined
above: we require at most 8 IFT clusters, an extrapolated radius
$r_{\mathrm{IFT}}$ of $90\text{\,}\mathrm{mm}$ to $95\text{\,}\mathrm{mm}$
with respect to the IFT center, and apply no selection on
$r_{\mathrm{veto}\,\nu}$. None of the selected sideband events have a momentum
larger than $100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. We thus extrapolate
to the signal region by using a linear fit to the momentum distribution. We
correct this estimate to account for the $r_{\mathrm{veto}\,\nu}$ selection by
using the ratio of events with a radius smaller than
$120\text{\,}\mathrm{m}\mathrm{m}$ over all sideband events in the fitted
range. As we observe no events with
$r_{\mathrm{veto}\,\nu}<$120\text{\,}\mathrm{m}\mathrm{m}$$ in the fitted
range, we use the 3$\sigma$ upper limit of the expectation value of a Poisson
process for an observation of zero events of 5.9. With this we find
$n_{\mathrm{geo}}=0.01\pm 0.23$ background events in the sideband, with the
uncertainty denoting the statistical error. We extract a scaling factor
between this sideband and the signal region from simulations, probing
different momenta, angles, and position ranges, and use the resulting
deviation from the nominal simulation scenario as an uncertainty. This results
in a scaling factor of $f_{\mathrm{geo}}=7.9\pm 2.4$ and a total geometric
background estimate of $0.08\pm 1.83$ events.

#### Results

Figure 2:  The selected signal region in extrapolated radius $r_{\mathrm{veto}\,\nu}$ and reconstructed track momentum $p_{\mu}$ is depicted. The region with lower momenta and larger radii is dominated by background events consisting of charged particles that miss the FASER$\nu$ scintillator station.  Category | Events | Expectation
---|---|---
Signal | 153 | $n_{\nu}+n_{b}\cdot p_{1}\cdot p_{2}+n_{\mathrm{had}}+n_{\mathrm{geo}}\cdot f_{\mathrm{geo}}$
$n_{10}$ | 4 | $n_{b}\cdot(1-p_{1})\cdot p_{2}$
$n_{01}$ | 6 | $n_{b}\cdot p_{1}\cdot(1-p_{2})$
$n_{2}$ | 64’014’695 | $n_{b}\cdot(1-p_{1})\cdot(1-p_{2})$
Table 1:  Observed event yields in $35.4\text{\,}{\mathrm{fb}}^{-1}$of
collision data and their relation to neutrino and background events.

Figure 2 shows the selected events, as well as the background-enriched regions
with lower momentum or $r_{\mathrm{veto}\,\nu}>$ $120\text{\,}\mathrm{mm}$. In
total we observe 153 events passing all selection steps. Using GENIE we study
the composition of neutrino events passing this selection and find that 99%
originate from muon neutrino CC interactions.

We group the selected events into four categories to estimate the number of
neutrino ($n_{\nu}$) and background events ($n_{b}$). The categorization is
determined by whether the events pass or fail the FASER$\nu$ veto scintillator
selection criteria. This allows us to determine in a simulation-independent
way the inefficiencies of the two layers of the FASER$\nu$ veto scintillator
($p_{1}$, $p_{2}$) under the assumption that they are uncorrelated.

Besides the signal category, we select:

  * $n_{10}$:

Events for which the first layer of the FASER$\nu$ scintillator produces a
charge of $>$$40\text{\,}\mathrm{pC}$ in the PMT, but no signal with
sufficient charge is seen in the second layer.

  * $n_{01}$:

Analogous events for which more than $40\text{\,}\mathrm{pC}$ in the PMT was
observed in the second layer, but not in the first layer.

  * $n_{2}$:

Events for which both layers observe more than $40\text{\,}\mathrm{pC}$ of
charge.

Table 1 lists the observed event yields and their relation to the expected
number of neutrino and background events and the FASER$\nu$ veto scintillator
inefficiencies.

We analyze the observed number of events using a binned extended maximum
likelihood fit, implemented using the iminuit package iminuit . We introduce
nuisance parameters to constrain the estimated background events to their
expectations using Gaussian priors. The likelihood is numerically maximized,
and we use a discovery test statistic Cowan:2010js  to determine the
significance of the observed signal over the background-only hypothesis. We
find

$\displaystyle n_{\nu}$
$\displaystyle=153\,{}^{+12}_{-13}\,(\mathrm{stat.})\,\,{}^{+2}_{-2}\,(\mathrm{bkg.})=153\,\,{}^{+12}_{-13}\,(\mathrm{tot.})\,$

with a significance of $16$ standard deviations over the background-only
hypothesis and based on the asymptotic distribution of the test statistic. The
excess is compatible with the expected number of neutrino events
$n_{\nu}^{\mathrm{exp}}=151\pm 41$, but note that its error does not include
any systematic uncertainties from simulating the detector response and
selection. The determined inefficiencies of the two FASER$\nu$ scintillators
are $p_{1}=(6^{+4}_{-3})\times 10^{-8}$ and $p_{2}=(9^{+4}_{-3})\times
10^{-8}$, showing values close to the expected performance FASER:2022hcn .

Figure 3:  Extrapolated transverse position of the reconstructed tracks of
neutrino-like events to the FASER$\nu$ scintillator station. The ATLAS LOS is
indicated with a red marker and shifted $59\text{\,}\mathrm{mm}$ in the
negative $y$ direction from the center of the scintillator station.

We expect that the identified neutrino candidates are distributed around the
ATLAS LOS and do not cluster at a specific point of origin. We test this by
using the extrapolated position to the FASER$\nu$ scintillator station from
the reconstructed tracks of the neutrino-like events in the signal category.
Figure 3 shows the extrapolated positions and we observe the expected
behaviour.

Figure 4:  The figures depict the number of reconstructed clusters in the IFT,
track polar angle $\theta_{\mu}$, $q/p_{\mu}$, and the reconstructed momentum
$p_{\mu}$ for events in the signal region (black markers) and compare them to
the expectation from GENIE (blue) and muon-like events (grey markers). The
muon-like events are from the $n_{2}$ category, for which both layers of the
FASER$\nu$ scintillator observed a signal, and show the expected distributions
for non-neutrino backgrounds. The blue bands correspond to the statistical
error of the simulated samples and are luminosity scaled for $q/p_{\mu}$ and
$p_{\mu}$. The other figures are normalized to unity.

Figure 4 summarizes additional properties of the signal category events. The
CC neutrino interactions produce on average a larger number of particles than
MIP interactions, which appear in the IFT as charge depositions. The number of
IFT clusters of the signal category is very distinct from background-like
($n_{2}$) events and agrees well with the expectation from GENIE. We also
examine the polar angles $\theta_{\mu}$ of the neutrino candidates and observe
distributions close to the simulated neutrino events and distinctively
different from muon backgrounds. We observe a clear charge separation in
$q/p_{\mu}$ for the reconstructed tracks, with $q$ denoting the assigned track
charge. In total 40 events with a positively-charged track candidate are
observed, showing the presence of anti-neutrinos in the analyzed data set. The
reconstructed momentum of the muon produced in a CC $\nu_{\mu}$ interaction is
a good proxy for the incident neutrino energy. Using the simulated CC neutrino
interactions, we estimate that with our analysis strategy we select neutrino
events for which on average $>80\%$ of the incident neutrino momentum is
transferred to the final state muon. This indicates that a large fraction of
the reconstructed neutrino candidates have energies significantly larger than
$200\text{\,}\mathrm{GeV}$. A detailed study of these properties, which
accounts for systematic effects, is left for future work.

#### Summary

We report the first direct detection of neutrinos produced at a collider
experiment using the active electronic components of the FASER detector. We
observe $153^{+12}_{-13}$ neutrino events from CC interactions from
$\nu_{\mu}$ and $\overline{\nu}_{\mu}$ taking place in the tungsten-emulsion
detector of FASER$\nu$. The spatial distribution and properties of the
observed signal events are consistent with neutrino interactions, and the
chosen analysis strategy does not depend on the quality of the modeling of
detector effects in the simulation. For the signal events, the reconstructed
charge shows the presence of anti-neutrinos, and the reconstructed momentum
implies that neutrino candidates have energies significantly above
$200\text{\,}\mathrm{GeV}$. This result marks the beginning of the field of
collider neutrino physics, opening up a wealth of new measurements with broad
implications across many physics domains Feng:2022inv .

#### Acknowledgments

We thank CERN for the very successful operation of the LHC during 2022. We
thank the technical and administrative staff members at all FASER institutions
for their contributions to the success of the FASER project. We thank the
ATLAS Collaboration for providing us with accurate luminosity estimates for
the used Run 3 LHC collision data. FASER gratefully acknowledges the donation
of spare ATLAS SCT modules and spare LHCb calorimeter modules, without which
the experiment would not have been possible. We also acknowledge the ATLAS
collaboration software, Athena, on which FASER’s offline software system is
based ATL-PHYS-PUB-2009-011  and the ACTS tracking software framework
Ai:2021ghi . Finally we thank the CERN STI group for providing detailed FLUKA
simulations of the muon fluence along the LOS, which have been used in this
analysis. This work was supported in part by Heising-Simons Foundation Grant
Nos. 2018-1135, 2019-1179, and 2020-1840, Simons Foundation Grant No. 623683,
U.S. National Science Foundation Grant Nos. PHY-2111427, PHY-2110929, and
PHY-2110648, JSPS KAKENHI Grants Nos. JP19H01909, JP20K23373, JP20H01919,
JP20K04004, and JP21H00082, BMBF Grant No. 05H20PDRC1, DFG EXC 2121 Quantum
Universe Grant No. 390833306, ERC Consolidator Grant No. 101002690, Royal
Society Grant No. URF\R1\201519, UK Science and Technology Funding Councils
Grant No. ST/ T505870/1, the National Natural Science Foundation of China,
Tsinghua University Initiative Scientific Research Program, and the Swiss
National Science Foundation.

## References

  * (1) C. L. Cowan, F. Reines, F. B. Harrison, H. W. Kruse, and A. D. McGuire, “Detection of the free neutrino: A Confirmation,” Science 124 (1956) 103–104.
  * (2) G. Danby, J. M. Gaillard, K. A. Goulianos, L. M. Lederman, N. B. Mistry, M. Schwartz, and J. Steinberger, “Observation of High-Energy Neutrino Reactions and the Existence of Two Kinds of Neutrinos,” Phys. Rev. Lett. 9 (1962) 36–44.
  * (3) C. V. Achar et al., “Detection of muons produced by cosmic ray neutrinos deep underground,” Phys. Lett. 18 (1965) 196–199.
  * (4) F. Reines, M. F. Crouch, T. L. Jenkins, W. R. Kropp, H. S. Gurr, G. R. Smith, J. P. F. Sellschop, and B. Meyer, “Evidence for high-energy cosmic ray neutrino interactions,” Phys. Rev. Lett. 15 (1965) 429–433.
  * (5) Super-Kamiokande Collaboration, Y. Fukuda et al., “Evidence for oscillation of atmospheric neutrinos,” Phys. Rev. Lett. 81 (1998) 1562–1567, arXiv:hep-ex/9807003 [hep-ex].
  * (6) R. Davis, Jr., D. S. Harmer, and K. C. Hoffman, “Search for neutrinos from the sun,” Phys. Rev. Lett. 20 (1968) 1205–1209.
  * (7) B. T. Cleveland, T. Daily, R. Davis, Jr., J. R. Distel, K. Lande, C. K. Lee, P. S. Wildenhain, and J. Ullman, “Measurement of the solar electron neutrino flux with the Homestake chlorine detector,” Astrophys. J. 496 (1998) 505–526.
  * (8) T. Araki et al., “Experimental investigation of geologically produced antineutrinos with KamLAND,” Nature 436 (2005) 499–503.
  * (9) Kamiokande-II Collaboration, K. Hirata et al., “Observation of a Neutrino Burst from the Supernova SN 1987a,” Phys. Rev. Lett. 58 (1987) 1490–1493.
  * (10) R. M. Bionta et al., “Observation of a Neutrino Burst in Coincidence with Supernova SN 1987a in the Large Magellanic Cloud,” Phys. Rev. Lett. 58 (1987) 1494.
  * (11) IceCube Collaboration, M. G. Aartsen et al., “Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data,” Phys. Rev. Lett. 113 (2014) 101101, arXiv:1405.5303 [astro-ph.HE].
  * (12) A. De Rujula and R. Ruckl, “Neutrino and muon physics in the collider mode of future accelerators,” in SSC Workshop: Superconducting Super Collider Fixed Target Physics, pp. 571–596.  5, 1984.
  * (13) A. De Rujula, “Neutrino Physics at Future Colliders,” in Prague 1984, Proceedings, Trends in Physics, Vol. 1, 236-245. 1984\.
  * (14) F. Vannucci, “Neutrino physics at LHC/SSC,” Tech. Rep. LPNHE-93-03, Univ. Paris 7, Aug, 1993.  https://cds.cern.ch/record/253670.
  * (15) A. De Rujula, E. Fernandez, and J. Gomez-Cadenas, “Neutrino fluxes at future hadron colliders,” Nucl. Phys. B 405 (1993) 80–108.
  * (16) H. Park, “The estimation of neutrino fluxes produced by proton-proton collisions at $\sqrt{s}=14$ TeV of the LHC,” JHEP 10 (2011) 092, arXiv:1110.1971 [hep-ex].
  * (17) J. L. Feng, I. Galon, F. Kling, and S. Trojanowski, “ForwArd Search ExpeRiment at the LHC,” Phys. Rev. D 97 (2018) no. 3, 035001, arXiv:1708.09389 [hep-ph].
  * (18) FASER Collaboration, “Detecting and Studying High-Energy Collider Neutrinos with FASER at the LHC,” Eur. Phys. J. C 80 (2020) no. 1, 61, 1908.02310 [hep-ex].
  * (19) FASER Collaboration, H. Abreu et al., “First neutrino interaction candidates at the LHC,” Phys. Rev. D 104 (2021) no. 9, L091101, arXiv:2105.06197 [hep-ex].
  * (20) N. Beni et al., “Physics Potential of an Experiment using LHC Neutrinos,” J. Phys. G 46 (2019) no. 11, 115008, arXiv:1903.06564 [hep-ex].
  * (21) SHiP Collaboration, C. Ahdida et al., “SND@LHC,” arXiv:2002.08722 [physics.ins-det].
  * (22) SND@LHC Collaboration, G. Acampora et al., “SND@LHC: The Scattering and Neutrino Detector at the LHC,” arXiv:2210.02784 [hep-ex].
  * (23) SND@LHC Collaboration, C. Ahdida et al., “Presentation ’Results from SND@LHC’ at 57th Rencontres de Moriond presented by Ettore Zaffaroni,”.
  * (24) J. L. Feng et al., “The Forward Physics Facility at the High-Luminosity LHC,” J. Phys. G 50 (2023) no. 3, 030501, arXiv:2203.05090 [hep-ex].
  * (25) Particle Data Group Collaboration, P. A. Zyla et al., “Review of Particle Physics,” PTEP 2020 (2020) no. 8, 083C01.
  * (26) IceCube Collaboration, M. Aartsen et al., “Measurement of the multi-TeV neutrino cross section with IceCube using Earth absorption,” Nature 551 (2017) 596–600, arXiv:1711.08119 [hep-ex].
  * (27) FASER Collaboration, H. Abreu et al., “The FASER Detector,” arXiv:2207.11427 [physics.ins-det].
  * (28) FASER Collaboration, A. Ariga et al., “Letter of Intent for FASER: ForwArd Search ExpeRiment at the LHC,” arXiv:1811.10243 [physics.ins-det].
  * (29) FASER Collaboration, A. Ariga et al., “Technical Proposal for FASER: ForwArd Search ExpeRiment at the LHC,” arXiv:1812.09139 [physics.ins-det].
  * (30) FASER Collaboration, H. Abreu et al., “The trigger and data acquisition system of the FASER experiment,” JINST 16 (2021) no. 12, P12028, arXiv:2110.15186 [physics.ins-det].
  * (31) J. Arakawa, J. L. Feng, A. Ismail, F. Kling, and M. Waterbury, “Neutrino detection without neutrino detectors: Discovering collider neutrinos at FASER with electronic signals only,” Phys. Rev. D 106 (2022) no. 5, 052011, arXiv:2206.09932 [hep-ph].
  * (32) FASER Collaboration, H. Abreu et al., “The tracking detector of the FASER experiment,” Nucl. Instrum. Meth. A 1034 (2022) 166825, arXiv:2112.01116 [physics.ins-det].
  * (33) A. Abdesselam et al., “The barrel modules of the ATLAS semiconductor tracker,” Nucl. Instrum. Meth. A 568 (2006) 642–671.
  * (34) ATLAS Collaboration, “Luminosity determination in $pp$ collisions at $\sqrt{s}=13$ TeV using the ATLAS detector at the LHC,” arXiv:2212.09379 [hep-ex].
  * (35) ATLAS Collaboration, M. Aaboud et al., “Luminosity determination in pp collisions at $\sqrt{s}$ = 8 TeV using the ATLAS detector at the LHC,” Eur. Phys. J. C 76 (2016) no. 12, 653, arXiv:1608.03953 [hep-ex].
  * (36) C. Andreopoulos et al., “The GENIE Neutrino Monte Carlo Generator,” Nucl. Instrum. Meth. A 614 (2010) 87–104, arXiv:0905.2517 [hep-ph].
  * (37) C. Andreopoulos, C. Barry, S. Dytman, H. Gallagher, T. Golan, R. Hatcher, G. Perdue, and J. Yarba, “The GENIE Neutrino Monte Carlo Generator: Physics and User Manual,” arXiv:1510.05494 [hep-ph].
  * (38) F. Kling and L. J. Nevay, “Forward neutrino fluxes at the LHC,” Phys. Rev. D 104 (2021) no. 11, 113008, arXiv:2105.08270 [hep-ph].
  * (39) S. Roesler, R. Engel, and J. Ranft, “The Monte Carlo event generator DPMJET-III,” in International Conference on Advanced Monte Carlo for Radiation Physics, Particle Transport Simulation and Applications (MC 2000), pp. 1033–1038.  12, 2000.  arXiv:hep-ph/0012252.
  * (40) A. Fedynitch, Cascade equations and hadronic interactions at very high energies.  PhD thesis, KIT, Karlsruhe, Dept. Phys., 11, 2015.
  * (41) F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev, “Hadronic interaction model Sibyll 2.3d and extensive air showers,” Phys. Rev. D 102 (2020) no. 6, 063002, arXiv:1912.03300 [hep-ph].
  * (42) GEANT4 Collaboration, S. Agostinelli et al., “GEANT4–a simulation toolkit,” Nucl. Instrum. Meth. A 506 (2003) 250–303.
  * (43) A. Ferrari et al., “FLUKA: A multi-particle transport code (Program version 2005),”.
  * (44) G. Battistoni et al., “Overview of the fluka code,” Annals of Nuclear Energy 82 (2015) 10–18. Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo 2013, SNA + MC 2013. Pluri- and Trans-disciplinarity, Towards New Modeling and Numerical Simulation Paradigms.
  * (45) H. Dembinski and P. O. et al., “scikit-hep/iminuit,”. https://doi.org/10.5281/zenodo.3949207.
  * (46) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, “Asymptotic formulae for likelihood-based tests of new physics,” Eur. Phys. J. C 71 (2011) 1554, arXiv:1007.1727 [physics.data-an]. [Erratum: Eur.Phys.J.C 73, 2501 (2013)].
  * (47) ATLAS Collaboration, “The ATLAS Collaboration Software and Firmware.” ATL-SOFT-PUB-2021-001, 2021.  https://cds.cern.ch/record/2767187.
  * (48) X. Ai et al., “A Common Tracking Software Project,” Comput. Softw. Big Sci. 6 (2022) no. 1, 8, arXiv:2106.13593 [physics.ins-det].

##  I Appendix

###  I.1 Appendix A: Geometric Sideband

Figure 5 depicts the sideband used to estimate the geometric backgrounds of
the analysis. Background events are required to be consistent with a muon
candidate by having $\leq 8$ IFT clusters and an extrapolated radius
$r_{\mathrm{IFT}}$ of $90\text{\,}\mathrm{mm}$ to $95\text{\,}\mathrm{mm}$
with respect to the IFT center. This selection is dominated by geometric
background events that do not pass the signal selection steps of the analysis.
No events with $p_{\mu}>$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ are
observed. To estimate the number of events within this momentum range, we
linearly extrapolate the events between
$30\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and
$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and find $0.2\pm 4.1$ events,
with the error denoting the statistical error. To account for the
$r_{\mathrm{veto}\,\nu}$ requirement of the signal selection, we further apply
a requirement of $r_{\mathrm{veto}\,\nu}<120\,\mathrm{mm}$ to the sideband
events (orange distribution). No events with $p_{\mu}>30\,\mathrm{GeV}$ are
observed. We thus use $5.9$ as the 3$\sigma$ upper limit and use this to
calculate the ratio with respect to the number of events without any
$r_{\mathrm{veto}\,\nu}$ selection, to correct the sideband background events
for the $r_{\mathrm{veto}\,\nu}$ requirement. With this factor we find
$n_{\mathrm{geo}}=0.01\pm 0.23$ geometric background events. To account for
the fact that this number corresponds to an annulus, the correction factor
$f_{\mathrm{geo}}=7.9\pm 2.4$, determined from simulation, is applied. It is
obtained from simulation with the uncertainty spanning different assumptions
about the angle, momenta, and positions of the geometric background events.

Figure 5:  Sideband for geometric background estimation.

###  I.2 Appendix B: Event Display

Figure 6 shows an event display of an example neutrino candidate event. The
event has a momentum of $p_{\mu}=$843.9\text{\,}\mathrm{GeV}$$, negative
charge, $\theta_{\mu}=2.5\,\text{mrad}$,
$r_{\mathrm{veto}\,\nu}=$57.2\text{\,}\mathrm{mm}$$,
$r_{\mathrm{IFT}}=$55.8\text{\,}\mathrm{mm}$$ and produced 57 clusters in the
IFT.

Figure 6:  Event display of a neutrino interaction candidate in which
secondary particles produced in the CC interaction produce activity in the
IFT.

###  I.3 Appendix C: Likelihood Fit

The used likelihood has the form

$\displaystyle\mathcal{L}=\prod_{i}\mathcal{P}(N_{i}|n_{i})\cdot\prod_{j}\mathcal{G}_{j}\,.$
(1)

Here $\mathcal{P}$ denotes a Poissonian with the index $i$ running over the
four event categories with observed event counts $N_{i}$ and expectation
values $n_{i}$. We introduce nuisance parameters to constrain the estimated
number of background events to their expectations using three Gaussian priors
$\mathcal{G}_{j}$ . The used test statistic has the form

$\displaystyle
q_{0}=\bigg{\\{}\begin{matrix}-2\ln\lambda(n_{\nu}=0)&\widehat{n}_{\nu}\geq
0\\\ 0&\widehat{n}_{\nu}<0\end{matrix}$ (2)

and the significance of the observed signal $\widehat{n}_{\nu}$ over the
background-only hypothesis is given by $\sqrt{q_{0}}$ in the asymptotic limit.
Further
$\lambda(n_{\nu}=0):=\mathcal{L}(n_{\nu})/\mathcal{L}(\widehat{n}_{\nu})$
denotes the ratio of the likelihood maximized with the condition of no signal,
$n_{\nu}=0$, to the unconditionally maximized likelihood. The log-likelihood
ratio is shown in Figure 7.

Figure 7:  The log-likelihood ratio of the estimated number of neutrinos is
shown in blue. The dashed orange contour fixes the parameters of
$\mathcal{G}_{i}$ to determine the statistical uncertainty of the neutrino
signal yield.

###  I.4 Appendix D: Momentum Resolution

Data-driven alignment corrections are applied to the positions and
orientations of the modules of the tracking spectrometer stations using a
sample of reconstructed muons. In the case of perfect alignment of the FASER
tracking detectors, we expect a momentum resolution of 2.1% at 100 GeV, 4.7%
at 300 GeV, and 16.4% at 1 TeV. The accuracy of the alignment is validated
using a photon conversion sample for momenta up to 250 GeV.

###  I.5 Appendix E: Expected Number of Neutrino Events

Volume | Type | $0<E_{\nu}<$500\text{\,}\mathrm{GeV}$$ | $500<E_{\nu}<$1000\text{\,}\mathrm{GeV}$$ | $E_{\nu}>$1000\text{\,}\mathrm{GeV}$$ | $\sum$ |  $\overline{E}_{\nu}$ [GeV]
---|---|---|---|---|---|---
FASER$\nu$ | $\nu_{\mu}$ | 359 / 379 | 239 / 273 | 291 / 790 | 890 / 1442 | 880 / 1376
FASER$\nu$ | $\overline{\nu}_{\mu}$ | 116 / 130 | 62 / 85 | 49 / 151 | 227 / 367 | 657 / 1028
$r<$95\text{\,}\mathrm{mm}$$ | $\nu_{\mu}$ | 147 / 154 | 105 / 118 | 141 / 375 | 394 / 647 | 943 / 1477
$r<$95\text{\,}\mathrm{mm}$$ | $\overline{\nu}_{\mu}$ | 48 / 53 | 28 / 37 | 23 / 67 | 99 / 157 | 687 / 1057
Table 2:  The expected numbers of neutrino and anti-neutrino events from
SIBYLL (first number) and DPMJET (second number) for an integrated luminosity
of $35.4\text{\,}{\mathrm{fb}}^{-1}$ i and different energy intervals, along
with the sum over all energy intervals, and the average neutrino energy
$\overline{E}_{\nu}$. Results are shown requiring the interactions to be (1)
in the FASER$\nu$ detector volume or (2) in the target region and within a
radius of $95\text{\,}\mathrm{mm}$ from the center of the FASER detector.

The predicted numbers of neutrino and anti-neutrino interactions from SIBYLL
and DPMJET are listed in Table 2. Results are shown requiring the interactions
to be (1) in the FASER$\nu$ detector volume or (2) in the target region and
within a radius of $95\text{\,}\mathrm{mm}$ from the center of the FASER
detector. Note that no additional acceptance and efficiency corrections are
applied and the second requirement approximates the fiducial volume used in
the analysis.


# Tomographic entanglement indicators in a coupled oscillator model

Sreelekshmi Pillai<EMAIL_ADDRESS>S. Ramanan  V. Balakrishnan
Department of Physics, Indian Institute of Technology Madras, Chennai, India
Center for Quantum Information, Communication and Computing (CQuICC), Indian
Institute of Technology Madras, Chennai, India  S. Lakshmibala  Center for
Quantum Information, Communication and Computing (CQuICC), Indian Institute of
Technology Madras, Chennai, India

###### Abstract

We study entanglement in a simple model comprising two coupled linear harmonic
oscillators of the same natural frequency. The system is separable in the
center of mass (COM) and relative coordinates into two oscillators of
frequency $\omega_{c}$ and $\omega_{r}$. We compute standard entanglement
measures (subsystem linear entropy and subsystem von Neumann entropy) as well
as several tomographic entanglement indicators (Bhattacharyya distance,
Kullback-Leibler divergence and inverse participation ratio) as functions of
the frequency ratio $\eta=\omega_{c}/\omega_{r}$, keeping the COM oscillator
in the ground state. We demonstrate that, overall, the entanglement indicators
reflect quite faithfully the variations in the standard measures. The
entanglement is shown to be minimum at $\eta=1$ and maximum as $\eta\to 0$ or
$\infty$.

Bipartite entanglement measures, tomographic entangelment indicators,
Bhattacharyya distance, Kullback-Leibler divergence, inverse participation
ratio, coupled harmonic oscillators

##  I Introduction

Entanglement is an important resource for several tasks pertaining to quantum
information and computation. Investigations of entanglement have been carried
out on discrete variable (DV) systems such as two-level atomic ensembles and
spin arrays, and on continuous variable (CV) and hybrid quantum (HQ) systems
as with atoms interacting with a radiation field. These investigations have
led to a fair degree of understanding of diverse aspects such as the sudden
death of entanglement during the time evolution of specific systems [1], the
efficacy of different entanglement indicators (see, for instance, [2]),
signatures of entanglement in avoided energy level crossings of eigenstates of
generic Hamiltonians [3] and their relation to quantum chaos [4], and the
connection between entanglement and uncertainty relations [5].

However, the links between entanglement and other properties of atomic and
nuclear systems such as the binding energy have only attracted sporadic
attention in the past, although entangled states of nucleons are well known
even in the context of the charge independence of nuclear forces. It is only
in recent years that there has been an upsurge of interest in analyzing the
importance of entanglement in bipartite atomic and nuclear systems from the
information-theoretic point of view. In this context, the roles played by the
entanglement between the spin and spatial wave functions, and by the
entanglement between the space-spin state and the isospin state in determining
the structure of a nucleus, merit detailed attention. These aspects are as yet
in a nascent stage of investigation.

Commonly used entanglement measures and monotones are defined through the
density matrix of the system concerned. The density matrix itself has to be
reconstructed from experimentally obtained tomograms which are histograms of
measured observables. The state (i.e., density matrix) reconstruction program
turns out to be a challenging task if the dimensionality of the Hilbert space
is large, as in CV and HQ systems and in certain DV systems. For instance, two
commonly used bipartite entanglement measures, the subsystem linear entropy
(SLE) and the subsystem von Neumann entropy (SVNE) which are defined in terms
of the density matrix, are not always easy to compute. In terms of the reduced
density matrix $\widetilde{\rho}$ of either subsystem of a bipartite system,
the SLE is $1-\mathrm{Tr}\,(\widetilde{\rho}^{\,2})$, while the SVNE is given
by $-\mathrm{Tr}\,(\widetilde{\rho}\ \mathrm{ln}\,\widetilde{\rho})$. Both
these quantities vanish identically for unentangled states. In general, the
Schmidt decomposition (i.e., the diagonalization of $\widetilde{\rho}$) that
is necessary for computing these measures becomes a challenging numerical
problem for 3-dimensional CV systems. This is because a nontrivial integral
equation has to be solved to obtain the eigenvalues (or Schmidt coefficients)
in a continuous basis. Moreover, in certain cases the SLE and the SVNE may not
be useful in quantifying entanglement. For instance, even in the simple
example of the electron-proton entanglement in the hydrogen atom, under
certain conditions, the SVNE diverges while SLE tends to unity (its maximum
value) for _all_ energy eigenstates. Hence these measures are unsuitable for
relating entanglement with the energy spectrum [6],[7], [8]. It is therefore
desirable and important to identify a wider set of entanglement indicators.
Early investigations on the ground state of the hydrogen atom, with the center
of mass (COM) delocalized [9], revealed that the spread in the Schmidt
coefficients mirrors the spread in the momentum of the subsystem. The latter,
which is easier to compute, could therefore be a possible indicator of
entanglement. It has been shown [8] that this spread scales inversely as the
principal quantum number for excited states. Further, if the COM is localized
as a Gaussian wave packet, the second moments of the position and momentum
could be used to identify entanglement.

Much work remains to be done in order to arrive at universally applicable good
entanglement indicators in 3-dimensional CV systems. The aim of the present
investigation is to use tomographic entanglement indicators (TEIs) that are
readily calculable from appropriate tomograms, to assess entanglement. While
such TEIs have been proposed in the literature and their efficacy examined in
specific CV and HQ systems, their utility in understanding the link between
entanglement and nuclear/atomic structure has not yet been explored in detail.
While in this work we only report on exploratory investigations carried out on
a toy model, this approach opens up the possibility of readily understanding
the role of entanglement in the formation of many-nucleon bound states, and of
examining the links between entanglement and the stability of aggregates of
nuclei, circumventing arduous state reconstruction programs.

In order to study the effects of factors such as the localization of the COM
on entanglement, we analyze a 1-dimensional model of two coupled harmonic
oscillators, which is analytically tractable, and provides information on the
performance of the TEIs. Earlier literature on such bipartite oscillators
([10], [11], [7]) primarily contains analytical results pertaining to the SVNE
as a measure of entanglement between the two individual oscillators in any of
the eigenstates of the composite system. These results corroborate the
findings in the case of the hydrogen atom, namely, that the SVNE diverges when
the COM is not localized. Our emphasis is on the TEIs, and wherever necessary
we will compare their performance with the SVNE and SLE.

The organization of this paper is as follows. Section (II) defines the TEIs
that are used in this work. Section (III) deals with the description of the
model of interest, the results obtained, and the inferences drawn. In Sec.
(IV), we summarize our main results and suggest avenue for further research.

##  II Tomographic Entanglement Indicators

A quantum tomogram is a histogram comprising measurement outcomes of
observables in different basis sets. For pure states of a given system, the
tomogram of relevance to the present work is simply a collection of the
expectation values of the density matrix in a quorum of basis sets. As a
consequence, a tomogram does not explicitly involve the off-diagonal elements
of the density operator. All the elements of the bipartite density matrix can
be obtained, however, by carrying out an inverse Radon transformation on the
tomogram. This is a consequence of the fact that, in general, the off-diagonal
elements of the density operator, and the negativity of the Wigner function
for nonclassical states of light, are captured by performing the
transformation on the collection of measurements in several basis sets [12].

The noteworthy feature in the tomographic approach is that the bipartite
entanglement in any system is quantified by computing the TEIs in different
basis sets, and then evaluating their net effect using averaging procedures.
In practice, only measurements in a system-dependent finite number of
judiciously chosen basis sets will be available. The identification of this
choice of basis sets for a given system is nontrivial. From the tomographic
point of view, the TEIs computed directly in such selected basis sets (_prior_
to state reconstruction) should display the gross properties of the SVNE or
SLE (which can be computed only _after_ state reconstruction). Reconstruction
of the wave function of a many-body system is a computationally challenging
task, even when suitable approximations are made. It is in this context that
the TEIs could be very useful, if their slice-dependent (i.e., basis-
dependent) signatures in different settings are examined in detail and
categorized. Ideal candidate systems for this purpose are readily available in
nuclear physics, as nuclear wave functions are commonly expressed solely in
the position or momentum basis (equivalently, in canonically conjugate basis
sets). Bearing these considerations in mind, we have carried out the
computation of TEIs from specific slices of the tomogram for _known_ states of
the coupled oscillators. That is, we shall use the known position and momentum
space oscillator probability distributions as marginals in the calculation of
the TEIs. The efficacy of the TEIs can then be assessed by comparison with the
computed SVNE and SLE corresponding to the states concerned. In this sense, we
are merely demonstrating the feasibility of the tomographic approach in this
work.

We recall that the tomogram is a collection of normalized probability
distributions computed in different basis sets. A TEI is expected to assess
the difference between a probability distribution corresponding to an
entangled state and the factored product of the distributions corresponding to
the unentangled subsystem states. This exercise can be carried out readily if
we express the TEIs in terms of ‘distances’ between these distributions [13].
Two examples of the latter are the Bhattacharyya distance
($\mathrm{\epsilon_{{}_{BD}}}$) and the Kullback-Liebler divergence
($\mathrm{\epsilon_{{}_{KL}}}$), commonly used in classical probability
theory. While neither qualifies to be a distance in the strict sense of the
term, they are useful as entanglement indicators. We define them briefly
below.

The Bhattacharyya distance quantifies the similarity between two probability
distributions by estimating their mutual information. For a given slice (i.e.,
in a given basis) we have

$\displaystyle\mathrm{\epsilon_{BD}}=-\mathrm{log}_{2}\\!\int\\!\\!d\bm{X}\\!\int\\!\\!d\bm{Y}\big{[}\bra{\bm{X},\bm{Y}}\rho\ket{\bm{X},\bm{Y}}\bra{\bm{X}}\widetilde{\rho}_{1}\ket{\bm{X}}\bra{\bm{Y}}\widetilde{\rho}_{2}\ket{\bm{Y}}\big{]}^{1/2},$
(1)

where $\bm{X}$ and $\bm{Y}$ are the basis labels of subsystems 1 and 2
respectively. In the example we will consider subsequently, these are the
position labels of the two subsystems. $\rho$ is the full density matrix while
$\widetilde{\rho}_{1}$, $\widetilde{\rho}_{2}$ are the reduced density
matrices. The Kullback-Leibler divergence is given by

$\displaystyle\mathrm{\epsilon_{KL}}=\int\\!\\!d\bm{X}\\!\int\\!\\!d\bm{Y}\\!\bra{\bm{X},\bm{Y}}\rho\ket{\bm{X},\bm{Y}}\mathrm{log}_{2}\frac{\bra{\bm{X},\bm{Y}}\rho\ket{\bm{X},\bm{Y}}}{\bra{\bm{X}}\widetilde{\rho}_{1}\ket{\bm{X}}\bra{\bm{Y}}\widetilde{\rho}_{2}\ket{\bm{Y}}}.$
(2)

Both $\mathrm{\epsilon_{BD}}$ and $\mathrm{\epsilon_{{}_{KL}}}$ lie in the
range [0, $\infty$), and vanish for unentangled states. Another entanglement
indicator $\mathrm{\epsilon_{{}_{IPR}}}$, based on the inverse participation
ratio (IPR), can also be computed directly from the tomogram, as it involves
elements of the density matrix in specific basis sets. The IPR was introduced
originally in the context of estimating the localization properties of states.
It indicates how ‘spread out’ the state is in a given basis, and is related to
a generalized definition of entanglement [14]. By scaling the variables and
the density matrix elements appropriately, we obtain the corresponding
dimensionless quantities. Denoting the latter with an overhead bar, the IPR-
based entanglement indicator is given by

$\mathrm{\epsilon_{{}_{IPR}}}=1+\int\\!\\!d\overline{\bm{X}}\\!\int\\!\\!d\overline{\bm{Y}}\bra{\overline{\bm{X}},\overline{\bm{Y}}}{\overline{\rho}}\ket{\overline{\bm{X}},\overline{\bm{Y}}}^{2}-\int\\!\\!d\overline{\bm{X}}\bra{\overline{\bm{X}}}{\overline{\widetilde{\rho}}}_{1}\ket{\overline{\bm{X}}}^{2}-\int\\!\\!d\overline{\bm{Y}}\bra{\overline{\bm{Y}}}{\overline{\widetilde{\rho}}}_{2}\ket{\overline{\bm{Y}}}^{2}.$
(3)

The precise transformation to the dimensionless quantities in Eq. 3 depends on
the model considered (Eqs. (26) and (27) in the case of the model we study).

##  III Entanglement between two harmonic oscillators

Consider two coupled oscillators, each with mass m and frequency $\omega$. The
Hamiltonian is given by

$\displaystyle
H=\frac{1}{2m}({p_{1}}^{2}+{{p_{2}}^{2}})+\frac{1}{2}m\omega^{2}({x_{1}}^{2}+{x_{2}}^{2})+\frac{1}{2}m\lambda{x_{1}}{x_{2}},$
(4)

where $x_{i}$ and $p_{i}$ ($i=1,2)$ are the position and momentum of the
$i^{\rm th}$ oscillator, and $\lambda$ is the coupling constant. The
Hamiltonian is positive definite if $|\lambda|<2\omega^{2}$. In terms of the
relative position $\displaystyle{x}={x_{1}}-{x_{2}}$, the relative momentum
$\displaystyle{p}=\tfrac{1}{2}({p_{1}}-{p_{2}})$, the COM coordinate
$\displaystyle{X}=\tfrac{1}{2}({x_{1}}+{x_{2}})$ and the total momentum
$\displaystyle{P}={p_{1}}+{p_{2}}$, $H$ has the separable form

$\displaystyle
H=\frac{{P}^{2}}{4m}+\frac{{p}^{2}}{m}+m\omega_{c}^{2}{X}^{2}+\frac{m}{4}\omega_{r}^{2}{x}^{2},$
(5)

where $\omega_{c}^{2}=\omega^{2}+\frac{1}{2}\lambda$ and
$\omega_{r}^{2}=\omega^{2}-\frac{1}{2}\lambda$. It is convenient to work in a
basis set of direct products
$\ket{n_{c},n_{r}}\equiv\ket{n_{c}}\otimes\ket{n_{r}}$ of the energy eigenkets
of the COM and the relative system, where $n_{c}$, $n_{r}$ $=0,1,2,\dotsc$ .
In terms of the raising operators

$a_{r}^{\dagger}=\Big{(}\frac{m\omega_{r}}{4\hbar}\Big{)}^{1/2}\Big{(}{x}-\frac{2ip}{m\omega_{r}}\Big{)},\;\;a_{c}^{\dagger}=\Big{(}\frac{m\omega_{c}}{\hbar}\Big{)}^{1/2}\Big{(}{X}-\frac{iP}{2m\omega_{c}}\Big{)}$
(6)

and the corresponding lowering operators, the Hamiltonian becomes

$H=\hbar\omega_{r}\big{(}a_{r}^{\dagger}a_{r}+\tfrac{1}{2}\big{)}+\hbar\omega_{c}\big{(}a_{c}^{\dagger}a_{c}+\tfrac{1}{2}\big{)},$
(7)

whereas in terms of the ladder operators $a_{i}$ and $a_{i}^{\dagger}$
$(i=1,2)$ corresponding to the original pair of coupled oscillators, the
Hamiltonian is

$H=\hbar\omega(a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2}+1)+\frac{\hbar\lambda}{\omega}(a_{1}^{\dagger}a_{2}^{\dagger}+a_{1}^{\dagger}a_{2}+a_{1}a_{2}^{\dagger}+a_{1}a_{2}).$
(8)

###  III.1 Uncoupled oscillators

In the special case $\lambda=0$, the oscillators $1$ and $2$ are uncoupled and
$\omega_{c}=\omega_{r}=\omega$. The eigenstates of the Hamiltonian can then be
written in terms of the factored products of the energy eigenbasis of the two
oscillators. For notational clarity, we will denote $\ket{n_{1},n_{2}}$ by
$\ket{\nu_{1},\nu_{2}}$ and $\ket{n_{c},n_{r}}$ by $\ket{\nu_{c},\nu_{r}}$ in
the case $\lambda=0$. Further, we have in this case
$\nu_{1}+\nu_{2}=\nu_{c}+\nu_{r}$, and
$a_{c}^{\dagger}=(a_{1}^{\dagger}+a_{2}^{\dagger})/\sqrt{2}$,
$a_{r}^{\dagger}=(a_{1}^{\dagger}-a_{2}^{\dagger})/\sqrt{2}$.

We now consider the bipartite density matrix

$\displaystyle\rho$ $\displaystyle=\ket{\nu_{c},\nu_{r}}\bra{\nu_{c},\nu_{r}}$
$\displaystyle=\sum_{\nu_{1},\nu_{2}}\sum_{\nu_{1}^{\prime},\nu_{2}^{\prime}}\braket{\nu_{1},\nu_{2}}{\nu_{c},\nu_{r}}\braket{\nu_{c},\nu_{r}}{\nu_{1}^{\prime},\nu_{2}^{\prime}}\ket{\nu_{1},\nu_{2}}\bra{\nu_{1}^{\prime},\nu_{2}^{\prime}},$
(9)

where

$\braket{\nu_{1},\nu_{2}}{\nu_{c},\nu_{r}}=\frac{1}{(2^{\nu_{c}+\nu_{r}}\nu_{c}!\nu_{r}!)^{1/2}}\bra{\nu_{1},\nu_{2}}{}(a_{1}^{\dagger}+a_{2}^{\dagger})^{\nu_{c}}(a_{1}^{\dagger}-a_{2}^{\dagger})^{\nu_{r}}\ket{\nu_{c}=0,\nu_{r}=0}.$
(10)

In explicit form,

$\displaystyle\braket{\nu_{1},\nu_{2}}{\nu_{c},\nu_{r}}=$
$\displaystyle\sum_{k=0}^{\nu_{r}}\sum_{k^{\prime}=0}^{\nu_{c}}\Big{\\{}\frac{(-1)^{k^{\prime}}[\nu_{c}!\,\nu_{r}!\,(\nu_{r}+\nu_{c}-k-k^{\prime})!\,(k+k^{\prime})!]^{1/2}}{2^{(\nu_{c}+\nu_{r})/2}(\nu_{r}-k)!\,k!\,(\nu_{c}-k^{\prime})!\,k^{\prime}!}\,\times$
$\displaystyle\delta_{\nu_{1},\nu_{r}+\nu_{c}-k-k^{\prime}}\,\delta_{\nu_{2},k+k^{\prime}}\Big{\\}}.$
(11)

Tracing out oscillator 2 in Eq. (III.1), the reduced density matrix becomes

$\displaystyle\widetilde{\rho}$
$\displaystyle=\sum_{\nu_{1},\nu_{1}^{\prime}}\braket{\nu_{1},\nu_{c}+\nu_{r}-\nu_{1}}{\nu_{c},\nu_{r}}\delta_{\nu_{1},\nu_{1}^{\prime}}\braket{\nu_{c},\nu_{r}}{\nu_{1}^{\prime},\nu_{c}+\nu_{r}-\nu_{1}^{\prime}}\ket{\nu_{1}}\bra{\nu_{1}^{\prime}}$
$\displaystyle=\sum_{\nu_{1}}|\braket{\nu_{1},\nu_{c}+\nu_{r}-\nu_{1}}{\nu_{c},\nu_{r}}|^{2}\ket{\nu_{1}}\bra{\nu_{1}}.$
(12)

Since the Schmidt coefficients
$|\braket{\nu_{1},\nu_{c}+\nu_{r}-\nu_{1}}{\nu_{c},\nu_{r}}|^{2}$ are
available explicitly from Eq.(III.1), the SLE and SVNE can be obtained
analytically when $\lambda=0$. These are shown for the first few values of
$\nu_{r}$, keeping $\nu_{c}=0$, in Table I.

In what follows we examine, in the general case $\lambda\neq 0$, both the
standard entanglement measures as well as the TEIs, and compare the two sets
of quantities.

###  III.2 Coupled oscillators: SLE and SVNE

We analyze two cases in detail, with the COM oscillator in its ground state
($n_{c}=0$) in both instances: Case (a) $n_{r}=0$, and Case (b)
$n_{r}\geqslant 1$. We first compute the SLE and SVNE in Case (a), for which
both can be found analytically. Noting that the COM and the relative
coordinate wave functions are Gaussians in $X$ and $x$ respectively, the
matrix elements of the density matrix in the position basis of the oscillators
1 and 2 work out to

$\braket{x_{1},x_{2}}{\rho}{x_{1}^{\prime},x_{2}^{\prime}}=(4/\pi)\,\sqrt{c_{1}c_{2}}\,e^{-c_{1}[({x_{1}+x_{2}})^{2}+(x_{1}^{\prime}+x_{2}^{\prime})^{2}]}e^{-c_{2}[(x_{1}-x_{2})^{2}+(x_{1}^{\prime}-x_{2}^{\prime})^{2}]},$
(13)

where $c_{1}=m\omega_{c}/4\hbar$ and $c_{2}=m\omega_{r}/4\hbar$. Tracing over
oscillator $2$, the reduced density matrix element becomes

$\braket{x_{1}}{\widetilde{\rho}}{x_{1}^{\prime}}=\Big{[}\frac{8c_{1}c_{2}}{\pi(c_{1}+c_{2})}\Big{]}^{1/2}\,\mathrm{exp}\,\Big{\\{}-(c_{1}+c_{2})(x_{1}^{2}+x_{1}^{\prime
2})+\frac{(c_{1}-c_{2})^{2}}{2(c_{1}+c_{2})}(x_{1}+x_{1}^{\prime})^{2}\Big{\\}}.$
(14)

The SLE can now be computed in a straightforward manner. In terms of the ratio
$\eta=\omega_{c}/\omega_{r}$, it is given by

$1-\mathrm{Tr}\,\widetilde{\rho}^{\,2}=\frac{(1-\sqrt{\eta})^{2}}{1+\eta}.$
(15)

As expected, the SLE vanishes when $\eta=1$. We also note that it is unchanged
under $\eta\to 1/\eta$, i.e., the interchange
$\omega_{c}\leftrightarrow\omega_{r}$ (See Fig.1(a)).

A similar procedure cannot be used to compute the SVNE because the matrix
elements of $\mathrm{ln}\ \widetilde{\rho}$ in the position basis are not
readily calculable analytically. However, for the case $n_{c}=n_{r}=0$ alone
(with which we are presently concerned), an analytic expression for the
Schmidt coefficients in terms of the system parameters has been obtained
starting from a trial ansatz [10]. Further, a mapping between
$\widetilde{\rho}$ and the density matrix of a single oscillator in a heat
bath has also been mentioned. We indicate below the computation of the SVNE
using such a mapping. The reduced density matrix of a single oscillator of
frequency $\varpi$ in a heat bath at inverse temperature $\beta$ is given by
$\widetilde{\rho}=e^{-\beta{\mathcal{H}}}/{\mathrm{Tr}\,e^{-\beta{\mathcal{H}}}}$,
where $\mathcal{H}$ is its Hamiltonian. The eigenvalues of $\widetilde{\rho}$
are given by

$\displaystyle\lambda_{n}=(1-e^{-\beta\hbar\varpi})\,e^{-n\beta\hbar\varpi},\quad
n=0,1,2,\dotsc$ (16)

It follows at once that the SLE is

$1-\mathrm{Tr}\,\widetilde{\rho}^{\,2}=1-\mathrm{tanh}\,\tfrac{1}{2}\beta\hbar\varpi,$
(17)

while the SVNE is

$-\mathrm{Tr}\,(\widetilde{\rho}\,\,\mathrm{ln}\,\widetilde{\rho})=-\mathrm{ln}\,(1-e^{-\beta\hbar\varpi})+\frac{\beta\hbar\varpi}{e^{\beta\hbar\varpi}-1}.$
(18)

The matrix elements of $\widetilde{\rho}$ in the position basis are given by

$\braket{x_{1}}{\widetilde{\rho}}{x_{1}^{\prime}}=\sum_{n=0}^{\infty}\lambda_{n}\,\chi_{n}(x_{1})\,\chi_{n}(x_{1}^{\prime}),$
(19)

where $\chi_{n}$ is the oscillator wave function given by

$\chi_{n}(x_{1})=\frac{(m\varpi/\pi\hbar)^{1/4}}{(2^{n}n!)^{1/2}}H_{n}(x_{1}\sqrt{m\varpi/\hbar})\,e^{-m\varpi
x_{1}^{2}/2\hbar}$ (20)

and $H_{n}$ is the Hermite polynomial of order n. The sum in Eq. (19) is
evaluated using Mehler’s formula

$\sum_{n=0}^{\infty}H_{n}(u)H_{n}(v)\frac{s^{n}}{2^{n}n!}=(1-s^{2})^{-1/2}\,\mathrm{exp}\,\Big{\\{}-\frac{s^{2}u^{2}+s^{2}v^{2}-2uvs}{1-s^{2}}\Big{\\}},\
|s|<1.$ (21)

The result is

$\displaystyle\braket{x_{1}}{\widetilde{\rho}}{x_{1}^{\prime}}=\Big{[}\frac{m\varpi\,\mathrm{tanh}\,\tfrac{1}{2}\beta\hbar\varpi}{\pi\hbar}\Big{]}^{1/2}\mathrm{exp}\,\Big{[}-\frac{m\varpi}{2\hbar}$
$\displaystyle\Big{\\{}(x_{1}^{2}+x_{1}^{\prime
2})\,\mathrm{coth}\,\beta\hbar\varpi$
$\displaystyle-2\,x_{1}\,x_{1}^{\prime}\,\mathrm{cosech}\,\beta\hbar\varpi\Big{\\}}\Big{]}.$
(22)

The expressions in Eqs. (III.2) and (14) are the same if we make the
identifications

$\varpi=(\omega_{c}\,\omega_{r})^{1/2}$ (23)

for the frequency of the oscillator and

$\beta=\frac{1}{\hbar\varpi}\,\displaystyle{\mathrm{ln}\,\Big{(}\frac{1+\sqrt{\eta}}{1-\sqrt{\eta}}\Big{)}^{2}}$
(24)

for the inverse temperature of the heat bath. Since we already have the
expressions for the SLE and SVNE for this system, we can now readily identify
the corresponding entanglement measures for the original problem of the two
coupled oscillators. Equation (15) for the SLE is immediately corroborated.
For the SVNE, inserting the correspondences of Eqs. (23) and (24) in Eq. (18)
we get

$-\mathrm{Tr}\,(\widetilde{\rho}\,\,\mathrm{ln}\,\widetilde{\rho})=-\mathrm{ln}\,\Big{(}\frac{4\sqrt{\eta}}{(1+\sqrt{\eta})^{2}}\Big{)}+\frac{(1-\sqrt{\eta})^{2}}{4\sqrt{\eta}}\mathrm{ln}\,\Big{(}\frac{1+\sqrt{\eta}}{1-\sqrt{\eta}}\Big{)}^{2}.$
(25)

The SVNE diverges as $\eta\to 0$ (or $\omega_{c}\to 0$), since the COM
oscillator is delocalized [8] in this case. It vanishes when $\eta=1$, as
expected, and is symmetric (like the SLE, Eq. (15)) under the interchange
$\eta\leftrightarrow 1/\eta$. This symmetry is manifest in the log-log plot of
the SVNE versus $\eta$ as shown in Figs. 1(a) and 1(b).

We now proceed to discuss Case (b), i.e., $n_{c}=0$, $n_{r}\geqslant 1$. The
Schmidt coefficients, and hence the extent of entanglement between the coupled
oscillators 1 and 2, have been obtained as functions of $\eta$ by numerically
diagonalizing the reduced density matrix. Plots of the SVNE and SLE versus
$\eta$, for a set of values of $n_{r}$, are shown in Figs. 2(a) and 2(b). It
is clear that both these entanglement measures increase with $n_{r}$, or
equivalently, with the energy of the oscillator in the relative coordinate.

###  III.3 Coupled oscillators: Tomographic entanglement indicators

We now calculate the TEIs (defined in Eqs. (1)–(3)) in our oscillator model.
Once again, in the limit $\eta\rightarrow 0$ it can be shown that both
$\mathrm{\epsilon_{{}_{BD}}}$ and $\mathrm{\epsilon_{{}_{KL}}}$ diverge as a
consequence of the delocalization of the COM oscillator. Setting $n_{c}=0$ as
before, plots of $\mathrm{\epsilon_{{}_{BD}}}$ versus $\eta$ for a range of
values of $n_{r}$ have been obtained numerically and displayed in Fig. 3. We
see that the minimum in $\mathrm{\epsilon_{{}_{BD}}}$ as a function of $\eta$
in Fig. 3a drifts to smaller values of $\eta$ with increasing $n_{r}$. This is
an artefact arising from our choice of the position slice of the tomogram.
Figure 3b corresponds to the momentum slice of the same tomogram, in which the
minimum drifts to larger values of $\eta$ with increasing $n_{r}$. Figure 3c
depicts the average of the two slices, showing that the minimum value of
$\mathrm{\epsilon_{{}_{BD}}}$ remains at $\eta=1$ independent of $n_{r}$, in
accordance with the behavior of the SLE and SVNE. The same property is
exhibited by the averaged TEI $\mathrm{\epsilon_{{}_{KL}}}$, as seen in Fig.
4. Overall, the TEIs $\mathrm{\epsilon_{{}_{BD}}}$ and
$\mathrm{\epsilon_{{}_{KL}}}$ do reflect quite faithfully the variations in
the SLE and SVNE.

We now proceed to compute $\mathrm{\epsilon_{{}_{IPR}}}$ (Eq. (3)). We
introduce dimensionless variables $\overline{x}=x/L_{r}$ and
$\overline{X}=X/L_{c}$ where $L_{r}=(2\hbar/m\omega_{r})^{1/2}$ and
$L_{c}=(\hbar/2m\omega_{c})^{1/2}$. Defining
$\overline{x}_{1,2}=\overline{X}\pm\tfrac{1}{2}\overline{x}$ we have, in terms
of the original oscillator coordinates $x_{1}$ and $x_{2}$,

$\overline{x}_{1}=\tfrac{1}{2}(L_{c}^{-1}+L_{r}^{-1})x_{1}+\tfrac{1}{2}(L_{c}^{-1}-L_{r}^{-1})x_{2}$
(26)

and

$\overline{x}_{2}=\tfrac{1}{2}(L_{c}^{-1}-L_{r}^{-1})x_{1}+\tfrac{1}{2}(L_{c}^{-1}+L_{r}^{-1})x_{2}.$
(27)

The further requirement that $\overline{x}_{1}$ and $\overline{x}_{2}$ be
proportional to $x_{1}$ and $x_{2}$, respectively, imposes the condition
$L_{c}=L_{r}$ ($=L$, say). This restricts $\eta$ to the value $\frac{1}{4}$.
The dimensionless density matrix elements are given by

$\displaystyle\bra{\overline{x}_{1},\overline{x}_{2}}{\overline{\rho}}\ket{\overline{x}_{1}^{\prime},\overline{x}_{2}^{\prime}}=L^{2}\bra{Lx_{1},Lx_{2}}{\rho}\ket{Lx_{1}^{\prime},Lx_{2}^{\prime}}$
(28)

and the corresponding TEI is given by

$\displaystyle\mathrm{\epsilon_{{}_{IPR}}}=1+\int\\!\\!d\overline{x}_{1}\\!\\!\int\\!\\!d\overline{x}_{2}\bra{\overline{x}_{1},\overline{x}_{2}}{\overline{\rho}}\ket{\overline{x}_{1},\overline{x}_{2}}^{2}-\int\\!\\!d\overline{x}_{1}\bra{\overline{x}_{1}}{\overline{\widetilde{\rho}}}_{1}\ket{\overline{x}_{1}}^{2}-\int\\!\\!d\overline{x}_{2}\bra{\overline{x}_{2}}{\overline{\widetilde{\rho}}}_{2}\ket{\overline{x}_{2}}^{2}.$
(29)

A numerically generated plot of $\mathrm{\epsilon_{{}_{IPR}}}$ versus $n_{r}$,
in the position slice, for $\eta=\frac{1}{4}$ and $n_{c}=0$, is shown in Fig.
5. It is again evident that $\mathrm{\epsilon_{{}_{IPR}}}$ increases with
increasing $n_{r}$, following the same general trend as all the other
indicators considered.

##  IV Concluding remarks

We have investigated the entanglement between two coupled harmonic oscillators
with the same natural frequency. The system is separable in the center of mass
(COM) and relative coordinates. The standard entanglement measures (SLE and
SVNE) have been found analytically when the oscillators in the COM and
relative coordinates are in the ground state ($n_{c}=n_{r}=0$), and
numerically for $n_{c}=0,n_{r}\geqslant 1$. In the special case when the
oscillators in the COM and relative coordinates have the same frequency,
Schmidt diagonalization has been carried out to yield the SLE and SVNE
analytically. Tomographic entanglement indicators (the Bhattacharya distance,
the Kullback-Liebler divergence, and the IPR based indicator) have been
computed and compared with the SLE and SVNE. As with the standard entanglement
measures, the TEIs also indicate that localization is crucial to ensure
meaningful estimates of entanglement. (When the COM is not localized, all
indicators reach their maximum values for _all_ eigenstates.) We have also
shown that the entanglement increases with the energy of the oscillator in the
relative coordinate.

It would be desirable to extend the preliminary understanding of the effects
of localization and binding on entanglement gained by our study of this toy
model to many-body CV systems in both atomic and nuclear physics. The success
of TEIs in mimicking the entanglement measures could be utilized in more
complicated systems where the computation of conventional entanglement
measures is not viable, and state reconstruction from tomograms is extremely
challenging.

## Acknowledgements

We are pleased to dedicate this work to Prof. A. Ravi P. Rau with whom we have
had many fruitful discussions over the years. We acknowledge partial support
through funds from Mphasis to the Center for Quantum Information, Computing
and Communication (CQuICC). SL thanks the Department of Physics, Indian
Institute of Technology Madras for infrastructural support.

## References

  * [1] T. Yu and J. Eberly, Science 323, 598 (2009).
  * [2] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * [3] B. Sharmila, S. Lakshmibala, and V. Balakrishnan, J. Phys. B 53, 245502 (2020).
  * [4] J. Karthik, A. Sharma, and A. Lakshminarayan, Phys. Rev. A 75, 022304 (2007).
  * [5] S. Paul, S. Lakshmibala, V. Balakrishnan, and S. Ramanan, Front. Quantum Sci. Technol. 2, 1131798 (2023).
  * [6] K. Mallayya, R. Tibrewala, S. Shankaranarayanan, and T. Padmanabhan, Phys. Rev. D 90, 044058 (2014).
  * [7] S. M. Chandran and S. Shankaranarayanan, Phys. Rev. D 99, 045010 (2019).
  * [8] S. Qvarfort, S. Bose, and A. Serafini, New J. Phys. 22, 093062 (2020).
  * [9] P. Tommasini, E. Timmermans, and A. F. R. de Toledo Piza, Am. J. Phys. 66, 881 (1998).
  * [10] M. Srednicki, Phys. Rev. Lett. 71, 666 (1993).
  * [11] D. N. Makarov, Phys. Rev. E 97, 042203 (2018).
  * [12] M. Man’Ko, V. Man’Ko, G. Marmo, A. Simoni, and F. Ventriglia, Nuovo Cim. C 36, 163 (2013).
  * [13] T. Kailath, IEEE Trans. Commun. Technol. 15, 52 (1967).
  * [14] L. Viola and W. G. Brown, J. Phys. A 40, 8109 (2007).

$\nu_{c}$ | $\nu_{r}$ | SLE | SVNE
---|---|---|---
0 | 0 | 0 | 0
0 | 1 | $\frac{1}{2}$ | $\mathrm{ln}\ 2$
0 | 2 | $\frac{5}{8}$ | $\frac{3}{2}\ \mathrm{ln}\ 2$
0 | 3 | $\frac{11}{16}$ | $\mathrm{ln}\ 8-\frac{3}{4}\ \mathrm{ln}\ 3$
0 | 4 | $\frac{93}{128}$ | $\frac{1}{8}(21\ \mathrm{ln}\ 2-3\ \mathrm{ln}\ 3)$
0 | 5 | $\frac{193}{256}$ | $\frac{1}{16}(70\ \mathrm{ln}\ 2-15\ \mathrm{ln}\ 5)$

Table 1: SLE and SVNE for different values of $\nu_{r}$ and $\nu_{c}$ computed
using the Schmidt coefficients from Eq.(III.1)

(a)

(b)

Figure 1: Log-log plots of the SLE (Eq.(15)) and SVNE (Eq.(25)) versus $\eta$
in the case $n_{c}=0$, $n_{r}=0$.

(a)

(b)

Figure 2: Log-log plots of the (a) SLE and (b) SVNE for $n_{c}=0$,
$n_{r}\geqslant 1$. Figure 3: Log-log plot of $\mathrm{\epsilon_{{}_{BD}}}$
versus $\eta$ for $n_{c}=0$, $n_{r}\geqslant 1$. Top left and right, position
and momentum slices respectively; bottom, average of position and momentum
slices. Figure 4: Log-log plot of $\mathrm{\epsilon_{{}_{KL}}}$ averaged over
the position and momentum slices versus $\eta$, for $n_{c}=0$ and
$n_{r}\geqslant 1$. Figure 5: $\mathrm{\epsilon_{{}_{IPR}}}$ (in the position
slice) versus $n_{r}$, with $n_{c}=0$ and $\eta=\tfrac{1}{4}$.