††thanks: Corresponding author

# Correlations of Conserved Quantities at Finite Baryon Density

Oleh Savchuk<EMAIL_ADDRESS>Facility for Rare Isotope Beams, Michigan
State University, East Lansing, MI 48824 USA  Bogolyubov Institute for
Theoretical Physics, 03680 Kyiv, Ukraine  Scott Pratt  Department of Physics
and Astronomy, Michigan State University, East Lansing, MI 48824 USA  Facility
for Rare Isotope Beams, Michigan State University, East Lansing, MI 48824 USA

(September 1, 2024)

###### Abstract

Correlations involving the seven conserved quantities, namely energy, baryon
number, electric charge, strangeness, and the three components of momentum,
give rise to correlations in heavy-ion collisions. Through the utilization of
a simple one-dimensional hydrodynamic model, we calculate the evolution of the
entire $7\times 7$ matrix of correlations as a function of relative spatial
rapidity. This comprehensive analysis accounts for finite baryon density,
which results in off-diagonal correlations between the charge-related
quantities and the energy-momentum quantities. These correlations in
coordinate space are subsequently transformed into correlations in momentum
space using statistical weighting. The entire matrix of correlations is
revealed to be highly sensitive to the equation of state (EoS), viscosity, and
diffusivity.

##  I Introduction

The properties of QCD matter have been a major focus of the nuclear physics
community [1, 2]. Nucleus-nucleus (A+A) collisions provide the majority of
experimental information about the bulk properties of the matter reflected in
its phase diagram. Particularly, the existence phase transitions (PTs) related
to the change of hadronic to quark and gluon degrees of freedom have become a
focal point of the field. Recently, the properties of cold and dense nuclear
matter reached in low energy collisions has become increasingly relevant in
the context of binary neutron star mergers. Observations related to the radii
of neutron stars point to a remarkably stiff equation of state for baryon
densities reaching a few times nuclear saturation density, with the speed of
sound possibly exceeding half the speed of light [3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15]. Fluctuation and correlation observables have always played a
prominent role in representing the bulk properties of matter, particularly
near phase transitions or near critical points [16, 17, 18, 19]. Even for
featureless equations of state fluctuations represent fundamental physical
quantities. For example, the speed of sound can be expressed in terms of
fluctuation observables [20].

Unfortunately, the matter created in A-A collisions evolves and disassociates
in a short time. The ephemeral nature of the produced matter impedes the
opportunity for longer-range correlations to equilibrate. Thus, in the context
of heavy-ion collisions, the dynamical evolution of correlations must be
understood. Whereas in a macroscopic system with infinite time to equilibrate,
the charge within any finite volume fluctuates, in heavy-ion collisions one
needs to carefully specify the volume over which the fluctuation is being
considered. For conserved charges, there is zero fluctuation in the limit that
one considers the entire collision volume.

To understand correlations and fluctuations from a dynamic perspective, one
must understand how correlations are seeded, and how they spread with time.
For conserved charges, such as baryon or electric charge, this is typically
quantified through charge balance correlations [21, 22]. In coordinate space
such correlations can be expressed as

$\displaystyle C(\vec{r}_{1},\vec{r}_{2},t)$ $\displaystyle=$
$\displaystyle\langle[\rho(\vec{r}_{1},t)-\bar{\rho}(\vec{r}_{1},t)][\rho(\vec{r}_{2},t)-\bar{\rho}(\vec{r}_{2},t)]\rangle.$
(1)

If the correlations appear equilibrated at short relative distance, one can
divide the correlation into two pieces,

$\displaystyle C(\vec{r}_{1},\vec{r}_{2},t)$ $\displaystyle=$
$\displaystyle\chi(\vec{r}_{1})\delta(\vec{r}_{1}-\vec{r}_{2})+C_{B}(\vec{r}_{1},\vec{r}_{2},t).$
(2)

Here, the delta function is not meant literally, but only to represent some
function of short range that integrates to unity. If a system is equilibrated,
and if those correlations are not longer-range, the strength of the local
correlation is set by the charge susceptibility, which is a fundamental bulk
property of the matter. The second term, the balancing correlation, would
spread over the entire available volume if given the opportunity. However, its
spread is constrained by the amount of time elapsed since the correlations
were seeded. The integrated strength of the balancing correlation must reflect
the fact that for the entire volume, charge does not fluctuate, and

$\displaystyle\int d^{3}r\leavevmode\nobreak\
C_{B}(\vec{R},\vec{R}+\vec{r},t)$ $\displaystyle=$
$\displaystyle-\chi(\vec{R},t).$ (3)

If the system has translational invariance, the $\vec{R}$ dependence can be
ignored. At the end of the collision, one can measure the balancing
correlation, i.e. it is the correlations of all particles with those besides
themselves. At high energy, where the average charge densities are near zero,
$C_{B}(\vec{r},t)$ typically spreads diffusively. Local charge conservation
demands that the source function for the balance correlations, in the matter
frame, is given by

$\displaystyle S(\vec{R},t)$ $\displaystyle=$
$\displaystyle-(\partial_{t}+\nabla\cdot\vec{v})\chi(\vec{R},t).$ (4)

Thus, if one measures the correlations at latter times, and if one understands
how charges diffuse with time, one can infer the source function for the
correlations. One can then infer, not only the susceptibility at a specific
time, but the evolution of the susceptibility over time. Hydrodynamic
calculations provide temperatures, flow velocities and densities as a function
of time. This then provides the susceptibilities, which are functions of
density and temperature, as a function of time, which then determines the
source functions for $C_{B}$, and if one were to know the diffusivity, one can
then calculate the evolution of $C_{B}$. Indeed, such calculations have been
performed, overlaid on a three-dimensional hydrodynamic background, using the
equation of state, charge susceptibility, and diffusivity from lattice
calculations [23, 24, 25, 26].

In heavy-ion collisions, measurement is confined to recording the momenta of
final-state particles. Fortunately, collisions are highly explosive, and the
large flows provide the means to map correlations in relative velocity to
those in coordinate space, albeit with some smearing due to mainly thermal
motion. The corresponding correlations in relative velocity have been
measured, and are known as charge-balance functions. Such analyses, which can
be binned by any measure of the relative velocity, such as the relative
momentum or relative azimuthal angle, have proven remarkably adept at
constraining both the evolution of the charge susceptibilities as a function
of time, and the diffusivity of QCD matter. Remarkably, the analyses support
the notion that matter chemically equilibrates early in the collision, and
that baryon number, strangeness and electric charge diffuse at rates
calculated in lattice simulations [27].

The first goal of this paper is to understand how the above paradigm can be
extended to include other conserved quantities, i.e. energy and momentum. The
second goal is to understand how at non-zero average charge density, which is
increasingly the case for lower-energy collisions, energy and momentum
correlations mix with those indexed by charge. Whereas the first goal was
considered in [24], the effects of non-zero density were ignored. One can
define a $7\times 7$ susceptibility matrix, involving energy, three components
of momentum, and the three charges, baryon number, electric charge, and
strangeness. That matrix has cross terms between energy and the three charge
densities. For example, the fluctuation between energy and baryon density,
$\langle\delta E\delta B\rangle$, is non-zero if there is a finite net baryon
density. To make matters more complicated, energy and momentum correlations do
not spread diffusively, but due to viscous hydrodynamics. Energy and momentum
correlations then mix through the hydrodynamic evolution equations, regardless
of whether the net baryon density is zero. The next goal is to show how
correlations involving energy and momentum in coordinate space can be
projected onto final-state hadrons at the decoupling by providing a thermal
sampling of particles, with each particle accompanied by an additional weight
factor. Finally, as a fourth goal, we wish to see how the corresponding
$7\times 7$ matrix of correlation functions might be sensitive to variations
of the equation of state.

The model to be applied here is simple, longitudinally-boost-invariant
hydrodynamics without transverse flow. Hereafter, this will be referred to as
the Bjorken limit [28]. The neglect of transverse flow is not warranted, at
any energy, and boost invariance is poorly justified at lower energies.
However, this simple picture should be sufficient for understanding the issues
mentioned above. Correlations involving energy and momentum are likely to be
strongly affected by inhomogeneities of the initial state, and by jets. Even
crude estimates of the strength of correlations due to local
charge/energy/momentum conservation are valuable in that one can then
ascertain whether further, more realistic, modeling is warranted, or if such
correlations are so small that they are overwhelmed by the competing effects
related to the initial state and jets.

The next section of the paper focuses on the underlying theory. It presents
the generalized $7\times 7$ susceptibilities, and the expressions for then
evolving the $7\times 7$ matrix of balancing correlations. The evolution is
based on small perturbations of the uniform hydrodynamic background, which
spread both diffusively and hydrodynamically. Green’s functions are pursued as
a means by which to evolve the correlations from the sources. Expressions for
the evolution on the simple boost-invariant background are presented. The
final subsection of the theory section shows how one can project the spatial
correlations onto a simulation of outgoing particles, which can then be used
to calculate how such spatial correlations might drive the asymptotic final-
state correlations measured in momentum space. The final-state correlations
can be binned by any measure of relative momentum, such as relative rapidity
or azimuthal angle, and might be indexed by such quantities as transverse
energy, transverse momentum, or charge.

In Sec. III, after showing an illustrative example of Green’s functions,
final-state correlations are presented for three different equations of state.
Correlations are presented both in coordinate space, as a function of relative
spatial rapidity, and momentum space as a function of relative rapidity.
Implications of the results are then summarized in Sec. IV.

##  II Theoretical Foundation

This section consists of three parts. First, the connection between the
$7\times 7$ susceptibilities and the correlation functions are elucidated. The
source function, which is driven by the evolution of the susceptibility, is
connected to the final-state correlation through Green’s functions which
describe the response of the medium to small perturbations of the seven
quantities. The following sub-section shows how Green’s functions are
calculated on the background of a boost-invariant hydrodynamics description
lacking transverse flow. Finally, a method is presented for projecting the
correlations onto the momenta of particles in a hadron gas.

###  II.1 General Theory of Conservation-Driven Correlations

Susceptibilities are related to fluctuations, as calculated in the Grand
Canonical Ensemble. In general, one can weight a partition function according
to any conserved quantity,

$\displaystyle Z$ $\displaystyle=$ $\displaystyle{\rm Tr}\leavevmode\nobreak\
\exp\\{-\beta_{i}Q_{i}\\}.$ (5)

Here, Lagrange multipliers, $\beta_{i}$, are associated with each conserved
quantity $Q_{i}$. The seven quantities, $Q_{i}$, are
$E,P_{x},P_{y},P_{z},B,Q,S$, the energy, momenta, baryon number, electric
charge and strangeness. For some large volume $\Omega$, one can define
fluctuations as

$\displaystyle\frac{1}{\Omega}\langle\delta Q_{i}\delta Q_{j}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\Omega}\frac{\partial^{2}}{\partial\beta_{i}\partial\beta_{j}}\ln(Z).$
(6)

The Lagrange multiplier related to the energy is the inverse temperature,
whereas for the three charges they are $-\mu_{i}/T$, where $\mu_{i}$ is the
chemical potential and $T$ is the temperature. For the momenta, the Lagrange
multipliers can be related to $u_{i}/T$, where $u_{i}$ are the flow
velocities. Here, the term susceptibilities is used in a most general sense.
The usual, but not always, definition of charge susceptibilities is to define
them as the derivative of the charge densities with respect to the chemical
potentials. That definition differs from what is used here by a factor of $T$.
In the usage here, all the “susceptibilities” correspond to the fluctuations
without additional factors,

$\displaystyle\chi_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{\Omega}\langle\delta Q_{i}\delta Q_{j}\rangle.$ (7)

For the energy, the susceptibility, $\langle\delta E\delta E\rangle/\Omega$,
is $T^{2}$ multiplied by the specific heat. The $3\times 3$ subset of
$\chi_{ij}$ that refer to baryon number, electric charge and strangeness, are
the usual charge fluctuations. At finite number baryon density, the
correlations between the energy and the three charges becomes non-zero. These
cross terms have not been considered as extensively as correlations involving
only the three charges or only the energy. Finally, one can have correlations
involving the three momentum components. These have no cross terms, and at
equilibrium

$\displaystyle\frac{1}{\Omega}\langle\delta P_{i}\delta P_{j}\rangle$
$\displaystyle=$ $\displaystyle(P+\epsilon)T\delta_{ij}.$ (8)

In equilibrium, the correlations integrate to the susceptibilities. If
$\rho_{i}$ are the corresponding densities,

$\displaystyle C_{ij}(\vec{r}_{i},\vec{r}_{j})$ $\displaystyle=$
$\displaystyle\langle\delta\rho_{i}(\vec{r}_{i})\delta\rho_{j}(\vec{r}_{j})\rangle$
$\displaystyle=$ $\displaystyle\chi_{ij}\delta(\vec{r}_{i}-\vec{r}_{j}).$

In a heavy ion collision, all the quantities are conserved, so for any
operator $A$,

$\displaystyle\langle A\delta Q_{j}\rangle$ $\displaystyle=$ $\displaystyle
0,$ (10) $\displaystyle\int d^{3}r_{j}\langle
A\delta\rho_{j}(\vec{r}_{j})\rangle$ $\displaystyle=$ $\displaystyle 0,$

and for any conserved density

$\displaystyle\int d^{3}r^{\prime}\leavevmode\nobreak\
C_{ij}(\vec{r},\vec{r}^{\prime})$ $\displaystyle=$ $\displaystyle 0.$ (11)

The correlation can be divided into a local, equilibrated, part and a
balancing part,

$\displaystyle C_{ij}(\vec{r}_{i},\vec{r}_{j})$ $\displaystyle=$
$\displaystyle\chi_{ij}\delta(\vec{r}_{i}-\vec{r}_{j})+C_{B,ij}(\vec{r}_{i}-\vec{r}_{j}),$
(12) $\displaystyle\int d^{3}r_{j}C_{B,ij}(\vec{r}_{i}-\vec{r}_{j})$
$\displaystyle=$ $\displaystyle-\chi_{ij}(\vec{r}_{i}).$

If the local part were not equilibrated, the equivalence would still hold
true, but in that case $\chi_{ij}$ would not be given by the equilibrated
susceptibility. Instead, it might be some dynamic quantity that would relax
toward the equilibrated value.

As the system cools into a gaseous state, correlations between particles
disappear. In that limit, the equilibrated local correlation is only that
within a single particle, and the susceptibilities become

$\displaystyle\chi_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{\Omega}\sum_{a\in\Omega}Q_{ai}Q_{aj},$ (13)

where $Q_{ai}$ is the charge of type $i$ on particle $a$. In terms of the
phase space density,

$\displaystyle\chi_{ij}(\vec{r})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3}}\sum_{h}\int d^{3}p\leavevmode\nobreak\
f_{h}(\vec{p},\vec{r})Q_{i}(h,\vec{p})Q_{j}(h,\vec{p}),$ (14)

where the sum extends over hadron species. For example, if $i$ refers to the
momentum component $p_{x}$, $Q_{i}(\vec{p})=p_{x}$ and if $i$ refers to
strangeness, $Q_{i}$ is simply the strangeness of species $h$. Thus, at
decoupling, the balancing correlation, $C_{B,ij}$, covers all correlations
between different particles, assuming the matter is a non-interacting gas at
the time of decoupling.

In the absence of a source function
$\langle\delta\rho_{i}(x_{1})\delta\rho_{j}(x_{2})\rangle$ evolves as
$\langle\delta\rho_{i}(x_{1})\rangle\langle\delta\rho_{j}(x_{2})\rangle$ with
each $\delta\rho$ behaving independently in terms of the space time
coordinates $x_{1}$ and $x_{2}$. To satisfy the sum rule, one can then express
the correlation as

$\displaystyle\langle\delta j^{\mu}_{i}(x_{1})\delta
j^{\nu}_{j}(x_{2})\rangle$ $\displaystyle=$ $\displaystyle\int
d^{4}X\leavevmode\nobreak\
G^{\mu\mu^{\prime}}_{ii^{\prime}}(X,x_{1})G^{\nu\nu^{\prime}}_{jj^{\prime}}(X,x_{2})u_{\mu^{\prime}}(X)S_{i^{\prime}j^{\prime}}(X)u_{\nu^{\prime}}(X).$
(15)

Here, the formalism has been expressed with the charge density replaced by the
four-current. This enables one to project the correlations through a
hypersurface that is not necessarily at constant Euclidean time. For baryon,
electric and strangeness currents, the quantities $j_{i}^{\mu}$ are simply the
usual three charge currents. The energy and momentum current densities are
elements of the stress-energy tensor, $T^{0\mu}$ and $T^{k\mu}$, respectively,
where $k$ denotes the $k^{\rm th}$ component of the momentum density. The
Green’s functions satisfy the normalization conditions,

$\displaystyle\int
d\Omega_{\mu}G^{\mu\mu^{\prime}}_{ii^{\prime}}(X,x)u_{\mu^{\prime}}(X)=\delta_{ii^{\prime}},$
(16)

where $x$ refers to the space-time point along the hyper-surface where the
matter disassociates and $d\Omega_{\mu}$ is the differential volume in the
reference frame where the disassociation is locally simultaneous. The extent
of this hyper surface is fully in the absolute future relative to $X$, and the
integral covers the entire hyper volume. If one defines a step function, such
that inside the hyper volume the value is unity, and outside it is zero, the
differential volume element can be defined as

$\displaystyle d\Omega_{\mu}=d^{4}y\leavevmode\nobreak\
\partial_{\mu}\Theta(C(y)),$ (17)

where $C(y)$ is positive outside the surface and negative inside. For example,
defining $C(y)=T_{0}-T$, defines a hyper-surface separating the regions with
temperature above and below $T_{0}$. In this study, where a simple Bjorken
expansion is considered, space-like hyper-surfaces, like one would see with
evaporation from a surface, are not encountered.

The source functions are given by the rate of change of the susceptibilities.
If the collective velocity is noted by $u^{\mu}(x)$,

$\displaystyle S_{ij}(X)$ $\displaystyle=$
$\displaystyle[u\cdot\partial+(\partial\cdot u)]\chi_{ij}(X).$ (18)

In a small hydrodynamic volume $\delta V$, the rate of change of the volume is
$\delta V\partial\cdot u$, and the source function is $(1/\delta
V)(u\cdot\partial)(\delta V\chi_{ij})$. Thus, if $\chi_{ij}$ falls inversely
with the volume there is no source. This would happen for baryon, electric and
strange charges for the isentropic expansion of a non-interacting gas. If
entropy is conserved, one can identify the source function with

$\displaystyle S_{ij}(x)\approx s(u\cdot\partial)\frac{\chi_{ij}}{s}.$ (19)

where $s$ is the entropy density, with this becoming exact for isentropic
expansions. For non-zero charge densities, one can similarly replace $s$ with
the baryon, electric or strange densities, although this would neglect charge
diffusion.

The principal challenge in calculating correlations is in finding the Green’s
functions. At each point in spacetime, $X$, one considers a small perturbation
of each of the seven quantities. One then writes equations of motion for the
perturbation. These equations mix the various quantities, i.e.
$G^{\mu}_{ii^{\prime}}(X,x)$ is diagonal for times $x_{0}$ immediately
following $X_{0}$ but off-diagonal components ensue thereafter. If the average
charge densities were all zero, the Green’s functions for energy and momentum
do not mix with those for baryon number, electric charge or strangeness. The
Green’s functions for charge then simply require solving the diffusion
equation for an initial charge perturbation at $X$. If the diffusivity matrix
is diagonal, as one would expect for an idealized quark-gluon plasma, the
charge evolutions would not mix the three charges if the were expressed in the
$u,d,s$ basis. For the energy and momentum perturbations, one must address the
evolution with hydrodynamics. Hydrodynamics mixes the energy and momentum
components, and if the average charge densities are non-zero, as is the case
for lower-energy heavy-ion collisions, the Green’s function also develops off-
diagonal components between the energy/momentum and charge components.

###  II.2 Hydrodynamic Response

The Green’s function represents the response of the medium to a small
localized fluctuation. In this study, this is calculated by assigning a small
perturbation of charge type $i$ to a specfic point in space-time. The
perturbation integrates to unity, but is treated in the linear approximation.
Its spatial extent is set to a small, but non-zero, value. In the limit that
the extent is zero, the initial perturbation would be a delta function, but
that would preclude evolution equations involving spatial derivatives. Here,
we first review the treatment of linearized hydrodynamic fluctuations.

The stress-energy tensor of a viscous fluid in the Landau frame is represented
as follows [29]:

$\displaystyle
T^{\mu\nu}=(\varepsilon+P(\varepsilon,\rho_{B},\rho_{Q},\rho_{S}))u^{\mu}u^{\nu}-g^{\mu\nu}P+T_{\eta}^{\mu\nu}.$
(20)

In this equation, $\varepsilon$ and $\rho$ stand for energy and particle
densities, respectively. In this study, the shear contribution
$T^{\mu\nu}_{\eta}$ is assigned according to the Navier-Stokes equation,

$\displaystyle
T^{\mu\nu}_{\eta}=-\eta\left(\nabla^{\mu}u^{\nu}+\nabla^{\nu}u^{\mu}-\frac{2}{3}\Delta^{\mu\nu}\partial\cdot
u\right),$ (21)

where $\eta$ is the shear viscosity coefficient, and
$\Delta^{\mu\nu}=u^{\mu}u^{\nu}-g^{\mu\nu}$ acts as a projector eliminating
the collective velocity $u^{\mu}$, and the covariant derivative
$\nabla^{\mu}=\Delta^{\mu\nu}\partial_{\nu}$ represents the spatial
derivatives in the frame of the fluid. If the charge number is allowed to
diffuse, the charge current can be expressed as:

$\displaystyle
j^{\mu}_{i}=\rho_{i}u^{\mu}-D_{ij}\Delta^{\mu\nu}\partial_{\nu}\rho_{j}.$ (22)

In this equation, $D$ represents the diffusion matrix. The hydrodynamic
equations are equivalent to energy-momentum conservation:

$\displaystyle\partial_{\nu}T^{\mu\nu}$ $\displaystyle=$ $\displaystyle 0,$
(23) $\displaystyle\partial_{\nu}j^{\nu}$ $\displaystyle=$ $\displaystyle 0.$

Assuming that the solution for $\rho(t,\vec{r}),\varepsilon(t,\vec{r})$ and
$u^{\mu}(t,\vec{r})$ can be expanded into a series around the given solution,
with respect to the perturbations $\delta\rho,\delta\varepsilon,\delta
u^{\mu}$, our focus will mainly be on $\delta T^{\mu\nu}$, which is linear in
these perturbations. These perturbations should adhere to the following
equations:

$\displaystyle\partial_{\nu}\delta j^{\nu}$ $\displaystyle=$ $\displaystyle
0,$ (24) $\displaystyle\partial_{\nu}\delta T^{\mu\nu}$ $\displaystyle=$
$\displaystyle 0.$ (25)

To generate the Green’s function, one considers a small very localized
perturbation at space-time point $X$ that evolves over time. One can assign a
small Gaussian perturbation. If the initial perturbation is at space-time
point $x$,

$\displaystyle\delta j^{\mu}_{i^{\prime}}(X)$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi\sigma^{2})^{3/2}}e^{-|Y^{2}|/2\sigma^{2}},$ (26)
$\displaystyle Y^{\mu}$ $\displaystyle=$ $\displaystyle X^{\mu}-(u\cdot
X)u^{\mu}.$

Here, $Y$ is the same as $X$, but with the temporal component, as defined in
the fluid frame, projected away. The width $\sigma$ should be chosen as small
as possible while still enabling differentiation. In the limit of
$\sigma\rightarrow 0$ the initial perturbation becomes a delta function in
coordinate space. For the initial perturbation, one then solves for the
perturbation at a space-time point with time $x$, and the Green’s function is

$\displaystyle G^{\mu}(X,x)_{ii^{\prime}}$ $\displaystyle=$
$\displaystyle\delta j^{\mu}_{i}(x).$ (27)

Given the existence of dissipation, all solutions should asymptotically
approach a uniform solution. This is particularly relevant for thermodynamic
fluctuations in baryon number, momentum, or energy. Following the theoretical
framework outlined in [24], we aim to investigate the correlations of
conserved charges at freeze-out, specifically examining the impact of non-
diagonal susceptibilities in the presence of finite baryon density.

###  II.3 Bjorken Expansion as a Background

To most simply illustrate how correlations evolve during the hydrodynamic
stage of strongly interacting matter, a boost-invariant Bjorken solution is
chosen as the background hydrodynamic model. In this case, the evolution is
simplified by transforming the time $t$ and longitudinal space-coordinate $z$
to what is known as the proper time $\tau$ and the spatial rapidty, $\eta$:

$\displaystyle\quad t=\tau\cosh\eta,\quad z=\tau\sinh\eta.$ (28)

Energy and $B,Q,S$ density is uniform through the system and its time
evolution is given by:

$\displaystyle\partial_{\tau}\varepsilon=-\frac{\varepsilon+P}{\tau}+\frac{4\eta_{s}}{3\tau^{2}},$
(29) $\displaystyle\rho=\rho_{0}\frac{\tau_{0}}{\tau},$ (30)

where $\tau_{0},\rho_{0}$ represents the starting time and density. The four
velocity is

$\displaystyle u^{0}=\cosh\eta,\,u^{z}=\sinh\eta.$ (31)

The system of linear response equations is given by:

$\displaystyle\partial_{\tau}\delta\varepsilon$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau}\delta\left(\varepsilon+P\right)-\partial_{\eta}\frac{\varepsilon+P-\frac{8\eta_{s}}{3\tau}}{\tau}\delta
u^{\eta},$ (32)

$\displaystyle\partial_{\tau}\left(\varepsilon+P-\frac{4\eta_{s}}{3\tau}\right)\delta
u^{\eta}$ $\displaystyle=$
$\displaystyle-\frac{2}{\tau}\left(\varepsilon+P-\frac{4\eta_{s}}{3\tau}\right)\delta
u^{\eta}-c_{\rho}^{2}\frac{\partial_{\eta}\delta\varepsilon}{\tau}-\partial_{\rho}P\frac{\partial_{\eta}\delta\rho}{\tau}+\frac{4\eta_{s}\partial_{\eta}^{2}\delta
u^{\eta}}{3\tau^{2}},$ (33)

$\displaystyle\partial_{\tau}\delta\rho$ $\displaystyle=$
$\displaystyle-\frac{\delta\rho}{\tau}+\frac{\rho+D\partial_{\tau}\rho}{\tau}\partial_{\eta}\delta
u^{\eta}+\frac{D}{\tau^{2}}\partial_{\eta}^{2}\delta\rho,$ (34)

$\displaystyle\partial_{\tau}\left(\varepsilon+P-\frac{4\eta_{s}}{3\tau}\right)\delta
u^{x,y}$ $\displaystyle=$
$\displaystyle-\frac{2}{\tau}\left(\varepsilon+P-\frac{4\eta_{s}}{3\tau}\right)\delta
u^{x,y}+\frac{4\eta_{s}\partial_{\eta}^{2}\delta u^{x,y}}{3\tau^{2}},$ (35)

where $\delta u^{\eta}=\frac{\delta u^{z}}{u^{0}}$ is introduced in order to
have boost-invariant system of equations (uniform in $\eta$). These equations
describe the hydrodynamic response to some perturbation given by the
stochastic current $j_{\varepsilon},\vec{j}_{\vec{P}},j_{B},j_{Q},j_{S}$. This
system also conserves laboratory frame energy and momentum:

$\displaystyle\partial_{t}T^{tt}+\partial_{z}T^{tz}$ $\displaystyle=$
$\displaystyle\left(\partial_{\tau}+\frac{1}{\tau}\right)(\cosh\,\eta
T^{tt}-\sinh\eta\,T^{tz})$
$\displaystyle+\frac{1}{\tau}\partial_{\eta}(-\sinh\eta\,T^{tt}+\cosh\eta\,T^{tz})=0,$
$\displaystyle\partial_{t}T^{tz}+\partial_{z}T^{zz}$ $\displaystyle=$
$\displaystyle\left(\partial_{\tau}+\frac{1}{\tau}\right)(\cosh\,\eta
T^{tz}-\sinh\eta\,T^{zz})$
$\displaystyle+\frac{1}{\tau}\partial_{\eta}(-\sinh\eta\,T^{tz}+\cosh\eta\,T^{zz})=0,$

Furthermore, $\rho_{\tau}$ and $\rho_{\eta}$ are defined as follows:

$\displaystyle\rho_{\tau}$ $\displaystyle=$ $\displaystyle\cosh\,\eta
T^{tt}-\sinh\eta\,T^{tz}$ $\displaystyle=$
$\displaystyle\delta\varepsilon\cosh\eta+\left(\varepsilon+P-\frac{4}{3\tau}\eta_{s}\right)\delta
u^{\eta}\sinh\eta,$ $\displaystyle\rho_{\eta}$ $\displaystyle=$
$\displaystyle\cosh\,\eta T^{tz}-\sinh\eta\,T^{zz}$ $\displaystyle=$
$\displaystyle\delta\varepsilon\sinh\eta+\left(\varepsilon+P-\frac{4}{3\tau}\eta_{s}\right)\delta
u^{\eta}\cosh\eta.$

Together with $\rho_{B},\rho_{Q},\rho_{S}$, they satisfy the conservation law:

$\displaystyle\int\tau d\eta\rho(\tau,\eta)=\text{const}.$ (38)

In the following quantities:

$\displaystyle\delta P_{\eta}$ $\displaystyle=$
$\displaystyle\left(\varepsilon+P-\frac{4\eta_{s}}{3\tau}\right)\delta
u^{\eta},$ (39) $\displaystyle\delta P_{x,y}$ $\displaystyle=$
$\displaystyle\left(\varepsilon+P-\frac{4\eta_{s}}{3\tau}\right)\delta
u^{x,y},$

will be called longitudinal and $x,y$ momentum respectively. Using Green’s
functions of these equations, correlators can be expressed as follows:

$\displaystyle c_{AB}$ $\displaystyle=$ $\displaystyle\int
d\tau_{j}\tau_{j}d\eta_{j}G_{AA^{`}}(\eta_{1}-\eta,\tau_{1},\tau_{j})G_{BB^{`}}(\eta_{2}-\eta,\tau_{2},\tau_{j})\left(\partial_{\tau_{j}}+\frac{1}{\tau_{j}}\right)\chi_{AB}(\tau_{j}),$
(40)

where $A,B$ can be any of the following: $\delta E,\delta P_{\eta},\delta
P_{x,y},\delta B,\delta Q,\delta S$. The term $c_{AB}$ in Equation (40) does
not encompass the portion of correlation that should exist in the equilibrated
hydrodynamic medium right from the start. This initial correlation can take on
any form but is nevertheless bound by the overall conservation of charge, as
outlined in Eqs. (10). In the subsequent discussion, we will exclude this
initial correlation, and discuss it as a distinct element that can be
independently studied.

###  II.4 Projecting Correlations onto Final-State Hadrons

The previous subsections describe how one might find the correlations
$\langle\delta j^{\mu}_{i}(x)\delta j^{\nu}_{k}(x^{\prime})\rangle$, which are
functions of space-time. These correlations must then be projected onto final-
state particles. Techniques for this have been applied for purely charge
fluctuations [30, 31], but for this problem we need to extend those ideas to
include correlations involving momentum and energy. As was done with the case
with charges, the techniques will be based on the Cooper-Frye equation, where
emission from the hyper-surface element for uncorrelated particles from a
small hyper-surface element $\delta\Omega_{\mu}$ is given by

$\displaystyle\delta dN_{h}$ $\displaystyle=$
$\displaystyle(p\cdot\delta\Omega)\frac{d^{3}p}{E_{p}}f_{h}(\vec{p},x),$ (41)

where $f_{h}(\vec{p},x)$ is the phase space density of hadrons of type $h$ and
momentum $\vec{p}$ at space-time point $x$, with $x$ being at the hyper-
surface element. An efficient algorithm for generating particles consistent
with a thermal phase space density, including viscous corrections, is
described in [32]. The approach here will be to first integrate over each
source point $X$. For each point $X$, one then generates completely
uncorrelated pairs of hadrons. The two hadrons are uncorrelated with the
source point $X$ and are uncorrelated with one another. The correlation
function in momentum space is constructed by incrementing the correlation
functions by a correlation weight.

To calculate the weight, one can consider a small amount of energy, momentum,
and charge, $\delta Q_{i}$, that passes through a hyper-surface element. For a
thermal distribution, the phase-space density is altered by modifying the
Lagrange multipliers associated with each conserved quantity,

$\displaystyle f_{h}(\vec{p},x)$ $\displaystyle=$ $\displaystyle
f_{h}^{(0)}(\vec{p},x)\exp(-\delta\beta_{i}q_{i}(h,\vec{p})).$ (42)

Here, $i$ refers to the seven generalized charges of the particle of species
$h$ and given energy and momentum, as described at the beginning of this
section. The quantity $q_{i}(h,\vec{p})$ are the seven generalized charges of
a single hadron. The seven components of $\delta\beta$ are determined by
fixing the seven quantities, $\delta Q_{i}$, which are the charges passing
through the hyper-surface volume $V$.

$\displaystyle\delta Q_{i}$ $\displaystyle=$
$\displaystyle\sum_{h}\int\frac{d^{3}p}{(2\pi)^{3}E_{p}}(p\cdot
d\Omega)q_{i}(h,\vec{p})\delta f_{h}(\vec{p},x)$ $\displaystyle=$
$\displaystyle-\sum_{h,j}\int\frac{d^{3}p}{(2\pi)^{3}E_{p}}(p\cdot
d\Omega)q_{i}(h,\vec{p})q_{j}(h,\vec{p})f_{h}^{(0)}(\vec{p},x)\delta\beta_{j}.$

If the frame of hyper-surface dissolution is the fluid frame, can identify
$p\cdot d\Omega/E_{p}$ as a volume $V$, and the r.h.s. can be identified with
the susceptibilities,

$\displaystyle\delta Q_{i}$ $\displaystyle=$
$\displaystyle-V\chi_{ij}\delta\beta_{j},$ (44) $\displaystyle\delta\beta_{i}$
$\displaystyle=$ $\displaystyle-\chi^{-1}_{ij}\delta
Q_{j}/V=-\chi^{-1}_{ij}\delta\rho_{j}.$

Inserting this into Eq. (42), one then has the alteration of the phase space
density

$\displaystyle\delta f_{h}(\vec{p},x)$ $\displaystyle=$ $\displaystyle
f_{h}^{(0)}\left[-q_{i}(h,\vec{p})\chi_{ij}^{-1}(x)\delta\rho_{j}\right].$
(45)

Everything in the square brackets can be considered as a weight due to the
additional generalized charge density $\delta\rho_{j}$.

For each differential space-time volume of the source, one can then generate
two independent particles according to the phases-space densities at the
hyper-surface, $f^{(0)}(\vec{p},x)$. One then assigns a weight to the pair,

$\displaystyle w(\vec{p}_{1},x_{1},\vec{p}_{2},x_{2})$ $\displaystyle=$
$\displaystyle\int d^{4}X\leavevmode\nobreak\
q_{i}(h_{1},\vec{p}_{1})\chi_{ij}^{-1}(x_{1})d\Omega_{1,\mu}G^{\mu\mu^{\prime}}_{j^{\prime}j}(X,x_{1})u_{\mu^{\prime}}(X)$
$\displaystyle
S_{j^{\prime}k^{\prime}}(X)u_{\nu^{\prime}}(X)G^{\nu\nu^{\prime}}_{k^{\prime}m}(X,x_{2})\chi_{mn}^{-1}(x_{2})d\Omega_{2,\nu}q_{n}(h_{2},\vec{p}_{2}).$

One must decide whether to perform the integral over $d^{4}X$ with Monte Carlo
methods. For the case considered in the next section, one can exploit the
symmetries, that the source function depends only on the proper time $\tau$
and that the Green’s functions depend on the relative rapidities, to simplify
the integral. But, if these approaches were to be applied more generally, the
integral above would require; a much more nuanced approach.

One advantage of this approach is that weighted pairs are generated by Monte
Carlo, so it is straight forward to incorporate the effects of decays.
Hadronic decays can be simulated with the same weight from the parents being
assigned to the daughter particles. All the pairs of daughters, with one
daughter chosen from each parent, would be assigned the same weight as the
parental pair. It is also rather straight-forward to consider a variety of
binnings of the correlation. For example if one were to calculate correlations
in transverse energy as a function of relative pseudo-rapidity, one would
increment the bin with that relative pseudo-rapidity by the weight multiplied
by the product of the two transverse energies.

##  III Results

###  III.1 Equations of State

As an equation of state (EoS), we have used the free energy model for a
hadronic mixture (details can be found in the appendix). The background
charges were fixed to reproduce $Q/B=0.4$, as in heavy lead or gold nuclei,
and $S=0$. All trajectories start at $\varepsilon=3\leavevmode\nobreak\
\text{GeV}\text{fm}^{-3}$ and evolve for $10\leavevmode\nobreak\ \text{fm}/c$.
The final temperature falls within $150\pm 2$ MeV. The initial baryon density
was chosen at $8n_{0}$, where $n_{0}=0.16\leavevmode\nobreak\ \text{fm}^{-3}$
represents normal nuclear density, and the viscosity lies in the range of
$(6-8)\frac{s}{4\pi}$, as expected for the hadronic medium. This keeps the
temperature of the medium below $190$ MeV, within the range where hadronic
type of the equation of state and transport can be hypothetically feasible.
The diffusion coefficient $D$ was chosen to be lower than the values predicted
from lattice QCD at vanishing baryon charge, aligning more with hadronic phase
simulations, with $D=\frac{1}{8\pi}T$. However, in principle, the non-diagonal
diffusion matrix should be applied to charges, as diffusion predicted from
different kinetic models is strongly influenced by particle mass [33].
Electric charge is expected to be carried mainly by light mesons, and baryon
charge by heavier particles. Strangeness is interesting as it includes heavy
mesons and baryons, and it is thus correlated with baryon number. Therefore,
one might expect a strong sensitivity of the diffusive properties to the
hadronic composition of the medium.

Figure 1:  The speed of sound is displayed for a non-interacting (ideal)
hadron resonance gas and a hadronic fluid with two illustrative cases for
interactions. The equation of state that includes interactions can be either
stiffer (Hard-EoS) or softer (Soft-EoS) compared to the non-interacting one
(Ideal HRG).  Figure 2: The susceptibilities of energy-energy (a), energy-
baryon density (b), baryon-baryon densities (c) and momentum-momentum (d) are
shown for the three equations of state that were considered. The difference
between different EoS is reflected in the susceptibilities of conserved
charges.

Using the formalism in Appendix B, two equations of state were constructed
using simple parametric forms for the density dependence of the free energy.
The two interactions were not meant to be particularly realistic, but were
constructed to provide the means to compare the difference between a soft and
stiff equation state. With that goal in mind, the equations of state were
constructed to be represent the outer limit of expectations, rather than to
represent some sort of best guess. The free-energy densities were chosen to
depend only on baryon density $\rho_{B}$, and to have the form

$\displaystyle f_{\rm int}$ $\displaystyle=$
$\displaystyle\sum_{n=1,2}A_{n}\left[y^{1/3}-1\right],$ (47) $\displaystyle y$
$\displaystyle=$ $\displaystyle
1+\left(\frac{\rho_{B}-\rho_{0}}{\rho_{n}}\right)^{3}.$

Here, $\rho_{0}$ was set to the break up density so that the equation of
state, including fluctuations, would return to those of a non-interacting
hadron gas at that point. The parameters $A_{1},A_{2},\rho_{1}$ and $\rho_{2}$
were then adjusted to set the stiffness of the equation of state in various
density ranges. For both the stiff and soft equations of state, $\rho_{1}$ and
$\rho_{2}$ were set to 0.2 and 0.5 fm-3 respectively. For the stiff equation
of state the stiffness parameters were set to $A_{1}=0.05,A_{2}=0.1$ GeV/fm3,
and for the soft equation of state $A_{1}=-0.02,A_{2}=0.01$ GeV/fm3. The soft
equation of state has a phase transition with a critical point at $T=100$ MeV
at a density near $2\rho_{s}$. This rather extreme example should should help
bracket the range of fluctuations one might wish to consider.

The initial baryon density was chosen so that the initial baryon density was
eight times nuclear saturation density, $\rho_{s}=0.16$ fm-3, which because
the baryon density fell as $1/\tau$ always yielded the same final baryon
density, $\rho_{0}=(8/11)\rho_{s}$, regardless of the equation of state. The
initial temperature, $T_{i}$ was adjusted for each equation of state, and for
each choice of viscosity, so that the final temperature, $T_{f}$, was 150 MeV.
By construction, the speed of sound matches that of a non-interacting hadronic
mixture at and below $\rho_{0}$. In Fig. 1 the speed of sound is displayed for
the three equations of state. For the non-interacting case the speed of sound
stays fairly steady, with $c_{s}^{2}\approx 0.14$. For the stiff equation of
state the maximum speed of sound for the evolutions considered here was
approximately, $c^{2}_{s,{\rm max}}\approx 0.24$, while for the soft equation
of state the speed of sound dropped to $c^{2}_{s,{\rm min}}\approx 0.11$.

Susceptibilities of a few selected conserved charges can be observed in Fig.
2. Susceptibilities are multiplied by $\tau$ as the volume of a hydrodynamic
cell increases proportionally to $\tau$ in the Bjorken model. Thus, in this
case, the source functions are precisely the rates at which the curves in Fig.
2 rise or fall. The susceptibilities involving energy or baryon number all
have strong peaks or valleys for the same conditions when the speed of sound
has a maximum or minimum.

xRather remarkably, the momentum-momentum susceptibility (multiplied by
$\tau$), $(P+\epsilon)T\tau$, is relatively featureless and is thus rather
insensitive to the equation of state. This probably owes itself to the fact
that the interaction is a function of baryon density, and is of course an
energy, but does not much alter the momenta of particles. It would be
interesting to understand whether a different class of interaction, e.g. one
that alters the degrees of freedom, might manifest itself by providing more
structure for this case. Even if the flatness of momentum-momentum
susceptibility results in little source function, there would still be a
contribution from the initial thermalization, $\tau<1$ fm/$c$. The correlation
of the transverse components, e.g. $\langle\delta P_{x}\delta P_{x}\rangle$,
is expected to spread diffusively [24, 34], with the shear viscosity playing
the role of the diffusion constant. Thus, the correlations of transverse
momentum might be less sensitive to phase structure, but more sensitive to the
viscosity. Equations of state featuring a critical point have an adjacent
mixed phase on the phase diagram ($(\rho,T)$). For the soft equation of state
shown here, phase separation can ensue, but the critical temperature is 100
MeV, and the density-temperature trajectories are well above the critical
point, and therefore outside the mixed-phase region.

Within the mixed phase, or coexistence, region, susceptibilities are
undefined, i.e. fluctuations no longer scale linearly with volume [35, 36].
One might then anticipate the occurrence of nucleation or cavitation,
representing the dynamic formation of a new phase. This process unfolds
relatively slowly and challenges the application of traditional equilibrium
thermodynamics methods like Maxwell constructions. In such cases, a metastable
state can be contemplated, corresponding to the gradual growth of nucleation
or cavitation [37, 38]. Within the spinodal region, mechanical instability
triggers a phenomenon known as spinodal decomposition, where two phases
rapidly segregate. Ideally, this can be incorporated into our analysis by
integrating the appropriate hydrodynamic background. Specifically, the
inclusion of surface tension and energy terms in hydrodynamics can provide an
accurate description of the mixed phase (for further insights, see [39]). The
approach here would then certainly not apply if the density-temperature
trajectory were to pass through the coexistence region of the phase diagram.
Nevertheless, it is worth noting that substantial structure in
susceptibilities are expected, even if merely passes the vicinity of the phase
transition region.

In equilibrated systems correlation lengths diverge near the critical point
[40]. In that limit the finite lifetime of the collision prevents the system
having equilibrium. The assumption here that the short-range correlations can
be assigned the strength of the equilibrium susceptibility would then need to
be adjusted. In principle, the current approach could be modified to account
for this by treating that strength as a dynamic quantity. Non-equilibrium
considerations of the local correlations and the accompanying hydrodynamic
response are intriguing, but they will not be considered in this study.

###  III.2 Green’s Functions

Figure 3:  Green’s functions, denoted as
$G_{YX}(\tau_{\text{source}}=1\,\text{fm}/c,\tau=11\,\text{fm}/c,\Delta\eta)$,
represent the response to a perturbation in $X$, initiated at $1\,\text{fm}/c$
at $\eta=0$, and observed in $Y$ at $\Delta\eta$. Panels $a,b,c,d$ display
$G_{YX}$ for $Y=E,P_{\eta},P_{x}$ and $B$ respectively. For each $Y$, each
panel shows the non-zero Green’s functions for the quantities described in the
legends. These functions are computed for the soft version of the equation of
state, and further details can be found in the accompanying text. Cross-
correlations between different conserved charges are considered, including
longitudinal energy-momentum, baryon, and electric charges. To enhance
visibility, some of the functions were scaled by large factors.

In order to obtain a solution, a set of Green’s functions for point sources in
$(\tau,\eta)$ must be generated. Taking advantage of the boost invariance of
the Bjorken flow, we consider only $(\tau,0)$ and derive a solution for a
source at any $\eta$ using Lorentz transformations. Green‘s functions were
computed by assigning initial condition at $\tau_{\rm start}$ in some charge
$X$:

$\displaystyle G_{XX}(\tau=\tau_{\rm
start},\eta)=\frac{1}{\sqrt{2\pi\sigma^{2}_{\eta}}\tau}\exp\left\\{-\frac{\eta^{2}}{2\sigma^{2}_{\eta}}\right\\},$
(48)

as $\sigma_{\eta}\to 0$, $G_{XX}(\tau=\tau_{\rm
start},\eta)\to{\delta(\eta)\over\tau}$ approaches solution for point source
(Green‘s function by definition). Using very small but finite sigma
regularizes our system by regularizing high momentum modes. Figure (3) depict
the response to a Gaussian source with $\sigma^{2}_{\eta}=0.01$ and $\tau_{\rm
start}=1\text{fm}/c$ at $\tau=11\text{fm}/c$. It is noticeable that a
perturbation in one charge leads to a response in another. This can be
comprehended from the expressions, Eq. (32-35). Equation (34), in its second
term, describes the convection of background charge by a perturbation of
velocity. This effectively couples the response in charges to the response in
energy-momentum. Equations (32) and (33) are influenced by the response to
charge fluctuations through the dependence on the equation of state
$P(\varepsilon,\rho_{B},\rho_{Q},\rho_{S})$.

One expects two forms of response: two sound wave modes linear in wave number
and one quadratic mode that corresponds to dissipation. Sound waves propagate
at the adiabatic speed of sound of the medium and correspond to the long-range
part of correlations. Parts that are due to dissipative processes in the
medium are relatively localized and extend over approximately 1-2 units of
space-time rapidity. This is on the order of the thermal spread of momentum
rapidity for matter at freeze-out. Therefore, it is interesting to observe to
what extent the signal of hydrodynamic evolution will survive particlization.

###  III.3 Correlations

Three equations of state are compared here: the ideal Hadron Resonance Gas
(iHRG) and interacting hadron fluids, with one being harder and the other
softer than the iHRG. Fluctuations involving energy, longitudinal momentum,
baryon density, electric charge, and strange number densities are separated
from those involving transverse momentum, as correlations between them vanish.
The longitudinal hydrodynamic correlations are presented in Figure (4), and
they display noticeable sensitivity to the chosen equation of state.

Figure 4:  The matrix of two-point correlations for energy, longitudinal
momentum, and baryon, electric, and strange charges is presented as a function
of space-time rapidity distance at the freeze-out time $\tau=11\text{fm}/c$.
The considered equations of state include the ideal hadron resonance gas (in
blue), hard hadronic fluid (in orange), and soft hadronic fluid (in green).
Transverse fluctuations of the system are decoupled from other observables and
can be found in a separate plot. $N_{ch}$ is charged particle multiplicity per
unit of $\eta$ as predicted from statistical hadronization at $T=150\,MeV$
with transverse area $S=100\,fm^{2}$. Figure 5: The matrix of two-point
correlations for energy, longitudinal momentum, and baryon, electric, and
strange charges is presented as a function of momentum rapidity at the freeze-
out time $\tau=11\text{fm}/c$, computed after sampling all charged particles.
The equations of state considered include the ideal hadron resonance gas (in
blue), interacting hadronic fluid (in orange), and a fluid with a phase
transition (in green). It’s important to note that transverse fluctuations of
the system do not affect the longitudinal observables, and momentum
fluctuations completely cancel out. $N_{ch}$ is charged particle multiplicity
per unit of $\eta$ as predicted from statistical hadronization at $T=150\,MeV$
with transverse area $S=100\,fm^{2}$.

Electric and strange charges (s),(y) are only weakly influenced by the baryon
mean-field used. However, some effects are observable around
$\Delta\eta\approx 0$. In this case, a harder equation of state produces a
stronger peak, whereas a softer equation of state suppresses the peak in a
small region. These differences disappear at the tails of the correlations.
This effect is more pronounced in the strange charge compared to the electric
charge. This can be attributed to the transition from a medium that includes a
large number of strange and charged mesons at high temperatures to a medium
dominated by cold baryons at freeze-out. Baryons, being more sensitive to the
equation of state, contribute to the central part of the correlation function,
while the tails are created much earlier in the collision. The cross-
correlation between electric charge (Q) and strange charge (S) exhibits a
similar behavior (t) and (s), but here the peak structure changes from a
single maximum for an ideal gas to two peaks for the Hard-EoS or a much
stronger peak for the Soft-EoS.

The baryon charge (m) exhibits more pronounced features. In particular, the
value at the peak, $\Delta\eta=0$, changes from negative to positive when
transitioning from the iHRG to the Soft-EoS to the Hard-EoS. Interactions also
amplify both the maximal and minimal values of the correlation function. This
can be understood from Figure (2), where the rates of change of
$\tau\chi(\tau)$ are much higher for the equation of state with interactions
and sometimes have a different sign compared to the iHRG. The spread of the
correlation here is around one unit in space-time rapidity. Correlations of
the baryon charge with electric or strangeness experience similar effects from
interactions (r),(w),(n),(o).

The longitudinal momentum correlation $\langle\delta P_{\eta}\delta
P_{\eta}\rangle$ (g) for the iHRG is much narrower compared to the interacting
version. It is evident that, in this case, the correlation spreads further for
a hard equation of state (higher speed of sound) and less for the soft
equation of state (smaller speed of sound). The momentum correlation with
other charges appears to be proportional to the rate at which susceptibility
changes.

Energy correlations (a) have a width of about half a unit of rapidity. The
value at the peak changes sign for the hard and ideal hadronic medium, for the
$\langle\delta E\delta E\rangle$ correlation. It’s worth noting that the
transport coefficients corresponding to viscosity for different equations of
state are almost identical. It seems that the sound wave part of energy gets
suppressed and transferred into momentum. It is interesting to see how this
transformation manifests in the correlations at particlization.

Balancing correlation of the transverse momentum, Figure (6) (a), is weakly
influenced by the equation of state. It is noteworthy that in this case, the
correlations spread quite far, covering nearly 5 units of rapidity. The Hard-
EoS has the largest width, while the Soft-EoS has the smallest.

###  III.4 Freeze-out

Figure 6:  Balancing part for the correlations of transverse momentum
$P_{x,y}$ as function of space time (left) or momentum (right) rapidity
differences. As predicted from hydrodynamic evolution with thermal sources
before (left) and after particle creation (right). $N_{ch}$ is charged
particle multiplicity per unit of $\eta$ as predicted from statistical
hadronization at $T=150\,MeV$ with transverse area $S=100\,fm^{2}$.

The final form of the balancing part, computed for all charged particles at
freeze-out, is depicted in Figure (5) and Figure (6) for $E$, $P_{\eta}$, $B$,
$Q$, $S$, $P_{x}$, and $P_{y}$, respectively.

Before particlization, the $P_{x}$ and $P_{y}$ correlations exhibit a strong
positive peak at mid-rapidity with negative dips around it, Figure (6) (a).
However, after particlization, the value at the peak is suppressed, while the
dips are enhanced, Figure (6) (b). The differences between the equations of
state become more pronounced, especially in the case of the Hard-EoS, where
the width is much larger than in the other two cases.

In $(\eta,\phi)$ coordinates, correlations in transverse momentum have
undergone experimental investigation, serving as a sensitive probe for
viscosity. One notable feature that has been observed is the “Ridge” in
rapidity. Correlations in transverse momenta are influenced by various
physical phenomena, in addition to momentum conservation. These include
resonance and cluster decays, radial flow effects, anisotropic flow effects,
initial state fluctuations, and modified jet fragmentation [41, 42]. All of
these effects, with the exception of jets, can be incorporated into the (3D+1)
hydrodynamic modeling of the system. Jets are formed early in the collision
and are not typically part of the collective motion. However, jets interact
with the bulk of the QCD medium, serving as sources of additional correlation.
Therefore, it is important to understand the extent to which jets contribute
to the $p_{t}$ correlation. As evidenced by the comparison of our results with
data, the hydrodynamic history and balancing correlation should account for a
significant portion of the total correlation. The investigation of bulk-jet
interaction is indeed an intriguing and important topic.

While the other correlations are not as wide as the transverse momentum
correlation, the width of the balancing parts increases to 2 units of
rapidity. This is primarily attributed to the thermal spread of particles in
momentum rapidity, denoted as $y=\eta+\delta y$, with $\delta y$ having values
concentrated around the range of $[-0.5,0.5]$. This smearing effect blurs the
features and results in much smoother correlation.

The balancing parts for strange and electric charges mostly preserve their
features, differing at mid-rapidity but converging towards the tails. They
exhibit a simple unimodal structure where a stronger hydrodynamic response
leads to the generation of more balancing correlation. Since we have only
considered balancing correlation created during the hydrodynamic stage, the
normalization appears to be proportional to the difference in susceptibilities
at $\tau=1$ and $\tau=11\text{fm}/c$, which, as can be inferred from Figure
(2), is larger for the Hard-EoS. This difference is notably more significant
for the baryon charge. The minima surrounding the peak in this case become
less pronounced when compared to the correlation before freeze-out. In the
case of the ideal Hadron Resonance Gas (iHRG), the balancing part seems to
remain almost constant around 0, completely losing any rapidity dependence
that was present before.

When longitudinal momentum becomes involved, the correlation vanishes. This is
a consequence translation invariance in $\eta$. If other than periodic
boundary conditions were used that would not be true anymore. If one were to
consider a system that is not exactly transnational invariant, Green’s
functions would become functions of both rapidities (source and response), and
not just their difference.

For energy, the width is now comparable to that of electric or strange
charges, and it is even larger than the width of the baryon charge.
Additionally, it appears that a hard equation of state becomes much more
spread compared to the other equations of state considered. This phenomenon
reflects the presence of hydrodynamic modes in response to the initial
perturbation.

##  IV Conclusions

Here, the entire $7\times 7$ matrix of correlations or conserved quantities
was considered as a function of the measure of separation along the beam axis.
For correlations in coordinate space, correlations were binned as a function
of relative spatial rapidity, $\Delta\eta$, whereas for the final-state
correlations, the momentum rapidity difference $\Delta y$ was chosen. An
extremely simple model was used here. Boost invariance was assumed, which is a
poor choice at lower beam energies, and transverse flow was ignored, which is
a poor assumption at any energy. The equations of state applied were also not
particularly physical. The correlations are seeded by the evolution of the
susceptibilities and evolve according to viscous hydrodynamics and diffusion.
At a fixed breakup time, correlations were projected onto those of free-
streaming hadrons and their decay products, ignoring any additional
contributions from the hadron phase aside from those from decays. For example,
baryon annihilation is ignored here. Only thermal correlations were
considered, which ignore those from the initial state. Initial state
correlations can be driven by jets or minijets, or any local energy-density
fluctuations due to the stopping mechanism. Though the integrated strength of
such correlations is also constrained by conservation laws, their initial
structure is not well understood and might well be much stronger than the
thermally driven correlations discussed here.

Nonetheless, these calculations were sufficient to illustrate the potential of
such correlations to constrain the equation of state and to answer some basic
questions. Given that overcoming the shortcomings enumerated in the previous
paragraphs requires a substantial commitment from both the theory and
experimental communities, this study should, hopefully, clarify whether such a
significant investment might be worthwhile. Here, we summarize the insights
gained and the lessons learned from this study.

  1. 1.

It is possible to consistently model the entire matrix of correlations. This
requires simultaneous consideration of all seven conserved quantities. Because
boost invariance was assumed, and transverse momentum was ignored, the
calculations here were greatly simplified, but the paradigm presented here is
certainly extendable to more realistic cases. More realistic calculations
would first require careful consideration of the algorithms to handle the
additional dimensionality, and adding transverse flow would then allow the
correlations involving $P_{x}$ and $P_{y}$ to mix with the other five
quantities.

  2. 2.

All the correlations were significantly sensitive to the equation of state and
transport coefficients, and such correlations survived the particlization of
being projected from the hydrodynamic stage to free-streaming, but decaying,
hadrons. Correlations often had non-trivial structures, and the widths in
rapidity were $\lesssim 1$ unit of rapidity. Correlations from the initial
state, which were ignored here, should have structures characterized by
$\gtrsim 1$ unit of rapidity, because of the greater time they have to spread.

The strength of thermal correlations for central collisions was found to be
significant, similar in magnitude to those seen experimentally in
$p_{t}-p_{t}$ correlations [43], or correlations often referred to as the
“Ridge” [44]. Thus, the thermally sourced correlations considered here are not
overwhelmed by those of the initial state. Given their narrower structure,
there is hope they might somewhat separable from those of the initial state,
but indications are that this would require significant modeling and analysis,
both experimentally and theoretically.

  3. 3.

Such correlations are highly intertwined, especially at lower beam energy
where the net baryon density significantly differs from zero. At the highest
energies, charge-charge correlations do not mix with those involving energy
and momentum. For collisions studied at the SPS at CERN, or for the Beam-
Energy-Scan at RHIC, correlations appear between charges, such as baryon
number, and energy. Energy and momentum correlations then mix from
hydrodynamics. Thus, it is crucial to simultaneously consider the entire
matrix of correlations, which would include such quantities as the correlation
between baryon number and transverse energy.

  4. 4.

Calculating correlations given the equation of state, transport coefficients,
and initial state effects already represents a significant challenge.
Performing the inverse calculation, where one constrains or infers the
equation of state from measurements, is even more daunting. Given the current
uncertainties in the initial state, especially for lower beam energies,
pursuing this class of correlations at the current time is probably premature
for lower beam energies. However, once all the modeling elements are in place
and better understood, these correlations have the potential to provide
tremendous insight into the evolution of bulk properties of matter throughout
a heavy-ion collision. These correlations should have especially penetrating
insight into any phase transitions.

###### Acknowledgements.

Supported by the Department of Energy Office of Science through Grant No. DE-
FG02-03ER41259.

##  Appendix A Sum Rules

The correlation functions in this study, along with the Green’s functions are
all constrained by energy, momentum and charge conservation. These constraints
were rather easily expressed in Euclidean space, but in Bjorken coordinates
they are more subtle. Energy and momentum transform into one another, because
Bjorken coordinates involve Lorentz boosts, so the simple sum rules, such as
the correlations integrating to the susceptibilities, need to be restated.
This is the goal of this appendix.

In Bjorken coordinates the proper time $\tau$ and the spatial rapidity, $\eta$
replace the Euclidean coordinates of time $t$ and the longitudinal position
$z$ through the transformation,

$\displaystyle\tau$ $\displaystyle\equiv$ $\displaystyle\sqrt{t^{2}-z^{2}},$
(49) $\displaystyle\eta$ $\displaystyle\equiv$ $\displaystyle\tanh^{-1}(z/t).$

The collective velocity is $v_{z}=\tanh\eta$. For boost invariant systems, all
local thermodynamic quantities depend only on $\tau$, including the criteria
for decoupling. Thus, it makes sense to focus on correlations at fixed $\tau$.
At fixed $\tau$ correlations, $\langle\delta A(\eta_{1})\delta
B(\eta_{2})\rangle$ depend only on the relative rapidity, assuming that $A$ or
$B$ refer to the energy, momentum and charge densities expressed in the local
Bjorken frame. The local Bjorken frame involves boosting with velocity
$\tanh\eta$ from the laboratory frame. In the absence of fluctuations this
frame is the same as the local fluid frame, though the fluid frame will differ
slightly once the perturbations are considered. The goal for this section is
to express the sum rules for the correlations, defined in the local rest
frames, in terms of integrals over relative spatial rapidity, $\Delta\eta$.

###  A.1 Constraints for Correlation Functions

First, we consider correlations involving the energy and the longitudinal
momentum densities. For any operator $A(\eta=0)$,

$\displaystyle\int d\eta\langle\delta\rho_{E}(\eta)A(0)\rangle$
$\displaystyle=$ $\displaystyle 0,$ (50) $\displaystyle\int d\eta\langle\delta
M_{z}(\eta)A(0)\rangle$ $\displaystyle=$ $\displaystyle 0$

due to energy and momentum conservation, and $\delta\rho_{E}$ and $\delta
M_{z}$ are the changes to the energy and momentum per some slice of spatial
rapidity and the energy and momenta are defined in the laboratory frame. One
can write $\delta\rho_{E}$ and $\delta M_{z}$ in terms of $\delta\epsilon$ and
$\delta\tilde{M}_{z}$, which are the energy and momentum densities in the
local Bjorken frame,

$\displaystyle\delta\rho_{E}$ $\displaystyle=$
$\displaystyle\delta\epsilon\cosh\eta+\delta\tilde{M}_{z}\sinh\eta,$ (51)
$\displaystyle\delta M_{z}$ $\displaystyle=$
$\displaystyle\delta\tilde{M}_{z}\cosh\eta+\delta\epsilon\sinh\eta.$

One can define a four-vector, $\delta M^{\mu}=T^{\mu\nu}d\Omega_{\mu}$, where
$d\Omega$ is a hyper-volume element. In the local Bjorken frame, for a cross-
sectional area $A$, $d\Omega_{\mu}=A(\tau d\eta,0,0,0)$, and in the laboratory
frame $d\Omega_{\mu}=A(\tau d\eta\cosh\eta,0,0,\tau d\eta\sinh\eta)$. Given
that the break hyper-surface element, $d\Omega_{\mu}$, is defined by constant
$\tau$, the differential four-momentum in the Bjorken frame is
$\delta\tilde{M}^{\mu}=(\delta\epsilon,0,0,\delta\tilde{M}_{z})$.

Energy and momentum conservation thus give the following constraints:

$\displaystyle\int
d\eta\langle(\delta\epsilon(\eta)\cosh\eta+\delta\tilde{M}_{z}(\eta)\sinh\eta)\delta\epsilon(0)\rangle$
$\displaystyle=$ $\displaystyle 0,$ (52) $\displaystyle\int
d\eta\langle(\delta\epsilon(\eta)\cosh\eta+\delta\tilde{M}_{z}(\eta)\sinh\eta)\delta\tilde{M}_{z}(0)\rangle$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\int
d\eta\langle(\delta\tilde{M}_{z}(\eta)\cosh\eta+\delta\epsilon(\eta)\sinh\eta)\delta\epsilon(0)\rangle$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\int
d\eta\langle(\delta\tilde{M}_{z}\cosh\eta+\delta\epsilon(\eta)\sinh\eta)\delta\tilde{M}_{z}(0)\rangle$
$\displaystyle=$ $\displaystyle 0.$

Next we define correlations, $C_{ab}(\eta)$ in such a way that they ignore
zero range correlations, and assume the zero range correlations (particles
with themselves mostly) are defined by the susceptibilities.

$\displaystyle\langle(\delta\epsilon(\eta)\delta\epsilon(0)\rangle$
$\displaystyle=$ $\displaystyle C_{EE}(\eta)+\chi_{EE}\delta(\eta),$ (53)
$\displaystyle\langle(\delta\epsilon(\eta)\delta\tilde{M}_{z}(0)\rangle$
$\displaystyle=$ $\displaystyle C_{EP}(\eta),$
$\displaystyle\langle(\delta\tilde{M}_{z}(\eta)\delta\epsilon(0)\rangle$
$\displaystyle=$ $\displaystyle C_{PE}(\eta),$
$\displaystyle\langle(\delta\tilde{M}_{z}(\eta)\delta\tilde{M}_{z}(0)\rangle$
$\displaystyle=$ $\displaystyle C_{PP}(\eta)+\chi_{PP}\delta(\eta).$

The susceptibility $C_{EE}$ is the specific heat (multiplied by $T^{2}$) and
$C_{PP}=(P+\epsilon)T$. From inserting these definitions into Eq.(s) (52),

$\displaystyle\int
d\eta\left[C_{EE}(\eta)\cosh(\eta)+C_{PE}(\eta)\sinh(\eta)\right]$
$\displaystyle=$ $\displaystyle-A\tau\chi_{EE},$ (54) $\displaystyle\int
d\eta\left[C_{PP}(\eta)\cosh(\eta)+C_{EP}(\eta)\sinh(\eta)\right]$
$\displaystyle=$ $\displaystyle-A\tau\chi_{PP}.$

The functions $C_{EE}$ and $C_{PP}$ are even functions in $\eta$ and $C_{PE}$
and $C_{EP}$ are odd functions, with $C_{PE}(\eta)=-C_{EP}(\eta)$. These
correlations were defined in terms of operators at positions $\eta$ and zero,
but all the correlations $C_{ab}(\eta_{1},\eta_{2})$ are functions of
$\eta_{1}-\eta_{2}$ only, thus the $\eta$ dependence above can be replaced
with the relative spatial rapidity, $\Delta\eta$.

Subtracting the two expressions in Eq. (54) from one another provides a
particularly useful expression as it ultimately leads to expressions for the
transverse energy,

$\displaystyle\int d\eta\left[C_{EE}(\eta)\cosh\eta+C_{PE}(\eta)\sinh\eta-
C_{EP}(\eta)\sinh\eta-C_{PP}(\eta)\cosh\eta\right]$ $\displaystyle=$
$\displaystyle-\chi_{EE}+\chi_{PP}.$ (55)

At large times, when the particles are free streaming, the rapidity and
spatial rapidity become equal. The right-hand side of Eq. (55), which
represents the local correlation of a particle with itself simply becomes the
density of tracks per rapidity multiplied by the average $E_{t}^{2}$ per
track. The left-hand correlations involving the $\delta\tilde{M}_{z}$ all
vanish because, asymptotically, the particles have no longitudinal momentum as
the rapidity $y$ approaches $\eta$. The only correlation remaining is thus the
one involving the transverse energies, and

$\displaystyle\int dy\leavevmode\nobreak\
\langle\delta\rho_{t}(y)\delta\rho_{t}(0)\rangle\cosh\eta$ $\displaystyle=$
$\displaystyle-\frac{dN}{dy}\langle E_{t}^{2}\rangle.$ (56)

Here $\rho_{t}(y)$ is the transverse energy per rapidity, and the averaging,
$\langle...\rangle$, ignores the correlation of a particle with itself.

One might ask how the other combination, $C_{EE}+C_{PP}$, might constrain the
final state. Summing the two expressions,

$\displaystyle\int d\eta\left[C_{EE}(\eta)+C_{PP}(\eta)\right]\cosh\eta$
$\displaystyle=$ $\displaystyle(-\chi_{EE}-\chi_{PP})A\tau.$ (57)

However, once the particles rapidity approach their spatial rapidities, one
ends up with the same constraint. This is not to say there exists only one
useful constraint. Both constraints provide stringent tests of a correlation
calculation.

Constraints for correlations involving the transverse momentum components, or
the baryon, electric and strangeness charge densities are easier to derive
because the boost from the Euclidean coordinates to Bjorken coordinates does
not affect $P_{x},P_{y}$ or the charge densities. The constraints are then

$\displaystyle\int
d\eta\langle\delta\tilde{M}_{x}(\eta)\delta\tilde{M_{x}}(0)\rangle$
$\displaystyle=$ $\displaystyle\int
d\eta\delta\tilde{M}_{y}\delta\tilde{M_{y}}$ (58) $\displaystyle=$
$\displaystyle-A\tau(P+\epsilon)T^{2},$ (59) $\displaystyle\int
d\eta\langle\delta\rho_{i}(\eta)\delta\rho_{j}(0)\rangle$ $\displaystyle=$
$\displaystyle-A\tau\chi_{ij},$ $\displaystyle\int
d\eta\langle\delta\rho_{i}(\eta)\delta\epsilon(0)$ $\displaystyle=$
$\displaystyle 0.$

Due to odd reflection symmetries,

$\displaystyle\int d\eta\langle\delta M_{z}(\eta)\delta A(0)\rangle$
$\displaystyle=$ $\displaystyle 0$ (60)

for any $A=\epsilon,\rho_{i},M_{x}$ or $M_{y}$. Finally, because transverse
expansion is neglected, there are no cross terms involving $M_{x}$ or $M_{y}$,
and

$\displaystyle\int d\eta\langle\delta M_{x}(\eta)\delta A(0)\rangle=\int
d\eta\langle\delta M_{y}(\eta)\delta A(0)\rangle=0$ (61)

for any $A$ unless $A=M_{x}$ or $M_{y}$ for the two expectations above.

###  A.2 Constraints for Green’s Functions due to Conservation Laws

If one has energy density (energy per spatial rapidity), $\delta\epsilon$, at
time $\tau_{0}$, the resulting energy, momentum and charge conservation can be
expressed in terms of Green’s functions,

$\displaystyle\delta\epsilon(\tau,\eta)$ $\displaystyle=$ $\displaystyle\int
d\eta_{0}G_{EE}(\tau,\tau_{0},\Delta\eta=\eta-\eta_{0})\delta\epsilon(\tau_{0},\eta_{0}),$
(62) $\displaystyle\delta\tilde{M}_{z}(\tau,\eta)$ $\displaystyle=$
$\displaystyle\int
d\eta_{0}G_{PE}(\tau,\tau_{0},\Delta\eta=\eta-\eta_{0})\delta\epsilon(\tau_{0},\eta_{0}),$
$\displaystyle\delta\tilde{\rho}(\tau,\eta)$ $\displaystyle=$
$\displaystyle\int
d\eta_{0}G_{QE}(\tau,\tau_{0},\Delta\eta=\eta-\eta_{0})\delta\epsilon(\tau_{0},\eta_{0}).$

Considering the case where $\delta\epsilon(\tau_{0},\eta_{0})$ is a delta
function at $\eta_{0}=0$,

$\displaystyle\delta\epsilon(\tau_{0},\epsilon)=\delta E_{\rm
tot}\delta(\eta_{0}),$ (63)

the total energy is $E_{\rm tot}$, which along with momentum and charge must
then be conserved. The energy, momentum and charge in the lab frame are
defined as

$\displaystyle E_{\rm tot}$ $\displaystyle=$ $\displaystyle\int
d\eta\left(\delta\epsilon(\eta,\tau)\cosh\eta+\delta\tilde{M}_{z}(\eta,\tau)\sinh\eta\right),$
(64) $\displaystyle P_{z{\rm tot}}$ $\displaystyle=$ $\displaystyle 0=\int
d\eta\left(\delta\tilde{M}_{z}(\eta,\tau)\cosh\eta+\delta\epsilon(\eta,\tau)\sinh\eta\right),$
$\displaystyle Q_{\rm tot}$ $\displaystyle=$ $\displaystyle 0=\int
d\eta\delta\rho_{E}(\eta,\tau).$

The momentum and energy expressions above can be derived by starting with
writing the energy and momentum in terms of the stress-energy tensor,

$\displaystyle\delta M^{\mu}=\int\leavevmode\nobreak\ \delta
T_{\mu\nu}d\Omega_{\nu},$ (65)

where $d\Omega_{\nu}$ is the differential hyper volume, $\sim\tau d\eta$ in
the local Bjorken frame. In that frame $\delta T_{00}=\delta\epsilon/\tau$ and
$T_{0z}=\delta\tilde{M}_{z}/\tau$. Boosting the resulting energy-momentum four
vector to the lab frame adds the factors, $\cosh\eta$ and $\sinh\eta$.

Energy, momentum and charge conservation the reguire the following constraints
on the Green’s functions,

$\displaystyle\int
d\Delta\eta\left[G_{EE}(\tau,\tau_{0},\Delta\eta)\cosh\eta+G_{PE}(\tau,\tau_{0},\Delta\eta)\sinh\eta\right]$
$\displaystyle=$ $\displaystyle 1,$ (66) $\displaystyle\int
d\Delta\eta\left[G_{PE}(\tau,\tau_{0},\Delta\eta)\cosh\eta+G_{EE}(\tau,\tau_{0},\Delta\eta)\sinh\eta\right]$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\int
d\Delta\eta\leavevmode\nobreak\ G_{QE}(\tau,\tau_{0},\Delta\eta)$
$\displaystyle=$ $\displaystyle 0.$

Similarly, one can consider initial fluctuations in the momentum density. The
resulting constraints are:

$\displaystyle\int
d\Delta\eta\left[G_{PP}(\tau,\tau_{0},\Delta\eta)\cosh\eta+G_{EP}(\tau,\tau_{0},\Delta\eta)\sinh\eta\right]$
$\displaystyle=$ $\displaystyle 1,$ (67) $\displaystyle\int
d\Delta\eta\left[G_{EP}(\tau,\tau_{0},\Delta\eta)\cosh\eta+G_{PP}(\tau,\tau_{0},\Delta\eta)\sinh\eta\right]$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\int
d\Delta\eta\leavevmode\nobreak\ G_{QP}(\tau,\tau_{0},\Delta\eta)$
$\displaystyle=$ $\displaystyle 0.$

For the charge,

$\displaystyle\int d\Delta\eta\leavevmode\nobreak\
G_{QQ}(\tau,\tau_{0},\Delta\eta)$ $\displaystyle=$ $\displaystyle 1,$ (68)
$\displaystyle\int
d\Delta\eta\left[G_{EQ}(\tau,\tau_{0},\Delta\eta)\cosh\eta+G_{PQ}(\tau,\tau_{0},\Delta\eta)\sinh\eta\right]$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\int
d\Delta\eta\left[G_{PQ}(\tau,\tau_{0},\Delta\eta)\cosh\eta+G_{EQ}(\tau,\tau_{0},\Delta\eta)\sinh\eta\right]$
$\displaystyle=$ $\displaystyle 0.$

Eq.s (66), (67) and (68) represent the constraints on Green’s functions due to
conservation laws.

Finally, the transverse momentum components are also conserved. Given that, in
the Bjorken limit with no transverse expansion, they do not mix with the other
components, and that they are invariant to longitudinal boosts, the
constraints are simple,

$\displaystyle\int d\Delta\eta\leavevmode\nobreak\
G_{P_{x}P_{x}}(\tau,\tau_{0},\Delta\eta)=\int d\Delta\eta\leavevmode\nobreak\
G_{P_{y}P_{y}}(\tau,\tau_{0},\Delta\eta)$ $\displaystyle=$ $\displaystyle 1.$
(69)

##  Appendix B Incorporating Interactions into the Equation of State and
Fluctuations

The equations of state in this study use the hadron gas as a base, then add
contributions to interactions. For a given equation of state the calculations
require knowledge of all fluctuations. For hydrodynamic calculations the input
variables are the energy density and the three charge densities, i.e.,
thermodynamics in the microcanonical ensemble. However, the pressure, and
especially the fluctuations are most clearly defined in the grand canonical
ensemble. Further, the interaction terms used here are functions of the
density and temperature, i.e. the canonical ensemble.

###  B.1 Fluctuations in the Grand Canonical Ensemble

For a free gas, one can relatively easily program the pressure and
fluctuations in terms of the temperature and chemical potentials. The grand
canonical partition function is

$\displaystyle Z_{GC}$ $\displaystyle=$
$\displaystyle\frac{PV}{T}=\sum_{h}(2S_{h}+1)\frac{V}{(2\pi)^{3}}\int
d^{3}p\leavevmode\nobreak\ e^{-(E_{p}-\mu_{h})/T}$ $\displaystyle=$
$\displaystyle\sum_{h}(2S_{h}+1)m_{h}^{2}Te^{\mu_{h}/T}K_{2}(m_{h}/T),$

where $m_{h},S_{h}$ and $\mu_{h}$ are the mass, spin and chemical potential
for a hadron species $h$. The energy density and fluctuations can then readily
be obtained by taking derivatives of the grand-canonical partition function.
If one wants such quantities in terms of the energy and charge densities, one
can first find the temperature and chemical potentials from the energy and
charge densities through a Newton’s method. Once one knows the chemical
potentials the fluctuations can be found quickly.

Susceptibilities and fluctuations are more conveniently defined in terms of
$T$ and $\mu/T$, which is defined here as $\tilde{\mu}_{a}=\mu_{a}/T$. First,
the charge fluctuations are

$\displaystyle\chi_{ab}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{V}\langle\delta Q_{a}\delta Q_{b}\rangle$
$\displaystyle=$
$\displaystyle\frac{\partial\rho_{a}}{\partial\tilde{\mu}_{b}}.$

Energy-energy and energy-charge fluctuations are

$\displaystyle\chi_{EQ_{a}}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{V}\langle\delta E\delta Q_{a}\rangle$ $\displaystyle=$
$\displaystyle\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}=T^{2}\frac{\partial\rho_{a}}{\partial
T},$ $\displaystyle\chi_{EE}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{V}\langle\delta E\delta E\rangle$ $\displaystyle=$
$\displaystyle T^{2}\frac{\partial\epsilon}{\partial T}.$

###  B.2 Fluctuations in the Canonical Ensemble

Adding interactions complicates the procedure. Interactions are most easily
expressed in terms of density, which suggests using the canonical partition
function. The Helmholtz free energy density, which is simply related to the
partition function is convenient when adding interaction terms as functions of
the charge density $\rho$ and the temperature $T$. However, expressing the
fluctuations requires some care.

In terms of the canonical partition function the Helmholtz free energy density
is

$\displaystyle f(\rho,T)$ $\displaystyle=$
$\displaystyle\frac{1}{V}(E-TS)=-\frac{T}{V}\ln Z_{C}.$ (73)

The differential free energy density is, after using the identity $\delta
V/V=-\delta\rho/\rho$,

$\displaystyle\delta f$ $\displaystyle=$
$\displaystyle\frac{\delta\rho}{\rho}f+\frac{1}{V}(-S\delta T-P\delta V)$
$\displaystyle=$ $\displaystyle\frac{\delta\rho}{\rho}f-s\delta
T+\frac{P\delta\rho}{\rho}$ $\displaystyle=$
$\displaystyle\mu\delta\rho-s\delta T.$

The last relation used the identity,

$\displaystyle s=\beta\epsilon+P/T-\mu\rho/T.$ (75)

Reading off the last expression,

$\displaystyle s$ $\displaystyle=$ $\displaystyle-\frac{\partial f}{\partial
T},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\mu=\frac{\partial f}{\partial\rho},$ (76)

where $f$ is assumed to be a function of $\rho$ and $T$. Using the definition
of $f$,

$\displaystyle\epsilon$ $\displaystyle=$ $\displaystyle f+Ts=f+T\frac{\partial
f}{\partial T}=-T^{2}\frac{\partial}{\partial T}\left(\frac{f}{T}\right).$
(77)

Using Eq. (75) one can write the pressure,

$\displaystyle P$ $\displaystyle=$ $\displaystyle-f+\mu\rho.$ (78)

Summarizing, to this point we have the basic thermodynamic quantities in terms
of $f(T,\rho)$:

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle-T^{2}\frac{\partial}{\partial T}\left(\frac{f}{T}\right),$ (79)
$\displaystyle s$ $\displaystyle=$ $\displaystyle-\frac{\partial f}{\partial
T},$ $\displaystyle\mu$ $\displaystyle=$ $\displaystyle\frac{\partial
f}{\partial\rho},$ $\displaystyle P$ $\displaystyle=$
$\displaystyle-f+\mu\rho.$

####  B.2.1 Charge Susceptibilities in the Canonical Ensemble

To see how to express such quantities in the canonical ensemble, one may
consider the matrix

$\displaystyle M_{ab}$
$\displaystyle\equiv\frac{\partial\tilde{\mu}_{a}}{\partial\rho_{b}}.$ (80)

The matrix is simply the inverse of $\chi_{ab}$,

$\displaystyle\delta\rho_{a}$ $\displaystyle=$
$\displaystyle\chi_{ab}\delta\tilde{\mu}_{b},$ (81)
$\displaystyle\delta\tilde{\mu}_{a}$ $\displaystyle=$ $\displaystyle
M_{ab}\delta\rho_{b}.$

If the lower expression is multiplied by the inverse of $M$,

$\displaystyle\delta\rho_{a}$ $\displaystyle=$ $\displaystyle
M^{-1}_{ab}\delta\tilde{\mu}_{b}.$ (82)

Thus, one can see that

$\displaystyle M$ $\displaystyle=$
$\displaystyle\chi^{-1},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \chi=M^{-1}.$ (83)

Thus, to get the charge fluctuations, which were defined in the grand-
canonical ensemble, from a formulation based on the $f(T,\rho)$, one can
proceed via,

$\displaystyle\tilde{\mu}_{a}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\frac{\partial f}{\partial\rho_{a}},$ (84)
$\displaystyle\chi^{-1}_{ab}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\frac{\partial^{2}f}{\partial\rho_{a}\partial\rho_{b}}.$

The charge fluctuation is then found by inverting $\chi^{-1}$.

####  B.2.2 The Specific Heat and Energy-Charge Fluctuation in the Canonical
Ensemble

One can also find the fluctuation of the energy density, at fixed
$\tilde{\mu}$, in terms of $f(T,\rho)$.

$\displaystyle\delta\epsilon$ $\displaystyle=$
$\displaystyle\frac{\partial\epsilon}{\partial\rho_{a}}\delta\rho_{a}+\left.\frac{\partial\epsilon}{\partial
T}\right|_{{\rm fixed}\leavevmode\nobreak\ \rho}\delta T,$ (85)
$\displaystyle\delta\tilde{\mu}_{a}$ $\displaystyle=$
$\displaystyle\frac{\partial\tilde{\mu_{a}}}{\partial\rho_{b}}\delta\rho_{b}+\left.\frac{\partial\tilde{\mu_{a}}}{\partial
T}\right|_{{\rm fixed}\leavevmode\nobreak\ \rho}\delta T=0,$
$\displaystyle\delta\rho_{a}$ $\displaystyle=$
$\displaystyle-\left(\frac{\partial\tilde{\mu}}{\partial\rho}\right)_{ab}^{-1}\left.\frac{\partial\tilde{\mu}_{b}}{\partial
T}\right|_{{\rm fixed}\leavevmode\nobreak\ \rho}\delta T,$ $\displaystyle=$
$\displaystyle-\chi_{ab}\frac{\partial\tilde{\mu}_{b}}{\partial T}\delta T,$
$\displaystyle\delta\epsilon$ $\displaystyle=$ $\displaystyle\delta
T\left\\{-\frac{\partial\epsilon}{\partial\rho_{a}}\chi_{ab}\left.\frac{\partial\tilde{\mu_{b}}}{\partial
T}\right|_{{\rm fixed}\leavevmode\nobreak\
\rho}+\left.\frac{\partial\epsilon}{\partial T}\right|_{{\rm
fixed}\leavevmode\nobreak\ \rho}\right\\}.$

This yields

$\displaystyle\chi_{EE}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon}{\partial T}\right|_{{\rm
fixed}\leavevmode\nobreak\
\rho}-\frac{\partial\epsilon}{\partial\rho_{a}}\chi_{ab}\left.\frac{\partial\tilde{\mu}_{b}}{\partial
T}\right|_{{\rm fixed}\leavevmode\nobreak\ \rho}.$ (86)

Here, even if one is calculationg a contribution from a subset of the free
energy, e.g. that due to interactions, $\chi_{ab}$ is the total charge
susceptibility.

Next, one can calculate $\chi_{EQ}$. We do similar steps above, except set
$\delta T=0$.

$\displaystyle\delta\epsilon$ $\displaystyle=$
$\displaystyle\frac{\partial\epsilon}{\partial\rho_{a}}\delta\rho_{a},$ (87)
$\displaystyle\delta\rho_{a}$ $\displaystyle=$
$\displaystyle\chi_{ab}\delta\tilde{\mu}_{b},$ $\displaystyle\delta\mu_{a}$
$\displaystyle=$ $\displaystyle\chi^{-1}_{ab}\delta\rho_{b},$
$\displaystyle\chi_{EQ}$ $\displaystyle=$
$\displaystyle\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}=\frac{\partial\epsilon}{\partial\rho_{b}}\chi_{ba}.$

Again, even if one is calculating the contribution to some subset of the free
energy, $\chi_{ab}$ in Eq. (B.2.2)must refer to the total charge
susceptibility.

###  B.3 Adding Contributions from the Non-Interacting Gas and from
Interactions

Because the interaction energy is mainly a function of the density, it is
natural to express the Helmholtz free-energy as two separate contributions,
which then add together. The first contribution to the free energy, that from
the non-interacting gas, can be calculated by beginning with the grand
canonical ensemble. The contribution from interactions, can then be added. For
non-momentum dependent interactions one expects the free-energy density due to
interactions to depend only on the density, not the temperature, though the
derivation here does not assume that. We then express how the charge
susceptibilites can then be generated from the Helmholtz free-energy density.
It should be emphasized that stating that the interaction contribution to
$f(\rho,T)$ adds as separate term, is crucial to this derivation. One could
have stated that the interactions provided a separate term to the pressure,
with both terms evaluated at the same chemical potential. But, that represents
a different physical assumption that seems more difficult to motivate given
that it is more natural to assume that the interaction contributes as a
function of the densities, not as a function of the chemical potentials. For
example, if one added an additional non-interacting species, the chemical
potentials would change, but the densities of the interacting particles would
not. For mean-field interactions, one might expect that the energy, $E$,
changes with density, mostly independent of temperature, and that the entropy
is unaffected. That implies that $F=E-TS$ would be altered by adding some
function of the density. Further, one may wish to understand the
susceptibilities in the mixed-phase region, i.e. densities at which the grand-
canonical ensemble gives undefined fluctuations.

For our aims, we need to have thermodynamic quantities as a function of
density, not chemical potential. Going from one to the other is not
particularly difficult, unless one has an unstable mixed-phase region. If one
wishes to model the dynamics of such a phase it is imperative to work with
quantities which are functions of density. Further, the interaction energy is
more straight-forward to describe as an addition to the free-energy density,
because it is more often more easily understood as a function of density. For
that reason one might think it reasonable to express the free-energy density
as

$\displaystyle f(\rho,T)$ $\displaystyle=$ $\displaystyle f_{\rm
free}(\vec{\rho},T)+f_{\rm int}(\vec{\rho},T).$ (88)

In fact, if the interaction is purely density dependent, i.e. independent of
the momenta of the underlying particles, then $f_{\rm int}(\rho,T)$ is simply
the potential energy density, and is independent of temperature.

Given the definition in Eq. (107), one can immediately see that

$\displaystyle\tilde{\mu}_{\rm tot}(\vec{\rho},T)$ $\displaystyle=$
$\displaystyle\tilde{\mu}_{\rm free}(\vec{\rho},T)+\tilde{\mu}_{\rm
int}(\vec{\rho},T),$ (89) $\displaystyle P_{\rm tot}(\vec{\rho},T)$
$\displaystyle=$ $\displaystyle P_{\rm free}(\vec{\rho},T)+P_{\rm
int}(\vec{\rho},T),$ $\displaystyle\epsilon_{\rm tot}(\vec{\rho},T)$
$\displaystyle=$ $\displaystyle\epsilon_{\rm free}(\vec{\rho},T)+\epsilon_{\rm
int}(\vec{\rho},T).$

Note that this definition of $P_{\rm int}$ differs from one where the total
pressure is defined by separating out the two contributions at the same
chemical potential. The latter choice would be natural if one were working in
the grand canonical ensemble, in terms of $\vec{\tilde{\mu}}$, rather than in
terms of the density.

Thus, give the choice of Eq. (107) one can see that

$\displaystyle\frac{\partial^{2}f_{\rm
tot}}{\partial\rho_{a}\partial\rho_{b}}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}f_{\rm
free}}{\partial\rho_{a}\partial\rho_{b}}+\frac{\partial^{2}f_{\rm
int}}{\partial\rho_{a}\partial\rho_{b}},$ (90) $\displaystyle\chi^{-1}_{\rm
tot}$ $\displaystyle=$ $\displaystyle\chi^{-1}_{\rm free}+\chi^{-1}_{\rm
int}.$

This if one wants the charge susceptibility of the entire system, one should
calculate $\chi^{-1}$ for both pieces, add them together, then take the
inverse of the matrix.

$\displaystyle\chi$ $\displaystyle=$ $\displaystyle(\chi^{-1}_{\rm
free}+\chi^{-1}_{\rm int})^{-1}.$ (91)

The susceptibilities which include energy can also be found in terms of
density-dependent (rather that $\tilde{\mu}-$dependent functions. Beginning by
expressing the grand partition function in terms of $T$ and $\tilde{\mu}$, one
can write the energy charge fluctuation,

$\displaystyle\zeta_{a}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{VT}\langle\delta E\delta
Q_{a}\rangle=\frac{1}{T}\frac{\partial\epsilon}{\partial\tilde{\tilde{\mu}}_{a}}=T\frac{\partial\rho_{a}}{\partial
T}.$ (92)

One can then write this as

$\displaystyle\zeta_{a}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\frac{\partial\epsilon}{\partial\rho_{b}}\frac{\partial\rho_{b}}{\partial\tilde{\mu}_{a}}$
$\displaystyle=$
$\displaystyle\frac{1}{T}\chi_{ab}\frac{\partial\epsilon}{\partial\rho_{b}}.$

Rewriting the specific heat requires remembering that derivatives w.r.t.
temperature require noting whether $\vec{\tilde{\mu}}$ or $\vec{\rho}$ are
being held constant. First, we define the specific heat so that it corresponds
to the energy fluctuation in the grand canonical ensemble,

$\displaystyle\xi$ $\displaystyle\equiv$
$\displaystyle\frac{1}{T^{2}V}\langle\delta E\delta E\rangle$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon}{\partial
T}\right|_{\vec{\tilde{\mu}}}.$

Next, one can determine how $\vec{\tilde{\mu}}$ must be changing in order to
maintain constant density,

$\displaystyle\delta\rho_{a}$ $\displaystyle=$ $\displaystyle
0=\left.\frac{\partial\rho_{a}}{\partial
T}\right|_{\vec{\tilde{\mu}}}+\frac{\partial\rho_{a}}{\partial\tilde{\mu}_{b}}\delta\tilde{\mu}_{b},$
(95) $\displaystyle\delta\tilde{\mu}_{a}$ $\displaystyle=$
$\displaystyle-\chi^{-1}_{ab}\left.\frac{\partial\rho_{a}}{\partial
T}\right|_{\vec{\tilde{\mu}}}.$

This then gives

$\displaystyle\delta\epsilon$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon}{\partial
T}\right|_{\vec{\tilde{\mu}}}+\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}\delta\tilde{\mu}_{a}$
(96) $\displaystyle\left.\frac{\partial\epsilon}{\partial
T}\right|_{\vec{\rho}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon}{\partial
T}\right|_{\vec{\tilde{\mu}}}-\left.\frac{\partial\rho_{a}}{\partial
T}\right|_{\vec{\tilde{\mu}}}\chi^{-1}_{ab}\frac{\partial\epsilon}{\partial\tilde{\mu}_{b}}$
$\displaystyle=$ $\displaystyle\left.\frac{\partial\epsilon}{\partial
T}\right|_{\vec{\tilde{\mu}}}-\zeta_{a}\chi^{-1}_{ab}\zeta_{b}.$

For our calculation, we have two components, that of a non-interacting gas,
plus an interaction. The thermodynamics of the non-interactiong gas are most
easily expressed in terms of chemical potentials, whereas the interaction term
is most easily expressed in terms of densities. To add quantities together, we
assume that the net free energy is the sum of the free energy of a gas and the
interaction part,

$\displaystyle f$ $\displaystyle=$ $\displaystyle f^{\rm(gas)}+f^{\rm(int)}.$
(97)

This then implies that

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\epsilon^{\rm(gas)}+\epsilon^{\rm(int)},$ (98)
$\displaystyle\tilde{\mu}_{a}$ $\displaystyle=$
$\displaystyle\tilde{\mu}_{a}^{\rm(gas)}+\tilde{\mu}_{b}^{\rm(int)},$
$\displaystyle\chi^{-1}$ $\displaystyle=$
$\displaystyle\chi^{-1{\rm(gas)}}+\chi^{-1{\rm(int)}},$
$\displaystyle\left.\frac{\partial\epsilon}{\partial T}\right|_{\vec{\rho}}$
$\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon^{\rm(gas)}}{\partial
T}\right|_{\vec{\rho}}+\left.\frac{\partial\epsilon^{\rm(int)}}{\partial
T}\right|_{\vec{\rho}}.$

Note that the chemical potentials of the two parts add (they are not the same)
and that the densities are the same.

The order of operations for adding the two contributions are:

  1. 1.

Calculate all quantities for the gas in the grand canonical ensemble. This
includes $\epsilon^{\rm(gas)}$, $\zeta_{a}^{\rm(gas)}$, and $\chi^{\rm(gas)}$.
Also take the inverse of $\chi^{\rm(gas)}$ to get $\chi^{-1{\rm(gas)}}$. Also
calculate

$\displaystyle\frac{\partial\epsilon^{\rm(gas)}}{\partial\rho_{a}}$
$\displaystyle=$
$\displaystyle\frac{1}{T}\chi^{-1{\rm(gas)}}_{ab}\zeta^{\rm(gas)}_{b},$ (99)
$\displaystyle\left.\frac{\partial\epsilon^{\rm(gas)}}{\partial
T}\right|_{\vec{\rho}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon^{\rm(gas)}}{\partial
T}\right|_{\vec{\tilde{\mu}}}-\zeta^{\rm(gas)}_{a}\chi^{-1{\rm(gas)}}_{ab}\zeta^{\rm(gas)}_{b},$

  2. 2.

Calculate the quantities for the interaction contribution:

$\displaystyle\epsilon^{\rm(int)}$ $\displaystyle=$
$\displaystyle-T^{2}\frac{\partial}{\partial T}(f^{\rm(int)}/T),$ (100)
$\displaystyle\tilde{\mu}_{a}^{\rm(int)}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\frac{\partial f^{\rm(int)}}{\partial\rho_{a}},$
$\displaystyle\chi^{-1{\rm(int)}}_{ab}$ $\displaystyle=$
$\displaystyle\frac{\partial\tilde{\mu}^{\rm(int)}_{a}}{\partial\rho_{b}}.$

along with $\partial\epsilon^{\rm(gas)}/\partial\rho_{a}$ and
$\partial\epsilon^{\rm(gas)}/\partial T|_{\vec{\rho}}$.

  3. 3.

Add contribution from gas and interactions to get combined quantities,

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\epsilon^{\rm(gas)}+\epsilon^{\rm(int)},$ (101) $\displaystyle
P$ $\displaystyle=$ $\displaystyle P^{\rm(gas)}+P^{\rm(int)},$
$\displaystyle\tilde{\mu}_{a}$ $\displaystyle=$
$\displaystyle\tilde{\mu}_{a}^{\rm(gas)}+\tilde{\mu}_{a}^{\rm(int)},$
$\displaystyle\chi^{-1}_{ab}$ $\displaystyle=$
$\displaystyle\chi^{-1{\rm(gas)}}_{ab}+\chi^{-1{\rm(int)}}_{ab},$
$\displaystyle\frac{\partial\epsilon}{\partial\rho_{a}}$ $\displaystyle=$
$\displaystyle\frac{\partial\epsilon^{\rm(gas)}}{\partial\rho_{a}}+\frac{\partial\epsilon^{\rm(int)}}{\partial\rho_{a}},$
$\displaystyle\left.\frac{\partial\epsilon}{\partial T}\right|_{\vec{\rho}}$
$\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon^{\rm(gas)}}{\partial
T}\right|_{\vec{\rho}}+\left.\frac{\partial\epsilon^{\rm(int)}}{\partial
T}\right|_{\vec{\rho}}.$

Now, from all the combined quantities, one can calculate all the net
susceptibilities,

$\displaystyle\chi_{ab}$ $\displaystyle=$
$\displaystyle\left(\chi^{-1}\right)^{-1}_{ab},$ (102)
$\displaystyle\zeta_{a}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\chi_{ab}\frac{\partial\epsilon}{\partial\rho_{b}},$
$\displaystyle\xi$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon}{\partial
T}\right|_{\vec{\rho}}+\zeta_{a}\chi^{-1}_{ab}\zeta_{b}.$

###  B.4 The speed of sound

When taking partial derivatives, assume functions are of $T$ and
$\tilde{\mu}_{a}\equiv\mu_{a}/T$.

$\displaystyle c_{s}^{2}$ $\displaystyle=$ $\displaystyle\left.\frac{\partial
P}{\partial\epsilon}\right|_{s/\rho},$ (103) $\displaystyle\delta P$
$\displaystyle=$ $\displaystyle\frac{\partial P}{\partial T}\delta
T+\frac{\partial P}{\partial\tilde{\mu}_{a}}\delta\tilde{\mu}_{a},$
$\displaystyle=$ $\displaystyle\frac{P+\epsilon}{T}\delta
T+\rho_{a}T\delta\tilde{\mu}_{a},$ $\displaystyle\delta\epsilon$
$\displaystyle=$ $\displaystyle\frac{\partial\epsilon}{\partial T}\delta
T+\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}\delta\tilde{\mu}_{a},$
$\displaystyle\delta\rho_{a}$ $\displaystyle=$
$\displaystyle\frac{\partial\rho_{a}}{\partial T}\delta
T+\frac{\partial\rho_{a}}{\partial\tilde{\mu}_{b}}\delta\tilde{\mu}_{b}$
$\displaystyle=$ $\displaystyle\frac{\partial\rho_{a}}{\partial T}\delta
T+\chi_{ab}\delta\tilde{\mu}_{b},$ $\displaystyle\delta s$ $\displaystyle=$
$\displaystyle\frac{\partial s}{\partial T}\delta T+\frac{\partial
s}{\partial\tilde{\mu}_{a}}\delta\tilde{\mu}_{a},$ $\displaystyle
s^{2}\delta\left(\frac{\rho_{a}}{s}\right)$ $\displaystyle=$ $\displaystyle
s\delta\rho_{a}-\rho_{a}\delta s=0,$ $\displaystyle 0$ $\displaystyle=$
$\displaystyle s\frac{\partial\rho_{a}}{\partial T}\delta
T+s\chi_{ab}\delta\tilde{\mu}_{b}-\rho_{b}\frac{\partial s}{\partial T}\delta
T-\rho_{a}\frac{\partial s}{\partial\tilde{\mu}_{b}}\delta\tilde{\mu}_{b},$
$\displaystyle\delta\tilde{\mu}_{a}$ $\displaystyle=$ $\displaystyle
M_{a}\frac{\delta T}{T},$ $\displaystyle M_{a}$ $\displaystyle=$
$\displaystyle T\left(s\chi-\rho\frac{\partial
s}{\partial\tilde{\mu}}\right)_{ab}^{-1}\left(\rho_{b}\frac{\partial
s}{\partial T}-s\frac{\partial\rho_{b}}{\partial T}\right).$

By substituting for $\delta\tilde{\mu}$ in the expressions for $\delta P$ and
$\delta\epsilon$, then taking the ratio $\delta P/\delta\epsilon$, the factors
of $\delta T$ drop out and

$\displaystyle c_{s}^{2}$ $\displaystyle=$
$\displaystyle\frac{(P+\epsilon)/T+\rho_{a}M_{a}}{\frac{\partial\epsilon}{\partial
T}+\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}M_{a}/T}$ (104)

Next, substitute to replace derivative of $s$, with derivatives of $\epsilon$
or $\rho$, so that all derivatives in the Eq. (B.4) are related to
fluctuations.

$\displaystyle\frac{\partial P}{\partial T}$ $\displaystyle=$
$\displaystyle\frac{P+\epsilon}{T},$ (105) $\displaystyle\frac{\partial
s}{\partial T}$ $\displaystyle=$ $\displaystyle\frac{\partial}{\partial
T}\left(\frac{P+\epsilon-\tilde{\mu}_{a}\rho_{a}T}{T}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{T}\left(\frac{P+\epsilon}{T}+\frac{\partial\epsilon}{\partial
T}-T\tilde{\mu_{a}}\frac{\partial\rho_{a}}{\partial
T}-\tilde{\mu}_{a}\rho_{a}\right)-\frac{s}{T}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\frac{\partial\epsilon}{\partial
T}-\tilde{\mu}_{a}\frac{\partial\rho_{a}}{\partial T},$
$\displaystyle\frac{\partial s}{\partial\tilde{\mu}_{a}}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\frac{\partial}{\partial\tilde{\mu}_{a}}(P+\epsilon-\tilde{\mu}_{b}\rho_{b}T),$
$\displaystyle=$
$\displaystyle\frac{1}{T}\left(\rho_{a}T+\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}-\rho_{a}T-\tilde{\mu_{b}}T\frac{\partial\tilde{\rho}_{b}}{\partial\tilde{\mu}_{a}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{T}\left(\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}-\tilde{\mu}_{b}T\frac{\partial\rho_{b}}{\partial\tilde{\mu}_{a}}\right)$
$\displaystyle M_{a}$ $\displaystyle=$ $\displaystyle
A^{-1}_{ab}\left(\rho_{b}\frac{\partial\epsilon}{\partial
T}-\rho_{b}T\tilde{\mu}_{c}\frac{\partial\rho_{c}}{\partial
T}-Ts\frac{\partial\rho_{b}}{\partial T}\right),$ $\displaystyle A_{ab}$
$\displaystyle=$ $\displaystyle
s\chi_{ab}-\rho_{a}\left(\frac{1}{T}\frac{\partial\epsilon}{\partial\tilde{\mu}_{b}}-\tilde{\mu}_{c}\frac{\partial\rho_{c}}{\partial\tilde{\mu}_{b}}\right)$
$\displaystyle=$ $\displaystyle
s\chi_{ab}-\frac{\rho_{a}}{T}\frac{\partial\epsilon}{\partial\tilde{\mu}_{b}}+\rho_{a}\chi_{bc}\tilde{\mu}_{c}.$

Now, let’s summarize the results for the speed of sound using the following
definitions for the susceptibilities,

$\displaystyle\xi$ $\displaystyle\equiv\frac{\partial\epsilon}{\partial
T}=\frac{1}{T^{2}V}\langle\delta E\delta E\rangle,$ (106)
$\displaystyle\zeta_{a}$
$\displaystyle\equiv\frac{1}{T}\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}=T\frac{\partial\rho_{a}}{\partial
T}=\frac{1}{TV}\langle\delta E\delta Q_{a}\rangle,$ $\displaystyle\chi_{ab}$
$\displaystyle\equiv\frac{\partial\rho_{a}}{\partial\tilde{\mu}_{b}}=\frac{1}{V}\langle\delta
Q_{a}\delta Q_{b}\rangle,$ $\displaystyle F_{p}$ $\displaystyle=$
$\displaystyle\frac{P+\epsilon}{T}=\frac{1}{T^{2}v}\langle\delta P_{x}\delta
P_{x}\rangle,$

where $P_{x}$ is a $x-$component of the momentum.

This gives

$\displaystyle c_{s}^{2}$ $\displaystyle=$
$\displaystyle\frac{(P+\epsilon)/T+\rho_{a}M_{a}}{\xi+\zeta_{a}M_{a}},$ (107)
$\displaystyle M_{a}$ $\displaystyle=$ $\displaystyle
A^{-1}_{ab}[(\xi-\tilde{\mu}_{c}\zeta_{c})\rho_{b}-s\zeta_{b}],$
$\displaystyle A_{ab}$ $\displaystyle=$ $\displaystyle
s\chi_{ab}-\rho_{a}\zeta_{b}.$

###  B.5 Derivatives at fixed energy density

Stepping away from the topic of fluctuations, hydrodynamic codes often require
taking derivatives of various thermodynamic quantities with respect to density
at fixed temperature. Here, we summarize the relations that were required
here. First, we review derivative of the pressure.

$\displaystyle\left.\frac{\partial P}{\partial\rho_{a}}\right|_{T}$
$\displaystyle=$ $\displaystyle\frac{\partial
P}{\partial\tilde{\mu}_{b}}\frac{\partial\tilde{\mu}_{b}}{\partial\rho_{a}}$
$\displaystyle=$ $\displaystyle\chi^{-1}_{ab}\rho_{b}T,$
$\displaystyle\left.\frac{\partial\epsilon}{\partial\rho_{a}}\right|_{T}$
$\displaystyle=$
$\displaystyle\chi^{-1}_{ab}\frac{\partial\epsilon}{\partial\tilde{\mu}_{b}}$
$\displaystyle=$ $\displaystyle\chi^{-1}_{ab}\zeta_{b}T=0,$
$\displaystyle\left.\delta P\right|_{\rho}$ $\displaystyle=$
$\displaystyle\left[\left.\frac{\partial P}{\partial
T}\right|_{\tilde{\mu}}-\left.\frac{\partial
P}{\partial\tilde{\mu}_{a}}\right|_{\tilde{\mu}}\frac{1}{T}\chi^{-1}_{ab}\zeta_{b}\right]\delta
T,$ $\displaystyle\left.\frac{\partial P}{\partial T}\right|_{\rho}$
$\displaystyle=$
$\displaystyle\frac{P+\epsilon}{T}-\rho_{a}\chi^{-1}_{ab}\zeta_{b}.$

Similarly, derivatives of $\epsilon$ are

$\displaystyle\left.\frac{\partial\epsilon}{\partial T}\right|_{\rho}$
$\displaystyle=$ $\displaystyle\left.\frac{\partial\epsilon}{\partial
T}\right|_{\tilde{\mu}}-\frac{\partial\epsilon}{\partial\tilde{\mu}_{a}}\frac{1}{T}\chi^{-1}_{ab}\zeta_{b}$
$\displaystyle=$ $\displaystyle\xi-\zeta_{a}\chi^{-1}_{ab}\zeta_{b}.$

Next, we need to find derivatives of $P$ and $\epsilon$ w.r.t $T$ at fixed
density,

$\displaystyle\delta P$ $\displaystyle=$ $\displaystyle\frac{\partial
P}{\partial T}\delta T+\frac{\partial
P}{\partial\tilde{\mu}_{a}}\delta\tilde{\mu}_{a},$ (110)
$\displaystyle\delta\rho_{a}$ $\displaystyle=$
$\displaystyle\frac{\partial\rho}{\partial T}\delta
T+\frac{\partial\rho_{a}}{\partial\tilde{\mu}_{b}}\delta\tilde{\mu}_{b}=0$
$\displaystyle=$ $\displaystyle\frac{1}{T}\zeta_{a}\delta
T+\chi_{ab}\delta\tilde{\mu}_{b}$

Now, to find $\partial P/\partial\rho_{a}|_{\epsilon}$,

$\displaystyle\delta P$ $\displaystyle=$ $\displaystyle\left.\frac{\partial
P}{\partial\rho_{a}}\right|_{T}\delta\rho_{a}+\left.\frac{\partial P}{\partial
T}\right|_{\rho}\delta T,$ (111) $\displaystyle\delta\epsilon$
$\displaystyle=$
$\displaystyle\left.\frac{\partial\epsilon}{\partial\rho_{a}}\right|_{T}\delta\rho_{a}+\left.\frac{\partial\epsilon}{\partial
T}\right|_{\rho}\delta T=0$ $\displaystyle\delta P$ $\displaystyle=$
$\displaystyle\left.\frac{\partial
P}{\partial\rho_{a}}\right|_{T}\delta\rho_{a}-\frac{\delta\rho_{a}}{\partial\epsilon/\partial
T|_{\rho}}\left.\frac{\partial P}{\partial
T}\right|_{\rho}\left.\frac{\partial\epsilon}{\partial\rho_{a}}\right|_{T},$
$\displaystyle\left.\frac{\partial P}{\partial\rho_{a}}\right|_{\epsilon}$
$\displaystyle=$ $\displaystyle\left.\frac{\partial
P}{\partial\rho_{a}}\right|_{T}-\frac{1}{\xi}\left.\frac{\partial P}{\partial
T}\right|_{\rho}\left.\frac{\partial\epsilon}{\partial\rho_{a}}\right|_{T}.$

Finally, to find $\partial P/\partial\epsilon$ at fixed density, one can apply
Eq.s (B.5) and (B.5),

$\displaystyle\left.\frac{\partial P}{\partial\epsilon}\right|_{\rho}$
$\displaystyle=$ $\displaystyle\frac{\left.\frac{\partial P}{\partial
T}\right|_{\rho}}{\left.\frac{\partial\epsilon}{\partial T}\right|_{\rho}}$
(112)

###  B.6 Algorithmic Design

The code used in this study has simple routines that calculate the free-
particle contributions to $\vec{\rho},P,\epsilon$ and $\chi$ in terms of
$\vec{\mu}$ and $T$. An optimized Newton’s method routine is used to find
$\vec{\mu}$ in terms of $\vec{\rho}$ and $T$. This results in an easily
accessible routine that gives all the thermodyanmic quantities in terms of
$\vec{\rho}$ and $T$. One can invert $\chi$ to find the free-energy
contribution to $\chi^{-1}$. This thus provides $P_{\rm free}$, $\epsilon_{\rm
free},\vec{\mu}_{\rm free}$ and $\chi^{-1}_{\rm free}$.

Next, one calls a function that provides the free energy, plus all its
derivatives with respect to $\vec{\rho}$. This provides the interaction
contributions, $f_{\rm int},P_{\rm int}$, $\epsilon_{\rm int},\vec{\mu}_{\rm
int}$ and $\chi^{-1}_{\rm int}$. If one has an analytic form for $f_{\rm
int}(\vec{\rho},T)$ it is usually simple to find all the interaction
contributions by taking derivatives of $f_{\rm int}$ with respect to
temperature and densities. Finally, one finds the thermodynamic quantities of
the entire system by simply adding the two pieces according to the steps
outlined earlier.

## References

  * Busza  _et al._ [2018] W. Busza, K. Rajagopal, and W. van der Schee, Heavy Ion Collisions: The Big Picture, and the Big Questions, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Bzdak  _et al._ [2020] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan, Phys. Rept. 853, 1 (2020), arXiv:1906.00936 [nucl-th] .
  * Dietrich  _et al._ [2020] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Multimessenger constraints on the neutron-star equation of state and the hubble constant, Science 370, 1450 (2020).
  * Beloin  _et al._ [2019] S. Beloin, S. Han, A. W. Steiner, and K. Odbadrakh, Simultaneous fitting of neutron star structure and cooling data, Phys. Rev. C 100, 055801 (2019).
  * Most  _et al._ [2020] E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, and H. Stoecker, Eur. Phys. J. A 56, 59 (2020).
  * Weih  _et al._ [2020] L. R. Weih, M. Hanauske, and L. Rezzolla, Phys. Rev. Lett. 124, 171103 (2020).
  * Most  _et al._ [2018] E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018).
  * Legred  _et al._ [2021] I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Impact of the psr $\mathrm{J}0740+6620$ radius constraint on the properties of high-density matter, Phys. Rev. D 104, 063003 (2021).
  * Drischler  _et al._ [2021] C. Drischler, S. Han, J. M. Lattimer, M. Prakash, S. Reddy, and T. Zhao, Limiting masses and radii of neutron stars and their implications, Phys. Rev. C 103, 045805 (2021).
  * Huth  _et al._ [2022] S. Huth, P. T. H. Pang, I. Tews, T. Dietrich, A. L. Fèvre, A. Schwenk, W. Trautmann, K. Agarwal, M. Bulla, M. W. Coughlin, and C. V. D. Broeck, Constraining neutron-star matter with microscopic and macroscopic collisions, Nature 606, 276 (2022).
  * [11] D. Neill, R. Preston, W. G. Newton, and D. Tsang, arXiv:2208.00994 [astro-ph.HE] .
  * Sorensen  _et al._ [2023a] A. Sorensen _et al._ , Dense nuclear matter equation of state from heavy-ion collisions (2023a), arXiv:2301.13253 [nucl-th] .
  * Hebeler  _et al._ [2010] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Constraints on neutron star radii based on chiral effective field theory interactions, Phys. Rev. Lett. 105, 161102 (2010).
  * Gandolfi  _et al._ [2012] S. Gandolfi, J. Carlson, and S. Reddy, Maximum mass and radius of neutron stars, and the nuclear symmetry energy, Phys. Rev. C 85, 032801 (2012).
  * Drischler  _et al._ [2014] C. Drischler, V. Somà, and A. Schwenk, Microscopic calculations and energy expansions for neutron-rich matter, Phys. Rev. C 89, 025806 (2014).
  * Steinheimer and Randrup [2012] J. Steinheimer and J. Randrup, Spinodal amplification of density fluctuations in fluid-dynamical simulations of relativistic nuclear collisions, Phys. Rev. Lett. 109, 212301 (2012), arXiv:1209.2462 [nucl-th] .
  * Steinheimer and Randrup [2013] J. Steinheimer and J. Randrup, Spinodal density enhancements in simulations of relativistic nuclear collisions, Phys. Rev. C 87, 054903 (2013), arXiv:1302.2956 [nucl-th] .
  * Savchuk  _et al._ [2023] O. Savchuk, R. V. Poberezhnyuk, A. Motornenko, J. Steinheimer, M. I. Gorenstein, and V. Vovchenko, Phase transition amplification of proton number fluctuations in nuclear collisions from a transport model approach, Phys. Rev. C 107, 024913 (2023), arXiv:2211.13200 [hep-ph] .
  * Vovchenko  _et al._ [2020] V. Vovchenko, O. Savchuk, R. V. Poberezhnyuk, M. I. Gorenstein, and V. Koch, Connecting fluctuation measurements in heavy-ion collisions with the grand-canonical susceptibilities, Phys. Lett. B 811, 135868 (2020), arXiv:2003.13905 [hep-ph] .
  * Sorensen  _et al._ [2023b] A. Sorensen, D. Oliinychenko, L. McLerran, and V. Koch, Measuring the Speed of Sound Using Cumulants of Baryon Number, Acta Phys. Polon. Supp. 16, 1 (2023b), arXiv:2209.04957 [nucl-th] .
  * Ling  _et al._ [2014] B. Ling, T. Springer, and M. Stephanov, Hydrodynamics of charge fluctuations and balance functions, Phys. Rev. C 89, 064901 (2014), arXiv:1310.6036 [nucl-th] .
  * Bass  _et al._ [2000] S. A. Bass, P. Danielewicz, and S. Pratt, Clocking hadronization in relativistic heavy ion collisions with balance functions, Phys. Rev. Lett. 85, 2689 (2000), arXiv:nucl-th/0005044 .
  * Pratt  _et al._ [2015] S. Pratt, W. P. McCormack, and C. Ratti, Production of Charge in Heavy Ion Collisions, Phys. Rev. C 92, 064905 (2015), arXiv:1508.07031 [nucl-th] .
  * Pratt and Young [2017] S. Pratt and C. Young, Relating Measurable Correlations in Heavy Ion Collisions to Bulk Properties of Equilibrated QCD Matter, Phys. Rev. C 95, 054901 (2017), arXiv:1612.08983 [nucl-th] .
  * Pratt and Plumberg [2019] S. Pratt and C. Plumberg, Evolving Charge Correlations in a Hybrid Model with both Hydrodynamics and Hadronic Boltzmann Descriptions, Phys. Rev. C 99, 044916 (2019), arXiv:1812.05649 [nucl-th] .
  * Pratt and Plumberg [2021] S. Pratt and C. Plumberg, Charge balance functions for heavy-ion collisions at energies available at the CERN Large Hadron Collider, Phys. Rev. C 104, 014906 (2021), arXiv:2104.00628 [nucl-th] .
  * Aarts  _et al._ [2015] G. Aarts, C. Allton, A. Amato, P. Giudice, S. Hands, and J.-I. Skullerud, Electrical conductivity and charge diffusion in thermal QCD from the lattice, JHEP 02, 186, arXiv:1412.6411 [hep-lat] .
  * Bjorken [1983] J. D. Bjorken, Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region, Phys. Rev. D 27, 140 (1983).
  * Romatschke [2010] P. Romatschke, New Developments in Relativistic Viscous Hydrodynamics, Int. J. Mod. Phys. E 19, 1 (2010), arXiv:0902.3663 [hep-ph] .
  * Pratt  _et al._ [2018] S. Pratt, J. Kim, and C. Plumberg, Evolution of Charge Fluctuations and Correlations in the Hydrodynamic Stage of Heavy Ion Collisions, Phys. Rev. C 98, 014904 (2018), arXiv:1712.09298 [nucl-th] .
  * Pradeep  _et al._ [2022] M. Pradeep, K. Rajagopal, M. Stephanov, and Y. Yin, Freezing out fluctuations in Hydro+ near the QCD critical point, Phys. Rev. D 106, 036017 (2022), arXiv:2204.00639 [hep-ph] .
  * Pratt and Torrieri [2010] S. Pratt and G. Torrieri, Coupling Relativistic Viscous Hydrodynamics to Boltzmann Descriptions, Phys. Rev. C 82, 044901 (2010), arXiv:1003.0413 [nucl-th] .
  * Fotakis  _et al._ [2020] J. A. Fotakis, M. Greif, C. Greiner, G. S. Denicol, and H. Niemi, Diffusion processes involving multiple conserved charges: A study from kinetic theory and implications to the fluid-dynamical modeling of heavy ion collisions, Phys. Rev. D 101, 076007 (2020), arXiv:1912.09103 [hep-ph] .
  * Pratt  _et al._ [2011] S. Pratt, S. Schlichting, and S. Gavin, Effects of Momentum Conservation and Flow on Angular Correlations at RHIC, Phys. Rev. C 84, 024909 (2011), arXiv:1011.6053 [nucl-th] .
  * Kuznietsov  _et al._ [2022] V. A. Kuznietsov, O. Savchuk, M. I. Gorenstein, V. Koch, and V. Vovchenko, Critical point particle number fluctuations from molecular dynamics, Phys. Rev. C 105, 044903 (2022), arXiv:2201.08486 [hep-ph] .
  * Poberezhnyuk  _et al._ [2020] R. V. Poberezhnyuk, O. Savchuk, M. I. Gorenstein, V. Vovchenko, K. Taradiy, V. V. Begun, L. Satarov, J. Steinheimer, and H. Stoecker, Critical point fluctuations: Finite size and global charge conservation effects, Phys. Rev. C 102, 024908 (2020).
  * Kuznietsov  _et al._ [2023] V. A. Kuznietsov, O. Savchuk, R. V. Poberezhnyuk, V. Vovchenko, M. I. Gorenstein, and H. Stoecker, Molecular dynamics analysis of particle number fluctuations in the mixed phase of a first-order phase transition, Phys. Rev. C 107, 055206 (2023), arXiv:2303.09193 [hep-ph] .
  * Sorensen and Koch [2021] A. Sorensen and V. Koch, Phase transitions and critical behavior in hadronic transport with a relativistic density functional equation of state, Phys. Rev. C 104, 034904 (2021), arXiv:2011.06635 [nucl-th] .
  * Pratt [2017] S. Pratt, Consistent implementation of non-zero-range terms into hydrodynamics, Phys. Rev. C 96, 044903 (2017).
  * Mukherjee  _et al._ [2015] S. Mukherjee, R. Venugopalan, and Y. Yin, Real-time evolution of non-gaussian cumulants in the QCD critical regime, Physical Review C 92, 10.1103/physrevc.92.034912 (2015).
  * Bozek and Broniowski [2013] P. Bozek and W. Broniowski, Charge balancing and the fall off of the ridge, Nucl. Phys. A 904-905, 479c (2013), arXiv:1210.4315 [nucl-th] .
  * Bozek and Broniowski [2012] P. Bozek and W. Broniowski, Charge conservation and the shape of the ridge of two-particle correlations in relativistic heavy-ion collisions, Phys. Rev. Lett. 109, 062301 (2012), arXiv:1204.3580 [nucl-th] .
  * Adam  _et al._ [2019] J. Adam _et al._ (STAR Collaboration), Collision-energy dependence of ${p}_{t}$ correlations in au + au collisions at energies available at the bnl relativistic heavy ion collider, Phys. Rev. C 99, 044918 (2019).
  * Abelev  _et al._ [2009] B. I. Abelev _et al._ (STAR), Long range rapidity correlations and jet production in high energy nuclear collisions, Phys. Rev. C 80, 064912 (2009), arXiv:0909.0191 [nucl-ex] .