# Dense Nuclear Matter Equation of State from Heavy-Ion Collisions

Agnieszka Sorensen<EMAIL_ADDRESS>Kshitij Agarwal  Kyle W.
Brown  Zbigniew Chajęcki  Paweł Danielewicz  Christian Drischler  Stefano
Gandolfi  Jeremy W. Holt  Matthias Kaminski  Che-Ming Ko  Rohit Kumar  Bao-An
Li  William G. Lynch  Alan B. McIntosh  William G. Newton  Scott Pratt  Oleh
Savchuk  Maria Stefaniak  Ingo Tews  ManYee Betty Tsang  Ramona Vogt  Hermann
Wolter  Hanna Zbroszczyk
Endorsing authors: Navid Abbasi  Jörg Aichelin  Anton Andronic  Steffen A.
Bass  Francesco Becattini  David Blaschke  Marcus Bleicher  Christoph Blume
Elena Bratkovskaya  B. Alex Brown  David A. Brown  Alberto Camaiani  Giovanni
Casini  Katerina Chatziioannou  Abdelouahad Chbihi  Maria Colonna  Mircea Dan
Cozma  Veronica Dexheimer  Xin Dong  Travis Dore  Lipei Du  José A. Dueñas
Hannah Elfner  Wojciech Florkowski  Yuki Fujimoto  Richard J. Furnstahl
Alexandra Gade  Tetyana Galatyuk  Charles Gale  Frank Geurts  Sašo Grozdanov
Kris Hagel  Steven P. Harris  Wick Haxton  Ulrich Heinz  Michal P. Heller  Or
Hen  Heiko Hergert  Norbert Herrmann  Huan Zhong Huang  Xu-Guang Huang
Natsumi Ikeno  Gabriele Inghirami  Jakub Jankowski  Jiangyong Jia  José C.
Jiménez  Joseph Kapusta  Behruz Kardan  Iurii Karpenko  Declan Keane  Dmitri
Kharzeev  Andrej Kugler  Arnaud Le Fèvre  Dean Lee  Hong Liu  Michael A. Lisa
William J. Llope  Ivano Lombardo  Manuel Lorenz  Tommaso Marchi  Larry
McLerran  Ulrich Mosel  Anton Motornenko  Berndt Müller  Paolo Napolitani
Joseph B. Natowitz  Witold Nazarewicz  Jorge Noronha  Jacquelyn Noronha-
Hostler  Grażyna Odyniec  Panagiota Papakonstantinou  Zuzana Paulínyová  Jorge
Piekarewicz  Robert D. Pisarski  Christopher Plumberg  Madappa Prakash  Jørgen
Randrup  Claudia Ratti  Peter Rau  Sanjay Reddy  Hans-Rudolf Schmidt  Paolo
Russotto  Radoslaw Ryblewski  Andreas Schäfer  Björn Schenke  Srimoyee Sen
Peter Senger  Richard Seto  Chun Shen  Bradley Sherrill  Mayank Singh
Vladimir Skokov  Michał Spaliński  Jan Steinheimer  Mikhail Stephanov  Joachim
Stroth  Christian Sturm  Kai-Jia Sun  Aihong Tang  Giorgio Torrieri  Wolfgang
Trautmann  Giuseppe Verde  Volodymyr Vovchenko  Ryoichi Wada  Fuqiang Wang
Gang Wang  Klaus Werner  Nu Xu  Zhangbu Xu  Ho-Ung Yee  Sherry Yennello  Yi
Yin  Institute for Nuclear Theory, University of Washington, Seattle, WA
98195, USA  Physikalisches Institut, Eberhard Karls Universität Tübingen,
D-72076 Tübingen, Germany  Facility for Rare Isotope Beams, Michigan State
University, East Lansing, MI 48824, USA  Department of Chemistry, Michigan
State University, East Lansing, MI 48824, USA  Department of Physics, Western
Michigan University, Kalamazoo, MI 49008, USA  Department of Physics and
Astronomy, Michigan State University, East Lansing, Michigan 48824, USA
Department of Physics and Astronomy, Ohio University, Athens, OH 45701, USA
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545,
USA  Department of Physics and Astronomy, Texas A&M University, College
Station, TX 77843, USA  Cyclotron Institute, Texas A&M University, College
Station, TX 77843, USA  Department of Physics and Astronomy, University of
Alabama, Tuscaloosa, AL 35487, USA  Department of Physics and Astronomy, Texas
A&M University-Commerce, Commerce, TX 75429, USA  Bogolyubov Institute for
Theoretical Physics, 03680 Kyiv, Ukraine  Department of Physics, Ohio State
University, Columbus, OH 43210, USA  GSI Helmholtz Centre for Heavy-ion
Research, 64291 Darmstadt, Germany  Nuclear and Chemical Sciences Division,
Lawrence Livermore National Laboratory, Livermore, CA 94551, USA  Department
of Physics and Astronomy, University of California Davis, Davis, CA 95616, USA
Faculty of Physics, University of Munich, D-85748 Garching, Germany  Faculty
of Physics, Warsaw University of Technology, 00-662 Warsaw, Poland  School of
Nuclear Science and Technology, Lanzhou University, Lanzhou 730000, China
SUBATECH, Nantes University – IN2P3/CNRS – IMT Atlantique, Nantes, France
Frankfurt Institute for Advanced Studies, D-60438 Frankfurt am Main, Germany
Institut für Kernphysik, Westfälische Wilhelms-Universität Münster, 48149
Münster, Germany  Department of Physics, Duke University, Durham, NC 27708,
USA  Dipartimento di Fisica e Astronomia, Università di Firenze, Sesto
Fiorentino I-50019, Firenze, Italy  INFN Sezione di Firenze, 50019 Sesto
Fiorentino, Firenze, Italy  Institute of Theoretical Physics, University of
Wrocław, 50-204 Wrocław, Poland  Center for Advanced Systems Understanding
(CASUS), D-02826 Görlitz, Germany  Helmholtz-Zentrum Dresden-Rossendorf
(HZDR), D-01328 Dresden, Germany  Institut für Theoretische Physik, Goethe
Universität Frankfurt, D-60438 Frankfurt am Main, Germany  Helmholtz Research
Academy Hesse for FAIR (HFHF), GSI Helmholtz Center, Campus Frankfurt, 60438
Frankfurt am Main, Germany  Institut für Kernphysik, Goethe Universität
Frankfurt, D-60438 Frankfurt am Main, Germany  Nuclear Science and Technology
Department, Brookhaven National Laboratory, Upton, NY 11973, USA  Department
of Physics, California Institute of Technology, Pasadena, CA 91125, USA  LIGO
Laboratory, California Institute of Technology, Pasadena, CA 91125, USA
GANIL, CEA/DRF-CNRS/IN2P3, F-14076 Caen Cedex, France  INFN-Laboratori
Nazionali del Sud, I-95123, Catania, Italy  Horia Hulubei National Institute
of Physics and Nuclear Engineering (IFIN-HH), 077125 Mǎgurele-Bucharest,
Romania  Department of Physics, Kent State University, Kent, OH 44242 USA
Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA
94720, USA  Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld,
Germany  Department of Physics, McGill University, Montreal, Quebec H3A 2T8,
Canada  Departamento de Ingeniería Eléctrica y Centro de Estudios Avanzados en
Física, Matemáticas y Computación, Universidad de Huelva, 21007 Huelva, Spain
Institute of Theoretical Physics, Jagiellonian University, 30-348 Kraków,
Poland  Department of Physics, Technische Universität Darmstadt, 64289
Darmstadt, Germany  Department of Physics and Astronomy, Rice University,
Houston TX 77005, USA  Higgs Centre for Theoretical Physics, University of
Edinburgh, Edinburgh, EH8 9YL, Scotland  Faculty of Mathematics and Physics,
University of Ljubljana, SI-1000 Ljubljana, Slovenia  Department of Physics,
University of California Berkeley, Berkeley, CA 94720, USA  Department of
Physics and Astronomy, Ghent University, 9000 Ghent, Belgium  Department of
Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Physikalisches Institut, Universität Heidelberg, Heidelberg, Germany
Department of Physics and Astronomy, University of California Los Angeles, Los
Angeles, CA 90095, USA  Physics Department & Center for Particle Physics and
Field Theory, Fudan University, Shanghai 200438, China  Key Laboratory of
Nuclear Physics and Ion-Beam Application (MOE), Fudan University, Shanghai
200433, China  Shanghai Research Center for Theoretical Nuclear Physics
(NSFC), Fudan University, Shanghai 200438, China  Department of Agricultural,
Life and Environmental Sciences, Tottori University, Tottori 680-8551, Japan
Department of Chemistry, Stony Brook University, Stony Brook, NY 11794, USA
Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA
Instituto de Física, Universidade de São Paulo, 05508–090 São Paulo-SP, Brazil
School of Physics and Astronomy, University of Minnesota, Minneapolis, MI
55455, USA  Faculty of Nuclear Sciences and Physical Engineering, Czech
Technical University in Prague, 11519 Prague 1, Czech Republic  Center for
Nuclear Theory, Department of Physics and Astronomy, Stony Brook University,
Stony Brook, NY 11794, USA  Nuclear Physics Institute, The Czech Academy of
Sciences, 25068 Rez, Czech Republic  Center for Theoretical Physics,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Department of
Physics and Astronomy, Wayne State University, Detroit, MI 48201, USA
Department of Physics and Astronomy, University of Catania, I-95123 Catania,
Italy  INFN, Laboratori Nazionali di Legnaro, I-35020 Legnaro, Italy  Institut
für Theoretische Physik, Universität Giessen, Giessen, Germany  Helmholtz
Research Academy Hesse for FAIR (HFHF), Campus Giessen, Giessen, Germany
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  Illinois
Center for Advanced Studies of the Universe & Department of Physics,
University of Illinois Urbana-Champaign, Urbana, IL 61801, USA  Rare Isotope
Science Project, Institute for Basic Science, Daejeon 34000, Korea
Prírodovedecká fakulta, Univerzita Pavla Jozefa Šafárika, Košice, Slovakia
Department of Physics, Florida State University, Tallahassee, FL 32306, USA
Natural Science Division, Pepperdine University, Malibu, CA 90263, USA
Physics Department, University of Houston, Box 351550, Houston, TX 77204, USA
Institute of Nuclear Physics Polish Academy of Sciences, PL-31-342 Kraków,
Poland  Institute for Theoretical Physics, Regensburg University, D-93040
Regensburg, Germany  Department of Physics and Astronomy, Iowa State
University, Ames, IA 50011, USA  Facility for Antiproton and Ion Research,
Darmstadt, Germany  Department of Physics and Astronomy, University of
California Riverside, Riverside, CA 92521, USA  RIKEN BNL Research Center,
Brookhaven National Laboratory, Upton, NY 11973, USA  Department of Physics,
North Carolina State University, Raleigh, NC 27695, USA  National Centre for
Nuclear Research, 02-093 Warsaw, Poland  Physics Department, University of
Białystok, 15-245 Białystok, Poland  Department of Physics, University of
Illinois, Chicago, IL 60607, USA  Institute of Modern Physics, Fudan
University, 200438, Shanghai, China  Instituto de Fisica Gleb Wataghin -
UNICAMP Universidade Estadual de Campinas, 13083-859, Campinas SP, Brazil
Institute of Physics, Jan Kochanowski University, 25-406 Kielce, Poland  INFN
Sezione di Catania, I-95123 Catania, Italy  Department of Physics and
Astronomy, Purdue University, West Lafayette, IN 47907, USA  Chemistry
Department, Texas A&M University, College Station, TX 77843, USA  Quark Matter
Research Center, Institute of Modern Physics, Chinese Academy of Sciences,
Lanzhou, Gansu, 073000, China

###### Abstract

The nuclear equation of state (EOS) is at the center of numerous theoretical
and experimental efforts in nuclear physics. With advances in microscopic
theories for nuclear interactions, the availability of experiments probing
nuclear matter under conditions not reached before, endeavors to develop
sophisticated and reliable transport simulations to interpret these
experiments, and the advent of multi-messenger astronomy, the next decade will
bring new opportunities for determining the nuclear matter EOS, elucidating
its dependence on density, temperature, and isospin asymmetry. Among
controlled terrestrial experiments, collisions of heavy nuclei at intermediate
beam energies (from a few tens of MeV/nucleon to about 25 GeV/nucleon in the
fixed-target frame) probe the widest ranges of baryon density and temperature,
enabling studies of nuclear matter from a few tenths to about 5 times the
nuclear saturation density and for temperatures from a few to well above a
hundred MeV, respectively. Collisions of neutron-rich isotopes further bring
the opportunity to probe effects due to the isospin asymmetry. However,
capitalizing on the enormous scientific effort aimed at uncovering the dense
nuclear matter EOS, both at RHIC and at FRIB as well as at other international
facilities, depends on the continued development of state-of-the-art hadronic
transport simulations. This white paper highlights the essential role that
heavy-ion collision experiments and hadronic transport simulations play in
understanding strong interactions in dense nuclear matter, with an emphasis on
how these efforts can be used together with microscopic approaches and neutron
star studies to uncover the nuclear EOS.

###### keywords:

heavy-ion collisions, hadronic transport, nuclear matter, equation of state,
symmetry energy

††journal: Progress in Particle and Nuclear Physics

###### Contents

  1. 1 Introduction
    1. 1.1 Constraining the nuclear matter EOS using heavy-ion collisions
    2. 1.2 Connections to fundamental questions in nuclear physics
    3. 1.3 Upcoming opportunities
    4. 1.4 Scientific needs
  2. 2 The equation of state from $0$ to $5n_{0}$
    1. 2.1 Transport model simulations of heavy-ion collisions
      1. 2.1.1 Transport theory
      2. 2.1.2 Selected constraints on the EOS obtained from heavy-ion collisions
      3. 2.1.3 Challenges and opportunities
    2. 2.2 Microscopic calculations of the EOS
      1. 2.2.1 Status
      2. 2.2.2 Challenges and opportunities
    3. 2.3 Neutron star theory
      1. 2.3.1 Status
      2. 2.3.2 Challenges and opportunities
  3. 3 Heavy-ion collision experiments
    1. 3.1 Experiments to extract the EOS of symmetric nuclear matter
      1. 3.1.1 Measurements sensitive to the EOS
      2. 3.1.2 Experiments probing densities between 1–2.5$n_{0}$
      3. 3.1.3 Experiments probing densities above 2.5$n_{0}$
      4. 3.1.4 Challenges and opportunities
    2. 3.2 Experiments to extract the symmetry energy
      1. 3.2.1 Experiments that probe low densities
      2. 3.2.2 Measurements to extract symmetry energy up to 1.5$n_{0}$
      3. 3.2.3 Selected constraints on the symmetry energy around 1.5$n_{0}$
      4. 3.2.4 Challenges and opportunities
  4. 4 The equation of state from combined constraints
    1. 4.1 Constraints
    2. 4.2 EOS obtained by combining various constraint sets
  5. 5 Connections to other areas of nuclear physics
    1. 5.1 Applications of hadronic transport
      1. 5.1.1 Detector design
      2. 5.1.2 Space exploration, radiation therapy, and nuclear data
    2. 5.2 Hydrodynamics
      1. 5.2.1 Status
      2. 5.2.2 Range of applicability
      3. 5.2.3 Challenges and opportunities
  6. 6 Exploratory directions
    1. 6.1  Dense nuclear matter EOS meeting extreme gravity and dark matter in supermassive neutron stars
    2. 6.2 Nuclear EOS with reduced spatial dimensions
    3. 6.3 Interplay between nucleonic and partonic degrees of freedom: SRC effects on nuclear EOS, heavy-ion reactions, and neutron stars
    4. 6.4 High-density symmetry energy above $2n_{0}$
    5. 6.5 Density-dependence of neutron-proton effective mass splitting in neutron-rich matter
  7. A Executive Summary

## 1 Introduction

The equation of state (EOS) is a fundamental property of nuclear matter,
describing its emergent macroscopic behavior originating from the underlying
strong interactions. For ordinary nuclear matter, the EOS controls the
structure of nuclei through the binding energy [1], the incompressibility [2,
3], or the neutron-skin thickness in neutron-rich nuclei [4, 5]. The EOS also
determines the properties of nuclear matter at extreme densities and/or
temperatures, corresponding to conditions produced in experiments colliding
heavy nuclei [6, 7, 8, 9, 10] or observed in neutron stars [11] and neutron
star mergers [12]. Far beyond describing the properties of matter composed of
only protons and neutrons, the EOS can also reflect the appearance of new
degrees of freedom, e.g., strange particles in the cores of neutron stars [13,
14, 15, 16, 17] or quarks and gluons in ultrarelativistic heavy-ion collisions
[18, 19, 20, 21], or the emergence of new states of matter, e.g., chirally-
restored matter [22, 23, 24], meson condensates [25, 26, 27], or quarkyonic
matter [28, 29, 30].

Figure 1.1: Schematic depiction of the ranges of density and temperature
probed in experiments and astronomical observations sensitive to the EOS of
nuclear matter (counterclockwise from bottom left): neutron star (NS) crust
physics, including nuclear pasta structures; properties of nuclei; structure
of neutron stars; dynamics of neutron star mergers; and outcomes of heavy-ion
collisions (HIC) which can probe both symmetric and asymmetric matter. Figures
adapted from [31, 32, 33, 34, 35].

In heavy-ion collision experiments, the EOS is studied by detecting particles
emerging from the collision zone and measuring observables sensitive to the
properties of nuclear matter. Interpretation of these observables, including
quantitative constraints on the EOS, requires comparisons of experimentally
measured observables to results obtained in dynamic simulations. This White
Paper highlights the essential role of hadronic transport simulations of
heavy-ion collisions in advancing our understanding of the EOS. It also
elucidates the many connections between inferences of the EOS from heavy-ion
collision data and other efforts aiming to describe and understand the
properties of nuclear matter.

### 1.1 Constraining the nuclear matter EOS using heavy-ion collisions

Figure 1.2: Constraints on the zeroth ($S_{v}$) and first ($L$) coefficient of
the symmetry energy expansion. Experimental constraints are derived from
heavy-ion collisions (HIC) [36], neutron-skin thicknesses of Sn isotopes [37],
giant dipole resonances (GDR) [38], the dipole polarizability of 208Pb [39,
40], nuclear masses [41], and isovector skins (IAS+$\Delta R$) [42]. Also
shown are constraints from $\chi$EFT (GP-B) [43], microscopic neutron-matter
calculations (H, G) [44, 45], and from the unitary gas limit (UG) [46]. Figure
from [43].

The last decade has brought tremendous progress in extracting the EOS as a
function of baryon density $n_{B}$, temperature $T$, and the isospin asymmetry
$\delta$ (or, equivalently, the proton fraction) from a variety of
experimental and astronomical data as well as theoretical calculations. Many-
body theory, based on sophisticated approaches with input from nucleon
scattering or nuclear structure data, can now state the EOS below and near the
saturation density $n_{0}$ with meaningful uncertainties [47, 48, 49, 50] (see
Section 2.2, “Microscopic calculations of the EOS”). New classes of
experiments have extracted the thickness of neutron skins in nuclei [51, 52,
5, 53, 54, 55], shedding light on the isospin-dependence of the EOS (or,
equivalently, the symmetry energy) near or below $n_{0}$. High-energy heavy-
ion collisions [56, 57, 58, 59] have constrained the EOS of the quark-gluon
plasma at high temperatures and small baryon densities [60], while ongoing
experimental efforts worldwide focus on the EOS of nearly-symmetric dense
baryonic matter, probed in collisions at intermediate energies. Meanwhile,
collisions at lower energies have led to experimental constraints on the
symmetry energy at sub- and suprasaturation densities [61, 62, 63, 64, 65].
Most remarkably, a revolution in the quality and breadth of astronomical
observations, highlighted by the first simultaneous detection of gravitational
waves and electromagnetic signals from a neutron-star merger [66], ushered in
a new era of multi-messenger astronomy (see Section 2.3, “Neutron star
theory”). Together with the newly available experimental capabilities at the
Facility for Rare Isotope Beams (FRIB), there are unprecedented opportunities
to probe the isospin-dependence of the EOS through astronomical and
terrestrial measurements.

Among the experimental efforts discussed above, heavy-ion collisions probe the
widest range of baryon densities and, moreover, represent the only means to
address the EOS away from $n_{0}$ in controlled terrestrial experiments, see
Fig. 1.1. Indeed, heavy-ion reactions at beam kinetic energies from a few tens
of MeV/nucleon to about 25 GeV/nucleon in the fixed-target frame probe the EOS
of hadronic matter at baryon densities from a few tenths to about 5 times
$n_{0}$. Controlling the properties of matter produced in these experiments is
possible by varying the beam energy, collision geometry, and isotopic
composition of the target and projectile. Insights and constraints obtained
from transport model analyses of these experiments are relevant both for our
understanding of nuclear matter as found on Earth and for our understanding of
neutron stars from crust to core.

Within ongoing efforts, the STAR experiment’s Beam Energy Scan (BES) fixed-
target (FXT) program [7] at the Relativistic Heavy Ion Collider (RHIC) at the
Brookhaven National Laboratory (BNL), which collided gold nuclei at
intermediate beam energies and which completed data taking in 2022, leads the
U.S. efforts to constrain the EOS of

Figure 1.3: Pressure in neutron star matter as a function of density from a
Bayesian analysis combining nuclear theory and data from multi-messenger
neutron-star observations and heavy-ion collisions [67]: the dark blue and
light blue region corresponds to the 68% and 95% credible interval,
respectively, while the gray dashed line shows the 95% bound obtained in
$\chi$EFT calculations and used as a prior. Figure from [67].

nearly-symmetric nuclear matter at high baryon densities up to around
$5n_{0}$, corresponding to densities present in the deep interiors of neutron
stars. Among comparable efforts in Europe, the HADES experiment [68] at GSI,
Germany, probes matter at densities up to $2.5n_{0}$. Preliminary results from
these contemporary efforts, as well as measurements from other heavy-ion
collision experiments in the past, have led to competitive constraints on the
EOS of symmetric nuclear matter [69, 70, 71, 72, 73], with future measurements
expected to shed more light on its high-density behavior. Detailed constraints
on the isospin-dependence of the EOS can be obtained by varying the isospin
content of the target and projectile nuclei. Here, the ability to use
radioactive isotopes, as in, e.g., intermediate-energy heavy-ion collision
experiments at RIKEN and FRIB, is crucial to resolve the subtle effects
arising from changes in the isospin asymmetry of colliding systems [9].

Above all, obtaining constraints on the EOS from heavy-ion measurements would
not have been possible if not for advances in theory, and in particular for
the collaborative effort to test the robustness and quantify the uncertainties
of hadronic transport simulations (see Section 2.1, “Transport model
simulations of heavy-ion collisions”). At the same time, much remains to be
learned, as tight constraints on both the symmetric and asymmetric EOS at
higher densities have so far remained elusive. This is predominantly due to
model uncertainties [74], which themselves are rooted in the inherent
complexity of nucleus-nucleus collisions and the challenging task of
describing all processes contributing to the final state observables.

### 1.2 Connections to fundamental questions in nuclear physics

The wealth of data from efforts conducted in recent years has brought forward
fascinating questions challenging our understanding of strong interactions.

Following the successful BES-I campaign at RHIC, questions remain about the
structure of the QCD phase diagram at large baryon densities, where the sign
problem prevents predictions from lattice QCD calculations [75] and
extrapolations of lattice QCD results become unreliable [21]. Surprisingly,
the expected disappearance of quark-gluon plasma signatures has not been
unequivocally observed in BES-I, with some observables suggesting that the QCD
first-order phase transition may be located in the statistically demanding
low-energy region probed by BES-II [76], including the region probed by the
currently analyzed BES-II FXT data [77, 78]. If this is the case, then
constraining the EOS at lower densities and describing the approach to the
transition from the hadronic side, which would manifest as a softening of the
EOS, will be crucial for a robust interpretation of BES-II measurements. Due
to the largely out-of-equilibrium evolution of collision systems probing that
region of the QCD phase diagram, hadronic transport simulations will play a
dominant role in describing the dynamics of the collisions, and therefore in
constraining the EOS of nearly-symmetric dense nuclear matter.

Understanding the physics of neutron-rich matter across a range of densities
is necessary not only to explain the properties of rare neutron-rich isotopes
and the structure of neutron stars, but also to constrain microscopic
interactions in isospin-asymmetric nuclear matter. At low densities, this
challenge is addressed by experimental and theoretical analyses

Figure 1.4: Compilation of the world data on the slope of the directed flow at
mid-rapidity ($dv_{1}/dy|_{y^{\prime}=0}$, top) and the elliptic flow
($v_{2}$, bottom) as functions of the reduced center-of-mass energy
$\sqrt{s_{NN}}-2m_{N}$ for protons, $Z=1$ nuclei, and inclusive charged
particles. Figure modified from [79].

of nuclear structure observables [80, 81, 82]. An important objective of
nuclear many-body theorists is to accurately calculate these observables and
reliably deduce the EOS using microscopic interactions derived within the
framework of chiral effective field theory ($\chi$EFT) [83, 84, 85]. Probing
the symmetry energy over a range of densities wider than found in nuclei is
possible through heavy-ion collisions and neutron star studies. Often,
knowledge of the isospin-asymmetric EOS is encoded in terms of constraints on
the Taylor expansion coefficients of the symmetry energy around $n_{0}$.
Numerous analyses yield consistent constraints on the first few expansion
coefficients [43] (see Fig. 1.2), although they rely on an assumption that the
expansion remains accurate away from $n_{0}$. The recent advent of Bayesian
inference techniques allows one to pursue a different approach, within which
the isospin-asymmetric EOS is described in terms of the dependence of the
pressure on baryon density [67] (see Fig. 1.3). Moreover, Bayesian analyses
can shed more light on densities at which measurements constrain the symmetry
energy and quantify the uncertainties of the extracted EOS. As a result,
combining diverse measurements and using advanced analysis techniques can lead
to significantly tighter constraints, especially on the high-density behavior
of the symmetry energy (or, equivalently, on the higher-order symmetry energy
expansion coefficients), so far poorly known.

Constraints on the EOS of neutron-rich matter at high densities have been
dramatically affected by discoveries of heavy neutron stars [86, 87]. Combined
with the properties of all known compact stars, these observations indicate
that while the EOS of neutron-rich matter is relatively soft around
(1–2)$n_{0}$, the pressure steeply rises with density for $n_{B}\gtrsim
2n_{0}$ [88, 89]. In fact, multiple analyses show that describing the known
population of neutron stars is only possible for EOSs in which the speed of
sound in neutron-star matter breaks the conformal limit at high densities,
that is exceeds $1/\sqrt{3}$ of the speed of light $c$ for $n_{B}\gtrsim
2n_{0}$ [90, 91, 92, 93, 94]. This striking behavior remains to be understood.
In particular, it is currently not known whether the speed of sound exceeds
$c/\sqrt{3}$ above certain densities at all isospin fractions of nuclear
matter or, alternatively, only in neutron-rich matter. Robust constraints on
the symmetric matter EOS at $n_{B}\gtrsim 2n_{0}$, obtained from heavy-ion
collisions at intermediate to high beam energies, would answer this question
as well as put constraints on the isospin-dependent part of the EOS through
comparisons with the EOS inferred from neutron star studies, thus uncovering
the magnitude of isospin-related effects at high baryon density.

.

### 1.3 Upcoming opportunities

The next decade will be an era of high-luminosity heavy-ion collision
experiments at high baryon density with modern detector and analysis
procedures, as well as detailed studies of the symmetry energy with collisions
of proton- and neutron-rich isotopes.

Many of the discoveries of the BES program in ultra-relativistic heavy-ion
collisions at RHIC, e.g., the discovery of the triangular flow [95, 96] and
elliptic flow fluctuations [97], illustrate that modern analyses of heavy-ion
collisions bring new quality to the understanding of the underlying processes.
Because of this, revisiting the intermediate to high beam energies, previously
explored at the AGS at BNL as well as at SIS18 at GSI and now explored by the
STAR FXT program and the HADES experiment, is imperative to enable putting
tighter constraints on the EOS of dense nuclear matter. Moreover, the future
CBM experiment [10, 98] at the Facility for Antiproton and Ion Research
(FAIR), Germany, will be able to measure interaction rates exceeding those
currently used by several orders of magnitude, allowing for exploration of
multiple high-statistics observables [99]. Furthermore, the explored beam
energy range is where lower-order flow observables, reflecting the collective
motion of the colliding system due to the underlying hadronic EOS, are
particularly

Figure 1.5: Predicted slope of the directed flow at mid-rapidity
($dv_{1}/dy|_{y^{\prime}=0}$, top) and elliptic flow ($v_{2}$, bottom) as
functions of the incompressibility and the in-medium nucleon-nucleon
scattering cross section modification factor, generated in simulations of
Au+Au reactions using the isospin-dependent BUU (IBUU) transport model [100,
101]. Figure from Ref. [102].

prominent (see Fig. 1.4). Therefore, the corresponding precision measurements
carry with them the opportunity to bring a richer perspective and a better
understanding of the physics underlying the complex dynamics of nuclear matter
at extreme conditions (see Section 3.1, “Experiments to extract the EOS of
symmetric nuclear matter”). This advancement can only occur provided a
simultaneous development of hadronic transport simulations, as only a detailed
understanding of various factors affecting the dynamics of heavy-ion
collisions can lead to meaningful descriptions of the experimental data, and,
consequently, more robust constraints on the EOS of nearly-symmetric nuclear
matter (see Section 2.1, “Transport model simulations of heavy-ion
collisions”). As an example of the sensitivity of observables to various
details of the underlying physics, Fig. 1.5 shows the dependence of the slope
of the directed flow (top panel) and of the elliptic flow at midrapidity
(bottom panel) on the stiffness of the EOS, parametrized by the
incompressibility, and on the in-medium nucleon-nucleon scattering cross-
section modification factor.

Unprecedented possibilities are on the horizon for studies of the isospin-
dependence of the EOS, which is critical for connecting heavy-ion collision
measurements to astrophysical observations. The difficulties in using nuclei
with significant variations in the isospin asymmetry, along with the paucity
of neutron measurements at midrapidity, have in the past greatly restricted
the capability to put tight constrains on the EOS of asymmetric nuclear
matter. Fortunately, at this time modern neutron detectors are available for
heavy-ion measurements in many facilities, including at accelerators
performing collisions at high beam energies such as GSI, while radioactive
beam measurements are entering a new era at RIKEN and FRIB. FRIB will provide
proton- and neutron-rich beams of not only the highest-intensity worldwide,
but also characterized by the widest currently accessible range of the isospin
asymmetry [103, 9]. Establishing a strong heavy-ion program at FRIB will
therefore enable previously inaccessible exploration of the symmetry energy
(see Section 3.2, “Experiments to extract the symmetry energy”). Moreover, the
proposed FRIB400 beam energy upgrade would not only allow exploration of
densities up to around $2n_{0}$, but it would also provide increased
resolution of the isospin-dependence of the EOS [9]. In particular, among
observables sensitive to the symmetry energy [104, 105], both charged pion
yields and the absolute magnitude of the elliptic flow (see Fig. 1.4)
significantly increase between the current top FRIB beam kinetic energy of 200
MeV/nucleon in the fixed-target frame and the proposed 400 MeV/nucleon [9].

The increase in available computing power and advances in statistical methods
make it possible to perform wide-ranging comparisons of heavy-ion collision
simulations with experimental data (e.g., using Bayesian analysis), allowing
one to vary multiple model assumptions at the same time as well as to put
robust uncertainties on the obtained constraints. Furthermore, given the
wealth of the upcoming independent data, e.g., from heavy-ion collision
experiments, neutron star observations, and microscopic nuclear theory
calculations, global analyses of complementary efforts have likewise a strong
potential for putting tight constraints on the EOS (see Section 4, “The EOS
from combined constraints”).

Beyond the much-needed interpretation of intermediate energy heavy-ion
collisions, advances in transport theory can lead to significant contributions
to other areas of nuclear physics. Recently, attention has been given to
cross-cutting opportunities for employing state-of-the-art hadronic transport
codes in studies supporting space exploration and advanced medical treatments
(see Section 5.1, “Applications of hadronic transport”). Transport theories
may also be used in tests of extensions of hydrodynamic approaches supporting
far-from-equilibrium evolution (see Section 5.2, “Hydrodynamics”), which are a
focus of intense studies due to their importance for modeling heavy-ion
collisions at high energies. Finally, constraining the dense nuclear matter
EOS through interpretations of heavy-ion collision measurements may have other
profound consequences, including helping to answer fundamental questions about
the possible existence of dark matter in the cores of neutron stars or
providing the impetus for studies of nuclear systems in fractional dimensions
(see Section 6, “Exploratory directions”).

### 1.4 Scientific needs

The next-generation experimental measurements of observables sensitive to the
nuclear matter EOS are imminent. Capitalizing on the enormous worldwide
scientific effort aimed at uncovering the dense nuclear matter EOS through
heavy-ion collision experiments is contingent on enhanced theory support. In
particular, the development of transport theories based on microscopic
hadronic degrees of freedom, which are the only means of interpreting
measurements from heavy-ion collisions at intermediate beam energies, must be
strengthened and expanded. Support for individual scientists, and in
particular creating opportunities for early career researchers, is imperative
to maintain the health of and diversify the U.S. hadronic transport community.
Collaborative research programs are needed to enable both a systematic
advancement of hadronic transport models as well as explorations of new
directions in microscopic descriptions of heavy-ion collisions. With this
support, innovative theory research will enable the exploration of the dense
nuclear matter EOS and help fully realize the potential of the U.S. nuclear
physics facilities.

## 2 The equation of state from $0$ to $5n_{0}$

Efforts to determine the equation of state (EOS) of nuclear matter are at the
forefront of nuclear physics. An EOS contains fundamental information about
the properties of a many-body system (see, e.g., Section 1.2), and is, in
essence, any nontrivial relation between the thermodynamic properties of a
given type of matter. In nuclear physics, the form of the EOS that is most
often pursued is the relation between energy per baryon or pressure and baryon
density $n_{B}$, isospin asymmetry $\delta$, and temperature $T$. For
symmetric matter, the isospin excess vanishes ($\delta=0$), and for asymmetric
matter the energy per baryon or pressure are commonly partitioned into a part
corresponding to symmetric matter and the remainder, which contains all
information about the isospin-dependence of the EOS. Due to the charge
invariance of strong interactions, the latter part is (to a good accuracy)
quadratic in $\delta$ at densities around the nuclear saturation density
$n_{0}$, relevant to nuclear experiments and astrophysical observations. The
quadratic coefficients in the expansion around $\delta=0$ are independent of
$\delta$, and are often referred to as the symmetry energy (denoted as
$S(n_{B})$ at $T=0$) or symmetry pressure, respectively. These, together with
the EOS of symmetric matter, are then sufficient to describe the EOS of
nuclear matter at any isospin asymmetry.

While many approaches to constraining the nuclear matter EOS are pursued, here
we describe three research areas which have the capability to constrain the
EOS over wide ranges of density: inferences of the EOS from comparisons of
experimental measurements to model simulations of heavy-ion collisions
(Section 2.1), microscopic calculations of the EOS using chiral effective
field theory (Section 2.2), and EOS inferences from neutron star studies
(Section 2.3). To fully utilize the opportunity behind these complementary
research directions, efforts must be made, both within these areas as well as
across the different approaches, to describe the underlying physics
consistently and to combine different data sets in a well-controlled way. This
is expanded upon in Section 4.

Given the wealth of data expected in the near future from heavy-ion collision
experiments, nuclear structure studies, and multi-messenger astronomical
observations, a concerted theoretical effort aimed at a consistent
interpretation of precise measurements across varying thermodynamic conditions
is needed. Some efforts of this nature are already underway. The Nuclear
Physics from Multi-Messenger Mergers (NP3M) NSF Focused Research Hub [106]
develops theoretical models and numerical simulations of dense and hot matter
to connect multi-messenger observations of neutron stars to the underlying
merger dynamics. The goals of the Network for Neutrinos, Nuclear Astrophysics,
and Symmetries (N3AS) NSF Physics Frontier Center [107] include developing
simulations of supernovae, mergers, kilonovae, and cooling neutron stars
incorporating the most advanced treatments of the underlying neutrino and
nuclear matter microphysics, enabling robust connections between observations
and fundamental neutrino, dense matter, and dark matter properties. The
Modular Unified Solver of the Equation of State (MUSES) NSF CSSI Framework
[108] is developing a cyberinfrastructure with novel tools to answer critical
interdisciplinary questions in nuclear physics, and will provide, e.g.,
modules for generating ensembles of EOSs relevant for heavy-ion collisions and
neutron stars or neutron star mergers. Finally, the Beam Energy Scan Theory
(BEST) DOE Topical Collaboration in Nuclear Theory [109], operating in the
years 2016–2020, supported the development of hadronic transport models for
the description of the final state of heavy-ion collision experiments
performed within the Beam Energy Scan (BES) program at the Relativistic Heavy
Ion Collider (RHIC) [110].

At this time, results from the fixed-target (FXT) campaign of the BES are
imminent, while future experiments at intermediate beam energies are being
planned worldwide. In this regime, hadronic transport simulations are
currently the only way to interpret these measurements and use them to
understand the properties of QCD interactions at high baryon densities and
finite temperatures. Consequently, there is an urgent need for a collaborative
theoretical research program aimed at a further development of hadronic
transport models as well as at explorations of new directions in microscopic
descriptions of heavy-ion collisions. Possible research directions that could
be addressed within such an effort are described in Section 2.1.3.

### 2.1 Transport model simulations of heavy-ion collisions

Heavy-ion collisions at low to intermediate beam energies provide the means to
probe nuclear matter at different baryon densities (from subsaturation to
several times the saturation density), temperatures (from a few MeV to well
above one hundred), and neutron to proton ratios (from nearly symmetric
nuclear matter, where $N_{n}/N_{p}\approx 1$ and $\delta\approx 0$, to very
neutron-rich matter, where $N_{n}/N_{p}\approx 2$ and $\delta\approx 0.25$).
An illustrative calculation of the beam-energy-dependence of heavy-ion
collision trajectories in the $T$-$n_{B}$ phase diagram, obtained from
simulations using two schematic EOSs, can be seen in Fig. 2.1 (note that the
trajectories are only evaluated at times when temperature extraction is fairly
well-defined). These wide ranges of system properties accessed in heavy-ion
collisions position them as a perfect tool to extract the nuclear matter EOS,
test predictions and extrapolations from regions of the QCD phase diagram
accessed by other approaches, and provide a valuable input to nuclear theory
and nuclear astrophysics calculations. For example, the density- and momentum-
dependence of the nuclear potential in both symmetric and asymmetric matter,
and thus of the corresponding EOS, can shed light on modeling effective
nuclear interactions in the medium [45, 111, 112, 113] or constrain approaches
using the density functional theory [114, 115, 116].

However, systems created in heavy-ion collisions are short-lived, and their
dynamics is out of equilibrium over significant fractions of the total
collision time. The evolution of a colliding system, which strongly depends on
both the energy

Figure 2.1: Phase diagram trajectories of the central region (cubic box of
volume $27\ \rm{fm}^{3}$) in Au+Au collisions at zero impact parameter,
obtained from UrQMD simulations with a soft or a hard (characterized by
$K_{0}=200$ or $K_{0}=380\ \rm{MeV}$, respectively) EOS [117, 118, 119]. The
trajectories only follow the evolution at times when temperature is fairly
well-defined, from the moment of the highest compression to densities around
$0.5n_{0}$.

and centrality of the collision, progresses through initial compression,
growth of the compression zone, development of flows, and overall
decompression with a gradual local equilibration throughout the process, see
Fig. 2.2. The inherent complexity of the evolution means that the
corresponding transport equations cannot be solved directly due to their high
non-linearity, and therefore any inferences from heavy-ion collision
experiments require comparisons to results of simulations. These are obtained
using transport models which are able to describe the non-equilibrium
evolution of nuclear matter over substantial ranges of density, as well as
naturally include baryon, strangeness, and charge diffusion, and describe
effects due to the interplay between the evolving collision zone and the
spectators, which are crucial for a correct description of, e.g., flow
observables. Beyond modeling the dynamics of the collisions, the dependence of
the evolution on single-particle interactions provides a connection allowing
one to use transport models for inferring equilibrium properties such as the
EOS [69, 73], transport coefficients [120], as well as the in-medium
properties and cross-sections of hadrons [121, 104, 122].

#### 2.1.1 Transport theory

At its core, transport theory aims to describe the time evolution a
dissipative system composed of a large number of particles (here, a system of
two heavy nuclei colliding at an energy per nucleon which is typically larger
than the Fermi energy) in terms of the one-body phase-space distribution
function in a semi-classical approximation. The theoretical foundations of
transport theory include the BBGKY hierarchy of coupled equations for reduced
density matrices [123] as well as the equations for nonequilibrium single-
particle Green’s functions, known as the Kadanoff–Baym equations [124, 125,
126]. In turn, the nonequilibrium Green’s function theory was built on the
equilibrium theory developed by Martin and Schwinger [127], retaining some of
its features, hence in this context the approach is often also called the
Martin–Schwinger (or Schwinger–Keldysh) formalism (see also Section 5.2).

Figure 2.2: Contour plots of the system-frame baryon density $n_{B}$ (top row)
and local excitation energy $E^{*}/A$ (bottom row) as a function of time
(columns from left to right), obtained from a transport simulation [128] of
the 124Sn+124Sn reaction at beam kinetic energy $E_{\rm{kin}}=800\ A\rm{MeV}$
($\sqrt{s_{{}_{\rm NN}}}=2.24\ \rm{GeV}$) and impact parameter $b$ = 5 fm. The
contour lines for the density use increments of $0.4n_{0}$, starting from
$0.1n_{0}$, while the contour lines for the local excitation energy correspond
to the values of $E^{*}/A=\\{5,20,40,80,120\\}\ \rm{MeV}$; the energy contour
plots have been suppressed for baryon densities $n_{B}<0.1n_{0}$. Figure
modified from [128].

To arrive at transport equations, one employs (among others) a Wigner
transformation and coarse-graining as well as a gradient expansion. The Wigner
transformation and coarse-graining lead to positive-definite phase-space
distributions [129] that can be efficiently sampled with Monte-Carlo
techniques, while the gradient expansion yields, for each particle species,
the force acting on a particle and the particle’s velocity as gradients of its
total energy with respect to the spatial position and momentum, respectively.
Knowledge of the kinematics of all particles, together with the elementary
collision rates, drives the evolution in the phase space. Finally, to arrive
at a set of Vlasov–Boltzmann-like equations, one employs the quasi-particle
approximation, neglecting details of the spectral functions and treating all
particles as on-shell (we note here that while there are some transport codes
with off-shell particle treatment, e.g., [130, 131, 132], this approach is
still an outstanding challenge in the transport theory, as will be discussed
further below). Alternative approaches to arriving at a transport theory for
heavy-ion collisions include using the relativistic Landau quasiparticle
theory [133] or, in approaches starting from a molecular picture, representing
the global wavefunction as a product (sometimes antisymmetrized) of single-
particle Gaussian wavepackets [134].

Figure 2.3: Comparison of rapidity distributions (top) and transverse flow of
nucleons (bottom) as functions of the scaled rapidity, obtained with different
transport codes (identified in the legend) within the TMEP initiative. The
results shown were obtained for 132Sn+124Sn collisions at beam kinetic energy
$E_{\rm{kin}}=270\ A\rm{MeV}$ ($\sqrt{s_{{}_{\rm NN}}}=2.01\ \rm{GeV}$) and
impact parameter $b=4\,\text{fm}$, using controlled input models for the EOS
and the cross sections as well as identically initialized nuclei. Figure from
[135].

The particle species considered in transport theory depend on the collision
energy and may range from nucleons, through pions and the delta resonances, to
higher resonances, kaons, and hyperons. Some transport formulations further
incorporate light clusters (e.g., deuterons, tritons, and 3He nuclei) as
independent degrees of freedom, with recent extensions also including alpha
particles [136] which appear abundantly in experiments and are of particular
importance for collisions at fixed-target beam energies on the order of
hundreds of MeV/nucleon. In some of these approaches, clusters are produced
through multi-particle reactions, as discussed further below. For the lowest
energy collisions, nonrelativistic formulations of the transport theory may be
employed, but the majority of the available codes are relativistic, with many
addressing collisions at energies from tens of MeV/nucleon to at least a few
GeV/nucleon (see [74, 137, 122] for reviews).

Transport approaches can be generally divided into those concentrating on a
single-particle characterization of the colliding system and those attempting
to describe many-particle correlations. Both types of approaches are highly
complex and nonlinear, and the relevant equations are solved by simulations.
The single-particle approaches typically solve a set of Boltzmann–Vlasov-type
equations [74, 138] (also known as the Boltzmann–Uehling–Uhlenbeck, or BUU
equations) in which the evolution of the system is governed by a mean-field
evolution of the phase space distribution (Vlasov equations) and a collision
term which drives the dissipation (the Boltzmann collision term). While, in
principle, the Boltzmann–Vlasov equation is deterministic, numerical solutions
contain numerically-induced fluctuations due to the fact that the evolution is
obtained using the method of test particles, in which the continuous
distribution function is represented by a large, but finite, number of test
particles sampling the phase space. To include fluctuations of a physical
origin, one can add a fluctuation term to the two-particle collision term,
thus arriving at the Boltzmann–Langevin formulation [122, 139].

In contrast, quantum molecular dynamics (QMD) approaches include classical
many-body correlations in the ansatz of the many-body wave function [74, 140],
which is postulated as a product of single-particle wave packets of a fixed
width, with the width regulating the amount of fluctuations and correlations
in QMD [141] (this width is usually fixed to reproduce realistic properties of
nuclei). In Anti-Symmetrized Molecular Dynamics (AMD) [134], the product wave
function is anti-symmetrized and the formulation includes Pauli correlations
in the propagation as well as in, to a certain extent, the collision term.

The fact that hadronic transport approaches are built on firm theoretical
foundations has been crucial for the continued development of simulation
frameworks. Reaching back to the roots of the nuclear transport theory has
made it possible to resolve ambiguities which would be otherwise hard to
tackle by purely phenomenological means, including descriptions of cluster
production [142], low relative-velocity correlations (Hanbury–Brown-Twiss
correlations) [143], and off-shell transport [130, 144, 145, 138, 146]. The
strong theoretical foundation of transport theory has also been effective in
ensuring covariance of the theory and preserving conservation laws in case of
interactions that stray beyond outcomes of field-theoretic models, in
particular interactions employing energy density functionals [133, 147, 148,
73] which are often needed for realistic descriptions of bulk properties of
nuclear matter. The many-body theoretical origin of transport equations
provides a connection between the mean fields, entering the drift terms, and
the collision integrals, both of which originate from consistent self-energies
in the many-body theory, and also provides expressions for spectral functions
describing the widths of all particles [126, 125, 149]. These theoretical
foundations both serve as a basis for currently used frameworks as well as
offer a means for future improvements and expansions.

An important effort to validate conclusions reached from comparing transport
model results to data has been recently intensified by the formation of the
Transport Model Evaluation Project (TMEP) [74]. Within this endeavor,
predictions from different models are compared in controlled settings (e.g.,
ensuring the same physical input such as the EOS, initial densities, and cross
sections), oftentimes with comparisons to known results that can be achieved
analytically or by other methods. Similar controlled comparisons of complex
simulations have been done in other fields of physics: from atomic traps,
through ultra-relativistic heavy-ion collisions, to core-collapse supernova
calculations [150, 151, 152, 153], and they are known to be very fruitful for
their respective fields. The TMEP analyses not only enable identifying models
that produce outlier predictions, but also determine details of implementation
or physical assumptions behind the diverging results. An example of such a
comparison of codes for simulations of heavy-ion collisions at lower energies,
with controlled input, can be seen in Fig. 2.3, showing results for rapidity
distributions (left) and the transverse flow (right) [135]. In general, the
codes agree with each other reasonably well, however, differences between the
codes are visible and, moreover, can be traced to specific model choices in
the simulations. For example, the generally lower values of the transverse
flow in the case of QMD codes are a result of an approximation used in the
evaluation of a non-linear term in the mean-fields, which becomes relevant
when density fluctuations become large, as often occurs in QMD [154, 155].
Beyond identifying this and similar problems, the Project has yielded
recommendations for optimal algorithms used in transport codes, e.g., for
ensuring obeying the Pauli principle in elementary two-body collisions [156]
or for integration of equations of motion with mean-fields [155]. Moreover,
the Project has identified a set of tests for transport codes that ensure
their credibility when addressing different heavy-ion collision observables.
Stringent tests of hadronic transport codes are especially important for
studies aimed at constraining the nuclear symmetry energy, which, compared to
other model parameters, has a comparatively weak effect on heavy-ion
observables and which therefore demands maximal precision from transport
simulations. Below, we will also discuss the role that such comparisons can
play in determining the uncertainty of transport model investigations.

#### 2.1.2 Selected constraints on the EOS obtained from heavy-ion collisions

A selection of important constraints on the EOS obtained from heavy-ion
collisions can be found in Fig. 2.4 for both symmetric matter (pressure as a
function of density, left panel) and asymmetric matter (symmetry energy as a
function of density, right panel). Additional constraints and measurements are
also discussed in Sections 3.1.2, 3.1.3, 3.2.3, and 4.1. All quoted beam
energies $E_{\rm{kin}}$ are the single-beam kinetic energies per nucleon, in
units of $A\rm{MeV}$ or $A\rm{GeV}$.

We note here that while many results are reported in terms of constraints on
the incompressibility $K_{0}$, in the context of heavy-ion collision studies
of the EOS, $K_{0}$ should be understood as a parameter which specifies the
behavior of the EOS in the range of densities probed by a given study. For
example, in the case of experiments probing mostly densities above $2n_{0}$,
constraints on $K_{0}$ are only indicative of the behavior of the EOS above
$2n_{0}$, and in particular do not constrain the behavior of the EOS around
$n_{0}$. This subtle, and often confusing, point is a consequence of simple
parametrizations of the EOS used in many transport codes, where the only
parameter controlling the behavior of the EOS both around $n_{0}$ and at
higher densities is $K_{0}$. Recently, flexible parametrizations of the EOS
have been developed (see, e.g., [148, 73]) and implemented (e.g., in hadronic
transport code SMASH [157, 158]) which allow one to vary the incompressibility
$K_{0}$ and the high-density behavior of the EOS independently, in turn
enabling description of non-trivial features at high density such as, e.g., a
phase transition.

The collective behavior of matter created in the collisions, especially the
directed and elliptic flow, has been shown to be a very sensitive probe of the
EOS [159, 160, 69, 72, 161]. In contrast to collisions at the Fermi energies,
where all nucleons within nuclei participate in the collisions, and unlike in
collisions at ultrarelativistic energies, where the evolution of the colliding
nuclei can be understood in terms of participant nucleons, at intermediate
energies the interplay between the expanding collision zone and the dynamics
of the spectators are key ingredients to understanding experimental results. A
seminal constraint on the symmetric nuclear matter EOS [69] in the density
range $(2$–$4.5)n_{0}$ was obtained by comparing measurements of collective
flow from heavy-ion collisions [162, 163, 164, 165] at beam energies
$E_{\rm{kin}}=0.15$–$10\ A\rm{GeV}$ (corresponding to nucleon-nucleon center-
of-mass energies $\sqrt{s_{{}_{\rm NN}}}=1.95$–$4.72\ \rm{GeV}$) with results
from hadronic transport simulations using EOSs with different values of the
incompressibility at saturation density $K_{0}$. The outcome of this study
suggests a symmetric-matter EOS to lie between those labeled with $K_{0}=210\
\text{MeV}$ and $K_{0}=300\ \text{MeV}$ (see the region with black horizontal
stripes in the left panel of Fig. 2.4). For densities in the range
$(1.0$–$2.5)\,n_{0}$, probed in collisions below $E_{\rm{kin}}\lesssim 1.5\
A\rm{GeV}$ ($\sqrt{s_{{}_{\rm NN}}}\lesssim 2.5\ \rm{GeV}$), the EOS may be
inferred from meson yields [166, 167, 168]. Indeed, subthreshold production of
strange mesons (specifically, $K^{+}$ and $K^{0}$), which interact weakly with
nuclear matter, depends on the highest densities sampled in the collision,
which in turn depend on the stiffness of the EOS [169]. In [70], ratios of
experimentally measured kaon yields in Au+Au and C+C collisions have been
reproduced in hadronic transport simulations with soft mean-field interactions
yielding $K_{0}=200\ \rm{MeV}$ and an EOS [71] consistent with the constraint
from [69] (see the region with red forward stripes in the left panel of Fig.
2.4). In [72], the elliptic flow data measured at $E_{\rm{kin}}=0.4$–$1.5\
A\rm{GeV}$ ($\sqrt{s_{{}_{\rm NN}}}=2.07$–$2.52\ \rm{GeV}$) by the FOPI
collaboration [170] were used together with simulations from Isospin Quantum
Molecular Dynamics (IQMD) [171, 172] to constrain the incompressibility at
$K_{0}=190\pm 30\rm{MeV}$, again indicating a rather soft EOS (see the region
with blue backward stripes in the left panel of Fig. 2.4). Recently, new
measurements by the STAR collaboration from the fixed target (FXT) program at
RHIC have become available, providing an opportunity to expand the set of
world data utilized to deduce the baryonic EOS. A Bayesian analysis study
[73], in which the speed of sound was independently varied in specified
intervals of baryon density (thus providing a more flexible EOS at higher
densities), suggests a tension between the E895 [173, 165, 174, 175] and STAR
[176, 77] data. Using only the STAR measurements, the study [73] further found
that EOSs which simultaneously describe the slope of the directed flow and the
elliptic flow, in the considered energy range of $E_{\rm{kin}}=2.9$–$8.9\
A\rm{GeV}$ ($\sqrt{s_{{}_{\rm NN}}}=3.0$–$4.5\ \rm{GeV}$), are relatively
stiff at lower densities and relatively soft at higher densities (see the
region with green vertical stripes in the left panel of Fig. 2.4). However,
the model used in that work did not include the momentum dependence of the
EOS, which has been to shown to result in a spuriously stiff EOS at
intermediate densities. As such, the study should be treated as a proof of
principle that a tight constraint on the EOS at high densities can be achieved
by using a combination of precise data, flexible forms of the EOS used in
simulations, state-of-the-art models, and advances in analysis techniques.

Figure 2.4: Left: Selected constraints on the symmetric EOS obtained from
comparisons of experimental data to hadronic transport simulations in [69]
(region with black horizontal stripes), [70, 71] (region with red forward
stripes), [72] (region with blue backward stripes), and [73] (region with
green vertical stripes); see text for more details. Also shown are results of
analytical calculations for the free Fermi gas (green dotted line) and in the
linear Walecka model (pink dashed line). Right: Selected constraints on the
symmetry energy obtained from comparisons of hadronic transport simulations to
experimental data in [36] (region with purple forward stripes), [177] (region
with green backward stripes), [178] (the solid orange region), and [65] (the
red circle, square, and triangle symbols). Also shown are symmetry energy
constraints obtained in [65] based on a novel interpretation of analyses of
dipole polarizability $\alpha_{D}$ [179] (green diamond), of nuclear masses in
DFTs [41, 180] (cyan dot symbol) and in Skyrme models [181] (cyan star
symbol), of Isobaric Analog States (IAS) energies [182] (magenta plus symbol),
and of PREX-II experiment [52] (blue inverted triangle symbol), as well as the
$68\%$ confidence region consistent with the best fit of experimental data
points presented in [65] (region with yellow vertical stripes).

The symmetry energy contribution to the EOS can be studied at low collision
energies $E_{\rm{kin}}\lesssim 1.0\ A\rm{GeV}$ ($\sqrt{s_{{}_{\rm
NN}}}\lesssim 2.32\ \rm{GeV}$), where in particular observables such as
charged pion yields [183] or neutron and proton flow [184, 185] have been
proposed as sensitive to the asymmetric contribution to the EOS. Some of the
constraints derived from such studies are shown in the right panel of Fig.
2.4, where, in addition to the usual EOS constraint bands, symbols with
uncertainty bars represent results from analyses in which the symmetry energy
has been determined for the most sensitive density of a given measurement. At
incident energies below $E_{\rm{kin}}=100\ A\rm{MeV}$ ($\sqrt{s_{{}_{\rm
NN}}}=1.93\ \rm{GeV}$), low densities are probed after the initial impact and
compression of the projectile and target [36, 186]. Since the symmetry
potentials for neutrons and protons have opposite signs, emission of a
particular nucleon type is enhanced or suppressed depending on the asymmetry.
A comparison of the experimental measurements of isospin diffusion and the
ratio of neutron and proton spectra in collisions of 112Sn+124Sn at
$E_{\rm{kin}}=50\ A\rm{MeV}$ ($\sqrt{s_{{}_{\rm NN}}}=1.90\ \rm{GeV}$) to
results from ImQMD [187] simulations produced a constraint on the symmetry
energy for densities $(0.3$–$1)n_{0}$ [36] (see the region with purple forward
stripes in the right panel of Fig. 2.4). Collisions at higher energies
($E_{\rm{kin}}>200\ A\rm{MeV}$, or $\sqrt{s_{{}_{\rm NN}}}>1.97\ \rm{GeV}$)
probe the EOS at $n>n_{0}$. In the FOPI-LAND experiment, constraints on the
symmetry energy were obtained from studies of the ratio of the elliptic flow
of neutrons and hydrogen nuclei in Au+Au collisions at $E_{\rm{kin}}=0.4\
A\rm{GeV}$ ($\sqrt{s_{{}_{\rm NN}}}=2.07\rm{GeV}$) [177], while the ASY-EOS
experiment used neutron to charged fragments ratios measured in Au+Au
collisions [178] (see the region with green backward stripes and the solid
orange region, respectively, in the right panel of Fig. 2.4). In [65], a
comprehensive analysis was performed with the goal of identifying the values
of the symmetry energy at densities to which given experiments are most
sensitive. Using the isospin diffusion in collision systems with different
proton to neutron ratios [61], neutron to proton energy spectra in Sn+Sn
systems [62], and spectral pion ratios measured by the S$\pi$RIT collaboration
in Sn+Sn collisions at $E_{\rm{kin}}=270\ A\rm{MeV}$ ($\sqrt{s_{{}_{\rm
NN}}}=2.01\ \rm{GeV}$) [63, 64], that work [65] put constraints on the values
of the symmetry energy at about $0.2n_{0}$, $0.4n_{0}$, and $1.5n_{0}$,
respectively (see the red circle, square, and triangle symbols in the right
panel of Fig. 2.4). Also shown in the right panel of Fig. 2.4 are symmetry
energy constraints obtained in [65] based on a novel interpretation of the
analyses of dipole polarizability $\alpha_{D}$ [179] (green diamond symbol),
of nuclear masses in DFTs [41, 180] (cyan dot symbol) and in Skyrme models
[181] (cyan star symbol), of the Isobaric Analog State (IAS) energies [182]
(magenta plus symbol), and of the PREX-II experiment result [52] (blue
inverted triangle symbol), as well as the $68\%$ confidence region consistent
with the best fit of experimental data points presented in [65] (region with
yellow vertical stripes), where identifying the specific densities at which
the measurements constrain the symmetry energy allowed for a comparatively
tight constraint at low densities.

#### 2.1.3 Challenges and opportunities

Selected results presented in Fig. 2.4 showcase significant achievements in
determining the EOS and, simultaneously, the need to develop improved
transport models to obtain tighter and more reliable constraints. Answering
this need will require support for a sustained collaborative effort within the
community to address remaining challenges in modeling collisions, in
particular in the intermediate energy range ($E_{\rm{kin}}\approx 0.05$–$25\
A\rm{GeV}$, or $\sqrt{s_{NN}}\approx 1.9$–$7.1\ \rm{GeV}$). In the following,
we will address selected areas where we see the need for such developments:
(1) comprehensive treatment of both mean-field potentials and the collision
term in transport codes, (2) use of microscopic information on mean fields and
in-medium cross sections, such as discussed in Section 2.2, in transport, (3)
better description of the initial state of heavy-ion collisions in hadronic
transport codes, (4) deeper understanding of fluctuations in transport
approaches, which affect many aspects of simulations, (5) inclusion of
correlations beyond the mean field into transport, which is crucial for a
realistic description of, e.g., light-cluster production, (6) treatment of
short-range-correlations in transport, which are tightly connected to multi-
particle collisions as well as off-shell transport, (7) sub-threshold particle
production, (8) connections between quantum many-body theory and semiclassical
transport theory, (9) investigations focused on extending the limits of
applicability of hadronic transport approaches, (10) studies of new
observables, e.g., azimuthally resolved spectra, to obtain tighter constraints
on the EOS, (11) the question of quantifying the uncertainty of results
obtained in transport simulations, and (12) the use of emulators and flexible
parametrizations for wide-ranging explorations of all possible EOSs.
Fortunately, advances in transport theory as well as the greater availability
of high-performance computing make many of these improvements possible.
Support for these developments will lead to a firm control and greater
understanding of multiple complex aspects of the collision dynamics, allowing
comparisons of transport model calculations and heavy-ion experiment
measurements to provide an important contribution to the determination of the
EOS of dense nuclear matter, which, in particular, cannot be determined by any
other method at intermediate densities (1–5)$n_{0}$.

Comprehensive treatment of mean-field potentials and the collision term

Over the last two decades, driven by specific experimental needs, the
refinement of hadronic transport codes has diverged into two complementary
branches: Codes which were applied to describing experiments at very low to
low energies ($E_{\rm{kin}}\lesssim 1.5\ A\text{GeV}$, or $\sqrt{s_{{}_{\rm
NN}}}\lesssim 2.5\ \text{GeV}$), such as IQMD [188, 172], AMD [134, 189] and
pBUU [190, 147], have become progressively better at describing the momentum-
and isospin-dependence of the interaction, while codes which were primarily
designed for simulations of relativistic and ultra-relativistic heavy-ion
collisions ($E_{\rm{kin}}\gtrsim 25\ A\text{GeV}$, or $\sqrt{s_{{}_{\rm
NN}}}\gtrsim 7\ \text{GeV}$), such as SMASH [157] or UrQMD [117, 118], were
developed to offer a fully relativistic evolution as well as scattering and
decay modes taking into account a large number of established particle and
resonance species. As heavy-ion collisions are entering an era of precision
data on symmetric nuclear matter at high densities (e.g., in experiments at
HADES, BES FXT, and future CBM, probing densities up to several times the
saturation density) and on asymmetric nuclear matter at normal and supranormal
densities (e.g., at FRIB and future FRIB400, probing densities up to twice the
saturation density), where features of both diverging branches of hadronic
transport codes are important, a vigorous development of transport models is
needed to incorporate all relevant physics.

For example, numerous studies show the importance of including the momentum-
dependence of the interactions, which is observed in elastic scattering of
hadrons off nuclei. Moreover, momentum-dependence naturally occurs in
microscopic effective interactions [126, 191] where it contributes to the
calculated mean fields, whether near or away from saturation density.
Incorporating single-particle energies with momentum dependence different than
that in free space, which is often quantified with effective masses, is
crucial in hadronic transport both for studies of symmetric nuclear matter
[192, 193, 69, 161] as well as studies of the symmetry energy and its relation
to effects such as the neutron-proton effective mass splitting [194, 195, 196]
(see also Section 6.5 for more discussion on effective masses and the nuclear
symmetry energy).

As another example, inclusion of high-mass resonant degrees of freedom, which
can be produced in substantial amounts at collision energies
$E_{\rm{kin}}\gtrsim 1\ A\text{GeV}$ ($\sqrt{s_{{}_{\rm NN}}}\gtrsim 2.32\
\text{GeV}$), can significantly affect the evolution of the system. In
particular, since some part of the collision energy can be used for resonance
production (so that a non-negligible fraction of the nuclear matter becomes
“resonance matter”), it may affect both the degree of initial compression and
the spectrum of mesons emitted and absorbed throughout the evolution, which
has consequences for the collective flow [197] and thus directly influences
extraction of the EOS. Therefore, inclusion of all relevant resonant species
and ensuring correct description of meson production and absorption is crucial
for a reliable inference of the EOS.

While some of the theoretical and implementation solutions needed to improve
the models in the intermediate energy region have already been established,
others will require devising new approaches. When possible, best practices
need to be carried over across the domains, as has been exemplified in, e.g.,
the development of the SMASH code, which combines multiple well-tested
implementation solutions from UrQMD, GiBUU [146], and pBUU.

Microscopic input to transport

One of the most prominent opportunities for improvement in transport models
concerns implementations of the EOS informed by state-of-the-art many-body
studies. Such efforts are especially timely given that sophisticated
microscopic calculations of the properties of nuclear matter are currently
becoming available for large ranges of baryon density, temperature, and
isospin fraction (see Section 2.2 for more details). To incorporate the
effects of the resulting EOSs in hadronic transport calculations, the
corresponding Lorentz-covariant single-particle potentials as well as the in-
medium interactions (both as functions of density, asymmetry, and momentum)
are needed. A particular challenge is to determine the connection between the
EOS inferred from a transport calculation and the zero-temperature EOS
obtained from microscopic calculations [198], or even the finite-temperature
EOSs that are becoming increasingly available [199, 200]. In a heavy-ion
collision, the medium progresses through a set of non-equilibrium states that
relax toward a local equilibrium, however, the nature of the local equilibrium
also evolves during the collision due to the system expansion, so that even if
the system approaches a local equilibrium at any given moment of the
evolution, that agreement is only temporary. Errors incurred due to
differences between non-equilibrium and equilibrium states of high-density
matter contribute to the systematic error in inferring the EOS when comparing
transport to experimental data (see Fig. 2.4 and [69]). Here, the availability
of state-of-the-art microscopic calculations at finite temperature could
reduce systematic errors in connecting the finite- and zero-temperature EOSs.
Moreover, the use of microscopic input would provide a consistency between the
effective in-medium cross sections in the collision term and the mean fields
used in the propagation of the phase space distribution. It could also help
address the question of the extent to which nonlocalities in the microscopic
theory should be reflected in the propagation and the collision term [201,
202] (where, in particular, departures from standard approaches modify the
entropy to take a form different than that obtained in the Landau
quasiparticle theory [133, 203]). Microscopic calculations could also inform
the treatment of strange degrees of freedom in hadronic transport through a
better understanding of nucleon-hyperon and hyperon-hyperon interactions,
which can be pursued within, among others, lattice QCD approaches (see, e.g.,
[204]).

To accelerate progress at the interface of the transport description of heavy-
ion collisions and microscopic nuclear matter theory, direct collaboration of
practitioners in the two research areas is required to assess how the needs of
transport simulations can be answered by what can be currently calculated in
microscopic theories. Conversely, the use of microscopic interactions in
transport could validate the many-body theory results in regions of density
and temperature which are only accessible by heavy-ion collisions [205].

Initial state

Numerous studies point toward the dependence of outcomes of heavy-ion
collision experiments on details of the initial conditions. In
ultrarelativistic heavy-ion collisions, understanding these effects have led
to the discovery of higher order flow harmonics [95, 96] and flow fluctuations
[206]. (Interestingly, the importance of the initial state for experimental
outcomes also positions heavy-ion collisions at high energies as an unusual,
but complementary probe of nuclear structure, see, e.g., a white paper on
Imaging the initial condition of heavy-ion collisions and nuclear structure
across the nuclide chart [207].) Given the high sensitivity of flow
observables to both the EOS and the initial state of collisions, the impact of
the initial conditions on outcomes of heavy-ion collisions needs to be
thoroughly understood in order to narrow the constraints on the EOS of both
symmetric and asymmetric matter. Aspects of initial conditions that need to be
considered include event-by-event fluctuations of the initial state [95, 96,
206], relative distributions of neutrons and protons and shell effects [208],
and correlations tied to deformation [209] or short-range correlations [210].
Some of these elements will be further discussed below in the context of the
dynamics of heavy-ion collisions.

Fluctuations

Fluctuations of the phase space distribution are an important ingredient of
transport simulations. In particular, fluctuations of the one-body density are
important for including the consequences of the dissipation-fluctuation
theorem in the reaction dynamics as well as for describing effects due to the
largely unknown, neglected many-body correlations, thus going beyond the mean-
field description. The question of how to include them properly and of their
consistency with the nucleon-nucleon correlations explicitly implemented in
transport theories, however, has not been completely clarified. As discussed
above, fluctuations are included in a different manner in the two families of
transport approaches. While in the BUU transport fluctuations can be
introduced by the Langevin extension of the Boltzmann–Vlasov equation, which
adds a fluctuation term to the collision term (and which is still rarely
implemented), in the molecular dynamics approach fluctuations are introduced
in a classical way by using finite-size particles, the width of which
regulates the amount of fluctuations. Fluctuations then affect the outcome of
simulations in many ways, including by regulating the formation of
intermediate-mass fragments (IMFs) which appear through the growth of
fluctuations in regions of spinodal instability. It was also shown in box
calculations that fluctuations have a strong influence on the efficiency of
Pauli-blocking [156] and even on the calculation of the force in the Vlasov
term for QMD codes in which non-linear parametrizations of the fields are used
[155].

Correlations

Correlations in transport simulations strive to address intermediate-range
correlations beyond the mean-field picture, that is beyond fluctuations of the
one-body density. Physically, such correlations also contribute to one-body
density fluctuations, but at the same time they have other additional impacts,
including, e.g., influencing the production of light clusters (LCs), that is
light nuclei up to the alpha particle which are copiously produced in heavy-
ion collisions. However, the mean-field models used in transport calculations
are usually not detailed enough to realistically describe very light nuclei
with their particular spin-isospin structure reflecting strong quantum
effects. An additional complication results from the fact that in a collision,
clusters often appear in the nuclear medium where their properties are
drastically changed (e.g., the binding energy of clusters is reduced with
increasing density until the Mott point, at which they dissolve). Currently,
most codes describe the production of clusters by using a cluster-finding
algorithm, based on particle proximity in coordinate and/or momentum space
(coalescence) toward the end of the evolution, which in more advanced versions
also takes into account criteria related to the binding energy of the produced
clusters [211]. However, these late-stage algorithms do not take into account
the dynamic role played by both correlations and LCs in the evolution of the
collision. One of the known approaches to this problem has been to consider
LCs as separate degrees of freedom, with their own distribution functions and
corresponding transport equations, where the collision terms can lead to
creation or destruction of clusters (pBUU, SMASH) and which in particular can
also take into account the in-medium modifications of clusters. However, this
approach becomes increasingly complex as heavier clusters are characterized by
more and more production channels, and consequently it is significantly
challenging to include, e.g., alpha particles. Another approach is to modify
the phase space of the correlated nucleons according to the Wigner function of
the cluster, but then to propagate them after the collision again as nucleons,
which still requires using a cluster-finding step at the end (this is done in,
e.g., the AMD code [129], within which it has also been demonstrated that the
clustering effects may influence pion production [189]). In both approaches,
the production and destruction of clusters necessarily requires multi-particle
collisions to ensure energy-momentum conservation. Finally, at lower incident
energies the LC production can also be described in terms of the catalyzing
effect of spectator nucleons in few-particle collisions [190, 136]. To explain
LC production in high-energy collisions, where LCs are produced in numbers
that cannot be obtained through nucleon catalysis due to the relatively few
nucleons present in the final stages of these collisions, a similar mechanism
of catalysis by pions [212, 142, 213] can be invoked.

Short-range correlations

A particular aspect of describing correlations in transport simulations is the
treatment of short-range-correlations (SRCs), which have been measured in
nucleon knock-out experiments [214, 215, 216, 217]. Along with the
experiments, microscopic many-body calculations show that SRCs introduce a
high-momentum tail (HMT) into the nucleon momentum distribution and, moreover,
reduce the kinetic symmetry energy relative to the Fermi gas kinetic energy,
which is a consequence of the fact that SRCs are more pronounced in symmetric
relative to asymmetric matter [218, 219, 220, 221, 222, 223, 224, 225] (see
also Section 6.3). Phenomenological methods have been used to include SRCs in
transport models, e.g., by initializing nuclei with a HMT, but such a
procedure does not take into account the dynamic role of SRCs in the initial
state, which in the case of the on-shell semiclassical equations of motion
results in obtaining non-stationary, excited states of nuclei. In on-shell
transport approaches, three- and many-body collisions, incorporated into
transport codes within varying approximations, have been suggested as a way of
treating SRCs. In particular, in an investigation [226] of three-body
collisions for pion production processes (e.g., $NNN\rightarrow NN\Delta$), it
was found that SRCs between two of the incident nucleons give a noticeable
contribution to pion yields. Another approach [227], based on a mean-free-path
approximation to the collision integral, observed large effects also on bulk
observables. The incorporation of $n$-body collisions in transport equations
within a schematic cluster approximation was also studied [228], however, the
effects were found to be rather small. So far, none of these methods have been
widely exploited in the description of heavy-ion reactions. Since HMTs are
tied to the tails of the nucleon spectral functions (away from the
quasiparticle peaks), a consistent description of SRCs should involve an off-
shell transport formulation.

Threshold effects

An important influence of mean-field potentials in heavy-ion transport appears
in the form of threshold shifts and the related subthreshold production of
particles. Thresholds of particle production are modified in a medium since
the mean-field potentials have to be taken into account in the energy-momentum
balance of a two-body collision. Specifically, when the mean-field potentials
are momentum-dependent and/or as a consequence of other model assumptions for
the mean-field potentials of the produced particles, the thresholds are
shifted away from their free-space values. This may strongly change the
production rates of particles. Moreover, the threshold shifts make it
necessary to involve other nucleons, besides the two collision partners in the
process, to ensure the energy-momentum conservation. Various schemes to
achieve this locally or globally have been in use [229, 205]. Indeed,
explaining recent heavy-ion collision subthreshold pion yields, measured by
the S$\pi$RIT Collaboration [64], required invoking many-body elementary
effects in the form of mean-field effects on thresholds in two-particle
collisions [63, 168]. However, because the physics invoked in describing the
threshold effects is similar to that invoked for other multi-particle effects,
alternative multi-particle options remain to be investigated, including
producing pion degrees of freedom in multi-particle collisions or in the
aftermath of an off-shell propagation between binary collisions. (We note here
that there is a physics overlap between these mechanisms and the impact of
SRCs on pion production [210, 226, 130].) Notably, theoretical explorations
find sequences of on-shell binary processes to dominate the production at
higher beam energies [226, 131, 145], and no comparable difficulties have been
encountered in describing the data [230, 231] by transport models without
multi-particle effects. The contrasting struggles of transport models which do
not include threshold or other multi-particle effects of this type [64],
together with expected further theoretical explorations and future
measurements of the subthreshold production in heavy-ion collisions, offer
exciting possibilities for gaining understanding of the more exotic in-medium
processes.

Interface between quantum many-body theory and semiclassical transport theory

In deriving the transport theory in the nonequilibrium Green’s functions
formalism, one considers the evolution of the one-body Green functions up to
1-body level in terms of self-energies. In the semiclassical limit, these
self-energies yield both the potentials and the in-medium cross sections in
the collision terms of the transport equations [126, 125, 149, 146]. In the
usual applications of the transport theory, these driving terms are modeled
with phenomenological density functionals and effective in-medium cross
sections, respectively. Within a more microscopic approach, the self-energies
would be calculated in the $T$-matrix approximation similar to that employed
in the Brueckner theory. The Hermitian part of the resulting self-energy, in a
given intrinsically consistent approximation, would yield the potential for
the drift term in transport. The anti-Hermitian part would yield widths of the
spectral functions as well as contribute to the feeding and depletion in the
transport collision integrals. By retaining additional terms in the expansion
leading to transport, one may arrive at equations where the phase-space
evolution is coupled to the evolution of the spectral functions. With some
additional approximations, the spectral functions for particles (including
those which are stable in free space such as, e.g., nucleons) have been
implemented in GiBUU [146], PHSD [232, 132], and PHQMD [233] codes. A study
[234], performed within the approach developed in GiBUU, showed that the
momentum distribution naturally develops a HMT in off-shell transport.

Further development of these aspects will become more important as data and
theory become more precise. Their impact on inferences of the symmetry energy
from heavy-ion collision data based on, e.g., charged pion subthreshold
production yields, can be particularly consequential, but has yet to be
systematically investigated. Fully quantum transport approaches with SRCs (or
equivalent content), without any semi-classical expansions as are present in
current off-shell transport approaches, remain a long-term goal, and progress
in this area has not ventured yet beyond schematic models [235, 236]. However,
increasing computational power combined with emulation techniques may make
such efforts more realistic and enable, e.g., a seamless integration of the
treatment of shell effects in the initial state and collision dynamics.

Limits of hadronic transport

In heavy-ion collisions at higher beam energies, increasing numbers of nucleon
resonances and heavier mesons are produced, which at a certain point start to
significantly affect the evolution of the system (this includes the emergence
of “resonance matter”, mentioned above). In some transport codes, e.g., UrQMD
or SMASH, a large number of such states is included. Eventually, however, as
one approaches energies at which quark-gluon plasma (QGP) is produced, sub-
nuclear or partonic degrees of freedom may be needed for a correct description
of systems created in heavy-ion collisions. Notably, some transport codes
include both hadronic and partonic degrees of freedom simultaneously (e.g.,
PHSD [232]).

The separation of scales underlying transport theory, which can be encompassed
in the condition that the duration of a collision is short compared to the
mean free flight time, is increasingly undermined when the systems probe large
densities. The latter also lead to fast equilibration and the associated
hydrodynamic behavior of the system. In principle, transport theory can be
used to model hydrodynamics, even in cases when the degrees of freedom
employed in transport are not the ones truly underlying the system. At some
point, however, modeling heavy-ion collisions in terms of hydrodynamics
becomes more straightforward, especially given the fact that the influence of
the degrees of freedom believed to be appropriate can be easily included both
in the EOS and the transport coefficients, which can be pursued within
microscopic approaches [237, 238]. At the same time, hydrodynamic modeling of
heavy-ion collisions may face certain challenges below $E_{\rm{kin}}\approx
25\ A\rm{GeV}$ ($\sqrt{s_{{}_{\rm NN}}}\approx 7\ \rm{GeV}$), where not only
is the assumption of near-equilibrium not very well satisfied, but where one
also needs to take into account, among others, the hydrodynamic evolution of
all conserved charges (e.g., baryon number $B$, strangeness $S$, and electric
charge $Q$), possible evolution through unstable regions of the phase diagram,
or the influence of the spectators. Extensions of hydrodynamics that may make
it possible to consistently apply it in this regime are discussed in Section
5.2.

New observables

Upcoming precision data will further bring unprecedented observables that
could be previously considered only in theory, such as triple-differential
spectra tied to a fixed orientation of the reaction plane [239, 240, 241, 79]
not only for protons (see Fig. 3.4) and most abundant mesons, but also for
deuterons, tritons, light nuclei, and hypernuclei. The potential of such
spectra for the determination of the EOS is still to be fully explored, but a
preliminary investigation [240] indicates a rich structure with spectra which
exhibit a maximum away from the beam direction, characterized by slopes
dependent on azimuthal angle and slope discontinuities. Models that might have
agreed with each other in describing low-order Fourier coefficients of flow
will likely find describing such detailed observables difficult. Challenges
remain even at the level of the low-order coefficients, as many models now
reproduce proton flow, but not Lambda or pion flow (see, e.g., Fig. 3.3).
Understanding the relations between observables for various particle species
will lead to constraints on the physics driving the evolution of heavy-ion
collisions in simulations and, through that, to understanding cluster
formation, hyperon yields, in-medium interactions with of strange hadrons, and
more (see also the white paper on QCD Phase Structure and Interactions at High
Baryon Density: Continuation of BES Physics Program with CBM at FAIR [99]).

Quantifying uncertainties of transport predictions

In the era of multi-messenger physics, where information on the EOS is derived
from different areas of physics such as nuclear structure, nuclear reactions,
and astrophysics, the ability to assess the uncertainty of a particular result
is of crucial importance. This problem is especially relevant for evaluations
of constraints on the EOS from transport simulations of heavy-ion reactions,
since it has been found that using different transport models to describe the
same data can lead to very different conclusions. As found in the TMEP
comparisons (see [74] for a review), even with controlled input the results
from different models may vary considerably due to different implementation
strategies which in themselves are not dictated by the underlying physics. In
such a situation, calculating the mean and variance of different model
predictions is not a reliable way of determining the uncertainties. An
approach currently considered for ensuring a robust quality control in
combining inferences from different models is to weigh the models with a
Bayesian weight which could be based, e.g., on the performance of a given
model in benchmark tests and/or its ability to reproduce all key observables
of a given reaction (for example, flow observables, particle multiplicities,
and spectra). Bayesian analysis can be also used for model selection through a
comparison of results from a list of available models with data, during which
one assigns to each model a probability of being correct based on the quality
of the fit. However, this approach implicitly assumes that among the
considered models there is at least one “true” model (also known as the
$\mathcal{M}$-closed assumption), which is often not fulfilled. Efforts have
been taken to analyze data with an $\mathcal{M}$-open assumption, where the
existence of a perfect model is not assumed. For nuclear physics efforts, this
is being attempted within the Bayesian Analysis of Nuclear Dynamics (BAND)
group [242] by using Bayesian model mixing, where information from different
models is combined for inference.

Emulators and flexible EOS parametrizations

Robust explorations of the possible physics underlying various observables
often necessitate repeating the calculations many times for different
combinations of physics parameters. When high event statistics is needed, the
computational task can easily overwhelm the available computational resources.
An additional computational strain often arises from assessing Bayesian
probability distributions for any conclusions. Increasingly, emulators are
going to be used for this task, with some steps having been already made [64,
243, 73]. Notably, similar issues emerge in the area of applications of
hadronic transport [244] (see also Section 5.1).

For explorations focused on the EOS, it may be of advantage to fit various
possible EOSs with flexible relativistic density functionals as suggested in
[148, 158]. This approach, given the freedom in varying both the functional
form of the EOS as well as the EOS parameters, is particularly amenable to
Bayesian analyses (see, e.g., [73] for a Bayesian analysis with a
parametrization of the EOS in terms of the functional dependence of the speed
of sound on density).

The above list of issues facing the application of transport theory to heavy-
ion collisions highlights the fact that this approach to putting tighter
constraints on the EOS rests on overcoming certain challenges. In simple
terms, one attempts here to use a very dynamic and complex non-equilibrium
process to obtain information describing a relatively simple and well-defined
system, namely the equilibrated EOS of nuclear matter for different densities,
temperatures, and isospin asymmetries. To achieve this in a reliable way,
multiple complex issues of many-body physics have to be well controlled. On
the other hand, several of the needed improvements are relatively well-
understood, and tackling some of the unresolved problems poses an exciting
intellectual challenge. As a reward for undertaking this effort, one gains the
opportunity to obtain information on the EOS in a region which cannot be
accessed through any other means: For densities below saturation, there is
strongly constraining information from nuclear structure, with significant
contributions coming also from low-energy heavy-ion collisions. Astrophysical
observations on neutron stars and neutron star mergers are mainly sensitive to
densities above about $3n_{0}$. The gap between these domains can only be
filled with intermediate energy heavy-ion collisions, and transport studies
are the essential tool to extract the information on the EOS from experimental
data.

### 2.2 Microscopic calculations of the EOS

Over the past decade, many-body nuclear theory has made significant progress
in deriving microscopic constraints on the nuclear EOS at low densities from
chiral effective field theory ($\chi$EFT) [47, 48, 49, 50]. The progress has
been driven by improved two-nucleon (NN) and three-nucleon (3N) interactions,
rigorous uncertainty quantification, and algorithmic and computational
advances in the frameworks used to solve the many-body Schrödinger equation
with these interactions (see also the recent white paper on Dense matter
theory for heavy-ion collisions and neutron stars [245]).

#### 2.2.1 Status

Chiral EFT [83, 84, 246, 247, 248] provides a systematic way to construct
nuclear interactions consistent with the low-energy symmetries of QCD, using
nucleons (N’s), pions ($\pi$’s), and (in the case of delta-full $\chi$EFT),
$\Delta$-resonances ($\Delta$’s) as the relevant effective degrees of
freedoms. Nuclear interactions in $\chi$EFT are expanded in powers of momenta
or the pion mass over a hard scale at which $\chi$EFT breaks down; this
breakdown scale is expected to be of the order of the $\rho$-meson mass,
$\Lambda_{b}\approx 600$ MeV. At each order in the EFT expansion, only a
finite number of diagrams enter the description of the interaction according
to a chosen power counting scheme, of which the Weinberg power counting has
been predominant. (Several alternative power counting schemes have been
developed, see, e.g., [249, 250, 251, 252]; however, these approaches are not
well-developed at higher expansion orders and in some cases fail to reproduce
properties of nuclear matter at leading- or next-to-leading order.) At the
leading-order (LO) in Weinberg’s power counting one includes contribution from
the one-$\pi$ exchange between two nucleons as well as momentum-independent
contact interactions, which allow one to describe key features of the nuclear
interaction already at the lowest order. At next-to-leading-order (NLO),
two-$\pi$ exchanges are included as well as momentum-dependent contact
interactions, and similarly, more involved terms appear at higher orders. The
various low-energy coupling constants are determined from fits to experimental
data, e.g., the $\pi$-N couplings are fit to $\pi$-N scattering, while those
describing NN short-range interactions are fit to NN scattering data. The
advantage of $\chi$EFT over phenomenological approaches is that multi-nucleon
interactions, such as the important 3N interactions, naturally emerge in the
EFT expansion and, moreover, are consistent with the NN sector. The
contribution of chiral EFT 3N forces to the neutron-matter EOS naturally
includes both the interaction as well as method uncertainties, typically
adding a repulsion of about 2–4 MeV at nuclear saturation density [253, 254,
255, 256]. Forces involving increasingly more nucleons are correspondingly
more suppressed, e.g., the leading contribution to 3N forces (four-nucleon
(4N) forces) appears at N2LO (at N3LO) in Weinberg’s power counting.
Furthermore, there are only two new low-energy couplings appearing in the
three- and four-body forces to N3LO, which govern the strengths of the
intermediate- and short-range contribution to the leading 3N forces,
respectively. Consequently, $\chi$EFT 3N and 4N interactions at N3LO are
completely determined by constraints on the coupling constants obtained from
NN and $\pi$-N scattering, usually resulting in tight constraints on very
neutron-rich matter from $\chi$EFT.

Another key feature of $\chi$EFT is that order-by-order calculations in the
$\chi$EFT expansion have enabled estimation of theoretical uncertainties due
to truncating the chiral expansion at a finite order [257, 43, 258, 50].
Quantifying and propagating these EFT truncation errors enables meaningful
comparisons between competing nuclear theory predictions, see Fig. 2.5, and/or
constraints from nuclear experiments and neutron-star observations in the
multi-messenger astronomy era [259]. Such comparisons are facilitated by
Bayesian methods in a statistically rigorous way [89, 50, 259] to take full
advantage of the great variety of empirical EOS constraints we anticipate in
the next decade.

Figure 2.5: Comparison of the energy per particle $E/N$ (left) and the
pressure $P$ (right) as functions of density for pure neutron matter in
different many-body calculations using interactions from $\chi$EFT. The left
panel also shows low-density QMC results of Ref. [260] and the conjectured
unitary-gas lower bound on the energy per particle of pure neutron matter from
Ref. [46]. Figure from Ref. [261].

Chiral EFT also provides nuclear Hamiltonians governing the interactions in
nuclear systems. However, to calculate properties of a many-body system,
computational methods able to solve the Schroedinger equation for this system
are necessary. Among various frameworks used to solve the nuclear many-body
problem in dense matter, quantum Monte Carlo (QMC) methods and many-body
perturbation theory (MBPT) have been the main tools employed to study the
physics of neutron-star matter in recent years. Both methods have recently
made tremendous advances in predicting properties of nuclei and calculating
the nuclear matter EOS [262, 48, 263, 264, 265, 266, 50].

QMC frameworks, such as the auxiliary field Diffusion Monte Carlo (AFDMC)
method, are based on imaginary-time propagation of a many-body wave function
and enable us to extract ground-state properties of a nuclear many-body system
with high statistical precision [262, 48]. Their nonperturbative nature also
allows for the treatment of nuclear interactions at high momentum cutoffs,
providing important insights into nuclear interactions at relatively short
distances that may help to improve the modeling of $\chi$EFT interactions. QMC
calculations of binding energies, radii, and electroweak transitions of nuclei
up to $A=16$ [267, 268, 269, 270, 271, 272, 273] using $\chi$EFT NN and 3N
interactions are in very good agreement with experimental data [274, 275, 276,
277]. QMC methods were also used to calculate the EOSs of matter up to about
twice the nuclear saturation density $n\approx 2n_{0}$ [278, 92, 279, 280,
281]. The calculated EOSs include estimates of systematic truncation
uncertainties, and are commonly used to constrain properties of neutron stars
[92, 282, 283].

The past decade has also seen a renaissance for many-body perturbation theory
(MBPT) calculations in nuclear physics [50, 266]. Key to this development has
been the discovery that nuclear potentials with momentum-space cutoffs in the
range $400\,{\rm MeV}\lesssim\Lambda\lesssim 500$ MeV (not to be confused with
the breakdown scale of $\chi$EFT, $\Lambda_{b}$) are sufficiently soft to
justify the use of perturbation theory methods [253] (see [284] for a Weinberg
eigenvalue analysis). Such low-momentum potentials can be obtained from
renormalization group methods [285] or by directly constructing chiral
effective field theory potentials at a coarse resolution scale. Furthermore,
recent advances in automatic diagram generation [286] combined with automatic
code generation [287] and high-performance computing have led to a fully
automated approach to MBPT calculations in nuclear physics [50], in which
chiral two- and multi-nucleon forces can be included to high orders in the
chiral and MBPT expansions. MBPT has been demonstrated to be a computationally
efficient and versatile tool for studying the nuclear EOS as a function of
baryon number density $n_{B}$, isospin asymmetry
$\delta=(n_{n}-n_{p})/(n_{n}+n_{p})$, and temperature $T$ [199, 288, 289, 200]
with implications for neutron star structure [50] and astrophysical
simulations [290]; here, $n_{n}$ and $n_{p}$ correspond to the neutron and
proton densities, respectively. In particular, MBPT allows us to compute the
EOS of neutron-star (i.e., $\beta$-equilibrated) matter explicitly, which can
help improve isospin asymmetry expansions of the low-density nuclear EOS such
as the standard quadratic expansion [291, 292, 293, 288, 294, 295]. MBPT also
allows us to study nuclear properties other than the nuclear EOS, including
the linear response and transport coefficients that could be used to inform
more accurate numerical simulations of supernovae and neutron-star mergers.
Furthermore, MBPT for (infinite) nuclear matter has been used to construct a
microscopic global optical potential with quantified uncertainties based on
$\chi$EFT NN and 3N interactions [296, 297]. Altogether, MBPT calculations of
nuclear matter properties can provide important constraints that enable
microscopic interpretations of future nuclear reaction experiments [298]
(e.g., at the Facility for Rare Isotope Beams) and neutron star observations.

To date, theoretical predictions for the nuclear EOS, optical potentials, and
in-medium NN scattering cross sections have been computed at finite
temperature at various levels of approximation starting from fundamental two-
and multi-nucleon forces. These quantities are inputs to transport model
simulations [299, 171] of heavy-ion collisions used to extract constraints on
the properties of hot and dense nuclear matter (see Section 2.1 for more
details). In transport simulations, the EOS, single-particle potentials, and
in-medium NN cross sections are usually obtained from effective
phenomenological interactions [300, 301] that are fitted to the properties of
finite nuclei and cold nuclear matter, and then extrapolated into the finite-
temperature regime. Recently, some effort has been devoted to benchmarking
[198] the temperature dependence of these effective interactions against
predictions from $\chi$EFT or directly using EFT constraints in fitting
effective interactions [302, 303, 304]. To enable such comparisons, the free
energy of homogeneous nuclear matter as a function of temperature, baryon
number density, and isospin asymmetry has been calculated using $\chi$EFT
interactions up to second order in many-body perturbation theory [199] and
within the Self-Consistent Green’s Function (SCGF) approach [305], which
resums particle-particle and hole-hole ladder diagrams to all orders. The
resulting EOS has been shown to be consistent with the critical endpoint of
the symmetric nuclear matter liquid-gas phase transition [199, 305] as well as
the low-density/high-temperature pure neutron matter EOS from the virial
expansion [293]. Furthermore, single-particle potentials have been computed at
finite temperature at the Hartree-Fock level [306], from G-matrix effective
interactions [307], and in SCGF theory [308, 290]. Of particular importance is
the associated nucleon effective mass, which is obtained from a momentum
derivative of the single-particle energy. The nucleon effective mass is
directly related to the density of states and hence governs entropy generation
at finite temperature, with consequences for the dynamical evolution of core-
collapse supernovae and neutron star mergers. Finally, in-medium NN scattering
cross sections have been computed at finite density and zero [309] as well as
at finite [307] temperature using high-precision nuclear forces. Overall,
improved effective field theory methods can provide much-needed microscopic
inputs for transport simulations of heavy-ion collisions and astrophysical
simulations, leading to a consistent description of microscopic processes
affecting their dynamical evolution.

#### 2.2.2 Challenges and opportunities

To fully capitalize on experimental and observational data and extract key
information on fundamental questions in nuclear physics, continued progress in
nuclear theory is crucial. The combination of $\chi$EFT with modern
computational approaches like machine learning, artificial intelligence,
emulators, and Bayesian inference have provided EOS results for a wide range
of densities, and at various proton-to-neutron asymmetries and temperatures,
with quantified uncertainties [200, 258] (see also [310] for a broader review
of machine learning in nuclear physics). Future progress in the development of
fundamental interactions, combined with these tools, will increase the
precision of the results and enable us to answer open problems in chiral EFT.
Among these, the most pressing is at which densities and how $\chi$EFT breaks
down [92, 258]. In particular, for studies of neutron-star mergers it is of
great importance to describe dense matter at finite temperatures [293, 290,
289], however, these might influence the breakdown of the theory in dense
matter. In the next decade, it will be crucial to reliably determine how far
one can push the $\chi$EFT approach in nucleonic matter.

While microscopic calculations have been very successful in calculating
properties of nuclei and homogeneous matter at densities up to 1–2 times the
nuclear saturation density, we need improved microscopic descriptions of
neutron-rich dense matter beyond that regime, at a few times nuclear
saturation density and finite temperatures, with quantified uncertainties.
This can be achieved by employing models derived within relativistic mean-
field or density functional theory that are firmly rooted in microscopic
theory at lower densities. Such models will be very important to connect
theoretical calculations within the framework of $\chi$EFT to heavy-ion
collision experiments at accelerator facilities around the world. Heavy-ion
collision experiments at intermediate beam energies bridge the low- and high-
density regimes of the EOS and provide complementary information to that
obtained from nuclear structure or neutron-star studies [67] (see Section
2.1). Robust inferences from the experimental data will require more accurate
predictions from transport theory, which strongly depend on, among others,
mean-field or density functional models. It will be imperative to test and
constrain such models for the EOS with more rigorous microscopic calculations.
Beyond their use in hadronic transport simulations, these models are also a
crucial input for calculations of properties of neutron star crusts (see
Section 2.3).

Additional theoretical constraints might be provided by high-density
calculations within the framework of perturbative QCD (pQCD) [311], which can
be applied at very high densities of the order of 40 times the nuclear
saturation density, where the strong interactions among quarks become
perturbative. Constraints on the EOS based on pQCD, together with assumptions
on causality and stability, have been used to constrain the EOS at lower
densities probed in the core of neutron stars [312, 313, 314, 315]. However,
it has been found that the constraining power of pQCD calculations is strongly
dependent on the way in which they are implemented [315, 316]. Future studies
have to establish to what extent pQCD constraints are robust at densities of
the order of several times nuclear saturation density, and how constraining
future higher-order calculations may become. In this regard, improved
microscopic calculations of the nuclear EOS using the functional
renormalization group [317, 318] will provide important insights.

### 2.3 Neutron star theory

#### 2.3.1 Status

Measurements of the EOS, masses of neutron-rich isotopes far from the band of
stability, and experimental constraints on nucleon effective masses provide
essential input into neutron star models, progressing our understanding of the
structure and dynamics of these astronomically important objects. Several
properties of neutron stars, including the mass-radius relation and their
tidal deformabilities, can be calculated once the EOS is provided. This, in
turn, enables us to constrain the EOS once those properties are observed
[319].

Nuclear EOSs for neutron stars can be constructed from, for example, ab initio
calculations and density functionals [320, 321, 322, 323, 324, 325, 326] or,
more schematically, from meta-models [327, 328, 329] parametrized by nuclear

Figure 2.6: Impact of nuclear physics theory and experiment, and different
astrophysical measurements on constraining the cold neutron-star EOS. Blue
lines show a family of EOS that are constrained by chiral EFT at low
densities. At higher densities, the EOS can then be constrained using GWs from
inspirals of neutron star mergers, data from radio and X-ray observations of
pulsars, and electromagnetic signals associated with neutron star mergers. The
indicated boundaries between regions affected by these measurements are not
strict and depend on the EOS and properties of the astrophysical system.
Figure from [330].

matter parameters, which can be used to make contact with heavy-ion collisions
[67]. Ab initio calculations take into account more fundamental properties of
the nuclear force (see Section 2.2), but prohibit the calculation of large
ensembles of EOSs spanning the nuclear parameter space. Meta-models allow
rapid computation of such large ensembles, but encode mainly bulk properties
of nuclear matter, which excludes them from being used to model finite nuclei.
Density functionals represent a compromise, allowing both rapid computation of
EOSs and use in finite nuclear models, and thus are more suited to combining
nuclear experimental and astrophysical information. Many of these models can
be smoothly extrapolated from the saturation-density to arbitrarily high
density, in which case astronomical observations can be used to constrain the
saturation-density nuclear matter parameters and their density dependence
[329, 331, 332]. This extrapolation, however, is model-dependent, as different
density functionals have different dependence on density. Additionally, this
extrapolation might not be physically well-founded.

As densities inside neutron stars can reach up to several times nuclear
saturation density, at some (as-yet not determined) density a description in
terms of purely nucleonic degrees of freedom is expected to break down. Heavy-
ion collisions can help us constrain that point, and the nature of any phase
transitions that occur above saturation density. The nuclear EOSs can be then
combined with models describing the EOS at higher densities. Models that
explicitly include a range of possible high-density degrees of freedom, such
as hyperons and quarks, can be constructed; the predicted neutron star
compositions are then dependent on the particular model used. Another approach
is to use more general models that give up the explicit dependence on the
underlying degrees of freedom, thus losing information on, e.g., appearance of
exotic particles at high densities, in favor of spanning the full space of
physically consistent EOSs, reducing the model dependence of inferences from
astrophysical observations [333]. These schemes include piecewise polytropes
[334, 335, 44, 336], line segments [337, 282], speed-of-sound models [92, 279,
338, 339, 340], spectral models [341] and non-parametric models generated from
Gaussian processes (GPs) [342, 343, 344, 89, 345] or machine learning
techniques [346]. If these more general approaches are used down to the
nuclear saturation density, extra modeling is required to connect them to the
microscopic nuclear EOS and nuclear observables [347].

Once the EOS is specified, the solution of the Tolman–Oppenheimer–Volkoff
equations and their extensions including rotation, determining the structure
of a neutron star through balancing the attractive force of gravity and the
repulsion coming from the EOS, provide predictions for bulk properties of the
neutron star such as radii, tidal deformabilities, moments of inertia, and
break-up frequencies of neutron stars as a function of their mass. All of
these properties can be compared with multi-messenger observations, including
gravitational waves and electromagnetic signals from neutron-star mergers and
isolated neutron stars [283].

The systematic construction of neutron star EOS models and statistical
inference of EOS parameters from data is an endeavor that is just over a
decade old [334, 335, 44, 336]. This effort has matured in the current era of
multi-messenger astronomy with a large push to explore the model-dependence of
EOS inferences [338, 348] and ways of connecting the EOS with astrophysical
and nuclear data [340, 283, 339, 88, 259, 67]. Different choices of which
observables to include or infer can be made. For example, astrophysical
observations can be used to infer the EOS, which can then be connected to
nuclear models to inform their parameters and predict nuclear observables.
Conversely, nuclear observables can be used to infer nuclear parameters, which
can then inform the neutron star models and predict astrophysical observables.
The future lies in combining more and more sets of data of both types to
understand nuclear and neutron star models better.

Exciting progress has been made in gathering astrophysical data to constrain
our dense matter theories (see Fig. 2.6 for an illustration of density regions
affected by different observables). Neutron-star data from the last 5 years
identified the heaviest neutron star known to date with a mass of
$2.08(7)M_{\odot}$ [86, 87] (where $M_{\odot}$ is the solar mass), excluding
EOSs which cannot reach that limit, while the kilonova AT2017gfo, associated
with GW170817, has placed an upper limit on the maximum mass to be on the
order of $2.3M_{\odot}$ [349, 350]. The detection of GW170817 by the LIGO-
Virgo Collaboration has enabled us to place constraints on the tidal
deformability of this system, $\tilde{\Lambda}_{\rm GW170817}\leq 720$ [351,
352]. Neutron Star Interior Composition Explorer Mission (NICER) has provided
two mass-radius measurements by observing X-ray emission from several hot
spots on the neutron star surface, finding a radius of
$13.02^{+1.24}_{-1.06}$km for a star with mass $1.44^{+0.15}_{-0.14}M_{\odot}$
(PSR J0030+0451) and $13.7^{+2.6}_{-1.5}$km for a star with mass
$2.08(7)M_{\odot}$ (PSR J0740+6620) in the analyses of Refs. [353, 88, 354,
355]. X-ray observations of the temperature of the neutron star in the Cas A
supernova remnant have revealed core cooling on the timescale of years,
hinting at the possible superfluid properties of the core [356].

These observations have enabled meaningful constraints on the EOS to be set
and have already allowed us fascinating glimpses into the possible properties
of high-density matter. For example, perturbative QCD predicts that the speed
of sound squared approaches the conformal limit of $c_{\rm s}^{2}=1/3$ from
below as the density becomes arbitrarily high. Meanwhile, inferences of the
neutron star EOS from observational data indicate that the speed of sound
rises in the core to significantly above $1/3$ [90, 91, 92, 93, 94].
Consequently, this suggests that the speed of sound has a non-trivial behavior
with increasing density [311, 92, 357], first rising above $1/3$ to reach a
maximum at moderate densities, then decreasing, and finally approaching $1/3$
from below at asymptotically high densities. Several microscopic mechanisms
for this behavior have been suggested, including the appearance of quarkyonic
matter [29, 358], hadron-quark continuity within a topological Skyrmion model
[359], or superconducting quark matter [360, 361, 362]. At the same time,
tentative evidence for a softening of the EOS in neutron star cores, based on
an analysis using a novel speed-of-sound interpolation, has likewise been
suggested [340].

If we want to leverage the substantial data we have on neutron star cooling
and dynamical evolution, additional EOS quantities need to be supplied
consistently for each EOS model, such as the effective masses (see also
Section 6.5) and superfluid neutron and proton gaps, essential for modeling
thermal and dynamical properties of neutron stars. For example, the mutual
friction of the core – the strength of the coupling between the charged
particles (electrons, protons) and superfluid neutrons – depends on the
effective neutron mass and the proton fraction [363], which both also
correlate with the symmetry energy [196]. A consistent extraction of both
symmetry energy parameters and effective masses from heavy-ion collision data
is therefore required.

In contrast to efforts devoted to systematic, statistically meaningful
inferences of the EOS in the cores of neutron stars, modeling the neutron star
_crust_ is still in its infancy: The first calculations of large ensembles of
systematically parametrized crust models and their use in statistical analysis
have only been carried out recently [364, 365, 366, 367, 368]. However, much
more nuclear experimental data can be brought to directly bear on crust
physics, and we have entered an era where we can access information about the
crust with unprecedented fidelity. For example, we have now observed the same
neutron-star crust as it first cooled, then became heated by accreted matter,
and then cooled again [369, 370, 371, 372, 373]. We have followed a pulsar
through a glitch – a sudden change in the rotation period of the pulsar – and
glitch recovery with a resolution of a few seconds [374]. These observations
have provided very strong evidence that the crust is solid, that there exist
superfluid neutrons in the inner crust which can be decoupled from the nuclei
in the crustal lattice, and that nuclear reactions from accreted material
sinking into the crust provide deep crustal heating [371, 375, 376].

Additionally, models of the neutron star crust predict that, prior to the
transition to homogeneous matter, isolated nuclei in the crust fuse to form
cylindrical, planar, and more exotic shapes, termed “nuclear pasta”, that can
affect neutron-star observations [377, 378, 379, 380]. This crust-core
boundary region, often referred to as the mantle, is likely a complex fluid.
Density functional theory and molecular dynamics calculations of these
structures reveal a complex energy landscape with many coexisting shapes, and
correspondingly complex mechanical and transport properties [381, 382, 383,
31, 384, 385, 386], which are strongly influenced by the EOS at around
0.5$n_{0}$ through the pressure, proton fraction, and surface energy of the
structures. These properties can also be studied in multifragmentation
reactions, which probe, among others, the competition between nuclear surface
energy and Coulomb energy at sub-saturation density [387, 388, 389].

Inhomogeneous matter in the crust of a neutron star, including the dripped
neutrons expected in the inner crust, can be modeled using a variety of
nuclear theory techniques. These usually involve calculations within a single,
repeating unit (Wigner–Seitz cell) of matter, typically containing a single
nucleus [390, 391, 392]. The compressible liquid drop model (CLDM) treats the
nuclear matter inside and outside of nuclei as homogeneous and described by
the bulk matter EOS, while the surface energy is specified by a separate
function with additional parameters [393, 381, 394, 395, 392]. The surface
parameters and those that define the dimensions of the cell and nucleus are
minimized to obtain the ground state. The Thomas-Fermi model employs the local
density approximation, modeling matter with a specified form of the
inhomogeneous nuclear matter density in the unit cell; here, the parameters of
the density distribution are varied to obtain the ground state configuration
[396]. Microscopic approaches to describing inhomogeneous nuclear matter, in
which individual neutrons and protons are the degrees of freedom, include
quantum Hartree-Fock or Relativistic Mean Field models [397, 398, 399, 400,
401, 402], and semi-classical molecular dynamics approaches [403, 384].

There is a great need for nuclear physics input into models of the neutron
star crust, which analyses of heavy-ion collision data can provide. For
example, the thickness, mass and moment of inertia of the crust depend on the
higher-order symmetry energy parameters $L$, $K_{\rm sym}$, and $Q_{\rm sym}$
[404, 364, 366]. Thus measurements of the symmetry energy parameters up to
third order in heavy-ion collision experiments are essential to understand the
properties of the crust. The symmetry energy, effective masses, and surface
energies of nuclear clusters strongly affect the proton fraction on either
side of the crust-core transition density, the extent of nuclear pasta near
the crust-core boundary, the mechanical and transport properties, the thermal
conductivity and specific heat, the electrical conductivity, and the shear
modulus of the crust [396, 390, 405, 401]. Nuclear experiment can thus
constrain neutron star crust models, and astrophysical observables associated
with the crust can measure nuclear observables as well as measurements of
neutron star bulk properties. For example, the symmetry energy can be
constrained by combining nuclear data with crust and core observables, e.g.,
through a potential multi-messenger measurement of the resonant frequency of
crust-core interface oscillations [368].

#### 2.3.2 Challenges and opportunities

The next decade will provide a wealth of new data on neutron stars, as the
LIGO-VIRGO-KAGRA detectors are expected to observe many new binary neutron-
star mergers, some of them with electromagnetic counterparts [406, 407, 408].
As NICER continues to measure more neutron star masses and radii, next-
generation X-ray timing missions such as Strobe-X [409] and radio telescopes
such as the Square-Kilometer Array will increase the number of pulsars we see
and are able to measure by an order of magnitude. Long-timescale observations
of individual pulsars (using radio timing) and persistent gravitational waves
from deformations of neutron stars will lead to measurements of their moments
of inertia. These new data points might enable us to pin down the nuclear
matter EOS, to discover or rule out the existence of phase transitions to
exotic forms of matter in the cores of neutron stars, and to reliably
constrain microscopic interactions between fundamental particles. Further
ahead, next-generation detectors, such as the Einstein Telescope [410] and the
Cosmic Explorer [411], will be able to observe the neutron-star inspiral phase
as well as the onset of tidal effects with high signal-to-noise ratio, with
direct consequences for the resolving the interior structure of neutron stars
and probing fundamental properties of matter at highest densities.

Although model-agnostic extrapolations to higher densities such as through the
use of polytropes [334, 253], speed of sound schemes [92, 338, 412], Gaussian
processes [344, 343] and spectral methods [341], combined with robust data
analysis, will eventually allow us to pin down the dense-matter EOS, they
cannot answer the question about the relevant microscopic degrees of freedom
at high densities. Hence, it is crucial to develop improved microscopic models
with well-quantified uncertainties in this regime. At the same time, creating
ensembles of outer core and crust models that allow for inclusion of
astrophysical and nuclear data requires underlying nuclear models to have
enough freedom to explore a large region of parameter space, and allow fast
computation of relevant quantities that also capture the essential physics.
Currently, it is energy density functionals like Skyrme, Gogny, and
Relativistic Mean Field models that provide these properties. Consequently,
progress could be made by making a stronger connection between these models
and microscopic approaches, e.g., connecting energy-density functionals to ab
initio calculations allowing a more direct link to $\chi$EFT [391, 413, 414].
In the same spirit, EFT calculations of the EOS can be used as a “low-density
limit” to calibrate higher-density models for neutron stars and heavy-ion
collisions.

The crust can be modeled consistently with nucleonic matter in the core using
density functional theory to model both. When choosing a model, a compromise
must be made between accurate modeling of microscopic quantum effects, such as
shell effects in the nucleus and surrounding neutron gas, and the
computational expediency required to construct large ensembles of crust models
needed for statistical inference. For example, quantum shell effects strongly
determine the evolution of the mass and charge number of nuclei with density,
alter the effective mass of dripped neutrons, and drive the complex energy
landscape of nuclear pasta. Fully microscopic quantum calculations include
shell effects self-consistently, but are computationally expensive. The CLDM
approach can be used to construct large numbers of crust models, but requires
shell effects to be added by hand. Future work needs to develop schemes of
incorporating such microscopic effects in large ensembles of crust models. The
method that may allow that is the Extended Thomas-Fermi method, incorporating
shell effects through the Strutinsky Integral (see, e.g., [415]).

Models should also incorporate nuclear pasta, as its extended structures may
contribute to the mechanical and thermal properties of matter at the crust-
core boundary. It is computationally demanding to model transport and
mechanical properties of the crust microscopically or in simulations [416],
particularly in the nuclear pasta phases, and it is unrealistic to include