text:
# Hydrodynamics Near a Chiral Critical Point

## I Introduction

The hydrodynamical model is frequently employed to describe multiparticle production processes in hadronic collisions . In particular, it predicts characteristic flow-signatures as “fingerprints” for non-trivial equations of state of hot and dense matter . Such equations of state can occur when the effective potential, obtained by integrating out some degrees of freedom, exhibits features characteristic of a phase transition in thermodynamics .

More specifically, if there exist two (or more) collective states with the same free energy but separated by a barrier, then behavior characteristic of a first-order phase transition may emerge. On the other hand, if no free-energy barrier exists, one might expect resemblence to a second-order phase transition. This analogy of interacting quantum field theories with thermodynamics is believed to have played an important role in the evolution of the early universe  and is currently being investigated in accelerator experiments by colliding beams of protons and heavy ions . Classical energy flow and hydrodynamic scaling behavior emerges in high-energy inclusive processes  and from the real-time evolution of some quantum field theories .

In this paper, we extend the hydrodynamical transport model such that phase transitions related to the restoration (or breaking) of some global symmetry can be studied dynamically. In particular, we focus on chiral symmetry breaking at finite temperature , for which we shall adopt a relatively simple and tractable phenomenological model, i.e. the Gell-Mann and Levy model .

It has been argued  that the chiral phase transition for two massless quark flavors is second-order at baryon-chemical potential $`\mu _B=0`$, which then becomes a smooth crossover for small quark masses. On the other hand, a first-order phase transition is predicted for small temperature $`T`$ and large $`\mu _B`$. If, indeed, there is a smooth crossover for $`\mu _B=0`$ and high $`T`$, and a first-order transition for small $`T`$ and high $`\mu _B`$, then the first-order phase transition line in the $`(\mu _B,T)`$ plane must end in a second-order critical point. This point was estimated  to be at $`T100`$ MeV and $`\mu _B600`$ MeV (see also ). More recently, it has been attempted to determine the endpoint of the line of first-order phase transitions from the lattice, using $`2+1`$ quark flavors on $`N_t=4`$ lattices  (see also ). Those authors locate the critical point at $`T=160\pm 3.5`$ MeV and $`\mu _B=725\pm 35`$ MeV. Note, however, that a reliable extrapolation to the continuum limit and to physical pion mass has not been attempted so far.

There is an ongoing experimental effort to detect that chiral critical point in heavy-ion collisions at high energies. Note that both, high temperature and high baryon density are required to have dynamical trajectories in heavy-ion collisions pass reasonably close by the critical point. Some dynamical computations  of the energy deposition and baryon stopping process during the initial stage of head-on collisions of large nuclei within semi-realistic multi-fluid dynamical models suggest that the required conditions may be reached in the central region of collisions at $`E_{lab}2080`$ AGeV on a fixed target, or in the fragmentation regions of collisions at higher energies. However, to our knowledge there has been so far no attempt to describe hydrodynamical expansion of the hot and dense droplet produced initially for dynamical trajectories close to the critical point. This paper represents an attempt in that direction.

## II The Model

In this section we shall present our model for the dynamics of a droplet of quarks and antiquarks, starting at high temperature in a state with (approximately) restored chiral symmetry, and evolving towards a state where the symmetry is spontaneously broken. The quarks will be described as a heat bath in local thermal equilibrium that evolves in 3+1 dimensions according to the conservation laws for energy and momentum, i.e. relativistic hydrodynamics. However, the “fluid” of quarks interacts locally with the chiral fields, that is, they can exchange energy and momentum. In turn, the (long wavelength modes of the) fields obey the classical equations of motion which follow from the underlying Lagrangian in the presence of the quarks and anti-quarks. Similar models for the dynamics of quarks coupled to chiral fields were considered in the past. In , a background of freely streaming quarks was assumed, and the classical evolution of the chiral fields was discussed. More realistic dynamical descriptions for the quark medium followed shortly, treating them as either a relativistic fluid , as also envisaged here, or within classical Vlasov transport theory . Those studies focused on a second-order chiral phase transition, or, in the presence of explicit symmetry breaking, on a smooth crossover. However, it turns out that one can also address first-order chiral phase transitions within the very same model, at least in a phenomenological fashion, by chosing larger quark-field coupling $`g`$  (see below). Integrating out the quarks then leads to an effective potential exhibiting two degenerate states around $`T_c`$.

Here, we extend the previous work mentioned above, and at the same time shift our focus somewhat. Namely, the early studies were mainly concerned with the dynamical evolution of the long-wavelength chiral fields, and of classical pion production; that is, they mainly adressed issues related to the possiblity of forming “Domains of Disoriented Chiral Condensates” (DCC), as suggested by Rajagopal and Wilczek, and others , see also . Our present work puts more emphasis on the dynamics of the heat-bath of quarks, rather than on that of the soft modes of the chiral field. We shall point out qualitative changes in the classical energy-momentum flow of the “fluid” of quarks in the proximity of a chiral critical point, rather than look for “rare phenomena” like DCC formation.

Moreover, refs.  all employed the mean field approximation for the chiral fields. Field fluctuations at the phase transition were not considered. As an example, for the first-order phase transition discussed in  dynamical bubble nucleation (“boiling”) could not be described, as it requires large coherent thermal field fluctuations from the symmetry restored phase, over the barrier and into the symmetry broken phase (see e.g.  for results of such dynamical simulations, and  for a computation of bubble nucleation rates from the linear sigma model). Thus, the main improvement here is that we do include a dynamical treatment of field fluctuations in the vicinity of a critical point, and their influence on the dynamical evolution of the quark fluid.

On the technical side, going beyond the mean field approximation requires us to introduce appropriate subtractions in all thermodynamical functions, as explained in appendix A. Moreover, the coupled system of non-linear partial differential equations has to be solved numerically in 3+1 dimensions, without imposing any space-time symmetry assumptions (while  all simplified the solution greatly by assuming special symmetries which essentially reduced the problem to 0+1 or 1+1 space-time dimensions). That is because fluctuations break any space-time symmetry that may be obeyed by the mean field, as for example spherical symmetry or symmetry under Lorentz boosts in a particular direction.

As mentioned in the introduction, physically the chiral critical point is expected to occur for some specific values of temperature $`T`$ and baryon-chemical potential $`\mu _B`$. To simplify the problem and its numerical solution, however, here we rather choose to consider only locally baryon symmetric matter, i.e. equal numbers of quarks and anti-quarks. Instead, we can “shift” the critical point by varying the quark-field coupling constant $`g`$, that is, by increasing or decreasing the vacuum mass of the constituent quarks. As explained in more detail below, large values for $`g`$ result in a first-order phase transition, while small $`g`$ leads to a crossover. For the observables discussed in section III, the qualitative difference between the two realizations of a chiral critical point mentioned above should not matter much***However, we note that there are indeed differences on the quantitative level. For example, for a first-order phase transition in baryon dense matter the isentropic speed of sound does not vanish in general at $`T_c`$ , as it does for zero net baryon charge, $`\mu _B=0`$.. On the other hand, our simplified treatment disables us from studying fluctuations of net baryon charge .

In section II A we discuss the effective potential “seen” by the long wavelength modes of the chiral fields in the presence of a heat bath of quarks and anti-quarks. Then, in II B we present our model for the non-equilibrium dynamical treatment of field and fluid evolutions. Some numerical algorithms and details are mentioned briefly in section II C. In section III we present our results and end with a summary and an outlook in section IV.

### A Effective Potential

As an effective theory of the chiral symmetry breaking dynamics, we assume the linear $`\sigma `$-model coupled to two flavors of quarks :

$``$ $`=`$ $`\overline{q}\left[i\gamma ^\mu _\mu g(\sigma +\gamma _5\stackrel{}{\tau }\stackrel{}{\pi })\right]q`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}\left(_\mu \sigma ^\mu \sigma +_\mu \stackrel{}{\pi }^\mu \stackrel{}{\pi }\right)U(\sigma ,\stackrel{}{\pi }).`$ (2)
The potential, which exhibits both spontaneously and explicitly broken chiral symmetry, is

$$U(\sigma ,\stackrel{}{\pi })=\frac{\lambda ^2}{4}(\sigma ^2+\stackrel{}{\pi }^2v^2)^2h_q\sigma U_0.$$
(3)
Here $`q`$ is the constituent-quark field $`q=(u,d)`$. The scalar field $`\sigma `$ and the pseudoscalar field $`\stackrel{}{\pi }=(\pi _1,\pi _2,\pi _3)`$ together form a chiral field $`\varphi _a=(\sigma ,\stackrel{}{\pi })`$. The parameters of the Lagrangian are chosen such that chiral $`SU_L(2)SU_R(2)`$ symmetry is spontaneously broken in the vacuum. The vacuum expectation values of the condensates are $`\sigma =f_\pi `$ and $`\stackrel{}{\pi }=0`$, where $`f_\pi =93`$ MeV is the pion decay constant. The explicit symmetry breaking term is due to the non-zero pion mass and is determined by the PCAC relation, which gives $`h_q=f_\pi m_\pi ^2`$, where $`m_\pi =138`$ MeV. This leads to $`v^2=f_\pi ^2m_\pi ^2/\lambda ^2`$. The value of $`\lambda ^2=20`$ leads to a $`\sigma `$-mass, $`m_\sigma ^2=2\lambda ^2f_\pi ^2+m_\pi ^2`$, approximately equal to 600 MeV. In mean field theory, the purely bosonic part of this Lagrangian exhibits a second-order phase transition  if the explicit symmetry breaking term, $`h_q`$, is dropped. For $`h_q0`$, the transition becomes a smooth crossover from the phase with restored symmetry to that of broken symmetry. The normalization constant $`U_0`$ is chosen such that the potential energy vanishes in the ground state, that is

$$U_0=\frac{m_\pi ^4}{4\lambda ^2}f_\pi ^2m_\pi ^2.$$
(4)
For $`g>0`$, the finite-temperature one-loop effective potential also includes a contribution from the quarks. In our phenomenological approach, we shall consider the quarks as a heat bath in (local) thermal equilibrium. Thus, it is possible to integrate them out to obtain the effective potential for the chiral fields in the presence of that bath of quarks. Consider a system of quarks and antiquarks in thermodynamical equilibrium at temperature $`T`$ and in a volume $`𝒱`$. The grand canonical partition function reads

$$𝒵=𝒟\overline{q}𝒟q𝒟\sigma 𝒟\stackrel{}{\pi }\mathrm{exp}\left[_0^{1/T}d(it)_𝒱d^3\stackrel{}{x}\right].$$
(5)
(Anti-) Periodic boundary conditions for the (fermion) boson fields are implied. In mean-field approximation the chiral fields in the Lagrangian are replaced by their expectation values, which we denote by $`\sigma `$ and $`\stackrel{}{\pi }`$. Then, up to an overall normalization factor,

$`𝒵`$ $`=`$ $`𝒩_U{\displaystyle }𝒟\overline{q}𝒟q\mathrm{exp}\{i{\displaystyle _0^{1/T}}dt{\displaystyle _𝒱}d^3\stackrel{}{x}`$ (7)
$`\overline{q}[i\gamma ^\mu _\mu g(\sigma +i\gamma _5\stackrel{}{\tau }\stackrel{}{\pi })]q\}`$
$`=`$ $`𝒩_U\mathrm{det}_p\left\{\left[p_\mu \gamma ^\mu g(\sigma +i\gamma _5\stackrel{}{\tau }\stackrel{}{\pi })\right]/T\right\},`$ (8)
where

$$𝒩_U=\mathrm{exp}\left(\frac{𝒱U(\sigma ,\stackrel{}{\pi })}{T}\right).$$
(9)
Taking the logarithm of $`𝒵`$, the determinant of the Dirac operator can be evaluated in the standard fashion , and we finally obtain the grand canonical potential

$$\frac{T}{𝒱}\mathrm{log}𝒵=U+\stackrel{~}{V}_q(T),$$
(10)
where

$$\stackrel{~}{V}_q(T)=d_q\frac{d^3\stackrel{}{p}}{(2\pi )^3}\left\{E+T\mathrm{log}\left[1+e^{E/T}\right]\right\}.$$
(11)
Here, $`d_q=24`$ denotes the color-spin-isospin-baryon charge degeneracy of the quarks. For our purposes, the zero-temperature contribution to $`\stackrel{~}{V}_q`$, i.e. the first term in the integral in (11), can be absorbed into $`U`$ via a standard renormalization of the bare parameters $`\lambda ^2`$ and $`v^2`$. The logarithmic dependence on the renormalization scale is ignored in the following.

Adding the $`T=0`$ and the finite-$`T`$ contributions defines our effective potential $`V_{\mathrm{eff}}`$:

$$V_{\mathrm{eff}}(\varphi _a,T)U(\varphi _a)d_qT\frac{d^3p}{(2\pi )^3}\mathrm{log}\left(1+e^{E/T}\right).$$
(12)
$`V_{\mathrm{eff}}`$ depends on the order parameter field through the effective mass of the quarks, $`m_q=g\sqrt{\varphi ^2}g\sqrt{_a\varphi _a\varphi _a}`$, which enters into the expression for the energy, $`E=\sqrt{p^2+g^2\varphi ^2}`$.

In principle, one should also integrate out short wavelength fluctuations of $`\varphi `$, which would lead to an additional contribution to the effective potential (12), see for example . For the present phenomenological analysis, however, the expression (12) is already sufficient in that it exhibits a critical endpoint for some particular value of the coupling constant $`g`$. Since we do not expect the simple model (2) to be quantitatively reliable anyway, it is not unreasonable to employ the effective potential (12) for a study of qualitative effects near a chiral critical point.

We now turn to a discussion of the shape of the effective potential, cf. Fig. 1. For sufficiently small $`g`$ one still finds the above-mentioned smooth transition between the two phases. At larger coupling to the quarks, however, the effective potential exhibits a first-order phase transition . Along the line of first-order transitions, for temperatures near the critical temperature, $`V_{\mathrm{eff}}`$ displays a local minimum $`\sigma =\sigma _1(T)0`$ which is separated by a barrier from another local minimum at $`\sigma =\sigma _2(T)>0`$. (There is another local minimum for negative $`\sigma `$ which is of higher energy and does not concern us.) These two minima are degenerate at $`T=T_c`$. For example, $`g=5.5`$ leads to a critical temperature of $`T_c123.3`$ MeV. Lowering the value of $`g`$ leads to a smaller barrier between the two degenerate states. Also, $`\sigma _1`$ approaches $`\sigma _2`$, i.e., the phase transition weakens, and moreover the spinodal temperature approaches $`T_c`$ . At $`g_c3.7`$, finally, the barrier disappears, and so the latent heat vanishes. This is the second-order critical point, where the potential about the minimum is flat.

A different possibility of making the $`\sigma `$-meson much lighter at $`T_c`$ than at $`T=0`$ is to reduce the self-coupling of the chiral fields , $`\lambda ^2`$, rather than that to the quarks. For example, one may choose $`\lambda ^22.2`$ with the pion decay constant $`f_\pi `$ and vacuum mass $`m_\pi `$ fixed, such that $`v^2=0`$. However, within the present model this also reduces $`T_c`$ significantly, to less than $`100`$ MeV . Such low $`T_c`$ appear to be excluded by present lattice QCD results , and moreover would generate too small thermal fluctuations in the heat bath. Therefore, we keep $`\lambda ^2=20`$ fixed throughout the manuscript.

### B Coupled Dynamics of Fields and Fluid

The classical equations of motion for the chiral fields are

$`_\mu ^\mu \sigma +{\displaystyle \frac{\delta U}{\delta \sigma }}`$ $`=`$ $`g\overline{q}q=g\rho _s,`$ (13)
$`_\mu ^\mu \stackrel{}{\pi }+{\displaystyle \frac{\delta U}{\delta \stackrel{}{\pi }}}`$ $`=`$ $`g\overline{q}\gamma _5\stackrel{}{\tau }q=g\stackrel{}{\rho }_{ps},`$ (14)
where

$`\rho _s`$ $`=`$ $`\overline{q}q=g\sigma d_q{\displaystyle \frac{\mathrm{d}^3p}{(2\pi )^3}\frac{1}{E}f(p)},`$ (15)
$`\stackrel{}{\rho }_{ps}`$ $`=`$ $`\overline{q}\gamma _5\stackrel{}{\tau }q=g\stackrel{}{\pi }d_q{\displaystyle \frac{\mathrm{d}^3p}{(2\pi )^3}\frac{1}{E}f(p)}`$ (16)
are the scalar and pseudoscalar densities generated by the heat bath of quarks and anti-quarks, respectively. The distribution of the quarks and anti-quarks in momentum space is given by the Fermi-Dirac distribution.

The form (14) for the coupling to the heat bath of quarks can be derived from the Lagrange density (2) in mean-field approximation. The energy density of the quarks is given by

$$_q=\overline{q}i/q+g\overline{q}\sigma q+g\overline{q}\gamma _5\stackrel{}{\tau }\stackrel{}{\pi }q.$$
(17)
The finite-$`T`$ contribution to the equation of motion for the $`\sigma `$, $`\stackrel{}{\pi }`$ fields is obtained from the variation of the effective potential with respect to $`\sigma `$ or $`\stackrel{}{\pi }`$, respectively:

$`g\overline{q}q`$ $`=`$ $`{\displaystyle \frac{\delta _q}{\delta \sigma }}={\displaystyle \frac{\delta V_{\mathrm{eff}}U}{\delta \sigma }},`$ (18)
$`g\overline{q}\gamma _5\stackrel{}{\tau }q`$ $`=`$ $`{\displaystyle \frac{\delta _q}{\delta \stackrel{}{\pi }}}={\displaystyle \frac{\delta V_{\mathrm{eff}}U}{\delta \stackrel{}{\pi }}}.`$ (19)
Applying this to the right hand side of equation (12) yields the expressions for the scalar density and the pseudoscalar density as given in (15).

As already mentioned above, we shall assume that the quarks constitute a heat bath in local thermal equilibrium. Thus, their dynamical evolution is determined by the local conservation laws for energy and momentum in relativistic hydrodynamics. For simplicity, we shall further assume that the stress-energy tensor of the quark fluid is of the “perfect fluid” form (corrections could in principle be taken into account in the future along the lines discussed in ),

$$T^{\mu \nu }=(e+p)u^\mu u^\nu pg^{\mu \nu }.$$
(20)
Here,

$$u^\mu \frac{T^{\mu \nu }u_\nu }{\sqrt{u_\sigma T^{\sigma \rho }T_{\rho \alpha }u_\alpha }}$$
(21)
is the local four-velocity of the fluid. $`g^{\mu \nu }=\mathrm{diag}(1,1,1,1)`$ is our metric tensor, and so the line element is $`ds^2=dt^2d\stackrel{}{x}^{\mathrm{\hspace{0.17em}2}}`$. The time $`t`$ is measured in the global restframe. Furthermore,

$`e(\varphi ,T)`$ $`=`$ $`_q,`$ (22)
$`p(\varphi ,T)`$ $`=`$ $`V_{\mathrm{eff}}(\varphi ,T)+U(\varphi _a),`$ (23)
denote the energy density and the pressure of the quarks at temperature $`T`$. Note that we do not assume that the chiral fields are in equilibrium with the heat bath of quarks. Hence, both $`e`$ and $`p`$ depend explicitly on $`\varphi =\sqrt{\sigma ^2+\stackrel{}{\pi }^2}`$. Given an initial condition on some space-like surface, $`T^{\mu \nu }`$ at any other causally connected space-time point is determined by

$$_\mu T^{\mu \nu }=S^\nu .$$
(24)
In the absence of interactions between the chiral fields and the quarks, the source term $`S^\nu `$ vanishes, and the energy and the momentum of the quarks,

$$(E,\stackrel{}{P})=𝑑\sigma _\mu T^{\mu \nu }$$
(25)
are conserved:

$$u_\mu ^\mu E=u_\mu ^\mu P^i=0.$$
(26)
With interactions turned on, this is obviously not true any longer. Rather, the total energy and momentum of fluid plus fields is the conserved quantity:

$$S^\nu =_\mu T_\varphi ^{\mu \nu }.$$
(27)
The stress-energy tensor for the fields can be computed from the Lagrange density in the standard fashion. The effective mass of the quarks, $`m_q=g\sqrt{_a\varphi _a^2}`$, is already accounted for in the EoS for the quark fluid, see eqs. (23). Thus, $`T_\varphi ^{\mu \nu }`$ is the stress-energy tensor of the chiral fields alone:

$`T_\varphi ^{\mu \nu }`$ $`=`$ $`{\displaystyle \underset{a}{}}{\displaystyle \frac{_\varphi }{\left(_\mu \varphi _a\right)}}^\nu \varphi _ag^{\mu \nu }_\varphi `$ (28)
$`_\varphi `$ $`=`$ $`{\displaystyle \underset{a}{}}{\displaystyle \frac{1}{2}}(_\mu \varphi _a)(^\mu \varphi _a)U.`$ (29)
Its divergence is given by

$`_\mu T_\varphi ^{\mu \nu }`$ $`=`$ $`{\displaystyle \underset{a}{}}\left\{_\mu ^\mu \varphi _a+{\displaystyle \frac{\delta U}{\delta \varphi _a}}\right\}^\nu \varphi _a`$ (30)
$`=`$ $`g\rho _s^\nu \sigma +g\stackrel{}{\rho }_{ps}^\nu \stackrel{}{\pi }.`$ (31)
In the second step we made use of the equations of motion (14). This is the source term in the continuity equation (24) for the stress energy tensor of the quark fluid. A different derivation based on moments of the classical Vlasov equation for the quarks was given in , see also .

We emphasize again that we employ eqs. (14) not only to propagate the mean field through the transition but fluctuations as well. The initial condition includes some generic “primordial” spectrum of fluctuations, see section III A, which then evolve in the effective potential generated by the matter fields, i.e. the quarks. Near the critical point, those fluctuations have small effective mass and “spread out” to probe the flat effective potential. Since all field modes are coupled, effectively the fluctuations act as a noise term in the equation of motion for the mean field, similar to the familiar Langevin dynamics ; near the phase transition, however, the noise is neither Gaussian (the effective potential is not parabolic) nor Markovian (zero correlation length in time) nor white (zero correlation length in space). Rather, correlation lengths and $`n`$-point functions of the “noise” are governed by the dynamics of the fluctuations in the effective potential generated by the fluid of quarks, including also effects from the finite size and the relativistic expansion of the system.

### C Numerics and Technical Details

We briefly describe how we solved the coupled system (14,24) of partial differential equations in 3+1 space-time dimensions.

We follow the evolution on $`t=const`$ hypersurfaces, and in a fixed spatial cube of volume $`L^3`$. We discretize three-space in that cube by introducing a $`160^3`$ grid with a spacing of $`\mathrm{\Delta }x=0.2`$ fm (thus, $`L=160\times 0.2`$ fm=32 fm). On that grid, we solve the hyperbolic continuity equations of fluid dynamics (24) using the so-called phoenical SHASTA flux-corrected transport algorithm with simplified source treatment. It is described and tested in detail in Ref. , and we refrain from a discussion here. The time step was chosen as $`\mathrm{\Delta }t=0.4\mathrm{\Delta }x`$, as appropriate for the SHASTA . We performed each time step in the standard fashion with $`S^\nu =0`$, then subtracted the sources $`S^\nu `$ from the energy and momentum density in the calculational frame (i.e. the global rest frame of the fluid). Finally, from $`T^{00}`$, $`T^{0i}`$ and the EoS $`p=p(e,\varphi _a)`$ we solved for the velocity of the local restframe (LRF) of each cell, and for the energy density $`e`$ of the fluid in the LRF. Such a treatment of sources in the continuity equations for energy and momentum proved to work well in relativistic multi-fluid dynamics , where one encounters a similar equation due to interactions between various fluids. The boundary conditions at the edge of the computational grid are such that the fluid simply streams out when reaching the boundary. This can be monitored by checking the conservation of the total energy as a function of time.

Regarding the fields, we solve the classical equations of motion using a staggered leap-frog algorithm with second-order accuracy in time, see for example Ref. . The non-linear wave-equations (14) are split into two coupled first-order equations (in time) by considering separately the field $`\varphi _a`$ and its canonically conjugate field $`d\varphi _a/dt`$. For this algorithm, the time step for propagation of the fields has to be chosen smaller than that for propagation of the fluid. In practice, we found that $`\mathrm{\Delta }t=0.1\mathrm{\Delta }x`$ gave reasonable accuracy. As the time step for propagation of the fluid was chosen to be four times larger, the fields where propagated for four consecutive steps, with the fluid kept “frozen”, followed by one single step for the fluid. Spatial gradients of the fields where set to zero at the edges of the computational grid, rather than employing periodic boundary conditions. The spatial derivatives in the equations of motion for the fields (14) where discretized to second order accuracy in the grid spacing to match the second order accuracy in time. We employed the same small grid spacing $`\mathrm{\Delta }x=0.3`$ fm as for the fluid. Although the mean fields vary on a larger scale, it is important to allow for non-linear amplification of harder fluctuations which can couple to the soft modes and affect the dynamical relaxation to the vacuum state with broken symmetry, see for example  and .

For the fluid, one needs to introduce nine 3d spatial gridsAt finite baryon density, one would need another grid for $`\rho _B`$. Also, inverting the Lorentz transformation to obtain the LRF densities $`e`$ and $`\rho _B`$ from $`T^{00}`$, $`T^{0i}`$ and the net baryon current $`J_B^\mu \rho _Bu^\mu `$ in each time step would require a three-dimensional root search for the function $`p=p(e,\rho _B,\varphi )`$., for $`e`$, $`p`$, $`T^{0\mu }`$, and $`\stackrel{}{u}`$. In addition, one needs two more 3d grids for each field component, namely for $`\varphi _a`$ and $`d\varphi _a/dt`$. To save computational resources, we have therefore decided to neglect the pion field altogether, i.e., to set it to zero everywhere in the forward light cone. Thus, only the $`\sigma `$-field is considered in the actual computations described below. However, this represents a minor restriction only, since we focus here on the bulk evolution of the quark fluid rather than on fluctuation observables or coherent pion production from the decay of a classical pion field. (Coherent pion production is known to contribute a small fraction of the total pion yield only . Incoherent particle production from the decay of the fluid by far dominates, if no kinematic cuts are applied.)

## III Results

Having formulated our model, we now proceed to show some specific examples of numerical solutions. In particular, we would like to examine whether the dynamical evolution changes as one crosses the chiral critical point.

### A Initial Condition

As an example, we employ the following initial conditions for the results described below. Of course, one could employ more refined initial conditions and “tune” them such as to reproduce various aspects of the final state, which in principle could be compared to experimental data. At present, our more modest goal is to illustrate qualitative effects originating from the phase transition. The initial conditions are meant to provide a semi-realistic parameterization of the hot fireball created in a high-energy heavy-ion collision.

At time $`t=0`$, the distribution of energy density for the quarks is taken to be uniform in the $`z`$-direction (with length $`2l_z=12`$ fm) and ellipsoidal in the $`xy`$-plane:

$$e(t=0,\stackrel{}{x})=\{\begin{array}{cc}\hfill e_{\mathrm{eq}}:& x^2b^2+y^2a^2<(ab)^2z<l_z\hfill \\ \hfill 0:& x^2b^2+y^2a^2>(ab)^2z>l_z\hfill \end{array},$$
(32)
where $`e_{\mathrm{eq}}`$ denotes the equilibrium value of the energy density taken at a temperature of $`T_i160`$ MeV. In our calculation we choose $`a=r_A\stackrel{~}{b}/2`$ and $`b=\sqrt{r_A^2\stackrel{~}{b}^2/4}`$, where $`\stackrel{~}{b}`$ denotes the impact parameter of two nuclei with radius $`r_A`$ (below, we choose $`r_A=6.5`$ fm and $`\stackrel{~}{b}=6`$ fm). Thus, the ellipsoidal shape resembles the almond shaped overlap of two colliding nuclei. The $`\mathrm{\Theta }`$-function distribution of the initial energy density is evidently somewhat unrealistic; a smoother distribution with non-zero surface thickness would be more realistic and perhaps affect the results somewhat. However, as already mentioned above, here we do not aim at quantitative fits to experimental data but at illustrating qualitative effects related to the shape of the effective potential in the transition region.

The collective longitudinal velocity of the fluid of quarks is assumed to rise linearly with $`z`$: $`v_z(t=0,\stackrel{}{x})z/l_zv_{max}`$, where $`v_{max}=0.2`$ . The transverse components of $`\stackrel{}{v}`$ are set to zero at $`t=0`$.

Our initial conditions for the chiral fields are

$`\sigma (t=0,\stackrel{}{x})`$ $`=`$ $`\delta \sigma (\stackrel{}{x})+f_\pi +`$ (35)
$`(f_\pi +\sigma _{\mathrm{eq}})\left[1+\mathrm{exp}\left({\displaystyle \frac{\stackrel{~}{r}\stackrel{~}{R}}{\stackrel{~}{a}}}\right)\right]^1`$
$`\left[1+\mathrm{exp}\left({\displaystyle \frac{|z|l_z}{\stackrel{~}{a}}}\right)\right]^1,`$
$`\stackrel{}{\pi }(t=0,\stackrel{}{x})`$ $`=`$ $`\delta \stackrel{}{\pi }(\stackrel{}{x}),`$ (36)
where $`\stackrel{~}{r}=\sqrt{x^2+y^2}`$,

$$\stackrel{~}{R}=\{\begin{array}{cc}\hfill \frac{ab\stackrel{~}{r}}{\sqrt{b^2x^2+a^2y^2}}:& \stackrel{~}{r}0\hfill \\ \hfill a:& \stackrel{~}{r}=0\hfill \end{array},$$
(37)
and $`\stackrel{~}{a}=0.3`$ fm is the surface thickness of this Woods-Saxon like distribution. Here $`\sigma _{\mathrm{eq}}0`$ is the value of the $`\sigma `$ field corresponding to $`e_{\mathrm{eq}}`$. Thus, the chiral condensate nearly vanishes at the center, where the energy density of the quarks is large, and then quickly interpolates to $`f_\pi `$ where the matter density is low.

$`\delta \sigma (\stackrel{}{x})`$ and $`\delta \stackrel{}{\pi }(\stackrel{}{x})`$ represent the initial random fluctuations of the fields. Our focus at this stage is on how those “primordial” fluctuations evolve through the phase transition (or crossover) and how they affect the hydrodynamic expansion of the thermalized matter fields (the fluid). Thus, for a first qualitative analysis we do not rely on additional physics inputFor example, one might assume thermalized primordial fluctuations, in which case their distribution depends on the effective potential at the temperature $`T_i`$. We have checked, however, that at high temperature $`V_{\mathrm{eff}}(\varphi ,T_i)`$ looks rather similar for all values of $`g`$ considered here, regardless of whether later on the evolution proceeds through the crossover, the critical endpoint or the first-order phase transition regimes. for the primordial fluctuations but rather chose a generic Gaussian distribution,

$$P[\delta \varphi _a]\mathrm{exp}\left(\delta \varphi _a^2/2\delta \varphi _a^2\right),$$
(38)
with the variance $`\delta \varphi _a^2`$ left as a free parameter. (Here, no summation over the index $`a`$ (internal space) is carried out.) We performed simulations with $`\sqrt{\delta \varphi _0^2}\sqrt{\delta \sigma ^2}=v/3`$. Such fluctuations are large enough to probe the barrier of the effective potential in case of a first order transition, or the flat region close to the critical endpoint (Fig. 1). On the other hand, they are sufficiently small to allow for a one-loop subtraction of their contribution to the effective potential, as described in appendix A.

The time derivatives of the fields where set to zero at $`t=0`$, with no fluctuations. As already mentioned in section II C, the actual numerical computations described below were performed with $`\delta \stackrel{}{\pi }0`$, i.e. $`\delta \varphi _{1,2,3}^2=0`$.

One must further take into account that the field fluctuations are correlated over some space-like distance $`\xi 1`$ fm. This is a physical scale which is present in the initial conditions; if one simply picks random fluctuations at each point of the grid, the correlation length will instead be given by the artificial numerical discretization scale $`\mathrm{\Delta }x`$.

In practice, a useful and simple procedure to implement physical correlations in the initial condition is as follows. First, at each point of the grid one samples the distribution (38) at random. Then, one smoothly sweeps a ”sliding window” of linear size $`\xi =n\mathrm{\Delta }x`$ over the grid and averages the fields:

$`\varphi _a^{}(\stackrel{}{x})`$ $`=`$ $`{\displaystyle \frac{1}{n^3}}`$ (40)
$`\times {\displaystyle \underset{i,j,k=0..n1}{}}\varphi _a(\stackrel{}{x}+i\mathrm{\Delta }x\stackrel{}{e}_1+j\mathrm{\Delta }x\stackrel{}{e}_2+k\mathrm{\Delta }x\stackrel{}{e}_3).`$
Here, $`\stackrel{}{e}_1`$, $`\stackrel{}{e}_2`$, $`\stackrel{}{e}_3`$ define a global cartesian orthonormal basis, and $`\stackrel{}{x}=1\stackrel{}{e}_1+1\stackrel{}{e}_2+1\stackrel{}{e}_3`$, $`1\stackrel{}{e}_1+1\stackrel{}{e}_2+2\stackrel{}{e}_3`$, $`\mathrm{}`$, is the sequence of grid points. Clearly, the above averaging procedure “cools” the fluctuations, in that $`\delta \varphi _a^2\delta \varphi _a^2`$. In particular, in the continuum limit $`\mathrm{\Delta }x0`$ with $`\xi =n\mathrm{\Delta }x`$ held constant, one of course finds that $`\delta \varphi _a^20`$ for the distribution (38). Therefore, in a final sweep over the grid one has to rescale the fields at each grid point by

$$\varphi _a^{\prime \prime }(\stackrel{}{x})=\varphi _a^{}(\stackrel{}{x})\sqrt{\frac{\delta \varphi _a^2}{\delta \varphi _a^2}}.$$
(41)
This procedure leads to an initial field configuration that exhibits both the proper physical correlation length and the desired fluctuations. Moreover, it prevents the initial energy density from field gradients to grow like $`1/\mathrm{\Delta }x^{\mathrm{\hspace{0.17em}2}}`$.

### B Numerical Solution

In the following we consider both a first order phase transition corresponding to $`g=5.5`$, as well as the critical point at $`g=3.7`$.

Figs. 2 and 3 depict the time evolution of the $`\sigma `$ field along the $`x`$-axis and the $`y`$-axis, respectively. At $`t=0`$ the field within the hot region has small amplitude, corresponding to the chiral symmetry restored phase. That region is surrounded by the physical vacuum with $`\sigma =f_\pi =93`$ MeV.

For the first order phase transition ($`g=5.5`$) a barrier separates the two degenerate minima of the effective potential (Fig. 1) at $`T_c`$. Figs. 2 and 3 show that this barrier leads to a rather well-defined surface in coordinate space, separating the vacuum from the symmetric phase. For the above-mentioned initial conditions, the fluctuations are not strong enough for the field to easily overcome the barrier. Nevertheless, one can observe dynamical fluctuations into the broken symmetry state, e.g. at $`x1`$ fm and $`t24`$ fm/c, which however collapse again. Our dynamical results agree with previous arguments that nucleation is a slow process on the time scale of heavy-ion collisions, and so the Gibbs phase equilibrium is not established dynamically . At time $`t9`$ fm the phase transition occurs spontaneously “in an instant”, that is, on a space-like surface which can clearly be seen in Figs. 2 and 3. A “bubble” of the symmetric phase survives at $`x3`$ fm for a long time.

The picture is rather different for the transition at the critical point, i.e. $`g=3.7`$. Here, the barrier between the degenerate minima vanishes and the potential is flat. As is evident, there are no clear surfaces separating either the vacuum from the center or high-density bubbles (or “droplets”) from their surrounding. Due to the flatness of the potential, near the center the field performs large-amplitude oscillations (from $`\sigma 0`$ to $`\sigma >f_\pi `$) for a long time; they extend in space over distances $`13`$ fm (e.g. at $`t9`$, 12, and 14 fm/c in Fig. 2), which is not much less than the initial size of the hot region. Figure 4 shows a histogram of the field distribution at the center ($`\sqrt{x^2+y^2+z^2}<2fm`$). One observes that the distribution broadens from $`t=4`$ fm to $`t=10`$ fm, and then narrows again at later times after the transition to the broken phase occured. Also, notice that the distributions at $`t=4`$ fm and $`t=14`$ fm are well described by Gaussians, i.e. the effective potential is essentially parabolic, while for $`t=10`$ fm there are visible deviations from a simple Gaussian.

The time evolution of the local rest-frame energy density $`e`$ of the quarks is shown in Figs. 56. Again we see large-scale structures for the first order phase transition ($`g=5.5`$), while the energy density is rather homogeneous on large time and distance scales if the expansion trajectory goes through the critical endpoint. For the first-order transition, quarks can be “trapped” in droplets with $`\sigma 0`$ (the minimum of the effective potential where the symmetry is restored) because the mass barrier can keep them from escaping. The droplet of high-density matter at $`x3`$ fm and $`t1216`$ fm/c can easily be associated to the region of nearly vanishing chiral scalar field from Fig. 2. Eventually, that region must perform the transition to the symmetry broken state, either by a strong thermal fluctuation or when reaching the spinodal point. At the spinodal, the system is as far from local thermal equilibrium as it can get, and the “roll-down” of the order parameter field to the global minimum of the potential can influence the collective expansion of the quark fluid.

Note that the energy density at the center drops more rapidly for the first-order transition than near the critical endpoint. This has consequences for the build-up of azimuthally asymmetric flow, as we shall discuss below.

Fig. 7 depicts the time evolution of the azimuthal momentum anisotropy

$$ϵ_p=\frac{T_{xx}T_{yy}}{T_{xx}+T_{yy}},$$
(42)
where the averages of the stress-energy tensor of the fluid are taken at fixed time:

$$T_{ij}(t)d^3xT_{ij}(t,\stackrel{}{x}).$$
(43)
Also, we average over a few initial field configurations, which gave similar results for $`ϵ_p(t)`$, though.

For the above-mentioned initial conditions the energy-momentum tensor is symmetric, and so $`ϵ_p=0`$ at $`t=0`$ (this might be different in more realistic treatments ). Pressure gradients in $`x`$ and $`y`$ directions are different, though. Therefore, the acceleration of the fluid is stronger in the reaction ($`xz`$) plane than out of plane, leading to a nonzero azimuthal asymmetry $`ϵ_p>0`$ at times $`t>0`$. The asymmetry first grows nearly linearly with time but saturates when the asymmetry of the energy density and of the pressure gradients becomes small. As explained above (Figs. 56), this happens earlier for a first-order transition than for trajectories near the chiral critical endpoint. This is then reflected in the final value of $`ϵ_p`$. We stress that the more rapid saturation of the azimuthal asymmetry in case of a first-order transition is not in contradiction to the fact that hot (high-energy density) “droplets” survive for rather long times, as seen in the figures. Rather, such “droplets” typically turn out to be more or less rotationally symmetric, or at least exhibit deformations which are uncorrelated to the reaction plane (the $`xz`$ plane in our case). Thus, they tend to reduce the average azimuthal asymmetry of the energy-momentum tensor.

For comparison, in Fig. 7 we also show the result for an equilibrium first-order phase transition. Here, the equations of motion for the chiral fields, eqs. (14) are not solved but rather the $`\sigma `$-field is required to populate the (global) minimum of the effective potential,

$`{\displaystyle \frac{\delta V_{\mathrm{eff}}}{\delta \sigma }}`$ $`=`$ $`0,`$ (44)
$`{\displaystyle \frac{\delta ^2V_{\mathrm{eff}}}{\delta \sigma ^2}}`$ $`>`$ $`0.`$ (45)
That is, the chiral field is in equilibrium with the quark-antiquark fluid and does not exhibit any explicit space-time dependence. At $`T_c`$, where two degenerate minima exist, one performs the usual Maxwell-Gibbs construction to determine the fractions of the total volume occupied by matter in the symmetric and the broken symmetry phases, respectively. Evidently, for the above-mentioned initial conditions the equilibrium phase transition leads to nearly the same azimuthal asymmetry of $`T^{ij}`$ as for the cross over. Therefore, it is indeed the non-equilibrium real-time dynamics (field fluctuations over the free-energy barrier) that is responsible for the observed reduction of $`ϵ_p`$ in the regime of first-order chiral phase transitions.

## IV Summary and Outlook

In summary, we have introduced a simple phenomenological description of the non-equilibrium real-time dynamics of the chiral phase transition in an expanding (relativistic) fluid of quarks. More precisely, we coupled the linear sigma model, which describes the dynamics of the long-wavelength modes of the chiral order parameter field, to the hydrodynamical evolution of a system of quarks. The chiral field(s) evolve according to the finite-temperature effective potential that is generated by integrating out the quarks from the Lagrangian; in turn, the field(s) determine the effective quark mass (i.e. the equation of state of the quark fluid) dynamically.

The above model exhibits a first order phase transition for large $`g`$, which is the quark-field coupling constant. The line of first order transitions ends in a critical point when $`g`$ is lowered, i.e. the transition turns into a cross over for smaller couplings. Thus, by varying $`g`$ one can qualitatively compare the hydrodynamic expansion pattern of the quark fluid for dynamical trajectories that cross the line of first order transitions to that obtained in the cross over regime.

We have obtained numerical solutions in 3+1 space-time dimensions, using simple initial conditions that might be appropriate for relativistic heavy-ion collisions. The hydrodynamical expansion pattern clearly depends on the structure of the effective potential. For trajectories in the cross over regime or near the critical endpoint the overall bulk dynamics is found to be rather “smooth”, in that the space-time distribution of the energy density of the fluid is not affected very much by the fluctuations of the order parameter field. In the absence of a latent heat, the energy density can not jump much between regions where the field amplitude is different.

In contrast, if the effective potential exhibits a barrier between the symmetry restored and broken phases, respectively, we do see that large-scale structures are formed dynamically, e.g. “droplets” of the symmetric phase may survive for rather long times before becoming mechanically unstable (at the spinodal). In that sense, the overall time scale is longer for trajectories crossing the line of first order transitions. Nevertheless, typically such structures are not correlated to the reaction plane; thus, the direct correspondence of spatial anisotropies in the initial condition to momentum-space anisotropies in the final state predicted by equilibrium hydrodynamics (that is, when the phase transition is not treated dynamically but modelled by a Maxwell-Gibbs construction) is weakened. For example, we find much smaller momentum-space anisotropy for a dynamical first-order transition than for a trajectory through the chiral critical endpoint (for the same initial condition). This could be a very useful prediction with regard to the experimental search for the chiral critical endpoint of QCD in heavy-ion collisions at the BNL-AGS, the CERN-SPS and the envisaged new GSI heavy-ion accelerator. Until now, experiments focused on fluctuation observables, but inclusive observables usually are much easier to analyze accurately.

In the future, we intend to scrutinize other inclusive observables as to their sensitivity to non-equilibrium effects from phase transitions. Of course, there is also plenty of room to improve on the model in order to obtain more quantitative predictions. The present paper represents a first step towards an actual real-time description of the chiral phase transition on either side of the critical endpoint in expanding relativistic fluids with realistic 3+1d geometries.
Note added: After this manuscript was submitted for publication the NA49 collaboration published the elliptic flow at $`E_{\mathrm{lab}}=40A`$ GeV . From Fig. 24 of that publication, the dependence of $`v_2`$ on $`\mathrm{log}\sqrt{s}`$ is approximately linear. However, the ”natural” scale for $`v_2`$ is set by $`p_\mathrm{t}`$, not $`\mathrm{log}\sqrt{s}`$, as pointed out by Snellings . Indeed, at high energies the differential $`v_2(p_\mathrm{t})`$ of charged hadrons is approximately proportional to $`p_\mathrm{t}`$, such that the averaged $`v_2p_\mathrm{t}`$. In fact, for mid-central collisions $`v_2`$ increases from $`3\%`$ at top SPS energy ($`\sqrt{s}=18A`$ GeV) to $`4.5\%`$ at RHIC energy ($`\sqrt{s}=130A`$ GeV). When scaled by the average transverse momentum, though, the elliptic flow in that energy regime is nearly constant . To scrutinze deviations from the ”natural” scaling $`v_2p_\mathrm{t}`$, we plot the excitation function of $`v_2/p_\mathrm{t}`$ in Fig. 8.

One observes that, as already mentioned above, the data is compatible with no energy (or $`p_\mathrm{t}`$) dependence above top SPS energy. Clearly, there is a systematic drop of $`v_2`$ relative to $`p_\mathrm{t}`$ towards lower energies. For instance, at $`E_{\mathrm{lab}}=40A`$ GeV corresponding to $`p_\mathrm{t}350`$ MeV, $`v_2/p_\mathrm{t}`$ is lower by about two standard deviations than at higher energies. Qualitatively, this interesting behavior is similar to the reduction of the azimuthal momentum asymmetry, predicted above, caused by crossing the second order critical point into the regime of first order phase transitions. Additional studies of “conventional” non-equilibrium effects unrelated to a phase transition are certainly required, however, before firm conclusions can be drawn.

###### Acknowledgements.

We thank T. Kodama, I. Mishustin, O. Scavenius and R. Snellings for helpful discussions, and D. Rischke for a critical reading of the manuscript as well as for communicating his 3+1d SHASTA code. K.P. and H.S. gratefully acknowledge support by GSI, DFG, DESY, the Bergen Computational Physics Laboratory (Project No. 43) and the Frankfurt Center for Scientific Computing. A.D. was partly supported by the U.S. Department of Energy under contract No. DE-AC02-98CH10886 and by the German Ministry for Education and Research (BMBF).

## A Subtracting Initial Fluctuations

In section II A we discussed our phenomenological ansatz for the effective potential for the long wavelength modes of the chiral fields, as generated by the heat bath of quarks. Formally, it is obtained from the Lagrangian by integrating out the quarks to one loop.

Our main objective here is to study dynamically fluctuations of the chiral order parameter (or effects generated by those fluctuations) in the vicinity of the chiral critical point as the system makes the transition to broken chiral symmetry at low temperature. Thus, we have to allow for primordial fluctuations of the chiral fields, also. However, those fluctuations of the fields at time $`t=0`$ will of course also contribute to the effective potential and “distort” its shape. In order to restore our original choice for $`V_{\mathrm{eff}}`$ from eq. (12), and thus ensure the correct dynamics for the long-wavelength modes, we have to introduce appropriate subtractions.

The procedure is as follows. The scalar and pseudoscalar densities are given by eqs. (16). They depend explicitly on the value of the fields $`\varphi _a`$, which we formally separate into short and long wavelengths:

$$\varphi _a(x)=\varphi _a+\delta \varphi _a(x).$$
(A1)
Here, $``$ denotes a spatial average over a volume large enough for the fluctuations to average out:

$$\delta \varphi _a=0.$$
(A2)
The linear dimension of that volume will be given by the wavelength of the soft modes of interest.

We now substitute (A1) into eqs. (16) and perform an expansion up to second order in $`\delta \varphi _a(x)`$. We then perform the averaging over the fluctuations with the distribution (38) and obtain

$`\rho _s`$ $`=`$ $`\rho _s(\varphi )`$ (A3)
$`+`$ $`{\displaystyle \frac{1}{2}}g\sigma d_q{\displaystyle }{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}}\{.`$ (A8)
$`{\displaystyle \underset{a}{}}\delta \varphi _a^2[{\displaystyle \frac{g^2f(k)}{E^2T}}({\displaystyle \frac{T}{E}}+f(k)e^{E/T})`$
$`+{\displaystyle \frac{g^4\varphi _a^2}{E^3T^2}}f(k)(2f^2(k)e^{2E/T}f(k)e^{E/T}`$
$`+3{\displaystyle \frac{T}{E}}f(k)e^{E/T}+3{\displaystyle \frac{T^2}{E^2}})]`$
$`2\delta \sigma ^2{\displaystyle \frac{g^2f(k)}{E^2T}}({\displaystyle \frac{T}{E}}+f(k)e^{E/T})\}`$
Here, $`\delta \varphi ^2_a\delta \varphi _a\delta \varphi _a`$ is the variance of the fluctuations in the initial condition, i.e. at the initial time $`t=0`$, summed over internal quantum numbers. We made use of the fact that the fluctuations (38) are diagonal in internal space, i.e. $`\delta \varphi _a\delta \varphi _b=0`$ if $`ab`$. The second term is the additional contribution seen by the long wavelength modes $`\varphi _a`$, which is due to the fluctuations. To restore the original effective potential, we have to subtract that term, i.e. redefine the scalar density as

$`\rho _s(\varphi ,T)`$ $`=`$ $`g\sigma d_q{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}\frac{1}{E}f(k)}`$ (A9)
$``$ $`{\displaystyle \frac{1}{2}}g\sigma d_q{\displaystyle }{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}}\{.`$ (A14)
$`{\displaystyle \underset{a}{}}\delta \varphi _a^2[{\displaystyle \frac{g^2f(k)}{E^2T}}({\displaystyle \frac{T}{E}}+f(k)e^{E/T})`$
$`+{\displaystyle \frac{g^4\varphi _a^2}{E^3T^2}}f(k)(2f^2(k)e^{2E/T}f(k)e^{E/T}`$
$`+3{\displaystyle \frac{T}{E}}f(k)e^{E/T}+3{\displaystyle \frac{T^2}{E^2}})]`$
$`2\delta \sigma ^2{\displaystyle \frac{g^2f(k)}{E^2T}}({\displaystyle \frac{T}{E}}+f(k)e^{E/T})\}`$
This expression has to be substituted for $`\rho _s`$ on the right-hand-side of the equation of motion (14). The subtracted term cancels the “distortion” of the scalar density caused by using the local values $`\sigma (x)`$, $`\stackrel{}{\pi }(x)`$ for the fields, rather than their long-wavelength components $`\sigma (x)`$, $`\stackrel{}{\pi }(x)`$.

Along the same lines one derives the following expressions for the fluctuation-subtracted pseudo-scalar density, and for the pressure of the quarks:

$`\stackrel{}{\rho }_{\mathrm{ps}}(\varphi ,T)`$ $`=`$ $`g\stackrel{}{\pi }d_q{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}\frac{1}{E}f(k)}`$ (A15)
$``$ $`{\displaystyle \frac{1}{2}}g\stackrel{}{\pi }d_q{\displaystyle }{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}}\{.`$ (A20)
$`{\displaystyle \underset{a}{}}\delta \varphi _a^2[{\displaystyle \frac{g^2f(k)}{E^2T}}({\displaystyle \frac{T}{E}}+f(k)e^{E/T})`$
$`+{\displaystyle \frac{g^4\varphi _a^2}{E^3T^2}}f(k)(2f^2(k)e^{2E/T}f(k)e^{E/T}`$
$`+3{\displaystyle \frac{T}{E}}f(k)e^{E/T}+3{\displaystyle \frac{T^2}{E^2}})]`$
$`2\delta \stackrel{}{\pi }^2{\displaystyle \frac{g^2f(k)}{E^2T}}({\displaystyle \frac{T}{E}}+f(k)e^{E/T})\}`$
$`p(\varphi ,T)`$ $`=`$ $`d_q{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}T\mathrm{log}\left(1+e^{E/T}\right)}`$ (A21)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{a}{}}\delta \varphi _a^2\{{\displaystyle \frac{g^2f(k)}{E}}`$ (A22)
$`+`$ $`{\displaystyle \frac{g^4\varphi _a^2}{E^3T}}f(k)[T+EEf(k)]\}.`$ (A23)
One can verify that the identities $`g\rho _s=\delta p/\delta \sigma `$ and $`g\rho _{ps}=\delta p/\delta \stackrel{}{\pi }`$ are satisfied, as it should be.

At fixed values for the fields, the energy density of the quarks at a temperature $`T`$ is given by

$$e(\varphi ,T)=T\frac{p(\varphi ,T)}{T}p(\varphi ,T).$$
(A24)
Using the expression for the fluctuation-subtracted pressure given above one obtains

$`e(\varphi ,T)`$ $`=`$ $`d_q{\displaystyle \frac{d^3k}{(2\pi )^3}Ef(k)}`$ (A25)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{a}{}}\delta \varphi _a^2\{{\displaystyle \frac{g^2f(k)}{ET}}[TEf(k)e^{E/T}]`$ (A26)
$``$ $`{\displaystyle \frac{g^4\varphi _a^2}{E^3T}}f(k)\left[T+EEf(k)\right]`$ (A27)
$`+`$ $`{\displaystyle \frac{g^4\varphi _a^2}{ET^2}}f^2(k)e^{E/T}\mathrm{tanh}\left({\displaystyle \frac{E}{2T}}\right)\}`$ (A28)
The source term (27) changes due to the fluctuations and one has to use the modified scalar (A9) and pseudoscalar (A15) densities, respectively.

To second order in the fluctuations, the self-interaction of the chiral fields is renormalized as

$`U(\varphi _a)`$ $`=`$ $`{\displaystyle \frac{\lambda ^2}{4}}\left(\sigma ^2+\stackrel{}{\pi }^2v^2\right)^2h_q\sigma U_0`$ (A29)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{a}{}}\delta \varphi _a^2\lambda ^2\left(2\varphi _a^2+\sigma ^2+\stackrel{}{\pi }^2v^2\right)`$ (A30)
The above expressions for $`U(\varphi _a)`$, $`\rho _s`$, $`\stackrel{}{\rho }_{ps}`$, $`e`$, and $`p`$ are to be used in the equations of motion for the chiral fields (14), in the stress-energy tensor of the quark fluid (20), and in the source term for its divergence $`S^\nu `$ (27,31). We point out that we subtract those quantities for the contribution from initial fluctuations of $`\varphi `$ only up to second order in $`\delta \varphi `$. We can therefore not employ initial conditions with very large local fluctuations about the mean field.