Dyson-Schwinger study of chiral density waves in QCD

The formation of inhomogeneous chiral condensates in QCD matter at nonzero density and temperature is investigated for the first time with Dyson-Schwinger equations. We consider two massless quark flavors in a so-called chiral density wave, where scalar and pseudoscalar quark condensates vary sinusoidally along one spatial dimension. We find that the inhomogeneous region covers the major part of the spinodal region of the first-order phase transition which is present when the analysis is restricted to homogeneous phases. The triple point where the inhomogeneous phase meets the homogeneous phases with broken and restored chiral symmetry, respectively, coincides, within numerical accuracy, with the critical point of the homogeneous calculation. At zero temperature, the inhomogeneous phase seems to extend to arbitrarily high chemical potentials, as long as pairing effects are not taken into account.

The properties of strong-interaction matter under extreme conditions, such as high temperature or density, and the corresponding phase structure of Quantum Chromodynamics (QCD) are subject of extensive theoretical and experimental investigations [1,2]. At vanishing net baryon density, first-principle lattice gauge calculations have revealed that the approximate chiral symmetry of QCD, which is spontaneously broken at low temperature, gets restored at high temperature in a cross-over transition [3,4]. At low temperature and high baryon density, where lattice calculations are inhibited by the sign problem, effective model studies typically predict that chiral symmetry is restored in a first-order phase transition, which weakens with increasing temperature and eventually ends at a critical point [5]. More recently, this picture was confirmed with Dyson-Schwinger equations (DSEs) applied to QCD [6,7,8,9].
A basic assumption in these investigations was that the phases are homogeneous, i.e., in particular, the chiral order parameter is constant in space. On the other hand, phases with non-uniform chiral order parameters have been proposed already long time ago, see Ref. [10] for a brief historical review. Starting with Migdal's pwave pion condensation [11], the idea was generalized to relativistic systems [12,13,14] and studied in high-density QCD with large number of colors, applying weak-coupling methods [15,16]. More recently, it gained new attention after it was found in effective models that the firstorder chiral phase boundary between homogeneous phases is covered completely by an inhomogeneous phase [17,18].
The aim of our work is to study the chiral phase transition with the possibility of inhomogeneous condensates in strong-coupling QCD with DSEs. We restrict ourselves to a so-called chiral density wave (CDW), where the chiral condensate rotates along the chiral circle, when moving into a fixed direction. Specifically, the scalar and pseudoscalar condensates behave like qq ∝ cos(Qz) , with Q being the modulus of a wave vector, which we have chosen to point into the z direction. Moreover, as indicated by the Pauli matrix τ 3 , we have chosen the third isospin component of the pseudoscalar condensate. Assuming isospin invariance of the QCD Lagrangian, this can be done without loss of generality. The DSE for the dressed quark propagator S is diagrammatically depicted in Fig. 1. In coordinate space with Euclidean metric, it is given by , depending on two space-time variables, x and x ′ . S 0 denotes the bare propagator, Σ the selfenergy and Z 2 is the wave-function renormalization constant of the quark field.
In homogeneous, i.e., translationally invariant matter, the propagator depends only on the relative coordinate x − x ′ . In momentum space, this translates into a dependence on a single 4-momentum, and the inverse propagator at temperature T and chemical potential µ can be parametrized as with p := ( p, p 4 = ω n + iµ), the Matsubara frequencies ω n = (2n + 1)πT , and three dressing functions A, B, and C. In vacuum, due to Lorentz covariance, these functions depend on p 2 only and A(p) = C(p). The wave-function renormalization constant Z 2 is then fixed by the condition that A(p)| p 2 =ν 2 = 1 at an arbitrary renormalization point ν. This prescription remains valid for our analysis of inhomogeneous phases since Z 2 is always fixed in vacuum, which is homogeneous. In inhomogeneous matter, the quark selfenergy and, hence, the propagator depend separately on both coordinates x and x ′ , or, equivalently, on the relative coordinate x − x ′ and the center of momentum coordinate (x + x ′ )/2. In momentum space they thus depend on two momenta, and the DSE reads Here p and p ′ correspond to the out-and ingoing momenta of the quark, which do not need to be identical. Physically, this means that the inhomogeneous condensates carry momentum, so that the quark can change its momentum by scattering off the condensate. The (inverse) propagator can be viewed as a continuous matrix in momentum space. For the formal manipulations to be discussed below it is useful to introduce a finite quantization volume V in 3-space and take periodic boundary conditions. Thus, together with the finite extent in the imaginary time direction, the 4-volume is finite as well, V = V × [0, 1/T ], and the 4-momenta take discrete values. At the end we will take the limit V → ∞, so that the 3momenta will be continuous variables again. Momentum sums and Kronecker symbols should then be replaced as where n and n ′ label the Matsubara frequencies.
In this Letter we consider two quark flavors with vanishing bare mass. Since the bare propagator is constructed on the (homogeneous) perturbative ground state, it stays diagonal in momentum space and keeps its familiar form, The quark selfenergy is given by (see Fig. 1) with the QCD coupling constant g, the bare and dressed quark-gluon vertices, gΓ a,0 µ and gΓ a ν , respectively, and the dressed gluon propagator D µν . As detailed below, we neglect possible modifications of these quantities with respect to the homogeneous case. As a consequence, the gluon propagator depends only on a single momentum variable, and 4-momentum is conserved at the vertices, i.e., k = p − q and q ′ = q + p ′ − p.
The bare vertex is given by Γ a,0 µ = Z 1F γ µ λ a /2 with a renormalization constant Z 1F and the Gell-Mann matrix λ a . The dressed gluon propagator and the dressed vertex are in principle given by their own DSEs. Since these depend on even higher n-point functions, truncations are necessary to get a closed set of equations. Here we adopt the truncation scheme described in Ref. [19], and we refer to that reference for details and parameters.
In this scheme the dressed vertex is taken to have the same structure as the bare vertex, Γ a µ (p, q) = Γ(p − q)γ µ λ a /2 with a dressing function Γ(k), which has the correct perturbative running in the ultraviolet and a phenomenological enhancement in the infrared. The gluon propagator is based on a parametrization of lattice data for the Yang-Mills system, which is corrected for quark effects by perturbatively adding a polarization loop in hard-thermal-loop-hard-dense-loop approximation. This accounts for Debye screening and Landau damping at high temperature or chemical potential, but neglects the dressing of the quarks in the polarization loops. As a consequence, the dressed gluon propagator remains diagonal in momentum space, as already mentioned above.
The task is now to generalize the structure Eq. (2) of the quark propagator in a homogeneous ground state to an inhomogeneous medium where the quark condensate takes the form of a CDW, Eq. (1). Starting from the definition Z 2 S(x, x ′ ) = T (q(x)q(x ′ )) of the Euclidean propagator, where T is denotes the imaginary time ordering operator, the condensates are related to the propagator as with the trace in Dirac, flavor and color space. Comparing this with Eq. (1), we find that the desired spatial behavior is obtained if with the wave vector Q ≡ Qe 3 being a 4-vector of length Q, pointing to the 3-direction. This suggests to generalize the dressing function B , which is the only chiral-symmetry breaking term in Eq. (2), in a similar way. Specifically, we make the ansatz that the inverse propagator S −1 (p, p ′ ) contains a term where B(p) is closely related to the B function of the homogeneous case. In fact, when we take Q = 0, the matrix B(p, p ′ ) becomes purely scalar and diagonal in momentum space, with the diagonal matrix elements essentially given by B(p). 1 It is also instructive to compare our ansatz with the Nambu-Jona-Lasinio model where the quark selfenergy is local in coordinate space. For homogeneous matter, this leads to a constant selfenergy in momentum space, while for inhomogeneous matter the selfenergy can only depend on the difference p − p ′ but not on the sum. For a CDW the selfenergy is then given by Eq. (9) with B(p) = const. In our case, B(p) is an unknown function, which must be determined through the DSE. To that end, the inverse propagator with the dressing function Eq. (9) must be inverted and inserted into Eq. (6). It turns out that this induces further structures, and we need additional dressing functions to achieve a self-consistent solution. For instance, since the wave vector Q = Qe 3 defines a preferred direction, the A function in Eq. (2) must be replaced by two independent functions, corresponding to the momentum component p 3 in Q direction and to the perpendicular part p ⊥ ≡ p 1 e 1 + p 2 e 2 . The complete ansatz contains in total 10 dressing functions and reads withB(p, p ′ ) = 1 2 (B(p) + B(p ′ )) and similar forF ,Ḡ and H.
For a general inhomogeneous ansatz, the inversion of S −1 , which is an infinite matrix in momentum space, is a highly non-trivial task. However, for the CDW it turns out that S −1 has a relatively simple block structure so that it can be inverted analytically. The resulting dressed propagator has the same tensor structure as S −1 .
Starting from Eq. (10) one can find self-consistent solutions of the DSE for arbitrary values of the wave number Q. We thus need an additional constraint to fix Q. This is provided by the requirement that the free energy of the system is minimal for the stable solution or, equivalently, the pressure is maximal. The latter corresponds to the effective action where the traces are over momentum, color, flavor and Dirac components. The last term denotes the two-particle irreducible interaction part. In our truncation scheme it corresponds to the diagram shown in Fig. 2 and is given by Tr Γ a,0 µ S(q, q ′ )Γ a ν (q ′ , p ′ )D µν (k)S(p ′ , p) , where we have written the momentum sums explicitly, so that the trace is only over internal degrees of freedom. The variation of the effective action with respect to the dressed propagator, We now take the infinite-volume limit as specified in Eq. (4) and solve Eqs. (3) and (13) numerically. Results for the mass amplitude M (0) ≡ B( 0, n = 0)/C( 0, n = 0) and the wave number Q are presented in Fig. 3. In the upper panel we show them for T = 0 as functions of µ. At low chemical potential, we only find a homogeneous solution, whereas above µ = 350 MeV, there is also an inhomogeneous solution. Comparing the pressure, we find that the inhomogeneous phase becomes favored above a critical chemical potential of about 365 MeV. At this point a firstorder phase transition takes place, where the wave number jumps from zero to a finite value, while the mass amplitude drops discontinuously, but remains nonzero. When we increase the chemical potential further, the latter continues to decrease, but with decreasing slope, suggesting that the inhomogeneous phase survives up to arbitrarily high chemical potentials. In fact, within numerical accuracy, we never find the chirally restored phase to be favored at T = 0.
The lower panel shows the same quantities at µ = 370 MeV as functions of temperature. At this chemical potential the system is inhomogeneous at low temperatures, as evident from the nonvanishing Q. With increasing temperature, Q increases further, while the mass continuously Collecting the results from different temperatures and chemical potentials we obtain the phase diagram displayed in Fig. 4. At high temperature there is a second-order phase transition between the homogeneous chirally broken phase and the chirally restored phase (purple dash-dotted line), which becomes first order at low temperatures (green dotted line) when the analysis is restricted to homogeneous phases. The corresponding spinodals are indicated by the red dashed lines.
The limits of the region where we find an inhomogeneous solution with our CDW ansatz are marked by the blue solid lines. It can be deduced from the behavior of dΓ dQ (Q) that in this region the inhomogeneous phase is always favored over the restored phase. The restored phase is eventually reached in a second-order phase transition. Moreover, from the fact that the homogeneous chirally broken phase is energetically degenerate with the restored phase along the green dotted line, we conclude that the phase transition from the homogeneous to the inhomogeneous chirally broken phase must be to the left of this line. As already seen for T = 0, this phase transition is first order. Hence, the phase boundary must be somewhere between the left solid and the dotted line. To locate it more precisely, we have to compare the pressure of the two solutions, which is numerically quite demanding. Within numerical precision we find that at T = 0 the critical chemical potential is about 10 MeV lower than in the homogeneous case.
Our most important result is that, within numerical resolution, the inhomogeneous phase covers the firstorder phase boundary of the homogeneous case completely. Moreover, the point where the inhomogeneous phase and the two homogeneous phases meet seems to coincide with the homogeneous critical point. In this respect, Fig. 4 has great similarities with the phase diagram in the Nambu-Jona-Lasinio (NJL) model [18]. A qualitative difference is that in the NJL model the inhomogeneous phase ends at some upper critical chemical potential, whereas in the present case it seems to extent to arbitrarily high µ at T = 0.
On the other hand, at high densities inhomogeneous chiral symmetry breaking should become disfavored against quark pairing (color superconductivity) [16], which we have neglected here. We have recently studied color superconductivity in a similar framework [19] and it should be a feasible task to extend the present analysis in this direction. Additionally it needs to be checked if the results of this work are robust under the improvement of the truncation. The consideration of more complicated inhomogeneous structures than the CDW ansatz would be interesting but extremely difficult as the propagator can no longer be inverted analytically.