Thermodynamics from the quark condensate

We present a method to compute thermodynamic quantities within functional continuum frameworks which is independent of the employed truncation. As a proof of principle, we first apply it to a Nambu-Jona-Lasinio model in mean-field approximation. Then, we use the method with solutions obtained from a coupled set of truncated Dyson-Schwinger equations for the quark and gluon propagators of (2+1)-flavor quantum chromodynamics in Landau gauge to obtain the pressure, entropy density, energy density, and interaction measure across the phase diagram of strong-interaction matter. We also discuss the limitation of the proposed method.


I. INTRODUCTION
The thermal properties of strong-interaction matter described by the theory of quantum chromodynamics (QCD) at nonvanishing temperature and density are subject to both experimental and theoretical efforts [1,2]. Future as well as already operating heavy-ion-collision facilities such as CBM at FAIR/GSI, NICA at JINR, and RHIC at BNL aim to probe the phase structure of QCD [3][4][5][6][7]. In order to understand and explain the various phases and thermodynamic properties of QCD one needs the pressure, entropy density, and energy density-in short the equation of state (EoS)-as a function of temperature and chemical potential. More generally, the knowledge of the EoS of matter under a given physical environment is important for numerous reasons ranging from hydrodynamic simulations in the context of heavy-ion collisions to astrophysical issues like supernovae or neutron stars; see, e.g., Refs. [8,9] for review articles.
On the theoretical side and at vanishing chemical potential, first-principle lattice-regularized QCD provides the following well-established picture: The confined hadronic low-temperature phase characterized by dynamical chiral symmetry breaking is connected through an analytic crossover to the deconfined high-temperature phase of the quark-gluon plasma with (partially) restored chiral symmetry [10][11][12][13][14]. Furthermore, the full continuumextrapolated EoS is available [15,16]. The situation is fundamentally different at nonvanishing (real) chemical potential. There, lattice QCD is hampered by the sign problem. Results for the EoS are limited to rather small chemical potential and usually obtained via extrapolation from imaginary chemical potential or using the Taylorexpansion technique [17][18][19][20]. Thus, other approaches are necessary to complement the lattice calculations.
In this work we aim to compute thermodynamic quantities within the DSE approach. To this end, one needs basically the thermodynamic potential. In contrast to the FRG where solving its flow equation yields directly the thermodynamic potential, accessing this quantity within the DSE framework is extremely difficult and limited to simple models [39][40][41]. Generally speaking, this is due to the fact that the DSE approach starts with the first derivative of the thermodynamic potential and an integration is needed to get hold of the potential itself. Unfortunately, this integration is only possible for certain truncations. It is thus desirable to develop a truncation-independent way to calculate thermodynamic quantities from DSEs. Purpose of this work is to introduce such a method.
The remainder of this paper is organized as follows: In Sec. II we detail the aforementioned method and apply it to a Nambu-Jona-Lasinio (NJL) model in Sec. III as a testing ground to gauge the method's effectiveness. In Sec. IV we summarize our DSE framework which solutions serve as input to obtain the pressure, entropy density, energy density, and interaction measure in (2 + 1)-flavor QCD. In Sec. V we present and discuss our results for these quantities and finally conclude in Sec. VI.

II. THERMODYNAMICS AND QUARK CONDENSATE
The fundamental quantity for QCD thermodynamics at nonzero temperature T and quark chemical potential µ is the thermodynamic potential with the grand-canonical partition function Z. Here, V is the volume of the system and we consider only one light flavor first. Thermodynamic quantities like pressure (p), arXiv:2012.04991v1 [hep-ph] 9 Dec 2020 entropy density (s), and number density (n) follow from the standard relations Furthermore, a Legendre transform of the pressure yields the energy density and with that one defines the interaction measure It is related to the trace of QCD's energy-momentum tensor (hence also referred to as the trace anomaly) and measures the deviation of the EoS from the one of an ideal gas given by ε = 3p. In addition to temperature and chemical potential, the current-quark mass m can be seen as an additional variable the thermodynamic potential depends on. It appears as an external source for the field bilinearψψ in the QCD action and the quark condensate is obtained via In principle, this relation can be inverted to obtain the thermodynamic potential as an integral of the quark condensate with respect to the current-quark mass, i.e., (6) Unfortunately, this relation is not suitable for an actual calculation since the thermodynamic potential and the quark condensate are both divergent. The divergence is caused by the vacuum contribution contained in Ω and is even present in the noninteracting theory [42]. Since the divergence is independent of temperature and chemical potential, suitable derivatives of the potential and condensate are expected to be finite. 1 Therefore, differentiating Eq. (6) with respect to T yields the welldefined (divergence-free) equation In order to use this relation in practical calculations, we have to specify the integral boundaries. We set the lower one to the physical current-quark mass, m 1 = m, and send the upper one to infinity, m 2 → ∞. An infinitely heavy quark freezes out of the system and does not contribute to thermodynamics. The corresponding entropy density is then simply the one of pure Yang-Mills theory: s(T, µ; m 2 → ∞) = s YM (T ). Thus, our final expression for the entropy density reads If gluons are no active degrees of freedom, like, e.g., in the NJL model, the Yang-Mills contribution is set to zero. For QCD, s YM is taken from the lattice [43,44]. Note that Eq. (8) implies ∂s/∂m = −∂ ψ ψ /∂T , which is nothing but the Maxwell-like relation Having the entropy density at hand, the pressure at vanishing chemical potential follows thermodynamically consistent from and an additional integration over the number density yields the pressure at nonvanishing chemical potential, (11) Analogously to s YM , the value p(T 0 , 0) at a reference temperature T 0 is treated as an input parameter and taken from the lattice [15,16].
In the case of 2 + 1 flavors we have to consider two degenerate light quarks with m u = m d , µ u = µ d and a heavier strange quark with mass m s m u . To calculate the entropy density we integrate within the Columbia plot as follows: first the up-quark mass from the physical point to the physical strange-quark mass and then both masses to infinity (see Fig. 1). With {µ} = (µ u , µ s ), the generalization of Eq. (8) reads We would like to emphasize that Eqs. (8) and (12) are obtained without any approximation and are therefore exact. Only the quark condensate as a function of the quark mass (at fixed T and {µ}) is needed to compute the entropy density. This renders Eqs. (8) and (12) quite general and not constraint to certain approaches-they are applicable as soon as the quark condensate is available. This is particularly useful within the framework of DSEs where accessing the thermodynamic potential is extremely difficult and limited to simple models with rainbow-ladder-like truncations [39][40][41]. The method presented here allows to bypass these limitations in order to compute thermodynamic quantities regardless of the chosen truncation. 2

III. NJL-MODEL STUDY
To show that the method described in the previous section works effectively, we use a two-flavor NJL model [45,46] in mean-field approximation where the thermodynamic potential can be computed analytically.
The two-flavor NJL Lagrangian in mean-field approximation is given by [47][48][49] where m is the current-quark mass, M the constituentquark mass, and G denotes the coupling constant. The thermodynamic potential is thus simply the noninteracting one [42] shifted by a field-independent term, i.e., with E k = √ k 2 + M 2 ; N f = 2 and N c = 3 denote the number of flavors and colors, respectively. We regularize the divergent vacuum integral with a sharp threemomentum cutoff, |k| ≤ Λ, but leave the convergent medium contribution unaltered. The physical constituentquark mass is obtained by minimizing the potential, and the quark condensate reads ψ ψ Finally, the model is complete once the parameters are fixed. We use m = 5.6 MeV, Λ = 587.9 MeV, and G = 2.44/Λ 2 . These values were determined in Ref. [49] to yield a pion mass and decay constant of m π = 135 MeV and f π = 92.4 MeV in vacuum. The resulting constituentquark mass in vacuum is M vac = 400 MeV.
We are now in a position to compute the entropy density first directly from the thermodynamic potential (14) using Eqs. (2) and second via the method described in Sec. II, i.e., by means of Eq. (8). As mentioned earlier, s YM = 0 since gluons are no active degrees of freedom in the NJL model. We find that both results cannot be distinguished by the eye and show the relative error between them in Fig. 2. It is smaller than 0.05 % across the whole covered temperature range. Thus, we are confident that our method to obtain the entropy density from the quark condensate is able to yield reliable results in the Dyson-Schwinger approach, too.

IV. DYSON-SCHWINGER EQUATIONS
In the following we briefly summarize our functional framework of Dyson-Schwinger equations used to determine the dressed quark propagator. From it the quark condensate and eventually the entropy density are obtained. We use the same setup as in our previous work [26] and solve a coupled set of Landau-gauge Dyson-Schwinger equations where the back reaction of the quarks onto the Yang-Mills sector is explicitly taken into account. The dressed quark propagator S f for a flavor f at nonzero temperature T and quark chemical potential µ f is the solution of the DSE 3 which is depicted in Fig. 3. Here, q = (ω q , q) is the fourmomentum withω q = ω q + iµ f and fermionic Matsubara frequencies ω q = (2 q + 1)πT , q ∈ Z. Furthermore, Z f 2 and Z f m denote the wave function and mass renormalization constants; m f is the renormalized current-quark mass. The inverse dressed quark propagator is given by with scalar dressing functions C f , A f , and B f . They carry the nonperturbative information, thus having a nontrivial momentum dependence, and depend on temperature and chemical potential as well. A further dressing function corresponding to the tensor structure γ 4 γ ·q is in principle possible but qualitatively negligible [50]. The explicit form of the quark self-energy appearing in Eq. (17) reads with the dressed gluon propagator D µν , dressed quarkgluon vertex Γ f σ , strong coupling constant g, ghost renormalization constantZ 3 , and loop momentum k = (ω k , k). The prefactor 4/3 stems from the color trace.
In order to solve the quark DSE self-consistently for the dressed quark propagator, we need to specify the dressed gluon propagator and quark-gluon vertex. The truncation used in this work evolved gradually (see Ref. [37] and references therein) and is characterized as follows. First, we use temperature-dependent lattice data for the quenched gluon propagator [51,52] as input and incorporate unquenching effects by explicitly evaluating the quark-loop diagram for each of the N f = 2 + 1 quark flavors considered here. This results in the DSE for the unquenched gluon propagator as shown in Fig. 4. Consequently, the quark and gluon DSEs are nontrivially coupled and need to be solved simultaneously. This construction allows for a dependence of the gluon on temperature and chemical potential controlled by QCD dynamics rather than modeling. Second, we use an ansatz for the dressed quark-gluon vertex motivated by its known perturbative running in the ultraviolet combined with an approximate form of the Slavnov-Taylor identity in the infrared based on the Ball-Chiu vertex construction [53]. Since our setup is identical to the one used recently, we shall not repeat explicit expressions regarding the truncation for the sake of brevity and refer the reader to Ref. [26].
The free parameters of the truncation are the infrared strength of the vertex ansatz and the quark masses. They are fixed to yield a pseudocritical chiral transition temperature at vanishing chemical potential of T c = 156 MeV, defined by the inflection point of the light quark condensate with temperature, in agreement with lattice results [12,13,18,20,54]. We work in the isospin-symmetric limit m u = m d , µ u = µ d ≡ µ and choose µ s = 0 for simplicity. This implies that the baryon chemical potential is given by µ B = 3µ.
Finally, after solving the coupled set of quark and gluon DSEs for the propagators, the corresponding quark condensate is obtained via

V. RESULTS AND DISCUSSION
For the sake of completeness, we first recapitulate the phase diagram obtained with the DSE setup described in Sec. IV. It is shown in Fig. 5 and we refer the reader to Ref. [26] for more details. The chiral crossover line (dashed black) starting at T c (µ = 0) = (156 ± 1) MeV becomes steeper with increasing chemical potential and terminates in a second-order CEP (solid black circle) located at (T CEP , µ CEP ) = (119 ± 2, 165 ± 2) MeV. This corresponds to a ratio of (µ B /T ) CEP ≈ 4.2. The errors are purely numerical. Beyond the CEP we find the coexistence region of a first-order transition (shaded area) bounded by spinodals (solid black).
We now focus on our thermodynamic results obtained with the method described in Sec. II used with condensate data computed from the (2 + 1)-flavor DSE framework summarized in the previous section. Starting point is the entropy density obtained via Eq. (12). Lattice simulations of pure Yang-Mills theory show that there is basically no sizable contribution to s YM /T 3 for T T YM c ≈ 270 MeV [43,44], i.e., for temperatures around and below the critical Yang-Mills temperature. Thus, to a good approximation we set s YM = 0 in Eq. (12) for the temperature range covered in this work.

A. Zero chemical potential
Our result for the entropy density (scaled to T 3 ) at vanishing chemical potential is shown in the upper diagram of Fig. 6 (solid black line) compared to results from lattice QCD [15,16] (colored symbols). Up to T ≈ 175 MeV, it is a monotonically increasing function of temperature and the agreement with lattice data is satisfying. Beyond that temperature the entropy density starts to decrease. This unphysical behavior can be traced back to a deficiency in our vertex ansatz and became apparent already in the calculation of quark and baryon number fluctuations [26]. At the moment we take only the leading Dirac structure γ σ of the dressed quark-gluon vertex into account. However, the full vertex Γ f σ contains twenty-four different Dirac tensor structures (in Landau gauge) with half of them reacting strongly to the (partial) restoration of chiral symmetry around and above the pseudocritical chiral transition temperature. These terms are missing in our current setup. Their inclusion would cause a continuous weakening of the quark-gluon interaction for T T c , thereby resolving the issue of a decreasing entropy density at high temperatures. 4 But since the focus of this article is on the method presented in Sec. II, we leave our setup unchanged. It yields, however, satisfying results in the temperature range 100-160 MeV, i.e., below and around the crossover temperature. FIG. 6. Entropy density (top), pressure (center), and energy density (bottom) at vanishing chemical potential. The lattice data is taken from Refs. [15,16].
From the entropy density we obtain the pressure via integration, see Eq. (10), and use p(T 0 )/T 4 0 = 0.242 at T 0 = 110 MeV [15]. The result depicted in the middle diagram of Fig. 6 is in good agreement with the lattice but starts to deviate for T 185 MeV. This is inherited from the erroneous high-temperature behavior of the entropy density. The pressure saturates at p/p SB ≈ 0.3, where is the Stefan-Boltzmann pressure of an ideal gas of massless quarks and gluons. Finally, combining the entropy density and pressure we obtain the energy density, Eq. (3), shown in the lower diagram of Fig. 6. Since it is a combination of s and p, the agreement with lattice results is reasonable for temperatures below and around T c while a decreasing behavior stemming from the entropy density is found at high temperatures. Ignoring the high-temperature artifacts, one can also define the pseudocritical chiral transition temperature as the inflection point of, e.g., the pressure with temperature. We find T (p) c = 157 MeV consistent with 156 MeV obtained from the light quark condensate.

B. Nonzero chemical potential
We now turn to nonvanishing chemical potential and show our results in Fig 7. The entropy density is depicted in the upper left diagram as a function of temperature for different chemical potentials starting from zero up to the critical-endpoint value µ CEP = 165 MeV. A bulge develops around the pseudocritical chiral transition temperature and becomes more pronounced with increasing chemical potential. Close to and across the CEP we find a strong increase of the entropy density with temperature and the slope becomes maximal at T = T CEP . The incorrect high-temperature behavior persists and becomes nonmonotonic.
The pressure follows again via integration according to Eq. (11) while the number density is computed as described in Ref. [26]. As seen in the upper right diagram, p gets larger with increasing chemical potential across the whole temperature range but the changes are less noticeable at low temperatures. For chemical potentials close to the CEP a kink starts to form at the corresponding transition temperature T c (µ). After that, the pressure rises stronger with a steeper slope as a function of T ; most pronounced and noticeable directly at the CEP. The pressure is, however, a smooth function of temperature for all µ up to µ CEP . These results are consistent with FRG results from the (Polyakov-loop enhanced) quark-meson model; see, e.g., Ref. [31].
The effect of a nonzero and increasing chemical potential is most prominent in the energy density (bottom left diagram) due to the additional number-density term µn(T, µ) from the Legendre transform; see Eq. (3). Its steep rise close to CEP indicates a rapid increase of degrees of freedom from hadrons to quarks and gluons. This behavior carries over to the interaction measure (bot-tom right diagram) which reacts strongly to chemical potential, too. It is shape-consistent with lattice results at small chemical potential and experiences a strong increase from intermediate chemical potentials onwards to the CEP. There, at µ = µ CEP , the slope becomes infinite at the corresponding critical temperature T CEP . The peak-like structure of I/ T 4 close to and at the CEP with a large magnitude indicates that nonperturbative effects are manifest in this region of the phase diagram.

C. The first-order phase boundary
With the pressure as a function of temperature and chemical potential at hand, the next obvious step is in principle the determination of the first-order phase boundary which lies between the spinodals (shaded area in Fig. 5). Unfortunately, as discussed in the following, there we hit a limitation of the method described in Sec. II.
In order to locate the first-order phase boundary, one considers the pressure difference between the chirally broken Nambu (N) and the chirally symmetric Wigner (W) phase for a fixed T < T CEP as a function of the chemical potential. Clearly, B(T, µ) is only defined up to the chemical potential µ N c = µ N c (T ) above which the Nambu solution does not exist anymore and only the Wigner solution can be found. The physically realized phase maximizes the pressure: B(T, µ) > 0 indicates that the Nambu phase is more stable than the Wigner phase and vice versa. Therefore, B(T, µ) = 0 defines the phase boundary. By finding the root of B(T, µ) with respect to µ for various T ∈ [0, T CEP ) one can draw the first-order phase boundary in the phase diagram.
In our approach, the pressure difference is explicitly given by with the functions and The first part (24) is evaluated at µ = 0 only, i.e., it does not depend on chemical potential. Moreover, ∆ µ=0 (T ) is positive for all T ∈ [0, T CEP ) since this region of the phase diagram is well within the hadronic phase where the Nambu solution is realized. Next, we find both in the NJL model as well as in our DSE setup that the number density as a function of µ (at fixed T < T CEP ) of the Wigner phase is generally larger than the number density of the Nambu phase: n W (T, µ) − n N (T, µ) > 0 and therefore ∆ µ =0 (T, µ) > 0 for all µ ∈ [0, µ N c ]. Furthermore, ∆ µ =0 (T, µ) is a monotonically increasing function of the chemical potential. It follows that there exists µ 1st c ∈ [0, µ N c ]-the location of the first-order phase boundary-where ∆ µ =0 (T, µ 1st c ) = ∆ µ=0 (T ) and consequently B(T, µ 1st c ) = 0. From the above discussion it becomes apparent that the correct location of the phase boundary depends crucially on the value of ∆ µ=0 . In particular, we need the pressure difference p N (T 0 , 0) − p W (T 0 , 0) as input while all other quantities in Eqs. (24) and (25) are computed from DSEs. Alas, lattice QCD provides to our knowledge only the physical Nambu pressure p N (T 0 , 0) but not the unphysical Wigner pressure p W (T 0 , 0) at a reference temperature T 0 . Thus, we are not able to obtain a reliable location of the phase boundary in the first-order region of the phase diagram. However, we would like to note that this is not an inherent flaw of our method but more that we hinge on the availability of an external input parameter.

VI. SUMMARY AND CONCLUSIONS
In this work we have studied the thermodynamics of strong-interaction matter within the DSE framework. We proposed a method to compute the entropy density solely from the quark condensate; a subsequent integration yields the pressure. Key feature of the method is that no approximation is used during its derivation and only the quark condensate is needed as input. This is particularly useful for DSEs where accessing the thermodynamic potential is extremely difficult and limited to truncations of the rainbow-ladder type. Even then, the proper removal of the quartic divergence contained in the potential is a nontrivial task. The proposed method provides a truncation-independent and straightforward way to compute thermodynamic quantities as soon as one gets hold of the quark condensate. The results, however, are truncation dependent since the quantitative behavior of the quark condensate depends on the chosen truncation.
That the method works effectively and yields reliable results was shown successfully using a NJL model. Then, we used condensate data obtained from a coupled set of Dyson-Schwinger equations for the quark and gluon propagators of (2 + 1)-flavor QCD within a truncation scheme used and discussed previously [26,37] to obtain the pressure, entropy density, energy density, and interaction measure from zero chemical potential up to the CEP as functions of temperature. These thermodynamic results are to our knowledge the first ones obtained from DSEs with a beyond-rainbow-ladder truncation-emphasizing the usefulness of the method presented in Sec. II.
At vanishing chemical potential, we find that our results for the pressure, entropy density, and energy density are in very good agreement with lattice QCD for temperatures below and around the pseudocritical chiral transition tem-perature. However, at high temperatures we observe an unphysical decrease of the entropy density with temperature. This erroneous behavior is rooted in the vertex ansatz of our DSE setup. Results at nonzero chemical potential show the expected behavior if one increases the chemical potential and are in qualitative agreement with results obtained within the FRG applied to the (Polyakovloop enhanced) quark-meson model. Unfortunately, we are not able to determine the firstorder phase boundary below the CEP. For that we need the pressure difference between the Nambu and Wigner phase at a reference temperature and at vanishing chemical potential as input. However, this quantity is as far as we know not provided by lattice-QCD calculations.
Finally, the results obtained in this work together with the ones on quark and baryon number fluctuations [26] made clear that an elaborate dressed quark-gluon vertex is needed for proper thermodynamics at high temperatures and/or densities. Especially terms in the vertex which react strongly to the (partial) restoration of chiral symmetry at high temperatures and/or chemical potential are of crucial importance. This extension of our current setup could be guided, e.g., by vacuum results for the dressed quark-gluon vertex or by an explicit calculation of (parts of) the vertex in medium [55,56]. A different yet complementary approach which takes more vertex structures into account is the FRG-assisted difference-DSE method proposed in Refs. [29,30] where nonzero temperature and chemical potential is treated as a fluctuation around the vacuum.