Freezing Out Fluctuations in Hydro+ Near the QCD Critical Point

We introduce a freeze-out procedure to convert the critical fluctuations in a droplet of quark-gluon plasma (QGP) that has, as it expanded and cooled, passed close to a posited critical point on the phase diagram into cumulants of hadron multiplicities that can subsequently be measured. The procedure connects the out-of-equilibrium critical fluctuations described in concert with the hydrodynamic evolution of the droplet of QGP by extended hydrodynamics, known as Hydro+, with the subsequent kinetic description in terms of observable hadrons. We introduce a critical scalar isoscalar field sigma whose fluctuations cause correlations between observed hadrons due to the couplings of the sigma field to the hadrons via their masses. We match the QGP fluctuations obtained by solving the Hydro+ equations describing the evolution of critical fluctuations before freeze-out to the correlations of the sigma field. In turn, these are imprinted onto correlations and fluctuations in the multiplicity of hadrons, most importantly protons, after freeze-out via the generalization of the familiar half-century-old Cooper-Frye freeze-out prescription which we introduce. The proposed framework allows us to study the effects of critical slowing down and the consequent deviation of the observable predictions from equilibrium expectations quantitatively. We also quantify the suppression of cumulants due to conservation of baryon number. We demonstrate the procedure in practice by freezing out a Hydro+ simulation in an azimuthally symmetric and boost invariant background that includes radial flow discussed in arXiv:1908.08539.

Understanding the physics of strongly interacting matter in extreme conditions and mapping the phase diagram of QCD has been a major goal of theoretical and experimental efforts from high-energy heavy-ion collisions to neutron star mergers [2][3][4]. One possible central feature of the phase diagram of quantum chromodynamics (QCD) -a QCD critical point at the baryon chemical potential µ B above which the crossover via which cooling quark-gluon plasma (QGP) becomes ordinary hadronic matter becomes a first order phase transition -still remains a theoretical conjecture. The challenge of its discovery is being taken up by the current Beam Energy Scan (BES) program at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory as well as at planned facilities worldwide. The intriguing hints observed in the first phase of the BES [2,5], in particular the deviations of certain measures of fluctuations from their non-critical baseline, deviations that vary nonmonotonically as a function of √ s, motivate the ongoing experimental efforts in its second, higher statistics, phase (BES II). BES II data has been taken by the STAR collaboration over the course of 2019-2021 in AuAu collisions at a sequence of collision energies, which corresponds to a scan in µ B [6,7]. The data is currently being analyzed and we look forward with considerable anticipation to learning much from these measurements.
On the theory side, there have been many studies of the observable consequences of critical fluctuations in heavy ion collisions that produce a droplet of QGP that cool close to a critical point upon making the (greatly simplifying, but unrealistic) assumption that these fluctuations stay in equilibrium [4,[6][7][8][9][10][11][12][13][14]. A part of the essence of critical fluctuations is that since their correlation length, ξ, grows near the critical point, the typical timescale for their evolution grows also -this is referred to as critical slowing down. This means that in the rapidly cooling droplets of QGP produced in heavy ion collisions critical fluctuations cannot be described by their equilibrium values slaved to hydrodynamic fields. Fortunately, much progress has also been achieved in describing the non-equilibrium evolution of hydrodynamic fluctuations [1, , whose distinctive behavior at long wavelengths is governed by the universality of critical behavior and thus can serve as a signature of the critical point.
The challenge which still needs to be addressed is establishing a connection between the hydrodynamic fluid -including its critical fluctuations -and the observed particle yieldsand their fluctuations. This work is aimed at closing this gap.
Traditionally one thinks of hydrodynamics as a deterministic theory of fluid evolution.
The hydrodynamic variables are the local fluid velocity u(x) as well as average densities of conserved quantities like the energy density ε and the baryon number density n, or, equivalently, the corresponding conjugate variables such as the temperature T and chemical potentials (like the chemical potential for baryon number, µ B ) characterizing the local equilibrium conditions. These variables evolve deterministically according to hydrodynamic equations. Heavy ion collision experiments, of course, do not measure these hydrodynamic or thermodynamic quantities directly. The conversion to experimentally observable hadron multiplicities occurs at freezeout -a point in the evolution of the expanding cooling droplet of matter where the density becomes low enough that the kinetic description in terms of a hadron gas becomes applicable and the scattering rates for processes that modify the particle species (i.e., chemical) composition is negligible. At that point one can convert the hydrodynamic variables into particle yields and momentum distributions, i.e. spectra. The well-known procedure known as Cooper-Frye freezeout [41] maps the local fluid velocity u(x) and hydrodynamic fields such as T (x) and µ B (x) on the freezeout surface (a hypersurface in spacetime) to a simplified hadronic description in terms of kinetic variables of an expanding ideal resonance gas of hadrons, namely a gas of particles with momenta distributed according to boosted Fermi-Dirac or Bose-Einstein distributions. The interactions are encoded in resonances, which later decay, modifying the distributions ultimately measured by experiment. The average densities of conserved quantities such as energy or baryon number are guaranteed to match as long as the hadron resonance gas provides a good description of the equation of state, including the relations between T and µ B and ε and n in that regime.
The Cooper-Frye freezeout procedure [41] has been successfully employed in the description of high energy heavy-ion collision data for more than four decades to describe average particle yields and spectra. The procedure ensures that the event-averaged baryon number and energy-momentum densities are matched between the hydrodynamic and kinetic theory descriptions. At sufficiently high collision energies √ s, the data from many experiments are in reasonable agreement with this description across a broad kinematic regime. The Cooper-Frye framework, however, does not describe fluctuations in either the hydrodynamic fluid or the kinetic theory particles. Our goal is to extend the Cooper-Frye procedure to the description of critical fluctuations. Such a description is crucial in the special case of heavy ion collisions that freeze out in the vicinity of a critical point. In this case, fluctuations are both enhanced and of considerable interest, since it is via detecting critical fluctuations that we hope to discern the presence of a critical point [6,7]. These fluctuations are due to thermal noise and their magnitude, or more importantly, their correlations are a sensitive signature of the proximity of a thermodynamic singularity, such as the critical point. Obviously, we cannot expect to match these critical fluctuations using a free gas of hadrons.
The effect of the critical correlations can be captured by a critical scalar field, which we call σ -a collective mode which becomes "soft", long-range correlated and slow, at the critical point, justifying its treatment as a collective field. One can then match the singular part of fluctuations of hydrodynamic variables by the fluctuation of the field σ which, via its coupling to the observed particles, causes their masses to fluctuate at the time of freezeout and consequently yields observable fluctuations in particle multiplicities.
In this paper we show how to implement such a procedure and demonstrate its application in a model of the hydrodynamic evolution near the critical point already studied in Ref. [1].
Our paper is organized as follows. In the remainder of this Introduction, in Section I A we review some foundational aspects of critical fluctuations in equilibrium and of critical slowing down. We also introduce the Hydro+ equations that describe the dynamics of out-of-equilibrium fluctuations near a critical point. In Section I B we review the standard Cooper-Frye freezeout procedure that neglects fluctuations. With this groundwork in place, in Section II, we derive and explain our freezeout procedure that extends the Cooper-Frye approach so as to match the critical fluctuations just before freezeout, as described by Hydro+, to observable fluctuations in particle multiplicities just after. In Section III we apply our freezeout procedure to the Bjorken scenario: a fluid that is undergoing boostinvariant longitudinal expansion, meaning that it is cooling, but that is translation-invariant and at rest in the transverse directions. In this simplified setting, we are able to push much of the calculation through analytically and in so doing gain intuition and elucidate general features that arise again in the next Section. In Section IV we illustrate the use of the freezeout procedure that we have introduced and fully exercise its salient features by obtaining the two-point correlations of particle multiplicities from the more realistic Hydro+ simulation of Ref. [1] in which the fluid is boost invariant and azimuthally symmetric but is finite in transverse extent and thus exhibits radial flow. We conclude in Section V with a summary of the main qualitative lessons that we draw from our results of Sections III and IV as well as a look ahead at important next steps.  [42][43][44][45], with a single scalar field becoming soft and slow at the critical point. 1 In QCD, the critical field is a linear combination of scalar operators such as the chiral condensate qq and the baryon number density. This means, in particular, that, in equilibrium near the critical point, the k th cumulants of the baryon number density, which are related to the the derivatives of pressure P with respect to baryon chemical potential µ B at fixed temperature T , diverge as certain powers of ξ [9]: where η ≈ 0.04 is the well-known Ising critical exponent [44,45]. The smallness of this exponent makes it negligible in practical applications we will be concerned with here. In this paper we shall only be concerned with the variance of particle multiplicities, which is to say we shall only need Eq. (1) with k = 2.
The singular behavior of fluctuations predicted by Eq. (1) points the way toward finding signatures of the presence of a critical point in the phase diagram of QCD in heavy-ion collision experiments. As the location of the freeze-out point on the QCD phase diagram, which moves in response to experimentally varying a parameter such as center of mass collision energy √ s, approaches and then passes the location of the QCD critical point, the magnitude and higher cumulants of fluctuations should show a characteristic non-monotonic dependence. There are, however, two essential and quite nontrivial steps which must be taken in order to connect the elegant scaling equation (1) to experimental observables.
First, due to the fact that the droplet of matter produced in a heavy-ion collision expands and cools rapidly it is far from being a static thermodynamic system, meaning that nonequilibrium effects on fluctuations must be considered [15,16] non-equilibrium effects is provided by an extension of hydrodynamics known as Hydro+ [18] that we shall review briefly in this Section.
The second necessary advance arises due to the fact that experimental measurements do not access hydrodynamic variables, such as the densities of conserved quantities arising in Eq.
(1), directly. Therefore, a connection needs to be made between such quantities, and their fluctuations, and the experimentally observable multiplicities of protons, pions, etc. In this paper we propose how this can be done via introducing a scalar field σ in the kinetic theory description of the observable hadrons whose fluctuations can be determined via matching to the critical fluctuations of hydrodynamic quantities, and illustrate it using simplified, but longest times scale [18]. Thus the non-equilibrium dynamics is most important in this mode.
Furthermore, near the critical point, we can imagine mapping the QCD energy density ε and baryon number density n to the 3D Ising entropy density S and Ising magnetization density M or, equivalently, mapping the QCD phase diagram variables T and µ B to the 3D Ising variables r (reduced Ising model temperature) and h (magnetic field) as set up explicitly in Refs. [48,49], and then determining which combination of ε and n corresponds to the most singular Ising model fluctuations. In practice, though, any combination of S and M has the same leading singular behavior in powers of the correlation length ξ as long as it involves fluctuations of M , and a similar statement applies to any combination of ε and n as long as it involves fluctuations of the entropy per baryonŝ ≡ s/n [18,20]. It therefore suffices to use the fluctuations ofŝ, which also happens to be the slowest hydrodynamic mode. We shall match these fluctuations to the leading singular contribution to the fluctuations of the σ field that we shall introduce in Section II.
In the local rest frame of the fluid, the equal-time correlation function ofŝ can be ex-pressed in terms of its Wigner transform: Here x = (x + + x − )/2 and ∆x = x + − x − = (0, ∆x) in the local rest frame of the fluid at the point x. The relaxation of this quantity to its local equilibrium valueφ Q is governed by the equation [18]: For the purposes of this paper, the equilibriumφ Q can be adequately approximated by the Ornstein-Zernike ansatz 2φ where Q ≡ |Q| and c p and n are the heat capacity at constant pressure and the baryon number density.
The Q-dependent relaxation rate, Γ, controls how slowly φ Q (x) evolves toward its equilibrium valueφ Q via Eq. (3). The leading critical behavior of Γ depends on the dynamic universality class. For the QCD critical point [50], it is the one of Model H (liquid-gas critical point) in Halperin and Hohenberg's classification [51], where the linear combination of δε and δn given by δŝ = (δε − (w/n)δn)/(T n) is the slow, diffusive, scalar mode nonlinearly coupled to diffusive (transverse) momentum modes. At the same level of approximation as in Eq. (4), the leading critical behavior of the relaxation rate in model H is given by [52]: where ξ 0 is a typical value of the correlation length well away from the critical point, D 0 is a constant with the dimensions of length that we shall take as a free parameter, and As we shall demonstrate, the most important property of the critical mode relaxation rate given by Eqs. (5)(6) is that it vanishes as Q → 0: 2 While the value ofφ Q at Q = 0 in Eq. (3) is determined by thermodynamics, the dependence on Q in this expression is an often used approximation which takes into account the nonzero correlation length.
A more sophisticated form can be found in Ref. [18].
This reflects the fact that φ Q measures the fluctuation of hydrodynamic variables, which are conserved. The relaxation rate of the 2-point correlator of fluctuations is twice the rate of the relaxation of the corresponding mode, whose relaxation is also diffusive with a diffusion coefficient given by which vanishes at the critical point, where ξ → ∞. We can think of the parameter D 0 which we introduced as the diffusion constant at some reference point well away from the critical point. A crude estimate for D 0 could be obtained by using [51,52] where η and s are the shear viscosity and entropy density, respectively, and where we have taken η ≈ s/(4π). Taking s(T ) =sT 3 withs ≈ 6 as is reasonable around T = T c [53,54] and choosing T c = 160 MeV and ξ 0 = 0.5 fm as we shall throughout, we estimate a critical contribution of D 0 ≈ 0.3 fm. Assuming that the non-critical contribution to D 0 is not too large, we expect D 0 > 0.3 fm but not D 0 0.3 fm. To bracket the uncertainty in this estimate, we shall typically illustrate our results by plotting the results obtained from calculations employing D 0 = 0.25 fm and D 0 = 1 fm.
To elucidate and emphasize the importance of the conservation laws in the dynamics of fluctuations and, consequently, in the experimental signatures of the critical point we shall compare and contrast results obtained using the model H universality class with those which one would have obtained using model A universality class. In the model A universality class, the critical order parameter is not a conserved quantity and the relaxation rate of the fluctuations does not vanish as Q → 0. To the same level of approximation as we have been using so far we can utilize the following ansatz for the relaxation rate in model A: where Γ 0 is a constant with the dimensions of rate (1/time) which we can think of as the relaxation rate at a point well away from the critical point where the correlation length is As the fireball expands, its temperature T drops and the point characterizing the thermodynamic state of the system on the QCD phase diagram moves past the critical point.
The correlation length ξ reaches the maximum value ξ max (see Fig. 1) which depends on how The fluctuation evolution equation (3) depends on the correlation length ξ via the dependence ofφ Q and Γ(Q) on ξ. In a realistic hydrodynamic simulation, ξ will be determined upon solving the hydrodynamic equations with a given equation of state. Since our purpose in this paper is to describe how to freeze out critical fluctuations in hydrodynamics and translate them into observables based on particle multiplicity fluctuations, we shall instead, for simplicity, choose a plausible parametrization of ξ along the expansion trajectory in terms of T .
Specifically, we shall adopt the parametrization of the correlation length along the trajectory of the expanding fireball on the QCD phase diagram in terms of temperature used previously in Ref. [1]: with ∆T = T c /5. We shall not attempt to refine this parametrization in this work. Alternate parametrizations for the correlation length are discussed, e.g., in Refs. [11,30]. The ansatz in Eq. (11) reflects the main features of ξ(T ) relevant for this work -the correlation length reaches a maximum value ξ max at a certain temperature T c (close to the crossover temperature) and then decreases as the system continues to cool on its way to freezeoutas shown in Fig are well described by a procedure which we shall summarize below known as Cooper-Frye freezeout [41] which starts from the output of a hydrodynamic simulation.
In the traditional Cooper-Frye procedure, the macroscopic evolution of the conserved charges and fluid velocity field obtained from a hydrodynamic calculation are converted into a microscopic description in terms of particles in a Hadron Resonance Gas model. The freeze-out hypersurface where this switching is done is determined based on some thermodynamic condition for eg., fixed temperature or energy density. The averages of the conserved densities are equated to those of a hadron resonance gas of particles via the Cooper-Frye formula. Let dS µ be the differential element pointing along the normal vector to the freeze-out surface. The mean multiplicity of particle species A ( N A ) according to the Cooper-Frye formula is given by, Here, d A is the degeneracy of particle species A and Dp A is the Lorentz invariant measure: f A is the momentum-dependent particle distribution function which is either taken to be Fermi-Dirac or Bose-Einstein based on the spin and statistics of hadron species A. For simplicity, throughout this paper we shall ignore the spin and statistics and consider f A to be the Boltzmann distribution where T (x), µ A (x) and u(x) are the temperature, the chemical potential of species A, and the local fluid velocity at a point x on the freeze-out hypersurface. For mesons, µ A = 0, while for baryons/antibaryons µ A = ±µ B , respectively. In addition to ignoring the modification of the distribution function due to spin and statistics, we also ignore further modifications to f A due to viscous effects [55][56][57][58][59][60] in this preliminary study.
In Section II we shall turn to describing our extension of the Cooper-Frye procedure that will enable us to translate the output of a Hydro+ simulation, with traditional hydrodynamic variables as well as fluctuations described by φ Q (x), into particles in a way that faithfully turns the critical fluctuations in the fluid into fluctuations and correlations of the hadrons.
We are pursuing this goal within what is often referred to as a deterministic framework for describing the fluctuations: Hydro+ adds new deterministic equations of motion to the equations of hydrodynamics, equations that describe the evolution of quantities that characterize the fluctuations starting with φ Q (x) that describes their two-point correlation function. Fluctuations can also be described stochastically, where one evolves an ensemble of configurations each with its own realization of the fluctuations [32][33][34][35][36][37][38][39][40]. It would be natural in a stochastic description to analyze freezeout via extending the Cooper-Frye procedure in a manner that follows an analogous logic to that we shall employ here, but we leave this to future work.

II. COOPER-FRYE FREEZE-OUT FOR CRITICAL FLUCTUATIONS
The Cooper-Frye freezeout procedure described in the previous section converted hydrodynamic variables into (event) mean multiplicities. In order to describe the signatures of the QCD critical point we need to be able to describe also the fluctuations and correlations of these multiplicities. We shall now describe a freezeout procedure which extends Cooper-Frye freezeout by connecting the fluctuations of hydrodynamic variables to the fluctuations of the particle multiplicities. We shall focus on the fluctuations of the slowest, and thus most out of equilibrium mode -ŝ. In Hydro+ the two point function of this mode is given by (its Fourier transform) φ Q . Our goal is to connect it to the two-point correlation function of the multiplicity fluctuation δf where f A (x, p) is given by Eq. (14). We shall use the model of critical correlations which incorporates critical fluctuations in the hadronic description via the interaction of the hadrons with a critical σ field. Such a description of critical fluctuations in a hadron gas has been used in equilibrium [7,9,[11][12][13][14][61][62][63][64] as well as with some out-of-equilibrium effects included [65]. In this approach the interaction with the σ field modifies the masses of the hadrons, to linear order in σ, as follows: where σ = 0. The proportionality constant g A plays the role of the coupling constant between the hadron species A and the σ field. The critical contribution to the fluctuations of f A is due to the dependence of the averaged particle distribution function f A on the mass, and is given by where f A is the Boltzmann distribution in Eq. (14). As a result, fluctuations of the σ field translate into fluctuations and correlations between particles, as in the effective action (or free energy) which can be written for small fluctuations at long wavelengths in an expansion in powers of the field σ and its gradients around its equilibrium value, σ = 0, as follows: The equilibrium two-point correlator can be then found from Eq. (19) and is given by where ∆x = x + − x − and ξ ≡ 1/m σ is the correlation length of the σ-field fluctuations. As we shall only be interested in the two-point correlator in this work and as we are neglecting the (small) nonzero value of the critical exponent η, we will be able to neglect the terms of order σ 3 and higher in the expansion (19). The Fourier/Wigner transform of the two-point correlator is then given by In the approximate equalities in Eqs. (20) and (21) we ignore loop corrections, which are known to be small in the 3D Ising universality class in which η is small.
We choose the units of length in Eq. (19) so that the value of ξ introduced in Eqs. (19) and (20) Equivalently, the Wigner transforms are related via the same proportionality constant Z: Using Eqs. (4) and (21) we find Note that, while both c p and ξ 2 diverge at the critical point, their ratio is finite in the approximation we are using. 3 We shall apply the relationship in Eq. preserves the information about important non-equilibrium effects, including the effects of conservation laws. 3 Our approximation sets the critical exponent to its mean-field value η = 0, which is a good approximation to make for a critical point in the 3D Ising universality class where η ∼ 0.04. If one uses a more sophisticated, non-mean-field equation of state as in, e.g., Ref. [48], and/or more sophisticated form of φ Q and χ Q as in Ref. [18], the value of the normalization constant will nevertheless be determined by the matching equation (23), which is more general than the approximation in which we have derived it.
We shall thus determine the correlation functions of σ at freezeout as follows: where Z is a normalization constant which can be obtained by matching the fluctuations obtained in the kinetic description to fluctuations (i.e., susceptibilities) obtained from the QCD equation of state using Eq. (23). Since in this paper our focus is entirely on developing and exploring the implementation of the freezeout prescription that we introduce to describe fluctuations, we shall take the constants Z in Eq. (25b) and g A in Eqs. (16)-(18) as given and postpone their determination by matching a particular QCD EoS to future work. We also note that we shall find ways to express our results that are independent of those unknown parameters. Note that there is a subtlety in defining Eq. (25b) relating to the choice of frame in which x + and x − are at equal time; we shall discuss this in Section II A.
Due to Eq. (25a), the mean number of particles is unmodified by critical fluctuations and is given by Eq. (12). Integrating the spatial correlations given by Eq. (18) over the full freeze-out hypersurface and using Eq. (25b), we can express the leading critical contribution to the correlator of particle multiplicities N A and N B as: with d A the spin (and/or isospin) degeneracy of the particle species A. The subscript σ in δN A δN B σ is there to remind us that this is the contribution due to critical fluctuations. The expressions (26) and (27) constitute the central result whose consequences we shall explore over the course of the rest of this paper by making them explicit in simplified settings. Note that the field σ has now done its job and has now disappeared; in (26) and (27) we have a relationship between the critical fluctuations of hydrodynamic variable on the right-hand side of (26) and the correlator of observable particle multiplicities on the left-hand side.
Although we shall not go beyond two-point correlators in our explicit explorations to come, we note with future work in mind that a straightforward generalization of Eq. (26) yields the form for the critical contribution to the k th cumulant of the multiplicity of particle species A. We have extended the Cooper-Frye procedure in a way that will allow us to see how the critical, i.e. most singular, contribution the two-point correlations ofŝ translates into the variance of observed particle multiplicities. We leave continuing onward to higher-point correlations and non-Gaussian cumulants for future work.
Finally, we note that the total variance of the particle multiplicity has an additional non-critical contribution which is usually taken as Poissonian: There can certainly be corrections to the non-critical contributions that we represent here by the Poisson distribution. These may arise from global charge conservation [66,67] or initial fluctuations [68], for example, but we shall not discuss them in this work. This work is intended only as a prescription for freezing out the fluctuations near the critical point that encode information about the leading singularity.

A. Toward explicit evaluation
The extended Cooper-Frye procedure that we have derived above involves expressions with a certain formality to them. Noting that sometimes the devil turns out to be found in the details, we shall now begin to take the steps need to turn these expressions into tools that can be used in explicit calculations.
The Wigner transform φ Q , as defined in Eq. [1] that we shall analyze in Section IV. In order to translate φ Q into the correlator δŝ(x + )δŝ(x − ) in such a case one needs to be able to evolve this correlator not only in time u·x (using Eq. (3)) but also in time difference u · ∆x. We shall show below that because this evolution is slow If the typical range of the correlator is of order * , then the typical value of the time difference u(x) · ∆x ∼ β * . The typical scale * can be determined by the condition that the relaxation rate Γ(Q * ) ∼ DQ 2 * of the corresponding modes Q * = 1/ * is of order the expansion rate 1/τ , where D was introduced in Eq. (8).
The evolution of the correlator δŝ(x + )δŝ(x − ) as a function of the time separation u · ∆x occurs with the same rate, also of order Γ(Q * ). As a result, the correction to the correlator is of order Γ(Q * ) u · ∆x ∼ β D/τ . This quantity is small already because τ is a macroscopic scale, while D is microscopic, i.e., τ D. Furthermore, near the critical point, D itself is vanishing: as seen in Eq. (8), it is smaller than the microscopic scale by another factor of More formally, let us define a projection of the separation four-vector ∆x onto the plane perpendicular to u(x): Then we can express the correlator δŝ( where the three-component vector ∆x ⊥ is the projection of the four-vector ∆x onto the equal-time hyperplane orthogonal to the vector u(x) as defined in Eq. (30) and illustrated in Fig. 3 Eq. (31) formalizes the qualitative argument from the preceding paragraph.
We shall be comparing our results with what one would obtain upon assuming the fluctuations are in equilibrium. Up to corrections suppressed by ratios of microscopic (e.g., correlation length ξ) to macroscopic scales (such as hydrodynamic gradient scales, e.g., τ ) we can replace the correlation function in Eq. (20) with a delta function: Substituting into Eq. (26) we find wheren(x) is the unit vector along the normal on the freeze-out hypersurface at x. This expression straightforwardly generalizes existing results for equilibrium fluctuations, see for example Ref. [11], to locally equilibrated fluctuations in a (more realistic) inhomogeneous fireball. We shall make comparisons between our full results and the equilibrium fluctuation predictions (34) in order to highlight the importance of non-equilibrium effects, especially the effects due to conservation laws.

B. Ratios of observables
We shall calculate the contribution of critical fluctuations to the variance of the particle multiplicity of species A (we shall consider protons, A = p, and pions, A = π) in a specified finite rapidity and transverse momentum acceptance window. To eliminate the dependence on the volume (i.e., the transverse size) of the droplet of QGP we shall introduce the intensive ratio which was referred to as ω A,σ in Ref. [11]. We note that ω A depends on the choice of acceptance window. (See, e.g., Ref. [13].) Since this acceptance dependence is not the main focus of the present study, while criticality and non-equilibrium effects are, we shall often illustrate our results by plotting the ratiõ where ω nc A is the ω A calculated upon assuming freezeout well away from the critical point, i.e., upon setting ξ max = ξ 0 . We have checked (for a few sets of parameters) that the acceptance dependence of the numerator and denominator in Eq. (36) is similar and, thus, largely cancels. In contrast, the numerator ω A is strongly enhanced by critical fluctuations (for example, in equilibriumω eqbm A = ξ 2 /ξ 2 0 ) and is sensitive to the non-equilibrium effects of critical slowing down, while the denominator ω nc A is, by construction, not affected by critical fluctuations. Although ω A defined in Eq. (36) depends on the unknown parameters g A and Z via the ratio g A / √ Z, all dependence on these unknowns cancels in the ratio of ratios defined in Eq. (36), within the approximations that we shall make. This is a second significant benefit, and we suggest employingω A (and its generalizations to higher cumulants) in the comparison between future theoretical calculations and experimental measurements.
In subsequent Sections, we shall computeω A in two different model hydrodynamic backgrounds. In Section III, we study an analytically solvable scenario with longitudinal boost invariance and no dependence on transverse coordinates, which is to say Hydro+ in Bjorken flow. Then, in Section IV we shall freeze out the numerical simulation of Hydro+ with azimuthal symmetry, radial expansion, and longitudinal boost invariance from Ref.
[1]. Before we conclude this Section, though, we need to establish some notation with which to describe two points on the freeze-out hypersurface and the separation between them.

C. Establishing notation for azimuthally symmetric boost invariant freezeout
To specify the shape of an azimuthally symmetric boost-invariant freezeout surface, it is convenient to use Bjorken coordinates defined in terms of the Cartesian coordinates The mutually orthogonal set of unit vectors corresponding to each of the Bjorken coordinates can be expressed in terms of the Cartesian coordinates aŝ The radial profile of a boost-invariant and azimuthally-symmetric freezeout surface can then be expressed in a parametric form using an arbitrary parameter α as in Ref. [69] τ so that the point on the freeze-out hypersurface corresponding to parameters α, η, ϕ is given by: Then, the volume vector normal to the freeze-out hypersurface can be written as d 3 S = n dατ dηrdϕ where the vector n is given by: The flow four-velocity on the freezeout surface is given by in the coordinates with which we are working.
In defining the two-point correlation function we shall need to specify two points on the freeze-out hypersurface. Let x ± ≡ x(α ± , η ± , ϕ ± ) be any two such points on with x ≡ (x + + x − )/2 being their midpoint and ∆x ≡ x + − x − being the separation vector between them.
Let us denote similarly τ = (τ where ∆x ⊥ was defined in Eq. (30). While the points x + and x − are on the freezeout surface by construction, the midpoint x, in general, is not. The displacement between the midpoint and the freezeout surface is, however, small when the typical range of the correlator is much shorter than the typical curvature radius of the freezeout surface. We can use an argument similar to the one preceding Eq. (30) to simplify the calculation by neglecting the difference between the correlator at the actual midpoint and the correlator at the point on the freezeout surface given by where α ≡ (α + + α − )/2. Henceforth, we shall use x to denote the on-hypersurface approximation (47) to the actual midpoint. Again, neglecting the effect of the curvature and linearizing in ∆η, ∆ϕ and ∆α = α + − α − , the projection of the separation vector ∆x ⊥ from Eq. (46d) onto the hyperplane normal to the four-vector u can be approximated as where u is the 4-velocity of the fluid at the point x, the vector n is defined in Eq. (44), and we have introduced a spacelike unit vector The vectorsα ⊥ ,η andφ form a basis in the hyperplane perpendicular to the four-vector u With all of this notation established, we can now take a step toward making the expression Eq. (26) for the squared variance of the multiplicity of species A that we have derived above as our central result more explicit, writing it as where ∆x ⊥ is a three-vector whose components in theα ⊥ ,η, ϕ basis are given by Eq. (48) and , where n is given by Eq. (44).
The integral in Eq. (27) expressed in Bjorken coordinates takes the form where we used Eq. (14). The Cartesian coordinates in the lab frame of the particle fourmomentum are given by in terms of the particle rapidity y and transverse mass m T ≡ p 2 T + m 2 .

III. FREEZING OUT A HYDRO+ SIMULATION WITH BJORKEN FLOW
In this Section, we apply our approach to the well-known Bjorken scenario: a hot fluid that is undergoing idealized boost-invariant longitudinal expansion, so that it cools as a function of proper time, but that is translation-invariant and at rest in the transverse directions [70].
We shall obtain some results in this simplified scenario in analytic form, thus allowing us to elucidate general features that we shall observe again in a more realistic scenario with transverse expansion in the next Section.  (3) for the fluctuation measure φ Q then takes the form of an ordinary differential equation: where Γ(Q) depends on τ through ξ(τ ) and is specified via Eqs. (5,6,11). Since our focus throughout is on the effects caused by fluctuations near a critical point, for simplicity we shall assume that the initial fluctuations are in equilibrium, i.e., This assumption could certainly be improved in future, but for our present purposes any choice in which the initial fluctuations are small compared to those that develop later will suffice. Since, in the Bjorken scenario, the temperature depends only on τ , the unit fourvector normal to the isothermal freeze-out hypersurface In Fig. 4 we plot φ Q obtained by solving Eq. (53)   The conservation laws keep φ Q stuck at Q = 0, which corresponds to keeping the integral of ∆x 2φ constant. This means that their effect in panel (b) of Fig. 5 is that if there is a (large, critical) correlation at small ∆x they produce a corresponding compensating anticorrelation at longer ∆x. One can also show generally that a peak in φ Q away from Q = 0 corresponds to the anticorrelation (i.e., negative tail) in its Fourier transformφ(∆x ⊥ ). Indeed, if there exist a value of Q at which φ Q > φ 0 , then where we used the fact thatφ(∆x) is an even function. Since | cos(Q · ∆x)| ≤ 1, the inequality (55) cannot be satisfied ifφ(∆x) is always positive.

B. Multiplicity fluctuations and their rapidity correlator
Upon substituting the solution to Eq. (53), or rather its inverse Fourier transformφ, into Eq. (50), we can now calculate the square variance of multiplicity fluctuations δN 2 A . In the simplified setup that we are employing in this Section, we can go one step farther and exploit Bjorken boost invariance to compute explicit results for the rapidity correlator defined as (In the next Section where we employ a more realistic hydrodynamic solution, we shall only compute δN 2 A , not C A .) The correlator C A measures the correlations between the multiplicity of particle species A at rapidities y + and y − and can be determined similarly to Eq. (50) in terms of φ Q or its inverse Fourier transformφ(∆x ⊥ ). For the idealized Bjorken scenario, some of the integrals in Eq. (50) (e.g., over transverse coordinates) can be taken analytically. In order to make even further analytical progress we shall consider the case of particles with mass much bigger than the temperature, m A T . This is an adequate approximation for protons and will allow us to perform an additional integral analytically in that case. We shall not use this approximation in the next Section, where we shall anyway be doing the analogous integrals numerically, but it will be helpful here to make the result and its general features more explicit. As described in detail in Appendix B, upon doing the integrals we obtain which is Eq. (B13), where Q ≡ Q ηη , ∆y ≡ y + − y − , and A ⊥ is the transverse area.
We see that in the simplified setup of this Section in which the fluid is translation invariant in the transverse directions, the modes that contribute in Eq. (57) are those whose Q is directed along theη direction. Also, the effect of the last Gaussian factor in the Q η integral in Eq. (57) is to limit the range of that integral to values of order The fact that the characteristic wavenumber Q of the fluctuations responsible for the correlations at freezeout is not zero (despite considering a volume of fluid that is infinite in extent in rapidity η in this idealized scenario) is ultimately due to the fact that in the laboratory frame the fireball is not spatially homogeneous: the fluid velocity varies over a longitudinal distance of order τ f due to the longitudinal expansion. One can check that taking τ f → ∞ results in only Q = 0 contributing. However, the characteristic Q η is not just 1/τ f , but rather depends on the mass of the particle. This is due to the thermal smearing, or "blurring", which translates spatial Bjorken rapidity η into kinematic particle rapidity y [13,71].
As we are assuming that m A T , the factor T f /(m A cosh η) in Eq. (58) is the typical thermal rapidity of the particles at temperature T f / cosh η, which can be understood as the freezeout temperature "red-shifted" by longitudinal expansion.
The final piece that we need in order to compute C A is the denominator in Eq. (56). By explicit calculation starting from Eq. (12), in the Bjorken scenario in which we are working dN A /dy is given by: Substituting Eqs. (57) and (59) into Eq. (56), we can evaluate C A which, because of boost invariance, is a function of ∆y only. In Fig. 6, we plot our results for C A , normalized by the its non-critical value at ∆y = 0, C nc A (0), which we estimate by substituting φ Q = Z T (Q 2 + ξ −2 0 ) −1 into Eq. (57), where ξ 0 is, as before, the correlation length away from critical point defined in Eq. (11). That is, we define the ratio that we have plotted in Fig. 6 asC Since ξ max > ξ 0 , critical fluctuations make the ratioC A (0) larger than unity.
In Fig. 6 we also illustrate the dependence of the rapidity correlations on the value of the diffusion parameter D 0 . Stronger diffusion (larger D 0 ) enhances the effects of the critical point in φ Q , as we saw in Fig. 5. This enhancement is reflected in the corresponding particle rapidity correlations, as seen in Fig. 6 at small ∆y. Due to conservation laws, anticorrelations at large ∆y are also enhanced. (For any value of D 0 , the consequence of conservation is that the integral under the curveC A over all separations ∆y is independent of D 0 and is determined by the initial fluctuations. This means that when we increase D 0 and seeC A (∆y) growing at small ∆y, it must also become more negative at large ∆y.) However, unlike the direct effect of diffusion on the range of the spatial correlatorφ in Fig. 5b, the effect on the range of C(∆y) in Fig. 6 is minor. This is due to the fact that this range is mostly determined by the effect of thermal rapidity smearing or "blurring" [13]. Now, the variance of the particle multiplicity δN 2 A σ that in the general case takes the form (50) can in this Bjorken scenario be obtained from the rapidity correlator C(∆y) by integration over the rapidity window y ± ∈ (−y max , y max ), i.e., Upon using our expression for C(∆y) from Eqs. (56), (57), and (59) and using the fact that boost invariance implies that N A = 2y max dN A /dy , we now obtain the result The ∆y dependence of C(∆y) translates into the rapidity acceptance window dependence of δN 2 A σ , which has been discussed in the literature, e.g., in Ref. [13], and will not be discussed here. (ii) in Model H dynamics,ω A is less sensitive to the freeze-out temperature than would be the case if the fluctuations were in equilibrium throughout. In the next Section, we shall investigate how radial flow modifies these and other observations.

IV. FREEZING OUT A SEMI-REALISTIC NUMERICAL HYDRO+ SIMULA-TION
In this Section, we demonstrate the use of the freezeout procedure introduced in Section. II and employed in a Bjorken scenario in Section III to obtain the two-point correlations of particle multiplicities from the Hydro+ simulation that was introduced and analyzed, but not frozen out, in Ref.
[1]. As in the previous Section, the system under consideration is boost invariant and has azimuthal symmetry. Unlike in the previous Section, the system we consider here is finite in transverse extent and thus exhibits radial flow. We give a brief description of the details of simulation here. For more details the reader may refer to Ref. [1].
The evolution of the energy density, ε(r, τ ) and the fluid four-velocities u(r, τ ) in our simulation is determined using the standard MIS second order hydrodynamic equations as implemented in the publicly available VH1 + 1 hydrodynamic code [72][73][74]. The equation of state p(ε) used in the simulation was introduced in Ref.
[1] and was already employed in Section III and, for convenience, is reviewed in Appendix A. We set the shear viscosity to entropy density ratio to η/s = 0.08 throughout, and solve the equations numerically using a spatial (radial) lattice with 1024 points spaced by 0.0123 fm and a time step of 0.005 fm. In this Section, we initialize the hydrodynamic simulation at τ i = 1 fm, with an initial central temperature of 330 MeV, following Ref. [72]. We set the radial flow v r and the viscous part of the stress-energy tensor Π µν to zero initially at τ = τ i . We employ the standard Glauber model radial profile corresponding to a central Au-Au collision at √ s = 200 GeV for ε(r) at τ = τ i , again following Ref. [72].
As in Ref.
[1], in our Hydro+ simulation the hydrodynamic densities ε(r, τ ) and Π µν (r, τ ) and the four-velocities u(r, τ ) provide the background for the evolution of the fluctuations described by φ Q according to Eq. (3). Again following Ref.
[1], we choose to initialize the fluctuations φ Q at τ = τ i by setting them to the corresponding equilibrium values determined by the local temperature at this initial time, i.e.
For the interested reader, the limitations of the various assumptions made in setting up this Hydro+ simulation, as well as possible future improvements to it, are detailed in Ref. [1].
We calculate the evolution of φ Q using the same code as in to highlight the effects of conservation laws on the dynamics of φ Q and on the resulting particle multiplicity fluctuations by comparing the results of this Section to those obtained by repeating the calculations of this Section using Model A evolution. We perform this comparison in Appendix C. The second change that we make relative to Ref.
[1] is that here we shall neglect the back-reaction of the fluctuations on the hydrodynamic densities.
The modifications to the bulk dynamics of the hydrodynamic fluid, in particular its entropy density s(r, τ ), introduced by the presence of the fluctuations was computed in Ref. [1,30] and in fact the fluctuations and the hydrodynamic densities were computed self-consistently.
However, these authors showed that including back reaction self-consistently introduces fractional changes to ε(r, τ ) and v r (r, τ ) that are small, rarely comparable to 1% and typically much smaller. For this reason we neglect these effects.
In the remainder of this Section, we demonstrate the implementation of the freezeout prescription introduced in Section II to freeze out the Hydro+ simulation with the back- In Fig. 9, we plot our results for the fluctuation measure φ Q in the hydrodynamic background illustrated in Fig. 8 at three different times τ along two hydrodynamic flow lines, one close to the center of the fireball (r(τ i ) = 0.7 fm) and one further out (r(τ i ) = 5 fm).
We display results from simulations performed with D 0 = 0.25 fm (slower diffusion) and (black curve) φ Q is given by its equilibrium value. In the upper (lower) four panels, the red curves at τ = 9.19 fm (τ = 5.14 fm) are drawn at the time when when the fluid cell moving along the flow line that started at r i = 0.7 fm (r i = 5 fm) has cooled to the temperature T = T c = 160 MeV and the blue curves at τ = 11.36 fm (τ = 6.72 fm) are drawn at the time when these fluid cells have cooled further to T = 140 MeV. Increasing ξ max , i.e., bringing the evolution trajectory closer to the critical point, causes the magnitude of equilibrium fluctuations to increase. However the relaxation to the equilibrium value becomes slower since its rate Γ(Q) ∝ DQ 2 and D = D 0 ξ 0 /ξ is proportional to 1/ξ. We find that the former effect "wins" in the sense that φ Q at the time the fluid cell trajectory cools to T = 140 MeV (i.e. the blue curves in Fig. 9), which is well after the fluid cell trajectory passes the point where ξ = ξ max , see Fig. 1, increases with ξ max , at least in the range of the parameters we have considered.
We see by comparing the left and right columns of Fig. 9 that increasing the diffusion parameter D 0 , which increases the relaxation rate, has the consequence that φ Q is closer to its instantaneous equilibrium formφ Q during the course of the evolution. The value of φ Q at Q = 0, however, remains invariant during the evolution due to the conservation laws inherent in model H: Γ(Q = 0) = 0.
The Q-dependence of φ Q is shaped by two competing effects. As a given hydrodynamic cell, represented by a point on the phase diagram (see Fig. 1) moving along the expansion trajectory, approaches the critical point, the "desired" equilibrium values ofφ Q , to which φ Q is forced to relax by Eq. (3), increases across all values of Q. However, while at larger Q, the relaxation is fast enough to effectively equilibrate φ Q to these larger equilibrium values, at lower Q conservation laws slow down the evolution, making the φ Q values lag behindφ Q more significantly. This produces a peak in φ Q at a characteristic value of Q denoted by Q peak in Ref. [30] which moves to lower values of Q as D 0 is increased. These features are evident in Fig. 9 across the range of parameters we have considered. It is also instructive to compare and contrast Fig. 9 with the results that would be obtained if the fluctuations followed model A dynamics where the relaxation rate of low-Q modes is not suppressed and, consequently, Q peak = 0. We perform this comparison in Appendix C; see Fig. 17 from that Appendix which is to be compared with Fig. 9.
In the simpler Bjorken scenario of Section III, we described "memory effects" and looked at their dependence on Q and the diffusion parameter D 0 . We can do the same here, for the φ Q obtained in this more realistic r-dependent calculation, by displaying our results as in Fig. 10, where suitably normalized φ Q andφ Q are plotted as a function of the local temperature along a fluid cell trajectory for three different values of Q. In accordance with Eq. (53), the value of φ Q increases when φ Q <φ Q and decreases when φ Q >φ Q , as φ Q "tries" to relax toward the rising, and later falling, equilibrium valueφ Q , as the critical point is approached and later passed. For larger values of D 0 (as in the bottom half of Fig. 10) the rate of relaxation is greater, meaning that φ Q rises more rapidly, and therefore higher, in the critical region. Although it also drops more rapidly as the temperature drops further, overall a larger D 0 yields larger fluctuations, at least within the reasonable range of values of D 0 that we explore. For small Q (see the left column in Fig. 10), the value of φ Q grows very slowly, and reaches values much lower that the equilibriumφ Q before it starts decreasing. However, for low Q, the rate at which φ Q decreases after the critical point has been passed is also slow, and as a result significant memory of the fluctuation magnitude near the critical point (albeit itself smaller than equilibrium magnitude) is retained at freezeout. This dynamics is qualitatively similar to the dynamics first described in Ref. [15] in a very simplified model of the out-of-equilibrium evolution of critical fluctuations with no spatial-or Q-dependence.

B. Fluctuations on the freezeout surface
As in the Bjorken scenario discussed in Section III, we consider two isothermal freezeout scenarios with T f = 140 MeV and T f = 156 MeV. The main difference here relative to Section III is that the temperature is now not only a function of τ but also of the radial coordinate r. An isothermal surface, therefore, is not simply τ = const for all r, as in the previous section. The surface T (τ, r) = T f can be parametrized according to the discussion in Section II C and we use the notations and approximations discussed in that section. For simplicity, we choose the parameter α introduced in Eq. (42) according to α = r.
The magnitude of fluctuations at T = 156 MeV in equilibrium is several times higher than that at T f = 140 MeV, since T = 156 MeV is closer to the critical temperature T c = 160 MeV.
However, the time that the system spends in the critical region before freezing out is shorter for T f = 156 MeV than for T f = 140 MeV. By comparing these two freeze-out scenarios, we can understand the sensitivity of out-of-equilibrium fluctuations to the proximity of the freeze-out temperature to the critical point.
In Fig. 11, suitably normalized plots of φ Q are shown for three points on the freeze-  Fig. 11 can be compared with the results obtained in the case of Model A dynamics in Fig. 18 in Appendix C.) At moderate Q the "memory" effect weakens and at large Q the modes closely track their equilibrium values, which rises and then falls as the critical point is approached and then passed. By comparing the plots in Fig. 11 for different T f , we can also see that at smaller, but not too small, Q, the "memory" causes the fluctuation measure φ Q to be larger than its equilibrium valueφ Q . This effect is more pronounced for lower T f , due to the fact that the equilibrium fluctuations are smaller there.
Having understood the effects of varying the parameter D 0 and the proximity to the critical point on the fluctuation measure φ Q , as in Section III the next step toward the calculation of observable particle multiplicity fluctuations is to computeφ(∆x ⊥ ), the inverse Fourier transform of φ Q defined in Eq. (32). In Fig. 12, we plot ∆x 2φ (∆x ⊥ ) as a function of the spatial separation ∆x ≡ |∆x ⊥ | between the two points in the correlator δŝ(x + )δŝ(x − ) , see Eq. (31). By comparing to Fig. 5(b), we see that the D 0 -dependence is qualitatively similar to that in the Bjorken scenario, discussed at length in Section III. The small ∆x (large Q) behavior of the fluctuations is not affected by changing D 0 , while at the same time the spatial correlator becomes longer ranged as D 0 is increased. The correlator goes negative at larger values of ∆x; this is a consequence of conservation as can be seen by comparing Fig. 12 to Fig. 19 in Appendix C and as explained in the context of the Bjorken scenario in Eq. (55). Finally, consistent with what we have already seen in Fig. 11

C. Variance of particle multiplicities
Because of the greater symmetry in the simpler Bjorken scenario that we were employing there, in Section III we were able to compute the rapidity correlator C(∆y) of the multiplicity fluctuations, before integrating over a rapidity window to obtain the variance of the particle multiplicity. Here, with nontrivial radial dependence and radial flow we shall instead obtain the variance of the multiplicity of particles of species A directly fromφ by employing the more general expressions in Eqs. (50) and (51). As we did in Section III, we shall compute ω A , the ratio of the variance of the multiplicity of species A to its mean, see Eq. (35). We can obtain the mean multiplicity of protons and pions for a rapidity acceptance window [−y max , y max ] and acceptance cuts in p T using the Cooper-Frye formula  In all calculations, we choose the acceptance cuts p T (0.4 GeV, 2 GeV) and y max = 1. As already discussed above,ω A , should not depend on the acceptance. This is explicitly the case in equilibrium and we have verified that this remains approximately the case in all of our simulations. Our approach generalizes the well-known and well-tested Cooper-Frye freezeout [41], which translates hydrodynamic degrees of freedom (but not their fluctuations) into particle distributions. The Cooper-Frye procedure only specifies event-averaged single-particle quantities (multiplicities and spectra of each hadron species) and as such is not sufficient to describe the freezeout of fluctuations or correlations. Our more general freezeout procedure allows us to perform such a translation, or matching, of Hydro+ (hydrodynamic and fluctuation) degrees of freedom to particle multiplicities and their fluctuations.
We demonstrate our generalized freezeout procedure in practice by freezing out a simpli- the radial flow, the typical thermal velocity spread of the produced particles, as well as the acceptance window in rapidity if this acceptance window is larger than the typical thermal spread in rapidity. This characteristic Q is small compared to microscopic scales as it is typically of order 1/τ f (see Eq. (58)). Since the fluctuations at small Q are suppressed by conservation laws (another aspect of the out-of-equilibrium dynamics that is in a sense also a "memory" effect), the smallness of the characteristic Q relevant for the freezeout contributes to the suppression of the fluctuations relative to the equilibrium values.
Our study focused on Gaussian measures of fluctuations. The higher, non-Gaussian, cumulants are more sensitive to the proximity of the critical point [9,11]. It is, therefore, important to generalize our freezeout procedure to higher-order cumulants. This can be done using Eq. (28) and we leave implementation and the analysis of the results to future work. While we demonstrated the application of our procedure to the variance of particle multiplicities, it will also be straightforward to generalize to the cross-correlation of different particle species, as was done in Ref. [11] for equilibrium fluctuations. We expect the conclusion from that work that cumulants involving protons are most sensitive to critical fluctuations to persist, but we leave an investigation of how best to combine measurements of different (cross-)correlations so as to optimize the sensitivity to critical fluctuations while reducing dilution of their effects by backgrounds to future work.
We have focused on the dependence of the observable fluctuations on the proximity of the critical point, either by varying ξ max , which corresponds to varying freezeout µ B via changing collision energy √ s (see Fig. 1), or by varying the freezeout temperature (for the same trajectory). We also studied the dependence on the (thus far unknown) value of the diffusion parameter D 0 . In order to illustrate these dependences, we chose to present our results using normalized variables which did not depend on the absolute magnitude of the effect. In order to predict the absolute magnitude one would have to know the equation of state, i.e., the parameter Z in Eq. (22) which determines the absolute strength of the singularity. Also, one would have to determine the couplings g A of hadrons to the critical collective field σ. These couplings are, of course, also related to the equation of state and it would be interesting to make this relation more explicit. We leave this to future work.
For the purpose of using our results to make crude estimates, one could take Z ∼ 1/T 2 c and follow Ref. [11] and choose g p ∼ 7 and g π ∼ 2.
It would be also interesting to consider going beyond the leading critical behavior of fluctuations to take into account less singular critical contributions and modes which are not critical, including fluctuations of pressure and flow velocity. Extending our work in this way would be necessary in order to attempt to develop a fluctuation freezeout procedure that could be applied away from the critical point as well as near it, as in our paper.
It is also important to realize that we considered fluctuations and freezeout on the crossover side of the critical point (see Fig. 1). It would also be interesting and important to understand what happens on the other side, where the first-order phase transition occurs. The challenge in this domain begins already at the level of hydrodynamics and is beyond the scope of this paper.
Although, as we have detailed in this Section, there is a scope for improvement and generalization of the freezeout procedure that we have introduced and explored in this paper, we believe that the procedure can already be integrated into the full numerical simulation of heavy-ion collisions relevant for the beam energy scan (BES) program aimed at the search for the QCD critical point. With first results from high-statistics BES data taken at RHIC in 2019-2021 anticipated soon, this represents a high priority next step.  [18]. However, it has been observed in Refs.
[1] and [30] that these backreaction effects are smaller than 1 % in most cases. For this reason, throughout the present work we neglect the feedback of Hydro+ modes on the equation of state. That is, in the notation of Refs. [1, 18,30] we approximate the Hydro+ equation of state p + (ε) by the standard pressure p(ε) given by where ε is the energy density, β is the inverse temperature, and where the entropy density s is given as a function of the local temperature by with c V (T ) being the specific heat capacity at fixed volume. In Ref.
[1], the equation of state is specified close to and away from a critical temperature T c by choosing with (T L , T H ) = (T c − ∆T, T c + ∆T ) and where ∆T , which parametrizes the width of the critical region, is the same parameter that arises in Eq. (11). Here as there, we take ∆T = T c /5. Following Ref.
[1], we take c crit V (T ), the critical part of c V that shows the leading singular behavior near the critical point, to have the form where the temperature dependence of the correlation length of critical fluctuations, ξ(T ), needs to be specified. Following Ref.
[1], we do so as in Eq. (11) with a L = 0.1 a QGP , a H = 0.8 a QGP , and a QGP ≡ 4π 2 (N 2 where N c = 3 and N f = 3 are the number of flavors and colors respectively, and with The specification of the equation of state is completed by choosing the six constant coefficients c n that appear in Eq. (A3) so as to enforce that c V (T )/T 3 and its first two derivatives are continuous at T = T L and at T = T H .

Appendix B: Analytical calculations in a Bjorken scenario
In this Appendix, we derive an explicit expression for the numerator of C A (∆y) defined in Eq. (56), namely δ dN A dy + δ dN A dy − σ , in the Bjorken symmetric background described in Section III. We shall begin from the somewhat formal expression for δN 2 A σ that we obtained in Eq. (50), develop an explicit expression for this measure of fluctuations in the Bjorken background, and then see that we can obtain the explicit expression (57) for the rapidity correlator δ dN A dy + δ dN A dy − σ essentially by inspection of the explicit form of δN 2 A σ . In the second half of the Appendix, as a bonus we shall use the analytic control that we have over the calculation of δN 2 A σ to show that the low-Q modes of φ Q make the dominant contribution to this observable.
We begin by noting that in the Bjorken scenario we have the following expressions for x and ∆x ⊥ at τ = τ f : Eqs. (B1) are exact for two points x + and x − on the freeze-out hypersurface and can be obtained by substituting u τ = 1, u r = 0, τ = τ f and ∆τ = 0 in Eqs. (46). We have used x T and ∆x T to denote the transverse parts of x and ∆x ⊥ , namely their projections onto the plane spanned byr andφ. Note that because points on the freeze-out surface, including x + and x − , all have the same τ in the Bjorken scenario, ∆x ⊥ = ∆x in this setting and Eq. (50) need not be corrected as described around Fig. 3.
Next, we note that theφ that arises in the expression (50) that we wish to evaluate is the inverse Wigner transform of φ Q , see Eq. (21). In the Bjorken scenario, this transform takes the formφ where Q η ≡ Q ·η and Q ⊥ ≡ Q − Q ηη . Note that since points on the freeze-out hypersurface all have the same τ in the Bjorken scenario, the expression for the two-point correlator of σ between two points on the freeze-out hypersurface given by Eq. (31) does not receive the correction described following that equation. Note that φ Q that enters Eq. (B2) is obtained by solving Eq. (53), as discussed in Section III. Later in this Appendix, we shall need the formal solution to Eq. (53) that satisfies the initial conditions (54); it is given by The functional form for the evolution of temperature T (τ ), which also determines ξ(τ ) through Eq. (11), can be obtained from the condition that follows from the isentropic nature and Bjorken symmetry of the flow and that must be satisfied for all τ i < τ < τ f . And, we employ the equation of state s(T ) from Ref.
[1] that we describe briefly in Appendix A.
With these preliminaries in place, we now substitute Eq. (B2) into Eq. (50) and integrate over dx T , d∆x T and dQ ⊥ , obtaining where we have defined Q ≡ Q ηη , where A ⊥ is the transverse area in the plane spanned bŷ r andφ, and where we have defined Upon explicit evaluation, this function is given by where Using the expressions above, Eq. (B5) can be evaluated directly, numerically. However, to elucidate its main features we ignore the subleading corrections due to the curvature of the freeze-out hypersurface and assume that m A T , in both cases as discussed in Section II C.
This allows us to make the following approximations: The assumption m A T has allowed us to simplify Eq. (B10b) by expanding only the exponential term in F A as a function of ∆η and not its prefactor. We have verified by explicit calculation that this assumption is well justified for protons. With the above simplifications, after defining ∆y ≡ y + − y − and redefining the variables η and ∆η according to η → η − (y + + y − )/2 and ∆η → ∆η − ∆y, Eq. (B5) becomes This explicit expression for the observable measure of the fluctuations δN 2 A σ is the first main result of this Appendix.
Upon inspection of the result (B11), we see that the two point rapidity space correlator occurring in C A (∆y) is given by Eq. (B13), our second main result of this Appendix, is reproduced in Section III as Eq. (57) and its implications are discussed there.
In the remainder of this Appendix, we shall demonstrate that the low Q modes of φ Q contribute the most to the variance of particle multiplicities, δN 2 A σ . We shall expand φ Q , given by Eq. (B3), in powers of Q to O(Q 2 ) and compare the result for δN 2 A σ that we obtain starting from this expansion to the result that we obtain starting from the full form of φ Q . We denote the polynomial expansion for φ Q to quadratic order by The expression (B14) is a good approximation to φ Q for its low-Q modes, as we illustrate in Fig. 15a.
Upon making the low-Q approximation and working to order Q 2 as in Eq. (B14), we can perform the Q η integral in Eq. (B13), obtaining In Fig. 15b, we compareC p (∆y) (defined via Eqs. (56,59,60)) obtained from δ dN A dy + δ dN A dy − σ computed without making a low-Q approximation, namely Eq. (B13), which is plotted as the solid curves in Fig. 15b, to that computed upon working only to order Q 2 , namely Eq. (B16), which is plotted as the dashed curves in Fig. 15b. The qualitative, even semi-quantitative, agreement between them indicates that the low-Q modes contribute significantly to the variance of particle multiplicities.

Appendix C: Contrasting with Model A evolution
In this Appendix we repeat the calculations of Section IV in a scenario in which the hydrodynamic background is the same (and hence also the same as in Ref. [1]) but in which the dynamical evolution of the fluctuations differs. In Section IV we take conservation laws Normalized fluctuation measure observable (rapidity space correlator) for protonsC p (∆y) obtained with the full form (solid) and truncated form (dashed) of φ Q . The qualitative and even semi-quantitative agreement between the same colored curves in the right plot indicates that the low-Q modes contribute significantly to the variance of particle multiplicities. In obtaining these plots, ξ max was set to 3 fm and the fluctuations at τ i = 1 fm were initialized to their equilibrium value at τ = τ i = 1 fm with T i (τ i ) = 235 MeV.
into account, following Model H dynamics. Here, in contrast, we shall consider the case where the fluctuating slow mode is not a conserved quantity, meaning that the appropriate dynamics for the relaxation of the two-point function is that of Model A in the classification of Halperin and Hohenberg [51], with the relaxation rate given by Eq. (10), which we repeat here: We have performed simulations with Γ 0 = 1 fm −1 and 8 fm −1 , which correspond to Γ 0 ξ 2 0 = 0.25 fm and 2 fm, respectively. Comparing the "shapes" of all the results plotted in this Appendix to those in the analogous Figures in Section IV provides us with another way of seeing the impact of conservation laws on the results from Section IV. The main difference between Model A and Model H evolution arises from the qualitatively different relaxation rate for the low Q modes, Qξ 1, which goes as Γ 0 ξ 2 0 /ξ 2 in Model A and as (D 0 ξ 0 /ξ)Q 2 in Model H, with the Q 2 -suppression being a manifestation of conservation. A second motivation for this Appendix is that, because Model A dynamics is simpler to implement, in their pioneering calculation the authors of Ref.
[1] used Model A dynamics, meaning that in this Appendix we shall be freezing out the calculations of Ref.

Evolution of φ Q
In Fig. 16, which can be compared to the analogous Model H results shown in Fig. 9, we illustrate the Model A dynamics of φ Q for fluid cells following two different hydrodynamic flow lines, with two choices of Γ 0 and two choices of ξ max . As in Section IV, varying ξ max corresponds to varying how close the cooling trajectory of the fluid cell comes to the critical point on the phase diagram. The central qualitative difference between the Model A results in Fig. 16 and the Model H results in Fig. 9 is that in Model H φ Q at Q = 0 is unchanging in time, because of conservation, which means that as critical flucutations develop we see that in Fig. 9 φ Q takes on a shape in which it first rises as a function of increasing Q and then falls whereas in the Model A dynamics of this Appendix the maximum value of φ Q is found at Q = 0, and this value is time dependent. In Model A, here, as in Model H in Fig. 9, the fluctuations φ Q fall out of equilibrium, lagging behind the equilibrium fluctuationsφ Q as the latter change with time.
In Fig. 17, which can be compared to the analogous Model H results shown in Fig. 10, for all three representative Q modes that are plotted we notice the φ Q 's lagging behind their respectiveφ Q s, with the degree to which they fall out of equilibrium greater for smaller Γ 0 , meaning slower relaxation toward equilibrium. For the values of Γ 0 that we have considered in Fig. 17, we can see that fluctuations do depend on whether we choose a freeze-out temperature of 156 MeV or 140 MeV. As we also observed in Fig. 10, φ Q has an inflection point at T = T c where the relaxation rate takes its minimum value and the growth of φ Q stops when φ Q equals the instantaneousφ Q .

Fluctuations on the freezeout surface
In Fig. 18, which can be compared to Fig. 11, suitably normalized plots of φ Q are shown for three points on the freeze-out hypersurface, characterized by radial coordinate r = 0, 3 As in Fig. 12, in Fig. 19 we have computedφ(∆x ⊥ ), the inverse Fourier transform of φ Q defined in Eq. (32), and plotted ∆x 2φ (∆x ⊥ ) as a function of the spatial separation ∆x between the two points in the correlator δŝ(x + )δŝ(x − ) . As in the Model H evolution of Model H dynamics in Fig. 12 is that hereφ(∆x ⊥ ) is positive at large ∆x: the fact that it becomes negative in the large ∆x region in Fig. 12 is a direct consequence of conservation in Model H.

Variance of particle multiplicities
As in Section IV, but here with Model A dynamics, we close by computing and plotting the normalized fluctuation measure for protons and pions,ω p andω π , in Figs. 20 and 21. Figs. 13 and 14, to which these figures can be compared, these results demonstrate that for trajectories passing closer to the critical point (i.e., for trajectories with larger ξ max ) the magnitude of fluctuations is larger. Again as in Section IV, the magnitude of the effect depends on the rate of the relaxation of the fluctuations, controlled here by parameter Γ 0 . We see that for large enough values of Γ 0 , eg. 8 fm −1 , the proton and pion fluctuations are , as a function of the maximum equilibrium correlation length along the system trajectory, which is to say as a function of how closely the trajectory passes the critical point. As Γ 0 → ∞, theω π 's approach their equilibrium values shown in panel (c).

As in
able to come reasonably close to their equilibrium values, which as an aside means that they depend quite sensitively on the freeze-out temperature. The differences that we discussed at length in Section IV originate from the effects of conservation.