A Equation of state

# Spinodal phase separation in relativistic nuclear collisions

## Abstract

The spinodal amplification of density fluctuations is treated perturbatively within dissipative fluid dynamics for the purpose of elucidating the prospects for this mechanism to cause a phase separation to occur during a relativistic nuclear collision. The present study includes not only viscosity but also heat conduction (whose effect on the growth rates is of comparable magnitude but opposite), as well as a gradient term in the local pressure, and the corresponding dispersion relation for collective modes in bulk matter is derived from relativistic fluid dynamics. A suitable two-phase equation of state is obtained by interpolation between a hadronic gas and a quark-gluon plasma, while the transport coefficients are approximated by simple parametrizations that are suitable at any degree of net baryon density. We calculate the degree of spinodal amplification occurring along specific dynamical phase trajectories characteristic of nuclear collision at various energies. The results bring out the important fact that the prospects for spinodal phase separation to occur can be greatly enhanced by careful tuning of the collision energy to ensure that the thermodynamic conditions associated with the maximum compression lie inside the region of spinodal instability.

###### pacs:
25.75.-q, 81.30.Dz, 64.75.Gh, 64.60.an

## I Introduction

It is expected that the confined and deconfined phases of strongly interacting matter may coexist at net baryon densities above a certain critical value and significant experimental efforts are underway to search for evidence of the associated first-order phase transition and its critical end point: a systematic beam-energy scan is currently being performed at RHIC (BNL) to look for the critical point (1); the CBM experiment at FAIR (GSI) will study baryon-dense matter and search for the phase transition (2); and the proposed NICA (JINR) aims at exploring the mixed phase (3).

These studies are rather challenging, not only because it is inherently difficult to extract the properties of equilibrium bulk matter from collision experiments, but also because there is currently no suitable transport model available for guiding these efforts. As a result, there is currently an urgent need for identifying experimentally observable signals of the phase structure.

The present paper focusses on the possibility that the mechanism of spinodal phase decomposition may have effects that could be exploited as signals of the phase transition. Spinodal decomposition is a well-known generic phenomenon associated with first-order phase transitions that has been studied in a variety of substances and also found industrial application (4). Furthermore, nuclear spinodal fragmentation (5) was observed in nuclear collisions at intermediate energies (6) several years ago. A preliminary study (7) found grounds for guarded optimism that spinodal separation between the confined and deconfined phases could in fact occur in relativistic collisions and we have therefore undertaken the present more refined analysis.

While that earlier study (7) employed a somewhat schematic equation of state based on a generalized classical gas, the present uses a more realistic equation of state obtained by interpolating between a hadron gas and a quark-gluon plasma. An advantage of this procedure is that it automatically incorporates the increase in the number of degrees of freedom in the dense (deconfined) phase, a peculiar but important characteristic of strongly interacting matter. Building on the developments in Ref. (7), we take account of finite-range effects by including a gradient term in the equation of state. This refinement is essential for obtaining a physically meaningful description since it ensures both that there is an interface tension between the two coexisting phases and that the spinodal growth rates subside at large wave numbers.

We again carry out our studies within the framework of fluid dynamics, because this type of transport description has the distinct advantage that the complicated and still poorly understood microstructure of the system enters only via the equation of state and the transport coefficients. A general discussion of the fluid-dynamical description of  first-order phase transitions was given recently in Ref. (8).

Although the dispersion equation in (7) was derived with both shear and bulk viscosity included, the actual calculations were done for ideal fluid dynamics. In the present study, the dynamical calculations include not only viscosity but also heat conductivity, which proves to be as important as viscosity while affecting the spinodal growth rates oppositely, as we demonstrate quantitatively. The associated cubic dispersion equation is derived directly from relativistic fluid dynamics. The medium dependence of the transport coefficients is expressed in terms of the local values of the enthalpy density and the particle spacing, an approximation that applies not only in the baryon-poor regime but also in the baryon-dense media of relevance to the phase transition. The strength of the transport coefficients for arbitrary density and temperature can thus be related to the values obtained from analyses of the RHIC data.

We seek to construct plausible dynamical phase trajectories by invoking results from earlier calculations with various transport models and we examine in particular the crucial importance of using a collision energy for which the maximum compression occurs inside the spinodal phase region. Once the phase trajectory of the collision system has been specified, we may integrate the spinodal growth rate along the dynamical history and thus calculate the resulting degree of amplification as a function of the wave number of the density perturbation. We do that for a range of dissipation strengths that bracket those expected from the analysis of the RHIC data.

The present, more refined, studies suggest that spinodal phase decomposition may indeed be triggered during nuclear collisions within a certain (likely relatively narrow) optimal energy range. This expectation is rather insensitive to the still poorly known strength of the transport coefficients. Such a spinodal phase separation would result in an assembly of plasma drops embedded in a hadron gas and our present analysis permits us to estimate the typical drop size.

The presentation is organized as follows. We first discuss the expected thermodynamic phase structure of strongly interacting matter within the framework of the specific equation of state that we have constructed. We then turn to dissipative fluid dynamics within which we derive the dispersion equation for the collective modes in bulk matter. Subsequently we develop expressions for the transport coefficients in baryon-rich matter and use those in calculations of the spinodal growth rates. Finally, we obtain quantitative results for the degree of spinodal amplification experienced by the bulk of the collision system as it evolves along various plausible phase trajectories. The construction of the equation of state and the associated spline procedure are described in appendices.

## Ii Thermodynamic phase structure

In order to make quantitative studies, we need to employ a specific equation of state that is plausibly realistic and, in particular, exhibits the expected phase structure.

Although significant progress has been made in understanding the thermodynamic properties of the confined and deconfined phases separately, our current understanding of the phase coexistence region is not yet on firm ground. We therefore employ a conceptually simple approximate equation of state in which the region of phase transformation is described by means of a suitable interpolation between an idealized hadron gas and an idealized quark-gluon plasma. The details of this construction are described in Appendix A, while the resulting phase structure is shown in Fig. 1.

For the present study, it is convenient to work in the canonical representation where the thermodynamic state of the system is characterized by the temperature and the (net) baryon density . (In the phase region of primary interest, the temperature is relatively low and the chemical potential relatively high, so the equilibrium population of antibaryons is relatively insignificant, .) The key thermodynamic function is then the free energy density,s , from which the other quantities may be obtained,

 Chemical potential: μT(ρ) = ∂ρfT(ρ) , (1) Pressure: pT(ρ) = ρ∂ρfT(ρ)−fT(ρ) , (2) Entropy density: σT(ρ) = −∂TfT(ρ) , (3) Energy density: εT(ρ) = fT(ρ)−T∂TfT(ρ) . (4)

We may also express the isothermal sound speed ,

 v2T ≡ ρhT(∂p∂ρ)T =ρhT∂ρpT(ρ) =\ ρ2hT∂2ρfT(ρ) . (5)

where is the enthalpy density, as well as the isentropic sound speed ,

 v2s ≡ ρhT(∂p∂ρ)s =\ v2T+ThT1cv(∂TpT(ρ))2 ≥ v2T , (6)

where the specific heat is .

Different manifestations of the system that have the same values of temperature, chemical potential, and pressure are in mutual thermodynamic equilibrium and may thus coexist. As Fig. 1 shows, such phase coexistence occurs for temperatures below the critical value,  MeV. The corresponding density values, and , are traced out (outer solid curve). Inside this boundary, it is thermodynamically favorable for uniform matter to separate into the two coexisting phases. However, in the part of the phase coexistence region that is close to the boundary any small deviation from uniformity is thermodynamically unfavorable and here the system is mechanically metastable, as signalled by the fact that in this phase region.

This situation changes radically at the spinodal boundary where the isothermal sound speed vanishes, (the dashed curve in Fig. 1): Inside this boundary even infinitesimal deviations from uniformity are thermodynamcially favorable and bulk matter is thus mechanically unstable. As a consequence, small density fluctuations may be amplified and thereby cause the system to undergo a spontaneous phase separation. The phase region of spinodal instability extends in temperature all the way from zero up to .

Finally, Fig. 1 also shows the boundary where the isentropic sound speed vanishes, . Inside this smaller phase region, the spinodal amplification process can occur without any entropy production. This special region does not extend all the way up to and it narrows considerably as the temperature approaches its upper bound.

The equation of state applies to idealized uniform matter. However, since we are interested in following the evolution of density disturbances, it is essential to include finite-range effects. Indeed, the physical coexistence between two different phases along a common interface could not be realized without taking proper account of the density gradients, nor could the associated interface tension be obtained. We shall therefore augment the bulk thermodynamics with a gradient term as proposed in Ref. (7) (see Sect. III.3.1), which then extends the validity of the equation of state to non-uniform systems. In particular, it is possible to describe the diffuse interface between two coexisting phases and the associated tension. Furthermore, as we shall see, the gradient term is essential for obtaining a physically reasonable dispersion relation because it stabilizes disturbances having a small spatial scale.

## Iii Dissipative fluid dynamics

We wish to employ dissipative fluid dynamics for our dynamical studies. Fluid dynamics is convenient because the specific miscroscopic structure of the matter under consideration enters only via the equation of state and a few transport coefficients. On the other hand, the treatment relies on the assumption of approximate local equilibrium which may generally be questionable in nuclear collisions. Fortunately, our applications are to collisions of relatively modest energies and, moreover, to the later stages in the evolution. Thus the conditions for applicability should be reasonably favorable.

We base our treatment on the relativistic formulation by Muronga (9). The four-velocity of the local flow is and the symmetric tensor projects onto the 3-space orthogonal to , .
The space-time derivative decomposes, , where the convective time derivative is while the gradient is .

The equations of motion simplify when expressed in a specific reference frame. For the present study, it is convenient to use the Eckart frame, which is defined in terms of the charge flow and is usually employed in non-relativistic scenarios. Then, since the local charge flow, , vanishes by definition, the charge four-current density is , where is the charge density in the local flow frame. However, it may generally be preferable to use the Landau frame because it is defined in terms of the energy flow and is thus meaningful also for chargeless fluids and fluids with several conserved charges for which there is no unique generalization of the Eckart frame.

### iii.1 Small disturbances

We consider the early evolution of small deviations from uniformity. We generally assume that these are planar and harmonic, with , and similarly for the other quantities. We may then ignore terms of order and higher. It follows that the associated flow velocities are small, since , and thus we have .

Generally the energy flow is given by , where is the heat flow and the enthalpy density. Since the local charge flow vanishes when we use the Eckart frame (see above), the energy flow equals the heat flow, , and is . Hence , and the energy-momentum tensor simplifies to

 Tμν = εuμuν−pΔμν+πμν−ΠΔμν . (7)

Here is the energy density in the local flow frame and is the sum of the local isotropic pressure and the pressure induced by the bulk viscosity which enters through the bulk pressure,

 Π = −ζ∇μuμ ≈ −ζ∇ivi = −ζ∂ivi = −ζ\boldmath∇⋅\boldmathv . (8)

Furthermore, the heat flow is , while the shear viscosity enters via the stress tensor

 πμν = η[ΔμσΔντ+ΔμτΔνσ−23ΔμνΔστ]∇σuτ (9) ≈ η[ΔμiΔνj+ΔμjΔνi−23ΔμνΔij]∇ivj , (10)

where we have used that only the spatial components of contribute to leading order in and, moreover, any derivatives of can be ignored so only the spatial components contribute. Furthermore, since , which is , only the spatial elements contribute. Hence only the spatial part is non-vanishing. It has the following elements,

 πij ≈ −η[∂ivj+∂jvi−23δij∂kvk] . (11)

Thus, for small deviations from uniformity, the spatial part of the energy-momentum tensor is given by

 Tij ≈ δijp−η[∂ivj+∂jvi−23δij∂kvk]−ζδij\boldmath∇⋅\boldmathv . (12)

If the spatial variation of the viscosity coefficients (shear) and (bulk) may be ignored, we have

 \boldmath∇⋅\boldmathT ≈ \boldmath∇p−η\boldmathΔ\boldmathv−[13η+ζ]\boldmath∇(\boldmath∇⋅% \boldmathv) , (13)

The gradient simplifies further in a semi-infinite geometry (where the only spatial variation is in the direction),

 \boldmath∇⋅\boldmathT ≍ ∂xTxx ≈ ∂xp−[\makebox$43$η+ζ]∂2xv . (14)

Thus only the effective viscosity enters. It follows that a uniform stretching, , is dissipation free. We also wish to point out that the shear viscosity contributes even though the flow has no shear.

It is interesting to note that the above result reflects a general feature of isotropic expansions in dimensions. To see this, assume that and . The viscous term in the Euler equation may then be evaluated by use of spherical coordinates,

 η\boldmathΔ\boldmathv+[13η+ζ]\boldmath∇(\boldmath∇⋅\boldmathv) = ξ^\boldmathr∂r1rN−1∂rrN−1v . (15)

Thus it is always the combination that enters. Furthermore, it follows that a Hubble-type expansion, , is dissipation free in any dimension.

### iii.2 Equations of motion

The fluid-dynamic equations of motion reflect the conservation of (baryon) charge, momentum, and energy. We are interested in the dynamics of small deviations from uniformity in a semi-infinite configuration and we focus on harmonic disturbances,

 ρ(\boldmathr,t) = ρ0+δρ(x,t) ≐ ρ0+ρkeikx−iωt , (16) ε(\boldmathr,t) = ε0+δε(x,t) ≐ ε0+εkeikx−iωt , (17) p(\boldmathr,t) = p0+δp(x,t) ≐ p0+pkeikx−iωt , (18)

and similarly for the other dynamical variables.

The conservation of charge is ensured by the continuity equation, , which here becomes

 C: ∂tρ ≐ −ρ0∂xv ⇒ ωρk ≐ ρ0kvk . (19)

It serves to eliminate the flow velocity, . The momentum equation simplifies considerably for the present scenario of small disturbances,

 M: h0∂tv ≐ −∂x[p−ζ∂xv]−∂xπxx−∂tq , (20)

where is the enthalpy density of the uniform system and the heat flow is (see below). The equation for energy conservation is similarly simplified,

 E: ∂tε ≐ −h0∂xv−∂xq . (21)

By combining these latter two equations, (20) and (21), one obtains the sound equation,

 ∂tE−∂xM: ∂2tε ≐ ∂2xΔ[p−ζ∂xv]+∂2xπxx , (22)

which amounts to where we recall that , see Eq. (14).

### iii.3 Dispersion equation

When heat conductivity is ignored (), the energy density tracks the charge density, as follows immediately from the energy equation, . Furthermore, in the absence of a gradient term in the equation of state (see later), we have with and where is the microcanonical equation of state. Since the isentropic sound speed is given by , we obtain the familiar viscous dispersion equation, (7).

To obtain the dispersion equation with heat conductivity included, we must invoke the form of the heat current,

 q≈−κ[∂xT+T0∂tv]: qk=−iκ[kTk−T0ρ0ω2kρk] . (23)

Insertion of this expression into the energy equation (21) yields a relationship between , , and ,

 εk = h0ρ0ρk+kωqk ≈ h0ρ0ρk−iκk2ωTk , (24)

where the term has been ignored because if is in comparison with . Thermodynamics enables us to express in terms of and ,

 εk=(∂ε∂T)ρTk+(∂ε∂ρ)Tρk =\ −cvTk−σερσεερk , (25)

where and are second derivatives of the entropy density . We have also used , the heat capacity at constant density. Using furthermore

 Missing or unrecognized delimiter for \left (26)

we may then obtain from the energy equation (21),

 Tk ≈ 11+iκk2/ωcv T0ρ0(∂p∂ε)ρρk . (27)

The canonical equation of state allows us to express the pressure variation in terms of the variations in temperature and density,

 pk=(∂p∂T)ρTk+(∂p∂ρ)Tρk=(∂p∂ε)ρcvTk+h0ρ0v2Tρk, (28)

where is the isothermal sound speed, see Eq. (5).

With these preparations, the dispersion equation can then be obtained by substituting the relations (25), (27), (28) into the sound equation (22) and using the relationship ,

 ω2 ≐ v2Tk2−iξωh0k2+v2s−v2T1+iκk2/ωcvk2 , (29)

retaining heat conduction terms only up to . This is recognized as the dispersion equation given in Ref. (10).

As noted above (end of Sect. II), it is essential to take account of finite-range effects, without which the spinodal growth rate would become ever larger as the wave number is increased (11). Following Ref. (7), we introduce a gradient correction in the equation of state. To leading order in the disturbance amplitudes, the effect of the gradient term on the local pressure is given by

 p(\boldmathr) ≈ p0(ε(\boldmathr),ρ(\boldmathr))−Cρ0∇2ρ(\boldmathr) , (30)

where is the microcanonical equation of state, i.e. the pressure in uniform matter having the specified energy and charge densities. The pressure amplitude is then modified accordingly,

 pk ⇝ pk+Cρ0k2ρk . (31)

Hence we should augment the term in Eq. (28),

 h0ρ0v2Tρk ⇝ [h0ρ0v2T+Cρ0k2]ρk . (32)

The full dispersion equation is then

 ω2 ≐ v2Tk2+Cρ20h0k4−iξωh0k2+v2s−v2T1+iκk2/ωcvk2 , (33)

i.e. the gradient term is simply added, just as when there is no heat conductivity (7).

#### Solution of the dispersion equation

When heat conductivity is included, , the dispersion equation is of third order and, consequently, there are three eigenmodes for each wave vector , as one would generally expect since the energy density is now no longer tied to the baryon density . This equation has one purely imaginary solution, , and a pair of generally complex solutions, , which are either both also imaginary, , or have the form . In the latter case it is easy to see that in the normal region where .

The two frequencies obtained for ideal (i.e. non-dissipative) fluid dynamics are either purely real (outside the isentropic spinodal phase region where ) or purely imaginary (inside the isentropic spinodal region where ). Thus, in ideal fluid dynamics, the region of spinodal instability is bounded by the isentropic spinodal, . The introduction of viscosity adds a negative imaginary amount to the frequency. We then have to first order in , where we have introduced the characteristic viscous length . But the inclusion of viscosity does not change the region of instability. In contrast, the inclusion of heat conductivity expands the region of instability from the isentropic spinodal to the isothermal spinodal, , and generally increass the spinodal growth rates.

For small distortions of uniform matter at given density and temperature, the dispersion relation yields the eigenfrequencies in terms of the equation of state , the strength of the gradient term , and the transport coefficients , , and . We describe below what we adopt for these key functions.

## Iv Transport coefficients

The deviation of the dynamical evolution from that of an ideal fluid is governed by three transport coefficients: the shear viscosity and the bulk viscosity (which here enter only through the effective viscosity ) as well as the heat conductivity . Neither their magnitudes nor their dependencies on the environment (through and ) are very well known. We shall therefore employ simple parametrizations of their functional form and introduce one adjustable overall strength parameter for each one, thus enabling us to conveniently explore a range of physical scenarios.

### iv.1 Viscosity

Using string theory methods, Kovtun et al. (12) made general arguments that the shear viscosity , in any relativistic quantum field theory at finite temperature and zero chemical potential, has a lower bound, , where is the associated entropy density.

It might therefore seem natural to use , i.e. simply scale the minimum value by the factor . However, the entropy density , while appropriate in the context of ultrarelativistic nuclear collisions where the medium has a high temperature and a very small net baryon density, is not suitable in the present context where the focus is on matter at high baryon density and relatively modest temperature (13). A more appropriate quantity, suitable in both limits, is the enthalpy density (13). When the net baron density vanishes it becomes , whereas approaches the mass density at high baryon density and low temperature, , where is density of particles and is their mass. (For cold nuclear matter we have , since the pions and anti-nucleons have negligible populations, and hence .)

Furthermore, one would expect the viscosity to be proportional to the interparticle spacing (13), which provides a convenient measure of the mean free path in a dense fluid. When plasma has no net baryon density, is inversely proportional to the temperature , . The conversion constant is given by

 c0 ≡ 14π[(gg+\makebox$32$gq)ζ(3)π2]13 ≈ 0.12779 . (34)

Therefore, based on these considerations, we shall make the following ansatz,

 η(ρ,T) ≐ η0c0cd(ρ,T)h(ρ,T) , η0≥1 , (35)

which amounts to in the baryon-free plasma, while it becomes in the non-relativistic gas. This latter expression is similar to the familiar expression from classical transport theory for the shear viscosity coefficient in a dilute one-component gas, , where is the mean particle speed and is its mean free path.

Since the bulk viscosity is generally expected to be significantly smaller than the shear viscosity , we shall take the effective viscosity to be . It is convenient to introduce the associated characteristic length which is proportional to the interparticle spacing,

 λvisc(ρ,T) ≡\ 1cξ(ρ,T)h(ρ,T)/c2 =\ \makebox$43$η0c0d(ρ,T) . (36)

### iv.2 Heat conductivity

The thermal conductivity is fundamentally related to the viscosity because it derives from the same microscopic transport processes. In a dilute classical gas it can be expressed as , where is the mean particle speed and is the heat capacity (equal to for a dilute classical gas). Since and , we see that . We therefore make the following ansatz,

 κ(ρ,T) ≐ κ0c0cd(ρ,T)cv(ρ,T) , κ0≥1 , (37)

with as given in Eq. (34) and assuming that the overall strength factor should be at least unity. Although the relation would suggest that the two normalization constants should be similar, , we prefer to leave them separately adjustable in order to make it possible to explore the effects of the two distinct types of dissipation. The characteristic length scale associated with heat conduction is then

 λheat(ρ,T) ≡\ 1cκ(ρ,T)cv(ρ,T) =\ κ0c0d(ρ,T) . (38)

While the approximate expressions for the transport coefficients, Eqs. (35) and (37), should not be expected to be accurate, they will serve well for our present purpose of exploring the effect of the dissipative mechanisms on the spinodal decomposition, since they enable us to easily control the overall dissipative effects.

### iv.3 Growth rates

The spinodal isothermal and isentropic boundaries determined by and , respectively, pertain to the thermodynamic limit of very long wave lengths, . When the wave number is increased, the region of instability steadily shrinks as the gradient term gives an ever larger contribution to the pressure. Thus, at a given temperature , the lower spinodal boundary density increases steadily with , while the upper spinodal boundary density decreases steadily. For a fixed value of , a contour plot of the growth rate in the phase plane therefore exhibits a ridge between those two boundaries. The height of of the ridge decreases steadily as is increased, first rather gently due to the dominance of the fermions at low temperatures. The local value of the maximum wave number for which instabilities exist, , will have a similar appearance.

We now illustrate, in Fig. 2, the resulting dispersion relations for thermodynamic scenarios relevant to the present study. Selecting a phase point in the central region of the phase coexistence region where both the isothermal and the isentropic sound velocities are imaginary, and (see Fig. 1), we consider the growth rate as a function of the wave number of the density undulation being amplified. The non-dissipative treatment with ideal finite-range fluid dynamics provides a convenient reference result. It yields a fastest growth time of about which occurs for wave numbers near , corresponding to an optimal wave length of .

Relative to this reference, the inclusion of viscosity slows the growth but does not change the domain of instability which is still delineated by the vanishing of the isentropic sound speed . We see that the inclusion of a minimal amount of viscosity () leads to a significant reduction in and also shifts the optimal length scale towards larger values.

On the other hand, relative to the ideal scenario, the inclusion of heat conductivity enlarges the domain of instability, the boundary being now determined by the vanishing of the isothermal sound speed . Thus, generally, the inclusion of heat conductivity increases the growth rates, particularly at the high end of the unstable range.

While the inclusion of both minimal viscosity and minimal heat conduction necessarily enlarges the unstable range, it does somewhat reduce the fastest growth rates. However, it hardly affects the scale of the fastest-growing modes, . As the strengths of the dissipative terms are further increased, the growth rate decreases steadily and, at the same time, the maximum in the dispersion relation moves gradually downwards in .

These features are present throughout the unstable region of the phase plane. In order to obtain an impression of how the growth rates change with and , we show in Fig. 3 the fastest growth rate as we move away from the phase point explored in Fig. 2, either at constant (upper panel) or constant (lower panel), as indicated in Fig. 1. In these illustrations we show the ideal scenario where no dissipation is present (), the scenario with minimal dissipation (), and a scenario with five times stronger dissipation (), thus covering the range suggested by the RHIC data (14); (15); (16); the currently favored estimate is . Since, as noted above, the heat conductivity is fundamentally related to the shear viscosity, it should be expected that their strengths vary in approximate unison; the effect of varying their relative strengths may be judged from the results shown in Fig. 2.

As expected from the discussion in the beginning of this section, the fastest growth at a given temperature occurs about midway between the corresponding spinodal boundaries. Furthermore, the temperature generally reduces the growth rates, though only weakly at small . The curves in the lower panel of Fig. 3 do not exhibit a monotonic decrease because the selected path along a constant density does not follow the ridgeline: the maxima in the curves shown in the upper panel shift steadily towards smaller densities as is increased, as would be expected from the general appearance of the phase diagram.

Finally, we note that a five-fold increase in the dissipation above the minimal value, i.e. changing from to , leads to a reduction in the growth rates by only about a factor of two.

## V Dynamical evolution

After the above preparations, we are now in a position to address the dynamical evolution of the unstable collective modes in the spinodal region of the phase diagram.

### v.1 Dynamical phase trajectories

We first specify dynamical phase trajectories, , that are representative of the bulk matter in the collision zone, making use of the results presented in Ref. (17). In that work, a number of different dynamical models were used to extract the time evolution of the net baryon density, , and the energy density, , in the center of a head-on gold-gold collision for the range of collision energies anticipated at FAIR. The resulting dynamical trajectories in the phase plane were remarkably independent of the specific model, in large part probably due to the robust nature of the mechanical densities which are subject to conservation laws, in contrast to the corresponding thermodynamical variables . Nevertheless, there were significant variations in the detailed behavior. This is illustrated in Fig. 4 which shows density evolutions obtained with the 3-fluid model (18) and UrQMD (19) at beam kinetic energies of and . (Note that these beams are bombarded onto stationary targets, as will be done at FAIR; the same collision energies can be obtained in a collider configuration by using total beam energies of and , respectively.) We utilize these results to construct the dynamical phase trajectories considered.

As expected, the degree of amplification achieved depends strongly on the length of time spent in the phase region of spinodal instability. In order to elucidate this key feature, we consider in some detail the phase trajectories depicted on Fig. 5. The density evolutions are those calculated with the two models for , while the temperature evolutions are (somewhat arbitrarily) prescribed with an eye towards the freeze-out conditions extracted from data (20). The essential difference between these two trajectories is that UrQMD yields compressions that reach all the way into the deconfined phase region, whereas the maximum compression achieved with the 3-fluid model lies inside the unstable region. Such a situation would be obtained with UrQMD as well at a somewhat lower collision energy (3-4 ).

When the turning point of the phase trajectory lies inside the unstable region the modes are exposed to the spinodal amplification for a longer time and, presumably, the resulting degree of amplification would be maximized. This key fact is illustrated in Fig. 6. For the penetrating (UrQMD) trajectory, the modes are exposed to amplification only during the two relatively brief periods when the phase point is traversing the spinodal region (and the amplification gained during the compressional traverse is largely lost during the high-density stage due to the equilibration process so only the amplification received during the expansion traverse is relevant). By contrast, the more optimal (3-fluid) trajectory provides amplification for a sustained period of time, thus making a phase separation more likley to occur.

While it is yet difficult to make specific predictions about the optimal collision energy, it seems evident that such an energy range exists, since the location of the phase point associated with the maximum compression (the turning point) moves steadily downwards as the collision energy is lowered. The experimentally known freeze-out temperatures (indicated on the left in Fig. 5) provide a lower bound on the phase region that could be accessed by nuclear collisions. The turning point is therefore likely to traverse the unstable phase region at fairly high temperatures where it is relatively narrow. Consequently, the optimal range of collision energy is most likely not so wide. Because of this generic expectation, we shall take the above phase trajectories as being reasonably representative of what may happen in the bulk region of an actual collision.

### v.2 Dissipative collective dynamics

As the bulk of the evolving system traverses the unstable phase region, its density fluctuations may be amplified. In order to make quantitative estimates of this effect, we employ the method developed in Ref. (21) for the evolution of collective modes subject to a dissipative coupling to the environment, i.e. to the (many more) non-collective modes in the system.

The treatment considers an ensemble of macroscopically similar systems that have all been prepared with the same average density relative to which individual perturbations are introduced,

 ρ(n)(\boldmathr) = ρ0+δρ(n)(\boldmathr) = ρ0+∑\boldmath\scriptsizek≠0ρ(n)% \boldmath\scriptsizekei\boldmath\scriptsizek⋅\boldmath\scriptsizer. (39)

Generally, such as during the expansion stage of a nuclear collision, the bulk density is time dependent. We shall assume that the effect of this overall evolution may be taken into account approximately by using the corresponding time-dependent transport coefficients in the formulas below. We should then think of the mode index as a mode number rather than a wave number, since the expansion primarily stretches the modes without introducing much mode mixing (22). In this way we may obtain the spatial size of the resulting fluctuations by scaling the initial wave length with the linear expansion factor, , where is the effective dimensionality of the expansion.

The primary object of study is the associated spatial correlation function, , where

 ≺δρ(\boldmathr1)δρ(\boldmathr2)∗≻ ≡1N∑nδρ(n)(\boldmathr1) δρ(n)(\boldmathr2)∗ (40)

is the ensemble average and . The Fourier transform of provides a convenient measure of the degree of density fluctuation at a given scale, since it is the ensemble average of ,

 σ2\boldmath\scriptsizek = ≺|ρ\boldmath% \scriptsizek|2≻ =\ ∫d\boldmathrVe−i\boldmath\scriptsizek⋅\boldmath\scriptsizerσ(\boldmathr) . (41)

We are generally interested in the evolution of the collective modes in the system. Adopting the treatment of Ref. (21), we assume that the amplitude of a given collective mode is governed by an equation of motion having the form

 ddtρ\boldmath\scriptsizek(t) = −iω% \boldmath\scriptsizek(t)ρ\boldmath\scriptsizek(t)+B% \boldmath\scriptsizek(t) . (42)

The complex eigenfrequencies are determined by the dispersion equation derived above, while the Brownian term describes the residual coupling between the collective mode and the reservoir, i.e. the non-collective modes of the system. It is fluctuating in nature and is assumed to be Markovian,

 ≺B\boldmath\scriptsizek(t)B\boldmath\scriptsizek% ′(t′)∗≻ = 2D\boldmath\scriptsizek\boldmath\scriptsizek′(t)δ(t−t′) . (43)

The equal-time correlation coefficient for two collective modes, , is then given by (21),

 σ\boldmath\scriptsizek\boldmath\scriptsizek′(t)=σ\boldmath\scriptsizek\boldmath% \scriptsizek′(ti)e−iω\boldmath% \scriptsizek\boldmath\scriptsizek′t+2∫ttidt′D\boldmath\scriptsizek\boldmath% \scriptsizek′(t′)eiω\boldmath% \scriptsizek\boldmath\scriptsizek′(t′−t), (44)

where . It is readily seen that it satisfies the following differential equation,

 ddtσ\boldmath\scriptsizek\boldmath\scriptsize% k′(t) = −iω\boldmath\scriptsizek% \boldmath\scriptsizek′(t)σ\boldmath\scriptsizek\boldmath\scriptsizek′(t)+2D\boldmath% \scriptsizek\boldmath\scriptsizek′(t) , (45)

which was dubbed the Lalime equation in Ref. (21).

For our present studies, we are particularly interested in the time evolution of the diagonal components of the covariance matrix, , which are equal to the fluctuation coefficients introduced in Eq. (41). Then

 σ2\boldmath\scriptsizek(t) = [σ2% \boldmath\scriptsizek(ti)+∫tti2D\boldmath% \scriptsizek(t′)e−2Γ\boldmath\scriptsizek% (t′)dt′]e2Γ\boldmath\scriptsizek(t) , (46)

where and we have introduced the following amplification coefficient,

 Γ\boldmath\scriptsizek(t) ≡ ∫ttiIm[ω\boldmath\scriptsizek(t′)]dt′ = ∫ttiγ\boldmath\scriptsizek(t′)dt′ . (47)

If the environment is stationary, then and remain constant in time, so , and we obtain a simple exponential time evolution,

 σ2\boldmath\scriptsizek(t) = D\boldmath% \scriptsizekγ\boldmath\scriptsizek[e2γ\boldmath\tinyk(t−ti)−1]+σ2% \boldmath\scriptsizek(ti)e2γ\boldmath\tinyk(t−ti) . (48)

When is negative, as is normally the case, the mode is stable and relaxes exponentially towards its equilibrium value , , the associated relaxation time being . More generally, within the stable regime, will seek to relax towards its instantaneous equilibrium value , in accordance with the Einstein relation. In the opposite case, when is positive, the mode is unstable and the fluctuations exhibit an exponential growth, with both the original mode fluctuations and the noise being amplified.

The diffusion coefficients represent the coupling of the collective modes to the residual system. They thus play two important roles in the dynamical evolution of the density fluctuations. In the stable regime (as mentioned above) the coupling determines the relaxation times for the relaxation of the fluctuation coefficients towards their appropriate equilibrium values , while in the unstable regime it continually produces additional fluctuations that are also amplified. Thus, even in the absense of initial fluctuations, the coupling of the collective mode to the residual system will continually create fluctuations that will subsequently become amplified or damped, as governed by the Lalime equation (45). In particular, the diffusion coefficients enable the stable modes to continually adjust their fluctuations as the equilibrium variances evolve.

### v.3 Onset of phase separation

The amplification coefficient defined in Eq. (47) is obtained by integrating over the entire time interval considered. When there is no dissipation, is always real so vanishes outside the unstable region. But in the presence of dissipation, is negative outside the unstable region, and increasingly so the stronger the dissipation. As a result, the local relaxation time is reduced which in turn ensures that the fluctuations remain close to their equilibrium value until shortly before the trajectory enters the unstable region. However, by the same token, any amplification acquired while the trajectory is inside the unstable region is correspondingly quickly lost after the trajectory has reentered the normal region and again equilibrates. In fact, the dissipative length tends to be larger at late times, due to the increased particle spacing in the more dilute system, thus shortening the relaxation time.

Our present study aims to clarify the prospects for the spinodal amplification to become sufficiently large to cause a phase-separating clumping of the system. Since our treatment is perturbative, it is only valid as long as the fluctuations remain relatively small. The occurrence of large fluctuations in the calculation should then be taken as a signal that clumping is likely. It should be recalled that the bulk of the system is thermodynamcially unstable so even a modest degree of fluctuation may cause a catastrophic breakup. The prospects for this to occur are enhanced by the fact that fluctuations created in the mechanically unstable (spinodal) region may become futher amplified during the traverse of the adjacent metastable region, just as impurities introduced in this region may trigger condensation.

Thus the key question is how much amplification may occur during the unstable era. To elucidate this issue in a quantitative manner, we extract the following amplification factor,

 G\boldmath\scriptsizek ≡ exp(∫>γ% \boldmath\scriptsizek(t)dt) , (49)

where the integral is only over those times during which the mode is unstable, .

Figure 7 shows values of obtained for the entire range of wave numbers, as obtained with various degrees of dissipation. While we concentrate on the phase trajectory that reaches its maximum compression inside the spinodal region (the one based on the 3-fluid model results for ), we also show the result of a more penetrating trajctory to bring out the importance of tuning the collision energy for optimal effect.

In addition to the result of ideal fluid dynamics, we show results for various degrees of dissipation, ranging from minimal, i.e.