Can dissipation prevent explosive decomposition in high-energy heavy ion collisions?

We discuss the role of dissipation in the explosive spinodal decomposition scenario of hadron production during the chiral transition after a high-energy heavy ion collision. We use a Langevin description inspired by microscopic nonequilibrium field theory results to perform real-time lattice simulations of the behavior of the chiral fields. We show that the effect of dissipation can be dramatic. Analytic results for the short-time dynamics are also presented.

Lattice QCD results [1] suggest that strongly interacting matter at sufficiently high temperature undergoes a phase transition (or a crossover) to a deconfined quark-gluon plasma (QGP). Despite the difficulties in identifying clear signatures of a phase transition in ultrarelativistic heavy ion collisions, recent data from the experiments at BNL-RHIC clearly point to the observation of a new state of matter [2].
Depending on the nature of the QCD phase transition, the process of phase conversion in the expanding QGP generated in a high-energy heavy-ion collision may proceed in a number of different ways. In general, there is a competition between the mechanisms of nucleation and spinodal decomposition [3]. Results from CERN-SPS and BNL-RHIC feature what has been called sudden hadronization [4] or explosive behavior [5] in the hadronization process of the expanding QGP and seem to favor a fast (explosive) spinodal decomposition.
Recently, possible signatures of this behavior in high-energy nuclear collisions were proposed [6].
Most theoretical attempts to understand this behavior focused on the rapid changes in the effective potential of QCD near the critical temperature such as predicted, for instance, by the Polyakov loop model [7], which would be followed by a very fast spinodal decomposition process [5]. In this case, one observes a nearly instantaneous decay of the Polyakov loop condensate, producing an exploding source of pions as well as large RMS spatial fluctuations in the chiral fields. Previous studies using low-energy chiral effective models to investigate how nucleation rates compare to the time scale of expansion of the plasma also found very likely that most of the system would reach the spinodal region and then undergo an explosive phase conversion [8,9,10,11,12,13]. Therefore, both experiment and theory seem to suggest that a QGP generated in a high-energy heavy ion collision expands and cools down so fast that the processes of chiral symmetry breakdown and hadronization can be described within the framework of an effective potential that is quenched into the spinodal region, where long-wavelength fluctuations grow with no barrier to overcome. This leads to what we will refer to as the explosive spinodal decomposition scenario. However, we will argue that even if the system quickly reaches this unstable region there is still no guarantee that it will explode. To asses the different possibilities one has to study the time evolution of the order parameter of the transition after the system is quenched beyond the mean-field analysis, in order to incorporate the effects from dissipation, and ask how they could modify the explosive decomposition picture. In particular, dissipation effects have proved to be important in the context of disoriented chiral condensate (DCC) formation in heavy ion collisions [14,15,16,17,18].
In this Letter we present an exploratory investigation of the effects of dissipation in the explosive spinodal decomposition scenario of hadron production during the QCD transition after a high-energy heavy ion collision. We use a Langevin description inspired by microscopic nonequilibrium field theory results to perform real-time lattice simulations for the behavior of the inhomogeneous chiral fields. We show that the effects of dissipation can be dramatic even for very conservative assumptions. Analytic results for the short-time dynamics are also presented and discussed.
For the sake of simplicity, we consider an infinite system that is quenched to the spinodal and then evolves with a fixed effective potential. We believe that the subsequent evolution of the effective potential should not bring deep modifications to our general conclusions once the system has already reached the unstable spinodal region. On the other hand, effects from the finite size of the plasma will most likely play a non-trivial role [19]. Both effects will be taken into account in a future publication [20].
In what follows, we consider the real-time dynamics of chiral symmetry breakdown of a QGP created in a high-energy heavy ion collision. We assume that the system is characterized by a coarse-grained free energy where U ef f (φ, T ) is an effective potential of the Landau-Ginzburg form whose coefficients depend on the temperature, and φ( x, t) is a real scalar field which plays the role of an order parameter that is not conserved, such as the chiral condensate. To model the mechanism of chiral symmetry breaking found in QCD, we adopt the linear σ-model coupled to quarks, whose standard Lagrangian can be found, for instance, in Ref. [11]. This approach is widely used in the literature and its specificities imply no major limitations to our main results.
Quarks can be treated as fast-moving modes and integrated out using a classical approximation for the chiral field, yielding the effective potential U ef f (φ, T ) that is shown in Fig. 1. The pion directions play no major role in the process of phase conversion we are considering, so we concentrate on the sigma direction represented by the field φ (see Ref. [11] for details). In Fig. 1, we show U ef f (φ, T ) for a few values of temperature. The critical temperature is T c ≈ 123 MeV and the spinodal is reached at T sp ≈ 108 MeV. In order to simplify the numerics and to compare with analytic results for the short-time behavior of the order parameter evolution, we work with a fit of U ef f (φ) by a polynomial of sixth degree whose coefficients depend on the temperature. Such a polynomial fit is an almost perfect representation of U ef f , and it is the fit that is displayed in Fig. 1.
In our analysis, the evolution of the order parameter φ( x, t) and its approach to equilibrium will be dictated by a Langevin equation of the form where Γ, which can be seen as a response coefficient that defines a time scale for the system and encodes the intensity of dissipation, will be taken to be a function of temperature only, Γ = Γ(T ). The function ξ( x, t) represents a stochastic (noise) force, assumed Gaussian and white, so that ξ( x, t) = 0 and ξ( Eq. (2) could, in principle, be obtained from a microscopic field-theoretic description of the real-time nonequilibrium dynamics of the chiral field at finite temperature. This procedure was implemented in the case of a λφ 4 scalar field theory in Refs. [21,15]. The noise and dissipation terms, which originate from quantum fluctuations, are engendered by either self-interactions of the chiral field or coupling to one or more different fields that play the role of a heat bath, provided one incorporates higher order terms in the computation of the effective equation of motion for φ( x, t). In fact, it is well-known that one has to go up to two-loop corrections in order to pick up imaginary parts in the self-energy associated with viscosity and dissipation. Self-interactions of the chiral field, as well as its interactions with quarks and anti-quarks, justify the inclusion of a dissipation and a noise term such as done in Eq. (2) in the framework adopted in this paper. The same is true for an approach that also includes the dynamics of the Polyakov loop condensate coupled to the chiral field and, as we will argue, the inclusion of this effect could dramatically modify the results of Refs. [5,7].
The main physics behind Γ is the decay of the φ field and as such one would think that Γ would be significant only in regions where φ oscillates rapidly. However, this needs not be the case in general. The procedure detailed in Refs. [21,15] leads to an effective equation of motion that is richer in structure and much more complicated than (2). There is no a priori reason for the dissipation function Γ be so simple and the noise be Gaussian and white. In general, one obtains a complicated dissipation kernel that simplifies to a multiplicative dissipation term which depends quadratically on the amplitude of the field as η 1 (T )φ 2 ( x, t)φ( x, t), where η 1 is determined by imaginary terms of the effective action for φ and depends weakly (logaritmically) on the coupling(s). The fluctuation-dissipation theorem implies, then, that the noise term will also contain a multiplicative contribution of the form φ( x, t)ξ 1 ( x, t), and be in general non-Markovian. The white-noise limit is reobtained only in the limit of very high temperature. Nevertheless, assuming Γ to be a linear function of the temperature is a reasonable first approximation, as can be seen from the results presented in [15,16]. We will comment further on this point later on when discussing our results.
For the sake of simplicity, we adopt the simple approximate form of Eq. (2) for a phenomenological description of the dissipative evolution of the expectation value of the sigma field. Although assumedly simple, this analysis allows for a clear distinction and comparison of the roles played by dissipation and the (explosive) spinodal instability in the spinodal decomposition scenario of hadron production during the QCD transition after a high-energy heavy ion collision. The simple form of Eq. (2) is also convenient for a comparison of our numerical results to (linear) analytic estimates in the region of short-time evolution to measure the effect of nonlinearities.
In our numerical simulations we solve Eq. (2) on a cubic space-like lattice with 64 3 sites under periodic boundary conditions, with a lattice spacing of a = 0.91 fm. We use a semi-implicit finite-difference scheme for the time evolution and a Fast Fourier Transform for the spatial dependence [22]. Temperature is fixed to the spinodal value T sp ≈ 108 MeV. We perform several runs starting from different random initial configurations around the inflexion point of U ef f which happens at φ 0 ≈ 0.162 T and then average the results from the different initial configurations. For time steps of ∆t = 0.001/T the results become independent of the lattice spacing once it is smaller than a ≃ 1 fm.
Before presenting the results of the simulations, it is instructive first to analyze the short-time behavior of the solution. One can linearize the equation around the inflexion point φ 0 substituting φ by φ = φ ′ +φ 0 in Eq. (2) and average over the noise. For short times, φ ′ is small and the cubic and higher-order terms in U ef f (φ ′ ) can be neglected, so that the equation for φ ′ becomes linear. Since the equation is linear, the average over the noise can be done formally. An analytic form for the short-time solution of the linear equation can be found by using the polynomial fit of U ef f , which leads to following equation for the average φ ′ where A = 2a 2 + 6a 3 φ 0 + 12a 4 φ 2 0 + 20a 5 φ 3 0 + 30a 6 φ 4 0 . Note that the constant term in U ′ ef f does not contribute to the equation Eq. (4) since the average over noise of a constant is zero. One can write the solution of Eq. (4) in terms of the Fourier transform φ ′ (k, t) of φ ′ (x, t) as where C 1 and C 2 are integration constants and λ 1 (k) and λ 2 (k) are the roots of the quadratic equation where we used the fact that A < 0 in our case. From this one sees that for short wavelengths, such that k 2 ≫ Γ 2 /4 + |A|, we have (complex conjugate) imaginary roots and the solution is oscillatory. For long wavelengths, such that k 2 < Γ 2 /4 + |A|, we have real roots and there is an exponential growth of the Fourier components. This exponential growth yields the explosive spinodal decomposition.
Two other limits are also instructive. One is the strong dissipation limit of large Γ, such that the first order time derivative dominates over the second order one. In this case, the short-time solution to Eq. (4) is given by In this case, one sees that short wavelengths with k 2 > |A| are absorbed by the system, while those with k 2 < |A| explode exponentially. Of course, as time increases, φ ′ increases and the linear equation is not valid anymore and the fully nonlinear equation has to be solved. The other interesting limit is the one with Γ = 0, i.e. no dissipation. We will discuss this limit in connection with the solution of the full nonlinear equation in the following.
We show results of simulations for three different values of the dissipation coefficient, namely Γ = 0, 2 T and 4 T . It can be argued that the response coefficient has the form Γ(T ) ≈ 2 T /b, where b is a number of order one to first approximation [23]. The cases considered provide a conservative band around the value Γ(T ) ≈ 2 T to illustrate the effect of dissipation.
In Fig. 2 we compare the solutions of the full, nonlinear equation with the solution of the linearized equation for Γ/T = 2. Fig. 2 (a) shows that the rolldown for the O(φ ′2 ) potential is slower than for the full U ef f potential. This is obviously due to the fact that the falloff of full U ef f is steeper that that of the O(φ ′2 ) for small φ ′ . At short times, smaller than t = 5 fm −1 , one sees that both solutions are very close to each other. One interesting aspect of the solutions shown in Fig. 2 (b) is that the exponential explosion of the linear solution happens much later than the explosion of the full nonlinear equation. This seems at first sight very counterintuitive since, from the discussion above, one would expect an early-time explosion of the solution of linearized solution. This does not happen here because the O(φ ′2 ) potential is much shallower than the full U ef f for small φ ′ . In Figure 3 we show the average value of φ in units of its vacuum value, φ vac , as a function of time for the three different values of Γ mentioned above. In the absence of dissipation, the solution (dotted line) is obviously oscillatory, with an explosive early-time behavior. The effect of dissipation is of retarding the exponential growth, as shown by the solid and dashed lines in the Figure. The results clearly show that even for a very conservative value of dissipation, Γ = 2 T , the effect can be dramatic. For this value of Γ, dissipation retards the time evolution of φ towards its vacuum value in ∼ 100% compared to the case with Γ = 0. The important point to be noted here is that for expansion times of the order of 5 fm −1 , which in of the order of the time scales for the RHIC collisions, there might be not enough time for the onset of the spinodal explosion [5].
Of course, effects brought about by the expansion of the plasma [5] and by its finite size [19], as well as a more realistic treatment of dissipation from the microscopic point of view, will bring corrections to this picture. For instance, the authors of Ref. [5] consider a Hubble expansion of the system which introduces a dissipation-like term to the evolution equation of the form H(∂φ/∂τ ), where τ is the proper time and H = 1/τ is the expansion rate. Therefore, for a very rapid expansion, corrections due to dissipation, such as discussed here, should play a comparatively less important role. Also, dissipation being mainly the result of the decay of the φ field, common wisdom would suggest that its effect should be less important at short times, when the field is slowly starting to roll down the potential. However, as discussed earlier, Γ depends, in general, on φ -in a λφ 4 model it is proportional to φ 2 -and it is not a priori clear what will be the effect of such a dependence on the short-time evolution of the system. These issues will be addressed in a future publication [20].