Convergence of Ginzburg–Landau Expansions:
Superconductivity in the Bardeen–Cooper–Schrieffer Theory and
Chiral Symmetry Breaking in the Nambu–Jona-Lasinio Model

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

$$\mathcal{L}_{\text{BCS}}=\psi^{\dagger}\left(i\partial_{t}+\frac{\bm{\nabla}^{2}}{2m}+\mu\right)\psi+\frac{G}{2}(\psi^{\dagger}\psi)^{2}\,.$$

(2)

Here $\psi=(\psi_{\uparrow},\psi_{\downarrow})^{\top}$ is a two-component fermion spinor, $G$ is a four-fermion coupling constant used to model the attractive interaction between fermions, and $\mu$ is the chemical potential, which is equal to the Fermi energy $\epsilon_{\rm F}=k_{\rm F}^{2}/(2m)$. Introducing the gap parameter $\Delta={G}\langle\psi^{\top}C\psi\rangle/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 $\Delta$) in the mean-field approximation is

$$F=\frac{\Delta^{2}}{G}-2T\int_{|\bm{k}|\approx k_{\rm F}}\frac{d^{3}k}{(2\pi)^{3}}\ln\bigg{[}2\cosh\left(\frac{\beta}{2}\sqrt{\Delta^{2}+\xi_{\bm{k}}^{2}}\right)\bigg{]}\,,$$

(3)

where $\beta$ is the inverse temperature $\beta=1/T$ and $\xi_{\bm{k}}=\epsilon_{\bm{k}}-\mu$ with $\epsilon_{\rm F}=|{\bm{k}}|^{2}/(2m)$. Here and below, we take $\Delta$ as real without loss of generality.

The integral should be taken near the Fermi surface over the region $k_{\rm F}-k_{\rm D}<|\bm{k}|<k_{\rm F}+k_{\rm D}$, where $k_{\rm D}$ is the Debye wavelength, which functions here as an ultraviolet (UV) cutoff. Near the Fermi surface, we can approximate $d^{3}k\approx 4\pi k_{\rm F}^{2}dk$ and $\xi_{\bm{k}}\approx v_{\rm F}(k-k_{\rm F})$ with $v_{\rm F}=k_{\rm F}/m$ being the Fermi velocity. Therefore, assuming $k_{\rm D}\ll k_{\rm F}$, we have

$$F=\frac{\Delta^{2}}{G}-4\rho T\int_{0}^{\omega_{\rm D}}d\xi\,\ln\bigg{[}2\cosh\left(\frac{\beta}{2}\sqrt{\Delta^{2}+\xi^{2}}\right)\bigg{]}\,,$$

(4)

where $\rho=k_{\rm F}^{2}/(2\pi^{2}v_{\rm F})$ is the density of states per spin at the Fermi surface and $\omega_{\rm D}=v_{\rm F}k_{\rm D}$ is the Debye frequency.

It is now clear that after an integral transformation $t=\beta\xi$, the integral in Eq. (4) is essentially $J_{0}(\beta\Delta,0)$, except with the infinite upper limit replaced by a finite cutoff. Defining

$\displaystyle J_{\ell}^{\text{cut}}(x,y;\lambda)$

$\displaystyle=\int_{0}^{\lambda}dt\,t^{\ell}\sum_{\zeta=\pm}\ln\bigg{[}2\cosh\bigg{(}\frac{\sqrt{x^{2}+t^{2}}+\zeta y}{2}\bigg{)}\bigg{]}\,,$

(5)

we have

$$F=\frac{\Delta^{2}}{G}-2\rho T^{2}J_{0}^{\text{cut}}(\beta\Delta,0;\beta\omega_{\rm D})\,.$$

(6)

Notice that $\mu$ enters into the expression through the density of states $\rho$ in Eq. (4)—hence, $\mu$ 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.

The GL expansion is a systematic expansion of the free energy $F$ in powers of the order parameter $\Delta$, in this case given by

$$F=\alpha_{2}\Delta^{2}+\alpha_{4}\Delta^{4}+\cdots\;.$$

(7)

Historically, GL theory was proposed as a phenomenological model of superconductivity, and the coefficients $\alpha_{2}$ and $\alpha_{4}$, often denoted $\alpha$ and $\beta$ (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 Gorkov1959 , 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_{0}^{\text{cut}}(x,0;\lambda)$ in powers of $x$. If we can find a general $n$th-order formula for these coefficients, then all the GL coefficients will immediately be known. Thus, our task is to calculate the coefficients appearing in

$$J_{0}^{\text{cut}}(x,0;\lambda)=c_{2}x^{2}+c_{4}x^{4}+\cdots\;,$$

(8)

via

$$c_{2n}=\frac{1}{n!}\partial_{x^{2}}^{n}J_{0}^{\text{cut}}(x,0;\lambda)\Big{|}_{x=0}\,.$$

(9)

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 $\lambda\to\infty$ on $c_{2n}$ that remain finite in this limit, which turn out to be those with $2n\geq 4$. In the weak-coupling limit, defined by $\rho G\ll 1$, $\omega_{\rm D}$ is large compared to other quantities characterizing the system, so taking $\lambda\to\infty$ is physically justified. Indeed, from the condition $\partial F/\partial\Delta=0$ for Eq. (4) in the limit $T\rightarrow 0$, one easily finds that the zero-temperature gap $\Delta_{0}$ is given by

$$\Delta_{0}=2\omega_{\rm D}e^{-1/(\rho G)}\,.$$

(10)

Therefore, since in general $\Delta\leq\Delta_{0}$, we have $\lambda=\beta\omega_{\rm D}\gg\beta\Delta=x$, so that $\lambda$ 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_{\ell}(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 $\ell=0$ and $y=0$, and converting the integrand to a Matsubara sum gives

	$\displaystyle\partial_{x^{2}}J_{0}^{\text{cut}}$	$\displaystyle=2\int_{0}^{\lambda}dt\,\frac{\tanh\left(\frac{1}{2}\sqrt{x^{2}+t^{2}}\,\right)}{4\sqrt{x^{2}+t^{2}}}$
		$\displaystyle=2\int_{0}^{\lambda}dt\,\sum_{k=0}^{\infty}\frac{1}{\bar{\omega}_{k}^{2}+x^{2}+t^{2}}\,,$		(11)

where $\bar{\omega}_{k}\equiv\beta\omega_{k}$ with $\omega_{k}=(2k+1)\pi T$ being the fermionic Matsubara frequencies. Using the series expansion in $x^{2}$ of the summand, taking $(n-1)$ more derivatives with respect to $x^{2}$ under the integral and sum, and evaluating at $x=0$, we find

$\displaystyle\partial_{x^{2}}^{n}J_{0}^{\text{cut}}\Big{|}_{x=0}$

$\displaystyle=2\int_{0}^{\lambda}dt\,\sum_{k=0}^{\infty}\frac{(-1)^{(n-1)}(n-1)!}{(\bar{\omega}_{k}^{2}+t^{2})^{n}}\,.$

(12)

This expansion of the summand in $x^{2}$ has a positive and finite radius of convergence for each choice of $\bar{\omega}_{k}^{2}+t^{2}$, the smallest being $\bar{\omega}_{1}^{2}=\pi^{2}$.²²2In the bosonic case, on the other hand, the smallest value of $\omega_{k}^{2}$ is zero. Integrating this term and expanding in $x$ gives a series with an odd term $|x|^{3}$, which is not analytic in $x^{2}$ (see, e.g. Ref. Laine2016 ). Such a nonanalytic term does not appear in the fermionic case treated here due to nonzero $\omega_{k}$, which acts as an infrared (IR) cutoff. Since we will eventually substitute $x=\beta\Delta$, the convergence radius $x^{2}=\pi^{2}$ corresponds to $\Delta=\pi T$. This foreshadows—but does not immediately imply—the final result, that the radius of convergence of the GL expansion is $\pi T$.

The expression in Eq. (12) is finite in the limit $\lambda\to\infty$ for $n\geq 2$. For these $n$, we can interchange the sum and integral, and then use the formula

$$\int_{0}^{\infty}\frac{dt}{(1+t^{2})^{n}}=\frac{(2n-3)!!}{(2n-2)!!}\frac{\pi}{2}\,,$$

(13)

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

$$c_{2n\geq 4}=\frac{(-1)^{n+1}}{n}\frac{(2n-3)!!}{(2n-2)!!}\frac{1}{\pi^{2n-2}}\big{(}1-2^{-(2n-1)}\big{)}\zeta(2n-1)\;.$$

(14)

It follows from Eq. (6) that the GL coefficients of order $2n\geq 4$ are given by

$$\alpha_{2n\geq 4}=-\frac{2\rho}{T^{2n-2}}c_{2n}\,.$$

(15)

In particular, this reproduces the well-known result for the GL coefficient of the $\Delta^{4}$ term Gorkov1959 :

$$\alpha_{4}=\frac{7\zeta(3)}{16}\frac{\rho}{(\pi T)^{2}}\,.$$

(16)

The result of the GL coefficient $c_{2}$ is also known Gorkov1959 , 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 $\lambda\to\infty$ (physically, $\omega_{\rm 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

$\displaystyle\partial_{x^{2}}J^{\text{cut}}_{0}\Big{|}_{x=0}$

$\displaystyle=\frac{1}{2}\left[\ln\left(\frac{\lambda}{2}\right)\tanh\left(\frac{\lambda}{2}\right)-\int_{0}^{\lambda/2}dt\,\ln(t)\operatorname{sech}^{2}(t)\right].$

(17)

Now we clearly see a logarithmic UV divergence in the first term, although we can still take $\lambda\to\infty$ in the argument of $\tanh$ and on the remaining integral, the latter of which leads to

$\displaystyle\int_{0}^{\infty}dt\,\ln(t)\operatorname{sech}^{2}(t)=-\ln\left(\frac{4e^{\gamma}}{\pi}\right)\,,$

(18)

where $\gamma$ is the Euler–Mascheroni constant. Thus we have

$$c_{2}=\frac{1}{2}\ln\left(\frac{2e^{\gamma}}{\pi}\lambda\right),$$

(19)

and the related result for BCS theory using Eq. (6) is

$$\alpha_{2}=\frac{1}{G}-\rho\ln\left(\frac{2e^{\gamma}}{\pi}\frac{\omega_{\rm D}}{T}\right)\,.$$

(20)

The above result can be re-expressed in terms of the gap $\Delta_{0}$ at $T=0$, using Eq. (10), giving

$$\alpha_{2}=\rho\ln\left(\frac{T}{T_{\rm c}}\right)\,,$$

(21)

where $T_{\rm c}=e^{\gamma}\Delta_{0}/\pi$ is the critical temperature of the superconducting state, at which $\alpha_{2}$ changes sign. Combining Eq. (10) with the formula for $T_{\rm c}$, we also find

$$\frac{\omega_{\rm D}}{T_{\rm c}}=\frac{\pi}{2}{e}^{-\gamma+1/(\rho G)}\,.$$

(22)

The radius of convergence of the expansion in Eq. (8) can essentially be read off of the coefficient formula (14). We have $(x^{2})_{\text{conv}}=\pi^{2}$ when viewed as a series in $x^{2}$, or equivalently, $x_{\text{conv}}=\pi$ when viewed as a series in $x$. To make the argument rigorous, one can apply the ratio test, $(x^{2})_{\text{conv}}=\lim_{n\to\infty}|c_{2n+2}/c_{2n}|$. Finally, Eq. (6) shows that the GL expansion is essentially given by the same series, with $x=\beta\Delta$, so we have $\Delta_{\text{conv}}=x_{\text{conv}}T=\pi T$.

Let us explain the exact meaning and implications of this result. The radius of convergence $\pi$ applies technically to the series whose coefficients $c_{2n}$ are given by Eq. (14), but these are computed in the limit $\lambda\to\infty$, in which the original quantity $J^{\text{cut}}_{0}$ becomes ill-defined. Nonetheless, it is easy to see that the true coefficients in the expansion of $J_{0}^{\text{cut}}$, keeping $\lambda$ 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 $\pi$, and the series for $J_{0}^{\text{cut}}$ with $\lambda$ finite has its coefficients bounded by the former, we can actually conclude that the radius of convergence is at least $\pi$ for the true expansion of $J_{0}^{\text{cut}}$, and hence the GL expansion has radius of convergence of at least $\pi T$.

When using the GL expansion to approximate and solve the gap equation, the key question is whether the true gap $\Delta_{\text{exact}}$ falls within this convergence radius, $\Delta_{\text{exact}}<\Delta_{\text{conv}}$. When this condition is satisfied, then the solution of the $N$th-order GL expansion, $\Delta_{\text{GL},N}$, will converge to $\Delta_{\text{exact}}$ as $N\to\infty$. 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\leq\Delta<\Delta_{\text{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 $\Delta_{\text{exact}}$ drops below $\Delta_{\text{conv}}$ over a range of parameter values for $\rho G$. For each $\rho G$, the value of $\omega_{\rm D}/T_{\rm c}$ is fixed by Eq. (22), so $\rho G$ is the only parameter of the theory when all quantities are expressed relative to $T_{\rm c}$. Note that the lower boundary of the convergence region is nearly constant over the given range of $\rho G$ (e.g. for $\rho G\leq 0.21$ the GL solution converges when $T/T_{\rm c}\geq 0.53$). The right panel of Fig. 1 shows $\Delta$ vs. $T$ computed from the gap equation and also from 6th- and 20th-order GL expansions. Both GL expansions are reliable near $T_{\rm 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 $\partial F/\partial\Delta=0$ for an $N$th-order GL expansion involves finding the roots of an $(N-1)$th-order polynomial in $\Delta$, 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 $-\infty$ as $\Delta\to\infty$, 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\leq\Delta\leq\Delta_{\text{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_{\rm 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 $\alpha_{20}$ is positive; by contrast, the red dashed curve vanishes for $T/T_{\rm c}<0.69$ because the relevant local minimum at 6th order truly vanishes at these temperatures.

The two-flavor NJL Lagrangian is given by Nambu:1961tp ; Klevansky:1992qe ; Hatsuda:1994pi ; Buballa:2003qv :

$$\mathcal{L}_{\text{NJL}}=\bar{\psi}(i\gamma^{\mu}\partial_{\mu}-\hat{m})\psi+G\big{[}(\bar{\psi}\psi)^{2}+(\bar{\psi}i\gamma^{5}\tau^{a}\psi)^{2}\big{]}.$$

(23)

Here $\psi=({\rm u,d})^{\top}$ is a two-flavor quark field, $\hat{m}=\operatorname{diag}(m_{\rm u},m_{d})$ is the bare quark mass matrix in flavor space, $G$ is a four-fermion coupling constant, and $\tau^{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 $\hat{m}\to 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 $\sigma=\langle\bar{\psi}\psi\rangle$, called the chiral condensate. To study the behavior of the condensate $\sigma$ as a function of temperature $T$ and quark chemical potential $\mu$, we expand Eq. (23) about the ansatz:

$$\langle\bar{\psi}\psi\rangle=\sigma,\qquad\langle\bar{\psi}i\gamma^{5}\tau^{a}\psi\rangle=0\quad(a=1,2,3).$$

(24)

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 Buballa:2014tba , the simplest being the so-called chiral density wave, given by $\langle\bar{\psi}\psi\rangle+i\langle\bar{\psi}i\gamma^{5}\tau^{3}\psi\rangle=\sigma 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 Hatsuda:1994pi ; Buballa:2003qv :

$$\Omega_{\text{pre-reg}}=\frac{M^{2}}{4G}-2N_{\rm f}N_{\rm c}\int\frac{d^{3}k}{(2\pi)^{3}}\left[E_{\bm{k}}+\sum_{\zeta=\pm 1}T\ln\left(1+e^{-\beta(E_{\bm{k}}+\zeta\mu)}\right)\right],$$

(25)

where $M=-2G\sigma$ can be interpreted as the dynamical mass of quarks due to the chiral condensate; $N_{\rm f}=2$ and $N_{\rm c}=3$ are the numbers of flavors and colors, respectively; and $E_{\bm{k}}=\sqrt{\bm{k}^{2}+M^{2}}$. Making the transformations $x=\beta k$ and $y=\beta\mu$, we find

$$\Omega_{\text{pre-reg}}=\frac{M^{2}}{4G}-\frac{N_{\rm f}N_{\rm c}}{\pi^{2}}T^{4}J_{2}(\beta M,\beta\mu)\,.$$

(26)

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 Klevansky:1992qe . For consistency with the last section, we will impose a simple momentum cutoff $\Lambda$, and take $\Lambda\to\infty$ 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 Buballa:2014tba . In that context, covariant regularization methods, such as the Schwinger proper time method Schwinger:1951nm , 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 $\Lambda=\infty$ 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 $\Lambda\to\infty$; 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 $|\bm{k}|<\Lambda$ in Eq. (25), we find

$$\Omega^{\text{cut}}=\frac{M^{2}}{4G}-\frac{N_{\rm f}N_{\rm c}}{\pi^{2}}T^{4}J_{2}^{\text{cut}}(\beta M,\beta\mu;\beta\Lambda)\,.$$

(27)

As in the previous case, the GL expansion takes the form

$$\Omega=\alpha_{2}M^{2}+\alpha_{4}M^{4}+\cdots.$$

(28)

The coefficient $\alpha_{4}$ was previously computed in Ref. Yamamoto:2007ah , but the expression for the generic coefficient $\alpha_{2n}$ has not been calculated explicitly, to the best of our knowledge. These coefficients are determined by the expansion coefficients of $J_{2}^{\text{cut}}(x,y;\lambda)$ in powers of $x$,

$$J_{2}^{\text{cut}}(x,y;\lambda)=c_{2}x^{2}+c_{4}x^{4}+\cdots,$$

(29)

via

$$c_{2n}=\frac{1}{n!}\partial_{x^{2}}^{n}J_{2}^{\text{cut}}(x,y;\lambda)\Big{|}_{x=0}\,,$$

(30)

which we will now derive. Finally, once the $c_{2n}$ coefficients are known, the $\alpha_{2n}$ coefficients follow immediately from Eq. (27),

$$\alpha_{2n}=\frac{\delta_{n1}}{4G}-\frac{N_{\text{f}}N_{\text{c}}}{\pi^{2}}\frac{c_{2n}}{T^{2n-4}}\Big{|}_{y=\beta\mu,\lambda=\beta\Lambda}\,.$$

(31)

The two differences from the BCS case are that we now have $J_{2}$ rather than $J_{0}$, and we now have $y\neq 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 $\lambda$ finite on the $c_{4}$ coefficient in addition to $c_{2}$. Having $y\neq 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

	$\displaystyle\partial_{x^{2}}J_{2}^{\text{cut}}$	$\displaystyle=\int_{0}^{\lambda}dt\sum_{\zeta=\pm 1}\frac{t^{2}}{4\sqrt{x^{2}+t^{2}}}\tanh\left(\frac{\sqrt{x^{2}+t^{2}}+\zeta y}{2}\,\right)$
		$\displaystyle=\int_{0}^{\lambda}dt\,\sum_{\zeta=\pm 1}\frac{t^{2}}{\sqrt{x^{2}+t^{2}}}\sum_{k=0}^{\infty}\frac{\sqrt{x^{2}+t^{2}}+\zeta y}{\bar{\omega}_{k}^{2}+(\sqrt{x^{2}+t^{2}}+\zeta y)^{2}}\,,$		(32)

where $\bar{\omega}_{k}=(2k+1)\pi$, as before. After applying the identity

$$\sum_{\zeta=\pm 1}\frac{1}{b}\cdot\frac{b+\zeta c}{a^{2}+(b+\zeta c)^{2}}=\sum_{\zeta=\pm 1}\frac{1}{b^{2}+(a+i\zeta c)^{2}}\,,$$

(33)

we obtain

$$\partial_{x^{2}}J_{2}^{\text{cut}}=2\operatorname{Re}\sum_{k=0}^{\infty}\int_{0}^{\lambda}dt\frac{t^{2}}{x^{2}+t^{2}+(\bar{\omega}_{k}+iy)^{2}}\,.$$

(34)

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_{\text{conv}}^{2}=|t^{2}+(\bar{\omega}_{k}+iy)^{2}|$. Taking $t\to 0$ then gives $x_{\text{conv}}^{2}=|\bar{\omega}_{k}+iy|^{2}$, the minimum of which is $x_{\text{conv}}^{2}=\pi^{2}+y^{2}$; recalling that $x=\beta M$ and $y=\beta\mu$, this corresponds to $M_{\text{conv}}=\sqrt{\mu^{2}+(\pi T)^{2}}$. However, this heuristic argument is fairly crude, partly because the quantity $|t^{2}+(\bar{\omega}_{k}+iy)^{2}|$ need not be minimized in the limit $t\to 0$. We therefore proceed to the rigorous proof.

Taking $n-1$ more derivatives of Eq. (34) with respect to $x^{2}$, introducing a variable $s=t/(\bar{\omega}_{k}+iy)$, and letting $\lambda\to\infty$, we find

$$\partial_{x^{2}}^{n}J_{2}^{\text{cut}}\Big{|}_{x=0}=2(-1)^{n-1}(n-1)!\operatorname{Re}\sum_{k=0}^{\infty}\frac{1}{(\bar{\omega}_{k}+iy)^{2n-3}}\int_{0}^{\infty}ds\frac{s^{2}}{(s^{2}+1)^{n}}\,.$$

(35)

The integral is easily performed by using Eq. (13), giving

$$\int_{0}^{\infty}ds\frac{s^{2}}{(1+s^{2})^{n}}=\frac{(2n-5)!!}{(2n-2)!!}\frac{\pi}{2}\,.$$

(36)

Notice that the sum in Eq. (35) converges when $n\geq 3$. Identifying the sum with the standard series representation of the $n$th-order polygamma function $\psi^{(n)}(z)=(-1)^{n+1}n!\sum_{k=0}^{\infty}1/(z+k)^{n+1}$, we finally find

$$c_{2n\geq 6}=\frac{(-1)^{n}}{n!2^{n}(2\pi)^{2n-4}(2n-4)!!}\operatorname{Re}\psi^{(2n-4)}\left(\frac{1}{2}+i\frac{y}{2\pi}\right)\,.$$

(37)

For $c_{2}$ and $c_{4}$, it is more convenient to avoid writing $J_{2}$ as a Matsubara sum, starting instead from the expression

$\displaystyle J_{2}^{\text{cut}}$

$\displaystyle=\int_{0}^{\lambda}dt\,t^{2}\left[\sqrt{x^{2}+t^{2}}+\sum_{\zeta=\pm 1}\ln\left(1+e^{-\sqrt{x^{2}+t^{2}}+\zeta y}\right)\right].$

(38)

For $c_{2}$, it is straightforward to calculate

$\displaystyle\partial_{x^{2}}J_{2}^{\text{cut}}\Big{|}_{x=0}$

$\displaystyle=\frac{\lambda^{2}}{4}+\frac{1}{2}\sum_{\zeta=\pm 1}\left[\lambda\ln\left(1+e^{-\lambda+\zeta y}\right)-\operatorname{Li}_{2}\left(-e^{-\lambda+\zeta y}\right)+\operatorname{Li}_{2}\left(-e^{\zeta y}\right)\right],$

(39)

where $\operatorname{Li}_{n}(x)=\sum_{k=1}^{\infty}x^{k}/k^{n}$ is the $n$th polylogarithm. The first two terms under the sum in Eq. (39) vanish as $\lambda\to\infty$, and for the last term we can apply the identity Apostol1976 ,

$$\operatorname{Li}_{n}(-e^{y})+(-1)^{n}\operatorname{Li}_{n}(-e^{-y})=-\frac{(2\pi i)^{n}}{n!}B_{n}\left(\frac{1}{2}+\frac{y}{2\pi i}\right),$$

(40)

where $B_{n}(x)$ is the $n$th Bernoulli polynomial. The result is

$$c_{2}=\frac{1}{4}\left(\lambda^{2}-y^{2}-\frac{\pi^{2}}{3}\right)\,.$$

(41)

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 $\partial_{x^{2}}^{2}$ with $t^{2}$ in the integrand, we can transform them as $\partial_{x^{2}}\to(2t)^{-1}\partial_{t}$. This allows us to take the limit $x\to 0$ and then take derivatives with respect to $t$. Performing only the rightmost derivative, we find

$$\partial_{x^{2}}^{2}J_{2}^{\text{cut}}\Big{|}_{x=0}=\frac{1}{4}\int_{0}^{\lambda}dt\,t\,\partial_{t}\frac{1}{t}\left(1-\sum_{\zeta=\pm 1}\frac{1}{1+e^{t+\zeta y}}\right).$$

(42)

Integrating by parts, taking $\lambda\to\infty$ on the boundary terms, and by applying the formula Gyory ,

$$\lim_{\lambda\to\infty}\left[\ln\left(\frac{\lambda}{2\pi}\right)-\int_{0}^{\lambda}\frac{dt}{t}\left(1-\sum_{\zeta=\pm 1}\frac{1}{1+e^{t+\zeta y}}\right)\right]=\operatorname{Re}\psi\left(\frac{1}{2}+i\frac{y}{2\pi}\right)\,,$$

(43)

where $\psi(z)$ is the digamma function defined by $\psi(z)={d}\ln\Gamma(z)/dz$ with $\Gamma(z)$ the gamma function, we obtain

$$\partial_{x^{2}}^{2}J_{2}^{\text{cut}}\Big{|}_{x=0}\approx\frac{1}{4}\left[1-\ln\left(\frac{\lambda}{2\pi}\right)+\operatorname{Re}\psi\left(\frac{1}{2}+i\frac{y}{2\pi}\right)\right]\,,$$

(44)

for large $\lambda$, and accordingly,

$\displaystyle c_{4}$

$\displaystyle=\frac{1}{8}\left[1-\ln\left(\frac{\lambda}{2\pi}\right)+\operatorname{Re}\psi\left(\frac{1}{2}+i\frac{y}{2\pi}\right)\right]\,.$

(45)

Finally, applying Eq. (31), we have

	$\displaystyle\alpha_{2}$	$\displaystyle=\frac{1}{4G}-\frac{N_{\text{f}}N_{\text{c}}}{\pi^{2}}\frac{1}{4}\left(\Lambda^{2}-\mu^{2}-\frac{\pi^{2}}{3}T^{2}\right)\,,$		(46)
	$\displaystyle\alpha_{4}$	$\displaystyle=-\frac{N_{\text{f}}N_{\text{c}}}{\pi^{2}}\frac{1}{8}\left[1-\ln\left(\frac{\Lambda}{2\pi T}\right)+\operatorname{Re}\psi\left(\frac{1}{2}+i\frac{\mu}{2\pi T}\right)\right]\,,$		(47)
	$\displaystyle\alpha_{2n\geq 6}$	$\displaystyle=-\frac{N_{\text{f}}N_{\text{c}}}{\pi^{2}}\frac{(-1)^{n}}{n!2^{n}(2\pi T)^{2n-4}(2n-4)!!}\operatorname{Re}\psi^{(2n-4)}\left(\frac{1}{2}+i\frac{\mu}{2\pi T}\right)\,.$		(48)

We will show that the expansion of $J_{2}^{\text{cut}}(x,y;\lambda)$ has radius of convergence $x_{\text{conv}}=\sqrt{\pi^{2}+y^{2}}$. Then, from Eq. (27), it will follow that the GL expansion has convergence radius $M_{\text{conv}}=(x_{\text{conv}}|_{y=\beta\mu})T=\sqrt{\mu^{2}+(\pi T)^{2}}$.

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\geq 6}$ given by Eq. (37). It is straightforward to show that the radius of convergence is at least $\sqrt{\pi^{2}+y^{2}}$. Using the inequality $|\operatorname{Re}(z)|\leq|z|$ for complex $z$, we have

$$\big{|}c_{2n}\big{|}\leq\frac{1}{n!\,2^{n}(2n-4)!!}\frac{1}{(2\pi)^{2n-4}}\bigg{|}\psi^{(2n-4)}\left(\frac{1}{2}+i\frac{y}{2\pi}\right)\bigg{|}\eqqcolon b_{2n}\,,$$

(49)

for $2n\geq 6$.

Using the identity $\psi^{(n)}(x)=(-1)^{n+1}n!\,\zeta(n+1,x)$, where $\zeta(s,a)=\sum_{k=0}^{\infty}1/(k+a)^{s}$ is the Hurwitz zeta function, one can show that

$$\frac{\psi^{(n)}(x)}{\psi^{(n+2)}(x)}\xrightarrow{n\to\infty}\frac{1}{(n+1)(n+2)}x^{2}$$

(50)

for $\operatorname{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\to\infty$. It follows that

$\displaystyle\lim_{n\to\infty}\Bigg{|}\frac{b_{2n}}{b_{2n+2}}\Bigg{|}$

$\displaystyle=\pi^{2}+y^{2}.$

(51)

Thus, the series for $J_{2}^{\text{cut}}(x,y;\lambda)$ in Eq. (29) is bounded in magnitude by a series with radius of convergence $x_{\text{conv}}=\sqrt{\pi^{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 $\sqrt{\pi^{2}+y^{2}}$, but the proof is more involved. In particular, if we do not take the upper bound of the coefficients using $|\operatorname{Re}(z)|\leq|z|$, then the ratio test does not work, because the quantity $[\operatorname{Re}\psi^{(2n-4)}(\frac{1}{2}+i\frac{y}{2\pi})]\allowbreak/[\operatorname{Re}\psi^{(2n-2)}(\frac{1}{2}+i\frac{y}{2\pi})]$ oscillates erratically. One can instead use the Cauchy–Hadamard theorem; we sketch this proof in Appendix B.1.

Although the results of this section were computed using a simple momentum cutoff $\Lambda$, the radius of convergence is unaffected when using Schwinger proper time regularization. Whereas the former cutoff scheme involves taking the limit $\Lambda\rightarrow\infty$ to obtain parts of $c_{2n}$ and $\alpha_{2n}$ analytically, the latter scheme allows deriving their analytic expressions even without taking such a limit for a cutoff parameter $\Lambda$. As shown in Appendix A, the $\alpha_{2n\geq 6}$ coefficients are almost the same as in Eq. (48), except with small corrections proportional to inverse powers of $\Lambda$. The key fact is that these correction terms decay faster than the $\alpha_{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 $\Lambda=\infty$ made when deriving Eqs. (46)–(48) is only a rough approximation in the NJL setting. The values for the two parameters $\Lambda$ and $G$ are determined by a choice of values for $M_{0}$ and $f_{\pi}$, where $M_{0}$ is the dynamical quark mass and $f_{\pi}$ is the pion decay constant at $T=\mu=0$. Choosing $M_{0}=300\,\text{MeV}$ and $f_{\pi}=88\,\text{MeV}$ in the chiral limit Nickel:2009wj ; Buballa:2014tba , we find $\Lambda=614\,\text{MeV}$ and $G\Lambda^{2}=2.15$ Klevansky:1992qe . With these values, one can numerically find $T_{\rm c}=182\,\text{MeV}$ at $\mu=0$, and $\mu_{\rm c}=311\,\text{MeV}$, where $\mu_{c}$ is the chemical potential at which $M$ vanishes in a first-order transition with $T=0$ fixed. Thus, $\Lambda$ is not especially large on the scale of other characteristic quantities of the system, so taking $\Lambda\to\infty$ is not fully justified. For instance, finding $T_{\rm c}$ by setting $\alpha_{2}=0$ in Eq. (46) gives $T_{\rm c}=165\,\text{MeV}$. Although the series with coefficients given by Eqs. (46)–(48) converges with radius $\sqrt{\mu^{2}+(\pi T)^{2}}$, the value does not converge exactly to $\Omega$, even within this radius.

These issues are avoided when using Schwinger proper time regularization, described in Appendix A, because there the finite size of $\Lambda$ 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\,\text{MeV}$ and $f_{\pi}=88\,\text{MeV}$, we find $\Lambda=634\,\text{MeV}$ and $G\Lambda^{2}=6.02$. We calculated $T_{\rm c}$ over the range $0\leq\mu\leq 312\,\text{MeV}$, and we determined the region in the $\mu$-$T$ plane below $T_{\rm c}$ where $M<M_{\text{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 $\mu=0$ (bottom left panel), the GL solution converges only when $T>92\,\text{MeV}$ (right of the vertical line); both the $6$th- and $20$th-order solutions are very accurate near $T_{\rm c}$, but the $6$th-order solution becomes noticeably inaccurate at $T\lesssim 130\,\text{MeV}$. On the other hand, when $\mu=300\,\text{MeV}$ (bottom right panel), the phase transition is first order, and the $20$th-order solution is a noticeable improvement over the $6$th-order solution across the entire range $0<T<T_{\rm c}$. Moreover, the GL solution converges over this entire range, and hence the $20$th-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_{\text{conv}}=\mu$, and noticing that $M\leq 300\,\text{MeV}$ for $\mu>300\,\text{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 $\sqrt{\mu^{2}+(\pi T)^{2}}$. One could also ask about the radius of convergence in the cutoff method if not using the approximation $\Lambda=\infty$, i.e. if instead one were to calculate the GL coefficients from Eq. (27) numerically at finite $\Lambda$. It is possible to show that in this case the radius of convergence is still given by at least $\sqrt{\mu^{2}+(\pi T)^{2}}$ when $M>0$; we sketch a proof in Appendix B.2.