Inhomogeneous phases in the quark-meson model with explicit chiral-symmetry breaking

We investigate the existence of inhomogeneous chiral phases in the quark-meson model with explicit chiral-symmetry breaking. We find that the inhomogeneous region shrinks with increasing pion masses but survives for the physical value of m_pi. The instability towards inhomogeneous matter occurs in the scalar channel, while pseudoscalar modes are disfavored.


Introduction
Mapping the phase diagram of QCD at nonvanishing temperature T and quark chemical potential µ is one of the major challenges in strong-interaction physics [1,2]. Lattice QCD calculations at µ = 0 revealed that chiral symmetry, which is spontaneously broken in vacuum, gets approximately restored in a smooth crossover at T ∼ 150 MeV [3], while in the regime of low T and nonzero µ, where standard lattice methods are not applicable, model studies as well as continuum approaches to QCD indicate the possibility of a first-order phase transition, terminating at a second-order critical endpoint (CEP) [4,5,6,7]. Most of these calculations rely however on the assumption that these phases are homogeneous, i.e., that the chiral order parameter does not vary in space. Allowing for spatially non-uniform order parameters, such inhomogeneous phases often turn out to be favored in some region of the phase diagram, typically covering parts of or even the entire first-order boundary between the homogeneous phases. Specifically this was found within the Nambu-Jona-Lasinio (NJL) [8,9] and the Quark-Meson (QM) model [9,10], but also in QCD using Dyson-Schwinger equations [11]. (For a review about inhomogeneous chiral phases, see Ref. [12].) However, most of these studies have been performed in the chiral limit, while the situation for the more realistic case with a small explicit breaking of chiral symmetry is less clear. For the NJL model it was found that the inhomogeneous phase shrinks when a nonvanishing bare quark mass is introduced but is still present for realistic masses [9]. More generally it was shown in Ref. [13] that the inhomogeneous phase always reaches up to the CEP in that model and thus survives as long as there is a first-order phase transition in the homogeneous case. For the QM model, on the other hand, it was found in Ref. [14] that the inhomogeneous phase becomes disfavored a e-mail: michael.buballa@tu-darmstadt.de arXiv:2006.02133v1 [hep-ph] 3 Jun 2020 already for a rather small amount of explicit symmetry breaking, corresponding to a pion mass of about one quarter of the physical value. The calculation was however done only for one specific spatial modulation of the order parameter, a so-called chiral density wave (CDW). This modulation is relatively simple to handle but is known not to be the most favored shape in most cases, and away from the chiral limit it is not even a self-consistent solution.
In the present work we therefore study the effect of explicit chiral-symmetry breaking on inhomogeneous phases in the QM model, starting with a stability analysis of the homogeneous phase. This method, which has already been employed in Ref. [13] to the analogous problem in the NJL model, has the advantage that one does not need to know the explicit shape of the spatial modulation. It only relies on the assumption that at the phase boundary the homogeneous phase becomes unstable against small inhomogeneous fluctuations, i.e., that the phase transition is of second order. This analysis will therefore yield a sufficient criterion for the inhomogeneous regime in the model, while the true inhomogeneous phase can be larger. We will then support the results of the stability analysis with a calculation of the full model phase diagram employing a specific self-consistent ansatz away from the chiral limit, as well as by a Ginzburg-Landau expansion close to the CEP.
The remainder of this article is organized as follows: We introduce our theoretical framework in Sec. 2, then in Sec. 3 we discuss the QM parameter fitting procedure away from the chiral limit. We show our numerical results for the phase diagram in Sec. 4 and discuss our conclusions in Sec. 5.

Theoretical framework
We consider the QM model defined by the Lagrangian where ψ is a quark spinor field with N f = 2 flavor and N c = 3 color degrees of freedom, coupled via a Yukawa interaction with coupling constant g to the scalar sigma meson σ and the pseudoscalar pion triplet π. Here τ = (τ 1 , τ 2 , τ 3 ) denotes the Pauli matrices in isospin space. The meson kinetic contributions read and is the meson potential. In the limit c = 0 it is symmetric unter O(4) transformations of the meson vector φ = (σ, τ ), which can be identified with the chiral SU (2) L × SU (2) R symmetry. For c = 0 the O(4) symmetry is broken explicitly down to O(3), corresponding to the SU (2) isospin symmetry. The model parameters, g, λ, v 2 , and c, will be fitted to vacuum properties as discussed in Sec. 3. The thermodynamic properties of the model are encoded in the grand potential per volume V , Ω(T, µ) = − T V log Z(T, µ), where Z(T, µ) denotes the grand canonical partition function, which depends on the temperature T and the quark chemical potential µ. In the following we perform the mean-field approximation, replacing the quantum fields σ and π by their expectation values, i.e., by classical fields. We assume that these fields are time independent but we retain their dependence on the spatial coordinate x in order to allow for inhomogeneous phases. Moreover, we assume that only the third isospin component of the pion field develops a nonvanishing expectation value, which we call π(x), while, for simplicity, we keep the name σ(x) for the classical sigma field. The mean-field grand potential per volume ("thermodynamic potential") is then given by with a purely mesonic part and a quark part where corresponds to the inverse dressed quark propagator at chemical potential µ in the presence of the sigma and pion mean fields, and Tr denotes a functional trace running over the Euclidean four volume V 4 = [0, 1 T ] × V as well as color, flavor and spinor degrees of freedom. These expressions are basically identical to those in Ref. [10], where the same model in the chiral limit was considered. The only exception is that we now have to take into account the explicitly symmetry-breaking term −cσ in the mesonic potential.

Stability analysis
In order to determine the ground state of the system at given T and µ, we must minimize the thermodynamic potential with respect to the mesonic fields σ and π. While this is standard for spatially constant mean fields, the functional minimization of Ω MFA with respect to arbitrary non-uniform fields is obviously a much harder problem, which has not yet been solved in full glory for 3 + 1 space-time dimensions. Instead of tackling the full problem, one possibility is to perform a stability analysis, applying the same methods, which have been used in Ref. [13] to analyze inhomogeneous phases in the NJL model. To this end we split the meson fields into spatially constant parts, corresponding to the lowest homogeneous state of the system, and small fluctuations with arbitrary spatial shapes. Since in homogeneous systems the pion field is disfavored against the sigma field due to the symmetry-breaking term in the potential, the constant part appears only in the sigma sector, i.e., we have whereσ corresponds to the (in general T and µ dependent) value of the sigma field in the homogenous ground state, and δσ and δπ are the fluctuations. Plugging this into Eqs. (4) -(7), the thermodynamic potential can be decomposed as with Ω (n) being of the nth order in the fluctuations. Specifically one finds for the contributions up to quadratic order where is the inverse quark propagator Eq. (7) without fluctuations, S 0 is its inverse, and is the quark selfenergy correction due to the fluctuating fields.
Noting that S 0 corresponds to the propagator of a free fermion with mass these expressions are evaluated most easily in momentum space. Assuming spatially periodic fields we perform the Fourier expansions where q k are the elements of the corresponding reciprocal lattice. Since the meson fields and, thus, their fluctuations are real fields in coordinate space, the Fourier coefficients obey the relations δσ −q k = δσ * q k and δπ −q k = δπ * q k . Taking the infinite-volume limit V → ∞ one then obtains for the linear contribution of the fluctuations to the thermodynamic potential. Here we have introduced the loop integral F 1 , where withM as defined in Eq. (15) and fermionic Matsubara frequencies ν m = (2m+1)πT . Note that only the spatially constant q k = 0 mode of the fluctuations in the sigma channel contributes to Ω (1) . However, since we have assumed that Ω (0) corresponds to the lowest homogeneous state, this contribution must vanish, leading to the gap equation Indeed, the same equation can be obtained from the stationary condition dΩ (0) dσ = 0.
Unlike the linear term, the quadratic corrections of the fluctuations to the thermodynamic potential get contributions from all Fourier modes. One finds where q k = (0, q k ) is the four-momentum vector with vanishing energy and threemomentum q k , and are the (unrenormalized) inverse dressed meson propagators at four-momentum q, temperature T and chemical potential µ. Here are the sigma and pion tree-level masses, and Π M (q) denote the corresponding quarkantiquark polarization loops (cf. Ref. [10] for further details). The explicit evaluation yields where we have used the gap equation (19) to eliminate terms proportional to the loop function F 1 . Taking again the infinite-volume limit, the loop function L 2 is given by where ν n are again fermionic Matsubara frequencies and ω m is a bosonic Matsubara frequency. As pointed out above, we only need L 2 at zero energy at this point, i.e., ω m = 0. In the infinite-volume limit the crystal can take any geometry and size, and therefore the momenta q k of the reciprocal lattice are not a priori restricted to certain values. As can be seen from Eq. (20), the free energy of the homogeneous ground state can thus be lowered by the formation of inhomogeneous modes if D −1 σ (q) > 0 or D −1 π (q) > 0 in some region of q = (0, q). Note that q 2 = −q 2 in this case, so that the inverse propagator of a free meson, D −1 M,free = q 2 − m 2 M is always negative. The instability is therefore a pure interaction effect, as also known, e.g., from P-wave pion condensation in nuclear matter [15] (see [16] for a review). In the present model one can distinguish between meson-meson interactions (encoded in the tree-level masses m M,t ) and quark-meson interactions (encoded in the polarization functions Π M ). The latter are identical to the polarization functions in the NJL model, and it was shown in Ref. [13] that they favor an instability in the sigma channel over an instability in the pion channel. 1 In the QM model, however, the situation is more complicated because of the tree-level masses. At least, if we naively assume the ordering m 2 σ,t > m 2 π,t > 0, as for the physical masses in vacuum, we would expect that both masses stabilize the homogeneous phase but less in the pion channel than in the sigma channel. In order to find out the overall effect we therefore have to evaluate Eqs. (23) and (24) explicitly. The resulting stability boundaries of the homogeneous phase will be presented in Sec. 4.1.

The Real-Kink Crystal ansatz
The stability analysis described above has the clear advantage to provide general results independent of the specific shape of the spatial modulation of the order parameter. However, since it relies on a small-amplitude expansion, it can only provide a sufficient condition for an inhomogeneous phase, while the true inhomogeneous region can be larger. Thus, in order to complement its results and obtain an estimate of the size of the inhomogeneous window, we will also compute the full thermodynamic potential of the QM model for a specific ansatz for the order parameter, the socalled "real-kink crystal" (RKC). Aside from the advantage of being a self-consistent ansatz away from the chiral limit [9], this RKC is also the energetically most favored modulation considered so far in the literature [12,17,18].
The order parameter is expressed in terms of the Jacobi elliptic functions sn, cn, dn, and is characterized by the three variational parameters ∆, ν and b, which are determined by minimizing the free energy of the system [9]. For this type of modulation, an analytical expression for the density of states ρ(E) has been computed, so that the free energy can be obtained without having to resort to numerical diagonalization of the inverse quark propagator [9,12]. One finds where = E/∆,ν = 1−ν,λ = arcsin( / √ν ), λ = arcsin(1/ ), C(ν) = E(ν)/K(ν)−1, F and K are incomplete and complete elliptic integrals of first kind, respectively, and E are the complete and incomplete elliptic integrals of second kind.
The thermodynamic potential is then given by with δ = 1/sn 2 (b, ν) − 1 and where the first term corresponds to the vacuum quark contribution.
The meson potential depends on the following spatial averages of the order parameter over a period: where Z is the Jacobi Zeta function. Before getting to our results for the model phase structure, let us now discuss how the model parameters are fixed.

Parameter fixing
As standard, we determine the model parameters by fitting masses and the pion decay constant in vacuum. Thereby, in order to systematically investigate the effect of the explicit chiral-symmetry breaking, we first set the coupling c equal to zero and fix the remaining parameters g, λ and v 2 in the chiral limit. After that, we consider c = 0 but keep the other parameters at their chiral-limit values.
For fixing g, λ and v 2 in the chiral limit we follow Ref. [10], where this was done by fitting the vacuum values of the pion decay constant f π , of the sigma-meson mass m σ , and of the constituent quark mass. For homogeneous matter we can identify the latter withM as defined in Eq. (15) withσ being the homogeneous sigma field which minimizes Ω (0) . In vacuum, i.e., at T = µ = 0, we expect that it also minimizes Ω MFA , since phenomenologically the vacuum is homogeneous. This turns out to be true in our model as well, at least up to quadratic-order fluctuations. For m σ and f π it was shown in Ref. [19] that it is crucial to fit the pole mass and to take into account the renormalization of the pion wave function, corresponding to the pole of D σ and the residue of D π , respectively. The resulting expressions are (see Refs. [10,19] for details) where the subscript 0 inM 0 , m σ,0 and f π,0 indicates that these quantities correspond to the vacuum values in the chiral limit. Likewise F (m σ,0 ) means that the function L 2 ((iω m , q)) is analytically continued to the real time-like momentum q = (m σ,0 , 0). The explicit expressions can be found in Ref. [19]. Having fixed g, λ and v 2 in this way, we turn on the chiral-symmetry breaking term by varying the parameter c. The most important consequence is that the pion, which is massless in the chiral limit in agreement with the Goldstone theorem, gets a non-vanishing mass. We can therefore relate the parameter c to the pion pole mass, implicitly given by D −1 π (q = (m π , 0)) = 0. We then get from Eq. (24) where L (vac) 2 (m π ) is the function L 2 ((iω m , q)) evaluated in vacuum and analytically continued to the real time-like momentum q = (m π , 0). Note, however, that the quark massM , which also enters the function L (vac) 2 , is not the chiral-limit valueM 0 , as in Eqs. (34) -(36), but related to the solution of the gap equation (19), including the constant c. For a fixed value of m π , Eq. (37) must therefore be solved self-consistently together with Eq. (19).
Finally, we note that the vacuum parts of loop integrals F 1 and L 2 , as well as their chiral-limit versions, are ultraviolet divergent and must be regularized in order to get meaningful results. 2 Again following Refs. [10,19], we use Pauli-Villars regularization with three regulators, controlled by the cutoff parameter Λ. As a consequence, the model parameters for fixed values ofM 0 , m σ,0 , f π,0 , and m π depend on Λ.
In the following, we will always fix our model in the chiral limit by choosingM 0 = 300 MeV, m σ,0 = 600 MeV, and f π,0 = 88 MeV. In particular we have m σ,0 = 2M 0 , in which case Eq. (35) simplifies to λ = 2g 2 . The corresponding values of λ and v 2 as functions of Λ are displayed in the first two panels of Fig. 1. The results agree with those in Refs. [10,19], where the same vacuum observables have been fitted. In addition, we show in Fig. 1 the parameter c, multiplied with g, for m π = 140 MeV.
Since the QM model is renormalizable, all observables should remain finite in the limit Λ → ∞. As demonstrated in Ref. [10], this is also true for the phase diagram. It was found that the results remain practically unchanged when Λ exceeds 2 GeV, so that in practice Λ = 5 GeV can be considered as the "renormalized limit". 3 In Fig renormalized limit. By construction, they take of course their fit values in the chiral limit, i.e., at m π = 0. With increasing m π they increase as well but stay finite, even for arbitrarily large values of Λ. We note that the value of f π for the physical pion mass m π ∼ 140 MeV is too small compared with the empirical value of 92.2 MeV [22]. This could be cured by slightly changing the fit values in the chiral limit (which are admittedly somewhat ad-hoc) but it is not our intention here to perform a precision fit. Moreover, in Fig. 2, f π has been calculated as [10] which corresponds to the quark-level Goldberger-Treiman relation and is strictly speaking only valid in the chiral limit (cf. Eq. (34)).
A more severe problem is that g 2 and λ diverge at the point when L (vac,0) 2 = −2f 2 π,0 /M 2 0 . Within our regularization scheme and for our parameters this happens at Λ = Λ * ≈ 757 MeV. Beyond this point, g 2 and λ even turn negative, see Fig. 1, and, related to this, Ω MFA is no longer bounded from below in this regime [10,20]. Moreover, a negative g 2 obviously means that the Yukawa coupling g is imaginary. Hence, forcing the constituent quark mass, Eq. (15), to stay real, the field expectation valueσ becomes imaginary as well, in contradiction to our original assumption of σ and π being real fields. Although it has been argued in Ref. [20] that the unbounded potential is a known one-loop artifact and should be cured at higher orders, this is clearly worrisome. On the other hand, the phase diagram changes smoothly when passing through Λ = Λ * , i.e., focusing only on the phase diagram, one would not even notice that the problem exists. In Sec. 4 we will therefore discuss results for the renormalized limit, ignoring the inconsistencies, as well as for Λ = Λ * , being the largest possible cutoff outside the problematic regime.

Phase structure
We are now ready to discuss our results for the model phase structure, starting from the stability analysis to determine the boundary where inhomogeneous phases become favored.

Stability analysis
In Fig. 3 we show the stability boundaries of the homogeneous phase with respect to inhomogeneous fluctuations. More precisely, we show the lines where D −1 M (q) just touches the zero-line at some value of q = (0, q = 0), both for M = σ or M = π, for different values of the PV regulator and the vacuum pion mass. We recall that this type of analysis relies on the assumption that the spatially modulated order parameters are small, and thus can only give reliable results for second-order phase boundaries. According to explicit calculations with certain modulations, this is typically the case at the right phase boundary of the inhomogeneous region, while the left boundary cannot reliably be determined by the stability analysis. We will confirm this below in Sec. 4.2.
As demonstrated in Ref. [10] for the chiral limit, incorporating vacuum fluctuations shrinks the size of the inhomogeneous phase, which nevertheless survives in the renormalized limit. As we move away from the chiral limit, the stability lines in the two channels split, with the sigma line becoming the only relevant one since it is the first to appear when coming from the stable homogeneous region at higher chemical potential. Moreover, the pion lines decrease rapidly with growing m π and eventually disappear from the phase diagram. On the other hand, albeit reduced, the instability in the sigma channel is still present for a physical pion mass in the renormalized limit, so that we still expect an inhomogeneous phase driven by the scalar condensate.
In Fig. 4 we show the extension of the whole instability region, i.e., the whole chemical-potential interval where D −1 σ (q) is positive for some q = (0, q = 0), at vanishing temperature and varying pion mass. We find that even in the renormalized limit a finite window of instability persists for all values of m π considered. While going from the chiral limit to a physical pion mass reduces the size of the instability region, when m π becomes very large its extension starts increasing again, a similar behavior to the one observed in the NJL model away from the chiral limit [13]. 5 It is worth emphasizing that the outcome of our stability analysis is not in discrepancy with the renormalized-limit results of Ref. [14], where it was found for a CDW ansatz that the inhomogeneous phase becomes disfavored against homogeneous solutions already at m π = 37 MeV: This is due to the fact that the CDW ansatz enforces equal amplitudes for the scalar and pseudoscalar channels, the latter being disfavored according to our stability analysis. A different ansatz which allows for inhomogeneous condensation only in the σ channel on the other hand should be thermodynamically favored over homogeneous matter in this region of the phase diagram. In the following section we will demonstrate this with the specific example of the RKC modulation.

Full phase diagram for the RKC ansatz
Having determined the behavior of the instability lines for a generic inhomogeneous order parameter away from the chiral limit, we now compute the full phase diagram for a specific ansatz, the RKC modulation introduced in Sec. 2.2. To be consistent with the stability analysis, we regularize the vacuum contribution of Eq. (29) using three Pauli-Villars counterterms [9].
In Fig. 5 we show the phase diagram for a physical pion mass m π = 140 MeV, both for Λ = 757 MeV and in the renormalized limit. As expected, we find an inhomogeneous phase whose right boundary coincides with the instability line for the sigma channel found in the previous section. For comparison, we show in the figure also the left edge of the instability region, which, as expected, falls inside the inhomogeneous phase. In fact, the left edge of the instability region coincides with the first-order phase boundary one finds when the model is restricted to homogeneous phases. In other words: While the chirally almost restored phase just to the right of the first-order boundary is unstable against small inhomogeneous fluctuations, the larger homogeneous condensates to the left make it at least metastable. Our results with the RKC ansatz show however, that it is still possible to lower the free energy of the system by large inhomogeneous fluctuations in this region. Furthermore, our numerical results suggest that the tip of the inhomogeneous phase, the so-called pseudo-Lifshitz point (PLP), 6 coincides with the location of the CEP obtained when restricting the analysis to homogeneous matter. This is similar to what happens in the NJL model, and it can be understood in a general way via a Ginzburg-Landau analysis, as discussed in the following section.

Ginzburg-Landau expansion
The Ginzburg-Landau (GL) expansion is a systematic expansion of the thermodynamic potential in powers of the order parameter and its gradients. It is a powerful tool which allows to determine precisely the locations of the CEP and the PLP where both the amplitude and the gradients of the spatially modulated order parameter approach zero. In the following we want to use this method, which has been applied to the NJL model in Ref. [13], to study the behavior of the CEP and PLP in the QM model away from the chiral limit.
For this, following the steps performed in Ref. [13] for the NJL model, neglecting pseudoscalar fluctuations we write again gσ(x) =M + δM (x) and get to where the GL coefficients α i depend on T , µ andM . As shown in Ref. [13], we can localize the CEP as the point where the GL coefficients α 1 = α 2 = α 3 = 0 (the condition α 1 = 0 simply enforces the gap equation for the backgroundM ), whereas the PLP is identified as the point where both the quadratic and the first non-vanishing gradient term become zero: α 1 = α 2 = α 4,b = 0. For the relevant coefficients we find Upon close inspection, we see that when m σ,0 = 2M 0 , and thus λ = 2g 2 , the coefficients α 3 and α 4,b are proportional to each other, like in the NJL model, and as a result, the CEP and the PLP coincide, supporting the numerical results of our previous section. In Fig. 3 we have indicated the positions of these points for various values of m π by black dots.

Conclusions
We have investigated inhomogeneous phases in the renormalized limit of the quarkmeson model away from the chiral limit. Both the effect of the vacuum quark fluctuations in the QM model [10] as well as the inclusion of an explicit chiral-symmetry breaking term [13] are known to shrink the size of the inhomogeneous window in the phase diagram, so it is natural to ask whether an inhomogeneous phase survives at all when both effects are taken into account. A first investigation in this direction found that if one restricts the analysis to a CDW modulation the inhomogeneous phase quickly disappears and is not present for physical pion masses [14]. On the other hand, it is known that other types of spatial modulations of the order parameter are usually thermodynamically more favored.
Thus, in order to obtain a modulation-agnostic answer, we first looked for the appearance of instabilities of homogeneous matter towards inhomogeneous phases for arbitrary shapes of the order parameter, and found that such an instability exists even for pion masses above the physical one. This instability occurs with respect to scalar fluctuations, whereas the instability in the pseudoscalar channel disappears quickly as m π increases. Explicit chiral-symmetry breaking thus strongly suppresses fluctuations in the pseudoscalar channel, explaining the rapid disappearance of the CDW modulation from the phase diagram.
We supported these findings with an explicit calculation of the model phase diagram considering a specific modulation of the order parameter involving only the scalar channel, the so-called real-kink crystal, which provides a self-consistent ansatz away from the chiral limit, and checked that indeed the inhomogeneous phase has a non-vanishing extension for a physical m π in the renormalized limit of the model.
The presence of inhomogeneous phases thus seems to be a robust model feature, even though the size of the inhomogeneous window found in this work is relatively small. In particular, our GL analysis revealed that for m σ,0 = 2M 0 the PLP, i.e., the tip of the inhomogeneous phase, coincides with the CEP of the first-order phase boundary in the homogeneous case and, hence, the inhomogeneous phase is as robust a feature of the model as the existence of a first-order phase transition if the analysis is restricted to homogeneous phases. We must keep in mind, however, that the present analysis has been performed in mean-field approximation. It is thus an interesting question, both in the chiral limit and away from it, whether these findings remain valid if fluctuation effects are taken into account. Investigations of such questions are presently subject of intese research [25,26], particularly within the functional renormalization-group approach [7,27,28] or by performing lattice simulations for lower-dimensional models [29,30].