Finite-volume effects in baryon number fluctuations around the QCD critical endpoint

We present results for the volume dependence of baryon number fluctuations in the vicinity of the (conjectured) critical endpoint of QCD. They are extracted from the nonperturbative quark propagator that is obtained as a solution to a set of truncated Dyson--Schwinger equations of ($2 + 1$)-flavor QCD in Landau gauge, which takes the backcoupling of quarks onto the Yang--Mills sector explicitly into account. This well-studied system predicts a critical endpoint at moderate temperatures and rather large chemical potential. We investigate this system at small and intermediate finite, three-dimensional, cubic volumes and study the resulting impact on baryon number fluctuations and ratios thereof up to fourth order in the region of the critical endpoint. Due to the limitations of our truncation, the results are quantitatively meaningful only outside the critical scaling region of the endpoint. We find that the fluctuations are visibly affected by the finite volume, particularly for antiperiodic boundary conditions, whereas their ratios are practically invariant.


Introduction
Fluctuations of the conserved quantities baryon number, electromagnetic charge, and strangeness are important quantities in our quest to locate a potential critical endpoint (CEP) in the phase diagram of strong-interaction matter [1]. It is one of the main goals of the experimental programs at the Relativistic Heavy-Ion Collider at Brookhaven National Laboratory to find such an endpoint in a dedicated beam energy scan [2] that probes a large baryon-chemical-potential area of the phase diagram, which is not accessible at the Large Hadron Collider. The search for the CEP is also one of the main motivations for the future Compressed Baryonic Matter experiment [3] at the Facility of Antiproton and Ion Research in Darmstadt and the Nuclotron-based Ion Collider Facility at the Joint Institute for Nuclear Research in Dubna.
While at small chemical potentials an analytic crossover from the low-temperature hadronic phase to a high-temperature quarkgluon plasma phase has been firmly established, see Refs. [4][5][6][7][8] and references therein, the location where the crossover turns into a first-order phase transition-marked by a second-order CEP-is up to now not unambiguously determined. Elaborate calculations within the functional approach to QCD using sets of coupled Dyson-Schwinger or functional renormalization group (FRG) equations [9][10][11][12][13][14] place the endpoint into the region ( CEP B , CEP ) = (495-654, 108-119) MeV, i.e., within a narrow temperature range with only moderate spread in chemical potential (see Fig. 1

below).
In experiments, signals for the CEP but also for the crossover and the first-order phase transition at large baryon chemical potential are expected to appear in ratios of event-by-event fluctuations, i.e., cumulants and related quantities. These may be particularly prominent in higher orders [15][16][17][18][19][20][21]. Ratios of cumulants have the advantage that the explicit volume dependence of the fluctuations drops out. However, as pointed out and investigated in Refs. [20,22], implicit volume dependences may remain. Whether these have to be taken into account when comparing theoretical calculations with experimental results is under debate. In any case, on systematic grounds it is interesting to study these implicit dependences. This is the topic of the present work. We study the explicit and implicit volume effects in baryon number fluctuations and ratios thereof in continuum QCD using the functional framework of Dyson-Schwinger equations (DSEs) at nonzero temperature and chemical potential. As a first step, we use the truncation scheme of Refs. [10,23] and investigate finite-volume effects for a range of cubic spatial volumes in the CEP region of the Bplane. Since this scheme does not take the effects of long-range degrees of freedom (e.g., those with the quantum numbers of the sigma meson) into account, we do not expect to be able to study the full volume dependence inside the critical scaling region. However, we believe that the scheme at hand provides a meaningful starting point for the study of fluctuations outside the critical scaling region. Furthermore, the truncation can be improved systematically in the future along the lines of the (more involved) scheme explored in Ref. [14].
We proceed as follows. First, our framework is detailed in Sec. 2. Then, we summarize the implementation of a finite spatial volume and the employed truncation and provide the reader with information how we determine fluctuations. Then, in Sec. 3, we discuss our results, and we conclude in Sec. 4.    [10][11][12][13][14] in comparison to an extrapolation from lattice QCD [24] and experimental freeze-out points [25][26][27][28][29][30]. Lower diagrams: volume dependence of the CEP of Ref. [10] as determined in Ref. [23] using antiperiodic (ABC) and periodic (PBC) spatial boundary conditions. We use a cubic volume = 3 , where denotes the edge length.

Finite-volume generalities
Introducing a finite, uniform, three-dimensional cubic volume with edge length amounts to the restriction of all appearing three-dimensional configuration-space integrals to the domain [0, ] 3 ⊂ ℝ 3 . Imposing either periodic or antiperiodic boundary conditions (PBC or ABC, respectively) in the spatial directions, each momentum-space integral is therefore replaced by a sum over discrete spatial modes according to 3 ] ⊤ and ∈ ℤ 3 , where the spatial Matsubara modes are given by ( = 1, 2, 3) When implementing a finite spatial volume in practical DSE calculations, it is (i) beneficial to rearrange the sum in Eq. (1) such that it resembles a spherical coordinate system [31]; (ii) vital to remove finite-size artifacts, which are caused by the mismatch of the discrete momentum grid and the O(3)symmetric continuum at large momenta; (iii) necessary to pay attention to the proper inclusion of the zero mode =0 = 0 that appears for PBC.
Since a detailed discussion of (i)-(iii) can be found in our previous work [23], we shall keep this section concise and refer the reader to that work for more details.

Dyson-Schwinger equations
Generally, the functional approach works with differential and/or integral equations, relating loops of nonperturbative correlation functions with each other. These can be formulated in configuration space or, more commonly, in momentum space.
As in our previous work [23], we use DSEs to investigate QCD in a finite, three-dimensional cubic volume and at nonzero temperature and quark chemical potential , where labels the flavor. The central quantity is the dressed quark propagator   that takes the form 2 with temporal Matsubara frequencies = (2 + 1) , ∈ ℤ, and = [ 1 , 2 , 3 ] ⊤ , ∈ ℤ 3 . The dressing functions , , and carry the nonperturbative information, which manifests in a nontrivial dependence on the Matsubara frequency and three-momentum.
The dressed quark propagator is the solution of its corresponding DSE that reads where 2 , , and denote the wave function renormalization constant, mass renormalization constant, and current quark mass, respectively. The former are determined in vacuum using a momentum-subtraction scheme. The self-energy is explicitly given by with the strong coupling constant , ghost renormalization con-stant˜3, and dressed gluon propagator ; two factors of i from the vertices appear explicitly such that Γ denotes the reduced dressed quark-gluon vertex. Color d.o.f. are already traced out resulting in the prefactor 4/3 for c = 3 colors. Equation (4) is shown diagrammatically in Fig. 2. In order to solve this self-consistent equation for the dressed quark propagator , we need to specify two yet unknown quantities: the dressed gluon propagator and the dressed quark-gluon vertex, which both fulfill their own DSEs. The truncation used here has been studied extensively in the past (see, e.g., Refs. [10,14,23,32,33] for 2 We work within the Matsubara (imaginary time) formalism, i.e., in Euclidean space-time with positive metric signature (++++). The gamma matrices are Hermitian and obey { , recent works and Ref. [34] for a review) and is characterized as follows.
First, in the full gluon DSE, we replace all pure Yang-Mills diagrams by fits to quenched, temperature-dependent lattice results [35,36] while the quark loop is evaluated explicitly, thereby unquenching the gluon. The resulting equation is shown in Fig. 3. As a consequence, the quark and gluon DSEs are nontrivially coupled. This results in a temperature and (implicit) chemicalpotential dependence of the gluon beyond modeling; moreover, the gluon becomes sensitive to the chiral dynamics of the quarks. Though this approximation misses implicit quark-loop effects in the Yang-Mills self-energies, these are subleading in a skeleton expansion in comparison to the quark loop. The impact of this approximate yet efficient way to compute the unquenched gluon can be estimated in vacuum within the framework of Ref. [37] and is found to be below the five-percent level. Furthermore, finite-volume effects in the quenched gluon propagator have been studied extensively at zero temperature both on the lattice and within DSEs. In both frameworks, substantial volume effects have only been found at very small momenta [38][39][40]. Since these are irrelevant for nonzero-temperature calculations due to the additional temperature scale, we work with the hypothesis that results of the present study will not be altered by the inclusion of volume effects in the gluon input in future studies.
Second, for the dressed quark-gluon vertex, we use an ansatz that supplements the leading term of the Ball-Chiu vertex [41] with a phenomenological dressing function. The latter is constructed based on a Slavnov-Taylor identity for the full vertex together with the requirement of the correct logarithmic running of the propagators in the ultraviolet momentum region. As a consequence, this construction also takes volume effects in the vertex via the quark dressing functions imposed by the Slavnov-Taylor identity into account.
For the sake of brevity and since our setup is identical to the one used in previous works, we do not show explicit expressions and refer the reader to Refs. [10,23,34] and references therein for more details. 3 The resulting set of truncated DSEs for the quark and gluon propagators is solved numerically. We use 2 + 1 quark flavors and work in the isospin-symmetric limit for the two light quarks, i.e., u = d and u = d , and choose s = 0. Thus, the baryon chemical potential is given by B = 3 u . The current quark masses are determined as follows (see Refs. [9,10] for more details): (i) the up/down quark mass is chosen such that the high-temperature behavior of the regularized quark condensate matches lattice results [43]; (ii) the strange quark mass follows from the ratio s / u = 25.7, which is obtained from in-vacuum pion and kaon masses using the Bethe-Salpeter framework of Ref. [44].
Finally, we would like to discuss the issue of the potential effects of (space-like) correlations with the quantum numbers of mesons in the quark DSE. Within the DSE framework, these correlations appear as part of a particular diagram in the DSE for the quark-gluon vertex and therefore feed back into the quark DSE. These correlations have been identified and explored in a number of works in the vacuum [45][46][47] as well as at finite temperature and chemical potential [14,48]. In Ref. [14], it was found that the inclusion of meson backcoupling effects onto the quark has only very little effect on the location of the CEP in the QCD phase diagram. On the other hand, however, these degrees of freedom will develop long-range correlations inside the critical regions around second-order transitions and will therefore be crucial when it comes to the calculation of anomalous dimensions and critical scaling.
Due to the substantial technical complications that arise when including these contributions, we ignore those in the present exploratory work. As already discussed in the introduction, we therefore do not expect to be able to study all (and likely not even the most important) contributions to volume effects inside the critical scaling region around the CEP. However, outside this region, these contributions are known to be subleading. Thus, our results should be meaningful and relevant-in particular since beyond-mean-field calculations indicate that the size of the critical region is rather small [49,50].

Quark and baryon number fluctuations
In the following, we first briefly summarize some general aspects of fluctuations (see, e.g., Refs. [1,51] for reviews) and then detail how we determine these from our solutions of the DSEs specified in the previous subsection.
In three-flavor QCD with quark chemical potentials u , d , and s , the quark number fluctuations are (dimensionless) derivatives of QCD's thermodynamic potential Ω with respect to these chemical potentials, viz.
Generally, fluctuations of conserved charges are sensitive to phase transitions and thus expected to be prime candidates to provide signatures of the (conjectured) CEP in experiments. Ratios are particularly interesting because they are equal to ratios of cumulants of the corresponding probability distributions that can be extracted from heavy-ion collisions by means of eventby-event analyses; see Refs. [1,2,51,52] and references therein for more details. Prominent ratios are the ones involving the skewness and kurtosis, namely where 2 B , B , and B denote the variance, skewness, and kurtosis of the net-baryon distribution, respectively. Analogous expressions hold for charge and strangeness.
In this work, we consider B 2 and the skewness and kurtosis ratios, which are determined via the quark number fluctuations [see Eq. (7)]. As in our previous work [10], we compute the latter from the quark number densities = − Ω/ ; for instance, In infinite volume, the regularized quark number density for a flavor is obtained from the dressed quark propagator via [10] The term sum contains a sum over the temporal Matsubara frequencies and is the one expected to yield the quark number density. However, it is divergent and needs to be regularized if-as in our approach-evaluated numerically where only a finite number of Matsubara frequencies is available. To this end, we subtract the term int that does not depend explicitly on temperature or chemical potential and is therefore known as a "vacuum contribution" [53]. This subtraction scheme is based on the Euclidean version of the contour-integration technique for Matsubara sums; see, e.g., Refs. [10,54,55] for more details. In a finite volume, the expression for turns into where we replaced the infinite-volume integral by the finitevolume sum with respect to spatial Matsubara sums in the first term sum but not in the subtraction term int . It turns out that this procedure is necessary to avoid artifacts in the resulting densities related to the high-momentum behavior of the subtraction term. The subscript ABC/PBC of the integral indicates that we treat the radial integral differently for different boundary conditions: whereas there is no momentum gap for PBC and the integral consequently starts at zero, for ABC we have ABC =0 ≠ 0 and therefore set the lower integration limit to the value of the smallest possible momentum magnitude.
Before presenting our results, we would also like to briefly comment on numerical stability. First, due to the necessary subtraction, the quark number density itself is numerically rather sensitive already in the infinite-volume setup and, as it turns  out, even more so in finite-volume calculations. Second, its derivatives, the fluctuations, are at present only accessible by means of finite-difference formulae. This leads to a limited accuracy with numerical uncertainties of the order of ten percent.

Results and discussion
Our starting point is the infinite-volume result for the location of the CEP from Ref. [10], which is displayed in the upper diagram of Fig. 1 together with recent results from other functional approaches [11][12][13][14] and an extrapolation obtained from lattice QCD [24]. As already discussed in the introduction, all results from these state-of-the-art functional calculations cluster in the region ( CEP B , CEP ) = (495-654, 108-119) MeV. The size of this region can be viewed as a measure of the lower bound for the remaining systematic error of the functional approach because the employed truncations differ in many technical details but are roughly on a similar level regarding the overall treatment of the tower of FRG equations and DSEs.
In the following, we discuss finite-volume baryon number fluctuations and their ratios around the CEP and compare them to their infinite-volume limit, which is explicitly calculated using the framework of Ref. [10] (black dot in Fig. 1). We study cubes with edge lengths of = 2.5, 3, 4, 5, 6, and 8 fm for both ABC and PBC.
We recall that the location of the CEP in the phase diagram is volume dependent-the corresponding calculations have been carried out in Ref. [23] and are shown again in the lower diagrams of Fig. 1. Sizable effects only occur for volumes ≲ (5 fm) 3 and are much larger for ABC than for PBC. For the sake of comparability, in the following we always show results obtained around the respective critical chemical potential for each system size. Additionally, we also normalize all temperatures to the corresponding critical temperatures. Since fluctuations around the CEP vary rapidly with temperature, a dense numerical grid is necessary to avoid misalignments. We accomplish this by using steps of one MeV in temperature.
In the interpretation of all results in this section, the reader should keep in mind that the present calculation is exploratory because effects from potential long-range and critical correlations (such as the sigma meson) are, as detailed in the previous section, not yet taken into account.

Baryon number fluctuations
We begin our discussion with the results for the baryon number fluctuations. In Fig. 4, we exemplarily show the second-order baryon number fluctuation B 2 for both ABC (left) and PBC (right) at finite volume as well as the infinite-volume result in black for comparison.
Before we discuss the details, a general comment is in order. As visible in both plots, the infinite-volume result does not show the expected divergence at the critical temperature. This is entirely due to our limited resolution in temperature. Given an unlimited amount of CPU time, we could determine the location of the CEP with arbitrary precision and perform calculations arbitrarily close to the CEP, thereby extracting the point of divergence exactly. In practice, for now, this is the best we could do and thus the finite width and height of the peak at infinite volume serves as our control quantity for the size of potential effects at finite volume.
Starting with the ABC results, we find clearly visible volume effects. First, we see a monotonous increase of B 2 with decreasing system size across the whole temperature range. Second, we find that the ratio of the peak height to the tails of the peaks, e.g., B 2 (1.00)/ B 2 (1.05), grows continuously with increasing system size (even within the resolution limits discussed above). The = 2.5 fm line especially stands out with the arsinh of its peak value nearly being doubled as compared to the infinitevolume result-its peak-to-tail ratio, however, is much smaller. For < CEP , we find a consistent infinite-volume limit, i.e., the results for ≳ 5 fm are very similar to one another, while the = 8 fm line coincides with the infinite-volume one. For > CEP , however, we see noticeable deviations between the finite-and infinite-volume results. We have not succeeded to  track down these deviations unambiguously, but we believe they originate from the finite-volume adjusted subtraction procedure of the density outlined in Sec. 2.3 either as remnants or as an overcompensation of the initial problem.
In contrast, the fluctuations using PBC shown in the right diagram of Fig. 4 are much less dependent on the volume of the system than the fluctuations with ABC. This resembles similar differences between PBC and ABC effects on the location of the CEP discussed above. In fact, the volume dependence of the fluctuations with PBC are within our margin of error. Nonetheless, for < CEP , one can observe a monotonous increase of B 2 with decreasing system size down to = 4 fm. For smaller box sizes, the values start decreasing again. This behavior is similar to the nonmonotonous volume dependence of the location of the CEP around = 3 fm that has been already noticed in Ref. [23]. It may be linked to the onset of the so-called epsilon regime at very small volumes [56].
In addition, the infinite-volume limit of the PBC fluctuations is much more consistent for > CEP as compared to ABC. That is, there are no substantial deviations between the = 8 fm and the → ∞ lines. As a consequence, one might conjecture that the numerical problems of ABC for > CEP are connected to infrared momentum modes because ABC introduce an effective infrared cutoff, ABC =0 ≠ 0, while PBC do not, PBC =0 = 0.
We also extracted the peak heights of the ABC and PBC results for B 2 as a function of volume and analyzed their behavior in terms of power laws. Whereas the PBC results are virtually independent of volume, i.e., roughly proportional to 0 , there is no clear power-law behavior visible for the ABC results if all volumes are considered. This may or may not be connected to the numerical error discussed at the end of Sec. 2.3 and needs to be revisited in the future in a more refined framework. Interestingly, as we will confirm in the next subsection, the volume dependence is similar for the higher-order fluctuations, too.

Ratios of baryon number fluctuations
Next, we turn to ratios of baryon number fluctuation. To this end, we display (the inverse hyperbolic sine of) both the skewness and kurtosis ratios for ABC and PBC in Fig. 5. For comparison, the infinite-volume result is again shown as a black line. In general, one can observe that both ratios for both boundary conditions and all system sizes qualitatively coincide very well for ≤ CEP . Additionally, they are also compatible with the respective infinite-volume results. Furthermore, the sign changes in both ratios are consistent with predictions based upon general grounds [17,18].
For ABC, however, we find again slight inconsistencies between the finite-and infinite-volume data for > CEP . Ne-glecting the obvious outlier at = 8 fm, the lines of the larger volumes deviate qualitatively from the → ∞ one, especially for larger temperatures. This deviation is more pronounced for the skewness ratio. Contrary to this, the PBC results exhibit once more no such behavior, and we observe a consistent infinitevolume limit also for > CEP . In fact, the = 8 fm and the infinite-volume curves are almost indistinguishable. This seems to corroborate our assumption that the infrared cutoff of ABC leads to some numerical problems for > CEP .
In addition to that, there are two notable outliers: the curves for ABC at = 8 fm and PBC at = 3 fm. Due to the randomness in their occurrence, they are most likely of purely numerical origin. We remark that the deviation from the rest of the curves in both cases occurs again for > CEP , which makes a connection to the subtraction procedure plausible. This also implies that PBC are not completely immune to these numerical problems.
Overall, we find the remarkable result that all individual volume dependences of the fluctuations cancel once ratios are studied. This not only true for large volumes but also for our smallest system sizes of = 5, 4, 3, and 2.5 fm. This is somewhat in contrast to the results of Ref. [22], where significant volume effects in the kurtosis ratio have been found for volumes ≲ (5 fm) 3 within an FRG treatment of the quark-meson model. Since the two approaches are rather different, e.g., our approach treats the gluonic sector explicitly but neglects a class of mesonic fluctuations (as discussed at the end of Sec. 2.2), and vice versa in their approach, it may be interesting to provide a systematic comparison in the future.

Summary and conclusions
In this work, we determined the volume dependence of the skewness and kurtosis ratios of baryon number fluctuations within a well-established functional approach to QCD. For a wide range of cubic spatial volumes with edge lengths between = 2.5 fm and = 8 fm and two different boundary conditions (ABC and PBC), we observe almost no volume dependence of these ratios. This is a highly nontrivial result because the individual results for the different fluctuations, both for ABC and for PBC, reveal a pattern that is at odds with the general expectation of linear dependence on volume: whereas the PBC results for B 2 do not change with volume, the ones for ABC are even inversely proportional to = 3 . Nevertheless, all these dependences cancel in the ratios, which are important when comparing with experimental results from heavy-ion collisions.
As explained in the main part of this exploratory work, our results need to be put into perspective. Mesonic degrees of freedom that are crucial to obtain the correct critical scaling inside the critical region around the CEP are not yet included. However, they have been taken into account in the infinite-volume calculation of the location of the CEP of Ref. [14]. Thus, their inclusion in our finite-volume study is in principle possible, although very demanding in terms of CPU time. This is left for future work.