Dynamics of inhomogeneous chiral condensates

We study the dynamics of the formation of inhomogeneous chirally broken phases in the final stages of a heavy-ion collision, with particular interest on the time scales involved in the formation process. The study is conducted within the framework of a Ginzburg-Landau time evolution, driven by a free energy functional motivated by the Nambu–Jona-Lasinio model. Expansion of the medium is modeled by one-dimensional Bjorken flow and its effect on the formation of inhomogeneous condensates is investigated. We also use a free energy functional from a nonlocal Nambu–Jona-Lasinio model which predicts metastable phases that lead to long-lived inhomogeneous condensates before reaching an equilibrium phase with homogeneous condensates.


Introduction
The study of the phase diagram of strongly interacting matter as a function of the temperature and baryon chemical potential has been an active field of research since the early days of nuclear physics.In the early days, the field was driven by studies of cold nuclear matter and aimed at understanding the structure of atomic nuclei and the central regions of neutron stars.With the discovery of the quark-gluon substructure of nucleons and the establishment of quantum chromodynamics (QCD) as the theory of the strong interactions in the 1970s, the field entered into a new era.The asymptotic freedom property of QCD predicts the possibility of creating through high-energy heavy-ion collisions a state of matter resembling the one that presumably existed in the early universe, in that quarks and gluons are no longer confined in the protons and neutrons as in cold nuclear matter [1].
In the last few decades, great experimental effort has been devoted to the extraction from heavy-ion collisions observables relevant for the phase diagram of QCD.For the high temperature region of the phase diagram, results coming from RHIC-BNL and LHC-CERN suggest that a new state of matter, the quark-gluon plasma (QGP), is indeed formed after the collision event [2].Yet, it is not clear what is the nature of the phase transition between ordinary nuclear matter and the QGP.The beam energy scan program at RHIC, and programs at the new accelerator facilities FAIR and NICA that are expected to become operational in a few years, aim at getting clues on the nature of QCD phase transition, in particular on the chiral critical point and on its vicinity.One exciting possibility is the search of signals of the formation of spatially varying chiral condensates which break translational invariance, a subject that has raised a great deal of interest in the last years -for a review and references, see Ref. [3].In addition to the expected tricritical endpoint of the first order chiral phase transition, model calculations suggest the existence of a so-called Lifshitz point, where two homogeneous phases and one inhomogeneous phase meet [4,5].Interestingly, even for ordinary cold nuclear matter there seems to be a possibility that an inhomogeneous phase can be realized at densities a few times larger than the ground-state density [6].
Possible observable consequences of inhomogeneous phases have been suggested in the context of neutron stars.Ref. [7] proposes that the existence of an inhomogeneous chiral phase leads to a novel cooling scenario for neutron stars, finding a neutrino emissivity with efficiency comparable with the usual quark cooling or pion cooling.Consequences of an inhomogeneous chiral phase for the equation of state of matter in neutron stars have been investigated in Refs.[8,9].On the other hand, to the best of our knowledge, there has been no investigation on the possibility of formation of an inhomogeneous chiral phase in the matter produced in the late stage of a heavy-ion collision.Motivated by this fact, in the present communication we present results of a study of the dynamics of the formation of inhomogeneous chirally broken phases, with the aim of getting insight on the possibility of formation of such phases at the final stages of a heavy-ion collision.We are particularly interested in the time scales involved in the formation process.Our study is conducted within the framework of the phenomenological time-dependent Ginzburg-Landau equation.

Dynamics of order parameters
It is well known that static phase transitions in many-particle systems with short-range interactions depend on the spatial dimensionality and symmetry properties of the order parameters characterizing the different phase transitions in the system.On the other hand, the dynamics of phase transitions depends, in addition, on whether there are conservation laws associated with the order parameters and also on the couplings among the order parameters and conserved densities.Like in the case of static phase transitions, dynamical phase transitions can be grouped into universality classes, or models, in the nomenclature of Hohemberg and Halpering [10].We are interested in the dynamics of the chiral order parameters.When neglecting the coupling of the order parameters to conserved densities, like the baryon density and the energy-momentum tensor, one is lead to assume that the dynamics belongs to the model A dynamical universality class [10][11][12][13][14].
Denoting by φ(x, t) the order parameter field, a protype equation of motion for φ(x, t) is the phenomenological time-dependent Ginzburg-Landau (GL) equation where H[φ] is the Ginzburg-Landau-Wilson (GLW) Hamiltonian, which is typically of the form where r 0 and u 0 are constants that can depend on the temperature, the latter being necessarily positive at this level of truncation.For r 0 < 0, U 0 (φ) has a double-well shape.In addition, ζ(x, t) is the noise field describing thermal fluctuations and Γ is an Onsager (dissipation) coefficient.For simplicity, we consider a simple white-noise correlation function for the noise field: where • • • ζ means average over noise realizations and T is the external temperature.
From the Fokker-Planck equation associated with the GL equation of Eq. ( 1), it is not difficult to show that the equilibrium probability distribution of field configurations φ(x) is the Boltzmann factor The equilibrium partition function is then given by the functional integral The meaning of this is that a correlation function, e.g. a two-point correlation function φ(x 1 )φ(x ) , defined as can also be calculated via the long-time solutions of the GL equation as where φ r (x, t) is a solution of the GL equation of Eq. ( 1), with the index r indicating r-th noise realization, and N r is the total number of realizations (N r is supposed to be large).Despite straightforward, the numerical solution of Eq. ( 1) requires some care with respect to the lattice spacing and time steps used in the discretization of the derivatives.One way to minimize the dependence on the lattice spacing is to use a renormalization procedure, as explained e.g. in Refs.[15][16][17].When considering the dynamics of a conserved density, e.g. the baryon number, one would be lead to consider a model B dynamics [18,19].Equations of the Ginzburg-Landau type can be derived from a microscopic model via the closedtime-path (CTP) effective action formalism [20].In general, the equations obtained with the CTP formalism contain nonlocal (non-Markovian) dissipation kernels and colored noise fields, as well as the possibility of field-dependent noise terms (multiplicative noise) [15,[21][22][23].In a phenomenological approach, which we employ in the present study, one takes the equilibrium Hamiltonian H[φ] and so the description is valid for situations not too far from equilibrium.In addition, the dissipation coefficient Γ that enters Eq. ( 1) must be obtained from independent estimates.
To illustrate a typical solution φ r (x, t) of Eq. ( 1) for a given noise realization, let us consider a double-well potential U 0 (φ).We take as an initial configuration φ(x, 0) that assumes a random value at each point x = (x, y, x), with amplitude φ 0 (x, 0)/φ 0 = 0.1 and zero average, where φ 0 = (|r 0 |/u 0 ) 1/2 is the positive minimum of U 0 (φ). Figure 1 shows density plots of a typical field configuration φ r (x, y, z = 0, t)/φ 0 at different time instants t/t 0 , where t 0 is a typical time scale.At t = 0, left panel of the figure, one has the initial random configuration with zero noise average.As t/t 0 increases, central e right panels, one starts seeing very clearly the inhomogeneous nature of the solutions.For t/t 0 very large, the system will eventually evolve to a nearly homogeneous equilibrium configuration for which φ r (x, t)/φ 0 is positive or negative, with |φ r (x, t)| φ 0 .The equilibrium value is not precisely the "bare vacuum" φ 0 because of the thermal fluctuations.

Chiral order parameters -statics
We consider NJL-type models with vanishing quark masses.The Lagrangian density defining the models is written generically as where ψ stands for the N f = 2 fermion doublet, and J a (x) are currents of the form where Γ a = (1 1, iγ 5 τ), and G(z) is a form factor that characterizes the effective interaction.For the traditional, local NJL model, teh form factor is given by G(z) = δ(z).One should be aware of the fact that NJL-type models are meant to model the chiral transition of QCD, but they miss two prominent features of QCD, asymptotic freedom and quark-gluon confinement, that certainly are very important at high densities.Models developed in the past that incorporate those features, as those of Refs.[25][26][27], have achieved relative success in describing vacuum hadron properties, and provide a good alternative for studying the confined phase.Still another alternative, although technically more elaborate, is the use of Dyson-Schwinger equations of QCD at finite density [28], like the recent study in Ref. [29].
To obtain an equilibrium GLW Hamiltonian as a functional of the order parameters to be used in Eq. ( 1), we perform a standard bosonization of the theory by integrating out the quark fields in favor of scalar and pseudoscalar fields.We employ the mean field approximation, in that the bosonic fields are replaced by their vacuum expectation values φ(x) = (σ(x), π(x)).Following Refs.[4,30], the mean field thermodynamic potential is expanded around the symmetric ground state in powers of the order parameters and their spatial gradients and define it as being the GLW Hamiltonian density where the α n , n = 2, 4, 4b, 6c, 6d depend upon the temperature T and baryon chemical potential µ; their expressions are given in Ref. [30] and their explicit behavior as a function of T and µ will be presented elsewhere [31].
Next, we consider the phase diagram from a nonlocal NJL model for a particular set of parameters [36].The model is completely determined by the form factor G(z) and the coupling constant G.The form factor is taken to be a Gaussian; its Fourier transform is g(p) = exp(−p 2 /Λ 2 ), where the parameter Λ is related to the range of the interaction.Those parameters are fixed by fitting the pion decay constant and the quark condensate in the chiral limit: f π = 86 MeV and qq = −(270 MeV) 3 .The numerical values of G and Λ are G = 14.668GeV −2 and Λ = 1.046GeV.Fig. 3 shows the corresponding phase diagram for a spatial modulation in the form of a dual chiral density wave, i.e. a plane wave for the scalar and pseudoscalar fields -see Ref. [36] for details.Although a stable inhomogeneous phase exists for very high densities only, there exists, however, a local stable minimum of the GLW Hamiltonian of Eq. ( 11) with an inhomogeneous condensate, but it is located within the extremely narrow band enclosed by the dotted and dash-dotted lines.This latter situation is very interesting, in that for a T and µ near the tricritical and Lifshitz points -e.g. the point (µ, T ) = (190 MeV, 50 MeV) identified in the figure -there appear inhomogeneous order parameters as solutions of the corresponding Ginzburg-Landau equation that are long-lived, with a life time much larger than Γ, as we discuss in the next section.It is important to note that the shape of the phase diagram and the location of the inhomogeneous minimum depend heavily on the value of the chiral condensate.

Chiral order parameters -dynamics
We have two order parameters, σ and π.The corresponding time-dependent GL equations are coupled due to the the nonlinear terms in H NJL [σ, π]; they are given by where ζ σ and ζ π are white-noise fields.We start presenting solutions for the local NJL model and a static medium -expansion of the medium is considered in the following section.We take π = 0 and, to avoid proliferation of different choices of parameters, we study the line α 4 = − √ 36/5 α 2 α 6 ; as discussed above, on this line the solution given in Eq. ( 12) becomes a hyperbolic tangent.Dimensionful quantities are rescaled with with the appropriate powers of Γ.We choose α 2 = α −1 6 = Γ −2 , then α 4 = − √ 36/5.The initial profile φ(z, 0) is chosen to be a random configuration generated with Gaussian noise, with the same parity of the inhomogeneous equilibrium solution -recall that on the line α 4 = − √ 36/5 α 2 α 6 , the homogeneous and inhomogeneous solutions coexist.As the inhomogeneous solutions occur at relatively low temperatures, the noise field ζ can be neglected.
t/Γ= 0.00 t/Γ= 10.00In Fig. 4 we plot the solutions φ(z, t) at different times t.Clearly, the long-time solution resembles the equilibrium hyperbolic tangent solution, achieved after a time t ∼ 10 Γ.Since we do not have a microscopic derivation of Eq. ( 1), the value of Γ must be taken from elsewhere.A good estimate for Γ, calculated from the decay σ → ππ, is Γ = 1/3 fm [22,37] at T = 0, which is the relevant temperature in the present case.For this value of Γ, equilibrium is achieved after a time interval of t = 30 fm.In the next section we compare this time scale with the expansion rate.It is important to mention that in principle there is a minimum size of the system for the decay σ → ππ to happen; if the size is too small, the decay cannot happen, there is no dissipation and fluctuation and, therefore, Γ would be zero.
We consider now the coupled system of GL equations, Eq. ( 13), for the non local NJL model.We evaluate the α n coefficients of the GLW Hamiltonian in Eq. ( 11) for (µ = 190 MeV, T = 50 MeV).This (µ, T ) point, identified on the phase diagram of Fig. (3), is located within the chirally broken homogeneous equilibrium phase, but it lies very close to the tricritical and Lifshitz points.We take as initial profiles σ(z, 0) and π(z, 0) dual chiral density waves superimposed by Gaussian noise.Fig. 5 shows the solutions for σ(z, t) and π(z, t) = π 3 (z, t).Clearly, the long-time solutions are the expected homogeneous solutions π = 0 and σ 0, the latter indicating that chiral symmetry is dynamically broken.The interesting feature here is that although the equilibrium order parameters are homogeneous, homogeneity is reached only after a very long time, t > 10 4 Γ.The investigation of the dependence of the lifetime of the inhomogeneous solutions on the size of the system is reserved for a future publication.

Chiral order parameters -dynamics with expansion of the medium
Here we address the evolution of the order parameter in an expanding medium using Bjorken flow [38].A more realistic simulation of the expansion would demand a full hydrodynamical simulation.The appropriate coordinates to describe longitudinal expansion are the so called Milne coordinates [39], the longitudinal proper time τ = (t 2 −z 2 ) 1/2 and space-time rapidity η = 1/2 log where t and z denote the time and longitudinal coordinates in the laboratory frame.We must transform the derivatives in the GL equation accordingly.The initial conditions are set at a time-like surface τ = τ 0 .The value of τ 0 sets how late the GL evolution starts after the beginning of the expansion.We consider the local NJL model only.Initially, we consider the homogeneous equilibrium phase and calculate the time evolution of the spatial average of the order parameter φ where L is the spatial extension of system.Again we simulate on the first critical line, α 4 = − √ 36/5 α 2 α 6 .The result for φ as a function of ∆τ ≡ τ − τ 0 is shown in the left panel of Fig. 6 for τ 0 = 10 Γ and τ 0 = Γ, as well as for the nonexpanding case.As expected, the sooner the GL dynamics starts with respect to the start of the expansion, the more delayed is the approach to the equilibrium configuration compared the to static case.In the right panel of the figure, we present the result for a delayed dynamics, i.e. for a case with τ 0 = 10 4 Γ.The equilibrium configuration is obtained for an elapsed time ∆τ = 12 Γ.

Conclusions and Perspectives
We have presented results of a study of the dynamics of the formation of inhomogeneous chirally broken phases.With the aim of getting insight on the possibility of formation of such phases at the final stages of a heavy-ion collision, we have investigated the time scales involved in the formation process.
Our study was based on the framework of the phenomenological time-dependent Ginzburg-Landau equation.The Ginzburg-Landau dynamics depends on two fundamental quantities: the energy functional that gives the equilibrium phase diagram and the Onsager coefficient Γ, a time scale related to fluctuation and dissipation processes.We employed NJL-type models to obtain the energy functional; we considered a local NJL model with point-like interactions [4] and a nonlocal NJL model with finite-range interactions [30].The results are qualitatively similar for both models, but there are interesting differences.In particular, we have shown that even in the case that the equilibrium phase is a homogeneous chirally broken phase, long-lived inhomogeneous phases can be formed when the model energy functional predicts metastable phases, like in the case of the nonlocal NJL model.For Γ = 1/3 fm, that is the value found in Ref. [37] from an estimate of the decay σ → 2π at zero temperature, we found that, for appropriate initial configurations, well-defined inhomogeneous domains are formed for a time interval of the order of ∆t = 10 fm, in both models.We have also considered expansion of the medium, a crucial feature of heavy-ion collisions that must be taken into account for a realistic description of the dynamics.We modeled the expansion by one-dimensional Bjorken flow and found, as expected, that the expansion delays the formation of inhomogeneous phases.Our work must be improved on several aspects.First of all, one needs to consider dynamics in three spatial dimensions.One needs also to investigate the role of confinement and asymptotic freedom in the phase diagram.For that, one needs models that incorporate such features and go beyond NJL-type models.In addition, possible observable signals linked to the decay of the inhomogeneous condensates must be found [40].Finally, one needs to go beyond the phenomenological Ginzburg-Landau approach and derive dynamical equations from microscopic models using, for example, the closed-time-path effective action formalism [20].

Figure 1 .
Figure 1.Typical spatial profiles of the order parameter -see text for explanations.Figure taken from Ref. [24].

Figure 3 .
Figure 3. T × µ phase diagram for a particular parametrization of the nonlocal NJL model [36].Solid (dashed) lines indicate first (second) order phase transitions.The dotted line corresponds to the homogeneous chiral restoration transition.The dash-dotted line is a boundary of a region in which there exists a local inhomogeneous minimum of the thermodynamical potential.HCB (HCR) and IH stand for homogeneous chirally broken (restored) and inhomogeneous phases.TCP and LP stand for tricritical and Lifshitz points.The values (µ = 190 MeV, T = 50 MeV) are those used in the evaluation of the α n in Eq. (11) for the dynamical simulation.Figure adapted from Ref. [36].

Figure 4 .
Figure 4. Results for the simulation on the first critical line, α 4 = − √ 36/5 α 2 α 6 .The initial condition is generated with a Gaussian noise and it has the parity of the equilibrium solution.

EPJ−Figure 5 .
Figure 5.Time dependence of spatial profiles of the pseudoscalar (left panel) and scalar (right panel) order parameters.The initial profiles are shown by the light dotted lines.