Time-dependent Ginzburg-Landau approach to the dynamics of inhomogeneous chiral condensates in a nonlocal Nambu--Jona-Lasinio model

We study the dynamics of inhomogeneous scalar and pseudoscalar chiral order parameters within the framework of the time-dependent Ginzburg-Landau equations. We utilize a nonlocal chiral quark model to obtain the phase diagram of the model as function of temperature and baryon chemical potential and study the formation of metastable spatial domains of matter where the order parameters acquire a spatial modulation in the course their dynamical evolution. We found that, before reaching the expected equilibrium homogeneous state, both scalar and pseudoscalar chiral condensates go through long-lived metastable inhomogeneous structures. For different initial configurations of the order parameters, the lifetimes of the inhomogeneous structures are compared to timescales in a relativistic heavy-ion collision.


I. INTRODUCTION
The phase diagram of strongly interacting matter at finite temperature (T ) and baryon chemical potential (µ) has been extensively studied along the last decades. Quantum Chromodynamics (QCD) predicts that at very high temperatures (T Λ QCD ) and low baryon densities this matter appears in the form of a plasma of quarks and gluons, and at very high baryon densities (µ Λ QCD ) and zero temperatures it is a color superconductor [1]. At such extremes, QCD is weakly coupled and first-principles perturbative calculations based on an expansion in the coupling constant can be used to explore the phase diagram. But there are regions of the phase diagram with temperatures and densities interpolating the extremes that remain poorly understood. An example is the region of low temperatures and densities close to or a few times larger than the nuclear saturation density, which is of great interest for the physics of compact stars [2][3][4]. Here, QCD is strongly coupled and coupling-constant expansions become inapplicable. First-principles nonperturbative lattice QCD methods based on large-scale Monte Carlo simulations are also not applicable in this case, because at finite µ QCD has a sign problem [5,6].
Lattice QCD has established, in particular, the existence of a finite-temperature crossover for chiral symmetry restoration at vanishing small densities [7][8][9]. In vacuum, the approximate chiral symmetry of the QCD Lagrangian in the light quark sector is dynamically broken, a feature that explains the lightness of the pion and is responsible for generating the bulk of the masses of the light hadrons, like protons and neutrons [10,11]. On the other hand, the behavior of chiral symmetry when moving from vacuum to densities of ordinary nuclear matter and higher is not well understood and all of what is presently known comes from model calculations. In this context, recent works [12][13][14][15][16][17] have revived the discussion on the possibility that chiral symmetry breaking in dense matter at low temperatures would drive the formation of nonuniform phases, i.e. the formation of spatially-varying chiral condensates which break translational invariance-Ref. [18] is a thorough recent review on the subject, with an account on earlier developments and a large list of references. One particularly interesting result suggests that, in addition to the expected tricritical endpoint of the first order chiral phase transition, there might exist a Lifshitz point, where two homogeneous phases and one inhomogeneous phase meet [15,17,19]. Interestingly, also for ordinary cold nuclear matter there seems to be the possibility that an inhomogeneous chiral phase appears at densities a few times larger than the normal density [20]. These are interesting features of dynamical chiral symmetry breaking in matter and it would be fascinating to find signals of their existence in real systems.
Recent studies [21][22][23] have suggested possible observable signals of the presence of inhomogeneous phases in hybrid stars. In Ref. [21] a novel cooling scenario was suggested, in that inhomogeneities induce modifications in momentum-conservation relations in quark beta decay, leading to neutrino emissivities with efficiencies comparable to those due to interacting quarks or due to the presence of a pion condensate. Refs. [22,23] investigated the consequences of an inhomogeneous chiral phase for the equation of state of matter in a hybrid star, finding substantial effects on the mass-radius relation of the arXiv:1804.05724v1 [hep-ph] 16 Apr 2018 star. Clues on such phases are also expected from experiments of heavy-ion collisions, like those ongoing within the beam energy scan program at Relativistic Heavy Ion Collider (RHIC) [24] and also from those planned to be conducted at the Facility for Antiproton and Ion Research (FAIR) [25], the Nuclotron-based Ion Collider Facility (NICA) [26] and the Japan Proton Accelerator Research Complex (J-PARC) [27]. If metastable inhomogeneous phases are present in the matter produced in a heavy-ion collision, the system would spend different times in the different phases during its evolution and therefore, this information could in principle be recovered employing for instance the freeze-out eccentricity which provides geometric information [24]. There has also been the suggestion that (thermal and quantum) fluctuations around a mean-field inhomogeneous condensate induce anomalies in thermodynamic quantities that could be revealed in particle production yields [28].
Parallel to the quest of observable signals in heavyion collisions, there is the important issue regarding the dynamical evolution of an inhomogeneous chiral configuration in this context. Particularly important is the time scale associated with the nonequilibrium evolution of such a configuration from its formation until its decay into freely steaming particles (mostly pions) that eventually will reach the detectors. Clearly, a firstprinciples theoretical study of this issue is currently out of reach. Nevertheless, likewise with the equilibrium situation, useful qualitative insight on the nonequilibrium evolution can be obtained through the Ginzburg-Landau (GL) mean-field approach. In this approach, given the GL functional, i.e. the thermodynamic potential, functional of the order parameters, the phase structure of the system is obtained by exploring the extrema of the potential with respect to the order parameters. When the system is brought out of equilibrium, the relaxation of the order parameters, i.e. their time evolution towards to an equilibrium state is described by the time-dependent Ginzburg-Landau (TDGL) equation [29]-a good review on the TDGL approach for condensed matter systems can be found in Ref. [30]. This approach and variants of it have been extensively used along the last decades to investigate different aspects of the relaxation dynamics of chiral order parameters and also order parameters associated with conserved charges and color confinement [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]. In the present paper, to get insight into the dynamics of inhomogeneous chiral condensates, we follow those lines and employ the TDGL approach. In general, to obtain the T − and µ−dependence of the GL functional, a model is required. Here we employ a nonlocal Nambu−Jona-Lasinio model (nlNJL) [51,52], in which quark fields interact through a non local chiral invariant four-fermion coupling, to obtain the thermodynamic energy functional. We focus on the time evolution of the the chiral order parameters at low values of the T , and different values of µ, both close to the tricritical and the Lifshitz points, as predicted by that model.
The paper is organized as follows. In the next section we define the nonlocal NJL model we use and briefly review the TDGL approach. Section III presents the results of the numerical simulations of the TDGL equations. Finally, Conclusions and Perspectives are presented in section IV.

II. TDGL APPROACH
As in Ref. [51], we consider the simplest version of a nonlocal SU(2) NJL model in the chiral limit. The corresponding Euclidean effective action is given by where ψ stands for the N f = 2 fermion doublet ψ = (u, d) T , and j a (x) for the nonlocal currents where we have defined Γ a = (1 1, iγ 5 τ ), and the function G(z) is a nonlocal form factor that characterizes the effective interaction.
To proceed, we perform a standard bosonization of the theory, in which bosonic fields are introduced and quark fields are integrated out. Within the GL approach, the bosonic fields are replaced by their (vacuum or thermodynamic) expectation (or mean-field) values φ( x) = (σ( x), π( x)). Since now parity is not necessarily an exact symmetry, one can get in general a nonzero value for the pseudoscalar field. The mean-field values are allowed to be inhomogeneous, hence the explicit dependence on spatial coordinates. The phase structure of the system is obtained by exploring the extrema of the GL functional, the thermodynamic potential Ω GL obtained from the action in Eq. (1), functional of the order parameter φ: This gives a coupled system of equations for σ and π, whose solutions reveal the presence, or absence, of inhomogeneous configurations of the fields σ and π in thermodynamic equilibrium at given values of T and µ. The question to be addressed is the time evolution of an initial configuration φ, which is not a fully-developed equilibrium configuration. Within the TDGL framework, the time evolution of φ is governed by the dynamical equation [29]: where Γ is a kinetic coefficient due to dissipation processes, like σ ↔ 2π and πσ → π; it is T − and µ−dependent and, in principle, different for σ and π.
The equilibrium configurations are the ∂φ/∂t = 0 stationary solutions for t → ∞, i.e. the mean-field equation Eq. (3). The equilibrium configurations are independent of Γ, but the time scale for them to become established is governed by Γ.
We note that thermal fluctuations are not taken into account in Eq. (4). Although they are small for low temperatures, as in the case we are going to study in the present paper, they play an important role at temperatures e.g. close to the crossover temperature. Fluctuations are usually taken into account phenomenologically by adding noise terms on the right hand side of Eq. (4), turning the TDGL equation into a stochastic equation [32,35,45,46,48]. The strength of the noise fields are constrained by the dissipation-fluctuation theorem. Equations of this sort can be derived from a microscopic model via influence functional techniques or the closed-time-path (CTP) effective action formalism [53]; in general, they contain nonlocal dissipation kernels and colored noise fields-see e.g. Refs. [54][55][56][57][58][59][60][61]. We also note that Eq. (4) is a TDGL equation for nonconserved order parameters, like for σ and π in the present case; for conserved order parameters, like the baryon density, different phenomenological equations are used [29].
As mentioned earlier, the situation of interest is the time evolution of the order parameters starting from initial configurations for which chiral symmetry is not fully restored. The physical picture behind such a scenario resembles the cooling stage in the course of the evolution of a heavy-ion collision toward a state of broken chiral symmetry. To make contact with the studies at equilibrium in the literature, we follow Refs. [15,51] and expand the mean-field thermodynamic potential in powers of the order parameters and their spatial gradients as: where ω GL (T, µ, φ) is the GL energy-density functional with the expansion coefficients α 2 , · · · , α 6d , which are functions of T and µ, given by where we have used the shorthand notation and The function g, evaluated at p 2 = p 2 n , is the Fourier transform of the non local form factor G(x) of the quark-antiquark currents, and g and g denote derivatives with respect to p 2 . We recall that the GL expansion in powers of the order parameter and its gradients, is expected to be valid only close to the second order transition to the chirally restored phase [62][63][64][65]. Moreover, at the Lifshitz point, both the order parameter and its gradient are expected to vanish. As in Refs. [15,51], we also restrict the analysis to one-dimensional modulations of the order parameters. In this case, the two TDGL equations for σ and π are given by Γ π ∂π( x, t) ∂t = −α 2 π + α 4 (π 2 + σ 2 ) − α 6 (π 2 + σ 2 ) 2 π + α 4b π + α 6b 3 where the primes denote spatial derivatives.

III. NUMERICAL RESULTS
In the chiral limit, the model is completely determined by the form factor g(p) and the coupling constant G. We choose the Gaussian form for g(p), considered in many previous studies [66][67][68][69][70][71][72]: where Λ is the range of the interaction in momentum space. It is usual to fix those parameters so as to get phenomenologically adequate values for the pion decay constant and the quark-antiquark condensate. Here, we use f π = 86 MeV and qq = −(270 MeV) 3 [73], to determine Λ and G. They are given by G = 14.668 GeV −2 and Λ = 1.046 GeV.
The explicit T and µ dependence of the GL coefficients α 2 , · · · , α 6d have not been presented previously for a nonlocal model, therefore in Fig. 1 we plot them as a function of temperature for different values of the chemical potential. The figure shows the dimensionless coefficients α n , which are ratios of α n to the appropriate powers of the vacuum quark condensate qq . Solid and dashed lines, enclosing the shaded areas, correspond to µ = 0 and 300 MeV, respectively.
Except for α 2 (upper panel), we see that the larger the temperature, the smaller are the magnitudes of the α n . Furthermore, the sign of the coefficients is strongly µ dependent. The energy density as a function of the order parameter for given values of T and µ is defined by the magnitude of the GL coefficients. Moreover, the analysis of the relative sign between those coefficients determines the different regions in the QCD phase diagram. In particular, the Lifshitz point (LP) is defined as the point where the inhomogeneous phase and the two homogeneous phases with broken and restored chiral symmetry meet, while the tricritical point (TCP) denotes the point where the second-order chiral phase transition turns into a first order one. The positions of the LP and the TCP can be determined by solving, as a function of T and µ, the following set of equations [18,51] For the set of parameters of the nonlocal NJL model given above, the coordinates in the (T, µ) plane for the LP and the TCP (in MeV) are ( Fig. 2 of Ref. [52]. In the following we explore the dynamics of the order parameters for values of T and µ close to the LP and TCP in the equilibrium phase diagram, namely (T, µ) = (50 MeV,190 MeV). At this point of the phase diagram, the σ(x) = σ ∼ 250 MeV and π(x) = 0.   Eqs. (9) and (10) are solved numerically using a finitedifference method. We define the dimensionless time coordinate τ = t/Γ and discretize it in time steps δτ = 10 −5 . The space coordinate is discretized into a onedimensional lattice with lattice spacing a = 0.1 fm. Since we do not have a microscopic derivation of Eq. (4), the value of Γ must be taken from independent sources. A good source are calculations using the influence functional in the linear sigma model [58][59][60], where Γ σ and Γ π as a function of T have been provided. For low values of T , within the range 0 ≤ T ≤ 50 MeV, which is the one of interest here, Ref. [58] finds in the chiral limit Γ σ 3.75 fm and Γ π 0 (at lowest order the coupling, Γ π = 0 at all temperatures in the chiral limit). The reason for the difference between the coefficients is due to the fact that the process σ → 2π is allowed for all temperatures (including T = 0), while the reverse process and πσ → π are strongly suppressed due to the large mass of the σ compared to the one of the π [58]. There is, however, one difficulty here, in that there have been no estimates of Γ as a function of T and µ. In the absence of such estimates, we explore the consequences of using either Γ π = Γ σ = 0 and Γ σ = 0 and Γ π = 0.
In the top panel of Fig. 2 we display the volume average of the scalar field σ(τ ): where x n = n a and N is the number of lattice sites. Results are shown for T = 50 MeV and three different values of chemical potential: µ = 180, 200 and 220 MeV. Whereas in the remaining panels we plot snapshots of the TDGL time evolution at τ = 0, 10, 10 3 , 2 × 10 5 for the scalar field σ(x, τ ). The initial profile for σ(x, 0) and π(x, 0) was set by imposing a unbiased white Gaussian noise for each position on the lattice, simulating a situation of an out-of-equilibrium state that has been quenched to a low-temperature phase. Here, and up to Fig. 5, we use Γ π = Γ σ . We also mention that the pseudoscalar field π(x, τ ) (not shown in the figure), after passing through inhomogeneous configurations, evolves to π(x, τ eq ) = 0, as it should. In addition, it is clear that chiral symmetry is restored for µ 220 MeV, which is the correct equilibrium state for this value of temperature.
Next, we investigate the effect of the initial configuration on the equilibration time. Fig. 3 shows results for initial profiles for the scalar and pseudoscalar fields given by dual chiral density waves (DCDW) superimposed with Gaussian noise. This kind of modulation is one of the few one-dimensional spatial dependences that can be derived analytically from the GL mean field equations [18]. As expected, for the point (T, µ) = (50 MeV, 190 MeV) of the phase diagram, the equilibrium profile is a homogeneous symmetry-broken configuration. One sees that the time to reach the homogeneous equilibrium configuration is longer than in the previous case. Here, one also sees an interesting feature of the volume average: for this initial configuration the system stays longer in metastable regions at early times than for the purely random configuration. It should be clear that if the homogeneous symmetry-broken solution would be reached without passing through intermediate inhomogeneous phases, the sigma field σ(x, τ ) would grow uniformly throughout the volume and, therefore, the volume average σ(τ ) would present a monotonic behavior. C learly, this is not the case. This oscillating behavior of the volume average signals inhomogeneities during the TDGL dynamics, compatible with a DCDW. For completeness, we show  in Fig. 4 the evolution of the pseudoscalar field when the initial condition is a DCDW. Before equilibrium, this field also goes through inhomogeneous metastable phases and reaches its equilibrium value π(x, t eq ) = 0 for long times. It should be noted that the comparable equilibration time for both fields is due to the fact that we are using Γ π = Γ σ here.
Similar metastable configurations appear in the course of the evolution for different initial conditions. The form of the intermediate-state configurations reflect the initial profile. This is emphasized in Fig. 5, where we show snapshots of inhomogeneous configurations of σ(x, τ ) for an initial antisymmetric random profile. The top panel of the figure reveals a hyperbolic tangent profile at long times, before reaching the expected homogeneous configuration (at τ ∼ 10 5 ). Next we show the effect of taking Γ σ = 0 and Γ π = 0: the bottom panel of the figure shows that the equilibrium configuration has a hyperbolic tangent shape-increasing the simulation time does not change this shape. This is not the expected equilibrium configuration, as discussed previously. The point is that  in this case of decoupled order parameters, and for this particular initial configuration, the σ field is driven into a local minimum of the energy functional. The coupling with the π field is essential to get the σ out of the local minimum. Of course, fluctuations, even at such low temperatures, might play a role here and change the picture.
Since the thermalization time depends on the size of the system, it is important to obtain a scaling law of τ eq with the length of the lattice, L = N a, where N is the number of lattice points. For all the three initial configurations employed along the work, we found that τ eq grows almost quadratically with L. A good fit to the data from simulations using several values of L is obtained by the formula τ eq = A L B , with A = 3.45 (69) and B = 2.154 (38), with a coefficient of determination of R 2 = 0.993. With such a scaling law, one can make a rough comparison of τ eq with time scales in typical heavyion collision. We emphasize that such a comparison is far from rigorous, and it might even not be entirely appropriate given the setting of the study: one-dimensional lattice, no expansion of the system and therefore fixed temperature, etc. Nevertheless, for orientation, but keeping this proviso in perspective, if one takes Γ σ = 3.75 fm [58], well-defined equilibrium configurations are reached after t eq ∼ 10 4 fm, for L = 20. That is, τ eq is three orders of magnitude larger than, say the decoupling time extracted in Pb-Pb collisions at the LHC [74], for which the kinetic freeze-out volume is L 3 ∼ 5 × 10 3 fm 3 . Now, it is important to mention that (T, µ) = (50 MeV,190 MeV) corresponds to a baryon density of the order of the nuclear matter density ρ 0 = 0.17 fm −3 . In a heavy-ion collision at FAIR, for example, after a time interval of the order of 15 fm after the collision, the density of the produced matter will be close to ρ 0 [25]. Therefore, if long-lived chiral inhomogeneities are produced from those high-density regions in such heavy-ion collisions, they should leave traces in observables, like e.g. in pion number fluctuations. As already mentioned, heavy-ion collisions producing large baryon densities at low temperatures are envisaged at NICA [26] and, more in the future, at J-PARC [27].

IV. SUMMARY AND CONCLUSIONS
In this work we have analyzed the dynamics of formation of inhomogeneous metastable chiral structures in the symmetry-broken phase of the phase diagram of strongly interacting matter as predicted by a chiral quark model. More specifically, we have employed a time-dependent Ginzburg-Landau (TDGL) equation to describe the dynamics of the scalar σ and pseudoscalar π chiral order parameters near the tricritical point (TCP) and the Lifshitz point (LP) of the equilibrium phase diagram. The former denotes the location in the phase diagram where a second order chiral phase transition turns into a first order transition, while the latter determines where one inhomogeneous phase and two homogeneous phases with broken and restored chiral symmetry meet. We have used the simplest nonlocal extension of a chiral SU(2) Nambu-Jona-Lasinio model. To the best of our knowledge, this work constitutes the first investigation of the dynamics of formation of inhomogeneous chiral phases QCD matter at low temperature and baryon density close to the saturation density of nuclear matter.
The solutions of the TDGL equations revealed the presence of long-lived metastable configurations of the order parameters in the course of the evolution. The time dynamic was studied in a region of the phase diagram where the equilibrium configurations of the σ and π order parameters have no spatial modulations. Initially, we verified our approach in regard to chiral restoration at low temperature as function of the chemical potential µ in the vicinity of the critical points. We verified that the TDGL equation leads an equilibrium solution for which the symmetry is restored for values of µ in good agreement with those obtained in Ref. [52]. We also found that when π is set equal to zero during the entire evolution, the σ field is not driven to the correct equilibrium configuration; it is driven into a local minimum of the energy functional. The coupling with the π field is essential to drive σ to the correct equilibrium configuration.
To finalize, the main conclusion of the paper is that inhomogeneous configurations of the chiral order parameters produced in the course of the evolution of matter can be long-lived, with lifetimes much larger than the typical lifetimes in a heavy-ion collision. This means that, even if at equilibrium the chiral parameters have no spatial modulation, the system decouples before reaching such a state and can leave traces of the inhomogeneities in observables. The investigation of such possible traces in observable is left for a future study, when several extensions of the present study will be made. First of all, it is underway the extension of the study to three-dimensions and inclusion of expansion of the system [75]. Finally, the derivation of a TDGL equation from the nonlocal NJL model used in the present work would be essential to avoid uncertainties regarding the kinetic coefficients.