Also at ]INFN-Roma Tre, 00146 Roma, ItalyNow at ]Department of Physics &
Astronomy, Bucknell University, Lewisburg, PA, USAAlso at ]Coimbra Polytechnic
- ISEC, 3030-199 Coimbra, PortugalAlso at ]Physikalisches Institut,
Universität Heidelberg, Heidelberg, GermanyXENON Collaboration

# XENONnT WIMP Search: Signal & Background Modeling and Statistical Inference

E. Aprile  0000-0001-6595-7098 Physics Department, Columbia University, New
York, NY 10027, USA  J. Aalbers  0000-0003-0030-0030 Nikhef and the University
of Groningen, Van Swinderen Institute, 9747AG Groningen, Netherlands  K. Abe
0009-0000-9620-788X Kamioka Observatory, Institute for Cosmic Ray Research,
and Kavli Institute for the Physics and Mathematics of the Universe (WPI),
University of Tokyo, Higashi-Mozumi, Kamioka, Hida, Gifu 506-1205, Japan  S.
Ahmed Maouloud  0000-0002-0844-4576 LPNHE, Sorbonne Université, CNRS/IN2P3,
75005 Paris, France  L. Althueser  0000-0002-5468-4298 Institut für
Kernphysik, Westfälische Wilhelms-Universität Münster, 48149 Münster, Germany
B. Andrieu  0009-0002-6485-4163 LPNHE, Sorbonne Université, CNRS/IN2P3, 75005
Paris, France  E. Angelino  0000-0002-6695-4355 INAF-Astrophysical Observatory
of Torino, Department of Physics, University of Torino and INFN-Torino, 10125
Torino, Italy  INFN-Laboratori Nazionali del Gran Sasso and Gran Sasso Science
Institute, 67100 L’Aquila, Italy  D. Antón Martin  0000-0001-7725-5552
Department of Physics & Kavli Institute for Cosmological Physics, University
of Chicago, Chicago, IL 60637, USA  F. Arneodo  0000-0002-1061-0510 New York
University Abu Dhabi - Center for Astro, Particle and Planetary Physics, Abu
Dhabi, United Arab Emirates  L. Baudis  0000-0003-4710-1768 Physik-Institut,
University of Zürich, 8057 Zürich, Switzerland  M. Bazyk  0009-0000-7986-153X
SUBATECH, IMT Atlantique, CNRS/IN2P3, Université de Nantes, Nantes 44307,
France  L. Bellagamba  0000-0001-7098-9393 Department of Physics and
Astronomy, University of Bologna and INFN-Bologna, 40126 Bologna, Italy  R.
Biondi  0000-0002-6622-8740 Max-Planck-Institut für Kernphysik, 69117
Heidelberg, Germany  A. Bismark  0000-0002-0574-4303 Physik-Institut,
University of Zürich, 8057 Zürich, Switzerland  K. Boese  0009-0007-0662-0920
Max-Planck-Institut für Kernphysik, 69117 Heidelberg, Germany  A. Brown
0000-0002-1623-8086 Physikalisches Institut, Universität Freiburg, 79104
Freiburg, Germany  G. Bruno  0000-0001-9005-2821 SUBATECH, IMT Atlantique,
CNRS/IN2P3, Université de Nantes, Nantes 44307, France  R. Budnik
0000-0002-1963-9408 Department of Particle Physics and Astrophysics, Weizmann
Institute of Science, Rehovot 7610001, Israel  J. M. R. Cardoso
0000-0002-8832-8208 LIBPhys, Department of Physics, University of Coimbra,
3004-516 Coimbra, Portugal  A. P. Cimental Chávez  0009-0004-9605-5985 Physik-
Institut, University of Zürich, 8057 Zürich, Switzerland  A. P. Colijn
0000-0002-3118-5197 Nikhef and the University of Amsterdam, Science Park,
1098XG Amsterdam, Netherlands  J. Conrad  0000-0001-9984-4411 Oskar Klein
Centre, Department of Physics, Stockholm University, AlbaNova, Stockholm
SE-10691, Sweden  J. J. Cuenca-García  0000-0002-3869-7398 Physik-Institut,
University of Zürich, 8057 Zürich, Switzerland  V. D’Andrea
0000-0003-2037-4133 [  INFN-Laboratori Nazionali del Gran Sasso and Gran Sasso
Science Institute, 67100 L’Aquila, Italy  L. C. Daniel Garcia
0009-0000-5813-9118 LPNHE, Sorbonne Université, CNRS/IN2P3, 75005 Paris,
France  M. P. Decowski  0000-0002-1577-6229 Nikhef and the University of
Amsterdam, Science Park, 1098XG Amsterdam, Netherlands  C. Di Donato
0009-0005-9268-6402 Department of Physics and Chemistry, University of
L’Aquila, 67100 L’Aquila, Italy  P. Di Gangi  0000-0003-4982-3748 Department
of Physics and Astronomy, University of Bologna and INFN-Bologna, 40126
Bologna, Italy  S. Diglio  0000-0002-9340-0534 SUBATECH, IMT Atlantique,
CNRS/IN2P3, Université de Nantes, Nantes 44307, France  K. Eitel
0000-0001-5900-0599 Institute for Astroparticle Physics, Karlsruhe Institute
of Technology, 76021 Karlsruhe, Germany  A. Elykov  0000-0002-2693-232X
Institute for Astroparticle Physics, Karlsruhe Institute of Technology, 76021
Karlsruhe, Germany  A. D. Ferella  0000-0002-6006-9160 Department of Physics
and Chemistry, University of L’Aquila, 67100 L’Aquila, Italy  INFN-Laboratori
Nazionali del Gran Sasso and Gran Sasso Science Institute, 67100 L’Aquila,
Italy  C. Ferrari  0000-0002-0838-2328 INFN-Laboratori Nazionali del Gran
Sasso and Gran Sasso Science Institute, 67100 L’Aquila, Italy  H. Fischer
0000-0002-9342-7665 Physikalisches Institut, Universität Freiburg, 79104
Freiburg, Germany  T. Flehmke  0009-0002-7944-2671 Oskar Klein Centre,
Department of Physics, Stockholm University, AlbaNova, Stockholm SE-10691,
Sweden  M. Flierman  0000-0002-3785-7871 Nikhef and the University of
Amsterdam, Science Park, 1098XG Amsterdam, Netherlands  W. Fulgione
0000-0002-2388-3809 INAF-Astrophysical Observatory of Torino, Department of
Physics, University of Torino and INFN-Torino, 10125 Torino, Italy  INFN-
Laboratori Nazionali del Gran Sasso and Gran Sasso Science Institute, 67100
L’Aquila, Italy  C. Fuselli  0000-0002-7517-8618 Nikhef and the University of
Amsterdam, Science Park, 1098XG Amsterdam, Netherlands  P. Gaemers
0009-0003-1108-1619 Nikhef and the University of Amsterdam, Science Park,
1098XG Amsterdam, Netherlands  R. Gaior  0009-0005-2488-5856 LPNHE, Sorbonne
Université, CNRS/IN2P3, 75005 Paris, France  M. Galloway  0000-0002-8323-9564
Physik-Institut, University of Zürich, 8057 Zürich, Switzerland  F. Gao
0000-0003-1376-677X Department of Physics & Center for High Energy Physics,
Tsinghua University, Beijing 100084, P.R. China  S. Ghosh  0000-0001-7785-9102
Department of Physics and Astronomy, Purdue University, West Lafayette, IN
47907, USA  R. Giacomobono  0000-0001-6162-1319 Department of Physics “Ettore
Pancini”, University of Napoli and INFN-Napoli, 80126 Napoli, Italy  R. Glade-
Beucke  0009-0006-5455-2232 Physikalisches Institut, Universität Freiburg,
79104 Freiburg, Germany  L. Grandi  0000-0003-0771-7568 Department of Physics
& Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL
60637, USA  J. Grigat  0009-0005-4775-0196 Physikalisches Institut,
Universität Freiburg, 79104 Freiburg, Germany  H. Guan  0009-0006-5049-0812
Department of Physics and Astronomy, Purdue University, West Lafayette, IN
47907, USA  M. Guida  0000-0001-5126-0337 Max-Planck-Institut für Kernphysik,
69117 Heidelberg, Germany  P. Gyoergy  0009-0005-7616-5762 Institut für Physik
& Exzellenzcluster PRISMA+, Johannes Gutenberg-Universität Mainz, 55099 Mainz,
Germany  R. Hammann  0000-0001-6149-9413<EMAIL_ADDRESS>Max-
Planck-Institut für Kernphysik, 69117 Heidelberg, Germany  A. Higuera
0000-0001-9310-2994 Department of Physics and Astronomy, Rice University,
Houston, TX 77005, USA  C. Hils  Institut für Physik & Exzellenzcluster
PRISMA+, Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany  L.
Hoetzsch  0000-0003-2572-477X<EMAIL_ADDRESS>Max-Planck-Institut
für Kernphysik, 69117 Heidelberg, Germany  N. F. Hood  0000-0003-2507-7656
Department of Physics, University of California San Diego, La Jolla, CA 92093,
USA  M. Iacovacci  0000-0002-3102-4721 Department of Physics “Ettore Pancini”,
University of Napoli and INFN-Napoli, 80126 Napoli, Italy  Y. Itow
0000-0002-8198-1968 Kobayashi-Maskawa Institute for the Origin of Particles
and the Universe, and Institute for Space-Earth Environmental Research, Nagoya
University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602, Japan  J. Jakob
0009-0000-2220-1418 Institut für Kernphysik, Westfälische Wilhelms-Universität
Münster, 48149 Münster, Germany  F. Joerg  0000-0003-1719-3294 Max-Planck-
Institut für Kernphysik, 69117 Heidelberg, Germany  Y. Kaminaga
0009-0006-5424-2867 Kamioka Observatory, Institute for Cosmic Ray Research,
and Kavli Institute for the Physics and Mathematics of the Universe (WPI),
University of Tokyo, Higashi-Mozumi, Kamioka, Hida, Gifu 506-1205, Japan  M.
Kara  0009-0004-5080-9446 Institute for Astroparticle Physics, Karlsruhe
Institute of Technology, 76021 Karlsruhe, Germany  P. Kavrigin
0009-0000-1339-2419 Department of Particle Physics and Astrophysics, Weizmann
Institute of Science, Rehovot 7610001, Israel  S. Kazama  0000-0002-6976-3693
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, and
Institute for Space-Earth Environmental Research, Nagoya University, Furo-cho,
Chikusa-ku, Nagoya, Aichi 464-8602, Japan  M. Kobayashi  0009-0006-7861-1284
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, and
Institute for Space-Earth Environmental Research, Nagoya University, Furo-cho,
Chikusa-ku, Nagoya, Aichi 464-8602, Japan  A. Kopec  0000-0001-6548-0963 [
Department of Physics, University of California San Diego, La Jolla, CA 92093,
USA  F. Kuger  0000-0001-9475-3916 Physikalisches Institut, Universität
Freiburg, 79104 Freiburg, Germany  H. Landsman  0000-0002-7570-5238 Department
of Particle Physics and Astrophysics, Weizmann Institute of Science, Rehovot
7610001, Israel  R. F. Lang  0000-0001-7594-2746 Department of Physics and
Astronomy, Purdue University, West Lafayette, IN 47907, USA  L. Levinson
0000-0003-4679-0485 Department of Particle Physics and Astrophysics, Weizmann
Institute of Science, Rehovot 7610001, Israel  I. Li  0000-0001-6655-3685
Department of Physics and Astronomy, Rice University, Houston, TX 77005, USA
S. Li  0000-0003-0379-1111 Department of Physics, School of Science, Westlake
University, Hangzhou 310030, P.R. China  S. Liang  0000-0003-0116-654X
Department of Physics and Astronomy, Rice University, Houston, TX 77005, USA
Y.-T. Lin  0000-0003-3631-1655 Max-Planck-Institut für Kernphysik, 69117
Heidelberg, Germany  S. Lindemann  0000-0002-4501-7231 Physikalisches
Institut, Universität Freiburg, 79104 Freiburg, Germany  M. Lindner
0000-0002-3704-6016 Max-Planck-Institut für Kernphysik, 69117 Heidelberg,
Germany  K. Liu  0009-0004-1437-5716 Department of Physics & Center for High
Energy Physics, Tsinghua University, Beijing 100084, P.R. China  J. Loizeau
0000-0001-6375-9768 SUBATECH, IMT Atlantique, CNRS/IN2P3, Université de
Nantes, Nantes 44307, France  F. Lombardi  0000-0003-0229-4391 Institut für
Physik & Exzellenzcluster PRISMA+, Johannes Gutenberg-Universität Mainz, 55099
Mainz, Germany  J. Long  0000-0002-5617-7337 Department of Physics & Kavli
Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637,
USA  J. A. M. Lopes  0000-0002-6366-2963 [  LIBPhys, Department of Physics,
University of Coimbra, 3004-516 Coimbra, Portugal  T. Luce
8561-4854-7251-585X Physikalisches Institut, Universität Freiburg, 79104
Freiburg, Germany  Y. Ma  0000-0002-5227-675X Department of Physics,
University of California San Diego, La Jolla, CA 92093, USA  C. Macolino
0000-0003-2517-6574 Department of Physics and Chemistry, University of
L’Aquila, 67100 L’Aquila, Italy  INFN-Laboratori Nazionali del Gran Sasso and
Gran Sasso Science Institute, 67100 L’Aquila, Italy  J. Mahlstedt
0000-0002-8514-2037 Oskar Klein Centre, Department of Physics, Stockholm
University, AlbaNova, Stockholm SE-10691, Sweden  A. Mancuso
0009-0002-2018-6095 Department of Physics and Astronomy, University of Bologna
and INFN-Bologna, 40126 Bologna, Italy  L. Manenti  0000-0001-7590-0175 New
York University Abu Dhabi - Center for Astro, Particle and Planetary Physics,
Abu Dhabi, United Arab Emirates  F. Marignetti  0000-0001-8776-4561 Department
of Physics “Ettore Pancini”, University of Napoli and INFN-Napoli, 80126
Napoli, Italy  T. Marrodán Undagoitia  0000-0001-9332-6074 Max-Planck-Institut
für Kernphysik, 69117 Heidelberg, Germany  K. Martens  0000-0002-5049-3339
Kamioka Observatory, Institute for Cosmic Ray Research, and Kavli Institute
for the Physics and Mathematics of the Universe (WPI), University of Tokyo,
Higashi-Mozumi, Kamioka, Hida, Gifu 506-1205, Japan  J. Masbou
0000-0001-8089-8639 SUBATECH, IMT Atlantique, CNRS/IN2P3, Université de
Nantes, Nantes 44307, France  E. Masson  0000-0002-5628-8926 LPNHE, Sorbonne
Université, CNRS/IN2P3, 75005 Paris, France  S. Mastroianni
0000-0002-9467-0851 Department of Physics “Ettore Pancini”, University of
Napoli and INFN-Napoli, 80126 Napoli, Italy  A. Melchiorre
0009-0006-0615-0204 Department of Physics and Chemistry, University of
L’Aquila, 67100 L’Aquila, Italy  M. Messina  0000-0002-6475-7649 INFN-
Laboratori Nazionali del Gran Sasso and Gran Sasso Science Institute, 67100
L’Aquila, Italy  A. Michael  Institut für Kernphysik, Westfälische Wilhelms-
Universität Münster, 48149 Münster, Germany  K. Miuchi  0000-0002-1546-7370
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan  A.
Molinario  0000-0002-5379-7290 INAF-Astrophysical Observatory of Torino,
Department of Physics, University of Torino and INFN-Torino, 10125 Torino,
Italy  S. Moriyama  0000-0001-7630-2839 Kamioka Observatory, Institute for
Cosmic Ray Research, and Kavli Institute for the Physics and Mathematics of
the Universe (WPI), University of Tokyo, Higashi-Mozumi, Kamioka, Hida, Gifu
506-1205, Japan  K. Morå  0000-0002-2011-1889 Physics Department, Columbia
University, New York, NY 10027, USA  Y. Mosbacher  Department of Particle
Physics and Astrophysics, Weizmann Institute of Science, Rehovot 7610001,
Israel  M. Murra  0009-0008-2608-4472 Physics Department, Columbia University,
New York, NY 10027, USA  J. Müller  0009-0007-4572-6146 Physikalisches
Institut, Universität Freiburg, 79104 Freiburg, Germany  K. Ni
0000-0003-2566-0091 Department of Physics, University of California San Diego,
La Jolla, CA 92093, USA  U. Oberlack  0000-0001-8160-5498 Institut für Physik
& Exzellenzcluster PRISMA+, Johannes Gutenberg-Universität Mainz, 55099 Mainz,
Germany  B. Paetsch  0000-0002-5025-3976 Department of Particle Physics and
Astrophysics, Weizmann Institute of Science, Rehovot 7610001, Israel  Y. Pan
0000-0002-0812-9007 LPNHE, Sorbonne Université, CNRS/IN2P3, 75005 Paris,
France  Q. Pellegrini  0009-0002-8692-6367 LPNHE, Sorbonne Université,
CNRS/IN2P3, 75005 Paris, France  R. Peres  0000-0001-5243-2268 Physik-
Institut, University of Zürich, 8057 Zürich, Switzerland  C. Peters
Department of Physics and Astronomy, Rice University, Houston, TX 77005, USA
J. Pienaar  0000-0001-5830-5454 Department of Physics & Kavli Institute for
Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
Department of Particle Physics and Astrophysics, Weizmann Institute of
Science, Rehovot 7610001, Israel  M. Pierre  0000-0002-9714-4929 Nikhef and
the University of Amsterdam, Science Park, 1098XG Amsterdam, Netherlands  G.
Plante  0000-0003-4381-674X Physics Department, Columbia University, New York,
NY 10027, USA  T. R. Pollmann  0000-0002-1249-6213 Nikhef and the University
of Amsterdam, Science Park, 1098XG Amsterdam, Netherlands  L. Principe
0000-0002-8752-7694 SUBATECH, IMT Atlantique, CNRS/IN2P3, Université de
Nantes, Nantes 44307, France  J. Qi  0000-0003-0078-0417 Department of
Physics, University of California San Diego, La Jolla, CA 92093, USA  J. Qin
0000-0001-8228-8949 Department of Physics and Astronomy, Rice University,
Houston, TX 77005, USA  D. Ramírez García  0000-0002-5896-2697
<EMAIL_ADDRESS>Physik-Institut, University of Zürich, 8057
Zürich, Switzerland  M. Rajado  0000-0002-7663-2915 Physik-Institut,
University of Zürich, 8057 Zürich, Switzerland  R. Singh  0000-0001-9564-7795
Department of Physics and Astronomy, Purdue University, West Lafayette, IN
47907, USA  L. Sanchez  0009-0000-4564-4705 Department of Physics and
Astronomy, Rice University, Houston, TX 77005, USA  J. M. F. dos Santos
0000-0002-8841-6523 LIBPhys, Department of Physics, University of Coimbra,
3004-516 Coimbra, Portugal  I. Sarnoff  0000-0002-4914-4991 New York
University Abu Dhabi - Center for Astro, Particle and Planetary Physics, Abu
Dhabi, United Arab Emirates  G. Sartorelli  0000-0003-1910-5948 Department of
Physics and Astronomy, University of Bologna and INFN-Bologna, 40126 Bologna,
Italy  J. Schreiner  Max-Planck-Institut für Kernphysik, 69117 Heidelberg,
Germany  D. Schulte  Institut für Kernphysik, Westfälische Wilhelms-
Universität Münster, 48149 Münster, Germany  P. Schulte  0009-0008-9029-3092
Institut für Kernphysik, Westfälische Wilhelms-Universität Münster, 48149
Münster, Germany  H. Schulze Eißing  0009-0005-9760-4234 Institut für
Kernphysik, Westfälische Wilhelms-Universität Münster, 48149 Münster, Germany
M. Schumann  0000-0002-5036-1256 Physikalisches Institut, Universität
Freiburg, 79104 Freiburg, Germany  L. Scotto Lavina  0000-0002-3483-8800
LPNHE, Sorbonne Université, CNRS/IN2P3, 75005 Paris, France  M. Selvi
0000-0003-0243-0840 Department of Physics and Astronomy, University of Bologna
and INFN-Bologna, 40126 Bologna, Italy  F. Semeria  0000-0002-4328-6454
Department of Physics and Astronomy, University of Bologna and INFN-Bologna,
40126 Bologna, Italy  P. Shagin  0009-0003-2423-4311 Institut für Physik &
Exzellenzcluster PRISMA+, Johannes Gutenberg-Universität Mainz, 55099 Mainz,
Germany  S. Shi  0000-0002-2445-6681 Physics Department, Columbia University,
New York, NY 10027, USA  J. Shi  Department of Physics & Center for High
Energy Physics, Tsinghua University, Beijing 100084, P.R. China  M. Silva
0000-0002-1554-9579 LIBPhys, Department of Physics, University of Coimbra,
3004-516 Coimbra, Portugal  H. Simgen  0000-0003-3074-0395 Max-Planck-Institut
für Kernphysik, 69117 Heidelberg, Germany  A. Takeda  0009-0003-6003-072X
Kamioka Observatory, Institute for Cosmic Ray Research, and Kavli Institute
for the Physics and Mathematics of the Universe (WPI), University of Tokyo,
Higashi-Mozumi, Kamioka, Hida, Gifu 506-1205, Japan  P.-L. Tan
0000-0002-5743-2520 Oskar Klein Centre, Department of Physics, Stockholm
University, AlbaNova, Stockholm SE-10691, Sweden  A. Terliuk  [  Max-Planck-
Institut für Kernphysik, 69117 Heidelberg, Germany  D. Thers
0000-0002-9052-9703 SUBATECH, IMT Atlantique, CNRS/IN2P3, Université de
Nantes, Nantes 44307, France  F. Toschi  0009-0007-8336-9207 Institute for
Astroparticle Physics, Karlsruhe Institute of Technology, 76021 Karlsruhe,
Germany  G. Trinchero  0000-0003-0866-6379 INAF-Astrophysical Observatory of
Torino, Department of Physics, University of Torino and INFN-Torino, 10125
Torino, Italy  C. D. Tunnell  0000-0001-8158-7795 Department of Physics and
Astronomy, Rice University, Houston, TX 77005, USA  F. Tönnies
0000-0002-2287-5815 Physikalisches Institut, Universität Freiburg, 79104
Freiburg, Germany  K. Valerius  0000-0001-7964-974X Institute for
Astroparticle Physics, Karlsruhe Institute of Technology, 76021 Karlsruhe,
Germany  S. Vecchi  0000-0002-4311-3166 INFN-Ferrara and Dip. di Fisica e
Scienze della Terra, Università di Ferrara, 44122 Ferrara, Italy  S. Vetter
0009-0001-2961-5274 Institute for Astroparticle Physics, Karlsruhe Institute
of Technology, 76021 Karlsruhe, Germany  F. I. Villazon Solar  Institut für
Physik & Exzellenzcluster PRISMA+, Johannes Gutenberg-Universität Mainz, 55099
Mainz, Germany  G. Volta  0000-0001-7351-1459 Max-Planck-Institut für
Kernphysik, 69117 Heidelberg, Germany  C. Weinheimer  0000-0002-4083-9068
Institut für Kernphysik, Westfälische Wilhelms-Universität Münster, 48149
Münster, Germany  M. Weiss  0009-0005-3996-3474 Department of Particle Physics
and Astrophysics, Weizmann Institute of Science, Rehovot 7610001, Israel  D.
Wenz  0009-0004-5242-3571 Institut für Kernphysik, Westfälische Wilhelms-
Universität Münster, 48149 Münster, Germany  C. Wittweg  0000-0001-8494-740X
Physik-Institut, University of Zürich, 8057 Zürich, Switzerland  V. H. S. Wu
0000-0002-8111-1532 Institute for Astroparticle Physics, Karlsruhe Institute
of Technology, 76021 Karlsruhe, Germany  Y. Xing  0000-0002-1866-5188
SUBATECH, IMT Atlantique, CNRS/IN2P3, Université de Nantes, Nantes 44307,
France  D. Xu  0000-0001-7361-9195 Physics Department, Columbia University,
New York, NY 10027, USA  Z. Xu  0000-0002-6720-3094<EMAIL_ADDRESS>Physics Department, Columbia University, New York, NY 10027, USA  M. Yamashita
0000-0001-9811-1929 Kamioka Observatory, Institute for Cosmic Ray Research,
and Kavli Institute for the Physics and Mathematics of the Universe (WPI),
University of Tokyo, Higashi-Mozumi, Kamioka, Hida, Gifu 506-1205, Japan  L.
Yang  0000-0001-5272-050X Department of Physics, University of California San
Diego, La Jolla, CA 92093, USA  J. Ye  0000-0002-6127-2582 School of Science
and Engineering, The Chinese University of Hong Kong, Shenzhen, Guangdong,
518172, P.R. China  L. Yuan  0000-0003-0024-8017 Department of Physics & Kavli
Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637,
USA  G. Zavattini  0000-0002-6089-7185 INFN-Ferrara and Dip. di Fisica e
Scienze della Terra, Università di Ferrara, 44122 Ferrara, Italy  M. Zhong
0009-0004-2968-6357 Department of Physics, University of California San Diego,
La Jolla, CA 92093, USA<EMAIL_ADDRESS>
(September 3, 2024)

###### Abstract

The XENONnT experiment searches for weakly-interacting massive particle (WIMP)
dark matter scattering off a xenon nucleus. In particular, XENONnT uses a
dual-phase time projection chamber with a 5.9-tonne liquid xenon target,
detecting both scintillation and ionization signals to reconstruct the energy,
position, and type of recoil. A blind search for nuclear recoil WIMPs with an
exposure of 1.1 tonne-years yielded no signal excess over background
expectations, from which competitive exclusion limits were derived on WIMP-
nucleon elastic scatter cross sections, for WIMP masses ranging from 6
GeV/$c^{2}$ up to the TeV/$c^{2}$ scale. This work details the modeling and
statistical methods employed in this search. By means of calibration data, we
model the detector response, which is then used to derive background and
signal models. The construction and validation of these models is discussed,
alongside additional purely data-driven backgrounds. We also describe the
statistical inference framework, including the definition of the likelihood
function and the construction of confidence intervals.

Dark Matter, Direct Detection, Xenon

##  I Introduction

Astrophysical and cosmological observations at various scales provide clear
evidence for the existence of dark matter (DM) as a fundamental building block
of our universe [1, 2]. Numerous extensions of the standard model of particle
physics predict additional fundamental particles as potential candidates for
DM [3, 4]. Measuring the direct interaction of particles from the DM halo of
the Milky Way with a suitable detector could provide conclusive evidence for
the particle DM hypothesis. The XENONnT experiment aims to detect a signal
from weakly interacting massive particles (WIMPs) elastically scattering off
xenon nuclei with a detector operated underground at the INFN Laboratori
Nazionali del Gran Sasso (LNGS) in Italy. The collaboration published the
first WIMP search results from a data-taking period between July
$6^{\mathrm{th}}$ and November $10^{\mathrm{th}}$ 2021 with a total live time
of $95.1$ days and a fiducial mass of $(4.18\pm 0.13)\;\mathrm{t}$, referred
to as science run 0 (SR0). We performed a blind analysis and observed no
significant excess above background expectations [5]. This led to a minimum
upper limit on the spin-independent WIMP-nucleon cross section of $2.58\times
10^{-47}\;\mathrm{cm}^{2}$ for a WIMP mass of
$28\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ at $90\%$ conﬁdence
level.

The detector used for the XENONnT experiment is a dual-phase xenon time
projection chamber (TPC). A particle interaction in the active
$5.9\text{\,}\mathrm{t}$ liquid xenon (LXe) target results in prompt vacuum
ultraviolet scintillation light as well as free ionization electrons. The
instantaneous light signal, called the S1 signal, is detected by two arrays of
photomultiplier tubes (PMTs) at the top and bottom of the approximately
cylindrical TPC. An electric field
($23\text{\,}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ average electric drift
field during SR0) is applied to the target volume to drift the ionization
electrons toward the liquid xenon surface. Here, they are extracted and
accelerated into a xenon gas volume above the liquid via a stronger electric
field ($2.9\text{\,}\mathrm{k}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$
electric extraction field during SR0). The kinetic energy of the accelerated
electrons in the gas phase is sufficient for the emission of
electroluminescent light, which is proportional to the number of extracted
electrons. This second light signal is referred to as the S2 signal. The
distribution of the S2 signal detected by the PMTs in the top array is used to
infer the position of the initial interaction in the horizontal plane
$\mathrm{(\mathrm{X},\mathrm{Y})}$ parallel to the liquid surface, which
defines the radial coordinate
$\mathrm{\mathrm{R}=\sqrt{\mathrm{X}^{2}+\mathrm{Y}^{2}}}$. The vertical
component $\mathrm{Z}$ follows from the drift time of electrons, which is
determined as the time difference between the S1 and S2 signals. This three-
dimensional position reconstruction of the interaction vertex enables the
discrimination of multi-site events as well as the selection of an inner
fiducial volume (FV) within the TPC. This volume benefits from a particularly
low background level thanks to the excellent self-shielding properties of LXe,
given its high density. Moreover, the relative magnitudes of the S1 and S2
signals can be used to determine if a particle interacted with the xenon
nucleus (expected from WIMPs) or xenon shell electrons (typical of the
dominant sources of background), i.e. nuclear recoil (NR) or electronic recoil
(ER) events. The electric drift and extraction fields are established by three
parallel-wire electrodes (cathode, gate, and anode from bottom to top). To
reduce wire sagging, two and four horizontal perpendicular wires support the
gate and anode electrodes, respectively.

The TPC is nested in an active neutron veto, which in turn is housed in an
active muon veto. Both veto systems are water Cherenkov detectors. More
details on the TPC, the veto systems, and other subsystems of the detector are
provided in [6, 7].

The data analysis chain of the WIMP DM search in XENONnT is subdivided into
two major parts, like in XENON1T [8, 9]. The first part is discussed in [10]
and comprises signal and event reconstruction, corrections, event selection,
and energy scale calibration. Here we present the second part of the analysis
chain: The detector response to the different interaction types is modeled
based on calibration data and discussed in Sec. II. The derived best-fit
models are important inputs for the definition of signal and background
models, detailed in Sec. III. This section also covers data-driven background
models that do not rely on the detector response. Finally, the statistical
methods used to compute constraints on DM given the XENONnT data are discussed
in Sec. IV.

##  II Modeling the detector response

Modeling the detector response to energy depositions allows for a powerful
discrimination between NR signal and ER background in S1–S2 space. NR events
are produced by particles scattering elastically off xenon nuclei, while in ER
events, particles scatter off xenon shell electrons. The different energy loss
processes involved after the two types of recoil cause ER and NR events to
form separate populations in S1–S2 space, which is exploited in the WIMP DM
search.

In this section, we describe the TPC response model to energy depositions in
LXe via ER or NR interactions, which is obtained from fits to calibration
data. This detector response model is separated into two parts. The first one
is the empirically parametrized LXe emission model, which describes the
production processes of the detectable quanta, i.e. the conversion of
deposited energy into scintillation photons and ionization electrons. The
second part is the detector reconstruction model, which covers the conversion
from the produced photons and electrons into the observed and spacetime-
dependence corrected S1 and S2 signals (cS1 and cS2). Due to its complexity,
it is unfeasible in the fit to simulate processes of photon and electron
propagation on a per-quantum basis. Instead, toy-Monte-Carlo (toy-MC)
simulations of the detector reconstruction model are used which sample from
effective maps provided either by data-driven analyses described in [10] or by
simulation-driven analyses using the XENONnT Monte Carlo (MC) framework [11,
12]. Then all the model parameters in the simulations are fit to both ER and
NR calibration data. The two parts of the detector response model together
with the toy-MC simulation workflow are described in Secs. II.1 and II.2, and
the fit to ER and NR calibration data is detailed in Sec. II.3.

###  II.1 Liquid Xenon Emission Model

The production of quanta from energy depositions in a xenon target is complex
and lacks a comprehensive description derived from first principles. Thus, a
semi-empirical parametrization of the emission model is commonly used. The
parametrization used in XENONnT SR0 is similar to XENON1T [9], which is based
on the Noble Element Simulation Technique (NEST) model described in [13, 14].
The simulation workflow of the emission model is described in detail in
Appendices A and B, and is summarized in the following.

In an ER event, the recoil energy is transmitted into the excitation and
ionization of xenon atoms. Recoiling xenon nuclei from an NR event, on the
other hand, lose a significant amount (roughly 80 %) of their kinetic energy
to atomic motion in collisions with other xenon atoms [13]. This
thermalization process is undetectable in LXe TPCs, resulting in an effective
signal loss. For both recoil types, the total number of detectable quanta
$N_{\mathrm{q}}$ equals the sum of the number of produced excitons
$N_{\mathrm{ex}}$ and ions $N_{\mathrm{i}}$. A fraction of ions and electrons
recombine, resulting in the production of additional excitons. An exciton can
combine with a ground-state xenon atom forming an excimer, which de-excites,
producing a scintillation photon. In an event, $N_{\gamma}$ photons are
generated, part of which are then detected by PMTs as the $\mathrm{S1}$
signal. The $N_{\mathrm{e}}$ electrons that do not recombine are drifted to
the liquid-gas interface and eventually form an $\mathrm{S2}$ signal.

In the region of interest (ROI) for the SR0 WIMP search, $\mathrm{cS1}$ in [0,
100] PE and $\mathrm{cS2}$ in [$10^{2.1}$, $10^{4.1}$] PE, NR interactions
have a larger recombination probability than ER interactions given the same
total energy deposition. This leads to smaller ratios of S2 to S1 areas, which
is the principle behind ER-NR discrimination in LXe TPCs. We parametrize the
mean recombination fraction for ER interactions using the Thomas-Imel box
model [15] with an additional Fermi-Dirac term, as in [9]. For NR
interactions, the emission model follows the parametrization used in the NEST
v1 model detailed in [14]. While there are newer versions of fits available in
the NEST framework that use a different parametrization for NR interactions
(which we refer to as NEST v2 in the following), we choose to stay with the
previous version (referred to as NEST v1). The model is simpler, fits our data
well, and the best-fit values and uncertainties of parameters of NEST v2 were
not published at the time of the analysis.

###  II.2 Detector Reconstruction Model

The detector reconstruction model covers all aspects of the signal
reconstruction, starting from the produced scintillation photons and
ionization electrons up to the measured S1 and S2 signal sizes. The XENONnT
detector reconstruction model is almost identical to the one presented in [9]
and is briefly outlined here. The detailed simulation workflow of the detector
reconstruction model is described in Appendix C.

For S1 signals, the spatial dependence of geometric effects during photon
propagation, the photon detection efficiency of the PMTs, and the effect of
double photoelectron emission (DPE) from the PMT photocathode [16, 17] are all
modeled in the simulations. For S2 signals, we model the loss of electrons
along their drift path due to attachment to electronegative impurities as well
as due to electric field effects close to the TPC wall, the efficiency of
electron extraction from the liquid into the gas phase, and the single-
electron gain of extracted electrons to detected photoelectrons (PE) via the
electroluminescence process in gas. In addition to these physical processes,
we also account for the influence of software reconstruction effects on the
signals. For both signal types, the software reconstruction process can
introduce biases and fluctuations in the recorded signal size. We also account
for the S1 signal reconstruction efficiency, which vanishes at about 3 PE, as
we require that at least 3 PMTs detect a photon within $\pm$50 ns around the
maximal amplitude. Finally, we model the uncertainty of the event position
reconstruction and the acceptances of data quality selections [10].

The modeling of ER interactions relies on calibration data from sources
dissolved in xenon, in particular 220Rn and 37Ar. The 220Rn progeny 212Pb
undergoes a $\beta$-decay with a Q-value at about 0.6 MeV [18], giving an
approximately flat ER energy spectrum in the WIMP ROI, i.e. below about $10$
keV. This dataset is used for the fit of the ER model. The 37Ar source
provides a mono-energetic ER peak when its K-shell electron capture process
leads to a fast cascade of X-rays and Auger-Meitner electrons, with a total
energy of about 2.8 keV [19]. Both ER sources provide single-scatter (SS)
events, producing one S1 and one S2 signal for the total deposited energy, and
the above-described detector reconstruction processes are simulated
accordingly.

The NR model is calibrated using an external 241AmBe neutron source (referred
to as AmBe source in the following). It emits a fast neutron via the
9Be$(\alpha,n)^{12}$C capture reaction. Neutrons often scatter multiple times
in the LXe target in so-called multi-scatter (MS) events, creating distinct
energy depositions at each interaction site, thus the topologies of neutron
events must be considered in the modeling. Since the neutrons travel at about
10 % of the speed of light and the S1 width is
$\sim$$\mathcal{O}(\mathrm{100\,ns})$, the S1 signals from the separate energy
depositions in a MS event get reconstructed into one merged S1 signal. The S2
signals might be distinguished based on the separation in spacetime between
the energy depositions of the individual scatters, but can get (partially)
merged if the interaction sites are close together. Due to the low drift field
of $23\,\mathrm{V/cm}$ in XENONnT SR0, the $\mathrm{Z}$ separation resolution
between two S2 signals is reduced with respect to the design because of the
increased drift time $\sim$$\mathcal{O}(\mathrm{ms})$ and consequently larger
S2 spread $\sim$$\mathcal{O}(\mathrm{\upmu s})$ due to the diffusion of
electrons. This effect is accounted for in the fit of the NR model to the
neutron calibration data: we first produce the photons and electrons for each
energy deposition in a MS event separately, based on the emission model. While
the primary scintillation photons are summed to produce the merged S1 signal,
we only sum S2 signals that will become part of the largest S2 signal after
going through the remaining detector reconstruction processes. This
information is provided by PMT waveform simulations of neutron scatters from
the AmBe source, using the Monte Carlo (MC) framework of XENONnT [11, 12] (see
also Sec. III.3 for details on neutron simulations). The results of these
waveform simulations are used as inputs to the toy-MC simulations and fitting
framework. Data selections that remove MS events (both resolved and
unresolved) are directly applied to these inputs, because their acceptances
can only be determined from full simulations, and are correspondingly dropped
from the selection acceptance curves for the NR fit.

Figure 1: AmBe neutron calibration events with (red) and without (gray) a
coincident signal in the neutron veto (NV). Selecting coincident events
ensures a clean nuclear recoil sample for detector response modeling. The
accidental coincidence population (cS1 below 5 PE) and misidentified single-
electron S1s (cS1 between 10 and 30 PE) are visible in the gray population.
After applying additional data quality selections, approximately 2000 events
(Fig. 2, right) are used for band fitting.

###  II.3 Fit to Calibration Data

Figure 2: Comparison between calibration data and the best-fit ER (left) and
NR (right) models. The equiprobable binning for the 2D binned Poisson
likelihood $\chi^{2}$ goodness-of-fit tests is shown. The color scale
indicates the deviation of the number of data points (overlaid as black dots)
in each bin from the best-fit model expectation $\mu_{\mathrm{bin}}$ in units
of sigma, which is 29.3 for the ER model and 36.9 for the NR model.

The parameters in the full detector response model are fitted to 220Rn, 37Ar,
and AmBe calibration data. The ROI selection is the same as the SR0 WIMP
search region defined in Sec. II.1. All data quality selections are applied to
the calibration data within the ROI. For the AmBe neutron data, an additional
selection is applied using coincident gammas detected in the neutron veto. In
9Be$(\alpha,n)^{12}$C capture reactions in the AmBe source, 12C is left in the
first excited state after neutron emission in about 50 % of the cases,
estimated with neutron veto data [20]. The 12C excited state promptly decays
via the emission of a 4.44 MeV gamma ray, coincident with the neutron
emission. The requirement of the time coincidence within $250\,\mathrm{\upmu
s}$ between the 4.44 MeV gamma ray detected in the neutron veto and the event
in the TPC, leads to a clean calibration dataset, with only about one expected
event from accidental coincidences between the two detectors in the dataset
($\sim$0.05 % of the resulting sample). The events recorded in the ROI that
pass (fail) the neutron veto coincidence selection are shown in red (gray) in
Fig. 1.

For the fit to the ER data, an additional background component is added to
account for accidental coincidence (AC, see subsection III.5 for details)
events formed by a wrong pairing of S1 and S2 signals in the calibration
dataset.

After all selections, the 220Rn and AmBe datasets consist of $\sim$2000 events
each. We use a downsampled dataset of 10 000 37Ar events to perform a combined
220Rn-37Ar ER fit. With the downsampling we avoid overconstraining the fit to
the narrow, very low energy range of 37Ar.

Figure 3: Photon (top) and charge (bottom) yields as functions of deposited
energy for ER events (left) and NR events (right). Dark red lines and shaded
bands show the XENONnT liquid xenon emission models with the best-fit
parameters and their uncertainties, respectively. The fit results from XENON1T
(blue shaded bands) [9, 21], and the models from NEST v1 (gray solid line) and
NEST v2 (black solid line), both evaluated at 23 V/cm, are shown as well.
Several measurements from literature performed at various drift fields as
labeled in the legends (ER yields from [9, 19, 22, 23, 24, 25], NR yields from
Aprile 2005 [26], Aprile 2006 [27], Aprile 2009 [28], XENON100 [29], Plante
2011 [30], Sorensen 2009 [31], Manzur 2010 [32], LUX 2016 [33], LUX 2022 [34])
are also shown. The gray shaded regions mark the energy ranges in which the
total detection + selection efficiency in the WIMP analysis is above 10 %.

For the fit, an affine-invariant Markov chain Monte Carlo (MCMC) algorithm
[35, 36] is used to sample from the high-dimensional posterior distribution of
the parameters. In each step of the sampling, a GPU-accelerated toy-MC
simulation is performed for every set of the ER and NR model parameters,
following the steps described in the previous two sections, producing a model
of the detector response in the space of corrected signals, i.e.
$\mathrm{cS1}$–$\mathrm{cS2}$. The model parameters are described in detail in
the appendix and are listed in Tab. 2. These toy-MC models are then compared
to calibration data by calculating the binned Poisson likelihood in
$\mathrm{cS1}$–$\mathrm{cS2}$ space. The likelihoods are multiplied by the
prior distributions of the parameters to yield the posteriors. After
sufficient iterations of the MCMC, the samples in the chain then converge to
the posterior distribution of the parameters.

In the ER emission model, the parametrization is the same as in XENON1T [9]
with very wide flat priors, referred to as free priors. For the NR emission
model, the parameter priors are taken from [14]. A Gaussian-like distribution
with different widths on either side is used as the prior of parameters with
asymmetric uncertainties. However, due to the low drift field of 23 V/cm,
where no literature data on NR yields exist, the validity of the field
dependence in the model is unverified. Therefore, the field dependence in the
emission model from [14] was modified. Two scaling parameters were freed
($\alpha$ and $\gamma$ in Eq. 24 and Eq. 27 in the appendix), and two field
dependence parameters ($\zeta$ and $\delta$) were fixed to the reported best-
fit values. The widths of all other parameter priors were doubled in order to
allow more freedom in the fit. Because the dependence of NR signals on the
drift field has been shown to be small (see also literature measurements shown
in Fig. 3), we consider this a reasonable choice.

The fit results of all emission model parameters are summarized in Tab. 2 in
the appendix. No significant tension between the posteriors and the priors of
any parameter is found. Suitable goodness-of-fit (GOF) tests are performed to
assess whether the best-fit models adequately describe the 220Rn and AmBe
calibration data. Specifically, we employ binned Poisson likelihood $\chi^{2}$
tests in the $\mathrm{cS1}$–$\mathrm{cS2}$ space. We adopt an equiprobable
binning scheme using the GOFevaluation package [37], ensuring that the number
of expected events under the best-fit model is the same in each bin,
$\mu_{\mathrm{bin}}$. In Fig. 2 the results of the 2D
$\mathrm{cS1}$–$\mathrm{cS2}$ Poisson likelihood $\chi^{2}$ test are shown for
both the ER and the NR fit, overlaid with the 220Rn and the AmBe calibration
data, respectively. The resulting p-values of the tests indicate no
significant discrepancy between the best-fit models and the calibration data.

Figure 3 shows the photon and charge yields as a function of deposited energy
for the ER and NR emission models using our best-fit parameters (dark red).
For comparison, NEST models calculated at a drift field of 23 V/cm and several
published measurements at different drift fields are shown. For the NR yields,
we also compare our results to the NEST v1 model at a drift field of 23 V/cm
(red dashed lines in Fig. 3) which uses the same parametrization as the model
in this work (see Sec. II.1), showing very good agreement.

Figure 4: Distribution of each background component (colored) in
$\mathrm{cS1}$–$\mathrm{cS2}$ space. The NR background includes both
radiogenic neutron background and NR events by neutrinos. The probability
density functions (PDFs) of backgrounds are shown as 1$\sigma$ (2$\sigma$)
dark (light) regions, containing 68% (95%) expected events in the ROI. For
reference, the contours are also shown for a spin-independent 200 GeV$/c^{2}$
WIMP as a solid (dashed) dark gray line. The light gray dash-dotted lines are
contours of constant NR-equivalent energy (and ER-equivalent energy in the top
left panel) in units of $\mathrm{keV}$.

##  III Signal and background modeling

The dominant sources of background in the WIMP search are ER events from
intrinsic $\beta$-decays from materials and neutrinos (Sec. III.2), NR events
from radiogenic neutrons from detector materials (Sec. III.3) and coherent
neutrino scattering (Sec. III.4), AC events (Sec. III.5), and surface
background (Sec. III.6). Besides their total expected rates, the distribution
of these background events in the analysis space cS1–cS2–$\mathrm{R}$ is
derived. In Fig. 4, the distributions of background models in the
$\mathrm{cS1}$–$\mathrm{cS2}$ space are shown, compared to the signal model
for elastic spin-independent scattering of a 200 GeV$/c^{2}$ mass WIMP.

###  III.1 WIMP Signal Model

For non-relativistic, spin-independent elastic WIMP-nucleus coherent
scattering, the event rate $R_{\mathrm{WIMP}}$ scales with the square of the
atomic mass number $A^{2}$ of the target and can be written as [38]

$\frac{\mathrm{d}R_{\mathrm{WIMP}}}{\mathrm{d}E}=\frac{\rho_{\chi}}{2m_{\chi}\mu_{\chi
N}^{2}}\left\langle\frac{1}{v}\right\rangle\sigma A^{2}F^{2}(q),$ (1)

where $E$ is the recoil energy, $\rho_{\chi}$ is the local DM density of 0.3
GeV/($c^{2}\,$cm3), $m_{\chi}$ is the WIMP mass, $\mu_{\chi N}$ is the WIMP-
nucleon reduced mass, $v$ is the WIMP velocity in the lab-frame, $\sigma$ is
the WIMP-nucleon cross section, and $F(q)$ is the nuclear form factor as a
function of the momentum transfer $q$ to the xenon nucleus [39]. The DM
velocity distribution is averaged using the parameter values of the standard
halo model with values from [40, 41, 42, 43, 44, 45], as recommended in [46].
Fig. 5 shows the cS1–cS2 distribution of WIMP-nucleus scattering for different
WIMP masses. For increasing masses, the $68\%$ contours extend to higher cS1
and cS2 values, up to about 200 GeV/$c^{2}$, where the shape does not change
significantly anymore due to kinematic effects. For this reason, confidence
intervals in $\sigma$ (see Sec. IV) for $m_{\chi}\gtrsim
200\,\mathrm{GeV}/c^{2}$ approximately scale proportional to the inverse of
the assumed WIMP mass $m_{\chi}$.

For the spin-dependent WIMP-nucleus interaction, assuming natural abundances
of xenon isotopes, only 129Xe and 131Xe can contribute since they are the only
two stable isotopes of xenon with non-zero nuclear spin $J$. The cross section
is usually written as

$\frac{\mathrm{d}\sigma_{\mathrm{SD}}}{\mathrm{d}q^{2}}=\frac{\sigma}{3\mu_{\chi
N}^{2}v^{2}}\frac{\pi}{2J+1}S_{N}(q).$ (2)

The axial-vector structure factor of xenon $S_{N}$ is taken from [47]. In the
XENONnT SR0 analysis, the mean of the structure factor is used and the
uncertainty is neglected since it is much smaller than the uncertainty from
the NR fit.

Uncertainties from the posterior distribution of the NR model parameters and
other efficiencies can be propagated to the NR rate uncertainty. For a 6
GeV$/c^{2}$ (50 GeV$/c^{2}$) WIMP, the relative rate uncertainty is $\sim$30%
(10%). The NR model shape can also be affected by the posterior distribution.
However, because of the low statistics of NR events in SR0 (for both
background and signal expectation), the impact of the shape uncertainty of the
NR model is negligible compared to the uncertainty of its absolute rate. Thus,
the NR model shape uncertainty is not propagated to the final likelihood.

Figure 5: Expected distribution in cS1–cS2 space from spin-independent WIMP
nuclear recoil for WIMP masses from 6 GeV/$c^{2}$ to 1000 GeV/$c^{2}$. The
contours contain 68% of expected events in the ROI. The blinded region is
shown for reference as a gray band and dash-dotted lines are constant NR-
equivalent energy contours.

###  III.2 Electronic Recoil Background Model

Figure 6: Shape uncertainty of the ER model constrained by the 220Rn
calibration data. The blue (gray) solid and dashed lines represent the
contours of ER background (200 GeV$/c^{2}$ WIMP) containing 68% and 95% of
expected events in the ROI. The transparent blue regions show the shape
uncertainty by varying the two shape nuisance parameters by $\pm 1\sigma$.
Dash-dotted lines are contours of constant ER-equivalent energy in units of
$\mathrm{keV}_{\mathrm{ER}}$.

ER events are one of the dominant background sources in the XENONnT SR0 WIMP
search, primarily originating from $\beta$-decays of intrinsic radioactive
contaminants such as 214Pb (a product of the 222Rn decay chain) and 85Kr.
Contributions from 136Xe double-$\beta$ decays, solar neutrino interactions,
and gamma events from the materials are also taken into account. In the ROI
for the WIMP search, the total ER energy spectrum is approximately flat. The
distribution of ER background events in $\mathrm{cS1}$–$\mathrm{cS2}$ space is
generated from the ER model fitted to 220Rn and 37Ar calibration data, as
discussed in Sec. II.3. The total ER rate is left unconstrained and is fitted
in the final WIMP likelihood, which will be discussed in Sec. IV.

In contrast to the NR model, the uncertainty in the shape of the ER model can
affect the sensitivity of the WIMP search. Ideally, the posterior distribution
of all parameters should be propagated to the ER model as nuisance parameters
in the WIMP search likelihood function. However, an excessive number of
correlated nuisance parameters becomes computationally challenging. In XENONnT
SR0, a principal component analysis (PCA) [48] is used to remove correlation
among parameters from the MCMC sampler. All ER parameters shown in Tab. 2,
together with $W$, $g_{1}$, and $g_{2}$ (see details in the appendix), are
included in the PCA. Different from the original PCA, the variance of

$t\equiv\sum_{i}\frac{s_{i}^{2}}{s_{i}+b_{i}}$ (3)

along each principal component is used to quantify its importance. Here,
$s_{i}$ and $b_{i}$ are the probability of a 50 GeV$/c^{2}$ WIMP and ER
background in the $i^{\text{th}}$ bin of cS1–cS2 space, respectively. A larger
variance of $t$ means that the uncertainty of that component can be more
impactful for the WIMP sensitivity. In the end, the two most important
principal components are kept as nuisance parameters. Their shape uncertainty,
constrained by the calibration data, is shown in Fig. 6.

###  III.3 Neutron Background Model

Figure 7: Spatial (left) and cS1–cS2 (right) distribution of events tagged by
the neutron veto (red markers). Black dashed and solid lines in the left plot
indicate the FV and physical TPC boundary, respectively. The regions
containing 68% (95%) of the total expected neutron background are indicated in
dark (light) yellow. The neutron veto tagging yields four SS events (circles)
and seven MS events, from which three are reconstructed with the main S2
signal inside the FV (squares) and the other four outside (diamonds). Solid
and hollow markers represent the reconstructed positions of the largest and
the second-largest S2 signal, respectively. The right plot only shows events
inside the FV; from these, the four events within the neutron background band
(large markers; one SS and three MS) were used to estimate the neutron
background rate. The other three SS events are in agreement with the
prediction of falsely tagged ER and surface background events (c.f. Fig. 4).
For reference, the events surviving all data selections, not tagged by the
neutron veto, are displayed as gray dots. Also here, the right panel displays
only those data points within the FV, which correspond to the unblinded WIMP
search events.

Radiogenic neutrons are mainly produced through spontaneous fission and
$(\alpha,n)$ reactions in the detector materials due to intrinsic traces of
radioactive impurities. The cosmogenic neutron background is subdominant [49],
hence it is neglected in XENONnT SR0 analysis. The neutron yield and energy
per isotope and material are calculated with the SOURCES-4A software [50],
using the radioactive impurity levels of the relevant detector components
obtained via the combined XENON1T and XENONnT radioassay campaigns [51, 52]. A
full-chain simulation pipeline [11] is used to estimate the neutron background
rate, the geometrical distribution, and the cS1–cS2 distributions for SR0.

The propagation of the neutrons in the XENONnT detector is simulated with the
Geant4 toolkit [53, 54], where the recoil type and energy deposition per
interaction in the target are recorded. We compute the number of photons and
electrons for a given interaction via the custom-developed epix package [55],
which utilizes the energy-dependent LXe response derived from the SR0
calibrations, shown in Fig. 3. This package also handles the clustering of the
individual energy-depositing steps at the LXe microphysics scale before the
quanta generation. The photons and electrons produced by epix are passed to
the waveform simulator (WFSim) [56], which computes the S1 and S2 signals up
to the waveform level, by means of a precise set of simulations- and data-
driven corrections which characterize the XENONnT detector response. The
event-by-event simulated waveforms share the same data structure as the
science data after applying PMT and DAQ effects, which allows us to process
them with the same software used for the real data (straxen [57]), as well as
to apply the same data selections. The final SS neutron rate arises from
weighting the rates obtained via the entire waveform simulation pipeline with
the specific activities of the corresponding material and isotopic neutron
yield.

When a neutron scatters in the $\sim$250 kg of LXe between the cathode and the
bottom PMT array, only the scintillation light for these events can be
detected. The electrons, in turn, are lost due to an electric field pointing
in the opposite direction to the active volume. Neutron events that consist of
scatters above and below the cathode are referred to as neutron-X events,
where “X” means additional S1 contribution from charge-insensitive scatters.
They are modeled as a separate background since they have a larger cS1 to cS2
ratio than normal neutron scatters.

The event parameters having discrimination power on MS interactions, such as
S2 pulse shape and PMT hit patterns of S1, are matched and validated with
calibrations, such that the relevant data quality selections can be applied to
the simulation outputs. Notably, a validation of the multiple-to-single
scatter ratio and total rate in the TPC of the AmBe calibration data was
conducted, from where we obtained the systematic uncertainty associated with
the accuracy of the full-chain simulations.

With this agreement between simulations and calibration data, we decided
before unblinding the WIMP ROI to proceed with the sideband unblinding of the
events tagged by the neutron veto. Initially, this confirmed the simulation-
driven neutron background prediction. However, a mistake in the definition of
the neutron veto time window was found after the unblinding of the WIMP ROI.
After fixing this issue, a mismatch was found, with the neutron background
rate being larger than predicted in the ROI [5]. We therefore decided to
constrain the neutron background rate in a purely data-driven way based on the
aforementioned sideband unblinding with the correct neutron veto tagging
window. The results of the sideband unblinding are shown in Fig. 7. Based on
the multiple-to-single scatter ratio of 2.2, the $53\%$ neutron veto tagging
efficiency, and the three observed MS events and one SS neutron event tagged
by the neutron veto in the fiducial volume, a data-driven prediction of
$1.1^{+0.6}_{-0.5}$ events was derived for the SR0 exposure. We used the
simulation-driven ratio between normal neutron background and neutron+X events
to estimate the neutron+X event rate.

No further modification was propagated into the analysis after the data-
simulations rate mismatch was identified: the $250\,\mathrm{\upmu s}$ time
veto window between the TPC and the neutron veto, chosen due to the reduced
background rate initially predicted, and the fiducial volume of the WIMP
search remained as defined prior to unblinding. An underestimated
contamination from some of the surrounding materials is considered as a
possible cause of the discrepancy between the expected rate and the observed
neutron background. Studies to constrain the material’s radiopurity by means
of the high-energy gamma rays for specific detector regions are ongoing.

###  III.4 CE$\nu$NS Background Model

Neutrinos can also contribute to the NR background via coherent elastic
neutrino-nucleus scattering (CE$\nu$NS). Signals from solar, atmospheric, and
diffuse supernova neutrinos (DSN) will be in the WIMP search ROI. Due to the
weak interaction between neutrinos and nuclei, CE$\nu$NS events are expected
to be single-scatter and spatially uniform.

The recoil energy spectrum of solar CE$\nu$NS is almost identical to that of a
6 GeV$/c^{2}$ spin-independent WIMP [58], and the flux is $(5.25\pm
0.20)\times 10^{6}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ [59, 60]. After applying
all selections and their corresponding efficiencies, the expected number of
events is $(0.19\pm 0.06)$ in SR0, also shown in Tab. 1. Atmospheric neutrinos
and DSN mainly affect the search for heavier WIMPs. Their recoil energy
spectra are taken from [49], and the SR0 expectation value is $(0.05\pm 0.02)$
events. Due to the low cross section and the similarity to WIMP interactions,
CE$\nu$NS background will be the major limitation to the WIMP search
sensitivity of the next generation of LXe experiments.

###  III.5 Accidental Coincidence Background Model

AC events consist of incorrectly paired S1 and S2 signals. These S1 and S2
signals can occur, for example, when either the S1 or the S2 signal of a
physical event is not reconstructed due to detector effects, or when a single
electron S2 signal is misclassified as an S1 signal. Such signals are referred
to as “isolated” S1 and S2 signals. If an isolated S1 and an isolated S2 fall
within the event-building time interval [8], they form an AC event.

The AC background is modeled with a data-driven approach. Isolated S1 and S2
signals as well as their surrounding S1 and S2 signals that are close in time
are sampled and paired into events. S1 signals $<$150 PE are selected as
isolated S1 signals, and the isolated S2 signals are taken from events whose
S1 area is $<$150 PE, together with all pulses in the event window. The AC
event rate is computed as

$R_{\mathrm{AC}}=R_{\mathrm{isoS1}}\times R_{\mathrm{isoS2}}\times\Delta t,$
(4)

where $R_{\mathrm{isoS1}}$ and $R_{\mathrm{isoS2}}$ are the isolated S1 and S2
rates, and $\Delta t$ is the event-building time interval, which is defined
according to the maximum drift time of 2.2 ms in XENONnT. Occasionally during
SR0, some localized high rates of S2 signals appeared in the TPC. Excluding
these periods from the analysis, $R_{\mathrm{isoS1}}$ and $R_{\mathrm{isoS2}}$
remained stable throughout SR0 at an average rate of 1.5 Hz and 80 mHz,
respectively. In XENONnT, the isolated S2 rate is an order of magnitude higher
compared to XENON1T [61]. This can be explained by the lower electron
extraction efficiency which causes an increased rate of delayed electrons. The
isolated S1 and S2 signals are fed into the data processing pipeline [62, 57]
to reconstruct the events. This provides the background distribution of all
relevant event properties, especially the distribution in cS1 and cS2, as
shown in Fig. 4 (lower right).

Validation performed on calibration data provides the rate and shape
uncertainty of the AC model. For this, we compared the predicted and observed
AC events in a dedicated dataset by performing both 1D and 2D GOF tests with
equiprobable bins in all relevant parameters. The AC prediction is provided by
the AC modeling method discussed above, and the AC dataset is selected by
inverting the data selections to avoid other physical events. Because of the
large statistics of 37Ar calibration data, it delivers the most stringent
constraint on the AC model. 37Ar events with S2 areas smaller than 400 PE are
selected to test the AC model. The predicted number of AC events was 731.6,
while 733 events were observed in the data. The statistics in the AC model is
very large so the uncertainty of the predicted AC rate is neglected in the
following GOF tests. The 1D GOF tests in S1, S2, $\mathrm{Z}$, $\mathrm{R}$,
and the 2D test in S1–S2 all yield p-values between 0.05 and 0.95. The result
of the 2D test with a p-value of 0.61 is shown in LABEL:fig:ac_${}^{37}$Ar.
The model was further validated with events removed by selections targeting
the AC background inside the WIMP ROI. Similar tests were also performed on
220Rn calibration data. All these tests show the model and data to be
compatible.

To suppress AC events, a selection criterion based on a gradient-boosted
decision tree (GBDT) classifier utilizing the S2 pulse shape and the drift
time was developed. Due to an insufficient model for the S2 pulse shape near
the perpendicular wires, the GBDT classifier is only applied in the region at
least 4.45 cm far from the wires (far-wire region), while in the region near
the wires (near-wire region), a data-driven S2 shape selection is applied
instead. The resulting AC background rate per ton-year in the near-wire region
is about 5.6 times higher compared to the rate in the far-wire region.
Consequently, the modeling of the TPC response is split into two parts. Due to
an over-fitting issue found in the GBDT training process, we conservatively
estimated the rate uncertainty to be 30% for both the near- and far-wire
regions. The expectation value of AC events in the WIMP search dataset in the
ROI is $(4.3\pm 0.9)$.

Figure 8: AC model validation with 37Ar calibration data in S2 versus S1 area.
The AC model was cross-checked against this high-statistics dataset using a 2D
GOF test evaluated in an equiprobable binning scheme. The color scale shows
the deviation from the predicted number of events in each bin
$\mu_{\text{bin}}=14.0$ in units of sigma. The projections to S1 and S2 are
shown above and to the right of the 2D plot. No significant deviations were
observed in the bulk of the distribution.

###  III.6 Surface Background Model

Figure 9: The radial surface background model is shown for four slices in S2
area. For each slice, a Student’s generalized skew-t distribution is fitted to
210Po $\alpha$-events. The outer radius of the physical TPC as well as of the
fiducial volume are shown as a solid and a dashed black line, respectively.

In XENONnT, the TPC wall is made of PTFE, which is known to accumulate
isotopes from the decay chain of atmospheric 222Rn, which decays down to
210Pb, an isotope with a half-life of 22.2 years. Radioactive decays on the
wall surface can result in events with reduced S2 signals due to charge
losses, which gives rise to a background that can leak into the WIMP signal
region. Three components contribute to the surface background in the WIMP ROI:
the two $\beta$-decays of 210Pb and 210Bi, and the recoiling 206Pb following
the $\alpha$-decay of 210Po. Due to the complexity of electric field
conditions near the surface and the loss of ionization electrons to the
detector walls, we employ a data-driven approach to model this background
component in the space of cS1–cS2–$\mathrm{R}$–$\mathrm{Z}$.

The surface background model in the space of $\mathrm{R}$–$\mathrm{Z}$ was
developed using 210Po $\alpha$-events. These events are mono-energetic (5.4
MeV) and thus easily identifiable through their characteristic S1 signal, yet
they present a wide range of S2 sizes due to variable charge loss to the
walls. The $\mathrm{R}$ and $\mathrm{Z}$ distributions of 210Po
$\alpha$-events were seen to match those of the lower-energy $\beta$-events.
The radial profiles were modeled using a Student’s generalized skew-t
distribution. The radial distribution was fitted independently for different
S2 sizes as shown in Fig. 9 to account for the S2-size dependent position
resolution [10]. The $\mathrm{Z}$ profile was modeled using a linear function,
to account for lower rates of surface events observed at the bottom of the
TPC, likely due to the increased charged insensitivity near the walls at the
bottom of the detector [7].

The cS1–cS2 distribution was modeled using a 2D Gaussian adaptive kernel
density estimation (aKDE), built using events reconstructed outside the
detector walls. The resulting model was then validated against the events
reconstructed between the walls and the fiducial volume, in order to rule out
any dependence between cS1, cS2, and R. Figure 4 (lower left) shows the
projection of the four-dimensional model on cS1–cS2 in the ROI. The absolute
rate of events in the blinded region was inferred from the radial
distributions of adjacent, non-blinded regions in cS1–cS2.

Uncertainties in the model were obtained from the $\mathrm{R}$ and
$\mathrm{Z}$ fit parameter uncertainties, as those are the parameters of
interest in the development of the FV. Uncertainties on the overall measured
surface event rate were also propagated. Uncertainties in the (cS1, cS2) aKDE
were neglected, as toy-MC tests were performed to show that they had little
impact on overall expectation. For the FV, a maximum radius of
$61.35\text{\,}\mathrm{cm}$ compared to the $63\text{\,}\mathrm{cm}$ used in
[63] was chosen to reduce the number of surface events to $(14\pm 3)$, which
improves the robustness against mismodeling.

##  IV Statistical Inference

In this section, we describe the statistical methods used to derive the WIMP
search results, which generally follow the recommendations formulated in [46].
We first describe the likelihood function used for the search in Sec. IV.1 and
its nuisance parameters in Sec. IV.2, followed by an illustration of the
procedure for constructing confidence intervals and computing the discovery
significance in Sec. IV.3. The GOF test to validate the best-fit models and
the blinding procedure are described in Sec. IV.4. We omit some details
already presented in [9], as the methods here presented build upon that
previous work.

###  IV.1  Definition of the Likelihood Function

The fundamental ingredient for the statistical analysis of our WIMP search
data is the likelihood function $\mathcal{L}(\sigma,\boldsymbol{\theta})$. It
depends on the WIMP-nucleon cross section $\sigma\geq 0$, which is our
parameter of interest, and a set of nuisance parameters $\boldsymbol{\theta}$.
The likelihood function is defined as a product of four terms: two for the
WIMP search dataset ($\mathcal{L}_{\text{near-wire}}$ and
$\mathcal{L}_{\text{far-wire}}$), one for the 220Rn calibration dataset
($\mathcal{L}_{\text{cal}}$), and one for ancillary measurements of nuisance
parameters ($\mathcal{L}_{\mathrm{anc}}$):

$\displaystyle\mathcal{L}(\sigma,\boldsymbol{\theta})={}$
$\displaystyle\mathcal{L}_{\text{near-wire}}(\sigma,\boldsymbol{\theta})$ (5)
$\displaystyle\times\mathcal{L}_{\text{far-wire}}(\sigma,\boldsymbol{\theta})$
$\displaystyle\times\mathcal{L}_{\mathrm{cal}}(\boldsymbol{\theta})$
$\displaystyle\times\mathcal{L}_{\mathrm{anc}}(\boldsymbol{\theta}).$

In the WIMP search data, we categorize events based on their proximity to the
perpendicular wires of the gate and anode electrodes, distinguishing between
those that are near ($\leq$$4.45\text{\,}\mathrm{c}\mathrm{m}$, corresponding
to about $17\%$ of the total FV) and far from the wires in the
$\mathrm{X}$–$\mathrm{Y}$ plane. This approach allows us to account for the
higher AC rate observed near the wires as discussed in Sec. III.5, without
introducing a full position-dependence in the likelihood. Other backgrounds,
in particular radiogenic neutron and surface background, exhibit a substantial
radial dependence. Thus, the likelihood is modeled and evaluated in
$\mathrm{R}$, in addition to cS1 and cS2. For the near-wire region, the radial
component is not included in the modeling. Each WIMP search term is an
extended unbinned likelihood function of the form

$\displaystyle\mathcal{L}_{\text{region}}(\sigma,\boldsymbol{\theta})={}$
$\displaystyle\mathrm{Pois}(N|\mu_{\mathrm{tot}}(\sigma,\boldsymbol{\theta}))$
(6)
$\displaystyle\times\prod_{i=1}^{N}\left[\sum_{c}\frac{\mu_{c}(\sigma,\boldsymbol{\theta})}{\mu_{\mathrm{tot}}(\sigma,\boldsymbol{\theta})}\times
f_{c}(\vec{x}_{i}|\boldsymbol{\theta})\right],$

where the index “region” runs over “near-wire” and “far-wire”. The index $i$
runs over all $N$ observed events $\vec{x}_{i}$ in the WIMP search far-wire
(near-wire) dataset, where $\vec{x}_{i}$ is a vector with entries cS1, cS2,
and $\mathrm{R}$ (cS1 and cS2). The PDFs $f_{c}$ of each signal and background
component $c$ with expectation values $\mu_{c}$ are evaluated for each event.
The total expectation value is given by
$\mu_{\mathrm{tot}}(\sigma,\boldsymbol{\theta})=\sum_{c}\mu_{c}(\sigma,\boldsymbol{\theta})$.

The term $\mathcal{L}_{\mathrm{cal}}(\boldsymbol{\theta})$ is the extended
unbinned likelihood function of the 220Rn calibration dataset. It depends on
the two shape parameters introduced in Sec. III.2 that parameterize the range
of models consistent with the ER calibration data, selected using the PCA
method. By incorporating this likelihood term, we simultaneously fit the shape
parameters to the 220Rn calibration dataset and the ER background model in the
WIMP search dataset. Through this procedure, the constraint from the
calibration data on the shape of the ER model is propagated to the final
inference.

Finally, the ancillary likelihood function
$\mathcal{L}_{\mathrm{anc}}(\boldsymbol{\theta})$ is a product of constraint
terms for some of the nuisance parameters $\boldsymbol{\theta}$, which are
detailed in the following section.

###  IV.2 Nuisance Parameters and Constraints

Table 1: Parameters of the XENONnT SR0 WIMP search likelihood function and their constraints. The right column shows best-fit values with an unconstrained $200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ WIMP signal component (spin-independent coupling). For the rate parameters, all values are given in units of “events in the SR0 exposure”. Near and far-wire refers to the region in the $\mathrm{X}$–$\mathrm{Y}$ plane near and far from the perpendicular wires. In the near-wire region, which constitutes approximately $17\%$ of the total FV, the AC expectation value is comparable to the one in the far-wire region due to the higher AC rate near the wires. Rate Parameter | Constraint | Nominal | Best Fit
---|---|---|---
ER | WIMP search data | $134.5$ | $135_{-11}^{+12}$
Neutron | Ancillary measurement | $0.8\pm 0.4$ | $0.8\pm 0.4$
Neutron-X | Ancillary measurement | $0.31\pm 0.16$ | $0.30\pm 0.15$
CE$\nu$NS (Solar $\nu$) | Ancillary measurement | $0.19\pm 0.06$ | $0.19\pm 0.06$
CE$\nu$NS (Atm + DSN) | Ancillary measurement | $0.05\pm 0.02$ | $0.04\pm 0.02$
AC (near-wire) | Ancillary measurement | $2.3\pm 0.7$ | $2.3\pm 0.6$
AC (far-wire) | Ancillary measurement | $2.0\pm 0.6$ | $2.1\pm 0.6$
Surface | Ancillary measurement | $14\pm 3$ | $12\pm 2$
220Rn calibration |  220Rn dataset | $2062\pm 210$ | $2052\pm 44$
Shape Parameter | Constraint | Nominal | Best Fit
ER shape parameter 1 |  220Rn dataset | $0.0$ | $0.4\pm 0.2$
ER shape parameter 2 |  220Rn dataset | $0.0$ | $-1.8_{-0.6}^{+0.7}$
Signal Parameter | Constraint | Nominal | Best Fit
Signal efficiency | NR model uncertainty | $1.0$ | $1.0$
WIMP cross section [$10^{-47}\mathrm{cm^{-2}}$] |  | Parameter of interest | 3.22

In addition to our parameter of interest $\sigma$, we parameterize the space
of background hypotheses using nuisance parameters $\boldsymbol{\theta}$ in
the likelihood function of Eq. 5. These parameters control both background
rates and shape as well as the signal efficiency. In total, twelve nuisance
parameters are defined, which are listed in Tab. 1. The uncertainties of the
nominal expectation values correspond to the width of a Gaussian constraint
term. The ancillary likelihood function
$\mathcal{L}_{\mathrm{anc}}(\boldsymbol{\theta})$ is the product of all
constraint terms. Most parameters are constrained via ancillary measurements,
which are obtained from sidebands (e.g. data outside the ROI) or external
measurements in combination with simulations. More details on the constraints
on neutron, CE$\nu$NS, AC, and surface background rates were discussed in Sec.
III. The ER rate is a free parameter in the fit and is fully constrained by
the WIMP search data. The ER shape parameters obtained with PCA are
constrained via the simultaneous fit of the ${}^{220}\mathrm{Rn}$ calibration
dataset via the term $\mathcal{L}_{\mathrm{cal}}(\boldsymbol{\theta})$ in Eq.
5.

The best-fit values from the XENONnT SR0 WIMP search data [5] for a fit
including an unconstrained
$200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ signal component
are given in the last column of Tab. 1. The corresponding uncertainties
represent the two-sided one-sigma confidence intervals derived from a profile
likelihood scan in the respective parameter. More details on the construction
of confidence intervals are discussed in the following.

###  IV.3 Confidence Intervals and Discovery Significance

The construction of confidence intervals is based on the profile likelihood
test statistic

$q(\sigma)\equiv-2\ln\frac{\mathcal{L}(\sigma,\boldsymbol{\hat{\hat{\theta}}})}{\mathcal{L}(\hat{\sigma},{\boldsymbol{\hat{\theta}}})}.$
(7)

Quantities with a single hat denote the global maximum likelihood estimator of
a parameter, while $\boldsymbol{\hat{\hat{\theta}}}$ denotes the set of
nuisance parameters that maximize the likelihood if the cross section is fixed
to a given $\sigma$. The likelihood is defined for a specific signal model,
for example a $\mathrm{WIMP}$ with a certain mass and interaction type (e.g.
spin-independent or spin-dependent coupling). In our statistical inference, we
consider signal models across a range of WIMP masses, from
$6\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ to
$500\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$. For higher WIMP
masses, the PDF remains constant and the flux for a given cross section
decreases approximately linearly with mass, as discussed in Sec. III.1.
According to this, limits can be extrapolated even beyond
$500\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$.

Knowing the distribution of $q(\sigma)$ under different hypotheses is
essential for calculating discovery significances and confidence intervals.
Due to the low rate of background resembling the signal, asymptotic formulae
as discussed in [64] are not applicable. Instead, toy-MCs are used to estimate
the distribution $g(q(\sigma)|\sigma)$ of the test statistic given that
$\sigma$ is the true cross section. We use the custom-developed blueice [65]
framework to define the likelihood function. This Python package provides an
efficient interpolation (“template morphing”) between PDFs evaluated for
different discrete values of shape nuisance parameters. This allows shape
parameters to be considered for which we have no analytical description.

For each signal model considered, we compute the discovery significance and
the confidence interval by comparing the measured test statistic $q(\sigma)$
with the distribution under each hypothesis $g(q(\sigma)|\sigma)$ [66, 67].
Testing $\sigma=0$ yields the (local) discovery significance, while confidence
intervals are computed by finding the region where $\sigma$ is rejected at a
$90\%$ confidence level (CL) given the data, which is illustrated for three
background-only toy-MCs (three black parabola-like lines) in Fig. 10 (top).
The minimum of each curve corresponds to the respective best-fit value
$\hat{\sigma}$. Upper limits (ULs) and lower limits are obtained by finding
the cross sections at which $q(\sigma)$ reaches the critical region. The
threshold of the 90 $\%$ CL critical region is precomputed as the 90th
percentile of the test statistic distribution from toy-MCs with injected
signal corresponding to the respective cross section. For low cross sections,
it deviates from the asymptotic two-sided threshold indicated as a horizontal
gray dotted line.

Repeating this procedure $\mathcal{O}(1000)$ times yields the expected
distribution of ULs that can be obtained in the absence of any signal, shown
in Fig. 10 (bottom). The experiment’s sensitivity is given by the median of
these ULs and the sensitivity band (“Brazil band”) marked by the $\pm$2-sigma
and $\pm$1-sigma quantiles indicates the spread. The distribution of lower
limits for background-only toy-MCs is also shown in Fig. 10 (bottom). As
recommended in [46], we decide on a discovery significance threshold for
reporting two-sided confidence intervals of $3$ sigma, following [9].

Figure 10: Illustration of the confidence interval construction and
distribution of limits on the WIMP-nucleon cross section $\sigma$ of a
$200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ WIMP for
background-only toy-MCs. Top: The test statistic $q(\sigma)$ as a function of
the cross section $\sigma$ is shown for three toy-MCs. The intersections with
the threshold of the critical region (gray line) yield the 90 $\%$ CL upper
(blue diamonds) and lower limits (purple diamonds). Bottom: (Inverse)
cumulative distribution function (CDF) of upper (lower) limits for background-
only toy-MCs. The bands containing 68 $\%$ (green) and 95 $\%$ (yellow) of
ULs, as well as the median UL (dashed line), are indicated. Note the shared
abscissa with the cross section (bottom) and the corresponding expected number
of signal events (top).

Statistical fluctuations as well as mismodeling, such as overestimated
background rates, can yield arbitrarily low ULs, which may result in the
spurious exclusion of models beyond the experiment’s sensitivity. To mitigate
this issue, various methods have been proposed [68, 69, 70]. The XENONnT SR0
WIMP search follows the recommendation of [46] to use power-constrained limits
(PCL) [70]. In this method, a signal size threshold is selected, at which the
experiment has a “rejection power” $M_{\mathrm{min}}$ – the probability of
excluding a given signal under the background-only hypothesis. If an UL falls
below this threshold, the threshold value is reported instead. These
thresholds correspond to the quantiles of the UL distribution used to compute
the sensitivity band illustrated in Fig. 10. For instance, choosing
$M_{\mathrm{min}}=0.16$ sets the threshold to the -1 sigma quantile of the
band, while $M_{\mathrm{min}}=0.5$ truncates ULs at the band’s median. In
[46], the fiducial choice of $M_{\mathrm{min}}=0.16$ was suggested. However,
the choice was erroneously based on the discovery power, which corresponds to
the probability of rejecting the background-only hypothesis given an
alternative hypothesis with a specific signal size. In the absence of an
updated recommendation, we chose the very conservative choice of
$M_{\mathrm{min}}=0.5$ for the first WIMP search of XENONnT.

###  IV.4 Goodness of Fit and Blinding Procedure

To verify that the best-fit models describe our WIMP search data well, we
defined a GOF test before unblinding, and after studies to optimize the power
to reject mismodeling. The test uses a binned Poisson likelihood $\chi^{2}$
with 15 equiprobable bins in the cS1–cS2 space, defined from the model being
tested. To compute a p-value, the distribution of the $\chi^{2}$ test
statistic under the best-fit hypothesis is derived through toy-MCs. Specifying
the test and acceptance threshold (90% CL) before unblinding ensured that the
results were not influenced by statistical fluctuations expected from the low-
statistic WIMP search dataset. For the background-only fit we found a p-value
of $0.67$ and for a fit with an additional
$200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ WIMP signal the
p-value is $0.63$. Both indicate no strong mismodeling in the predefined
parameter space. Using GOF tests with well-studied power to discover signal-
like mismodeling was also designed to replace the “safeguard” ER shape
parameter defined in [71], due to its computational cost and susceptibility to
some kinds of mismodeling.

When searching for new signals in data, it is crucial to avoid the
experimenter’s bias on the result [72]. In our WIMP search, we adopt a common
strategy for bias mitigation: blinding the signal region until all steps of
the data analysis are finalized. Initially, events in the ER and NR bands with
a reconstructed ER energy below $20\,\mathrm{keV_{ER}}$ were blinded. This
involved defining the ER and NR bands in cS1 and cS2 based on quantiles in cS2
for slices in cS1 of calibration data. In the first step, all events above
$10\,\mathrm{keV_{ER}}$ and those above the $-2$ sigma quantile of the ER band
were unblinded for the analysis presented in [63]. The signal region for the
WIMP search (indicated in Fig. 5) remained blinded until the analysis
procedure was finalized. The unblinded SR0 WIMP search data showed no
significant excess for any of the tested WIMP masses with local discovery
p-values $\geq 0.2$. For this reason, we reported new ULs on the spin-
independent WIMP-nucleon cross section across WIMP masses ranging from
$6\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ to
$500\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$, with a minimum of
$2.58\times 10^{-47}\;\mathrm{cm}^{2}$ at
$28\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$.

##  V Summary

The XENONnT WIMP dark matter search relies on a detector response model as
well as simulation- and data-driven background models. These were combined to
construct a statistical model in cS1, cS2, and $\mathrm{R}$, which was used to
infer limits on WIMP-nucleus scattering cross sections.

The full detector response model for electronic recoil (ER) and nuclear recoil
(NR) interactions, including both the xenon emission and detector
reconstruction model, was successfully fitted to calibration data. Accurate
simulations of particle interactions up to the data acquisition waveform level
made this possible, in particular, to correctly model the S2 multiplicity of
events with several, potentially unresolved energy deposits. Using these
models, we derived the distributions for ER and NR backgrounds, as well as the
signals, in our analysis space. Except for ER, the background rates were
constrained with ancillary measurements. The radiogenic neutron background
rate was constrained by first matching the simulated ratio of multiple-to-
single scatter interactions and the neutron veto tagging efficiency with NR
calibration data. After unblinding the neutron veto-tagged events, these three
inputs were combined to derive a prediction for the remaining non-vetoed
single-scatter neutron background. Two shape parameters were propagated to the
final inference to account for the uncertainty of the ER model in cS1–cS2
space. Accidental coincidence and surface background were modeled with data-
driven approaches. The validity of the models was confirmed in calibration
data as well as science data outside the region of interest of the WIMP
search.

A blind analysis was performed for the first science data from XENONnT. The
statistical methods largely follow the previous XENON1T approach and community
recommendations, by using a toy-MC-calibrated profile log-likelihood ratio
test statistic. One departure from these recommendations was raising the
minimal required power of the power-constrained limit (PCL) threshold from
0.15 to 0.5, corresponding to placing upper limits only at or above the median
unconstrained upper limit. Both the calibration fits and the final best-fit
model were assessed and found acceptable with goodness-of-fit tests that were
chosen based on their mismodeling rejection power, and defined prior to
unblinding the data. Analysis of the data with an exposure of 1.1 tonne-years
revealed no signal excess over backgrounds. Therefore, new upper limits on the
spin-independent WIMP-nucleon cross section were derived.

###### Acknowledgements.

We gratefully acknowledge support from the National Science Foundation, Swiss
National Science Foundation, German Ministry for Education and Research, Max
Planck Gesellschaft, Deutsche Forschungsgemeinschaft, Helmholtz Association,
Dutch Research Council (NWO), Fundacao para a Ciencia e Tecnologia, Weizmann
Institute of Science, Binational Science Foundation, Région des Pays de la
Loire, Knut and Alice Wallenberg Foundation, Kavli Foundation, JSPS Kakenhi
and JST FOREST Program ERAN in Japan, Tsinghua University Initiative
Scientific Research Program, DIM-ACAV+ Région Ile-de-France, and Istituto
Nazionale di Fisica Nucleare. This project has received funding/support from
the European Union’s Horizon 2020 research and innovation program under the
Marie Skłodowska-Curie grant agreement No 860881-HIDDeN. We gratefully
acknowledge support for providing computing and data-processing resources of
the Open Science Pool and the European Grid Initiative, in the following
computing centers: the CNRS/IN2P3 (Lyon - France), the Dutch national
e-infrastructure with the support of SURF Cooperative, the Nikhef Data-
Processing Facility (Amsterdam - Netherlands), the INFN-CNAF (Bologna -
Italy), the San Diego Supercomputer Center (San Diego - USA) and the Enrico
Fermi Institute (Chicago - USA). We acknowledge the support of the Research
Computing Center (RCC) at The University of Chicago for providing computing
resources for data analysis. We are grateful to Laboratori Nazionali del Gran
Sasso for hosting and supporting the XENON project.

## References

  * Bertone and Hooper [2018] G. Bertone and D. Hooper, History of dark matter, Rev. Mod. Phys. 90, 045002 (2018).
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck collaboration), Planck 2018 results: Vi. cosmological parameters, A&A 641, A6 (2020).
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Particle dark matter: evidence, candidates and constraints, Phys. Rept. 405, 279 (2005).
  * Roszkowski  _et al._ [2018] L. Roszkowski, E. M. Sessolo, and S. Trojanowski, WIMP dark matter candidates and searches—current status and future prospects, Rep. Prog. Phys. 81, 066201 (2018).
  * E. Aprile  _et al._ [2023a] E. Aprile _et al._ (XENON collaboration), First dark matter search with nuclear recoils from the XENONnT experiment, Phys. Rev. Lett. 131, 041003 (2023a).
  * E. Aprile  _et al._ [2024a] E. Aprile _et al._ (XENON collaboration), The XENONnT Dark Matter Experiment, preprint (2024a), arXiv:2402.10446 [physics.ins-det] .
  * Aprile  _et al._ [2024] E. Aprile _et al._ (XENON collaboration), Design and performance of the field cage for the XENONnT experiment, Eur. Phys. J. C 84, 138 (2024).
  * E. Aprile  _et al._ [2019a] E. Aprile _et al._ (XENON collaboration), XENON1T dark matter data analysis: Signal reconstruction, calibration, and event selection, Phys. Rev. D 100, 052014 (2019a).
  * E. Aprile  _et al._ [2019b] E. Aprile _et al._ (XENON collaboration), XENON1T dark matter data analysis: Signal and background models and statistical inference, Phys. Rev. D 99, 112009 (2019b).
  * E. Aprile  _et al._ [2024b] E. Aprile _et al._ (XENON collaboration), XENONnT Analysis: Signal Reconstruction, Calibration and Event Selection; In preparation (2024b).
  * Ramírez García [2022] D. Ramírez García, Ph.D. thesis, Albert-Ludwigs-Universität Freiburg (2022).
  * Zhu [2022] T. Zhu, Ph.D. thesis, Columbia University (2022).
  * Szydagis  _et al._ [2011] M. Szydagis, N. Barry, K. Kazkaz, J. Mock, D. Stolp, M. Sweany, M. Tripathi, S. Uvarov, N. Walsh, and M. Woods, NEST: a comprehensive model for scintillation yield in liquid xenon, JINST 6, P10002 (2011).
  * Lenardo  _et al._ [2015] B. Lenardo, K. Kazkaz, A. Manalaysay, J. Mock, M. Szydagis, and M. Tripathi, A global analysis of light and charge yields in liquid xenon, IEEE TNS 62, 3387 (2015).
  * Thomas and D. A. Imel [1987] J. Thomas and D. A. Imel, Recombination of electron-ion pairs in liquid argon and liquid xenon, Phys. Rev. A 36, 614 (1987).
  * López Paredes  _et al._ [2018] B. López Paredes, H. Araújo, F. Froborg, N. Marangou, I. Olcina, T. Sumner, R. Taylor, A. Tomás, and A. Vacheret, Response of photomultiplier tubes to xenon scintillation light, Astropart. Phys. 102, 56 (2018).
  * Faham  _et al._ [2015] C. Faham, V. Gehman, A. Currie, A. Dobi, P. Sorensen, and R. Gaitskell, Measurements of wavelength-dependent double photoelectron emission from single photons in VUV-sensitive photomultiplier tubes, JINST 10, P09010 (2015).
  * Jörg  _et al._ [2023] F. Jörg, S. Li, J. Schreiner, H. Simgen, and R. F. Lang, Characterization of a ${}^{{\textbf{220}}}$Rn source for low-energy electronic recoil calibration of the XENONnT detector, JINST 18, P11009 (2023).
  * E. Aprile  _et al._ [2023b] E. Aprile _et al._ (XENON collaboration), Low-energy calibration of XENON1T with an internal ${}^{{\textbf{37}}}$Ar source, Eur. Phys. J. C 83, 1 (2023b).
  * D. Wenz [2023] D. Wenz, Ph.D. thesis, Johannes Gutenberg-Universität Mainz (2023).
  * E. Aprile  _et al._ [2024c] E. Aprile _et al._ (XENON collaboration), Erratum to XENON1T Dark Matter Data Analysis: Signal & Background Models, and Statistical Inference, Phys. Rev. D (2024c), in preparation.
  * Baudis  _et al._ [2021] L. Baudis, P. Sanchez-Lucas, and K. Thieme, A measurement of the mean electronic excitation energy of liquid xenon, Eur. Phys. J. C 81, 1060 (2021).
  * Boulton  _et al._ [2017] E. M. Boulton _et al._ , Calibration of a two-phase xenon time projection chamber with a 37Ar source, JINST 12, P08004 (2017).
  * D. S. Akerib  _et al._ [2017] D. S. Akerib _et al._ (LUX collaboration), Ultralow energy calibration of LUX detector using 127Xe electron capture, Phys. Rev. D 96, 112011 (2017).
  * D. S. Akerib  _et al._ [2016a] D. S. Akerib _et al._ (LUX collaboration), Tritium calibration of the LUX dark matter experiment, Phys. Rev. D 93, 072009 (2016a).
  * Aprile  _et al._ [2005] E. Aprile, K. L. Giboni, P. Majewski, K. Ni, M. Yamashita, R. Hasty, A. Manzur, and D. N. McKinsey, Scintillation response of liquid xenon to low energy nuclear recoils, Phys. Rev. D 72, 072006 (2005).
  * E. Aprile  _et al._ [2006] E. Aprile _et al._ (XENON collaboration), Simultaneous measurement of ionization and scintillation from nuclear recoils in liquid xenon as target for a dark matter experiment, Phys. Rev. Lett. 97, 081302 (2006).
  * E. Aprile  _et al._ [2009] E. Aprile _et al._ (XENON collaboration), New Measurement of the Relative Scintillation Efficiency of Xenon Nuclear Recoils Below 10 keV, Phys. Rev. C 79, 045807 (2009).
  * E. Aprile  _et al._ [2013] E. Aprile _et al._ (XENON collaboration), Response of the XENON100 Dark Matter Detector to Nuclear Recoils, Phys. Rev. D 88, 012006 (2013).
  * Plante  _et al._ [2011] G. Plante, E. Aprile, R. Budnik, B. Choi, K.-L. Giboni, L. W. Goetzke, R. F. Lang, K. E. Lim, and A. J. Melgarejo Fernandez, New Measurement of the Scintillation Efficiency of Low-Energy Nuclear Recoils in Liquid Xenon, Phys. Rev. C 84, 045805 (2011).
  * P. Sorensen  _et al._ [2009] P. Sorensen _et al._ (XENON collaboration), The scintillation and ionization yield of liquid xenon for nuclear recoils, Nucl. Instrum. Meth. A 601, 339 (2009).
  * Manzur  _et al._ [2010] A. Manzur, A. Curioni, L. Kastens, D. N. McKinsey, K. Ni, and T. Wongjirad, Scintillation efficiency and ionization yield of liquid xenon for mono-energetic nuclear recoils down to 4 keV, Phys. Rev. C 81, 025808 (2010).
  * D. S. Akerib  _et al._ [2016b] D. S. Akerib _et al._ (LUX collaboration), Low-energy (0.7-74 keV) nuclear recoil calibration of the LUX dark matter experiment using D-D neutron scattering kinematics, preprint (2016b), arXiv:1608.05381 [physics.ins-det] .
  * D. S. Akerib  _et al._ [2022] D. S. Akerib _et al._ (LUX collaboration), Improved Dark Matter Search Sensitivity Resulting from LUX Low-Energy Nuclear Recoil Calibration, preprint (2022), arXiv:2210.05859 [physics.ins-det] .
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC hammer, PASP 125, 306 (2013).
  * Foreman-Mackey, Daniel and others [2019] Foreman-Mackey, Daniel and others, emcee v3: A Python ensemble sampling toolkit for affine-invariant MCMC, J. Open Source Softw. 4, 1864 (2019), arXiv:1911.07688 [astro-ph.IM] .
  * XENON collaboration [2022] XENON collaboration, XENONnT/GOFevaluation: v0.1.1 (2022).
  * Lewin and Smith [1996a] J. Lewin and P. Smith, Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil, Astropart. Phys. 6, 87 (1996a).
  * R. H. Helm [1956] R. H. Helm, Inelastic and elastic scattering of 187-MeV electrons from selected even-even nuclei, Phys. Rev. 104, 1466 (1956).
  * Lewin and Smith [1996b] J. D. Lewin and P. F. Smith, Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil, Astropart. Phys. 6, 87 (1996b).
  * Smith  _et al._ [2007] M. C. Smith _et al._ , The RAVE Survey: Constraining the Local Galactic Escape Speed, Mon. Not. Roy. Astron. Soc. 379, 755 (2007).
  * McCabe [2014] C. McCabe, The Earth’s velocity for direct detection experiments, JCAP 02, 027 (2014).
  * Schoenrich  _et al._ [2010] R. Schoenrich, J. Binney, and W. Dehnen, Local Kinematics and the Local Standard of Rest, Mon. Not. Roy. Astron. Soc. 403, 1829 (2010).
  * Bland-Hawthorn and Gerhard [2016] J. Bland-Hawthorn and O. Gerhard, The galaxy in context: Structural, kinematic, and integrated properties, Annu. Rev. Astron. Astrophys. 54, 529 (2016).
  * Abuter, R.  _et al._ [2021] Abuter, R. _et al._ (GRAVITY Collaboration), Improved gravity astrometric accuracy from modeling optical aberrations, A&A 647, A59 (2021).
  * Baxter  _et al._ [2021] D. Baxter _et al._ , Recommended conventions for reporting results from direct dark matter searches, Eur. Phys. J. C 81, 907 (2021).
  * Klos  _et al._ [2013] P. Klos, J. Menéndez, D. Gazit, and A. Schwenk, Large-scale nuclear structure calculations for spin-dependent wimp scattering with chiral effective field theory currents, Phys. Rev. D 88, 083516 (2013).
  * Pearson [1901] K. Pearson, LIII. On lines and planes of closest fit to systems of points in space, London Edinburgh Dublin Philos. Mag. & J. Sci. 2, 559 (1901).
  * Aprile  _et al._ [2020] E. Aprile _et al._ (XENON collaboration), Projected WIMP sensitivity of the XENONnT dark matter experiment, JCAP 11, P031 (2020).
  * Madland  _et al._ [1999] D. G. Madland _et al._ , SOURCES 4A: A code for calculating (alpha,n), spontaneous fission, and delayed neutron sources and spectra 10.2172/15215 (1999).
  * E. Aprile  _et al._ [2017] E. Aprile _et al._ (XENON collaboration), Material radioassay and selection for the XENON1T dark matter experiment, Eur. Phys. J. C 77, 890 (2017).
  * E. Aprile  _et al._ [2022a] E. Aprile _et al._ (XENON collaboration), Material radiopurity control in the XENONnT experiment, Eur. Phys. J. C 82, 599 (2022a).
  * S. Agostinelli  _et al._ [2003] S. Agostinelli _et al._ (Geant4 collaboration), Geant4: A Simulation toolkit, Nucl. Instrum. Meth. A 506, 250 (2003).
  * J. Allison  _et al._ [2016] J. Allison _et al._ (Geant4 collaboration), Recent developments in Geant4, Nucl. Instrum. Meth. A 835, 186 (2016).
  * XENON collaboration [2023] XENON collaboration, XENONnT/epix: v0.3.4 (2023).
  * XENON collaboration [2022a] XENON collaboration, XENONnT/WFSim: v1.0.2 (2022a).
  * XENON collaboration [2022b] XENON collaboration, XENONnT/straxen: v1.8.3 (2022b).
  * E. Aprile  _et al._ [2021] E. Aprile _et al._ (XENON collaboration), Search for Coherent Elastic Scattering of Solar 8B Neutrinos in the XENON1T Dark Matter Experiment, Phys. Rev. Lett. 126, 091301 (2021).
  * B. Aharmim  _et al._ [2013] B. Aharmim _et al._ (SNO collaboration), Combined Analysis of all Three Phases of Solar Neutrino Data from the Sudbury Neutrino Observatory, Phys. Rev. C 88, 025501 (2013).
  * M. Agostini  _et al._ [2018] M. Agostini _et al._ (BOREXINO collaboration), Comprehensive measurement of $pp$-chain solar neutrinos, Nature 562, 505 (2018).
  * E. Aprile  _et al._ [2018] E. Aprile _et al._ (XENON collaboration), Dark Matter Search Results from a One Tonne$\times$Year Exposure of XENON1T, Phys. Rev. Lett. 121, 111302 (2018).
  * Aalbers  _et al._ [2022] J. Aalbers _et al._ , AxFoundation/strax: v1.2.3 (2022).
  * E. Aprile  _et al._ [2022b] E. Aprile _et al._ (XENON collaboration), Search for new physics in electronic recoil data from XENONnT, Phys. Rev. Lett. 129, 161805 (2022b).
  * Cowan  _et al._ [2011a] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Asymptotic formulae for likelihood-based tests of new physics, Eur. Phys. J. C 71, 1554 (2011a).
  * Aalbers  _et al._ [2021] J. Aalbers, K. D. Morå, and B. Pelssers, JelleAalbers/blueice: v1.1.0 (2021).
  * G. J. Feldman and R. D. Cousins [1998] G. J. Feldman and R. D. Cousins, Unified approach to the classical statistical analysis of small signals, Phys. Rev. D 57, 3873 (1998).
  * Patrignani [2016] C. Patrignani (Particle Data Group collaboration), Review of Particle Physics, Chin. Phys. C40, 100001 (2016).
  * A. L. Read [2002] A. L. Read, Presentation of search results: the cls technique, J. Phys. G: Nucl. Part. Phys. 28, 2693 (2002).
  * Kashyap  _et al._ [2010] V. L. Kashyap, D. A. van Dyk, A. Connors, P. E. Freeman, A. Siemiginowska, J. Xu, and A. Zezas, On computing upper limits to source intensities, ApJ 719, 900 (2010).
  * Cowan  _et al._ [2011b] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Power-constrained limits, preprint (2011b), arXiv:1105.3166 [physics.data-an] .
  * Priel  _et al._ [2017] N. Priel, L. Rauch, H. Landsman, A. Manfredini, and R. Budnik, A model independent safeguard against background mismodeling for statistical inference, JCAP 2017, 013 (2017).
  * J. R. Klein and Roodman [2005] J. R. Klein and A. Roodman, Blind analysis in nuclear and particle physics, Annu. Rev. Nucl. Part. Sci. 55, 141 (2005).
  * Aprile and Doke [2010] E. Aprile and T. Doke, Liquid Xenon Detectors for Particle Physics and Astrophysics, Rev. Mod. Phys. 82, 2053 (2010).
  * Doke  _et al._ [1976] T. Doke, A. Hitachi, S. Kubota, A. Nakamoto, and T. Takahashi, Estimation of fano factors in liquid argon, krypton, xenon and xenon-doped liquid argon, Nucl. Instr. and Meth. 134, 353 (1976).
  * Lindhard  _et al._ [1963] J. Lindhard, M. Scharff, and H. E. Schioett, Range concepts and heavy ion ranges (notes on atomic collisions, II), Kgl. Danske Videnskab. Selskab. Mat. Fys. Medd. 33,  (1963).
  * Thomas and Imel [1987] J. Thomas and D. A. Imel, Recombination of electron-ion pairs in liquid argon and liquid xenon, Phys. Rev. A 36, 614 (1987).
  * D. M. Mei  _et al._ [2008] D. M. Mei, Z. B. Yin, L. C. Stonehill, and A. Hime, A model of nuclear recoil scintillation efficiency in noble liquids, Astropart. Phys. 30, 12 (2008).

##  Appendix A Details on the ER Emission Model

Let $E$ be the recoil energy of an ER event. The total number of detectable
quanta $N_{\mathrm{q}}$ is sampled from a normal distribution

$N_{\mathrm{q}}\sim\mathrm{Norm}(E/W,\sqrt{fE/W}).$ (8)

$W$ is the mean energy to generate a quantum and $f$ is the Fano factor
$\leavevmode\nobreak\ 0.059$ [73, 74]. The ER energy deposit will produce both
excitons and electron-ion pairs. Using the mean ratio between number of
excitons and ions $\left\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\right\rangle$,
we simulate their numbers from a binomial distribution:

$N_{\mathrm{i}}\sim\mathrm{Binom}\left(N_{\mathrm{q}},\frac{1}{1+\langle
N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle}\right),$ (9)
$N_{\mathrm{ex}}=N_{\mathrm{q}}-N_{\mathrm{i}}.$ (10)

A fraction $r$ of ions recombine with electrons, depending on the electric
drift field $F$. We parametrize the mean value of $r$ in the same way as in
XENON1T [9]:

$\langle
r\rangle=\frac{1}{\mathrm{e}^{-(E-q_{0})/q_{1}}+1}\left(1-\frac{\log(1+\langle
N_{\mathrm{i}}\rangle\varsigma)}{\langle
N_{\mathrm{i}}\rangle\varsigma}\right),$ (11)

where

$\varsigma=\frac{1}{4}\gamma\mathrm{e}^{-E/\omega}F^{-\delta}.$ (12)

The fluctuation of $r$ is parametrized via

$\Delta r=q_{2}(1-\mathrm{e}^{-E/q_{3}}),$ (13)

and the true recombination fraction $r$ is then sampled from

$r\sim\mathrm{Norm}(\langle r\rangle,\Delta r).$ (14)

All fitted ER parameters are listed in Tab. 2. Finally, the numbers of
produced electrons and photons are given by

$N_{\mathrm{e}}\sim\mathrm{Binom}(N_{\mathrm{i}},1-r),$ (15)
$N_{\gamma}=N_{\mathrm{q}}-N_{\mathrm{e}}.$ (16)

##  Appendix B Details on the NR Emission Model

Let $E$ be the recoil energy of an NR event. The total number of produced
quanta $N$ is

$N\sim\mathrm{Norm}(E/W,\sqrt{fE/W}).$ (17)

In contrast to electron recoils, recoiling xenon nuclei lose a significant
amount of their recoil energy via atomic motion and collisions with other
xenon atoms, which are processes that are undetectable in a LXe TPC. This
quanta loss is modeled following the Lindhard theory [75], with the so-called
Lindhard quenching factor $L$, such that the number of detectable quanta
becomes

$N_{\mathrm{q}}\sim\mathrm{Binom}(N,L).$ (18)

Following the parametrization in [14], the Lindhard factor is given by

$L=\frac{\kappa\,g(\epsilon)}{1+\kappa\,g(\epsilon)},$ (19)

where $g(\epsilon)$ is a function of deposited energy via

$g(\epsilon)=3\epsilon^{0.15}+0.7\epsilon^{0.6}+\epsilon,$ (20)
$\epsilon=11.5\,Z^{-7/3}(E/\mathrm{keV}),$ (21)

with the nuclear charge number $Z=54$ for xenon. The numbers of produced ions
and excitons are then simulated by

$N_{\mathrm{i}}\sim\mathrm{Binom}\left(N_{\mathrm{q}},\frac{1}{1+\langle
N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle}\right),$ (22)
$N_{\mathrm{ex}}=N_{\mathrm{q}}-N_{\mathrm{i}},$ (23)

with the exciton-to-ion ratio parametrized as

$\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle=\alpha
F^{-\zeta}(1-\mathrm{e}^{-\beta\epsilon}).$ (24)

Unlike ERs, the recombination fluctuation in NRs is usually small, thus the
number of photons produced from the recombination of electron-ion pairs is

$N_{\gamma}^{\mathrm{re}}\sim\mathrm{Binom}(N_{\mathrm{i}},r),$ (25)

where $r$ follows the Thomas-Imel box model [76]

$r=1-\frac{\log(1+N_{i}\varsigma)}{N_{i}\varsigma},$ (26) $\varsigma=\gamma
F^{-\delta}.$ (27)

The number of excitons is further reduced by bi-excitonic and Penning
quenching effects [77], such that the number of photons produced from de-
excitation becomes

$N_{\gamma}^{\mathrm{de}}\sim\mathrm{Binom}(N_{\mathrm{ex}},f_{l}),$ (28)

with the scintillation quenching factor

$f_{l}=\frac{1}{1+\eta\epsilon^{\lambda}}.$ (29)

Then the final numbers of photons and electrons are

$N_{\gamma}=N_{\gamma}^{\mathrm{de}}+N_{\gamma}^{\mathrm{re}},$ (30)
$N_{\mathrm{e}}=N_{\mathrm{i}}-N_{\gamma}^{\mathrm{re}}.$ (31)

##  Appendix C Details on the Detector Response Model

The simulations of $\mathrm{S1}$ and $\mathrm{S2}$ signals from $N_{\gamma}$
and $N_{\mathrm{e}}$ are almost independent of each other.

To simulate $\mathrm{S1}$ signals, we first introduce the spatially dependent
scintillation gain $\tilde{g_{1}}$ with DPE excluded:

$\tilde{g_{1}}(\mathrm{X},\mathrm{Y},\mathrm{Z})=\frac{g_{1}}{1+p_{\mathrm{DPE}}}\cdot\mathrm{S1corr}^{-1}(\mathrm{X},\mathrm{Y},\mathrm{Z}),$
(32)

where $g_{1}$ is the averaged scintillation gain, $\mathrm{S1corr}$ is the
relative spatial correction factor, and $p_{\mathrm{DPE}}$ is the probability
of the double photoelectron emission effect [17]. Then, for a given number of
photons $N_{\gamma}$ and the position of the recoil $\mathrm{X}$,
$\mathrm{Y}$, $\mathrm{Z}$, the number of photons detected by the PMTs is

$N_{\mathrm{PhD}}\sim\mathrm{Binom}(N_{\gamma},\tilde{g_{1}}(\mathrm{X},\mathrm{Y},\mathrm{Z})).$
(33)

Now accounting for the DPE effect, the expected number of photoelectrons (PE)
in the S1 signal is simulated via

$N_{\mathrm{DPE}}\sim\mathrm{Binom}(N_{\mathrm{PhD}},p_{\mathrm{DPE}}),$ (34)
$N_{\mathrm{S1,PE}}\sim N_{\mathrm{PhD}}+N_{\mathrm{DPE}}.$ (35)

However, the $\mathrm{S1}$ is not always equal to the number of PEs due to
reconstruction bias. This is modeled as

$\mathrm{S1}/N_{\mathrm{S1,PE}}-1\sim\mathrm{Norm}(\delta_{\mathrm{S1}},\Delta_{\mathrm{S1}}).$
(36)

Here, both $\delta_{\mathrm{S1}}$, $\Delta_{\mathrm{S1}}$ are functions of
$N_{\mathrm{PhD}}$ obtained from the XENONnT MC. In the end, the relative
spatial correction is applied back to the $\mathrm{S1}$ signal, yielding

$\mathrm{cS1}=\mathrm{S1}\cdot\mathrm{S1corr}(\mathrm{X}^{\prime},\mathrm{Y}^{\prime},\mathrm{Z}),$
(37)

where $\mathrm{X}^{\prime},\mathrm{Y}^{\prime}$ are the event positions as
reconstructed in data, which are smeared by the S2-size-dependent resolution
of the position reconstruction.

To get the $\mathrm{S2}$ signal from $N_{\mathrm{e}}$, the first step is to
simulate the electron loss during the drift. The survival probability due to
attachment to electronegative impurities in LXe is given by

$p_{\mathrm{loss}}(\mathrm{Z})=\mathrm{e}^{-\mathrm{Z}/(\tau v)},$ (38)

where $\tau$ is the electron drift survival time (“electron lifetime”) and $v$
is the drift velocity of electrons in LXe. Inside the FV the electric field is
uniform enough to approximate $v$ as a constant. Electron losses due to drift
field effects close to the wall are accounted for via a spatially dependent
charge-insensitive-volume (CIV) probability function
$p_{\mathrm{CIV}}(\mathrm{R},\mathrm{Z})$ from field simulations of XENONnT
[7]. The number of surviving electrons is then given by

$N_{\mathrm{surv}}\sim\mathrm{Binom}\left(N_{\mathrm{e}},\>p_{\mathrm{loss}}(\mathrm{Z})\cdot
p_{\mathrm{CIV}}(\mathrm{R},\mathrm{Z})\right).$ (39)

At the liquid-gas interface, a fraction of electrons are extracted,

$N_{\mathrm{extr}}\sim\mathrm{Binom}(N_{\mathrm{surv}},\epsilon(\mathrm{X},\mathrm{Y})).$
(40)

The extraction efficiency is $\mathrm{X}$-$\mathrm{Y}$ dependent and can be
calculated by

$\epsilon(\mathrm{X},\mathrm{Y})=g_{2}\cdot\mathrm{S2corr}^{-1}(\mathrm{X},\mathrm{Y})/G,$
(41)

where $g_{2}$ is the averaged ionization gain, $\mathrm{S2corr}$ is the
relative spatial correction factor, and $G\sim 31\,\mathrm{PE/e^{-}}$ is the
single electron gain. Assuming the fluctuation of the secondary scintillation
process is Poisson-like, the expected number of $\mathrm{S2}$ PEs is

$N_{\mathrm{S2,PE}}\sim\mathrm{Norm}(N_{\mathrm{extr}}G,\sqrt{N_{\mathrm{extr}}G}).$
(42)

Similar to $\mathrm{S1}$, the $\mathrm{S2}$ reconstruction bias is modeled by

$\mathrm{S2}/N_{\mathrm{S2,PE}}-1\sim\mathrm{Norm}(\delta_{\mathrm{S2}},\Delta_{\mathrm{S2}}).$
(43)

Finally, the correction is applied to the $\mathrm{S2}$,

$\mathrm{cS2}=\mathrm{S2}\cdot\mathrm{S2corr}(\mathrm{X}^{\prime},\mathrm{Y}^{\prime})\cdot\mathrm{e}^{\mathrm{Z}/\tau^{\prime}v}.$
(44)

Here, $\tau^{\prime}$ is the electron lifetime used in the correction, which
could be slightly off from the true value. The difference is one of the
parameters to fit in the band fitting.

For both S1 and S2 signals, data quality selection efficiencies are applied
which depend on the respective signal sizes. For S1 signals, we additionally
apply the S1 reconstruction efficiency which is a function of the number of
detected photons $N_{\mathrm{PhD}}$.

For NR multi-scatter events, the simulation is slightly altered. Here the
total number of detected photons is defined as the sum of the detected photons
for each single energy deposition $i$,

$N_{\mathrm{PhD}}=\sum_{i}N_{\mathrm{PhD}}^{(i)}.$ (45)

Similarly, we sum over the number of S2 PEs of only those scatters that
contribute to the S2,

$N_{\mathrm{S2,PE}}=\sum_{i_{\mathrm{S2}}}N_{\mathrm{S2,PE}}^{(i_{\mathrm{S2}})}.$
(46)

Reconstruction biases, signal corrections, and selection efficiencies are
correspondingly applied to the merged signals.

##  Appendix D Details on the Calibration Fit Results

The parameters in the ER and NR band fits are summarized in Tab. 2, showing
both the applied prior and the marginal posterior values. Parameters where
prior values are given without uncertainties are fixed in the fit. The prior
for $\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle$ is a uniform distribution
in the stated range. For the posteriors, the stated values and uncertainties
correspond to the median and the central 68% of the marginal posterior
distributions.

Table 2: Prior and marginal posterior distributions of each parameter in the
ER and NR emission models.

Parameter Prior Marginal posterior Unit $W$ $13.7\pm 0.2$ $13.7^{+0.2}_{-0.2}$
eV $f$ 0.059 0.059 - ER parameters $\langle
N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle$ $0.06$ \- $0.20$ $0.13^{+0.04}_{-0.04}$
- $\gamma$ free $0.13^{+0.03}_{-0.02}$ - $\delta$ free $0.34^{+0.07}_{-0.07}$
- $\omega$ free $57^{+15}_{-12}$ keV $q_{0}$ free $1.32^{+0.17}_{-0.20}$ keV
$q_{1}$ free $0.47^{+0.07}_{-0.05}$ keV $q_{2}$ free $0.030^{+0.002}_{-0.002}$
- $q_{3}$ free $0.47^{+0.40}_{-0.31}$ keV NR parameters $\zeta$ 0.047 0.047 -
$\delta$ 0.062 0.062 - $\alpha$ free $0.92^{+0.07}_{-0.06}$ - $\gamma$ free
$0.016^{+0.001}_{-0.001}$ - $\beta$ $239^{+56}_{-18}$ $334^{+40}_{-43}$ -
$\kappa$ $0.139^{+0.006}_{-0.005}$ $0.138^{+0.005}_{-0.006}$ - $\eta$
$3.3^{+10.6}_{-1.4}$ $10.0^{+6.8}_{-5.9}$ - $\lambda$ $1.14^{+0.90}_{-0.18}$
$1.40^{+0.61}_{-0.38}$ -