Convergence of Ginzburg-Landau expansions: superconductivity in the BCS theory and chiral symmetry breaking in the NJL model

We study the convergence of the Ginzburg-Landau (GL) expansion in the context of the Bardeen-Cooper-Schrieffer (BCS) theory for superconductivity and the Nambu–Jona-Lasinio (NJL) model for chiral symmetry breaking at finite temperature T and chemical potential µ . We present derivations of the all-order formulas for the coefficients of the GL expansions in both systems under the mean-field approximation. We show that the convergence radii for the BCS gap ∆ and dynamical quark mass M are given by ∆ conv = πT and M conv = (cid:112) µ 2 + ( πT ) 2 , respectively. We also discuss the implications of these results and the quantitative reliability of the GL expansion near the first-order chiral phase transition.


Introduction
Power series expansions are ubiquitous in physics.Some examples can be found in the perturbative expansions appearing throughout quantum mechanics and quantum field theory (QFT).Another example is the paradigm of effective field theories (EFT), based on a systematic expansion at certain scales.Hydrodynamics, for instance, is an EFT based on a gradient expansion.Ginzburg-Landau (GL) theory, originally introduced in 1950 as a phenomenological model of superconductivity, is also an EFT based on the expansion of the free energy of a system in powers of an order parameter near a phase transition.
It is well known that these power series do not always converge, and that even divergent series can generate useful results.Indeed, although the perturbative expansion in quantum electrodynamics (QED) should be a divergent asymptotic series with zero radius of convergence, as shown by Dyson [1], certain quantities calculated using perturbation theory, such as the anomalous magnetic dipole moment of the electron, are among the most precise predictions in physics.Several arguments have been put forward for the divergence of perturbative expansions in other types of QFTs, such as for scalar λϕ 3 theory [2], and by 't Hooft for quantum chromodynamics (QCD) [3].While convergent and divergent series can both be useful, there remain reasons-both practical and theoretical-for studying their convergence properties.On the practical side, it is important to know whether adding more terms to an expansion will produce better results, and if so, over what region of parameter space.On the theoretical side, at least in some situations, one can gain physical insight from an expansion's convergence or lack thereof.For recent developments on asymptotic series in QFTs, see, e.g. the review in Ref. [4].On the other hand, convergence properties of the systematic expansions for EFTs have yet to be understood generally.In this context, it was recently shown, using holographic duality methods, that the derivative expansions for specific physical quantities (shear and sound frequencies) in hydrodynamics have a finite radius of convergence in an N = 4 supersymmetric Yang-Mills plasma [5]; see also Ref. [6] for its extension to rotating plasmas.
In this paper, we study the convergence of the GL expansion.As paradigmatic examples, we consider the Bardeen-Cooper-Schrieffer (BCS) theory for superconductivity and the Nambu-Jona-Lasinio (NJL) model for chiral symmetry breaking [7] at finite temperature T and chemical potential µ.Despite the physical differences between the BCS superconductivity, characterized by the BCS gap ∆, and the chiral symmetry breaking in the NJL model, characterized by the dynamical quark mass M , their free energies in the mean-field approximation are closely related to a common mathematical expression, namely, where either ℓ = 0 or ℓ = 2.This expression must be regarded only schematically, because the first term in Eq. ( 1) is divergent and requires regularization.We will explain how J ℓ appears in each context, and then we will show that the GL expansion essentially reduces to expanding Eq. ( 1) in powers of x, subject to some form of regularization.We will then derive the nth-order coefficients of the GL expansions for arbitrary n in both systems. 1 We will show that the radius of convergence in each case is given by ∆ conv = πT and M conv = µ 2 + (πT ) 2 , respectively, and we will clarify the physical origin of the difference between these two formulas.The finite convergence radii show that results calculated using an nthorder GL expansion eventually improve with increasing n for sufficiently small values of the order parameter.
This paper is organized as follows.In Sects. 2 and 3, we derive the all-order formulas for the coefficients of the GL expansions and convergence radii in the BCS theory and NJL model, respectively.We discuss our results and give concluding remarks in Sect. 4. The technical details of the analysis are provided in appendices.
Throughout the paper, we use natural units ℏ = c = k B = 1.

GL expansion in BCS theory
Let us first review the basics of BCS theory and then derive the formula that gives the GL coefficients to all orders.Although most of the results in this section are known in the literature (see, e.g.Ref. [9]), we review these for completeness and to discuss the similarities and differences with the results of the NJL model later.

Microscopic theory and free energy
As we will not be interested in the electromagnetic responses in this paper, we can turn off the gauge fields.The BCS Lagrangian is then given by Here ψ = (ψ ↑ , ψ ↓ ) ⊤ is a two-component fermion spinor, G is a four-fermion coupling constant used to model the attractive interaction between fermions, and µ is the chemical potential, which is equal to the Fermi energy ϵ F = k 2 F /(2m).Introducing the gap parameter ∆ = G⟨ψ ⊤ Cψ⟩/2, which is assumed to be homogeneous for simplicity, with C the charge conjugation matrix, one can show that the free energy (except for the terms that do not depend on ∆) in the mean-field approximation is where β is the inverse temperature β = 1/T and ξ Here and below, we take ∆ as real without loss of generality.The integral should be taken near the Fermi surface over the region where k D is the Debye wavelength, which functions here as an ultraviolet (UV) cutoff.Near the Fermi surface, we can approximate where ρ = k 2 F /(2π 2 v F ) is the density of states per spin at the Fermi surface and ω D = v F k D is the Debye frequency.
It is now clear that after an integral transformation t = βξ, the integral in Eq. ( 4) is essentially J 0 (β∆, 0), except with the infinite upper limit replaced by a finite cutoff.Defining we have Notice that µ enters into the expression through the density of states ρ in Eq. ( 4)-hence, µ does not enter into the position of y in J 0 (x, y), unlike the case of the dynamical quark mass in the NJL model that we study in the next section. 3/21

nth-order coefficient formulas
The GL expansion is a systematic expansion of the free energy F in powers of the order parameter ∆, in this case given by Historically, GL theory was proposed as a phenomenological model of superconductivity, and the coefficients α 2 and α 4 , often denoted α and β (not to be confused with the inverse temperature) in the literature, were certain parameters.In the context of BCS theory, however, the coefficients-not only at the second and fourth orders [8], but at all orders-are determined by the microscopic theory defined in Eq. ( 2).From Eq. ( 4) we see that the GL coefficients depend on the expansion of J cut 0 (x, 0; λ) in powers of x.If we can find a general nth-order formula for these coefficients, then all the GL coefficients will immediately be known.Thus, our task is to calculate the coefficients appearing in via Note that in every expansion considered here, we ignore the constant term c 0 , because the physics is unchanged by an overall shift in the free energy.The trick to finding the (approximate) coefficients in Eq. ( 8) is to take λ → ∞ on c 2n that remain finite in this limit, which turn out to be those with 2n ≥ 4. In the weak-coupling limit, defined by ρG ≪ 1, ω D is large compared to other quantities characterizing the system, so taking λ → ∞ is physically justified.Indeed, from the condition ∂F/∂∆ = 0 for Eq.(4) in the limit T → 0, one easily finds that the zero-temperature gap ∆ 0 is given by Therefore, since in general ∆ ≤ ∆ 0 , we have λ = βω D ≫ β∆ = x, so that λ is large relative to the other variables in Eq. (5).We also comment later on how this approximation affects the radius of convergence.In this context, we could interpret J ℓ (x, y) not as divergent and requiring regularization, but as an abbreviated form of a more complicated expression in which the infinite upper limit applies only to certain terms.Taking a single derivative with respect to x 2 in Eq. ( 5), letting ℓ = 0 and y = 0, and converting the integrand to a Matsubara sum gives where ωk ≡ βω k with ω k = (2k + 1)πT being the fermionic Matsubara frequencies.Using the series expansion in x 2 of the summand, taking (n − 1) more derivatives with respect to 4/21 x2 under the integral and sum, and evaluating at x = 0, we find This expansion of the summand in x 2 has a positive and finite radius of convergence for each choice of ω2 k + t 2 , the smallest being ω2 1 = π 2 . 2 Since we will eventually substitute x = β∆, the convergence radius x 2 = π 2 corresponds to ∆ = πT .This foreshadows-but does not immediately imply-the final result, that the radius of convergence of the GL expansion is πT .
The expression in Eq. ( 12) is finite in the limit λ → ∞ for n ≥ 2. For these n, we can interchange the sum and integral, and then use the formula which follows from a recursive relation that can be obtained using integration by parts.The remaining Matsubara sum can then be related to the Riemann zeta function, yielding It follows from Eq. ( 6) that the GL coefficients of order 2n ≥ 4 are given by In particular, this reproduces the well-known result for the GL coefficient of the ∆ 4 term [8]: The result of the GL coefficient c 2 is also known [8], but we also provide its derivation to make the paper self-contained.For the coefficient c 2 , we must proceed more carefully, because we cannot simply take λ → ∞ (physically, ω D functions as a necessary UV cutoff).Using the expression in the first line of Eq. ( 11) at x = 0, and then performing the integral gives Now we clearly see a logarithmic UV divergence in the first term, although we can still take λ → ∞ in the argument of tanh and on the remaining integral, the latter of which leads to where γ is the Euler-Mascheroni constant.Thus we have and the related result for BCS theory using Eq. ( 6) is The above result can be re-expressed in terms of the gap ∆ 0 at T = 0, using Eq. ( 10), giving where T c = e γ ∆ 0 /π is the critical temperature of the superconducting state, at which α 2 changes sign.Combining Eq. ( 10) with the formula for T c , we also find

Radius of convergence
The radius of convergence of the expansion in Eq. ( 8) can essentially be read off of the coefficient formula (14).We have (x 2 ) conv = π 2 when viewed as a series in x 2 , or equivalently, x conv = π when viewed as a series in x.To make the argument rigorous, one can apply the ratio test, (x 2 ) conv = lim n→∞ |c 2n+2 /c 2n |.Finally, Eq. ( 6) shows that the GL expansion is essentially given by the same series, with x = β∆, so we have ∆ conv = x conv T = πT .Let us explain the exact meaning and implications of this result.The radius of convergence π applies technically to the series whose coefficients c 2n are given by Eq. ( 14), but these are computed in the limit λ → ∞, in which the original quantity J cut 0 becomes ill-defined.Nonetheless, it is easy to see that the true coefficients in the expansion of J cut 0 , keeping λ finite, are bounded from above in magnitude by c 2n in Eq. ( 14).This follows because the integrand of Eq. ( 12) is either entirely positive or entirely negative, so increasing the upper limit of integration must strictly increase its absolute value.Since c 2n of Eq. ( 14) lead to a series with radius of convergence π, and the series for J cut 0 with λ finite has its coefficients bounded by the former, we can actually conclude that the radius of convergence is at least π for the true expansion of J cut 0 , and hence the GL expansion has radius of convergence of at least πT .
When using the GL expansion to approximate and solve the gap equation, the key question is whether the true gap ∆ exact falls within this convergence radius, ∆ exact < ∆ conv .When this condition is satisfied, then the solution of the N th-order GL expansion, ∆ GL,N , will converge to ∆ exact as N → ∞.This is the convergence of the GL solution.Note that this should be distinguished from the convergence of the GL expansion itself over some (possibly vanishing) range 0 ≤ ∆ < ∆ conv for each fixed T .
The range of temperatures over which the GL solution converges is indicated in the left panel of Fig. 1 (orange shaded region).We searched numerically for the temperature at which ∆ exact drops below ∆ conv over a range of parameter values for ρG.For each ρG, the value of ω D /T c is fixed by Eq. ( 22), so ρG is the only parameter of the theory when all quantities are expressed relative to T c .Note that the lower boundary of the convergence region is nearly constant over the given range of ρG (e.g. for ρG ≤ 0.21 the GL solution converges when T /T c ≥ 0.53).The right panel of Fig. 1 shows ∆ vs. T computed from the gap equation and also from 6th-and 20th-order GL expansions.Both GL expansions are reliable near T c , as expected, and the 20th-order GL solutions remain reliable over a wider range of T .Let us also comment on how the GL solutions in the right panel of Fig. 1 are calculated.Since solving the stationary equation ∂F/∂∆ = 0 for an N th-order GL expansion involves finding the roots of an (N − 1)th-order polynomial in ∆, one should expect in general to find up to N − 1 potential solutions.Moreover, if the largest coefficient in such an N th-order GL expansion is negative, then the free energy must eventually decrease toward −∞ as ∆ → ∞, so one cannot determine the GL solutions simply by finding the global minimum.Instead, we essentially just search for the local minimum with the smallest free energy over 0 ≤ ∆ ≤ ∆ conv .Although many more real roots of the stationary equation can appear with increasing N , in practice these occur far outside the radius of convergence, and hence they are easily neglected.The only additional complication is that the minimum in the 20th-order case gradually shifts and moves outside the convergence region as T /T c drops below 0.53, and we continue to plot this minimum in the right panel of Fig. 1 (blue dotted curve left of vertical line).We treat this as the solution because it remains close to the convergence region and is continuously connected to the solution at higher T .Note that this minimum remains well defined at all T because α 20 is positive; by contrast, the red dashed curve vanishes for T /T c < 0.69 because the relevant local minimum at 6th order truly vanishes at these temperatures.

GL expansion in the NJL Model
In this section, we show how the expression J 2 (x, y) appears in the GL expansion for the dynamical quark mass in the NJL model and how it can be regularized in this setting.We then solve for the nth-order GL coefficients and convergence radius.

Microscopic theory and free energy
The two-flavor NJL Lagrangian is given by [7,[12][13][14]: Here ψ = (u, d) ⊤ is a two-flavor quark field, m = diag(m u , m d ) is the bare quark mass matrix in flavor space, G is a four-fermion coupling constant, and τ a are the Pauli matrices acting in flavor space.This model retains the approximate chiral symmetry of QCD, which becomes exact in the chiral limit m → 0. We will assume the chiral limit from now on.In vacuum, chiral symmetry is spontaneously broken by the formation of a quark-antiquark pairing σ = ⟨ ψψ⟩, called the chiral condensate.To study the behavior of the condensate σ as a function of temperature T and quark chemical potential µ, we expand Eq. ( 23) about the ansatz: We note that at sufficiently large densities, or even small finite densities and in the presence of a magnetic field, this ansatz is known to be disfavored over various spatially inhomogeneous condensates within the mean-field approximation [15], the simplest being the so-called chiral density wave, given by ⟨ ψψ⟩ + i⟨ ψiγ 5 τ 3 ψ⟩ = σe iqz .In this paper, we restrict to the homogeneous case for simplicity.
In the mean-field approximation, one can show that the free energy is given by [13,14]: where M = −2Gσ can be interpreted as the dynamical mass of quarks due to the chiral condensate; N f = 2 and N c = 3 are the numbers of flavors and colors, respectively; and Making the transformations x = βk and y = βµ, we find Here, the divergence of J 2 reflects the zero-point vacuum energy, and this divergence must be regularized in order to carry out numerical calculations.Several schemes are in common use [12].For consistency with the last section, we will impose a simple momentum cutoff Λ, and take Λ → ∞ on any terms that remain finite in this limit.Such a cutoff clearly violates Lorentz invariance, which can lead to undesired consequences when the condensate under consideration is inhomogeneous [15].In that context, covariant regularization methods, such as the Schwinger proper time method [16], are typically used instead.Since we consider only homogeneous condensates in this paper, the following analytical results based on a simple cutoff are physically meaningful.However, the numerical accuracy of the approximation Λ = ∞ is somewhat limited in this setting, as we will discuss in Sect.3.3.It is also possible to solve for the GL coefficients using the proper time method, and in fact, that allows for analytical formulas for the GL coefficients even without taking the limit Λ → ∞; these results are given in Appendix A. For the sake of accuracy, we use the proper time method in the numerical results presented in Fig. 2. Imposing the momentum cutoff |k| < Λ in Eq. ( 25), we find 8/21

nth-order coefficient formulas
As in the previous case, the GL expansion takes the form The coefficient α 4 was previously computed in Ref. [17], but the expression for the generic coefficient α 2n has not been calculated explicitly, to the best of our knowledge.These coefficients are determined by the expansion coefficients of J cut 2 (x, y; λ) in powers of x, via which we will now derive.Finally, once the c 2n coefficients are known, the α 2n coefficients follow immediately from Eq. ( 27), The two differences from the BCS case are that we now have J 2 rather than J 0 , and we now have y ̸ = 0.The appearance of J 2 adds a factor of t 2 to the integrand, which strengthens the UV divergence, and as a result, we will need to keep λ finite on the c 4 coefficient in addition to c 2 .Having y ̸ = 0 complicates the algebra, but does not significantly affect the overall approach.
Following the same procedure as in the BCS case, we start by taking a single derivative with respect to x 2 , giving where ωk = (2k + 1)π, as before.After applying the identity we obtain At this point, in close analogy to the comment after Eq. ( 12), we may see a hint of the final result for the radius of convergence.If one expands the integrand in powers of x 2 , the resulting series has radius of convergence x 2 conv = |t 2 + (ω k + iy) 2 |.Taking t → 0 then gives x 2 conv = |ω k + iy| 2 , the minimum of which is x 2 conv = π 2 + y 2 ; recalling that x = βM and y = βµ, this corresponds to M conv = µ 2 + (πT ) 2 .However, this heuristic argument is fairly crude, partly because the quantity |t 2 + (ω k + iy) 2 | need not be minimized in the limit t → 0. We therefore proceed to the rigorous proof.

9/21
Taking n − 1 more derivatives of Eq. (34) with respect to x 2 , introducing a variable s = t/(ω k + iy), and letting λ → ∞, we find The integral is easily performed by using Eq. ( 13), giving Notice that the sum in Eq. ( 35) converges when n ≥ 3. Identifying the sum with the standard series representation of the nth-order polygamma function For c 2 and c 4 , it is more convenient to avoid writing J 2 as a Matsubara sum, starting instead from the expression For c 2 , it is straightforward to calculate where Li n (x) = ∞ k=1 x k /k n is the nth polylogarithm.The first two terms under the sum in Eq. (39) vanish as λ → ∞, and for the last term we can apply the identity [18], where B n (x) is the nth Bernoulli polynomial.The result is The c 4 coefficient is more subtle because a naive separation of Eq. ( 38) into vacuum and medium contributions leads to IR divergences.An efficient way to compute this coefficient is to note that after commuting the derivatives ∂ 2 x 2 with t 2 in the integrand, we can transform them as ∂ x 2 → (2t) −1 ∂ t .This allows us to take the limit x → 0 and then take derivatives with respect to t. Performing only the rightmost derivative, we find Integrating by parts, taking λ → ∞ on the boundary terms, and by applying the formula [19], where ψ(z) is the digamma function defined by ψ(z) = d ln Γ(z)/dz with Γ(z) the gamma function, we obtain for large λ, and accordingly, Finally, applying Eq. ( 31), we have

Radius of convergence
We will show that the expansion of J cut 2 (x, y; λ) has radius of convergence x conv = π 2 + y 2 .Then, from Eq. ( 27), it will follow that the GL expansion has convergence radius First, since the convergence of any series only depends on its infinite tail, and not on any finite initial segment, we can focus on the coefficients c 2n≥6 given by Eq. (37).It is straightforward to show that the radius of convergence is at least π 2 + y 2 .Using the inequality | Re(z)| ≤ |z| for complex z, we have for 2n ≥ 6.
Using the identity ψ (n) (x) = (−1) n+1 n! ζ(n + 1, x), where ζ(s, a) = ∞ k=0 1/(k + a) s is the Hurwitz zeta function, one can show that for Re(x) > 0, in the sense that the ratio of the quantities on the left and right sides of the arrow approaches unity as n → ∞.It follows that Thus, the series for J cut 2 (x, y; λ) in Eq. ( 29) is bounded in magnitude by a series with radius of convergence x conv = π 2 + y 2 , so the original series has a convergence radius of at least this value.In fact, one can show the exact convergence radius is indeed π 2 + y 2 , but the proof is more involved.In particular, if we do not take the upper bound of the coefficients using | Re(z)| ≤ |z|, then the ratio test does not work, because the quantity [Re ψ (2n−4) ( Although the results of this section were computed using a simple momentum cutoff Λ, the radius of convergence is unaffected when using Schwinger proper time regularization.Whereas the former cutoff scheme involves taking the limit Λ → ∞ to obtain parts of c 2n and α 2n analytically, the latter scheme allows deriving their analytic expressions even without taking such a limit for a cutoff parameter Λ.As shown in Appendix A, the α 2n≥6 coefficients are almost the same as in Eq. ( 48), except with small corrections proportional to inverse powers of Λ.The key fact is that these correction terms decay faster than the α 2n coefficients calculated in this section, so they do not increase the convergence radius.We prove these statements in Appendix A.
Let us also mention that the assumption of Λ = ∞ made when deriving Eqs. ( 46)-( 48) is only a rough approximation in the NJL setting.The values for the two parameters Λ and G are determined by a choice of values for M 0 and f π , where M 0 is the dynamical quark mass and f π is the pion decay constant at T = µ = 0. Choosing M 0 = 300 MeV and f π = 88 MeV in the chiral limit [15,20], we find Λ = 614 MeV and GΛ 2 = 2.15 [12].With these values, one can numerically find T c = 182 MeV at µ = 0, and µ c = 311 MeV, where µ c is the chemical potential at which M vanishes in a first-order transition with T = 0 fixed.Thus, Λ is not especially large on the scale of other characteristic quantities of the system, so taking Λ → ∞ is not fully justified.For instance, finding T c by setting α 2 = 0 in Eq. ( 46) gives T c = 165 MeV.Although the series with coefficients given by Eqs. ( 46)-(48) converges with radius µ 2 + (πT ) 2 , the value does not converge exactly to Ω, even within this radius.
These issues are avoided when using Schwinger proper time regularization, described in Appendix A, because there the finite size of Λ is captured in the analytical formulas for the coefficients.Therefore, the numerical results for this section, presented in Fig. 2, are computed using the Schwinger method.In this regularization scheme, again choosing M 0 = 300 MeV and f π = 88 MeV, we find Λ = 634 MeV and GΛ 2 = 6.02.We calculated T c over the range 0 ≤ µ ≤ 312 MeV, and we determined the region in the µ-T plane below T c where M < M conv , i.e.where the GL solution converges (top panel, orange shaded region).The bottom panels of Fig. 2 show M vs. T computed from the gap equation and also from the GL expansion at orders 6 and 20.For µ = 0 (bottom left panel), the GL solution converges only when T > 92 MeV (right of the vertical line); both the 6th-and 20th-order solutions are very accurate near T c , but the 6th-order solution becomes noticeably inaccurate at T ≲ 130 MeV.On the other hand, when µ = 300 MeV (bottom right panel), the phase transition is first order, and the 20th-order solution is a noticeable improvement over the 6th-order solution across the entire range 0 < T < T c .Moreover, the GL solution converges over this entire range, and hence the 20th-order solutions are very reliable all the way down to T = 0.This can be understood by considering the radius of convergence formula at T = 0, which reduces to M conv = µ, and noticing that M ≤ 300 MeV for µ > 300 MeV.
We also note that the method used for determining the solutions of the 6th-and 20th-order GL expansions, plotted in the lower panels of Fig. 2, is exactly the same as in the BCS case, as discussed in the last paragraph of Sect.2.3.
As mentioned above and shown in Appendix A, the radius of convergence when using the proper time method is also µ 2 + (πT ) 2 .One could also ask about the radius of convergence in the cutoff method if not using the approximation Λ = ∞, i.e. if instead one were to calculate the GL coefficients from Eq. ( 27 μ = 300 MeV Fig. 2. Top: Region in which the GL solution for the dynamical quark mass converges in the NJL model using Schwinger proper time regularization.The black curve shows the critical temperature, with solid (dashed) lines indicating the first-order (second-order) phase transitions.The orange shaded region indicates where the GL solution converges when M > 0. Bottom: Dynamical quark mass versus temperature at µ = 0 (left) and µ = 300 MeV (right), computed by solving the gap equation using the exact free energy (black), a 6thorder GL expansion (pink, dashed), and a 20th-order GL expansion (blue, dotted).that in this case the radius of convergence is still given by at least µ 2 + (πT ) 2 when M > 0; we sketch a proof in Appendix B.2.

Discussions and outlook
In this paper, we studied the convergence of the GL expansion for two prominent examplesthe BCS theory for superconductivity and the NJL model for chiral symmetry breaking at finite temperature T and chemical potential µ.We have shown that the convergence radii are given by ∆ conv = πT and M conv = µ 2 + (πT ) 2 , respectively.The difference between these two expressions can be understood physically as follows.In the BCS theory, Cooper pairing 13/21 occurs close to the Fermi surface, so µ dependence enters only through the density of states, which is universal for all the GL coefficients.In the NJL model, on the other hand, the chiral condensate is a pairing between quarks and antiquarks inside the Dirac sea, so the GL coefficients depend rather nontrivially on µ, as shown in Eqs. ( 46)-(48).We also note that this difference leads to the known fact that the BCS instability appears for an infinitesimally small attractive interaction between fermions, whereas generation of the chiral condensate requires a sufficiently strong attractive interaction between quarks and antiquarks [7].
In this paper, we focused on the two-flavor NJL model, where only even powers of M appear in the GL expansion.In the three-flavor case, on the other hand, odd powers in M can also appear in the expansion [21] due to the instanton-induced interaction (or Kobayashi-Maskawa-'t Hooft interaction) [22][23][24].It would be interesting to see how the radius of convergence may be affected by such terms.One can also ask about the radius of convergence of the GL theory for the chiral phase transition in QCD per se.For the mean-field approximation to be well justified, one can take the large-N c limit [25,26].Because the radius of convergence for the dynamical quark mass in the NJL model above does not depend on the details of the model, such as the four-fermion coupling constant G and cutoff Λ, 3 one may conjecture that large-N c QCD has the same radius of convergence, M conv = µ 2 + (πT ) 2 . 4 It would be interesting to extend our approach to the diquark condensate for color superconductivity [28], the interplay between chiral and diquark condensates [29], the inhomogeneous condensates like Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) pairing [30,31] or inhomogeneous chiral condensates [15], and other orders in condensed matter systems.We defer these questions to future work.

Fig. 1 .
Fig.1.Left: Region in which the GL solution converges for BCS theory.The black dashed line denotes the second-order phase transition, and the orange shaded region indicates where the GL solution converges with ∆ > 0. Right: BCS gap versus temperature, computed by solving the gap equation using the exact free energy (black), a 6th-order GL expansion (pink, dashed), and a 20th-order GL expansion (blue, dotted).GL solution should converge to the right of the vertical gray line, i.e.T /T c ≥ 0.53.These data are computed for ρG = 0.18.
) numerically at finite Λ.It is possible to show 12