1 Introduction

There is compelling evidence that Quantum Chromodynamics (QCD) is the correct theory describing the strong interaction between quarks and gluons. QCD is a non-Abelian gauge theory, which has a nontrivial vacuum state and provides a rich phase structure and exciting physical phenomena in various environments [1,2,3,4]. At low energy densities quarks are confined inside hadrons, while they form a weakly coupled plasma in the very high energy limit. Therefore, the strongly interacting QCD matter is expected to undergo a hadron-quark phase transition at high temperature and/or chemical potential. Moreover, the chiral symmetry is spontaneously broken at low temperature/baryon density due to the nonvanishing chiral condensate, while it is expected to restore at high temperature and high density [5, 6]. Lattice calculations show that at small baryon chemical potential the transition is actually a smooth crossover [7,8,9,10,11]. On the other hand, many studies [12,13,14] suggest that at large chemical potential the smooth crossover meets a first order transition at a particular point in the \((T,\mu )\) plane named the critical endpoint (CEP), where the crossover becomes a second order phase transition. The quest for the possible existence of CEP as well as its precise location are still on-going, and have attracted a lot of attention, see e.g. [15,16,17,18,19,20,21].

In QCD, fluctuations of the conserved charges like the baryon and isospin numbers, exhibit critical behavior around a phase transition; these fluctuations are encoded into the corresponding susceptibilities as well as into higher order moments [22,23,24,25]. For example, the chiral susceptibility, which represents the modification of the chiral condensate to a small perturbation of the current quark mass, is often used to describe the chiral phase transition [26,27,28,29,30,31,32]. Besides, it has been suggested that the moments of the baryon number can be used to identify the CEP in the QCD phase diagram [24, 33,34,35,36,37,38,39,40,41,42,43,44].

In the two-flavor case the quark chemical potentials \(\mu _f\) can be expressed in term of the baryon chemical potential \(\mu _B=N_c(\mu _u+\mu _d)/2\), and the isospin chemical potential \(\mu _I=\mu _u-\mu _d\)Footnote 1 with \(N_c=3\) the number of colors. Along the \(\mu _B\) axis, the knowledge from first principle lattice simulation on the thermodynamics of strongly interacting matter is limited, due to the notorious sign problem [47]. Thus, effective field theories of QCD as well as phenomenological models are necessary, and they might provide powerful tools for us to have a better understanding of QCD physics in the nonperturbative regime. On the other hand, lattice simulations are feasible for \(\mu _B=0\) and \(\mu _I\ne 0\), therefore it is possible to study QCD at finite \(\mu _I\) by means of first principle calculations and compare these with the predictions of effective models. This is an useful exercise to test the effectiveness of QCD models. For example, an analysis based on the leading order (LO) [48, 49] and the next-to-leading order (NLO) [50, 51] chiral perturbation theory (CHPT) shows that a second order phase transition from the normal to pion superfluid phase is expected to happen at \(\mu _I^c(T=0)=m_\pi \), where \(m_\pi \) is the pion mass. This has been confirmed by lattice simulations [52,53,54,55] as well as by the Nambu−Jona-Lasinio (NJL) model analytically [56,57,58].

Many aspects of pion condensate have been extensively studied in the effective field theories [59,60,61,62,63,64,65,66], lattice simulations [45, 67,68,69,70] and effective models [71,72,73,74,75,76,77,78,79,80,81,82] (for a review, see [83]). We also mention that a new type of compact star made of a charged pion condensate has been proposed recently [63, 84,85,86]. While we refer to the original articles for the detailed picture for both two-flavor and three-flavor quark matter, we remind here the well established picture for the two-flavor case with finite bare quark masses. At zero temperature there is a second order phase transition from the normal phase, in which the pion condensate is vanishing, to the phase with nonzero pion condensate: this phase transition happens when \(|\mu _I|=m_\pi \); increasing \(|\mu _I|\) there is eventually another second order phase transition back to the normal phase. At finite temperature the gap between the two critical \(|\mu _I|\) becomes smaller and eventually the window for the pion condensate phase closes: at very high temperature there is room for the normal phase only, see for example Fig. 4 of [78].

In this paper, we study several susceptibilities of isospin imbalanced QCD matter at vanishing \(\mu _B\) and zero and/or finite temperature, which to the best of our knowledge has not yet been explicitly discussed in the literature. As mentioned above, doing this study is useful to test how effective models predictions compare with first principle calculations. The thermodynamic properties of pion superfluid phase at vanishing temperature has been studied at the lowest order of CHPT at zero temperature [64] and it has been found a good agreement with the lattice QCD simulations for \(\mu _I\lesssim 2m_\pi \), but also that the LO approximation breaks down at higher \(\mu _I\). The NJL model can be employed in a larger range of chemical potential and temperature with respect to LO CHPT and it is the purpose of this work to show how this model describes the thermodynamics and the fluctuations for the transition to the pion superfluid phase. We will use the two-flavor NJL model in this article for simplicity, leaving the three-flavor case to a future study.

The article is organized as follows. In Sect. 2, we give a brief introduction for the theoretical framework of NJL model at finite \(\mu _I\) and temperature. In Sect. 3, after comparing the NJL model results for several thermodynamic quantities with LO CHPT calculation and/or lattice data at vanishing temperature, we give the numerical results for several susceptibilities of our interest at finite \(\mu _I\) and temperature. Finally, we draw our conclusions in Sect. 4

2 The NJL model

The standard NJL model Lagrangian density for two flavors of quarks is given by [87,88,89,90,91,92]

$$\begin{aligned} \mathcal {L}=\bar{q}(i\gamma ^\mu \partial _\mu -m+\hat{\mu }\gamma _0)q+G[(\bar{q}q)^2+(\bar{q}i\gamma _5\tau q)^2], \end{aligned}$$
(1)

where q represent the quark fields, \(\tau \) the Pauli matrices in flavor space, G the four-fermion interaction coupling, m the degenerate quark mass of up and down quarks and \(\hat{\mu }=\mathrm {diag}\{\mu _u,\mu _d\}\) denotes the diagonal matrix for quark chemical potentials.

In the grand canonical ensemble the thermodynamics of strongly interacting matter can be considered as a function of the isospin and baryon chemical potentials and temperature. In this article we limit ourselves to the case \(\mu _B=0\) which leads to \(\mu _u=\mu _I/2\), \(\mu _d = -\mu _I/2\). For simplicity, we perform the calculations in the mean field approximation, i.e.,

$$\begin{aligned}&(\bar{q} q)^2 \approx 2(\bar{q} q)\langle \bar{q} q\rangle -\langle \bar{q} q\rangle ^2, \end{aligned}$$
(2)
$$\begin{aligned}&(\bar{q}i\gamma _5\tau q)^2 \approx 2(\bar{q}i\gamma _5\tau q)\langle \bar{q} i\gamma _5\tau q\rangle -\langle \bar{q} i\gamma _5\tau q\rangle ^2. \end{aligned}$$
(3)

We introduce the chiral condensates as

$$\begin{aligned} \langle \bar{q}q\rangle =\sigma =\sigma _u+\sigma _d, \end{aligned}$$
(4)

with \(\langle \bar{u}u\rangle =\sigma _u\) and \(\langle \bar{d}d\rangle =\sigma _d\). At \(\mu _I\ne 0\) a charged pion condensate is also expected: in order to account for this we introduce

$$\begin{aligned} \pi ^\pm = \langle \bar{q}i\gamma _5\tau _\pm q \rangle \equiv \Pi e^{\pm i \theta }, \end{aligned}$$
(5)

with \(\theta \) a real number and \(\tau _\pm =\tau _1 \pm i\tau _2\); the effective potential does not depend on \(\theta \) so we choose \(\theta =0\) in the above equation and we are left with the pion condensates, \(\Pi \), in the \(\tau _1\)-direction of isospin, namely

$$\begin{aligned} \Pi = \langle \bar{q}i\gamma _5\tau _1 q \rangle . \end{aligned}$$
(6)

The Lagrangian density can thus be written as

$$\begin{aligned} \mathcal {L}= & {} \bar{q}\big [i\gamma ^\mu \partial _\mu -M_q+\hat{\mu }\gamma ^0 \nonumber \\&+2G \Pi i\gamma ^5\tau _1 \big ]q +G(\sigma ^2+\Pi ^2) \end{aligned}$$
(7)

with the effective quark mass defined as \(M_q=m-2G\sigma \). The thermodynamic potential is

$$\begin{aligned} \Omega= & {} G(\sigma ^2+\Pi ^2)-\frac{2N_c}{\beta }\int \frac{d^3p}{(2\pi )^3} \nonumber \\&\times \Big \{\ln \left[ (1+e^{-\beta E_p^+})(1+e^{\beta E_p^+})\right] \nonumber \\&+\ln \left[ (1+e^{-\beta E_p^-})(1+e^{\beta E_p^-})\right] \Big \}, \end{aligned}$$
(8)

where \(\beta =1/T\) and

$$\begin{aligned} E_p= & {} \sqrt{p^2+M^2}, \end{aligned}$$
(9)
$$\begin{aligned} E_p^{\pm }= & {} \sqrt{\big (E_p\pm \mu _I/2\big )^2+4G^2\Pi ^2}. \end{aligned}$$
(10)

The ground state is given by the values of \(\sigma \), \(\Pi \) that minimize \(\Omega \): this is equivalent to require that \(\sigma \) and \(\Pi \) are solution of the gap equations,

$$\begin{aligned} \frac{\partial \Omega }{\partial \sigma }=\frac{\partial \Omega }{\partial \Pi }=0. \end{aligned}$$
(11)

Given \(\Omega \) at the ground state it is straightforward to obtain the isospin number density, the pressure and the energy density as

$$\begin{aligned} n_I= & {} -\frac{\partial \Omega }{\partial \mu _I}, \end{aligned}$$
(12)
$$\begin{aligned} P= & {} -(\Omega -\Omega _0), \end{aligned}$$
(13)
$$\begin{aligned} E= & {} -P-T\frac{\partial \Omega }{\partial T}-\mu _I\frac{\partial \Omega }{\partial \mu _I}, \end{aligned}$$
(14)

where \(\Omega _0\) is the thermodynamic potential at \(T=\mu _I=0\).

The two-flavor NJL model has three parameters: the current quark mass \(m=0.005\) GeV, the four-fermion coupling strength \(G=5.01\) GeV\(^{-2}\), and the hard three-momentum cutoff \(\Lambda =0.653\) GeV, which are fixed by reproducing the empirical values of the pion mass \(m_\pi =0.134\) GeV, the physical pion decay constant \(f_\pi =0.093\) GeV, and the chiral condensate \(\sigma _0=2(-0.25~\text {GeV})^3\) in vacuum [93].

Before going on it is useful to remind the results of LO CHPT for the pressure and the energy density in the pion superfluid phase, namely

$$\begin{aligned} P^{\text {LO}}= & {} \frac{f_\pi ^2}{2\mu _I^2}(\mu _I^2-m_\pi ^2)^2,\end{aligned}$$
(15)
$$\begin{aligned} E^{\text {LO}}= & {} \frac{f_\pi ^2}{2\mu _I^2}(\mu _I^4+2\mu _I^2m_\pi ^2-3m_\pi ^4), \end{aligned}$$
(16)

where \(f_\pi \) is the pion decay constant. The result for isospin density will be also useful in the following section:

$$\begin{aligned} n_I=\frac{\partial P^{\text {LO}}}{\partial \mu _I}=\frac{f_\pi ^2}{\mu _I^3}(\mu _I^4-m_\pi ^4). \end{aligned}$$
(17)

3 Results and discussion

3.1 Isospin number density, energy density, and trace anomaly

Fig. 1
figure 1

The normalized isospin number density as a function of \(\mu _I/m_\pi \) at fixed \(T=0\). The blue dashed line corresponds to the result obtained from the NJL model. The result from the CHPT (purple dotted line), i.e., Eq. (17), and recent lattice data [85] (magenta circles) are also included for comparison

Full size image

In Fig. 1, we plot the isospin number density normalized to \(m_\pi ^3\) at zero temperature as a function of \(\mu _I/m_\pi \). For comparison we also include the results from LO CHPTFootnote 2 (purple dotted line) as well as recent lattice data (magenta circles) [85]. The results obtained from these methods agree quantitatively with each other in the considered range of \(\mu _I\). When \(\mu _I<m_\pi \) the isospin number density is zero, which corresponds to the normal phase of the system. In the pion superfluid phase with \(\mu _I> m_\pi \), the isospin number density becomes nonzero, which increases monotonically with \(\mu _I\).

Fig. 2
figure 2

The normalized energy density \(E/\mu _I^4\) and trace anomaly \((E-3P)/\mu _I^4\) as functions of \(\mu _I/m_\pi \) at zero temperature. The result from CHPT at zero temperature is also included for comparison

Full size image

In Fig. 2, we plot the energy density normalized to \(\mu _I^4\) as a function of \(\mu _I/m_\pi \) at \(T=0\). The results from the LO CHPT (purple dotted line) in Eq. (16) are also shown for comparison. The result from the NJL model agrees with the one from CHPT at zero temperature though the energy density from the former is slightly larger than that from the latter at high \(\mu _I\). The quantity \(E/\mu _I^4\) becomes nonzero and positive when \(\mu _I\) exceeds the critical value \(\mu _I^c(T=0)=m_\pi \). In particular, it develops a peak at \(\mu _I^{\text {peak}}\simeq 1.274m_\pi \). Recently, a calculation of this peak based on the LO CHPT has been done in Ref. [64]; moreover, this quantity has also been computed on the lattice:

$$\begin{aligned} \mu _I^{\text {peak}}\simeq {\left\{ \begin{array}{ll} 1.274m_\pi , &{}\text {NJL},\\ 1.276m_\pi , &{}\text {CHPT}~\text {[16]},\\ \{1.20,1.25,1.275\}m_\pi , &{}\text {Lattice data}~\text {[55]}. \end{array}\right. } \end{aligned}$$
(18)

In the above equation the three lattice results quoted have been obtained with three different lattice volumes, namely \(L^3=\{16^3,~20^3,~24^3\}\), respectively. The extrapolation of the lattice results to the continuum limit gives \(\mu _I^{\text {peak}}=1.30(7)m_\pi \). We notice the nice agreement of the NJL model with LO CHPT and lattice data.

In Fig. 2, we also show the normalized trace anomaly, \((E-3P)/\mu _I^4\), as a function of \(\mu _I/m_\pi \) at \(T=0\). Again we notice the agreement between the NJL model and the LO CHPT results in the entire isospin chemical potential range considered. In more detail, we observe that the blue dashed line (lower one), which represents the result calculated in the NJL model at zero temperature, keeps zero at small \(\mu _I\), and becomes nonzero and positive when \(\mu _I\) slightly larger than \(m_\pi \). Similar to the behavior of \(E/\mu _I^4\), the normalized trace anomaly also develops a peak at an intermediate \(\mu _I\) although the peak of the former appears at a slightly smaller value of \(\mu _I\) with respect to the latter. As \(\mu _I\) increases inside the domain of the pion condensed phase, the normalized trace anomaly decreases and becomes negative. In fact, the point \(\bar{\mu }_I\) fulfilling the conformal relation \(E-3P=0\) can be evaluated in the NJL model at zero temperature as

$$\begin{aligned} \bar{\mu }_I \simeq 1.754m_\pi , \end{aligned}$$
(19)

which is in good agreement with the LO CHPT result given by \(\bar{\mu }_I=\sqrt{3}m_\pi \) [64]. This point separates the \(E>3P\) and \(E<3P\) regions. In fact, this value \(\bar{\mu }_I \simeq 1.754m_\pi \) is very close to the value of \(\mu _I\) corresponding to the BEC-BCS crossover, which is estimated to be \(1.702m_\pi \) at zero temperature according to the analyses in Refs. [95,96,97,98]. However, although these two values almost coincide with each other, whether \(\bar{\mu }_I\) can be identified with the isospin chemical potential corresponding to the BEC-BCS crossover or not need to be investigated in more detail in the future.

3.2 Chiral, pion, and isospin susceptibilities

Fig. 3
figure 3

The variation behaviors of chiral susceptibility \(\chi _\sigma \), pion susceptibility \(\chi _\pi \), and isospin number susceptibility \(\chi _I\) with respect to \(\mu _I/m_\pi \) for several values of the temperature

Full size image

The transition to the pion condensate phase can be well identified by means of the susceptibilities. In principle, they can also be measured on the lattice since there is no sign problem at finite \(\mu _I\) due to the real and positive fermionic determinant. Chiral condensate is the order parameter of spontaneous chiral symmetry breaking, while pion condensate indicates the spontaneous isospin symmetry breaking. The chiral susceptibility, \(\chi _\sigma \), which corresponds to the zero-momentum projection of the scalar propagator and encodes all fluctuations of the order parameter, is defined as [28, 30, 31, 99]

$$\begin{aligned} \chi _\sigma =-\frac{\partial \sigma }{\partial m}. \end{aligned}$$
(20)

We can also define the pion susceptibility, which is obtained as the first derivative of pion condensate with respect to the current quark mass, i.e.,

$$\begin{aligned} \chi _\pi =\frac{\partial \Pi }{\partial m}. \end{aligned}$$
(21)

We also define the isospin number susceptibility as

$$\begin{aligned} \chi _I=\frac{d n_I}{d\mu _I}. \end{aligned}$$
(22)

Note that we write the derivative operator above as a total derivative instead of a partial one because one might take into account the fact that the condensates may have a dependence on \(\mu _I\).

The chiral, pion, and isospin susceptibilities as functions of \(\mu _I/m_\pi \) for four values of the temperature are illustrated in Fig. 3. The temperatures are chosen as \(T=0\) (black short-dashed line), \(T=0.1\) GeV (blue dashed line), \(T=0.15\) GeV (red dotted line), and \(T=0.22\) GeV (magenta dot-dashed line). We notice that the zero temperature susceptibilities exhibit a discontinuity at \(\mu _I=\mu _I^c(T)\), signaling the boundary of a second order phase transition. At nonzero temperatures, e.g., \(T=0.1\) GeV and 0.15 GeV, the situation is similar except that the discontinuity shifts to larger \(\mu _I\). For \(T=0.22\) GeV instead the susceptibilities are continuous in the whole range of \(\mu _I\), which can be easily understood with the absence of a transition to the pion condensed phase in agreement with [56]. We notice that the chiral susceptibility stays finite in the pion condensed phase, while the pion susceptibility vanishes in the normal phase: this is obviously related to the fact that the bare quark mass is finite, so chiral symmetry is always broken explicitly and this leads to a finite chiral condensate. In the chiral limit this would not happen and the chiral susceptibility would vanish in the pion condensed phase.

Fig. 4
figure 4

The variation behaviors of the chiral susceptibility with respect to temperature at several values of \(\mu _I\)

Full size image

In Fig. 4, we show the thermal behavior of the chiral susceptibility at several values of \(\mu _I\). The chiral susceptibility at \(\mu _I=0\), which is denoted by the black short-dashed line, firstly grows up smoothly at low temperature, then it develops a smooth peak at \(T=0.2\) GeV in correspondence of the chiral crossover. At higher temperatures, the chiral susceptibility decreases with temperature. We also include several cases with \(\mu _I>\mu _I^c(T)\) for comparison. For \(\mu _I=0.3\) GeV, which is represented by the blue dashed line, the chiral susceptibility has a much larger value at \(T=0\) compared to the other cases. This line first slowly increases with temperature and reaches a smooth peak at \(\mu _I/m_\pi =0.175\), signaling a crossover. After the peak the line drops and show a kink at \(\mu _I/m_\pi =0.183\), in correspondence of the transition to the normal phase. Similarly, the lines for \(\mu _I=0.5\) GeV and 0.7 GeV, which are represented by red dotted and orange dot-dashed lines respectively, also show a modest kink peak at the critical temperature.

Fig. 5
figure 5

The variation behaviors of the pseudoscalar susceptibility against temperature for several values of \(\mu _I\)

Full size image

In Fig. 5, we plot \(\chi _\pi \) versus temperature for several values of \(\mu _I\), i.e., \(\mu _I=0.2,~0.3,~0.5,~0.7\) GeV. We do not illustrate the result for \(\chi _\pi \) with \(\mu _I\) smaller than the critical \(\mu _I^c(T)\) because in the normal phase \(\chi _\pi =0\). For \(\mu _I=0.2\) GeV, which is labeled by the black short-dashed line, \(\chi _\pi \) increases very smoothly at low temperature, and then grows up rapidly when T is close to the critical temperature; clearly it vanishes above this temperature because the pion condensate is zero there.

Fig. 6
figure 6

The isospin number susceptibility \(\chi _I\) as a function of the temperature for several values of \(\mu _I\)

Full size image

In Fig. 6 We plot \(\chi _I\) versus the temperature for several values of \(\mu _I\) . For the black short-dashed line, which represents the result obtained at \(\mu _I=0\), the isospin number susceptibility increases smoothly and monotonously with temperature. At higher \(\mu _I\), e.g., \(\mu _I=0.3\) GeV (blue dashed line), \(\mu _I=0.5\) GeV (red dotted line), and \(\mu _I=0.7\) GeV (magenta dot-dashed line), \(\chi _I\) keep almost unchanged in the entire considered temperature range expect for the cusp at the phase transition point.

4 Conclusions

In this work, we have studied the isospin chemical potential and temperature dependence of several thermodynamic quantities, in particular of several susceptibilities, of isospin-imbalanced QCD matter at vanishing baryon chemical potential within the two-flavor NJL model. We have compared the isospin number density at zero temperature with recent lattice data as well as with LO CHPT, finding a good agreement with the latter results, see Fig. 1. In addition, the result for the normalized energy density \(E/\mu _I^4\) as a function of \(\mu _I/m_\pi \) is also in fair agreement with that from LO CHPT at zero temperature. In particular, the location of the peak of \(E/\mu _I^4\) as well as the trace anomaly computed within the NJL model agree with LO CHPT and/or lattice simulations: we have found \(\mu _I^{\text {peak}}\simeq 1.274m_\pi \) for the peak of the energy density and \(\bar{\mu }_I \simeq 1.754m_\pi \) for the condition \(E-3P=0\).

We have then considered the chiral, pion and isospin susceptibilities at finite temperature. These are interesting quantities because they allow to define the critical lines in the \((T,\mu _I)\) plane and can also be computed on the lattice, so they can be used as a test of the effective models of QCD. The qualitative behavior of these susceptibilities is in agreement with naive expectations: in particular, the chiral susceptibility is always finite because of the nonvanishing bare quark mass that induces a finite chiral condensate even at large \(T/\mu _I\). On the other hand, the pion susceptibility vanishes in the normal phase and is nonzero in the pion condensed phase. The transition from the normal phase to the pion condensed phase is of the second order with a divergent pion susceptibility. The isospin number susceptibility also keep finite and almost unchanged in the normal phase except that it is zero for the zero temperature case, and shows a discontinuity at the phase transition.

In this study we have ignored the role of baryon and strangeness chemical potentials [64]. These external sources are known to play an important role and we leave the study of the more complete problem to future works.