Sensitivity of finite size effects to the boundary conditions and the vacuum term

Finite volume effects are studied both with low-momentum cutoff and with momentum discretization in the framework of an (axial)vector meson extended quark-meson model with Polyakov-loop variables. In the momentum cutoff scenario, the critical endpoint (CEP) moves to lower temperatures and larger quark chemical potentials as the characteristic system size is reduced, however, the treatment of the vacuum term significantly affects its trajectory. The size dependence of the baryon fluctuations is also studied by the kurtosis and the skewness, both of which show moderate dependence on temperature and some dependence on quark chemical potential. The order of the phase transition is also studied near the chiral limit at finite system size and found to be second order only at vanishing explicit breaking. The implementation of the finite size effect with momentum discretization is more complicated and shows peculiar behavior due to the different modes dropping below the Fermi surface and strong dependence on the type of the boundary condition chosen. We found that both the different boundary conditions and the treatment of the vacuum term cause significant changes in the trajectory of the CEP as the characteristic system size is changed.


I. INTRODUCTION
The phase diagram of the strong interaction and, in particular, the existence and location of the critical endpoint (CEP) has been the subject of intense study both in theory and in experiment.However, heavy-ion collisions always suffer from the effect of finite system size, in contrast to field theoretical calculations in the thermodynamic limit.It is expected that a sufficiently small volume can affect the phase diagram and the CEP.
To study the effects of finite size in theoretical models, it is common to consider the constraints in momentum space imposed by the finite spatial extension.Such constraints can be a discretization or a simple lowmomentum cutoff using the phenomenological observation that the lowest modes are the most relevant for the phase transition.In the former case, the momentum integrals are changed to a summation over the modes determined by the boundary conditions, most commonly periodic boundary condition (PBC) and antiperiodic boundary condition (APBC).In the present work, we will focus on these momentum space constraints, but we note that there are other possible implementations of the finite size effects, e.g., the use of different distributions in the thermodynamics [1].
There have been several attempts in the literature to study the finite volume effects with both discretization and low-momentum cutoff schemes .The treat- * kovacs.gyozo@wigner.hument of the vacuum term and its finite volume correction varies among the different approaches employed in the analysis.In linear sigma model (LSM) studies, the vacuum and matter contributions are usually separated after performing the Matsubara sum.However, in the previous works with either discretization [2,3] or lowmomentum cutoff [4,5], the vacuum contribution is not present, and therefore its size dependence is not studied.
However, when the functional renormalization group approach is applied to quark-meson models, there is a renormalized fermionic vacuum integral, which is modified by decreasing size to study the physical quantities at T = 0 [6,7] or the phase transition at zero [8] and at finite quark chemical potential (µ q ) [9][10][11].
The Nambu-Jona-Lasinio (NJL) model is nonrenormalizable and therefore requires the specification of an UV regularization scheme to define the model.When examining the finite size effects, an additional complexity arises on the infrared side.To address this, the introduction of a low-momentum cutoff can be employed [12][13][14] or via implementing a discretization scheme [15,16,19].
The finite size effects are also studied with the Dyson-Schwinger (DS) equations with a low-momentum cutoff [25] and with discretization [26][27][28][29].In the DS approach, the modified integrals are more complicated compared to the other discussed models.Because of the different handling of the Matsubara summation, they have renormalized four-momentum integrals.
Combining the results of the different studies, one could conclude that, in general, the CEP is likely to move to lower T and higher µ q with decreasing size.However, not only the details of the criticality, e.g., the path and existence of the CEP for very small sizes, but also the size dependence of the phase boundary show different behavior in the different works with different approximations.By implementing several scenarios within a given model, one can investigate whether these differences are due to the approximations used, which will be the goal of the present work.
First, we use a low-momentum cutoff to implement the finite size effects and study the phase transition and thermodynamics at different system sizes.This scenario is also motivated by hadron resonance gas calculations [30,31], where it was found that the volume correction can be correctly reproduced by implementing a lowmomentum cutoff in the thermodynamic limit.We will also investigate a less studied problem, the size dependence of baryon fluctuations around the CEP as well, which has been briefly studied in the DS approach [28].
To make a comparison between the different implementations and boundary conditions, we also consider the momentum discretization with PBC and APBC, which is widely used, especially in the case of the NJL model and the DS approach.Furthermore, in contrast to other effective models, most previous studies on the quark-meson model have neglected the vacuum contribution to the fermion determinant and hence the size dependence of this term.Therefore we also study how the treatment of the vacuum term-i.e., finite or infinite size-affects the size dependence.
For this purpose, we use a vector and axial vector meson extended Polyakov quark-meson model (ePQM) with 2 + 1 flavors [32], which is based on [33] with a fermion one-loop contribution, including its-properly renormalized-vacuum and matter part as well.This model gives a good description of the meson phenomenology at T = 0, µ q = 0, and also of the chiral phase transition at finite temperatures, in good agreement with lattice results, while predicting a CEP and a first-order region for finite µ q .The ePQM model has been already used for many studies, including the investigation of inmedium vector and axial vector masses [34] and neutron star properties at T = 0 as the high-density part of a hybrid approach [35,36], as well as special limits of the phase diagram such as N c → ∞ [37].
The paper is organized as follows.In Sec.II, we briefly introduce the ePQM model.Section III is dedicated to the low-momentum cutoff scenario of the finite volume effects.After its implementation in the ePQM model, the modification of the phase diagram is discussed.In addition to the importance of the treatment of the vacuum term, we study the size dependence of thermodynamic quantities such as the equation of state (EOS) at T = 0 and the baryon fluctuations in the vicinity of the CEP.Furthermore, we show the size dependence of the transition and its order in the chiral limit in Appendix A. In Sec.IV, we turn to the discretization scenario.We discuss separately how the size of the vacuum part affects the phase transition and show the size dependence of the CEP with APBC and PBC.Finally, we conclude in Sec.V.

II. THE EXTENDED POLYAKOV QUARK-MESON MODEL
A vector and axial vector meson extended linear sigma model was developed in [33] with A global symmetry.The model was further extended with constituent quarks and Polyakov-loop variables in [32] (ePQM model), where the main focus was on the thermodynamics and the T − µ B phase diagram.The Lagrangian of the ePQM model is based on a scalar and a pseudoscalar nonet as well as left-and right-handed vector nonet fields, where T a = λ a /2 are the generators of the U (3) group, with λ 0 = 2 3 1 3 , and λ i (i = 1, . . ., 8) are the Gell-Mann matrices.Following [32], the Lagrangian can be written as where ψ = (q u , q d , q s ) T are the constituent quarks, and the quark mass matrix is defined as M = S + iγ 5 P .
The covariant derivatives and the field strength tensors appearing in (2) can be written with the help of the leftand right-handed vector fields and the electromagnetic field A µ e as The external fields corresponding to the scalar and vector fields are defined as Although in nature the isospin symmetry is broken, since its effect is small here we consider the isospin symmetric case (i.e., m u = m d ).As a usual procedure, we assume nonzero vacuum expectation values to scalar fields with zero quantum numbers, that is, ϕ N = ⟨σ N ⟩ and ϕ S = ⟨σ S ⟩, which are called the nonstrange and strange condensates.The corresponding scalar fields are subsequently shifted by their expectation values, σ N/S → σ N/S + ϕ N/S .To somewhat take into account the effects of confinement, it is usual to introduce Polyakov loops, which can be used to define the Polyakov-loop variables that signal the breaking of the center symmetry (for details, see [32]).Then the grand potential is calculated in a hybrid approximation, where the mesons are treated at tree level and the fermions at one-loop level, Ω(T, µ q ) = U (⟨M ⟩) + Ω (0) qq (T, µ q ) + U (⟨Φ⟩, ⟨ Φ⟩), (6) where the first term is the tree-level mesonic potential, and the second term is the fermionic contribution that includes the effects of the Polyakov loops and consists of a vacuum and a matter part, The third term is the Polyakov-loop potential, for which we use the improved U = U glue , being defined in [32].Within the current level of mean-field approximation, the exact choice of the Polyakov potential does not change our qualitative results on the phase transition and its finite size dependence as long as the thermodynamics is correctly reproduced.Alternative potentials are available [38][39][40] if one aims at analyzing the fluctuations of the Polyakov loops and matching the curvature masses to lattice QCD data [39].In the current model, it is necessary to tune the parameter T glue c = 182 MeV in the U glue in order to obtain a reasonable chiral pseudocritical temperature at vanishing chemical potential.Even this ad hoc procedure can be elevated by taking into account the dressing of the four-quark coupling with quark loops [41], which naturally resolves the problem and provides a natural explanation for the phenomenon of inverse magnetic catalysis.However, the study of fluctuations lies beyond the purpose of this work and we defer such a study for future endeavors.
The model parameters are determined with a χ 2minimization method using physical quantities like meson masses and decay widths listed in Table V of [32].The fitted parameters are the bare masses m 2 0 and m 2 1 , the meson-meson couplings g 1 , g 2 , λ 1 , λ 2 , h 1 , h 2 , and h 3 , the meson condensates ϕ N/S , the external field δ S , the U (1) A anomaly parameter c, and the fermion coupling g F .The set of parameters can also be found in [32].
After parametrization of the model, the thermodynamics is given by the field equations (FEs), The explicit form of these in the infinite volume limit is where the quark-antiquark condensate is given by while the matter part of the tadpole integral reads as Here are the modified Fermi-Dirac distribution functions with Solving the field equations above one can study thermodynamics at finite T and/or µ q and the phase structure of the strongly interacting matter.To determine the phase boundary between the chirally broken and restored phases, it is common to use the so-called subtracted condensate, which is also calculated on the lattice [42] and can be written in our case as This quantity can be generalized for finite µ q by choosing a proper reference point of the denominator.The phase boundary is given by the set of inflection points of ∆ µq,0 (T ) [or ∆ T0 (µ q ) in the µ q direction].

III. THE EPQM MODEL WITH LOW-MOMENTUM CUTOFF A. Implementation of low-momentum cutoff
The low-momentum cutoff λ = π/L can be included by simply introducing a Heavyside function for momentum integrals as The fermionic vacuum part with the low-momentum cutoff λ and the UV regularization Λ reads as, where the m f tree-level constituent quark masses are After the same renormalization as in the case of the infinite volume model, the fermionic vacuum contribution reads If λ is removed, the resulting expression agrees with the formula in [32].Using Eq. ( 18), the fermionic matter contribution with low-momentum cutoff is given by where g ± f (p) is defined in Eq. ( 16).Accordingly, the modified quark-antiquark condensate reads Now, different linear sizes are chosen, L ∈ (2 fm, . . .6 fm), which corresponds to different low-momentum cutoffs according to λ = π/L.This can be also expressed in MeV by using 1 fm −1 = 197.3MeV.Then the system of field equations in Eqs. ( 9)-( 12) is solved with all the integrals replaced according to Eqs. ( 18) and (23).For this we also need a parameter set, which is taken from [32].
The curvature meson masses for the scalars and pseudoscalars-which include one-loop fermion corrections-are also modified in case of low-momentum cutoff.Details can be found in Sec.B of [32]; here only the modifications are highlighted.The curvature meson masses are defined as where the tree-level part m 2 i,ab is unchanged and can be found in Table I of [32].In the matter part δm 2 ij , the formulas of Table I of [32] can be used by changing only the tadpole T f and bubble B f integrals 1 to their volume-dependent version T λ f in Eq. ( 24), and B λ f = −dT λ f /(dm 2 f ). 2 Finally, the fermionic vacuum contribution ∆m 2 ij can be written as , where m 2 f,i = ∂m 2 f /∂φ i and m 2 f,ij = ∂ 2 m 2 f /∂φ i ∂φ j are shorthand notations for the meson field derivatives of the constituent quark masses and can be found in Table II of [32].Whereas the derivatives of Ω

B. The phase diagram
If we change the system size L as described above, the solution of the field equation at T = µ q = 0 gives decreasing values for the meson condensates ϕ N and ϕ S for increasing L, which is shown in the top panel of Fig. 1 for ϕ N .It can be seen that the pseudocritical temperature T c at µ q = 0-defined as the temperature at the inflection point of ϕ N (T ) [or ∆(T )]-decreases with decreasing L. It is worth noting that on the top figure the curve belonging to L = 2 fm is already in the region where the pion and sigma masses have already become degenerate (see Fig. 2).It is known that finite size effects are similar to thermal effects if we consider different quantities as a function of 1/L.In Fig. 2 the tree-level masses of the u and d constituent quarks, and the curvature masses of sigma (or f 0 ) and pion are depicted as a function of 1/L in the vacuum (T = µ q = 0).The pion and sigma masses become degenerate around L ≈ 2 fm, which signals chiral symmetry restoration.The transition is of crossover type.We also note that the monotonic increase of the pion mass at small sizes is consistent with results obtained with the Lüscher formalism in matter part of the tadpole integral.Therefore, there is a −1 factor difference in the definition of T f compared to the cited work.chiral perturbation theory [43][44][45], although in our case the increase turns out to be too large.In the top panel of Fig. 3, the chiral phase boundary is shown for different characteristic sizes.As it was already pointed out, the pseudocritical temperature T c at µ q = 0 decreases with the decreasing L, behavior that has already been seen in the literature [9,12,28]; however, opposite behavior was also seen in the case of linear sigma models [2,4].This difference comes from the fact that, in our case, the fermionic vacuum contribution is taken into account (see also Sec.III C).The second-order CEP moves to smaller temperatures and slightly higher chemical potentials 3 and even disappears around L ≈ 2.5 fm.This shrinking of the first-order transition line and the absence of the criticality at small volumes is consistent with previous studies [2,4,9,12,28].When the first-order transition at T = 0 is present, the critical quark chemical potential µ q,c also decreases with the decreasing system size; thus the chirally broken region on the T − µ q plane becomes smaller with decreasing L. We note that a similar picture can also be obtained in the case of the chiral limit studied in Appendix A. Instead of a crossover, one finds a second-order phase transition at µ q = 0 ending in a tricritical point, which moves to lower temperatures and chemical potentials with the decreasing system size, similar to the CEP in the case of a nonvanishing explicit symmetry breaking.

C. Importance of the vacuum contribution
In previous works on finite volume effects within the framework of the linear sigma model, usually the fermionic vacuum contribution was not taken into account (see, e.g., [2,4,5]).To show the importance of the vacuum contribution, we also calculated the L dependence of the condensates and the phase diagram by keeping the vacuum part of the fermion integral at infinite size, i.e., setting λ = 0 in Eqs. ( 21) and (23).It is important to note that it is not the mere presence of the vacuum contribution that matters, but the way of its handling, i.e., its finite or infinite size. 4In the bottom panel of Fig. 1, the same curves are shown as in the top panel, but for the case where the vacuum is kept at infinite size.It can be seen that the T = 0 value of ϕ N remains naturally fixed, while the L dependence of the finite temperature integral pushes the pseudocritical temperature to higher temperatures, as can be seen in Fig. 4, the opposite of our previous case, but consistent with results in [2,4,5] Fig. 3-moves to smaller temperatures and higher chemical potentials again.However, with no finite size effects in the vacuum contribution, the path of the CEP has a qualitatively different curvature, and consequently, the critical quark chemical potential µ q,c at T = 0 increases with the decreasing L. We note that in this scenario the CEP does not disappear until L ≈ 0.5 fm. 5enerally, it can be said that size dependence of the vacuum pushes the system toward the chirally symmetric phase, which is obvious already from the transition in Fig. 2. The broken phase is squeezed down for decreasing size and even disappears completely around L = 2 fm, if the full fermion one-loop contribution is size dependent.Contrarily, for the vacuum part being fixed at L = ∞ the whole phase boundary moves further from the origin in the T − µ q plane, i.e., the chirally broken phase even extends.Consequently, at very small sizes the two ap- proaches give qualitatively different results regarding the existence of a chirally broken phase.

D. Thermodynamics and baryon fluctuations
The pressure is defined as while other thermodynamic quantities like the entropy density s = ∂p/∂T , the quark number density ρ q = ∂p/∂µ q , and the energy density ϵ = −p + T s + µ q ρ q can be derived from it.The equation of state, i.e., the pressure as a function of energy density p(ϵ) is an important quantity, which is shown at T = 0 in Fig. 5.It can be seen that the handling of the vacuum contribution affects the EOS.The difference partially comes from the increasing width of the unstable region between the spinodals in µ q at T = 0 and hence from the increasing ϵ at the transition point for decreasing L in the case of an infinite size vacuum, which pushes downward the p(ϵ) curve for small sizes.
To study the baryon number fluctuations, the generalized susceptibilities of the baryon number can be defined, while the so-called cumulants are related to these susceptibilities as Since the cumulants are explicitly volume dependent usually their ratios are used, which are equal to the corresponding susceptibility ratios.The most commonly used quantities are the skewness Sσ = χ 3 /χ 2 and the (excess) kurtosis κσ 2 = χ 4 /χ 2 , where σ is the variance, for which σ = χ 2 holds.Although explicit volume dependence is canceled in these ratios, they can still have implicit dependence as shown in [10,46], which was already investigated in effective approaches at vanishing chemical potential in [5,13] and with functional methods around the critical endpoint in [29].The cumulants and their ratios are proportional to higher powers of the correlation length, and therefore they could play an important role in the identification of the CEP in experiments [47][48][49].
Baryon fluctuations for infinite volume were already studied in the ePQM model in [50].When system size is decreased, the kurtosis at µ q = 0 smoothens as shown in Fig. 6, which is consistent with previous results of other effective model calculations [5,13].This behavior is in FIG. 6.The temperature dependence of the kurtosis κσ 2 for different system sizes L at µq = 0. line with the behavior of the ϕ N condensate that also smoothens (see Fig. 1).The temperature dependence of the skewness and the kurtosis is also calculated at µ q = µ CEP q , which can be seen in Fig. 7. Similarly, the quark chemical potential dependence of the same quantities at .The insets show the unscaled T dependence T = T CEP are shown in Fig. 8.It should be noted that the arsinh function is used to squeeze the function in the vertical direction for better visibility and to be directly comparable to results in the literature [29].It can be seen in Fig. 7 that there is a slight increase in the kurtosis for temperatures away from T CEP , which becomes larger with decreasing volume.It should be noted that this increase is somewhat artificial in the sense that the T CEP values used for the rescale change substantially with L. For instance, for L = 3 fm T CEP ≈ 10 MeV, which is ≈ T CEP ∞ /5.This can be also seen in the inset, where a slight decrease can be observed for smaller and smaller L. On the other hand, in Fig. 8, in the case of µ q dependence, both the rescaled curves and the inset show a decrease away from the transition, which is due to the fact that the change in µ q,CEP is moderate with decreasing L.
In conclusion, kurtosis shows a slight change for finite volumes but this effect is most probably originated from the CEP movement toward the chemical potential axis.The increase of κσ 2 is not remarkable until the CEP comes very close to the chemical potential axis.It is important to note that our result does not contradict what  ) was found in a Dyson-Schwinger approach in [29]. 6There the authors concluded that the kurtosis can be size independent around the CEP, since only a small deviation was found and it had not followed a clear trend in the system size.We also note that their CEP is located at a lower T CEP /µ q,CEP and significantly further away from the chemical potential axis.Note that the temperature dependence of the kurtosis near the critical endpoint has also been studied in a Polyakov-Nambu-Jona-Lasignio (PNJL) model with momentum discretization [17].However, due to the unusual phase transition and the multiple CEPs arising in finite volumes when the momentum is discretized (see Sec. IV C below), these results cannot be directly compared.

IV. THE EPQM MODEL WITH MOMENTUM DISCRETIZATION
Fourier transformation naturally implies discretization in momentum space if the system has finite spatial size.The modes that have to be taken into account are determined by the boundary condition imposed on the surface of the system's volume, which therefore also set the lowest momentum mode.In most cases, the finite system is assumed to be a cube with side length L and with periodic or antiperiodic boundary condition.In these cases the momentum integral is substituted with a sum that runs over the modes p i = 2n i π/L = n i ∆p or p i = (2n i + 1)π/L = (n i + 1/2)∆p, respectively, with n i ∈ Z, i ∈ (x, y, z), while ∆p ≡ 2π/L is the size of the momentum grid.

A. Implementation of momentum discretization
In the present work, we implement the periodic boundary condition, the antiperiodic boundary condition, and the periodic boundary condition without the zero mode (PBC-0).It should be noted that other prescriptions are also possible for both the boundary condition (e.g., Dirichlet) and shape of the system (spherical or spheroid, which would also require different BCs) that would give a better description of a real fireball produced in heavy-ion collisions.However, our purpose is to investigate the differences between the results of the different approaches used in the literature, for which the most commonly used PBC and APBC are the best choices.
Summation over the momentum modes can be numerically very expensive for large volumes.To reduce the computational costs, following [28,51], the summation is rearranged as follows: for any integral kernel K(p).Here the original grid points of the cube are rearranged in spheres and p j 's are the radii of larger and larger spherical shells in the momentum space as j increases, while m is indexing the grid points on a given shell.This is illustrated in two dimensions for APBC in Fig. 9, where it can be seen that, for example, for j = 2 one has p 2 = 1/4 + 9/4∆p and a multiplicity of 8 from the trivial sum over m.It can also be seen in the case of the two-dimensional sketch already, that the distance δp j = p j − p j−1 , j ≥ 1 between two consecutive spheres shows a decreasing trend.Moreover, if the momentum is large enough, an approximate δp j ∝ 1/p j behavior can be found, which allows an additional approximation.For integration in spherical coordinates, the integrand goes to zero in the matter part and diverges as p 3 in the vacuum part of the fermion integral for p → ∞.Therefore, in the matter part, the so-called UV improvement [27,28] can be used.In the UV improvement, an intermediate momentum scale λ Σ is introduced and the summation is changed to integration for momenta above this scale.This means that for p < λ Σ there is a summation over the discrete modesordered now in spheres-while for p ≥ λ Σ there is an integration with a low-momentum cutoff, similar to our previous scenario.It is worth clarifying the fitting of the boundary of integration to the boundary of summation.Since summation goes over for cubic cells of size ∆p 3 and the integration starts from a sphere, there are some parts of cubes around the boundary that penetrate into the integration region and some other parts that are missingboth for PBC and APBC (see, e.g., Fig. 9 in case of 2D).Thus, there is no common boundary between the summation and the integration.However, for the matter part, the error made by this mismatch is suppressed by the δp j ∝ 1/p j factor and the convergent behavior of the integrand for large momentum.λ Σ has to be set for each volume separately, by finding a region where the results are already insensitive to the change of the value of this cutoff.Therefore, the continuum integration can be changed for finite volumes to where j max = max{j|p j ≤ λ Σ }.Consequently, the fermionic matter contribution can be written as Unfortunately, in the vacuum part, the integrand is divergent for p → ∞ strong enough to overcome the δp j ∝ 1/p j behavior.Consequently, the error made by changing from a cubic-based summation to a spherical integration at λ Σ increases with λ Σ .Note that this error purely comes from the mismatch of the boundaries of the summation and the integration.This can be seen for APBC as follows: let us keep the original summation starting from some large p 0 value up to p 0 + n∆p, i.e., p 0 ≤ p i ≤ p 0 + n∆p, i ∈ (x, y, z).This corresponds to an integration in Descartes coordinates for the interval p 0 − ∆p/2 < p i < p 0 + (n + 1/2)∆p, i ∈ (x, y, z).The summation and the integration are for exactly the same volume, thus the summation is nothing else but a Riemann sum of the integral.For arbitrary n and p 0 ≫ ∆p it can be shown that with increasing p 0 the error made by the summation with midpoint rule compared to the integration decreases, as can be seen in Fig. 10.In this case, K(p) = E(p) and E(p) ∝ p for large p.Consequently, considering a well-defined cubic-based boundary between the summation and the integration at a momentum p 0 , the error, made by changing from summation to integration, can be made arbitrarily small by setting p 0 large enough.Therefore, integration can be used in the large momentum part of the fermionic vacuum contribution too.In addition to making the calculation numerically less expensive, this also enables the renormalization of the vacuum part, similar to the infinite volume case.The simplest way to use the UV improvement is to change from summation to integration at λ Σ in each direction; that is, to sum up to a cube of side length 2λ Σ .Note that the value of λ Σ has to be chosen to be half-integer or an integer multiple of the cell size, i.e., λ Σ = (n + 1/2)∆p for PBC and λ Σ = (n + 1)∆p for APBC.To perform the integration that starts from the cube (of side length 2λ Σ ), the integration is divided into two parts, the first one from the cube to a sphere of radius √ 3λ Σ , i.e., the smallest sphere that includes the cube, and the second one which is a usual spherical integral that starts from p ≡ |p| = √ 3λ Σ .Since the integrand has rotational invariance, in the first part of the integral in spherical coordinates the solid angle integral can be calculated explicitly giving a p-dependent solid angle as with and where there is a change in the expression for the solid angle at p = √ 2λ Σ , when the sphere touches every middle point of the cube edges.The formula above is derived in Appendix B. Thus, the fermionic vacuum contribution in Eq. ( 7) can be written as three terms, a sum up to p i = ±λ Σ , an integration from the cube to the sphere with radius √ 3λ Σ , and a renormalized vacuum integral of the form of Eq. ( 21) with lower cut-off λ We emphasize that the need for such a complicated adjustment of the intermediate scale for the UV improvement arises from the fact that the discretization does not respect the symmetry of the integration and the renormalization schemes, which has rotational invariance.To overcome this problem, it would be possible to study the finite volume effect with a judicious choice of scheme, so that it is consistent with the symmetry chosen in the renormalization scheme, e.g., with discretization in the radial direction similar to [52], where an O(3) scheme was used for the Matsubara sum in 2 + 1-dimensional spacetime.This would also be in line with the considerations to implement a rotation invariant, spherical shape of the finite volume system.
To illustrate the stability of the UV improvement with the cubic boundary and the cubic-spherical integration, we show the λ Σ dependence of the tree-level constituent quark masses-which are proportional to the meson condensates-at fixed L in Fig. 11.The curves start at 1, since for λ Σ = 0 the summation is absent and integration is performed in the whole momentum range, which is equivalent to the infinite size case.In the calculations, we used λ Σ /∆p > 50 to ensure the stability.
The field equations and the curvature meson masses should be modified according to the grand potential in a similar fashion to the low-momentum cutoff case discussed in Sec.III A. Every integral derived from the fermionic matter part has to be changed to the UVimproved version given in Eq. (33).It is straightforward for Eqs. ( 9) and (10), while the matter part of the tadpole integral in Eq (14) should be modified as with T λΣ f being the tadpole with low-momentum cutoff [Eq.(24)].Moreover, the fermionic vacuum contribution to the quark-antiquark condensate [first part of Eq. ( 13)] can be derived from Eq. (38), where the last term is the renormalized contribution with a low-momentum cutoff λ = √ 3λ Σ defined in Eq. ( 23).The tree-level parts of the meson masses remain unchanged, while for the fermionic matter contribution, δm 2 ij in Table I of [32] can be used by replacing the tadpole and bubble integrals with T Σ f and B Σ f = −dT Σ f /(dm 2 f ), respectively.Moreover, the fermionic vac-uum contribution is given similarly as in Eq. ( 26), , where m 2 f,i and m 2 f,ij can be found in Table II of [32].

B. Vacuum contribution
Similar as in Sec.III C, we investigated the role of handling of the vacuum contribution, i.e., its infinite or finite size.As discussed in the previous section at T = µ q = 0 the finite size effect to the field equations comes from the discretization of the vacuum part of the chiral condensate ⟨qq⟩ Σ vac [Eq.(38)].Unfortunately, as the system size is decreased, around L ⪅ 4.5, the solution to the field equations ceases to exist.Contrary to the low-momentum cutoff scenario, here the magnitude of the vacuum contribution to the FEs increases with decreasing L, which subsequently causes the disappearance of the common solution to the two FEs [Eqs.(11), (12)] at T = 0. Note that this kind of instability for large meson condensates is a common feature of the quark-meson models at meanfield level with fermionic vacuum contribution.This is due to the fact that the vacuum term gives a leading contribution to the grand potential, which is ∝ −σ 3 log σ with the general order parameter σ, and hence the grand potential is not bounded from below as a function of σ.When the size decreases in the discretized scenario, the local minimum moves to larger values of the order parameter and collapses with the maximum, hence leaving no extremum of the grand potential.Nevertheless, we can still look at the effect of different boundary conditions as a function of 1/L if we fix m u (and consequently ϕ N ) and calculate ⟨qq⟩ Σ vac (L).This can be seen in Fig 12 .The most significant contribution to ⟨qq⟩ Σ vac is at p = 0, which is taken into account in the case of PBC and omitted in the case of APBC.Consequently, it would be expected that the absolute value of the chiral condensate and thus (via the field equations) ϕ N and ϕ S increases for PBC and decreases for APBC with decreasing L. In terms of chiral symmetry, it means that the APBC drives the system to chiral restoration while the PBC increases the spontaneous symmetry breaking.In NJL model studies with discretization of the regularized chiral condensate, such behavior was found in [15,16,19].A qualitatively similar result was obtained in a parity doublet model in [23].We note that in these cases the vacuum contribution is only regularized with special functions, and no renormalization is performed.However, in the ePQM model for L > 1 fm, the main difference between the integral at L = ∞ and the summation at finite L comes from higher modes, and both boundary conditions give an increasing absolute value for the chiral condensate, as shown in Fig. 12.Therefore, as L decreases from infinity and the finite size effect appears, the chirally broken phase extends to higher T (and µ q ) not only for PBC but also for APBC.
It is worth noting that if the vacuum size is kept infinite the solution exists, which will be investigated in the next section.On the other hand, it is still possible to take into account the finite size effect of the vacuum in a restricted way, as it was done in [2] for PBC, 7 if we consider only the lowest mode [i.e., ⃗ p = (0, 0, 0) for PBC and ⃗ p = (1/2, 1/2, 1/2) for APBC].Practically, this means that only the first term is kept from the sum in Eq. (40).In this way, the solution to the FEs will exist up to L ≈ 2.5 fm.
Finally, in Fig. 12, the PBC-0 case is also shown (see, e.g., in [28]).Here, the low-momentum cutoff caused by the exclusion of the zero mode gives the main modification.Consequently, the absolute value of the chiral condensate decreases as soon as the finite size effect appears around L = 10 fm, similar to the low-momentum cutoff scenario.this scenario, the solution to the FEs exists.

C. Finite temperature and chemical potential
As a consequence of the analysis presented in the previous section, from now on the fermionic vacuum contribution is either kept at infinite volume, or the discretization is taken into account only in the lowest mode, and we implement the summation only for the matter part.To get a better understanding of the size dependence at finite T and µ q in addition to the ePQM model, we also considered a simpler linear sigma model presented in [2] (called LSM A) and another LSM model with two different parametrizations from [53].In the latter case, the Lagrangian is similar to Eq. ( 2) without the (axial) vector meson sector and without fermionic vacuum contribution.The two sets of parameters used here are the one with m σ = 600 MeV (LSM B) and with m σ = 800 MeV (LSM C).The parameter sets for LSM B and C can be found in the 11th and 13th row [bottom panel, with U (1) A breaking] of Table II in [53].
Besides the absence of the fermion vacuum contribution in LSM A, LSM B, and LSM C, the location of the CEP is also very different.The ePQM model and LSM C both have a "low-lying" CEP with µ q /T = 5 − 6 at (µ q , T ) = (295, 53) and (328, 63) MeV, respectively, whereas LSM A and B have a "high-lying" CEP with µ q /T = 2 − 2.5 at (207, 100) and (220, 92) MeV, respectively. 8efore turning to the effects of finite size on the CEP in the different model,s we discuss the µ q = 0 and T ≈ 0 cases in the ePQM model.At µ q = 0 with the vacuum contribution of infinite size, the pseudocritical temperature increases for APBC (after a transient decrease from 1/L = 0.1 to 1/L = 0.25 fm −1 ), PBC and also for PBC-0 as shown in Fig. 13, labeled as T A pc , T P pc , and T P −0 pc , respectively.As can be seen in the case of PBC, the solu-0.9 is that it has no continuation at higher temperatures.This is illustrated in Fig. 14.Therefore, the crossover is suddenly replaced by a nonphysical first-order transition.By including the effect of the finite size, at least in the zero mode of the vacuum contribution, the zero mode of the matter part is compensated, and the above problem does not exist.On the other hand, the solvability of the model at T = 0, µ q = 0 is still limited to L > 2.5 fm, due to the increasing condensates (caused by the increasing vacuum contribution) discussed in the previous section.The resulting pseudocritical temperature, denoted as T P +v0 pc , is also shown in Fig. 13.For PBC, contrary to the low-momentum cutoff scenario, the finite size effects on the vacuum contribution (implemented only in the zero mode, but also if the full discretization is applied) enhance the chiral symmetry breaking, 9 hence the pseudocritical temperature T P +v0 pc increases with the decreasing system size.
In the T → 0 limit, the Fermi-Dirac distribution (as a function of µ q ) becomes a step function, which is one below the Fermi surface, i.e., E(p) < µ q , and zero above.When the momentum space is discretized, only the modes below the Fermi surface contribute.With increasing µ q , the lowest modes enter below the Fermi surface, causing jumps in the fermionic matter contribution in the FEs and thus leading to drops in the meson condensates. 10Moreover, the decrease of the condensates, and hence of the fermion masses in E(p), enhances the expansion of the Fermi surface (generally driven by the growth of µ q ), which gets closer to or even grows beyond the next modes.This leads to a staircaselike decrease of the condensates with multiple steps, hence multiple inflection points, which look like first-order transitions on of each other. 11The structure of ϕ N (µ q ) is shown in Fig. 15, where L = 6 fm and APBC is applied for the ePQM model at vanishing and nonzero temperatures.At T = 0 the horizontal steps are connected by quasisolutions (marked with dashed lines), where ∂Ω/∂ϕ N/S = 0 has a solution, but ∂Ω/∂ϕ S/N has a discontinuity.These quasisolutions change to real solutions for nonzero temperatures.Since the transition at T ≈ 0 is of first order new second-order points appear on the new quasi transition lines.A similar staircaselike phase transition has also been observed in PNJL model calculations [17], where it was named "quantized first-order phase transition."In these models, the gap equation, which is solved self-consistently to find the value of the chiral condensate at a given T and µ q , formally contains the same integral as the field equations in Eqs.(11) and ( 12) and thus suffers from the problem of discrete modes entering below the Fermi surface in the same way.For PBC, similar behavior to APBC can be seen when the unphysical first-order transition is not present.From this behavior for both boundary conditions it can be predicted that the low-T and high-µ q region of the phase diagram, where the CEP is located, is strongly influenced by the T ≈ 0 structure.The behavior of the CEP is considered separately for APBC and PBC in the following subsections.
We note that we are working in the grand canonical ensemble, where the chemical potential is an external 11 Instead of one phase transition from the chirally broken to the symmetric phase, in the present case there are smaller steps (marked generally by the lines of inflection points, but for a firstorder transition also by having multiple solutions at a given T and µq).Therefore, the chiral restoration happens in multiple steps, which we will call "quasitransitions" since they separate phases with only quantitatively different spontaneous symmetry breaking.
parameter and the Fermi surface is independent of L.
Meanwhile, the distance between the modes increases with the decreasing size, and hence some of them fall out for a given µ q .Alternatively, for finite size, one could first formulate the canonical picture with the discretized momentum space and then derive the Fermi surface from the chemical potential defined by the change of the free energy F with the particle number N of the discretized system µ q = dF/dN .To see whether such a treatment changes the results or not, further investigations are needed.

Size dependence of the CEP with APBC
The trajectories of the CEP can be seen in Fig. 16 as L is decreased in the case of APBC for the models discussed in the previous section.In each case, the CEP belonging to L = ∞ (marked with a black dot) starts to move to lower temperatures and higher quark chemical potentials.This can be seen with the solid line for the LSM A and LSM B models and with the dashed line for the LSM C and ePQM models.However, as discussed in the previous section, the ϕ N (µ q ) has multiple inflection points (Fig. 15).For the ePQM model and LSM C, the inflection point with higher µ q and lower ϕ N is connected to the CEP belonging to L = ∞ (dashed trajectory), but at some L value this part of the curve becomes a crossover.On the other hand, the second inflection point (with lower µ q and higher ϕ N ), which coexist in some region with the first one, separates the crossover and first-order regions for L ⪅ 4.5 fm; therefore it is reasonable to identify it as the CEP of the chiral phase transition (solid trajectory) below the given L. The second inflection point is induced by the first mode when entering below the Fermi surface.For sufficiently small L, this new CEP also turns to lower T and higher µ q (like LSM A and LSM B), while the original CEP is continued toward the µ q axis generated by the second lowest modes when entering below the Fermi surface.Similar nonmonotonic behavior was also found in [10] within a functional renormalization group approach.For LSM A and LSM B, the original CEPs start from a much lower µ q and higher T , so their trajectories continuously connect the CEP at L = ∞ to the CEP with the smallest possible L. The quasi-CEPs (second-order points at lower ϕ N values) of the transitions generated by the second, third, etc. lowest modes (when entering below the Fermi surface) are also present for LSM A and B, but their presence remains hidden.The 1/L dependence of µ CEP q and T CEP is also shown separately for each model in Fig. 17.It can be seen that for 1/L < 0. the curves are flat.For 0.2 fm −1 < 1/L < 0.5 fm −1 (2 fm < L < 5 fm) for the ePQM and LSM C models, the original CEP disappears, while a newly generated CEP emerges and becomes dominant.In this range the T CEP (µ CEP q ) is continuously decreasing (increasing) for LSM A and B. Finally, for 1/L > 0.5 fm −1 (L < 2 fm) both µ CEP q and T CEP show a very similar trend for each model.This also shows that for sufficiently small sizes the CEP corresponds to the second-order point belonging to the highest ϕ N (first step in Fig. 15) of the transition, which is generated by the first mode when entering below the Fermi surface.

Size dependence of the CEP with PBC
For PBC with a vacuum term of infinite size, the CEP moves to higher temperatures and lower quark chemical potentials with decreasing system size, leading to an expansion of the first-order region for each model, as shown (with solid lines) in Fig. 18.However, the appearance of the new solution of the field equations discussed in Sec.IV C turns the phase transition, either crossover or first order, into an unphysical first-order transition in each model at roughly the same size around L = 5.5 fm (at L = 5.46 for the ePQM), thus we cannot continue the calculation to smaller sizes.For LSM A and B models, the CEP reaches the T axis, leaving a first-order phase transition in the full phase boundary well before the nonphysical effect becomes dominant.For the LSM C and the ePQM models, there is still a crossover region where this unphysical transition occurs and the CEP is still very far from the T axis.This CEP is marked with a black cross in the figure.
As mentioned in Sec.IV C, by including the effect of the discretization in the vacuum part only in the zero mode, the nonphysical transition disappears.Furthermore, this zero-mode correction of the vacuum contribution 12 is large enough to even reverse the trend of the CEP trajectory as L is decreased.Therefore, it moves (for ePQM after some transient oscillations) to higher temperatures and lower chemical potentials with decreasing size, resulting in the reduction of the first-order region.The trajectory of the CEPs in this scenario is also shown in Fig. 18 with the dash-dotted lines.For LSM A, this approximation reproduces the results in [2] for PBC.We note that for low T and high µ q there are multiple solutions of the FEs [Eqs.(11) and (12)] due to the structure generated by the modes entering below the Fermi surface, which was discussed in the previous section (see Fig 15).However, in this case, the structure of ϕ N/S (T, µ q )-and hence the structure of the order parameter ∆(T, µ q )-is more complicated than in the case of APBC with infinite vacuum, and it is difficult to identify well-defined phase transitions for low T and high µ q .Therefore, when the T = 0 behavior begins to affect the T = T CEP region at sufficiently small sizes-below L =4-5 fm for LSM C and the ePQM models, and below L = 2 fm for the LSM A and B models-the CEP will not be tracked further with decreasing size.
For PBC-0, the influence of the missing zero mode of the fermionic thermal contribution-implying a cubic low-momentum cutoff-will be dominant.It overwhelms any effect of summation and pushes the CEP toward the chemical potential axis.However, the appearance of the new CEP arising from the T = 0 behavior influences the path of the critical endpoint also for this boundary condition.The location of the CEP for very small system sizes is determined only by the first mode-following the absent zero mode-entering below the Fermi surface.

V. CONCLUSION
In this paper we have investigated the effects of finite system size in the framework of the vector and axial vector meson extended Polyakov linear sigma model by using two types of implementations: a low-momentum cutoff of the momentum integrals and a discretization of the momentum space implied by the finite spatial extent of the system.The effect of these modifications on the phase diagram of strongly interacting matter, the role of the vacuum term and the different boundary conditions have been studied in detail.
First, the case of low-momentum cutoff was investigated.In particular, the effect of the size of the vacuum, either finite or left unchanged (infinite size), was studied.The latter leads to similar results as in previous studies where the fermionic vacuum corrections were not considered.It was found that the modification of the vacuum integrals strongly pushes the system toward the chirally symmetric phase.Therefore, the whole phase boundary moves to lower temperatures until the chirally broken phase vanishes around L ≈ 2 fm, where a crossoverlike transition can be observed in the vacuum physical quantities.However, if the vacuum is kept in its infinite form, the chirally broken phase even extends as the boundary moves to higher T and µ q .In both cases, the critical endpoint moves to lower temperatures and higher chemical potentials.For a vacuum term of finite size, it turns toward the chemical potential axis and even disappears around L ≈ 2.5 fm.On the other hand, for an infinite size vacuum, the CEP turns toward high µ q , and although the first-order line shrinks with decreasing size, the CEP does not reach the axis until at least L ≈ 0.5 fm.The volume dependence of the phase diagram and the order of the chiral phase transition were also studied in and near the chiral limit, i.e., the m u,d = 0, m s ̸ = 0 axis of the Columbia plot.It was found that the phase structure changes with decreasing size in a very similar way to the physical case, while a second-order phase transition is found only at vanishing explicit breaking and crossover elsewhere.
We also studied the modification of the baryon fluctuations at finite volume in the low-momentum cutoff scenario.At µ q = 0 we found a smoothening of the kurtosis as a function of temperature, which was expected.The kurtosis and the skewness were also calculated through the critical endpoint (µ q = µ CEP q or T = T CEP ) as a function of T /T CEP or µ q /µ CEP q .κ(T /T CEP ) µq=µ CEP q showed an increasing behavior on both sides, not very close to the CEP.However, we concluded that this trend might be a consequence of the rescaling by T CEP , since for small volumes the critical endpoint moves to lower T values.For κ(µ q /µ CEP q ) T =T CEP the rescaling by µ CEP q did not alter the results and a slight decrease was observed near the critical endpoint, which might be also connected to the significant shift of the CEP.Therefore, we claimed that the kurtosis probably does not possess a strong size dependence, which is consistent with what has been found in the literature using a DS approach.
As an alternative to the low-momentum cutoff, the discretization of the momentum space was also studied.On the one hand, it turned out that the implementation of this approach is more complicated, and on the other hand, in several cases, the solution of the field equations was lost.Moreover, the µ q dependence of the condensates showed a more complicated, staircaselike structure than in the low-momentum cutoff case, due to the different modes falling below the Fermi surface-which may be an artifact of our assumption due to the size-independent Fermi surface, but could also be realized in certain physical systems, e.g., in solid state physics.This behavior makes it more difficult to define a unique transition point and causes some unexpected behavior in the trajectory of the CEP as a function of the characteristic size L, and the resulting finite size effects are not conclusive.The choice of different boundary conditions, such as PBC, APBC, and PBC-0, also causes significant changes in the trajectory of the CEP.We have also compared the behavior of our model with other models in the literature.
Comparing the predictions for the size dependence of the CEP for the low-momentum cutoff with and without modification of the vacuum part (see Fig. 3) and for the discretization with APBC (Fig. 16) and PBC (Fig. 18), it can be seen that completely different behavior can be found.The role of the treatment of the vacuum term and the different boundary conditions is paramount in the exact trajectory of the CEP as L decreases.The only similarity seems to be that for small L the CEP moves toward the µ q axis, as expected.
Since the present results do not provide a consistent picture of the details of the finite size effects on the phase diagram, we need to find better approaches to this problem.For this purpose, in addition to improving the models used, we need more sophisticated methods to implement these effects, e.g., by considering more physically motivated shape and boundary conditions for the system.Furthermore, besides the momentum space constraints, one could consider the finite size dependence in direct space as well.However, this would lead to a nonhomogeneous potential and hence to nonhomogeneous meson condensates, which would require overcoming new obstacles. form where Z π is the pion wave function renormalization constant that is connected with the pion decay constant by f π = ϕ N /Z π .For h N → 0, and nonzero spontaneous symmetry breaking (i.e., ϕ N > 0), the pion mass has to vanish as expected in the chiral limit (Goldstone theorem).We note that to connect the parameter h N to the current quark mass, which does not appear naturally in a quark-meson model, a possible way is provided by the Gell-Mann-Oakes-Renner (GMOR) relation, m 2 π f 2 π = −2m 0 ⟨qq⟩, where m 0 is the current quark mass, while ⟨qq⟩ is the chiral condensate.According to Eq. (A1) the GMOR relation can be rewritten as h N ϕ N = −2m 0 ⟨qq⟩, which directly connects the product of the explicit and spontaneous symmetry breaking parameters in the linear sigma model to the parameters in a quark-based model, like the NJL.Although h N ϕ N gives the expected physical value for the product, to obtain m 0 (≡ m u ) the identification of the ⟨qq⟩ chiral condensate would be needed in the quark-meson model separately.However, in the ePQM model, the quark-antiquark condensate is included in the meson masses via the fermion one-loop correction.Therefore, the parametrization that is based on the meson masses fixes the value of the chiral condensate, which turns out to be too low in our case.
In Fig. 19, the chiral phase boundary is shown in the chiral limit for the N f = 2 (light) + 1 (heavy) case, i.e., for h N → 0, h S = h phys moves to lower temperature with the decreasing h N and the crossover becomes a second-order phase transition at h N = 0, which is consistent with the behavior found in previous studies [57][58][59].The CEP moves to lower quark chemical potentials and higher temperatures and ends up in a tricritical point, as it is expected.The universality class of the second-order phase transition is of the mean-field theory, which can be seen in Fig. 20, where the scaling of the order parameter ∆ (subtracted condensate) is presented in the chiral limit.In Fig. 21, the volume dependence of the phase boundary can be seen in the chiral limit.The solid and dashed lines correspond to first-and second-order phase transitions, respectively, while the black dots now mark the tricritical points.It is clear that the qualitative picture is the same as in the case of physical explicit breaking shown in the top panel of Fig. 3.A possible method to find the order of a phase transition is to look at the chiral susceptibility, i.e., the derivative of the order parameter with respect to the temperature χ ch = ∂∆/∂T .χ ch diverges in the case of a secondorder phase transition at the transition point.Although in a numerical calculation the divergence can never be  22 for L = 4 fm.It seems from the log-log plot that there will be no divergence for finite h N .The points can be fitted with the form a(x+c) b for x = h N /h phys N .However, the c parameter is getting smaller and smaller when the number of points is increased or the accuracy of the computation is improved.The best fit can be obtained by setting c = 0, in which case a = 16.99 and b = −0.33 is found.Furthermore, from the field equation in Eq. ( 11) it is clear that for any T < ∞ temperature and for a nonvanishing explicit breaking ϕ N = 0 cannot be a solution.Therefore, the order parameter will not vanish for finite temperatures, resulting in an increasingly sharper but still crossover transition when h N is decreased.On the other hand, we also emphasize that at h N = 0 one finds a second-order phase transition even for finite volumes.
Appendix B: Solid angle integration from a cube to a sphere In this section, Eqs. ( 35)-( 37) are derived, which are the result of a solid angle integral.Imagine a cube with side length 2λ Σ and draw a sphere with a radius R = p which is in the interval λ Σ < p < √ 3λ Σ .If p = λ Σ , the sphere is inside the cube and touches its sides, while if p = √ 3λ Σ , the sphere includes the cube and touches its vertices.In between, part of the sphere is outside the cube, while part is inside.Our goal is to calculate -for a general p in the interval mentioned -the area of the sphere that is outside of the cube.If we divide the resulting area by p 2 , we get the desired solid angle Ω = S/p 2 .This is illustrated in Fig. 23.
For λ Σ < p ≤ √ 2λ Σ there is a cap on each side of the cube.The surface of this spherical zone is S zone = 2πRh At p = √ 2λ Σ the caps touch each other.As the momentum increases, they overlap, so either the double counting must be removed or another method is needed.In both cases, one has to calculate spherical triangles and circular sectors of a spherical cap.We will calculate the contribution of modified caps that avoid double counting by dividing them into four congruent spherical triangles △ S (shown as ordinary triangles when drawn in a plane for simplicity) and four congruent sectors, as can be seen in angle covered by a spherical triangle is equal to its spherical excess, which in our case is E △ S = 2δ + β − π.The angle β can be calculated from the right triangle △ 1 on the side of the cube with vertices at the center of the side, at half of the edge of the cube, and at the point on the edge where the sphere intersects it.Its unknown leg with length g/2 can be obtained from the right triangle △ 2 , which has the same vertices on the edge of the cube, but the third one is in the center of the cube.

FIG. 2 .
FIG.2.The 1/L dependence of tree-level u and d constituent quark masses and one-loop level sigma (or f0) and pion curvature masses in vacuum (T = µq = 0).

FIG. 4 .
FIG.4.The normalized pseudocritical temperature Tc/T c,inf for characteristic sizes L > 2.5 fm with the vacuum contribution treated as finite or infinite.

FIG. 5 .
FIG.5.The equation of state, i.e., p(ε) for infinite and finite volumes of the vacuum term.

FIG. 7 .
FIG. 7. The T dependence of skewness (top) and kurtosis (bottom) for different L values at µq = µ CEP q

FIG. 8 .
FIG.8.The µq dependence of the skewness (top) and the kurtosis (bottom) at several volume sizes through the critical endpoint.The insets show the unscaled temperature dependence of these quantities.

FIG. 9 .
FIG.9.Sketch of the rearrangement of summation for APBC in two dimensions.

FIG. 10 .
FIG.10.The error made by replacing the summation with integration for the kernel E(p) with m = 300 MeV.Its decreasing trend allows the use of an UV improvement.

FIG. 11 .
FIG.11.The normalized fermion masses at L = 6 fm with APBC as a function of the grid size measured by λΣ/∆p in the case of UV improvement.

FIG. 15 .
FIG.15.ϕN as a function of the µq quark chemical potential for different temperatures for APBC.

FIG. 16 .
FIG. 16.Trajectories of the CEPs on the phase diagram with decreasing L for the LSM A, B, and C and the ePQM models with APBC.Arrows point to lower values of L.

FIG. 17 .
FIG. 17.The size dependence of T CEP and µ CEP q for the considered quark-meson models.

FIG. 18 .
FIG. 18.The path of the critical endpoints over the phase diagram for the different models with PBC.Arrows point to lower values of L

SFIG. 19 .
FIG.19.The phase diagram in the physical case and in the chiral limit, which is given by hN = 0 in the ePQM model.

FIG. 20 .
FIG.20.The critical scaling of the order parameter for the chiral symmetry in the second-order phase transition at = 0 for infinite and L = 4 fm sizes.The gray band shows scaling for the mean-field theory with critical exponent β = 1/2.

FIG. 21 .
FIG.21.The phase diagram for several finite sizes in the chiral limit.Here the dashed lines correspond to second-order phase boundary

FIG. 22 .
FIG.22.The maximum of the chiral susceptibility as a function of the explicit breaking hN /h phys N
FIG. 14. ϕN versus temperature for different values of L in the case of PBC with vacuum of infinite size.