Finite size effects in strongly interacting matter at zero chemical potential from Polyakov loop Nambu-Jona-Lasinio model in the light of lattice data

We study finite volume effects within the Polyakov loop Nambu-Jona-Lasinio model for two light and one heavy quarks at vanishing baryon chemical potential and finite temperatures. We include three different Polyakov loop potentials and ensure that the predictions of our effective model in bulk are compatible with lattice QCD results. Finite size effects are taken into account by means of the Multiple Reflection Expansion formalism. We analyze several thermodynamic quantities including the interaction measure, the speed of sound, the surface tension, and the curvature energy and find that they are sensitive to finite volume effects, specially for systems with radii below $\sim 10$ fm and temperatures around the crossover one. For all sizes, the system undergoes a smooth crossover. The chiral critical temperature decreases by around $5 \%$ and the deconfinement temperature by less than a $2 \%$ when the radius goes from infinity to 3 fm. Thus, as the drop's size decreases, both temperatures become closer. The surface tension is dominated by the contribution of strange quarks and the curvature energy by $u$ and $d$ quarks. At large temperatures both quantities grow proportionally to $T^{3/2}$ but saturate to a constant value at low $T$.


Introduction
Understanding the hadron-quark phase transition is still a challenge from both the theoretical and experimental points of view. The framework for describing it is provided by Quantum Chromodynamics (QCD), which is the fundamental theory of strong interactions. However, the nonperturbative character of QCD at low energies makes extremely difficult to solve it in the regime of intermediate temperatures and chemical potentials, although lattice methods had a huge progress in the last years [1][2][3][4][5]. In this context, effective models such as the Nambu-Jona-Lasinio (NJL) model [6][7][8][9] are very useful because they can address many aspects of the QCD phase diagram without computational shortcomings at finite chemical potentials. The NJL model has many similarities with the full QCD theory but does not take into account the property of confinement, since quarks interact each other via pointlike interactions without exchanged gluons. Thus, in order to obtain a more realistic description, taking into account the quark confinement at low energies, the Polyakov loop was introduced in the NJL model [10], leading to the so called Polyakov loop NJL (PNJL) model (see also [11]). From this widely studied effective QCD model, many properties of strongly interacting matter can be obtained, such as its phase diagram [12][13][14].
On the other hand, a comprehension of finite size properties is very important for situations where the deconfinement transition occurs over a finite volume as in relativistic heavy ion collisions and neutron stars. The strongly interacting matter formed in a heavy-ion collision is finite in volume, and its size depends on the size of the colliding nuclei, the collision center of mass energy, and the centrality of the collision. In neut-ron stars, a deconfinement transition to quark matter is possible and a hybrid star or a strange quark star can be formed. The conversion of the star is expected to start with the nucleation of small quark matter drops [15][16][17][18] which subsequently grow at the expenses of the gravitational energy extracted from the contraction of the object and/or through a strongly exothermic combustion process. Quark matter droplets with a variety of geometrical forms can also arise within the mixed hadronquark phase that is expected to form inside hybrid stars if global charge neutrality is allowed [19]. Also, the most external layers of a strange star may fragment into a charge-separated mixture, involving positivelycharged strange droplets (strangelets) immersed in a negatively charged sea of electrons, forming a crystalline solid crust [20].
In the past years, many theoretical studies of finitevolume effects have been performed based on the NJL model [21][22][23]. However, studies within the PNJL model are more recent [24][25][26]. To incorporate finite-size effects different procedures have been employed, such as Monte Carlo simulations [24], a renormalization group approach [27], and the implementation of a low momentum cutoff Λ on the integration of the thermodynamic potential density of the PNJL model [28].
In the present work we use a different approach for the inclusion of finite size effects, known as Multiple Reflection Expansion (MRE) formalism [29]. First, different thermodynamic quantities calculated within our effective model in the bulk (including three different Polyakov loop potentials) are compared to the corresponding lattice QCD results. This is necessary as a starting point to check the validity of our model. Then, we study the relevance of finite size effects on many properties of strongly interacting matter and analyze how they deviate from the bulk case.
A comparison with lattice QCD is always important to calibrate effective models, that can be later extrapolated to a higher density regime. For example, the effective model can be used to explore finite-volume effects in a regime where they are known to be essential, such as in relativistic heavy ion collisions. Additionally, some results could be of interest for the analysis of the cosmological quark-hadron transition, which occurred in the early Universe about 10 µs after the Big Bang, when a hot unconfined quark-gluon plasma was converted, as the Universe expanded and cooled, into a confined hadronic phase.
The paper is organized as follows. In Sec. 2 we review the PNJL model in bulk for different Polyakov loop potentials and in Sec. 3 we introduce finite size effects through the MRE formalism. Our results are presented in Sec. 4 where we analyze the behavior of several thermodynamic quantities such as the chiral critical temperature, the deconfinement temperature, the constituent masses, the interaction measure, the pressure, the energy density, the entropy density, the speed of sound, the surface tension, and the curvature energy for different system sizes. Finally, we present our conclusions in Sec. 5.

The PNJL model in the bulk
The Lagrangian of the Polyakov loop extended SU (3) f NJL model including the six-quark 't Hooft interaction reads where q = (u, d, s) represents the three flavor quark field with three colors andm = diag(m u , m d , m s ) stands for the current quark mass matrix. We assume the SU (2) V isospin symmetry limit in which m u = m d . The covariant derivative in the fermion kinetic term couples a temporal background gauge field, the Polyakov loop, to the quark fields through D µ = ∂ µ −iA µ with A µ = δ 0 µ A 0 in Polyakov gauge, and A 0 = −iA 4 . Here, we used the notation A µ = gA a µ λ a /2 with g the SU (3) c gauge coupling. The λ a s stand for the Gell-Mann matrices with λ 0 = 2/3 1 in flavor space. The four-quark interaction coupling in the (pseudo)scalar channel is denoted by g S and the six-quark 't Hooft interaction coupling, induced by instantons, is labeled by g D . The latter one breaks the axial U A (1) symmetry. Finally, the above Lagrangian includes an effective potential U(l,l; T ) that accounts for gauge field self-interactions and is a function of the temperature T and the normalized colortraced Polyakov loop expectation value and its Hermitean conjugate, defined by where the Polyakov loop L is an N c × N c matrix in color space, as a function of A 4 . The explicit form of the Polyakov loop potential U(l,l; T ) will be discussed in Sec. 2.2.

The thermodynamic potential
Different thermodynamic properties of our model can be obtained from the thermodynamic potential in the mean-field approximation (MFA). The thermodynamic grand potential Ω(T, µ) of the PNJL model in the MFA has been largely considered in the literature, see e.g. [30,31]. Based on [30] we write the thermodynamic grand potential per unit volume as follows The first term is the condensation energy, that contains the contribution of the scalar four-quark interaction proportional to g S plus the six-quark 't Hooft interaction, proportional to g D . In the MFA this term depends on the three condensates ūu , d d and ss as follows The zero point energy is clearly divergent. Since the PNJL model is non-renormalizable, the zero-point energy contribution requires an ultraviolet cutoff Λ. The quark quasiparticle energies are denoted by i (p) = p 2 + M 2 i where the constituent quark masses M i for flavors i = u, d, s are: being m i the current quark masses.
The Ω quark term is ultraviolet finite and hence no momentum cutoff is imposed on it. It contains the coupling between the chiral condensates and the Polyakov loop L, and reads [30]: As shown in [30], taking an average of the 3 × 3 determinant we obtain: where The fourth contribution in Eq. (3) is a constant Ω vac ≡ −P vac , which is usually introduced in order to obtain a vanishing pressure at vanishing temperature and chemical potential. We will discuss the procedure for fixing P vac and its effect on the thermodynamic quantities in the next section.
Finally, the term U(l,l; T ) in Eq. (3), represents the pure gluonic effective potential in terms of the Polyakov loop variables, which will be presented below in detail. Notice that the U(l,l; T ) potential and Ω quark are invariant under the simultaneous exchange of l ↔l together with −µ i ↔ +µ i . Let us remark that for three quark flavors the thermodynamic grand potential Ω(T, µ i ) generally depends on three independent quark chemical potentials µ i . As a consequence of the isospin symmetry, the light quark chemical potentials are also degenerated. In the present work, we consider quark matter to be symmetric and define a common chemical potential µ ≡ µ u = µ d = µ s . Moreover, since we want to compare our results in the bulk with lattice QCD results we will work at finite temperature and vanishing chemical potential.
In order to obtain the dependence of the order parameters on the temperature and the chemical potential, one has to solve the following set of coupled equations: These conditions are consequences from the fact that the thermodynamically consistent solutions correspond to the stationary points of Ω P N JL with respect to ūu , d d , ss , l andl.

Polyakov loop potentials
The choice of the effective Polyakov loop potential U is not unique. In general, it can be constructed from the center symmetry of the pure-gauge sector. The required parameters can be extracted from pure gauge lattice data at µ = 0 [32]. Among several possible choices, see e.g. [33], we will use the following effective Polyakov loop potentials: (i) Logarithmic potential : the logarithmic ansatz presented in [33] is: where a(T ) and b(T ) are defined by [34]: with a 0 = 3.51, a 1 = −2.47, a 2 = 15.2 and b 3 = −1.75. (ii) Polynomial potential : Another choice is [31]: where b 2 (T ) = a 0 + a 1 (T 0 /T ) + a 2 (T 0 /T ) 2 with a 0 = 6.76, a 1 = −1.95, a 2 = 2.625, a 3 = −7.44, b 3 = 0.75 and b 4 = 7.5. In the absence of dynamical quarks, in a pure gauge sector, one expects a deconfinement temperature T 0 = 270 MeV. Nevertheless, in [35] it has been shown that in the presence of two light dynamical quarks and a massive strange one, this temperature is rescaled to about 187 MeV, with an uncertainty of about 30 MeV. In fact, for N f = 2 + 1, T 0 = 187 MeV and T 0 = 190 MeV have been used in [33] and in [25] respectively. Here we use T 0 = 185 MeV. (iii) Fukushima potential : Finally, we will use the strongcoupling inspired version of the effective Polyakov potential with only two parameters a and b proposed by Fukushima [30] The first term (proportional to ll) reminds the nearest neighbor interaction in the effective action at strong coupling and its temperature-dependent coefficient controls the deconfinement phase transition temperature. The logarithmic term comes from the Haar measure of the group integration with respect to the SU(3) Polyakov loop matrix. The parameters a and b are independent of the temperature, the chemical potential and the number of quark flavors N f . The parameter a controls only the deconfinement transition temperature and can be determined by the condition that the first-order phase transition in pure gluodynamics takes place at T = 270 MeV, which results in a = 664 MeV. On the other hand, the parameter b can be used to control the relative value of the deconfinement and chiral restoration crossover temperatures. Since there is no established prescription for fixing b, we shall adopt here two different values. First, we consider b = (196.2 MeV) 3 as suggested in [30,33] leading to an almost simultaneous crossover for deconfinement and chiral restoration at a temperature of T 200 MeV (we call this case U F 1 ). The second choice is b = (115 MeV) 3 (we call this case U F 2 ) which gives lower critical temperatures as we will see below.

Parametrization
In order to fully specify the non-local model under consideration we fix the model parameters following Ref. [36]. For comparison with some recent results [37], we have considered the parameters in [9], m u = m d = 5.5 MeV, m s = 135.7 MeV, Λ = 631.4 MeV, g S · Λ 2 = 3.67 and g D · Λ 5 = −9.29.

Finite size effects within the MRE formalism
Now we are ready to introduce the effects of finite size in the thermodynamic potential. For doing so we consider the MRE formalism (see Refs. [29,[38][39][40] and references therein) which takes into account the modification in the density of states resulting when the system is restricted to a finite domain. For the case of a finite spherical droplet the density of states reads: where the surface contribution to the density of states is and the curvature contribution is given by Madsen's ansatz [38] f which takes into account the finite quark mass contribution. The MRE density of states for massive quarks is reduced compared with the bulk one, and for a range of small momenta becomes negative. This non-physical negative values are removed by introducing an infrared (IR) cutoff in momentum space [40]. Thus, we have to perform the following replacement in order to obtain the thermodynamic quantities The upper integration limit is either infinity or given by a cutoff Λ. The IR cut-off Λ i,IR is the largest solution of the equation ρ i,MRE (p, m i , R) = 0 with respect to the momentum p.
After the above replacement, the full thermodynamic potential Ω MRE for a finite size spherical droplet reads: Multiplying on both sides of the last equation by the volume of the quark matter drop, replacing the area S = 4πR 2 and the curvature C = 8πR for a spherical drop, and rearranging terms we arrive to the following form for Ω MRE where the pressure P , the surface tension α and the curvature energy density γ, are defined as in Ref. [41]: As we previously mentioned, the value of Λ IR is the largest root when solving ρ i,MRE (p, m i , R) = 0 with respect to the momentum p, i.e. Λ IR changes with m i and with the drop's radius R.
Finally, we will address some aspects of the present model that deserve a more detailed discussion: (i) In the present treatment finite-size effects enter the fermion loop integral only; i.e. these effects are not considered in the pure Yang-Mills sector. As a consequence, the Polyakov loop potential remains unchanged and feels volume effects only implicitly through the saddle point equations. A more detailed analysis is left for future work. (ii) The conventional procedure for fixing P vac is to impose that the grand thermodynamic potential Ω must vanish at zero temperature and vanishing chemical potential for matter in bulk. For the above quoted parametrization, this assumption leads to the value P vac = 5080 MeV fm −3 . Nevertheless, it has been emphasized in previous works [42][43][44] that this prescription is no more than an arbitrary way to uniquely determine the EOS of the NJL model without any further assumptions. A change in the value of P vac has no influence on the fittings of the vacuum values for the meson masses and decay constants and thus the standard prescription for P vac is not related to experimental values. In fact, different prescriptions for determining P vac have been adopted [43], including the alternative of taking it as a free parameter [44] as it is usually done within the MIT bag model for the bag constant. When studying finite size systems, the standard choice for P vac has an additional issue.
If Ω vanishes at T = µ = 0 for matter in bulk it will not do so for a finite size, due to the contribution of surface and curvature effects (as can be seen from Eq. (26)).
As in previous works [45,46], we will fix P vac in the standard way, i.e. setting Ω = 0 at T = µ = 0 for matter in bulk, and will use this value for any system's size. Nonetheless, it must be emphasized that most of the thermodynamic quantities of relevance here (such as the critical temperatures, the entropy density, the sound speed, the specific heat, the surface tension and the curvature energy) are independent of the choice P vac since they are related to derivatives of the grand thermodynamic potential Ω. The influence of the P vac choice on other thermodynamic quantities will be discussed below.

Results
In this section we present our numerical results for some thermodynamic properties of bulk and finite size quark matter systems. We will show the dependence of our results on the size of the system as well as for different choices of the Polyakov loop potential. We work at zero chemical potential to compare our numerical results for the bulk with those from lattice QCD for (2+1)-flavors    using the highly improved staggered quark action extrapolated to the continuum limit [47] (see also [48]). Then we describe our predictions for finite size systems.

Chiral and deconfinement transitions
Here we will focus on the order parameters for both chiral and deconfinement transitions showing that, as the temperature is increased at zero baryon chemical potential, the PNJL model presents a smooth crossover transition at T ∼ 150 − 200 MeV depending on the size. Our results for the bulk are compatible with lattice QCD ones for N f = 2 + 1, as shown in Ref. [47] where the authors find a critical temperature of 154 ± 9 MeV (see also [48]). The chiral condensate is an order parameter for the spontaneous breaking of chiral symmetry [33]. The corresponding crossover transition can be established, for instance, by looking the temperature slope from χ u T ≡ ∂ ūu /∂T and χ l T ≡ ∂l/∂T . The peak positions give the inflexion points at the chiral critical temperature T χ and the critical deconfinement temperature T d of the condensates and the Polyakov loop expectation value respectively. As discussed in [30], it is convenient to take the crossover temperature in the u−sector, because the crossover temperature in the s−sector is larger and would be far from the deconfinement transition. As also discussed in [30], the chiral and deconfinement transitions do not take place at the same temperature as long as we treat the chiral condensates and the Polyakov loop as independent variables. Anyway, the idea is exploring different parameters in the Polyakov loop potential to force as much as possible the proximity of both critical temperatures. Also, as shown in [30,49,50] the peak of χ l T occurs at a lower temperature than the one of χ u T , in coincidence with our results presented in Tables 1 − 4.
For the different choices of the Polyakov loop potential introduced in Sect. 2.2, we present in Tables 1, 2, 3 and 4 the critical temperatures T χ and T d for different radii of the system. We also display the temperature T * below which the drop's pressure P would become negative for the standard prescription of P vac . Below T * these results would be unphysical if that P vac is adopted. If another choice of P vac is used the curves in Fig. 1 would shift upwards or downwards and a better coincidence of the model curves with lattice results could be achieved for the bulk case. From Tables 1, 2, 3 and 4 we see that (except for the cases U L and U F 2 with R = 3 fm) T * is always below T χ . Since we are interested in the physics around the critical temperature (relevant for heavy ion collisions and the early universe) we will not show our results for T < T * . We emphasize that the values of T χ and T d are independent of the choice of P vac .
In Table 1, we show our results for the polynomial Polyakov loop potential. The chiral critical temperature T χ has a significant dependence on the system size; it varies from 186 MeV to 177 MeV as the radius shrinks from infinity to 3 fm. This effect is also apparent in the left panel of Fig. 2 where we see that the peaks of ∂ ūu /∂T move towards smaller temperatures as the radius reduces. As seen in the right panel of Fig.  2, ∂l/∂T is less sensitive to finite size effects. Thus, the critical deconfinement temperature T d varies over a narrower range than T χ , as can be verified in Table  1. This behavior could have been anticipated because l feels volume changes only indirectly through the gap equations, and the Polyakov loop potential does not depend explicitly on the size of the system. As a consequence, T χ and T d get closer to each other as the drop's size decreases. The gray band are the results for the equation of state in (2+1)-flavor QCD using the highly improved staggered quark action extrapolated to the continuum limit [47] (see also [48]).
The critical temperatures for the model with a logarithmic Polyakov loop potential can be seen in Table  2. In this case T χ varies from 192 MeV to 181 MeV as R decreases. Here the deconfinement temperatures are slightly smaller than in the previous case, and the chiral ones, larger. For R = 3 fm, we find that T d lies in the negative pressure region for the standard choice of P vac .
In the cases with U L and U P , the choice of T 0 affects both, the deconfinement and the chiral critical temperatures. Here we use T 0 = 185 MeV in agreement with the values used in [25,33,35] for N f = 2 + 1. For larger T 0 , T χ and T d approach each other but both values increase, spoiling the coincidence with lattice results. On the other hand, for smaller T 0 , T χ and T d are closer to lattice data but there is larger separation between them.
Finally, we show the critical temperatures for the Polyakov loop potential of Fukushima. Here we considered two different examples, as discussed in the previous section. In Table 3 we show the results for b = (196.2 MeV) 3 , and in Table 4 for b = (115 MeV) 3 . The first case gives higher T χ and T d but both temperatures are closer to each other. In the second case we obtain smaller critical temperatures (closer to lattice results for 2+1 flavors) but there is a larger separation between them.
Summing up, in all cases discussed above, we see as a common behavior that T χ decreases with the size of the system by around 5% when the radius goes from R = ∞ to R = 3 fm. We also note that T d varies by less than 2% in the same size interval. As a consequence, as the drop's size decreases, T χ becomes closer to T d . This behavior is in agreement with the results presented in [24,25].
In Fig. 3 we show the temperature variation of the constituent masses M u , M d and M s for different drop sizes and for all the Polyakov loop potentials presented in Sect. 2.2. In the chirally broken phase, we find that the constituent quark masses are somewhat smaller for drops with smaller radii. In this region, the volume dependence of the effective masses is stronger than in the chirally restored phase. Also, M u and M d show a steep slope around the crossover temperature while for M s the slope is smoother. As shown in in Tables 1−4, the chiral critical temperature T χ shifts to smaller values as the volume decreases. Such behavior is also apparent in Fig. 3 where we see that, for smaller systems, the constituent mass tends to the current value at lower temperatures. A similar behavior has been reported in [25].

Interaction measure
A thermodynamic quantity of special interest is the thermal expectation value of the trace of the energy momentum tensor: This quantity is known as trace anomaly or equivalently as interaction measure ∆(T ) ≡ (T ) − 3P (T ) since it is very sensitive to the non-perturbative effects in the quark-gluon plasma. Specifically, it measures the deviation from the equation of state of an ideal gas = 3P due to interactions and/or finite quark masses.
Here we focus on the quantity which allows a straightforward assessment of deviations from the Stefan-Boltzmann behavior. Within the present model, the energy density (T ) at zero chemical potential is given by where the entropy density is given by: The interaction measure is sensitive to the finite drop's volume, because the energy density has an explicit dependence on the surface tension and the curvature energy:  In addition, as apparent from Eqs. (28), (29) and (30), there is an additional dependence on finite size effects through the infrared cutoff Λ i,IR in the integrals for P , α and γ.
In Fig. 4 we show our results for the bulk and for finite size systems together with lattice QCD simulation data in the continuum limit [47]. In general, we observe that the predictions of our effective model in bulk are in qualitative agreement with lattice QCD results. The peak heights are somewhat larger that in lattice QCD; nonetheless, the peak positions are in good coincidence with lattice. As a global feature, common to all finite sizes models that include different Polyakov loop potentials U, the interaction measure presents a peak that moves towards decreasing temperatures as the radii decrease. Note that, even though the interaction measure is explicitly dependent on P vac , the temperatures at which the peaks take place are not affected by the P vac choice.
For the chirally broken phase, i.e. for temperatures below the one in the peak, the curves for the bulk case are in qualitative agreement with lattice data. Close to T * , for finite sizes, the curves have a local minimum and start to increase at lower temperatures due to the contribution of the surface tension and the curvature energy. Now let us concentrate on the peak of the curves. For R = ∞, the peak position of the curves with U L , U P and U F 2 are in better coincidence with lattice results. The one with U F 1 is shifted to higher temperatures. For finite systems the position of the peaks is shifted to lower temperatures in all cases. As a global feature, the peak heights with R = ∞ for all models are larger than lattice results. We get a better agreement for U F 2 Figure 4 We show ∆/T 4 as a function of temperature for different drop sizes and different Polyakov loop potentials. We also include lattice QCD simulations data from [47] (gray band).
which is a ∼ 25% higher than lattice. For finite systems we see that the height of the peaks increase as the radii decrease, and they shift to smaller temperatures.
For high enough temperatures, in the chirally restored phase, there is a reasonable agreement between the bulk models and lattice results, specially for the U L and U P Polyakov potentials. Results for the U F 1 and U F 2 potentials, are somewhat below the lattice data. For finite sizes our results are superposed with the corresponding bulk case.

Energy density and entropy density
In the bulk case, our results for the energy density and the entropy density are in qualitative agreement with lattice QCD results (see Fig. 5) and with Ref. [33], as can be seen from their Fig. 3. As previously mentioned for the results of Fig. 1, different choices of P vac would lead only to a vertical shift of the curves for the energy density but will not change the temperature of the inflexion points. One could take advantage from this feature and introduce a different procedure to fix P vac in such a way that our predictions for the bulk case are as close as possible to lattice data. Since our focus here is not centered on the equation of state we shall not explore such strategy in the present work.
For high enough temperatures, our curves for all thermodynamic quantities approach to the Stefan-Boltzmann limit. The Stefan-Boltzmann limit for the pressure is given by where N c and N f are the number of colors and flavors. The first term represents the gluonic contribution and the second, the quark's contribution. For N c = 3 and N f = 3 we have: which results in SB /T 4 = 5.26 + 10.36 = 15.62 and s SB /T 3 = 7 + 13.8 = 20.8. From Fig. 5 we see that, at high enough temperatures, models with the Fukushima potentials U F 1 and U F 2 tend to the Stefan-Boltzmann limit for quarks only (no gluons) while models with the potentials U P and U L tend to the Stefan-Boltzmann limit including quarks and gluons.
This behavior is already known from previous works [25,30,33]. In the case of the U P and U L potentials, both the unconfined transverse gluons as well as the Polyakov loop, which corresponds to longitudinal gluons, contribute to the thermodynamic quantities [30]. But, since the Polyakov-loop potentials are fitted to pure gauge lattice data, they thus reproduce the total pressure, energy density, and entropy density of both the longitudinal and the transverse gluons, overcounting the degrees of freedom in the chirally symmetric phase [30,33]. However, the potential ansatz by Fukushima excludes these transverse gluon contributions at high temperatures leading to the differences found in Fig.  5. Nonetheless, at temperatures around and below the transition temperature such differences tend to disappear.
It is worth to remark that in Fig. 5 there is a wide range of temperatures in which our results for U F 2 are in a quantitatively good agreement with lattice results.
For finite systems, we see that in all cases the curves converge to the bulk ones at high temperatures. Close to the transition region, the curves for different radii start Figure 5 Energy density and entropy density as functions of temperature. Lattice data (gray band) are taken from [47].
to separate each other as the temperature decreases. In coincidence with Ref. [25], we find that the smaller the radius the higher the temperature at the inflexion point. Nonetheless, in Ref. [25] the results for R = 5 fm and for the bulk case are coincident for all temperatures but in our case are not. In the chirally broken phase, the energy density and the entropy density change very little with the drop's size.

Specific heat and speed of sound
The specific heat at constant volume is given by and the corresponding results are summarized in Fig.  6. At low temperatures c V grows with T , then shows a peak at the transition temperature, and approaches the corresponding Stefan-Boltzmann limit for high enough T . For the bulk case our results are in agreement with Ref. [33] and with lattice data [47]. In fact, in the left panel of Fig. 6 we note that lattice data show a soft undulation around the critical temperature, whose position is close to the peaks considering U F 2 and U P . For U F 1 and U L the peaks are shifted to higher temperatures. In general, the best agreement with lattice data (up to the critical temperature and somewhat above it) is obtained with U F 2 . For finite size drops, we find that c V doesn't change significantly with the change in volume, except in the crossover region. In fact, we find that the height of the peaks decreases as the volume shrinks, in agreement with [25]. Also, the peak position shifts to smaller temperatures as the volume decreases, as in [25]. As for other thermodynamic quantities, we find that the specific heat for the models with U F 1 and U F 2 tend to the Stefan-Boltzmann limit for quarks while the models with U P and U L tend to the Stefan-Boltzmann limit for quarks and gluons, due to the differences in the contributing gluon degrees of freedom.
In Fig. 6 we show our results for the speed of sound [33] The behavior of c 2 s is associated directly with the role of interactions in the system. The strength of interactions can be quantified through the interaction measure ∆ calculated in Sec. 4.2. A comparison between ∆ presented in Fig. 4 and c 2 s depicted in Fig. 6, shows that these quantities are correlated. At large temperatures, as the value of ∆ goes to zero, the speed of sound tends to the ultarelativistic limit of an ideal gas, c s = 1/ √ 3. At lower temperatures, interactions become relevant and therefore ∆ grows and c 2 s decreases significantly. Except for U F 1 , all minima positions lie close to the lattice QCD one. In the chirally restored phase our results for all Polyakov potentials are in good agreement with lattice QCD data, except for the U F 1 case that approaches lattice at higher temperatures.
Contrary to previous findings [25], our results show that the speed of sound doesn't depend too much on the system's size. In fact, small variations are observed only in the transition region.

Surface tension and curvature energy
In Fig. 7 we show the total surface tension α TOT and the total curvature energy γ TOT for drops with different  sizes. α TOT = i α i and γ TOT = i γ i include the contribution of u, d and s quarks. We have checked that α s is more than 10 times larger than α u and α d , in qualitative accordance with results for cold quark matter at very high densities [45,46] that show that the total surface tension α TOT is largely dominated by the contribution of strange quarks. On the contrary, γ u and γ d are typically ∼ 1 − 2 times γ s and thus the behavior of γ TOT is controlled mainly by u and d quarks.
Both, α TOT and γ TOT show a significant variation with R at all temperatures, specially for small drops with radii below 10 fm. There is also a considerable dependence on the Polyakov loop potential.
At large temperatures α TOT and γ TOT are monotonically increasing functions of T . Moreover, for T 250 MeV the surface tension grows approximately as being C α ≈ 0.029−0.034 MeV −1/2 fm −2 , while the curvature energy grows as being C γ ≈ 0.030 − 0.035 MeV −1/2 fm −1 . At lower temperatures both α TOT and γ TOT have local minimums. In the case of the total curvature energy the minimum falls around the chiral critical temperature T χ of the u and d condensates, which evidences the fact that γ TOT is controlled mostly by up and down quarks and is sensitive to their chiral transition. On the other hand, the total surface tension is sensitive to the chiral transition of strange quarks and therefore its minimum falls at a larger temperature.
At temperatures below that of the minimum there is a narrow interval where α TOT and γ TOT are decreasing functions of T . For even smaller temperatures, α TOT and γ TOT tend to a constant value which is of the same order of the values obtained within the NJL model for cold quark matter (T = 0) at finite chemical potentials (µ = 0 − 450 MeV) [46]. In some cases such constant value is not shown in the figures because the pressure becomes negative for the standard choice of P vac .

Summary and conclusions
In this work we studied the thermodynamic properties of finite systems composed by quark matter containing two light and one heavy quark within the frame of the PNJL model. We have considered vanishing baryon chemical potential and finite temperatures. We compared our numerical results for the bulk case with those from lattice QCD simulations, and then we studied the finite size deviations from the bulk case. We included finite size effects through the Multiple Reflection Expansion formalism and explored the effect of using different Polyakov loop potentials. Finite size effects were incorporated in the fermion integrals but not in the Polyakov loop potentials. However, if the pure Yang-Mills theory were formulated with a finite radius, the deconfinement phase transition could be affected and presumably the first-order transition would turn into a smooth crossover for small enough radii. This is beyond the scope of the present work.
As the temperature is increased at zero baryon chemical potential, the order parameters for both chiral and deconfinement transitions indicate that the PNJL model presents a smooth crossover transition, in accordance with lattice QCD results. For different radii of the system and different choices of the Polyakov loop potential, we determined the chiral critical temperature T χ of the u and d condensates and the critical deconfinement temperature T d of the Polyakov loop expectation value (see Table 1). In general, T χ depends on the system's size, decreasing by around 5% when the radius goes from infinity to 3 fm, while T d varies by less than 2% in the same interval. Thus, as the drop's size decreases, T χ becomes closer to T d , in accordance with [24,25].
Then we focused on the interaction measure ∆(T ) ≡ (T ) − 3P (T ), which evaluates the deviation from an ideal gas behavior ( = 3P ) due to interactions and/or finite quark masses. ∆/T 4 goes to zero at low and large temperatures and presents a peak around the transition density. In the bulk case, our results for ∆/T 4 are in qualitative agreement with lattice QCD results. Moreover, for U F 2 we obtain a good quantitative agreement with lattice data up to temperatures around 250 MeV. For different Polyakov loop potentials U, we find that as the radii decrease the peak moves towards lower temperatures and its height increases. At temperatures below that of the peak the results show a stronger dependence on the system's size and on the choice of the Polyakov loop potential.
In the bulk case, our results for the energy density , the entropy density s and the specific heat c V are in qualitative agreement with previous calculations presented in Ref. [33] and with lattice QCD results [47]. At high temperatures, the curves for all system's radii converge to the bulk ones and approach to the Stefan-Boltzmann limit. However, models with the Polyakov loop potentials of Fukushima tend to the Stefan-Boltzmann limit for quarks only (without gluons) while models with the polynomial and logarithmic Polyakov loop potentials tend to the Stefan-Boltzmann limit including quarks and gluons, as already known from previous works [25,30,33]. In general, , s and c V don't change significantly with the change in volume, except for c V in the transition region and for at temperatures below the transition region.
At high temperatures the speed of sound tends to the ultarelativistic limit of an ideal gas, c s = 1/ √ 3 but at lower temperatures, interactions become relevant and c 2 s decreases significantly. Again, for the bulk case we find a qualitative agreement with lattice QCD results. Notwithstanding, contrary to previous findings [25], our results show that the speed of sound doesn't depend too much on the system's size. In fact, small variations are observed only in the transition region.
Two very relevant quantities for finite systems are the surface tension and the curvature energy which have been calculated for drops with different sizes. We find that α TOT is largely dominated by the contribution of strange quarks (in coincidence with previous results for cold quark matter at very high densities [45,46]), while γ TOT is controlled mainly by the behavior of u and d quarks. Both, α TOT and γ TOT change significantly with R at all temperatures, specially for small drops with radii below 10 fm. There is also a considerable dependence on the Polyakov loop potential. For T 250 MeV, α TOT and γ TOT grow proportionally to T 3/2 . At lower temperatures α TOT has a minimum related to the chiral transition of s quarks and γ TOT has a minimum associated with the u and d quarks chiral transition. For smaller temperatures, α TOT and γ TOT tend to constant values of the same order of the ones obtained for very dense cold quark matter [46].
In summary, our main conclusion is that several thermodynamic quantities are sensitive to finite size effects, particularly for temperatures around the crossover transition and for systems with radii below ∼ 10 fm. These results can be potentially relevant for the study of the QCD transition at the early Universe [51,52] and should be extended to other regions of the QCD phase diagram, specially the region of high temperatures and moderate baryon chemical potentials where heavy ion collisions take place.