Crystalline phases by an improved gradient expansion technique

We develop an innovative technique for studying inhomogeneous phases with a spontaneous broken symmetry. The method relies on the knowledge of the exact form of the free energy in the homogeneous phase and on a specific gradient expansion of the order parameter. We apply this method to quark matter at vanishing temperature and large chemical potential, which is expected to be relevant for astrophysical considerations. The method is remarkably reliable and fast as compared to performing the full numerical diagonalization of the quark Hamiltonian in momentum space, and is designed to improve the standard Ginzburg-Landau expansion close to the phase transition points. For definiteness we focus on inhomogeneous chiral symmetry breaking, accurately reproducing known results for 1D and 2D modulations and examining novel crystalline structures, as well. Consistently with previous results, we find that the energetically favored modulation is the so-called 1D real kink crystal. We propose a qualitative description of the pairing mechanism to motivate this result.


I. INTRODUCTION
Quantum Chromodynamics (QCD) at nonzero baryonic densities is expected to exhibit a rich variety of phases [1]. At sufficiently large values of the quark chemical potential, µ, the chirally broken and confined phase should turn into a chirally restored deconfined phase. This phase transition is accompanied by the melting of the chiral condensate and the possible formation of a color superconducting condensate [2]. These interesting phenomena are expected to occur in a regime where perturbative QCD computations are insufficient and ab-initio lattice simulations are currently unavailable due to the sign problem, see for example [3]. Thus, effective models such as the Nambu-Jona Lasinio (NJL) model are commonly used to describe this region of the phase diagram (see [4,5] for reviews).
Remarkably, model calculations indicate that various inhomogeneous phases may arise in quark matter at high density.
Two notable examples are the crystalline color superconducting phase and the inhomogeneous chiral symmetry broken (χSB) phase (see [6,7] for reviews). The former is probably located at the onset of the deconfined phase, for neutral and beta equilibrated matter, while the latter is expected to arise between the homogeneous χSB phase and the chirally restored phase if the quark-antiquark coupling strength is sufficiently large.
More specifically, the inhomogeneous χSB island extends from zero temperature to the chiral critical point, which then turns into a Lifshitz point where three phases (homogeneous χSB, inhomogeneous χSB and the chirally restored phase) coexist [8]. The onset of this region, separating it from the traditional homogeneous χSB phase on the left, can be characterized by either a first-or second-order phase transition, depending on the crystalline shape assumed by the chiral condensate [7]. On the other hand, the chiral restoration transition is always found (at least in the chiral limit) to be a secondorder one where the chiral condensate gradually melts to zero. A remarkable consequence of this is that the position of the chiral restoration transition is independent from the type of crystalline structure considered. Regarding the color superconducting phases, it is known that at asymptotic densities the (spatially homogeneous) color-flavor locked phase [9] is favored. At densities relevant for compact stars this homogeneous phase could nevertheless be superseded by a crystalline color superconducting phase, especially when the constraints of charge neutrality and beta equilibrium are considered. However, whether this phase transition is of the first or second order is not yet established. The transition to the normal phase should be of the second order, although some modulations seem to indicate a first order phase transition. The form of the crystalline pattern has only been semiquantitatively established in [10] by a Ginzburg-Landau (GL) expansion.
Having determined the existence of an inhomogeneous island in the phase diagram, it is natural to ask which crystalline structure will be the most favored one in this region. For definiteness we focus in this work on the inhomogeneous χSB phase (but we will comment on applications to the crystalline color superconducting phase, as well). Contrary to the case of crystalline color-superconductivity, mean-field model calculations on the inhomogeneous χSB phase seem to indicate that the favored type of modulation is a one-dimensional real structure called a "real-kink crystal", which can be expressed in terms of Jacobi elliptic functions [7,11,12]. In particular, two-dimensional structures are found to be strongly disfavored compared to the one-dimensional ones, both in the vicinity of the Lifshitz point [13] as well as at zero temperature [14]. The latter result has been obtained through a computationally expensive analysis. Indeed, even within the NJL model in the mean-field approximation, the evaluation of the free energy for inhomogeneous phases is a nontrivial task. One way to obtain it is by performing a full diagonalization of the quark Hamiltonian and sum over its eigenvalues; alternatively GL expansions including gradient terms of the order parameter can be used. Both methods, however, have their own limitations: the diagonalization of the Hamiltonian can be performed analytically only in very special cases, while for most types of modulations one has to resort to a brute-force numerical computation in momentum space [14,15]. On the other hand, the GL expansion of the free energy in the order parameter and its gradients is expected to be valid only close to the second order transition to the chirally restored phase [10,13,16]. But, gradients of the order parameter may not vanish sufficiently fast even close to the second order phase transition point, making the GL expansion power counting nontrivial. One notable exception is the Lifshitz point, where both the order parameter and its gradient are expected to vanish. At any rate, while the determination of the GL coefficients is in principle straightforward, in presence of an inhomogeneous order parameter the actual derivation of the GL energy functional becomes extremely tedious at higher orders. The number of possible terms steadily increases and no automated procedure for the derivation of the coefficients has yet been developed. So far, for inhomogeneous chiral condensates expansions up to the eight order have been derived [13].
In this work, with the aim of going beyond these limitations, we build a controlled framework for investigating inhomogeneous phases away from the Lifshitz point (for the astrophysical relevant case of matter at vanishing temperature) without having to resort to a brute-force numerical diagonalization of the quark Hamiltonian in momentum space. For this, we devise an improved Ginzburg-Landau (IGL) approximation which can correctly describe both phase transitions delimiting the inhomogeneous phase to (at least in principle) arbitrary precision. This is done on one hand by implementing "by construction" a correct description of the homogeneous phase which also contains information on long wavelength modulations of the chiral condensate, and on the other hand by incorporating a sufficiently large number of appropriate gradient terms. The latter can be determined straightforwardly without having to resort to the full computation of the GL co-efficients thanks to our analytical knowledge of the eigenvalue spectrum of a simple modulation of the condensate, namely a single plane wave. This paper is organized as follows. In Sec. II we introduce the IGL approximation to describe the inhomogeneous phases, focusing on the χSB case. In Sec. III we extract the coefficients of the IGL expansion governing the transition to the chirally restored phase starting from the single plane wave modulation. In Sec. IV we compare the results of the GL and IGL approximation with the numerical results of the diagonalization of the full quark Hamiltonian. In Sec. V we analyze 2D modulations including a novel ansatz. A qualitative discussion of the pairing mechanism is given in Sec. VI. We finally draw our conclusions in Sec. VII.

II. IMPROVED GINZBURG-LANDAU EXPANSION
To develop our formalism we focus on the phenomenon of inhomogeneous χSB breaking starting from a GL expansion for the free energy in the NJL model within the mean-field approximation [4,5,17].
The idea behind the GL expansion is that close to the Lifshitz point the thermodynamic potential can be written as an expansion in powers of the chiral order parameter M (z) = −2G[S(z) + iP (z)] and its spatial derivatives (here S(z) = ψ ψ and P (z) = ψ iγ 5 τ 3 ψ are the scalar and pseudoscalar chiral condensates, respectively, and G the scalar coupling constant in the NJL Lagrangian). More specifically, for a real modulation one can write [8,13] where α n are some coefficients depending on the microscopic model. The reasoning behind this expansion is that terms with the same α n are equally important. In other words, close to the Lifshitz point both M and ∇M are small, thus M n and ∇ m M n−m , with n > m, can be comparable.
However, this is a very special case, because approaching the second order phase transition M is expected to vanish, but (∇M )/M may be nonzero. In particular, at T = 0 and close to the phase transition to the chirally restored phase the power counting underlying Eq. (1) is expected to be incorrect [7], more specifically ∇ m M n−m terms can be larger than the M n terms. In principle, we expect this to happen for any sufficiently small temperature away from the Lifshitz point, but in the following we will focus for simplicity on the T = 0 case, which is relevant for sufficiently old compact stars. Furthermore, Eq. (1) is insufficient to provide a reasonable description of the homogeneous χSB phase. This calls for a different scheme and/or a different approach.
From a technical point of view, the GL coefficients α n in the NJL model can be easily determined, either from their general expression [8], or more pragmatically by evaluating the free energy for a homogeneous order parameter and performing an expansion in powers of M , isolating the coefficients multiplying the M n terms. The knowledge of the α n is, however, insufficient to build a GL functional for inhomogeneous condensates, as different gradient terms of the same order carry different relative prefactors compared to the M n term (see Eq. (1)). The evaluation of these terms in the inhomogeneous χSB phase is already quite tedious at order α 6 , and becomes increasingly more involved at higher orders [13].
To improve the GL scheme and to circumvent the above technical difficulties, let us take one step back and inspect the structure of Eq. (1). There we have ordered the terms in such a way that order by order the first one is M n and the last one (∇ n 2 −1 M ) 2 . As already noted, these two sets of terms are particularly relevant: indeed, they are the dominant contributions close to the phase transitions. In particular, The (∇ n 2 −1 M ) 2 terms are the dominant gradient contributions close to the chirally restored phase, because terms with higher power of M are suppressed. The M n terms on the other hand are particularly relevant close to to the transition to the homogeneous χSB phase (indeed these are the only terms present in the homogeneous phase, where gradients vanish). However, the free energy of the homogeneous χSB phase is known in an analytical form. These consideration motivate the following "improved Ginzburg-Landau" (IGL) expansion, which for a real order parameter reads where the first and the last terms in the square bracket characterize this novel expansion technique.
In particular, Ω hom (M (z) 2 ), is the free energy for an homogeneous order parameter, evaluated point by point for a moving average of the mass function defined as where, as we will see, the relevant wavelength scale λ is determined by the chemical potential µ. If M is spatially constant, this term gives by construction the free energy of the homogeneous χSB phase. On the other hand, for a general oscillation it captures the long-wavelength modulation of the condensate: from a point of view of an effective field theory, this can be seen as the dominant term for long-wavelength fluctuations while high frequency components have been integrated out. The (∇ n M ) 2 term is instead the dominant one at high frequencies, of the same order or higher than µ. Indeed, the last term in Eq. (2) includes the leading gradient contributions close to the second-order transition to the chirally restored phase. The requirement for this term to be dominant is the vanishing of the amplitude of the chiral condensate, which is what happens close to the second order transition to the normal phase. We labelled the coefficients multiplying these gradient terms asα n , as they will be equal to the α n up to a numerical prefactor. By inspecting Eq. (1) we can see immediately thatα 4 = α 4 ,α 6 = α 6 /2 andα 8 = α 8 /5, whereas for higher orders these relations have not been determined until now. The first and the last term in the expansion in Eq. (2) thus guarantee the agreement with the exact result close to the two phase transitions. The other terms, which are taken from the traditional GL expansion, are expected to be relevant in the region between the two phase transitions, when ∇M ∼ M 2 , or more precisely when the modulation wavelength, λ, satisfies |M | ∼ 1/λ. As we will see in the following, the α 6 terms will be sufficient to provide a good qualitative agreement with the full numerical results, and including the α 8 terms will allow us to obtain an excellent quantitative agreement throughout the whole inhomogeneous χSB phase.
The IGL prescription can of course be generalized to complex modulations: for our novel terms this amounts to simply replacing M 2 by |M | 2 in the moving average Eq. (3) and (∇ n M ) 2 by |∇ n M | 2 in the leading gradient terms.
One important aspect is that when derived from an NJL model, the α n coefficients do not only depend on µ, but are also sensitive to the regularization scale, Λ. However, once this scale is fixed, the coefficients themselves do not depend on the considered shape of the modulation of the condensate. At vanishing temperature and for a Pauli-Villars regularization with 3 counterterms, a regulator Λ and coefficients c 0 = 1, c 1 = −3, c 2 = 3, c 3 = −1 (see [12]) we find where N f and N c are the numbers of quark flavors and colors, respectively.

III. HIGHER ORDER GRADIENTS FROM A CDW MODULATION
The only missing ingredients in the IGL expansion of Eq. (2) are theα n coefficients for n > 8. We compute them by exploiting the analytical knowledge of the eigenvalue spectrum of the quark Hamiltonian for the so-called chiral density wave (CDW) ansatz corresponding to a static single plane wave modulation, chosen without loss of generality along the z-direction, with amplitude ∆ and wave number 2q. For this simple case the quasi-particle dispersion law is known [18,19]: where = ±1, p ⊥ = p 2 x + p 2 y and E z = p 2 z + ∆ 2 . At vanishing temperature the free energy for this modulation is given by where again the Pauli-Villars regularization scheme has been adopted, see [4,7,12]. We now expand the free energy in a Taylor-like series, starting from Ω 0 = Ω(∆ = 0). Each term can be separated in a vacuum contribution Ω V n and a medium contribution Ω µ n which (at T = 0) depends on the quark chemical potential µ.
For example, the zero-th order is which is minus the pressure of a free Fermi gas of massless particles and is q-independent, as it should. The first nontrivial term is proportional to ∆ 2 , and is given by where the first term is simply a constant due to the condensation energy, and we explicited the Pauli-Villars regularization of the vacuum part. Furthermore, we isolated the medium and vacuum q-dependent contributions to Ω 2 , which can be evaluated analytically: and Upon closer inspection, we note that both contributions carry a log(q) term which could possibly spoil any expansion. However, by adding them up these log pieces cancel out. Expanding in q/µ we obtain which is a controlled expansion as long as q < µ. This typically turns out to be the case, as we will see in the following sections. Furthermore, from this expansion we can obtain all theα n terms required for the IGL expansion. Indeed, by in-specting Eq. (2) it is clear that for the CDW all the terms in the form |∇ n 2 M | 2 turn into q n ∆ 2 terms and are therefore all contained in Ω 2 . Thus, by expanding Ω 2 in powers of q as in Eq. (13) we can extract these terms to arbitrarily high order in an extremely simple way. Comparing the lower-order coefficients with the expressions in Eq. (4) we see that they agree, and pushing our expansion to higher orders we find for exampleα 10 From the above expansion it is clear that the relevant frequencies are of the order of µ, suggesting that the scale to be employed in the moving average Eq. (3) should be of the order of 1/µ. This scale is also comparable to the radius for the single soliton introduced in [20] for the real kink crystal modulation (see later), R sol = π/( √ 12M vac ) (M vac being the vacuum constituent quark mass), since the inhomogeneous phase is typically realized for µ ∼ M vac . In the following we will choose λ = 1/µ, although any other choice in the same ballpark leads to similar results.

IV. BENCHMARKS AND APPLICATIONS OF THE IGL
We are now ready to evaluate and minimize the IGL free energy and compare it with the standard GL approximation (including up to O(α 8 ) terms, see Eq. (1)) and the full numerical result. We begin by considering two different 1D modulations of the condensate, in order to explore the reliability of the GL and IGL expansions. We work in the chiral limit using a Pauli-Villars regulator Λ = 757.048 MeV and a scalar coupling G = 6.002/Λ 2 , corresponding to a vacuum constituent quark mass M vac = 300 MeV and a pion decay constant f π = 88 MeV [12]. We remind that all of our results are obtained at zero temperature.

A. Chiral Density Wave
We start by considering the CDW modulation, see Eq. (5). In Fig. 1 we report the results obtained for the variational parameters ∆, q (top panel) and the free energy at the minimum (bottom panel). For the single plane wave, the numerical results (solid black line) are extremely reliable and can be used as a benchmark to test the GL (dashed blue line) and the IGL (dotted red line) approximations. As a first step, we truncated the IGL approximation to orderα 10  function of µ in the inhomogeneous phase, eventually vanishing at the second-order transition to the chirally restored phase.
The first remarkable result visible by inspecting Fig. 1 is that even at zero temperature the standard GL expansion provides a good quantitative agreement with the results of the full numerical computation. It fails, however, to properly reproduce the numerical results in two key regions: Close to the transition to the chirally restored phase, where it overshoots the transition point, and at the transition between the homogeneous and inhomogeneous χSB phases, failing to correctly reproduce the value of ∆ in the homogeneous χSB phase and the transition point. On the other hand, the IGL exactly does what it is designed for: it improves the description of these two regions. Most notably, it exactly reproduces the free energy in the homogeneous phase and the transition to the in-homogeneous χSB phase. Moreover, it shifts the transition to the chirally restored phase closer to the numerical result.
It is important to stress that, close to the chiral restoration transition the IGL is designed for giving a systematic controlled improvement over the standard GL by the inclusion of higher order |∇ n M | 2 terms, which, as already discussed, can be straightforwardly extracted from Eq. (13). This is shown in Fig. 2, where we can see how the second-order transition is increasingly better described as we include higher order terms. In particular, we can see that to get a good qualitative agreement we need at least the O(α 8 ) terms, otherwise the inhomogeneous χSB phase extends to arbitrarily high chemical potentials. We can interpret this result by inspecting Eq. (13), or equivalently the expressions for the GL coefficients: close to the chiral restoration transition and for reasonable values of µ/Λ, the leading O(α 4 ) coefficient is negative and provides an energy gain in the formation of an inhomogeneous phase, whereas higher order terms constitute energy costs. Indeed, while we find that in principle within our regularization scheme the coefficients α 4n for n > 1 might flip sign (see Eq. (4)) and actually favor large q solutions, in practice this would require chemical potentials too close to the regulator Λ for us to trust the model results in that regime 1 . Therefore, we find that higher order gradient terms provide (increasingly smaller) energy costs which gradually push the phase transition towards lower chemical potentials, gradually approaching the full numerical result. In practice the convergence of this sum turns out to be very rapid: by includingα 10 corrections we are off the full numerical result for the transition chemical potential by only 2 MeV, and the IGL results from order α 12 on become practically indistinguishable from the numerical result. In light of this, in the following we will consider for simplicity the truncated IGL expansion at orderα 10 , with the understanding that more refined results can be straightforwardly obtained by simply adding higher order gradient terms.

B. Real Kink Crystal
The results with the CDW ansatz suggest that the IGL approximation works extremely well. As a second check we compare the GL and IGL results with the numerical ones for the modulation which has been found to be the most favored in the inhomogeneous χSB window, namely the real kink crystal (RKC) [7,12] where sn(∆z|ν) is a Jacobi elliptic function, whose shape is characterized by ∆ and by the so-called elliptic modulus ν.
After computing the cell averages over M (z) and plugging them in the GL and IGL expression, we obtain the results shown in Fig. 3. Again, we find a good agreement of the GL result with the full computation, while the IGL provides a significant quantitative improvement, reproducing the full numerical results within a few percent error. In this case the effect of the first term in the IGL expansion in Eq. (2) is more evident. It forces the average value of the condensate to match the homogeneous value, sensibly improving the agreement with the numerical results. This effect is due to the fact that the RKC ansatz can be seen as a superposition on many different waves with different amplitudes. The longwavelength amplitudes dominate close to the phase transition with the homogeneous χSB phase. On the other hand, for a CDW, there is one single spatial frequency q, which is large. Therefore in that case the IGL does not improve much with respect to the GL approximation close to this phase transition. It might be interesting at this point to compare the spatially modulated quark number density of the system obtained within the GL and IGL approximations with the numerical results of [23]. In our case, it is simply obtained by differentiating the integrand of Eq. (2) with respect to µ. This basically amounts to an improvement over a "local Fermi-gas" approximation (amounting to simply considering the first term in Eq. (2)), which has already been found to reproduce very well the behavior of the density of the system [24]. Our results are shown in Fig. 4. There we can see that once again the IGL provides a better agreement with the full result as compared to the GL, although in this case the results do not match perfectly, especially close to the phase transition to the homogeneous broken phase. This is due to the fact that ∆ and ν give the amplitude and frequency of the density oscillations. Since both are slightly different from the ones obtained by the full numerical calculation, the resulting density has amplitude and oscillation period different from the numerical ones.
The RKC modulation is somehow similar to the Lasagna phase in nuclear matter, that is a type of Pasta phase [25] expected to be realized in the inner crust of compact stars. In this phase nuclei form sheets immersed in a liquid of nuclear matter. However, with changing densities the Lasagna phase is supposed to be superseded by different modulations, possibly giving rise to higher-dimensional structures. Following this analogy one might expect that higher dimensional modulations can become favored at different values of the quark chemical potential. Another argument in favor of higherdimensional modulations comes from quarkyonic matter studies, where it is expected that increasingly complex crystalline structures can be formed by the chiral condensate as the density increases [26].

V. TWO-DIMENSIONAL STRUCTURES
Since the IGL provides very accurate results for the order parameters and free energies of one-dimensional modulations with minimal computational effort, let us now move on and consider two-dimensional structures. Close to the Lifshitz point, a systematic GL analysis of different types of higher dimensional modulations has been performed in [13], while a complementary numerical analysis for the astrophysical relevant T = 0 case can be found in [14]. Comparing the IGL results with the numerical ones of [14] for a two-dimensional square lattice with a sinusoidal ansatz, that is, we obtain the order parameters and the free energies reported in Fig. 5. It is clear that the agreement is again extremely good 2 , and we recall that the IGL result can be computed with very limited numerical effort (basically amounting to the evaluation of Ω hom (M 2 ) , as all the other terms can be computed analytically).
Using the IGL method we are in a position to easily test different 2D modulations. First we consider a square lattice with two RKC-type modulations along the x and y directions, that is M (x, y) = ∆νsn(∆x, ν)sn(∆y, ν) .
The practical implementation of this modulation in the numerical framework of [14] would be extremely complicated, as it would in principle require an expansion of the order parameter in a large number of Fourier harmonics and a minimization of the free energy with respect to all of their amplitudes. Instead, within the IGL approximation it can be straightforwardly implemented in the same way as with the 2D cosine. The minimization of the IGL free energy with respect to ∆ and ν yields qualitatively similar results to the onedimensional RKC for the order parameters. When computing the free energy associated with this modulation we find, similarly to what happens with the one-dimensional modulations, that the RKC-type solution is almost degenerate with the cosine one with the exception of the region close to the onset of the inhomogeneous phase, as shown in Fig. 6. In that figure we also see that this type of modulation is also disfavored with respect to its one-dimensional counterpart. We performed a further check in this direction by considering the ansatz M (x, y) = ∆ √ ν x sn(∆x, ν x ) + √ ν y sn(∆y, ν y ) , (18) which can interpolate between a one-dimensional RKC mod- ulation and a more involved two-dimensional structure. Consistent with our other results, we find that the minimum solution always corresponds to one of the two ν being zero, while the other reduces to the value obtained when minimizing with the one-dimensional ansatz Eq. (15).
Thus, as it was already found in [14], we can confirm within our novel approach that 2D modulations are disfavored with respect to 1D modulations at vanishing temperatures. We therefore expect that the same "hierarchy" found in [13] close to the Lifshitz point holds also at vanishing temperatures, and that 3D modulations will thus be even further disfavored compared to two-dimensional ones.

VI. QUALITATIVE ANALYSIS OF PAIRING
The comparison between the considered 2D modulations and the 1D modulations suggests that the 1D RKC is always favored. This result is in contrast with what is expected to occur in crystalline color superconductors, where a crystalline 3D pattern seems to be favored [16,27]. It is believed that in color superconductors the occurrence of the crystalline phase is due to the maximization of pairing at the Fermi surface, indeed the presence of a collective Fermi surface phenomenon seems to be the key-point for obtaining a crystalline phase.
Quite generally, a certain modulation is energetically favored if the energy gain due to pairing is larger than the energy cost of having pairs with nonvanishing total momentum. Let us examine in detail what happens in an inhomogeneous χSB condensate. For a qualitative understanding of the phenomenon we consider first the effect of a nonvanishing momentum and then we allow for pairing. For understanding whether multidimensional pairing is favored, we consider what happens for a plane wave ansatz. As discussed in [6,16], one way of representing the Fermi surface effects is to inspect the integrand of the free energy, corresponding, in our case, to the integrand appearing in Eq. (7). In Fig. 7 we plot this function at µ = 335 MeV, that is within the inhomogeneous χSB window. The left panel corresponds to the free case, that is, q = 0 and ∆ = 0. The integrand is peaked at p = µ, meaning that the larger contribution comes from the Fermi surface, corresponding to the lighter region in Fig. 7. This is the so-called pairing region, while the parts well inside the Fermi sphere or well outside it correspond to the blocking regions (see the discussion in [16] about pairing and blocking regions in color superconductors). In other words, pairing well inside/outside the Fermi sphere has a large free-energy cost, because particles should climb to the tip of the free energy (integrand), which is at the Fermi sphere. On the other hand, particle and hole excitations at the Fermi sphere are already at the tip of the mountain, that is to the largest possible energy, and they can eventually pair at no cost to form a chiral condensate [28]. Now we consider a momentum shift of the fermions. When pairs have nonvanishing total momentum, one can imagine first to displace fermions by q and then to turn on pairing. This is exactly what one does when diagonalizes the full Hamiltonian for the single plane wave to obtain the free energy in Eq. (7). This procedure is discussed in detail in [16] for crystalline color superconductors, where it is shown how by a proper momentum shift the quark propagator becomes diagonal (see also [7] for an analogous discussion for inhomogeneous chiral condensates). This momentum shift has the effect of separating the Fermi spheres, as shown in the central panel of Fig. 7 for q = 241 MeV. Now, the only pairing region corresponds to the ribbon where the two Fermi spheres touch. This picture also explains why q < µ, indeed if this were not the case the two Fermi spheres would have no overlapping regions. Exciting quasiparticles and/or holes in the pairing ribbon has no free energy cost, whereas particles from all other regions should climb an energy barrier. Indeed, we see that this is exactly what happens for ∆ = 44 MeV, right panel, where the smearing of the touching regions is exactly due to pairing. It is not possible now to add pairing in different regions of the Fermi spheres, say in the region p ⊥ ∼ 0, because these regions are too far apart and therefore the free energy cost for exciting particles and/or holes would be too high. Therefore, multidimensional modulations are dis- favored in the χSB phase because q is too large.
This does not happen in color superconductors. Indeed, one important difference between the inhomogeneous χSB phase and the crystalline color superconductors, regards the magnitude of q. Broadly speaking, q has to be proportional to the stress exerted on the pairing mechanism. In the χSB phase the stress is proportional to µ, because pairing is related to the formation of a chiral condensate. On the other hand, in color superconductors q ∝ δµ, where δµ µ is the mismatch between the Fermi spheres due to an unbalance between quarks of different flavors. For illustrative purposes, let us consider the non-energetically-favored χSB configuration corresponding to small q and ∆, somehow mimicking what happens in color superconductors. We show in Fig. 8 the integrand of Eq. (7), with q = 30 MeV and ∆ = 5 MeV. Again, pairing can happen in the ribbon where the two Fermi spheres overlap. However, it is clear now that it would be possible to slightly modify the Fermi sphere for allowing pairing, say in the p ⊥ ∼ 0 plane, at a small free-energy cost.

VII. CONCLUSIONS
We have presented a novel approach to study spatially inhomogeneous pairing by an improved Ginzburg-Landau expansion, Eq. (2). This approach relies on a scale separation between long-wavelength fluctuations, dominating the transition to the homogeneous phase, and rapid fluctuations governing the transition to the chirally restored phase. The IGL reproduces correctly by construction the homogeneous limit and allows for a description of the chiral restoration transition from the inhomogeneous phase with arbitrarily high precision by a controlled gradient expansion.
We have applied the IGL to the study of the inhomogeneous χSB phase at T = 0, reproducing the results obtained by numerical methods and extending the analysis to novel structures. These structures can hardly be studied by the numerical method, because of the complicated Fourier expansion technique underlying these methods. On the other hand, the IGL expansion turns out to be an extremely powerful tool, allowing us to quickly examine various crystalline structures to an arbitrarily accurate approximation. In this way we checked that various 2D modulations are disfavored with respect to the 1D RKC one in Eq. (15), confirming and extending previous results obtained via brute-force numerical computations [14].
It is worth emphasizing that no approximate method so far has been used to analyze the T = 0 case, probably because the standard GL approximation was assumed to be unreliable. Actually, we find that the GL approximation at the O(α 8 ) gives a surprisingly good qualitative agreement with the results of the full numerical computations. However, to make a quantitative comparison of the free energies of different struc-tures a refined approach must be used, and the IGL devised here performs this task excellently. In particular, we showed that it is able to give an accurate description of both the second order phase transition to the chirally restored phase and of the phase transition to the homogeneous χSB phase. As it turns out, a small number of additional specific gradient terms is enough to provide an excellent agreement with the numerical data.
Finally, it is worth recalling that fluctuations are expected to have a strong effect on the formation of inhomogeneous condensates [29][30][31], particularly in the case of lowerdimensional modulations [32] (see also [7] for a discussion). The inclusion of fluctuations in the IGL framework would lead to a systematic improvement beyond the mean-field approximation.
The present work can be extended in many different ways. The IGL free-energy can be used to rapidly evaluate the free energy of various crystalline color superconducting configurations, as the ones considered in [10], and to extend the analysis to novel modulations. The only modifications needed in Eq. (2) are the replacement of the free energy of the homogeneous phase with the 2SC one (for two flavors) or of the color flavor locked one (for three flavors) and to replace the α n coefficients with the pertinent ones, which can in principle be obtained by considering modulations for which the eigenvalue spectrum is known, such as a simple Fulde-Ferrell type plane wave [6,33,34]. In this case one can also compare the IGL results with those obtained by the numerical method in [15]. We will shortly present results on this topic.
Moreover, the IGL can be modified to simultaneously include the chiral and diquark condensates for examining the coexistence of the inhomogeneous χSB and of the crystalline color superconducting phase. In this case, the colorsuperconducting phase is expected to arise where the chiral condensate is small, or equivalently, where the density is large. Since 1D chiral modulations are favored, we expect that a cosine modulation, see for example [6], could be favored.